Solving the where problem in neuroanatomy: a generative framework with learned mappings to register multimodal, incomplete data into a reference brain

Mapping information from different brains gathered using different modalities into a common coordinate space corresponding to a reference brain is an aspirational goal in modern neuroscience, analogous in importance to mapping genomic data to a reference genome. While brain-atlas mapping workﬂows exist for single-modality data (3D MRI or STPT image volumes), generally speaking data sets need to be combined across modalities with different contrast mechanisms and scale, in the presence of missing data as well as signals not present in the reference. This has so far been an unsolved problem. We have solved this problem in its full generality by developing and implementing a rigorous, non-parametric generative framework, that learns unknown mappings between contrast mechanisms from data and infers missing data. Our methodology permits rigorous quantiﬁcation of the local sale changes between different individual brains, which has so far been neglected. We are also able to quantitatively characterize the individual variation in shape. Our work establishes a quantitative, scalable and streamlined workﬂow for unifying a broad spectrum of multi-modal whole-brain light microscopic data volumes into a coordinate-based atlas framework, a step that is a prerequisite for large scale integration of whole brain data sets in modern neuroscience.


Summary
A current focus of research in neuroscience is to enumerate, map and annotate neuronal cell types in whole vertebrate brains using different modalities of data acquisition. A key challenge remains: can the large multiplicities of molecular anatomical data sets from many different modalities, and at widely different scales, be all assembled into a common reference space? Solving this problem is as important for modern neuroscience as mapping to reference genomes was for molecular biology. While workable brain-to-atlas mapping workflows exist for single modalities (e.g. mapping serial two photon (STP) brains to STP references) and largely for clean data, this is generally not a solved problem for mapping across contrast modalities, where data sets can be partial, and often carry signal not present in the reference brain (e.g. tracer injections). Presenting these types of anatomical data into a common reference frame for all to use is an aspirational goal for the neuroscience community. However so far this goal has been elusive due to the difficulties pointed to above and real integration is lacking.
We have solved this problem in its full generality by developing and implementing a rigorous, generative framework, that learns unknown mappings between contrast mechanisms from data and infers missing data. The key idea in the framework is to minimize the difference between synthetic image volumes and real data over function classes of non-parametric mappings, including a diffeomorphic mapping, the contrast map and locations and types of missing data/non-reference signals. The non-parametric mappings are instantiated as regularized but over-parameterized functional forms over spatial grids. A final, manual refinement step is included to ensure scientific quality of the results.
Our framework permits rigorous quantification of the local metric distortions between different individual brains, which is important for quantitative joint analysis of data gathered in multiple animals. Existing methods for atlas mapping do not provide metric quantifications and analyses of the resulting individual variations. We apply this pipeline to data modalities including various combinations of in-vivo and ex-vivo MRI, 3D STP and fMOST data sets, 2D serial histology sections including a 3D reassembly step, and brains processed for snRNAseq with tissue partially removed. Median local linear scale change with respect to a histologically processed Nissl reference brain, as measured using the Jacobian of the diffeomorphic transformations, was found to be 0.93 for STPT imaged brains (7% shrinkage) and 0.84 for fMOST imaged brains (16% shrinkage between reference brains and imaged volumes). Shrinkage between in-vivo and ex-vivo MRI for a mouse brain was found to be 0.96, and the distortion between the perfused brain and tape-cut digital sections was shown to be minimal (1.02 for Nissl histology sections). We were able to quantitatively characterize the individual variation in shape across individuals by studying variations in the tangent space of the diffeomorphic transformation around the reference brain. Based on this work we are able to establish co-variation patterns in metric distortions across the entire Introduction A current focus of research in neuroscience is to enumerate, map and annotate neuronal cell types in entire vertebrate brains, particularly mouse. To facilitate this, neuroanatomical atlases have been developed to establish predefined coordinate systems and corresponding images [2][3][4][5] . Mapping data from new experimental observations to these reference coordinates enables statistical ensembles to be built from multiple samples and multiple laboratories, as illustrated in Figure 1 (a-d). This mapping problem is generally approached through im-  Figure 1. The generative framework is based on comparing synthetic image data to real data, and minimizing the difference between these two using the Expectation Maximization (EM) algorithm. b) This includes changes to shape using the diffeomorphism model of Computational Anatomy, 3D to 2D slice estimation where necessary, intensity or color differences via nonmonotonic polynomial mappings, and modeling of signals not present in the reference image including fluorescent traces, missing tissue, or artifacts. We register data to common coordinates by calculating a penalized maximum likelihood estimator of all unknown transformations simultaneously, using a combination of the Large Deformation Diffeomorphic Metric Mapping algorithm and the EM Algorithm.
age registration, which is well established for the case of high quality samples and images acquired from the same modality 6 . In modern neuroscience heterogeneous images are produced leading to four main challenges to registration as shown in Figure 1e). Images may be 2D or 3D, and tissue may have dramatically different shapes associated to different preparations. New imaging modalities are continually being developed, and include different contrast mechanisms from reference atlases, and neuroimages include substantially different signals from reference atlases, such as experiment-specific fluorescence signals, damaged or missing tissue, or other artifacts.
We have solved this general mapping problem by developing an algorithm and computational implementation to jointly estimate unknown changes in shape, image contrast or color profile, and non-reference signal locations. Our approach is a semiparametric generative framework for atlas mapping, where a synthetic dataset is generated for any possible set of parameters. Images from one such generated dataset are shown in the left side of Figure 2a). This is complemented by an inference procedure that learns parameters to minimize the discrepancy between synthetic and observed data. The right hand side of Figure 2a) shows a comparison between synthetic and observed data. Inference is performed using a combination of the Expectation Maximization (EM) algorithm 7 and the Large Deformation Diffeomorphic Mapping (LDDMM) framework 8 , and applies to 2D serial sections or 3D volumes. In this approach, the majority of unknown parameters correspond to functions sampled on a smooth grid. These are . Registered data in the brain architecture portal. The results of the atlas-mapping framework developed in this paper were applied at scale to a large collection of multimodal whole-brain data sets of mouse. Three example datasets are shown at coarse millimeter (a,c,e) and fine (b,d,f) spatial scales, with anatomical partitions from the Allen Reference Atlas overlayed. These datasets illustrate accurate mapping despite missing tissue (a,b), non reference signals (c,d), and multiple modalities (e,f). Images in (e,f) show two modalities combined with transparency. All datasets are publicly available through the http://brainarchitecture.org web portal.
shown in Figure 2b) and include 3D change in shape, 2D slice positioning, transformations of image intensity or RGB values, and locations of any extra signal. This method allows us to combine the work of multiple different laboratories and techniques, enabling the community to benefit from the unique strengths of each.
The position of data in a reference coordinate system is not the whole story: to compute the true distribution of cells, knowledge of scale change is necessary. An illustrative example is to consider cells appearing densely packed when mapped to motor cortex of our reference image. This could be due to either a high density or a large motor cortex in our observed data. In addition to positioning experimental data into reference coordinates, our work presents important information about brain morphology derived from these mappings. We quantify local scale changes including patterns of tissue expansion or contraction characteristic of different imaging preparations. Despite being missing or not well defined in alternative approaches to mapping, we argue they are necessary for accurate density estimates in a cell census.
The large population brain mapping studies we are undertaking is providing an opportunity for a multivariate study of anatomical variability. For the first time in the mouse brain, we quantify patterns of individual variation by computing a distribution on anatomical shapes, as opposed to a univariate approach such as computing probabilities that a given structure is located at each voxel 9 . Our technique uses principal component analysis on a linear parameterization of diffeomorphisms 10,11 (see diffeomorphometry supplement for details). We use this analysis to quantify the magnitude of potential errors in cell density if local scale change is not accounted for, and to provide insight into impact of individual variation which is observed for even genetically identical mice.

Results
An end-to-end processing pipeline was built around our mapping algorithm. The pipeline was written with Matlab and Python and performs data collection, registration, quality control (QC), transformation of high resolution data, and online posting. When a brain is added to the pipeline, the program collects and downsamples all high resolution images into 15um small TIFs, which are used as inputs to our mapping algorithm. This results in deformation vector fields saved as VTK files (https://vtk.org/wp-content/ uploads/2015/04/file-formats.pdf), linear transformations as matrices, and low resolution overlay images which are used for QC. The procedure was deployed on a shared supercomputer cluster at Cold Spring Harbor Laboratory, requiring about 16 CPU threads and 24G memory for one brain to finish in 8 hours. Using parallelization, we could typically process 30 brains in 8 hours. If a registered brain passes QC, the outputs and high resolution images are sent to a multi-node custom built computer cluster to apply transformations at high resolution. This cluster contains 8 Figure 4. Local scale change is quantified from Jacobian of mappings. a) We illustrate how distortion is calculated based on scale change factors computed from diffeomorphic mappings. While a mapping transforms a reference volume to match the appearance of an observed volume (left), its Jacobian quantifies the expansion or contraction of local neighborhoods, well visualized through deforming grid elements (middle). The Jacobian's determinant quantifies local volume changes at every point, and therefore its cube root quantifies local length changes (right). Typical distortion is shown for four imaging technologies as a function of space and as histograms of all voxels inside the brain: b) Fluorescence Micro-Optical Sectioning Tomography, c) Nissl, d) Magnetic Resonance Imaging, and e) Serial Two Photon Tomography. nodes with 72 CPU threads per node, and 188G memory, and 2 Nvidia GTX 2080TI GPUs. The transformation code is tuned and highly parallel on the compute node, and can transform a terabyte size brain in 2 hours. The transformation step produces full resolution sections and deformed Atlas annotations, which are subsequently posted on brainarchitecture.org. Details of this procedure are included in the workflow supplementary material. The brainarchitecture.org portal utilizes a JPEG2000 viewer incorporating an Openlayer 4.0 (https://openlayers.org/) frontend image service and IIP backend image tiling server (https://iipimage.sourceforge.io/). This viewer provides a fully interactive zoom and pan, supports online adjustment of RGB dynamic range and contrast, as well as gamma adjustment. The viewer can simultaneously display multiple registered sections, e.g. Nissls overlaid with fluorescent tracer-labelled sections or brightfield immunihistochemically labelled sectons. The viewer shows atlas overlays at full resolution using a geojson structure, and can simultaneously show overlays of detected somata or process fragments.

Atlas mapping
Several hundred datasets registered to a reference brain and annotated using the Allen common coordinate framework (CCF) atlas using the methods in this manuscript are currently available through brainarchitecture.org. We illustrate three examples in Figure 3. corresponding to a) Nissl images with tissue removed for snRNA-seq 12, 13 , c) two photon tomography, e) viral tracing. Viewing higher resolution resolution data (b,d,f) shows experimental data aligned to annotations, and the web interface provides name and contours (red) upon mouse over. One important application of this method was to register the Allen atlas to in vivo MRI. This allowed us to identify the bregma location on the skull, which is absent from the Allen atlas, and provided as standardized origin to the coordinate space (see Supplemental data).

Scale change
Our procedure for calculating scale change is illustrated in Figure 4a), which shows that as diffeomorphic mappings change the position of points on a gridded coordinate system, they also change the volume of each region in the gridded coor- To build a multivariate model across voxels this matrix is factored e) to yield uncorrelated and orthogonal modes of variation. f) Diffeomorphisms are reconstructed from initial velocities, and their local scale change is computed for visualization and analysis. Bottom: g) The first few principal modes describe the majority of individual variability in mouse anatomy beyond that explained by age. h) The mean change from Allen CCF, effect of ageing (for 7 weeks from the mean), and largest mode of variability (sampled at 1.5 standard deviations from the mean) are illustrated in terms of scale change of the reconstructed diffeomorphism.
dinate system. The change is characterized by the Jacobian (the matrix of partial derivatives) of the mapping between the individual brain and the reference brain coordinates. The corresponding change in volume is characterized by its determinant. To make numbers more interpretable, we work with the cube root of the Jacobian determinant, corresponding to the geometric average of local scale changes at a point. Details are provided in the diffeomorphometry supplement.
We calculated scale change at every voxel for typical images from several different modalities, slices of which are shown in Figure 4b). These correspond to tissue distortion introduced by processing and imaging, and must be understood correctly to report accurate densities of neurons or other markers. For STPT data, we found a scale change relative to the CCF reference of 0.93 (7% shrinkage). It is to be noted that the Allen CCF, which is also gathered using the STPT technique, was originally artificially scaled up by about 7% to account for tissue shrinkage compared to the in-vivo brain. For fluorescence micro optical sectioning tomography, we found a scale change relative to the CCF of 0.84. For our Nissl histology sections, the scale change relative to CCF was 1.02. Using in vivo and ex vivo MRI images, we quantified the scale change between in-vivo and ex-vivo MRI for a mouse brain to be 0.96. While these overall scale changes are informative, note that we obtain the local scale factor at each image location, and these scale changes vary significantly in brain space (see histograms in Fig. 4, bottom row), depending on the data acquisition technique, and also on individual biological variation. The spatial non-uniformity of these scale changes points to the importance of the quantification of local 6/13 scale changes in doing quantitative neuroanatomical analysis. Our methodology for the first time permits such quantitative analysis in a mathematically rigorous as well as algorithmically robust manner.

Individual variation
Our multivariate approach to quantifying individual variation is illustrated in Figure 5 a-f). Because deformations (a) do not lie in a linear space (the sum or differences of two invertible mappings is not necessarily invertible), we work with velocity fields (b) which are a linear parameterization from which deformations can be reconstructed. Note that the metric tensor is positive definite in our case and prima facie cannot be modeled as having Gaussian fluctuations -in other words we cannot directly subject the length changes to PCA style analysis in a statistically meaningful manner. Each of the x, y, z components of velocity vectors at each voxel are stacked into a data matrix (c), from which univariate or bivariate variation can be quantified (d). The histograms show univariate distributions of individual components of the velocity field at fixed locations, whereas the density plots show bivariate distributions across two spatial locations. It can be seen that in general the velocity fields are not independent, so that two brain regions may either be correlated or anticorrelated in their fluctuations. Kernel principal component analysis is applied (e) to uncover the uncorrelated modes of variability in a multivariate Gaussian distribution. Reconstructing deformations from these modes (f) yields a non-Gaussian distribution on the nonlinear space of scale changes.
We examined a data set consisting of a population of 113 STPT volumes of C57BL/6 mice, with 47% female and mean age 10.3 weeks (standard deviation 3.28 weeks). We performed linear regression between age, sex and initial velocity at every voxel, modeling the observed initial velocity at voxel i for brain j by v j 0 (x i ) = a i + b i age j + c i sex j + noise i, j where age is measured in weeks, sex = 1 for female mice and 0 for male. The noise is assumed 0 mean Gaussian and colored based on the smoothness of velocity fields (i.e. Lv is modeled as white noise for L defined in methods section (4)).
The unknown parameters a i , b i , c i are estimated by maximum likelihood at each voxel location x i . We performed permutation testing with 10,000 permutations, using the norm of the residuals across all voxels to generate an F test statistic for comparing nested models. We found a significant relationship with age (we reject b i = 0 for all i with p < 1e − 4) but not sex (we fail to reject c i = 0 for all i with p = 0.0673 ), and proceeded with kernel principal component analysis on the residuals after regressing for age. As illustrated in Figure 5g), we found that the majority of signal power in our dataset is explained by the mean (parameter a), corresponding to the average difference in shape between the Allen atlas and an STPT image in our dataset. The effect of age (parameter b) accounts for 24.2% of the variance (i.e. expected norm squared of the residual) in our dataset. We found that a small number of modes of variability accounted for the majority of variance not explained by age, with the first principal component accounting for 9.47%, and the first 3 for 23.6%.
We illustrate the effect of these components by reconstructing a diffeomorphism from our linear parameterization, and visualizing the corresponding scale change superimposed upon the Allen atlas. The mean term (parameter a) is illustrated in the left column of Fig. 5h, and shows considerable differences between our samples and the Allen CCF. For example, the hippocampus, thalamus, and caudoputamen tend to be relatively larger in our samples. The effect of age (parameter b) is illustrated in the center column of Fig. 5h, and shows that some structures like the midbrain grow relative to others like the thalamus. In the supplementary materials we quantify this pattern throughout the brain by computing mean values of log scale change, as well as mean displacement magnitude, in brain regions of the Allen atlas ontology. For non cortical structures we use level 5 of the Allen atlas ontology (83 structures), and for cortical structures we consider only those in level 7 with "cortical plate" as a parent (31 structures).
We illustrate the first mode of variability at 1.5 standard deviations in the right column of Fig. 5h. Individual variation includes scale changes typically as high as 1.25x, associated with displacements inside the brain with a median of 28.2 microns, and 95th percentile of 139 microns. Visually, the first mode can be seen to be left right asymmetric, and includes strong correlations between thalamic, striatal, and superficial cortical structures, and anticorrelations with midbrain. We quantify this pattern throughout the brain in supplementary material as described above.

Discussion
The methods presented in this work represent three major contributions to the field of modern neuroanatomy. First, we solve the problem of mapping across imaging modalities, in the presence of partial data or non-reference signals, using a generative atlas mapping framework. We make all our mapped data publicly available through a web portal (http://brainarchitecture.org). Second, we present estimates of local scale change which we emphasize are necessary for calculating densities in reference coordinates. Third, we present for the first time a rigorous multivariate statistical framework that can deal with the non-Gaussian, non-linear statistical deformation fields underlying the diffeomorphic neuroanatomical model, and apply it to a collection of mouse brains to derive biologically meaningful results.
This work builds upon a long history of brain atlasing. The notion and utility of parcellating the brain into regions of similar cytoarchitectonic features was established by Brodmann 14 Economo 15 and others. For reasons of reproducibility and stereotactic surgery, these regions were augmented by a standard coordinate system based on the anterior and posterior commissure 16,17 , which has since been refined and modified.

7/13
For mouse neuroanatomical studies, coordinate systems indexed to the skull's bregma and lambda positions have been developed. The histology based Allen Reference Atlas 2 , later augmented to include 3D serial two photon imaging from multiple animals is used in this work and includes 671 annotated structures at the leaves of its ontology. Other atlases are commonly used by the community, including the Waxholm space MRI atlas 3 , which has one advantage of being in vivo rather than excised. One challenge in modern brain atlasing is the representation of variability, and while univariate methods such as probabilistic atlas construction have been well established 9 , multivariate methods are less so.
The nonrigid alignment of brain imaging data to common coordinate systems is a widely studied problem, with many established algorithms and software packages. In a head to head comparison of 14 methods 18 , those with a large number of degrees of freedom in the spatial transformation (often referred to as "freeform", or "fluid", or "elastic") such as SyN 19 , IRTK 20 , ART 21 , outperformed those with fewer degrees of freedom (such as polynomials or low order splines). These algorithms function by maximizing a similarity measure between a deformed reference image and a target observation with some imposed regularity, and those with a similarity that assumed a correspondence between atlas and target voxel intensities (i.e. mean square error) performed more poorly than those which did not. The approach we present here exploits both these findings, building fluid deformations using the LDDMM framework 8 , and iteratively estimating a correspondence between voxel intensities.
Modern research in brain to atlas registration has focussed on registering much more challenging datasets than classically treated. Our approach addresses three common issues: registration between 3D images and 2D slices, and in the presence of extra signals or artifacts, and between images from different modalities. A review of slice to volume image registration methods was recently compiled 22 , which identified a paucity of nonrigid approaches to the problem. Some approaches 23,24 first identify a slice plane using rigid transformations and follow this with nonrigid 2D registration on each slice, avoiding modeling potential out of plane deformations. In contrast, the technique we present here allows for arbitrary 3D or 2D deformations.
Working with histology sections in particular leads to specific challenges as reviewed in 25 . Signal not present in reference images, such as tears or holes, are described as important, but are not addressed adequately by the surveyed methods. In contexts other than histology, extra signal such as fluorescence from tracers presents a similar problem. Previous authors have approached these issues by masking out these regions 26,27 or filling them with specific intensities or textures 28 before registration, often manually. Other approaches have attempted to jointly estimate these signals during registration, such as with excised tissue 29 , tuberculosis 30 , or tumor 31 .
When these extra signals are absent, registration between different modalities has been applied successfully by using appropriate similarity functions such as local cross correlation 19 , mutual information 32, 33 , or others 34,35 . We follow a generative approach, where the appearance of target images are synthesized as a transformation of the atlas. This idea has been applied with high quality data 24,36 , but the multi modality registration problem is more challenging when extra signals are present. Related work 37,38 has approached this challenge, but have been restricted to simple deformations. The method we develop jointly addresses each of these challenges: estimating arbitrary 3D and 2D deformations, predicting locations of extra signal by treating them as missing data in an Expectation Maximization setting as in [38][39][40] , and modeling different contrasts by synthesizing them in a generative framework.
The registration results presented illustrate an important use of our method: associating experimental observations in a mouse brain to a particular region or layer in the Allen atlas. We analyzed datasets important to the ongoing Brain Initiative Cell Census Network, where registration to atlas coordinates could not be accurately performed using existing methods. For single neuron RNA sequencing data, 41 large regions of tissue were excised for analysis. These slices often contained different contrast profiles and large spatial distortions as described in 42 . Our registration methods are enabling accurate positioning of the snRNAseq results into the reference brain. For the serial two photon tomography data shown 43 , strong fluorescence signals have traditionally limited registration accuracy in exactly the regions where it is most important. Our Nissl and fluorescence data shows registration of two preparations simultaneously, and demonstrates accurate positioning for tracing and connectivity studies.
The study of local scale changes within the brain using the derivatives of mappings have become known as tensor based, or voxel based morphometry 44 . These techniques have been used extensively in studying the human brain, where contraction can be associated with tissue loss in neurodegenerative disease. Our local scale change results demonstrate a large difference in tissue for different imaging preparations. Despite a significant amount of handling, serial Nissl sections showed the least amount of distortion, whereas FMOST showed the most, relative to the Allen common coordinate system. Our in vivo to ex vivo MRI registration demonstrated tissue shrinkage of 4%, a factor that is relevant to other ex vivo imaging modalities.
The random orbit model of Computational Anatomy 45 provides a powerful approach to the quantitative study of shape and form. In this paradigm, the spatial transformations generated from registration approaches are treated as objects of study. Focusing modeling efforts on transformations enables a multivariate description of individual variation, rather than univariate probabilistic atlases 9 . Because diffeomorphisms and their derivatives belong to a nonlinear space (i.e. the sum or difference of two diffeomorphisms need not be a diffeomorphism) linear representations 46 have been used for Gaussian modeling 10 . In human data, we have used these models to build low dimensional empirical priors similarly to Active Shapes 47 that improve registration in the presence of noise 11,48 , and have quantified complexity in large populations in terms of these priors 49 . In this work we apply them to mouse neuroanatomy for the first time.
While human brains show marked variation in terms of gyral and sulcal patterns, our individual variation results show a surprising degree of variability in a standard population of mouse brains despite their lack of gyrification. Typical fluctuations in scale as high as 1.25x reveal the importance of scale change calculations even within the same imaging modality. Importantly, change in scale associated to individual variation is often larger than differences between different modalities.
One limitation of our approach is the required computation time. Estimating, deformations by integrating a time varying velocity field increases computation time relative to a single displacement field proportional to the number of timesteps in the discretization, and Expectation Maximization algorithms may also be slow to converge. As such we have put effort into into building an infrastructure for high performance and parallel computing (see workflow supplement). For example, when studying large populations, throughput can be increased by analyzing each brain on a different cluster node.
An important open question in the image registration field is one of standardizing accuracy. Images transformed tangent to their level appear identical 50 , so there are generally an infinite number of transformations with the same accuracy and thus regularity must be considered as well. The evaluation in 18 presented accuracy in terms of segmented structures, and provided 5 different measures defined for both 3D segmentations and 2D surfaces: target overlap, mean overlap, union overlap, false negative and false positive errors, as well as similarity in the volume of segmented structures and distance between their boundaries. The evaluation in 51 defined accuracy in terms of feature points identified in both atlas and target images, and presented several possible weighted averages or medians to combine accuracy of each point within a dataset. For 3D to 2D registration, accuracy is more difficult to define, as many standard overlap measures require sets of the same dimension and the same feature points can rarely be accurately placed in both image sets. When data is missing, quantifying accuracy becomes impossible even in principal. More effort is required within the registration community to standardize accuracy measurements for specific imaging tasks. When the task is annotation of observed images, accuracy can often be judged by eye as acceptable or unacceptable. Traditionally, any errors that occur would be corrected manually by an anatomist. However, keeping corrected annotations consistent with mappings is essential for accurate scale change measurements and variability quantification, and we have developed a refinement approach to update mappings based on manual labeling as described in supplemental material. Our code and data format specification are also made available in supplementary material.
A key implication of the method we have developed is that it provides automatic registration capabilities for even low quality or damaged images. This is particularly important for experimental imaging methods whose protocols have not been fine tuned, and recognizes the value of each tissue specimen even if damaged. Our multi modality registration approach is capitalizing on the trend of generative models to understand complex relationships in imaging. While the approach described in 24 synthesizes image contrasts, it adopts a discriminative objective function, rather than the purely generative log likelihood used here. As diverse data from more laboratories and experimental modalities are registered to reference coordinates and combined within community brain atlases, their utility and impact on neuroscience will continue to grow -and this increasing utility can be facilitated using the techniques presented in this paper.

Generation of synthetic datasets
Our method performs registration by producing synthetic datasets with shape and intensity or color profile that closely matches that of observed images. Specifically for this paper, synthetic data was generated as a sequence of transformations of the Allen Institute's common coordinate framework reference images for the mouse brain gathered using STPT. These transformations are illustrated in Figure 2 and can be grouped as (i) 3D shape changes (ii) 3D location, scale, or orientation changes (iii) 2D shape changes (if applicable) (iv) 2D location orientation changes (if applicable) (v) image intensity changes.
Transformations that encode shape are formulated using the Large Deformation Diffeomorphic Metric Mapping paradigm 8 (i: ϕ 0 , iii: ϕ i for the i-th slice), and parameterized in the tangent space of the diffeomorphism group via flows of smooth velocity fields v t (x) specified on a voxel grid from time t = 0 to 1: Other spatial transformations are formulated as affine transformation matrices (ii: affine A, iv: rigid R i for the i-th slice), parameterized in their tangent space using the matrix exponential (e.g. the rigid transformation group is parameterized by the antisymmetric matrices). Image intensity changes are modeled as 3rd degree polynomials of the pixel intensities or colors (v: f i (I) ), which is the only parametric component of our method, and provides enough degrees of freedom to permute the brightness of gray matter, white matter, and background into any observed order.
Given the i-th observed image J i , we transform our atlas image I to synthesize its appearance: J i (x) = f i (I(T i,−1 (x))) for T i,−1 the inverse of T i .

9/13
For 3D images, R i , ϕ i are set to identity and i = 1, and for 2D images with rigid motions only, such as acquired through tape transfer techniques 12,13 , ϕ i are set to identity. The transformations ϕ 0 and A contain global information about shape to be analyzed, whereas R i and ϕ i contain nuisance information that is inconsistent from slice to slice.

Inference and the EM-LDDMM likelihood function
Our generative formulation produces an expected target image, and differences from this expectation are well modeled as Gaussian noise. Estimation of the unknown transformation parameters are therefore posed as penalized maximum likelihood. For an atlas image I, a family of observed images J i we learn transformation parameters by solving the optimization problem: arg min The Reg terms provides necessary regularization for the infinite dimensional diffeomorphism group, and are parameterized by a spatial smoothness scale α, and a magnitude σ 2 R .
Regularization is defined through an inner product and thus metric on the tangent space of the diffeomorphism group. This inner product is large when velocity fields contain quickly varying (high spatial frequency) components, encouraging solutions to be smooth, and guaranteeing that the integral of Eq. (1) is a diffeomorphism. Any minimizer of (3) defines a shortest length geodesic curve 46 . The square error term corresponds to a negative log likelihood of Gaussian noise with meanĴ i and variance σ 2 .
This statistical interpretation allows us to accommodate images with non reference signals, such as missing tissue, tracer injection sites, or other anomalies. At each pixel, the identity of the signal type is modeled as missing data, and maximum likelihood estimators are computed using an Expectation Maximization algorithm, which alternates between the E step: compute posterior probability π i (x) that each pixel corresponds to the reference image rather than one of the nonreference types, and the M step: update parameters by solving a posterior weighted version of the above: arg min As an EM algorithm, this approach is guaranteed to be monotonically increasing in likelihood. An example of posterior weights are shown in the right hand column of Figure 2b.

Nonconvex optimization with low to high dimensional subgroups and resolutions
This registration problem is highly nonconvex, and allows for many local minima. To provide robustness in our solution, we solve a sequence of lower dimensional subproblems, initializing the next with the solution to the previous. (i) 2D slice translation by center of mass detection (ii) 2D slice to slice rigid alignment maximizing similarity to neighbors 52 (iii) 3D affine only alignment, registration using the full model at (iv) low (200 micron), (v) medium (100 micron), and (vi) high (50 micron) resolution. Time varying velocity fields are discretized into 5 timesteps and integrated using the Semi Lagrangian method 53 . For most subproblems, spatial transformation parameters are estimated by gradient descent, and intensity transformation parameters are updated by solving a weighted least squares solution at each iteration. For subproblems that include linear registration only, parameters are estimated using a second order Gauss-Newton optimization scheme.

Local scale change is necessary for cell density estimation
Because shape change may be expansive or contractive, local scale change must be calculated and stored to accurately estimate densities in reference coordinates. To correctly account for these scale changes we must estimate the determinant of Jacobian (matrix of partial derivatives) of our 3D mappings at each voxel. The diffeomorphism theory of our approach guarantees nondegeneracy 54 unlike spline based methods, and also leads to favorable performance in the presence of noise 50 . For more interpretable units, we work with the cube root of the determinant of Jacobian, which corresponds to length change rather than volume change. This data is calculated at every pixel and stored in reference coordinates, and also in coordinates of observed images. An illustration of this approach is shown Figure 4a).

Individual anatomical variation
The tangent space parameterization of the diffeomorphism group provides a vector space representation of anatomical shape, complete with an appropriate inner product that allows for linear modeling of nonlinear shape data 10,11 . Because solutions of Equation (5) are geodesics, the entire trajectory v t can be reconstructed from v 0 46 , a tangent vector to the identity element of the diffeomorphism group. Our method enables us to calculate and present the orthogonal modes of variation that describe individual variability of mouse neuroanatomy, using kernel principal component analysis.
For a reference image with M voxels, and N observations, we let X be the 3M × N matrix of initial velocity vector fields at each voxel, X 0 be centered based on a mean and age regression, and let K be the 3M × 3M kernel matrix defining 10/13 inner products between smooth vector fields consistent with regularization in Eq. (4). We compute the eigendecomposition of the inner product matrix (Gram matrix) G = X T 0 KX 0 . Its eigenvectors are linearly related to modes of variability, and eigenvalues are proportional to their respective variances. Importantly, each mode is orthogonal with respect to an inner product on the space of smooth functions appropriate for our application. An illustration of this approach is shown Figure  5a-f).