# Cross-modality mapping using image varifolds to align tissue-scale atlases to molecular-scale measures with application to 2D brain sections

Kaitlin M Stouffer, Alain Trouvé, Laurent Younes, Michael Kunst, Lydia Ng, Hongkui Zeng, Manjari Anant, Jean Fan, Yongsoo Kim, Xiaoyin Chen, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Kaitlin M Stouffer, Alain Trouvé, Laurent Younes, Michael Kunst, Lydia Ng, et al. Cross-modality mapping using image varifolds to align tissue-scale atlases to molecular-scale measures with application to 2D brain sections. Nature Communications, 2024, 15 (1), pp.3530. 10.1038/s41467-024-47883-4. hal-04583184

HAL Id: hal-04583184

https://u-paris.hal.science/hal-04583184

Submitted on 22 May 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Cross-modality mapping using image varifolds to align tissue-scale atlases to molecular-scale measures with application to 2D brain sections

Received: 13 April 2023

Accepted: 10 April 2024

Published online: 25 April 2024

<figure><img src="image_2.png" /></figure>

Check for updates

Kaitlin M. Stouffer Ⓔ<fnref n="1" />,²,³✉, Alain Trouvé Ⓔ³, Laurent Younes Ⓔ⁴, Michael Kunst⁵, Lydia Ng⁵, Hongkui Zeng Ⓔ⁵, Manjari Anant<fnref n="1" />, Jean Fan Ⓔ<fnref n="1" />, Yongsoo Kim Ⓔ⁶, Xiaoyin Chen Ⓔ⁵, Mara Rue Ⓔ⁵ & Michael I. Miller Ⓔ<fnref n="1" />,²✉

This paper explicates a solution to building correspondences between molecular-scale transcriptomics and tissue-scale atlases. This problem arises in atlas construction and cross-specimen/technology alignment where specimens per emerging technology remain sparse and conventional image representations cannot efficiently model the high dimensions from subcellular detection of thousands of genes. We address these challenges by representing spatial transcriptomics data as generalized functions encoding position and high-dimensional feature (gene, cell type) identity. We map onto low-dimensional atlas ontologies by modeling regions as homogeneous random fields with unknown transcriptomic feature distribution. We solve simultaneously for the minimizing geodesic diffeomorphism of coordinates through LDDMM and for these latent feature densities. We map tissue-scale mouse brain atlases to gene-based and cell-based transcriptomics data from MERFISH and BARseq technologies and to histopathology and cross-species atlases to illustrate integration of diverse molecular and cellular datasets into a single coordinate system as a means of comparison and further atlas construction.

Since the 17th century, scientists have seen living organisms as a hierarchy of biological mechanisms at work across scales. To understand the interplay of these mechanisms, reference atlases that incorporate genetic, cellular, and connectivity measures into a single coordinate space have been constructed and which aim to summarize the mass of data across scales through a set of discrete partitions. An instance of the more general segmentation problem in computer vision, atlas construction relies on the underlying assumption of homogeneity within each region. The optimal partitioning assigns a label to each region based on this homogeneity and the presence of sharp changes at the boundaries between regions. In biology, this label

frequently reflects behavior or function, as seen in two of the most common mouse brain atlases: the Allen Common Coordinate Framework (CCFv3)<fnref n="1" /> and the Franklin and Paxinos Atlas². Together, the common coordinate framework an atlas provides in addition to its ontology facilitates the comparison of different types or replicates of data in a single coordinate system and hones efforts of study to particular regions relevant to each unique investigation.

Hence, atlas construction and segmentation of new data with defined atlases comprise two classes of widespread problems across domains of biology. These are particularly relevant in the context of emerging datasets in spatial transcriptomics, where atlases have still

<footnote n="1">¹Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD, USA. ²Kavli Neuroscience Discovery Institute, Johns Hopkins University, Baltimore, MD, USA. ³Centre Borelli, ENS Paris-Saclay, Gif-sur-Yvette, France. ⁴Department of Applied Mathematics and Statistics, Johns Hopkins University, Baltimore, MD, USA. ⁵Allen Institute for Brain Science, Seattle, WA, USA. ⁶Department of Neural and Behavioral Sciences, Penn State University, College of Medicine, State College, PA, USA. ✉e-mail: kstouff4@jhmi.edu; mim@jhu.edu</footnote>
<figure><img src="image_3.png" /></figure>

# Cross-modality mapping using image varifolds to align tissue-scale atlases to molecular-scale measures with application to 2D brain sections

Received: 13 April 2023

Accepted: 10 April 2024

Published online: 25 April 2024

<figure><img src="image_4.png" /></figure>

Check for updates

Kaitlin M. Stouffer Ⓔ<fnref n="1" />,²,³✉, Alain Trouvé Ⓔ³, Laurent Younes Ⓔ⁴, Michael Kunst⁵, Lydia Ng⁵, Hongkui Zeng Ⓔ⁵, Manjari Anant<fnref n="1" />, Jean Fan Ⓔ<fnref n="1" />, Yongsoo Kim Ⓔ⁶, Xiaoyin Chen Ⓔ⁵, Mara Rue Ⓔ⁵ & Michael I. Miller Ⓔ<fnref n="1" />,²✉

This paper explicates a solution to building correspondences between molecular-scale transcriptomics and tissue-scale atlases. This problem arises in atlas construction and cross-specimen/technology alignment where specimens per emerging technology remain sparse and conventional image representations cannot efficiently model the high dimensions from subcellular detection of thousands of genes. We address these challenges by representing spatial transcriptomics data as generalized functions encoding position and high-dimensional feature (gene, cell type) identity. We map onto low-dimensional atlas ontologies by modeling regions as homogeneous random fields with unknown transcriptomic feature distribution. We solve simultaneously for the minimizing geodesic diffeomorphism of coordinates through LDDMM and for these latent feature densities. We map tissue-scale mouse brain atlases to gene-based and cell-based transcriptomics data from MERFISH and BARseq technologies and to histopathology and cross-species atlases to illustrate integration of diverse molecular and cellular datasets into a single coordinate system as a means of comparison and further atlas construction.

Since the 17th century, scientists have seen living organisms as a hierarchy of biological mechanisms at work across scales. To understand the interplay of these mechanisms, reference atlases that incorporate genetic, cellular, and connectivity measures into a single coordinate space have been constructed and which aim to summarize the mass of data across scales through a set of discrete partitions. An instance of the more general segmentation problem in computer vision, atlas construction relies on the underlying assumption of homogeneity within each region. The optimal partitioning assigns a label to each region based on this homogeneity and the presence of sharp changes at the boundaries between regions. In biology, this label

frequently reflects behavior or function, as seen in two of the most common mouse brain atlases: the Allen Common Coordinate Framework (CCFv3)<fnref n="1" /> and the Franklin and Paxinos Atlas². Together, the common coordinate framework an atlas provides in addition to its ontology facilitates the comparison of different types or replicates of data in a single coordinate system and hones efforts of study to particular regions relevant to each unique investigation.

Hence, atlas construction and segmentation of new data with defined atlases comprise two classes of widespread problems across domains of biology. These are particularly relevant in the context of emerging datasets in spatial transcriptomics, where atlases have still

<footnote n="1">¹Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD, USA. ²Kavli Neuroscience Discovery Institute, Johns Hopkins University, Baltimore, MD, USA. ³Centre Borelli, ENS Paris-Saclay, Gif-sur-Yvette, France. ⁴Department of Applied Mathematics and Statistics, Johns Hopkins University, Baltimore, MD, USA. ⁵Allen Institute for Brain Science, Seattle, WA, USA. ⁶Department of Neural and Behavioral Sciences, Penn State University, College of Medicine, State College, PA, USA. ✉e-mail: kstouff4@jhmi.edu; mim@jhu.edu</footnote>
yet to be constructed³ and where the different technologies used number quite high, thus emphasizing the need for mapping data across samples and technologies into a common coordinate framework for comparison. A central challenge in both atlas construction and segmentation, however, is in aligning in a single coordinate space data that often exists at different scales (e.g. coarse scale atlases to fine-scale datasets) and is fundamentally of different functional modalities (e.g. tissue regions to molecular images). These differences preclude an obvious definition of similarity for consequently optimizing alignment.

In the classical image setting, a large family of diffeomorphism-based methods⁴ has been developed within the field of Computational Anatomy (CA)⁵,⁶ for transforming coordinate systems at the tissue scales. The smoothness, invertibility, and non-rigid nature of diffeomorphisms are particularly relevant in medical image alignment as they reflect the mechanics/dynamics of soft tissue⁷. Consequently, methods rooted in diffeomorphic registration come particularly from multiple labs in the magnetic resonance imaging (MRI) community, with variations developed in the contexts of both atlas construction and segmentation at these tissue scales⁷⁻¹⁶. Large Deformation Diffeomorphic Metric Mapping (LDDMM)¹⁷, on which many such methods are built, specifically equips diffeomorphisms with a metric that allows them to be reduced to a low-dimensional representation that conveniently identifies shape. Consequently, it has been harnessed in the setting of image classification and the study of atrophy in diseased versus control populations¹⁷⁻¹⁹. While successful at aligning images of the same modality, most of these diffeomorphism-based methods focus exclusively on generating geometric transformations (e.g. diffeomorphisms). The challenge of crossing functional modalities has instead been addressed in a classical imaging setting by coupling such methods to variations in matching cost or additional transformations. These include matching based on analytical methods using cross-correlation¹³ or localized texture features²⁰, and methods for transforming one range space to another in crossing modalities and scales based on polynomial transformations²¹, scattering transforms²², graphical models²³, and neural networks²⁴⁻²⁶.

In the molecular setting, however, both the diversity and magnitude of data measured by spatial transcriptomics technologies²⁷ often prohibit the representation of such data as classical continuous images discretized as regular grids, and consequently, the direct use of these image alignment methods at the molecular scales. As seen in repositories generated in the BRAIN Initiative Cell Census Network (BICCN) and archived at the Brain Image Library (BIL), these datasets are already on the order of terabytes and will only continue to increase as technologies shift from mouse to human measurements. Hence, while some methods in deep learning have been applied to align single-cell datasets, modeled as regular grid images, both to atlases at the histological scale²⁸ as well as reference transcriptional atlases that are beginning to emerge²⁹, most image-based tactics are limited in their ability to represent this high-dimensional and memory-intensive data. Furthermore, many learning-based schemes in the context of classical images have relied on the use of extensive numbers (typically on the order of hundreds) of images from different subjects for extracting features useful for cross-modalities and mapping¹⁸,²⁶. In the setting of spatial transcriptomics technologies that have been developed only over the last few years, the acquisition of such large training datasets is typically prohibitive for some of these learning-based approaches.

Consequently, an independent class of methods has been developing for aligning spatially resolved transcriptional datasets at these molecular scales. With influence from image-based methods, some of these including GPSA³ and PASTE³⁰ focus exclusively on spot-based data in which gene expression is measured in a neighborhood of each spot for a regular array of spots, analogous to a voxel grid in an image. In contrast, image-based technologies such as STARmap, BARseq, and MERFISH, which take point measures of individual mRNA molecules or

cells, cover measurements irregularly sampled over space according to tissue architecture and dynamics. Furthermore, natural fluctuation in gene expression over time and space coupled to the dynamics of each spatial transcriptomics technology leads each tissue section, at the molecular/cellular (0.3–100 micron) scales, to have a varying number of such particles with no natural ordering of particles consistently apparent between sections. Consequently, landmark-based methods³¹ that assume direct permutation correspondence between particles are not applicable. Methods aimed at generalizing to allow alignment within and across these additional technologies are typically rooted in different data representations, such as graphs, in the case of SLAT³² that aims to find correspondence between cells or groups of cells in an atlas and target or CAST³³, a deep graph-based neural network (GNN) algorithm which learns a graph representation of single-cell resolution omics data and subsequently aims to align datasets via these graph representations. Additionally, while many methods have showcased success at aligning different replicates within a single technology for both human and mouse samples, cross-technology alignment has typically relied on a nonempty intersection of feature sets across these technologies²⁸,³²,³³, limiting extension of these methods to the type of cross-modality and cross-scale mapping required for integrating tissue-scale atlases and molecular scale data.

Here, we build on both of these lineages in presenting a model equally equipped at representing tissue-scale atlases and molecular-scale data, and an associated cross-modality registration mechanism that addresses both challenges of crossing scales and functional modalities by estimating simultaneously geometric and functional transformations to align each dataset to the other. Specifically, we harness the generalizability of image varifolds, which have emerged in molecular CA³⁴, for simultaneously modeling molecular and tissue-scale data with both irregular and gridlike sampling schemes in a common framework (addressed in SLAT and CAST, for instance, by a graph-based representation in the setting of molecular and cellular data³²,³³). A subproblem covered by the image-varifold theory outlined in ref. 36 is the mapping of molecular scale data to atlas coordinate systems. As specified there, we estimate minimal energy diffeomorphic transformations as in image-based LDDMM¹⁷, but where the action of diffeomorphisms on images has been adapted to the setting, here, of image varifolds in a consistent manner³⁴. Additionally, we estimate in tandem a latent distribution over the molecular functional space for each atlas partition, without the need for the large datasets often used in deep learning approaches, but by instead relying on an assumption of spatial homogeneity in distribution across each atlas partition. We use this latent distribution to transform the functional space of each atlas (i.e., its ontology) to the molecular functional space. We refer to this method of jointly estimating diffeomorphisms and latent feature distributions, as Cross Image-varifold LDDMM (xIV-LDDMM), where the cross emphasizes its ability to map across scales and modalities³⁵. We demonstrate this methodology with two common implementations of image varifolds (triangled meshes and point clouds). First, we introduce the methodology through a demonstration of mapping both gene-based and cell-based MERFISH datasets to corresponding sections of the Allen Common Coordinate Framework (CCFv3)¹. Second, we quantify accuracy in mapping specifically characterized cell types in BARseq data³⁶ to the CCFv3. Finally, we illustrate the method's generalizability to additional modalities and offer diverse examples of its use across the spans of biology in atlas construction, cross-replicate and cross-species comparison, and data segmentation.

## Results

### Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms

To accommodate the high feature dimensionality and spatial irregularity of molecular datasets, as described in the Introduction, we harness the recent work that extends the theory of diffeomorphisms and
yet to be constructed³ and where the different technologies used number quite high, thus emphasizing the need for mapping data across samples and technologies into a common coordinate framework for comparison. A central challenge in both atlas construction and segmentation, however, is in aligning in a single coordinate space data that often exists at different scales (e.g. coarse scale atlases to fine-scale datasets) and is fundamentally of different functional modalities (e.g. tissue regions to molecular images). These differences preclude an obvious definition of similarity for consequently optimizing alignment.

In the classical image setting, a large family of diffeomorphism-based methods⁴ has been developed within the field of Computational Anatomy (CA)⁵,⁶ for transforming coordinate systems at the tissue scales. The smoothness, invertibility, and non-rigid nature of diffeomorphisms are particularly relevant in medical image alignment as they reflect the mechanics/dynamics of soft tissue⁷. Consequently, methods rooted in diffeomorphic registration come particularly from multiple labs in the magnetic resonance imaging (MRI) community, with variations developed in the contexts of both atlas construction and segmentation at these tissue scales⁷⁻¹⁶. Large Deformation Diffeomorphic Metric Mapping (LDDMM)¹⁷, on which many such methods are built, specifically equips diffeomorphisms with a metric that allows them to be reduced to a low-dimensional representation that conveniently identifies shape. Consequently, it has been harnessed in the setting of image classification and the study of atrophy in diseased versus control populations¹⁷⁻¹⁹. While successful at aligning images of the same modality, most of these diffeomorphism-based methods focus exclusively on generating geometric transformations (e.g. diffeomorphisms). The challenge of crossing functional modalities has instead been addressed in a classical imaging setting by coupling such methods to variations in matching cost or additional transformations. These include matching based on analytical methods using cross-correlation¹³ or localized texture features²⁰, and methods for transforming one range space to another in crossing modalities and scales based on polynomial transformations²¹, scattering transforms²², graphical models²³, and neural networks²⁴⁻²⁶.

In the molecular setting, however, both the diversity and magnitude of data measured by spatial transcriptomics technologies²⁷ often prohibit the representation of such data as classical continuous images discretized as regular grids, and consequently, the direct use of these image alignment methods at the molecular scales. As seen in repositories generated in the BRAIN Initiative Cell Census Network (BICCN) and archived at the Brain Image Library (BIL), these datasets are already on the order of terabytes and will only continue to increase as technologies shift from mouse to human measurements. Hence, while some methods in deep learning have been applied to align single-cell datasets, modeled as regular grid images, both to atlases at the histological scale²⁸ as well as reference transcriptional atlases that are beginning to emerge²⁹, most image-based tactics are limited in their ability to represent this high-dimensional and memory-intensive data. Furthermore, many learning-based schemes in the context of classical images have relied on the use of extensive numbers (typically on the order of hundreds) of images from different subjects for extracting features useful for cross-modalities and mapping¹⁸,²⁶. In the setting of spatial transcriptomics technologies that have been developed only over the last few years, the acquisition of such large training datasets is typically prohibitive for some of these learning-based approaches.

Consequently, an independent class of methods has been developing for aligning spatially resolved transcriptional datasets at these molecular scales. With influence from image-based methods, some of these including GPSA³ and PASTE³⁰ focus exclusively on spot-based data in which gene expression is measured in a neighborhood of each spot for a regular array of spots, analogous to a voxel grid in an image. In contrast, image-based technologies such as STARmap, BARseq, and MERFISH, which take point measures of individual mRNA molecules or

cells, cover measurements irregularly sampled over space according to tissue architecture and dynamics. Furthermore, natural fluctuation in gene expression over time and space coupled to the dynamics of each spatial transcriptomics technology leads each tissue section, at the molecular/cellular (0.3–100 micron) scales, to have a varying number of such particles with no natural ordering of particles consistently apparent between sections. Consequently, landmark-based methods³¹ that assume direct permutation correspondence between particles are not applicable. Methods aimed at generalizing to allow alignment within and across these additional technologies are typically rooted in different data representations, such as graphs, in the case of SLAT³² that aims to find correspondence between cells or groups of cells in an atlas and target or CAST³³, a deep graph-based neural network (GNN) algorithm which learns a graph representation of single-cell resolution omics data and subsequently aims to align datasets via these graph representations. Additionally, while many methods have showcased success at aligning different replicates within a single technology for both human and mouse samples, cross-technology alignment has typically relied on a nonempty intersection of feature sets across these technologies²⁸,³²,³³, limiting extension of these methods to the type of cross-modality and cross-scale mapping required for integrating tissue-scale atlases and molecular scale data.

Here, we build on both of these lineages in presenting a model equally equipped at representing tissue-scale atlases and molecular-scale data, and an associated cross-modality registration mechanism that addresses both challenges of crossing scales and functional modalities by estimating simultaneously geometric and functional transformations to align each dataset to the other. Specifically, we harness the generalizability of image varifolds, which have emerged in molecular CA³⁴, for simultaneously modeling molecular and tissue-scale data with both irregular and gridlike sampling schemes in a common framework (addressed in SLAT and CAST, for instance, by a graph-based representation in the setting of molecular and cellular data³²,³³). A subproblem covered by the image-varifold theory outlined in ref. 36 is the mapping of molecular scale data to atlas coordinate systems. As specified there, we estimate minimal energy diffeomorphic transformations as in image-based LDDMM¹⁷, but where the action of diffeomorphisms on images has been adapted to the setting, here, of image varifolds in a consistent manner³⁴. Additionally, we estimate in tandem a latent distribution over the molecular functional space for each atlas partition, without the need for the large datasets often used in deep learning approaches, but by instead relying on an assumption of spatial homogeneity in distribution across each atlas partition. We use this latent distribution to transform the functional space of each atlas (i.e., its ontology) to the molecular functional space. We refer to this method of jointly estimating diffeomorphisms and latent feature distributions, as Cross Image-varifold LDDMM (xIV-LDDMM), where the cross emphasizes its ability to map across scales and modalities³⁵. We demonstrate this methodology with two common implementations of image varifolds (triangled meshes and point clouds). First, we introduce the methodology through a demonstration of mapping both gene-based and cell-based MERFISH datasets to corresponding sections of the Allen Common Coordinate Framework (CCFv3)¹. Second, we quantify accuracy in mapping specifically characterized cell types in BARseq data³⁶ to the CCFv3. Finally, we illustrate the method's generalizability to additional modalities and offer diverse examples of its use across the spans of biology in atlas construction, cross-replicate and cross-species comparison, and data segmentation.

## Results

### Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms

To accommodate the high feature dimensionality and spatial irregularity of molecular datasets, as described in the Introduction, we harness the recent work that extends the theory of diffeomorphisms and
image-based LDDMM¹⁷ to the setting of image varifolds and estimation of correspondences between them³⁴,³⁵. We first describe this framework of image varifolds to emphasize their capacity for modeling diverse types of tissue-, cellular-, and molecular-scale data. Image varifolds are geometric measures over the product space, $R^2 \times \mathcal{F}$, therefore encompassing measures over both physical and feature spaces, $R^2$ and $\mathcal{F}$, respectively. At the finest scale of capture in both classical imaging and image-based spatial transcriptomics technologies, we might model a set of pixels or detections (e.g. mRNA reads) as a discrete set of point measures (particles) with an elementary Dirac measure, centered over physical location and feature value for each detection. At coarser (e.g. cellular) scales or in spot-based technologies, however, point measurements taken at single physical locations may capture a range of individual detections. To generalize to the range of technologies, we model both tissue-scale atlases and molecular data with semi-discrete image varifolds where the physical arrangement of measurements is captured by a collection of discrete point (Dirac) measures, but where each such point measure is associated with a full distribution over feature values measured in its neighborhood. This collection of point measures is indexed by $i \in I$, with discrete measures, $\delta_{x_i}$, evaluating to 1 at locations $x_i \in R^2$, and distributions over feature values modeled as weighted probability distributions: $w_i p_i$. Weight, $w_i$, is representative of total mass (e.g. total mRNA reads or total cells) measured at location $x_i$, thus enabling an estimate of density over physical space, and probability distribution, $p_i$, captures the proportions of each feature type within that mass.

The measure of the complete collection is denoted by the sum:

$$ \mu \doteq \sum_{i \in I} w_i \delta_{x_i} \otimes p_i. \qquad (1) $$

Feature spaces considered here are finite, yielding probability distributions over 10s–100s of different elements depending on the application:

* MERFISH gene sections: $w_i$ is the total mRNA at location $x_i$, and $p_i$ is the probability distribution on gene (-700).

* BARSeq cell sections: $w_i$ is the total cells at location $x_i$ and $p_i$ is the probability distribution on cell type (-30).

* CCFv3 tissue sections: $w_i = 1$ for location $x_i$ in foreground tissue and $p_i$ is the probability distribution on ontology label (-700).

Within this established framework, we align image varifolds capturing tissue-scale atlases to those capturing molecular-scale data via the estimation of two types of correspondence: one between physical coordinates ($\varphi$) and the other between feature spaces ($\pi$). These correspondences act independently and in parallel on the physical and feature measure components of the image-varifold object. Consequently, while they are both applied in the setting of optimization to evaluate the alignment of the atlas to target (top panel, Fig. 1) they can be applied individually as relevant to specific applications including data segmentation, atlas construction, and cross-specimen comparison, as depicted in the middle panel in Fig. 1.

Regarding physical correspondence, in Computational Anatomy (CA), correspondence between tissue sections is computed using coordinate transformation between the sections by solving an optimization problem characterized by the set of possible transformations to optimize the image similarity function that specifies the alignment of the sections. These transformations are modeled as affine motions and diffeomorphisms, $\varphi$, which act to generate the space of all configurations. We model physical transformations in the setting of image varifolds similarly to diffeomorphisms, optimally solving for them with LDDMM¹⁷. However, while for classical images such as for MRI, LDDMM¹⁷ uses the action of diffeomorphisms on images $I$ as classical

functions using function composition on the right with the inverse of the diffeomorphism: $\varphi \cdot (I(x) = I \varphi^{-1}(x)$ for $x \in R^d$, the action of diffeomorphisms on image varifolds is defined, as in³⁴,³⁵:

$$ \varphi \cdot \mu \doteq \sum_{i \in I} |D\varphi|_{x_i} w_i \delta_{\varphi(x_i)} \otimes p_i, \qquad (2) $$

for the form of the image varifold defined in (1). The determinant of the Jacobian, $|D\varphi|$, capturing the local expansion/contraction of physical space, is introduced to retain the given spatial density of the original object following physical transformation. Notably, the estimated diffeomorphism can be applied in both the forward and inverse directions, taking tissue-scale atlas to molecular coordinates or vice versa, as shown in the middle panel in Fig. 1. The bottom panel depicts the determinant of the Jacobian of the diffeomorphism taking CCFv3 section to MERFISH section (B,G) with areas of local expansion in red and local contraction in blue.

To define a similarity metric between the deformed atlas and the target, we also need to carry the feature component of tissue-scale atlases to the feature space of the molecular target. For this, we associate to each feature value in the atlas ontology ($\ell \in \mathcal{L}$), a distribution (denoted $\pi_\ell$) over target feature values, capturing the assumption of spatial homogeneity we make within each atlas region. Importantly, these distributions are not normalized, enabling the generation of both a measure of target mass (e.g. mRNA or cells) over physical space (given by $w'_i$) and conditional probability distribution ($p'_i$) over target features (e.g. gene or cell type). Both molecular mass and conditional probability distribution are associated with each point measured in the atlas through the mixture distribution:

$$ w'_i p'_i = w_i \sum_{\ell \in \mathcal{L}} p_i(\ell) \pi_\ell, \qquad (3) $$

with $w_i, p_i$ denoting the initial weight and probability distribution of the point measure over the atlas feature space and with

$$ w'_i p'_i(f) = \begin{cases} \text{mass of cell type } f \text{ at location } x_i \text{ for } f \in \mathcal{F} = \{\text{cell types}\} \\ \text{mass of gene } f \text{ at location } x_i \text{ for } f \in \mathcal{F} = \{\text{genes}\}. \end{cases} $$

The bottom panel in Fig. 1c, h exhibits the estimated mRNA density ($w'_i$) for each location in the corresponding CCFv3 section while the middle column of the middle panel summarizes the estimated probability distribution over genes ($p'_i$) with the gene with the highest probability denoted for each location.

Finally, the estimation and application of both geometric transformations and latent feature distributions, $\varphi, \pi$, to the atlas image varifold take it to both the physical and feature space of the target as necessary to evaluate the similarity function in the optimization scheme (top panel, Fig. 1). In the image setting, the similarity function used is often a norm on functions, and solving the problem of minimization of the norm in the space of diffeomorphisms gives the metric theory of LDDMM for generating geometric matching between exemplar anatomies³⁷,³⁸. Here, the similarity function is a norm on image varifolds (see “Molecular scale varifold norm” section), capturing proximity in both physical and feature spaces of the target. xIV-LDDMM jointly estimates optimal $\varphi$ and $\pi$ (top panel, Fig. 1) to minimize this normed difference in the space of diffeomorphisms through either simultaneous or alternating optimization algorithms using LDDMM (see “Alternating LDDMM and quadratic program algorithm for joint optimization” section), with additional regularization imposed in both settings on the estimated $\pi$ to ensure, for instance, positive values (see “Variational problems” section for explicit variational problems). Note that throughout we highlight mapping examples using
image-based LDDMM¹⁷ to the setting of image varifolds and estimation of correspondences between them³⁴,³⁵. We first describe this framework of image varifolds to emphasize their capacity for modeling diverse types of tissue-, cellular-, and molecular-scale data. Image varifolds are geometric measures over the product space, $R^2 \times \mathcal{F}$, therefore encompassing measures over both physical and feature spaces, $R^2$ and $\mathcal{F}$, respectively. At the finest scale of capture in both classical imaging and image-based spatial transcriptomics technologies, we might model a set of pixels or detections (e.g. mRNA reads) as a discrete set of point measures (particles) with an elementary Dirac measure, centered over physical location and feature value for each detection. At coarser (e.g. cellular) scales or in spot-based technologies, however, point measurements taken at single physical locations may capture a range of individual detections. To generalize to the range of technologies, we model both tissue-scale atlases and molecular data with semi-discrete image varifolds where the physical arrangement of measurements is captured by a collection of discrete point (Dirac) measures, but where each such point measure is associated with a full distribution over feature values measured in its neighborhood. This collection of point measures is indexed by $i \in I$, with discrete measures, $\delta_{x_i}$, evaluating to 1 at locations $x_i \in R^2$, and distributions over feature values modeled as weighted probability distributions: $w_i p_i$. Weight, $w_i$, is representative of total mass (e.g. total mRNA reads or total cells) measured at location $x_i$, thus enabling an estimate of density over physical space, and probability distribution, $p_i$, captures the proportions of each feature type within that mass.

The measure of the complete collection is denoted by the sum:

$$ \mu \doteq \sum_{i \in I} w_i \delta_{x_i} \otimes p_i. \qquad (1) $$

Feature spaces considered here are finite, yielding probability distributions over 10s-100s of different elements depending on the application:

* MERFISH gene sections: $w_i$ is the total mRNA at location $x_i$, and $p_i$ is the probability distribution on gene (-700).

* BARSeq cell sections: $w_i$ is the total cells at location $x_i$ and $p_i$ is the probability distribution on cell type (-30).

* CCFv3 tissue sections: $w_i = 1$ for location $x_i$ in foreground tissue and $p_i$ is the probability distribution on ontology label (-700).

Within this established framework, we align image varifolds capturing tissue-scale atlases to those capturing molecular-scale data via the estimation of two types of correspondence: one between physical coordinates ($\varphi$) and the other between feature spaces ($\pi$). These correspondences act independently and in parallel on the physical and feature measure components of the image-varifold object. Consequently, while they are both applied in the setting of optimization to evaluate the alignment of the atlas to target (top panel, Fig. 1) they can be applied individually as relevant to specific applications including data segmentation, atlas construction, and cross-specimen comparison, as depicted in the middle panel in Fig. 1.

Regarding physical correspondence, in Computational Anatomy (CA), correspondence between tissue sections is computed using coordinate transformation between the sections by solving an optimization problem characterized by the set of possible transformations to optimize the image similarity function that specifies the alignment of the sections. These transformations are modeled as affine motions and diffeomorphisms, $\varphi$, which act to generate the space of all configurations. We model physical transformations in the setting of image varifolds similarly to diffeomorphisms, optimally solving for them with LDDMM¹⁷. However, while for classical images such as for MRI, LDDMM¹⁷ uses the action of diffeomorphisms on images $I$ as classical

functions using function composition on the right with the inverse of the diffeomorphism: $\varphi \cdot (I(x) = I \varphi^{-1}(x)$ for $x \in R^d$, the action of diffeomorphisms on image varifolds is defined, as in³⁴,³⁵:

$$ \varphi \cdot \mu \doteq \sum_{i \in I} |D\varphi|_{x_i} w_i \delta_{\varphi(x_i)} \otimes p_i, \qquad (2) $$

for the form of the image varifold defined in (1). The determinant of the Jacobian, $|D\varphi|$, capturing the local expansion/contraction of physical space, is introduced to retain the given spatial density of the original object following physical transformation. Notably, the estimated diffeomorphism can be applied in both the forward and inverse directions, taking tissue-scale atlas to molecular coordinates or vice versa, as shown in the middle panel in Fig. 1. The bottom panel depicts the determinant of the Jacobian of the diffeomorphism taking CCFv3 section to MERFISH section (B,G) with areas of local expansion in red and local contraction in blue.

To define a similarity metric between the deformed atlas and the target, we also need to carry the feature component of tissue-scale atlases to the feature space of the molecular target. For this, we associate to each feature value in the atlas ontology ($\ell \in \mathcal{L}$), a distribution (denoted $\pi_\ell$) over target feature values, capturing the assumption of spatial homogeneity we make within each atlas region. Importantly, these distributions are not normalized, enabling the generation of both a measure of target mass (e.g. mRNA or cells) over physical space (given by $w'_i$) and conditional probability distribution ($p'_i$) over target features (e.g. gene or cell type). Both molecular mass and conditional probability distribution are associated with each point measured in the atlas through the mixture distribution:

$$ w'_i p'_i = w_i \sum_{\ell \in \mathcal{L}} p_i(\ell) \pi_\ell, \qquad (3) $$

with $w_i, p_i$ denoting the initial weight and probability distribution of the point measure over the atlas feature space and with

$$ w'_i p'_i(f) = \begin{cases} \text{mass of cell type } f \text{ at location } x_i \text{ for } f \in \mathcal{F} = \{\text{cell types}\} \\ \text{mass of gene } f \text{ at location } x_i \text{ for } f \in \mathcal{F} = \{\text{genes}\}. \end{cases} $$

The bottom panel in Fig. 1c, h exhibits the estimated mRNA density ($w'_i$) for each location in the corresponding CCFv3 section while the middle column of the middle panel summarizes the estimated probability distribution over genes ($p'_i$) with the gene with the highest probability denoted for each location.

Finally, the estimation and application of both geometric transformations and latent feature distributions, $\varphi, \pi$, to the atlas image varifold take it to both the physical and feature space of the target as necessary to evaluate the similarity function in the optimization scheme (top panel, Fig. 1). In the image setting, the similarity function used is often a norm on functions, and solving the problem of minimization of the norm in the space of diffeomorphisms gives the metric theory of LDDMM for generating geometric matching between exemplar anatomies³⁷,³⁸. Here, the similarity function is a norm on image varifolds (see “Molecular scale varifold norm” section), capturing proximity in both physical and feature spaces of the target. xIV-LDDMM jointly estimates optimal $\varphi$ and $\pi$ (top panel, Fig. 1) to minimize this normed difference in the space of diffeomorphisms through either simultaneous or alternating optimization algorithms using LDDMM (see “Alternating LDDMM and quadratic program algorithm for joint optimization” section), with additional regularization imposed in both settings on the estimated $\pi$ to ensure, for instance, positive values (see “Variational problems” section for explicit variational problems). Note that throughout we highlight mapping examples using
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

both mesh-based³⁵ and particle-based (point cloud)³⁴ implementations of image varifolds, with detailed construction, similarities, and differences of each covered in “Construction of image-varifold representation for different modalities” section. Supplementary Fig. 1 shows equivalence in these two implementations and in the corresponding mappings estimated in each case for the CCFv3 and MERFISH sections shown in Fig. 1.

## Generalizability of xIV-LDDMM to settings of complete and incomplete gene-based and cell-based data

Here, we emphasize the generalizability of xIV-LDDMM for computing mappings from tissue-scale atlases to datasets spanning a range of both gene-based and cell-based features. The bottom panel of Fig. 1 specifically summarizes the estimated latent distribution over genes for each atlas region with the sum over gene counts, exhibiting the
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

both mesh-based³⁵ and particle-based (point cloud)³⁴ implementations of image varifolds, with detailed construction, similarities, and differences of each covered in “Construction of image-varifold representation for different modalities” section. Supplementary Fig. 1 shows equivalence in these two implementations and in the corresponding mappings estimated in each case for the CCFv3 and MERFISH sections shown in Fig. 1.

## Generalizability of xIV-LDDMM to settings of complete and incomplete gene-based and cell-based data

Here, we emphasize the generalizability of xIV-LDDMM for computing mappings from tissue-scale atlases to datasets spanning a range of both gene-based and cell-based features. The bottom panel of Fig. 1 specifically summarizes the estimated latent distribution over genes for each atlas region with the sum over gene counts, exhibiting the
<caption>**Fig. 1 | Methodology and results of xIV-LDDMM for transforming CCFv3 sections to MERFISH spatial transcriptomic counts of 20 selected genes.**</caption>

Top panel shows an iterative optimization scheme to estimate geometric transformation (φ) and latent feature distribution (π) by minimizing the normed difference (error) between the geometric and feature-transformed CCFv3 section to target MERFISH. Middle panel illustrates the application of estimated geometric transformation (φ) to deform the CCFv3 atlas to MERFISH coordinates (left); the application of latent feature distribution (π) to generate gene distributions on initial CCFv3 geometry (middle); and the application of inverse geometric transformation (φ⁻¹) to deform MERFISH genes to CCFv3 coordinates (right). Gene with the highest probability of expression at each location is shown as a MERFISH

feature. Bottom panel illustrates the results of mapping CCFv3 sections to corresponding MERFISH sections. **a, f** 10 µm atlas sections at Z=385 and Z=485 out of 1320 visually chosen to match MERFISH architecture (**e, j**) rendered as meshes at 100 µm. **e, j** MERFISH sections rendered as meshes at 50 µm, with mRNA density depicted as a feature. **b, g** Geometric mappings (φ) of CCFv3 sections to MERFISH coordinates with the approximate determinant of the Jacobian showing areas of contraction (blue) and expansion (red). **c, h** Estimated mRNA density per atlas region (wᵢ' in (3)), as given by π shown in CCFv3 coordinates. **d, i** Estimated mRNA density shown per atlas region following geometric deformation to MERFISH coordinates.

<figure><img src="image_6.png" /><figcaption>**Fig. 2 | Measured and predicted spatial patterns of gene expression for a subset of 20 genes measured by MERFISH.** **a–c** Relative expression on target section of MERFISH transcriptomics data for three genes (*Gfap* (*a*), *Trp53i11* (*b*), *Wipf3* (*c*)) out of a set of twenty (shown at the right), with demonstrated spatial variability according to a computed mutual information score. **e–g** Predicted expression for the same three genes (*Gfap* (*e*), *Trp53i11* (*f*), *Wipf3* (*g*)) in each region of the</figcaption></figure>

CCFv3 section Z=485 out of 1320, as part of the latent distribution over genes estimated in tandem with a geometric transformation to align the CCFv3 section to MERFISH section. **d** Gene with the highest probability in MERFISH target section. **h** Predicted gene with the highest probability in estimated latent distribution for the CCFv3 section.

estimated density of mRNA over physical space. However, as described in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, a full weighted probability distribution over feature values (e.g. gene type) is estimated via $\pi$ for each atlas location that not only contributes to achieving geometric alignment but can be visualized independently as reflective of mean distributions in target data. Figure 2a–c, for instance, shows the relative expression of three genes, (*Gfap*, *Trp53i11*, *Wipf3*, respectively) over space out of a total set of 20 in the target MERFISH section shown in Fig. 1. Figure 2e–g show the corresponding expression of each of these genes in each region of the CCFv3 section, as estimated through $\pi$, with notable similarity in probability magnitude and variation over space. Figure 2d, h showcase another comparative summary of the estimated feature distributions for the CCFv3 section versus MERFISH expression in depicting at each location the most probable gene type. The assumption of spatial homogeneity in distribution is evidenced here, particularly in large areas of the CCFv3, such as the striatum, where each gene carries a single probability of expression across the entire region and *Kirrel3* is uniformly the most probable gene type.

The previous result solved the mapping problem between CCFv3 and MERFISH sections based on the mRNA reads directly. Often, these raw mRNA reads are segmented into discrete cells as a mode of data reduction followed by downstream analyses clustering the cells into discrete cell types. We emphasize that there are various methods for solving the segmentation to cells and thereby dimension reduction as

determined by the specific imaging technology. Some of the methods are rooted in image-based segmentation schemes such as the Watershed algorithm, operating jointly on transcriptional data and immunofluorescence images such as DAPI stains³⁹, while others utilize learning-based methods⁴⁰ for accommodating often a wider diversity of cell shapes and sizes. In either case, the assignment of mRNA reads to specific cells introduces a layer of functional information at the micron scale, which can now be modeled in lieu of or in tandem with the functional information at the nanometer scale (e.g. raw mRNA reads) as the feature space of a target image varifold to which we wish to map sections of an atlas.

We demonstrate the efficacy of xIV-LDDMM for mapping tissue-scale atlases to cellular-scale data in mapping CCFv3 section Z=675 to a section of cell-segmented MERFISH transcriptional data (courtesy of the JEFworks Lab, Johns Hopkins University) (Fig. 3). The total gene set measured is ≈500 genes, with each transcript assigned to a single cell. Transcriptional profiles per cell are clustered into 33 distinct clusters using Leiden graph-based clustering⁴¹ and annotated as cell types based on known marker genes. This gives a cell-based dataset analogous to the transcript-based dataset used in the above where we now capture the spatial density of cells, and conditional probability distributions over cell types, which can be summarized via the depiction of the cell type with the highest probability at each location in space (Fig. 3a). The latent cell type distribution per each CCFv3 region, estimated in tandem with the geometric transformation taking CCFv3 section to MERFISH spatial coordinates exhibits similarity to the
<caption>**Fig. 1 | Methodology and results of xIV-LDDMM for transforming CCFv3 sections to MERFISH spatial transcriptomic counts of 20 selected genes.**</caption>

Top panel shows an iterative optimization scheme to estimate geometric transformation (φ) and latent feature distribution (π) by minimizing the normed difference (error) between the geometric and feature-transformed CCFv3 section to target MERFISH. Middle panel illustrates the application of estimated geometric transformation (φ) to deform the CCFv3 atlas to MERFISH coordinates (left); the application of latent feature distribution (π) to generate gene distributions on initial CCFv3 geometry (middle); and the application of inverse geometric transformation (φ⁻¹) to deform MERFISH genes to CCFv3 coordinates (right). Gene with the highest probability of expression at each location is shown as a MERFISH

feature. Bottom panel illustrates the results of mapping CCFv3 sections to corresponding MERFISH sections. **a, f** 10 µm atlas sections at Z=385 and Z=485 out of 1320 visually chosen to match MERFISH architecture (**e, j**) rendered as meshes at 100 µm. **e, j** MERFISH sections rendered as meshes at 50 µm, with mRNA density depicted as a feature. **b, g** Geometric mappings (φ) of CCFv3 sections to MERFISH coordinates with the approximate determinant of the Jacobian showing areas of contraction (blue) and expansion (red). **c, h** Estimated mRNA density per atlas region (wᵢ' in (3)), as given by π shown in CCFv3 coordinates. **d, i** Estimated mRNA density shown per atlas region following geometric deformation to MERFISH coordinates.

<figure><img src="image_8.png" /><figcaption>**Fig. 2 | Measured and predicted spatial patterns of gene expression for a subset of 20 genes measured by MERFISH.** **a–c** Relative expression on target section of MERFISH transcriptomics data for three genes (*Gfap* (*a*), *Trp53i11* (*b*), *Wipf3* (*c*)) out of a set of twenty (shown at the right), with demonstrated spatial variability according to a computed mutual information score. **e–g** Predicted expression for the same three genes (*Gfap* (*e*), *Trp53i11* (*f*), *Wipf3* (*g*)) in each region of the</figcaption></figure>

CCFv3 section Z=485 out of 1320, as part of the latent distribution over genes estimated in tandem with a geometric transformation to align the CCFv3 section to MERFISH section. **d** Gene with the highest probability in MERFISH target section. **h** Predicted gene with the highest probability in estimated latent distribution for the CCFv3 section.

estimated density of mRNA over physical space. However, as described in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, a full weighted probability distribution over feature values (e.g. gene type) is estimated via $\pi$ for each atlas location that not only contributes to achieving geometric alignment but can be visualized independently as reflective of mean distributions in target data. Figure 2a–c, for instance, shows the relative expression of three genes, (*Gfap*, *Trp53i11*, *Wipf3*, respectively) over space out of a total set of 20 in the target MERFISH section shown in Fig. 1. Figure 2e–g show the corresponding expression of each of these genes in each region of the CCFv3 section, as estimated through $\pi$, with notable similarity in probability magnitude and variation over space. Figure 2d, h showcase another comparative summary of the estimated feature distributions for the CCFv3 section versus MERFISH expression in depicting at each location the most probable gene type. The assumption of spatial homogeneity in distribution is evidenced here, particularly in large areas of the CCFv3, such as the striatum, where each gene carries a single probability of expression across the entire region and *Kirrel3* is uniformly the most probable gene type.

The previous result solved the mapping problem between CCFv3 and MERFISH sections based on the mRNA reads directly. Often, these raw mRNA reads are segmented into discrete cells as a mode of data reduction followed by downstream analyses clustering the cells into discrete cell types. We emphasize that there are various methods for solving the segmentation to cells and thereby dimension reduction as

determined by the specific imaging technology. Some of the methods are rooted in image-based segmentation schemes such as the Watershed algorithm, operating jointly on transcriptional data and immunofluorescence images such as DAPI stains³⁹, while others utilize learning-based methods⁴⁰ for accommodating often a wider diversity of cell shapes and sizes. In either case, the assignment of mRNA reads to specific cells introduces a layer of functional information at the micron scale, which can now be modeled in lieu of or in tandem with the functional information at the nanometer scale (e.g. raw mRNA reads) as the feature space of a target image varifold to which we wish to map sections of an atlas.

We demonstrate the efficacy of xIV-LDDMM for mapping tissue-scale atlases to cellular-scale data in mapping CCFv3 section Z=675 to a section of cell-segmented MERFISH transcriptional data (courtesy of the JEFworks Lab, Johns Hopkins University) (Fig. 3). The total gene set measured is ≈500 genes, with each transcript assigned to a single cell. Transcriptional profiles per cell are clustered into 33 distinct clusters using Leiden graph-based clustering⁴¹ and annotated as cell types based on known marker genes. This gives a cell-based dataset analogous to the transcript-based dataset used in the above where we now capture the spatial density of cells, and conditional probability distributions over cell types, which can be summarized via the depiction of the cell type with the highest probability at each location in space (Fig. 3a). The latent cell type distribution per each CCFv3 region, estimated in tandem with the geometric transformation taking CCFv3 section to MERFISH spatial coordinates exhibits similarity to the
<figure><img src="image_7.png" /><figcaption>**Fig. 3** | Mapping of CCFv3 section Z = 675 out of 1320 to single MERFISH cell-based section with either cell type or gene feature spaces. a Top row depicts MERFISH target rendered as a mesh with cell density (left) and cell type with the highest probability (right) used to summarize cell type distributions. Bottom shows predicted cell density (left) and cell type with the highest probability (right) for the latent feature distribution estimated for each CCFv3 region in the native CCFv3 coordinates.</figcaption></figure>

b, c Top row depicts the probability of expression for each gene out of a subset of 6 selected from a total measured set of ~500 as those with high spatial variance. Bottom row depicts the estimated probability of expression for each gene for the latent feature distribution estimated for each CCFv3 region in the native CCFv3 coordinates.

MERFISH target in both total cell density and cell type with highest probability predicted per CCFv3 region (Fig. 3a). Alternatively, each MERFISH section can be characterized with a feature space over genes, and a corresponding latent distribution over genes can be estimated via xIV-LDDMM in mapping the same CCFv3 section to MERFISH gene

target (Fig. 3b, c). Here, we subsample the entire gene space by selecting a subset of six genes (*Baiap2*, *Slc17a6*, *Adora2a*, *Gpr151*, *Gabbr2*, *Cckar*) as those most spatially varying according to Moran's I score⁴². We combine both cell density and gene expression by aggregating the individual mRNA transcripts into an average gene
<figure><img src="image_9.png" /><figcaption>**Fig. 3** | Mapping of CCFv3 section Z = 675 out of 1320 to single MERFISH cell-based section with either cell type or gene feature spaces. a Top row depicts MERFISH target rendered as a mesh with cell density (left) and cell type with the highest probability (right) used to summarize cell type distributions. Bottom shows predicted cell density (left) and cell type with the highest probability (right) for the latent feature distribution estimated for each CCFv3 region in the native CCFv3 coordinates.</figcaption></figure>

b, c Top row depicts the probability of expression for each gene out of a subset of 6 selected from a total measured set of ~500 as those with high spatial variance. Bottom row depicts the estimated probability of expression for each gene for the latent feature distribution estimated for each CCFv3 region in the native CCFv3 coordinates.

MERFISH target in both total cell density and cell type with highest probability predicted per CCFv3 region (Fig. 3a). Alternatively, each MERFISH section can be characterized with a feature space over genes, and a corresponding latent distribution over genes can be estimated via xIV-LDDMM in mapping the same CCFv3 section to MERFISH gene

target (Fig. 3b, c). Here, we subsample the entire gene space by selecting a subset of six genes (*Baiap2*, *Slc17a6*, *Adora2a*, *Gpr151*, *Gabbr2*, *Cckar*) as those most spatially varying according to Moran's I score⁴². We combine both cell density and gene expression by aggregating the individual mRNA transcripts into an average gene
expression feature per cell and normalizing the total mRNA per cell to 1. The target and estimated probabilities for each gene out of this subset of six are shown following the estimation of an alternative geometric transformation and latent distribution, now, over genes rather than cell types. Correspondence is seen in both the absolute magnitude as well as the relative probability of each gene between those in the target MERFISH section and those resulting from the estimated feature laws, $\pi$. Two examples include the area of high Baiap2 expression adjacent to the hippocampus (probability ≈ 0.95) and the area of high Slc17a6 expression in the rhomboid nucleus (probability ≈ 0.6), both depicted in Fig. 3b. This correspondence serves to reinforce the validity in the estimated feature laws in xIV-LDDMM.

Notably, the assumed biological correlation between cell type and pattern of gene expression implies that signals of variation across cell types at the scale of microns should also exist across gene types at the scale of nanometers. Consequently, we might expect similar spatial deformation of a tissue-scale atlas in mapping onto the same geometric target, but with conditional feature distributions defined over either gene or cell types, with partition boundaries deforming to match regions of homogeneity that would be roughly consistent across genes and cells. Supplementary Fig. 2 shows the diffeomorphism estimated in the two cases (6 genes and cell type) in Fig. 3 together with that for a third feature space comprised of a second set of 7 different genes selected for high spatial variance using mutual information. We observe the global similarity in spatial pattern and magnitude of contraction and expansion occurring but with nuanced differences in areas where these feature spaces likely emphasize different boundaries. Hence, the manifest stability in the geometric mappings jointly estimated with the feature laws over three different feature spaces supports the stability of our method in the face of different numbers and types of features but also speaks to the stability of the biological organization across tissue, cellular, and molecular scales.

Finally, we emphasize the generalizability of our mapping methodology across particularly image-based spatial transcriptomics technologies by showcasing the alignment of a cell-segmented BARseq³⁶,⁴³ partial coronal section to a corresponding full coronal CCFv3 section (Supplementary Fig. 3). Note that cells were typed according to the same procedure used for the full coronal section highlighted in “Quantitative and comparative evaluation of xIV-LDDMM” section. The costs involved with many of these emerging technologies coupled with their use in studying particular subregions or brain circuits of interest has caused the partial measurement of tissue areas to be incredibly prevalent across both image-based and spot-based technologies⁴⁴. Consequently, this generates even greater variety in the scope and shape of tissue sections measured, only further emphasizing the need for aligning such partial captures to a common scaffold, such as an atlas coordinate system, where information can be merged across the intersection of these captures. The image-varifold representation, particularly in its semi-discrete form as used here, is amenable to representing data with regions of missing or disrupted capture as an extension of the irregularity in sampling that is assumed in the general case. To estimate mappings in these cases, variation in both the forms and sizes of kernels governing the varifold norm (see “Molecular scale varifold norm” section) can control the granularity to which matches between atlas and target should be evaluated, with coarser scales more appropriate to noisier tissue sections with higher numbers of artifacts. Additionally, the varifold normed difference in the cost function, as presented in (10) in “Variational problems” section can be appended with spatially varying weights, as used in settings of mapping digital pathology to MRI²² to prioritize matching amongst certain intact regions over others. We utilize this strategy to align a partial BARseq section to a corresponding full CCFv3 section where initial geometric offset (Supplementary Fig. 3A) requires accurate estimation of scale,

rigid, and diffeomorphic transformations to position the partial capture within the scope of the correct CCFv3 hemisphere (Supplementary Fig. 3B).

## Quantitative and comparative evaluation of xIV-LDDMM

We specifically evaluated the efficacy of xIV-LDDMM in mapping CCFv3 sections to corresponding coronal sections of BARseq cell-segmented and subsequently cell-typed data. As described in ref. 37, cells were segmented in the BARseq data using Cellpose⁴⁵. Gene reads were assigned to cells, and cells containing fewer than 5 unique genes and 20 total gene counts were excluded. Cells were clustered using an iterative clustering approach based on Louvain clustering to achieve a similar resolution as a subclass in recent single-cell RNAseq studies⁴⁶.

Given the difference in atlas and target modalities mapped with xIV-LDDMM, traditional measures of accuracy, such as Dice overlap score, as are often used to evaluate image registration tools cannot be directly applied. Furthermore, landmarks common to both atlas and target modalities are not typically readily available or necessarily identifiable given the diversity in measures taken by molecular modalities³. Consequently, to evaluate the accuracy, particularly in the geometric alignment achieved with xIV-LDDMM, we matched corresponding atlas feature values (regions) with target feature values (e.g. cell types) and quantified the resulting distance individual cells of a given type from the atlas region in which we expect to find them.

We quantified accuracy by computing the set distance between cells of types specific to hippocampal regions (e.g. CA1, CA2, and CA3 pyramidal cells and DG granule cells) and particles of these regions in the CCFv3. We specifically defined this distance for each cell, indexed by $c \in C_R$, for region R, as:

$$d_c = \min_{a \in A_R} \| x_a - y_c \|_2^2, \quad (4)$$

where location $y_c$ was the given 2D coordinate of the cell in the native BARseq coordinate system, and the set $A_R$ refers to the set of indexed locations in the atlas section with feature value (region label) of R (e.g. with $p_a(R) = 1$, and $p_a(f) = 0$ for $f \neq R$). We compared these distances to those after deformation, by replacing $y_c$ with the mapped position of each cell marker to the CCFv3 section: $\phi^{-1}(y_c)$, via the inverse diffeomorphism estimated in our joint image-varifold-based method. Figure 4 depicts these measures for two separate CCFv3 sections mapped to two corresponding sections of BARseq cell data. Initial positions (Fig. 4a, c) of CCFv3 regions versus cells give distances (Fig. 4e, g) on the order of 1 mm, where notably, neighboring particles in the CCFv3 are at 10 µm, thus giving a lower bound to distance metrics we might expect for cells to be neighboring CCFv3 particles. In contrast, the positioning of CCFv3 regions versus cells following geometric transformation reflects the effects of rotations, translations, scale, and diffeomorphisms effectively estimated to bring them into correspondence (Fig. 4b, d). Kernel size governing diffeomorphism regularization and varifold norm matching were both on the order of 100 µm. Median distances for all three regions in the first section are on the order of 20–30 µm (Fig. 4f) whereas those in the second section are on the order of 5 µm for CA1 and DG and 300 µm for CA3 (Fig. 4h), which we assume is coming from the separate group of cells labeled as CA3 at the edge of CA1 (indicated by the white star) in contrast to the group overlaying the CCFv3 region, which falls within the first 10 µm bin. Note that this discrepancy likely stems in part from the approximation of each BARseq section as a strictly coronal section of the CCFv3, where the cutting plane is instead slightly offset from this coronal plane. In any case, we observe in both sections an accuracy on the order of 10–30 µm, analogous in the image setting to cells being mapped to within 1–3 pixels (at a resolution of 10 µm) to the appropriate atlas region.
expression feature per cell and normalizing the total mRNA per cell to 1. The target and estimated probabilities for each gene out of this subset of six are shown following the estimation of an alternative geometric transformation and latent distribution, now, over genes rather than cell types. Correspondence is seen in both the absolute magnitude as well as the relative probability of each gene between those in the target MERFISH section and those resulting from the estimated feature laws, $\pi$. Two examples include the area of high Baiap2 expression adjacent to the hippocampus (probability ≈ 0.95) and the area of high Slc17a6 expression in the rhomboid nucleus (probability ≈ 0.6), both depicted in Fig. 3b. This correspondence serves to reinforce the validity in the estimated feature laws in xIV-LDDMM.

Notably, the assumed biological correlation between cell type and pattern of gene expression implies that signals of variation across cell types at the scale of microns should also exist across gene types at the scale of nanometers. Consequently, we might expect similar spatial deformation of a tissue-scale atlas in mapping onto the same geometric target, but with conditional feature distributions defined over either gene or cell types, with partition boundaries deforming to match regions of homogeneity that would be roughly consistent across genes and cells. Supplementary Fig. 2 shows the diffeomorphism estimated in the two cases (6 genes and cell type) in Fig. 3 together with that for a third feature space comprised of a second set of 7 different genes selected for high spatial variance using mutual information. We observe the global similarity in spatial pattern and magnitude of contraction and expansion occurring but with nuanced differences in areas where these feature spaces likely emphasize different boundaries. Hence, the manifest stability in the geometric mappings jointly estimated with the feature laws over three different feature spaces supports the stability of our method in the face of different numbers and types of features but also speaks to the stability of the biological organization across tissue, cellular, and molecular scales.

Finally, we emphasize the generalizability of our mapping methodology across particularly image-based spatial transcriptomics technologies by showcasing the alignment of a cell-segmented BARseq³⁶,⁴³ partial coronal section to a corresponding full coronal CCFv3 section (Supplementary Fig. 3). Note that cells were typed according to the same procedure used for the full coronal section highlighted in “Quantitative and comparative evaluation of xIV-LDDMM” section. The costs involved with many of these emerging technologies coupled with their use in studying particular subregions or brain circuits of interest has caused the partial measurement of tissue areas to be incredibly prevalent across both image-based and spot-based technologies⁴⁴. Consequently, this generates even greater variety in the scope and shape of tissue sections measured, only further emphasizing the need for aligning such partial captures to a common scaffold, such as an atlas coordinate system, where information can be merged across the intersection of these captures. The image-varifold representation, particularly in its semi-discrete form as used here, is amenable to representing data with regions of missing or disrupted capture as an extension of the irregularity in sampling that is assumed in the general case. To estimate mappings in these cases, variation in both the forms and sizes of kernels governing the varifold norm (see “Molecular scale varifold norm” section) can control the granularity to which matches between atlas and target should be evaluated, with coarser scales more appropriate to noisier tissue sections with higher numbers of artifacts. Additionally, the varifold normed difference in the cost function, as presented in (10) in “Variational problems” section can be appended with spatially varying weights, as used in settings of mapping digital pathology to MRI²² to prioritize matching amongst certain intact regions over others. We utilize this strategy to align a partial BARseq section to a corresponding full CCFv3 section where initial geometric offset (Supplementary Fig. 3A) requires accurate estimation of scale,

rigid, and diffeomorphic transformations to position the partial capture within the scope of the correct CCFv3 hemisphere (Supplementary Fig. 3B).

## Quantitative and comparative evaluation of xIV-LDDMM

We specifically evaluated the efficacy of xIV-LDDMM in mapping CCFv3 sections to corresponding coronal sections of BARseq cell-segmented and subsequently cell-typed data. As described in ref. 37, cells were segmented in the BARseq data using Cellpose⁴⁵. Gene reads were assigned to cells, and cells containing fewer than 5 unique genes and 20 total gene counts were excluded. Cells were clustered using an iterative clustering approach based on Louvain clustering to achieve a similar resolution as a subclass in recent single-cell RNAseq studies⁴⁶.

Given the difference in atlas and target modalities mapped with xIV-LDDMM, traditional measures of accuracy, such as Dice overlap score, as are often used to evaluate image registration tools cannot be directly applied. Furthermore, landmarks common to both atlas and target modalities are not typically readily available or necessarily identifiable given the diversity in measures taken by molecular modalities³. Consequently, to evaluate the accuracy, particularly in the geometric alignment achieved with xIV-LDDMM, we matched corresponding atlas feature values (regions) with target feature values (e.g. cell types) and quantified the resulting distance individual cells of a given type from the atlas region in which we expect to find them.

We quantified accuracy by computing the set distance between cells of types specific to hippocampal regions (e.g. CA1, CA2, and CA3 pyramidal cells and DG granule cells) and particles of these regions in the CCFv3. We specifically defined this distance for each cell, indexed by $c \in C_R$, for region R, as:

$$d_c = \min_{a \in A_R} \| x_a - y_c \|_2^2, \quad (4)$$

where location $y_c$ was the given 2D coordinate of the cell in the native BARseq coordinate system, and the set $A_R$ refers to the set of indexed locations in the atlas section with feature value (region label) of R (e.g. with $p_a(R) = 1$, and $p_a(f) = 0$ for $f \neq R$). We compared these distances to those after deformation, by replacing $y_c$ with the mapped position of each cell marker to the CCFv3 section: $\phi^{-1}(y_c)$, via the inverse diffeomorphism estimated in our joint image-varifold-based method. Figure 4 depicts these measures for two separate CCFv3 sections mapped to two corresponding sections of BARseq cell data. Initial positions (Fig. 4a, c) of CCFv3 regions versus cells give distances (Fig. 4e, g) on the order of 1 mm, where notably, neighboring particles in the CCFv3 are at 10 µm, thus giving a lower bound to distance metrics we might expect for cells to be neighboring CCFv3 particles. In contrast, the positioning of CCFv3 regions versus cells following geometric transformation reflects the effects of rotations, translations, scale, and diffeomorphisms effectively estimated to bring them into correspondence (Fig. 4b, d). Kernel size governing diffeomorphism regularization and varifold norm matching were both on the order of 100 µm. Median distances for all three regions in the first section are on the order of 20–30 µm (Fig. 4f) whereas those in the second section are on the order of 5 µm for CA1 and DG and 300 µm for CA3 (Fig. 4h), which we assume is coming from the separate group of cells labeled as CA3 at the edge of CA1 (indicated by the white star) in contrast to the group overlaying the CCFv3 region, which falls within the first 10 µm bin. Note that this discrepancy likely stems in part from the approximation of each BARseq section as a strictly coronal section of the CCFv3, where the cutting plane is instead slightly offset from this coronal plane. In any case, we observe in both sections an accuracy on the order of 10–30 µm, analogous in the image setting to cells being mapped to within 1–3 pixels (at a resolution of 10 µm) to the appropriate atlas region.
<figure><img src="image_8.png" /><figcaption>**Fig. 4** | Quantification of accuracy in two coronal slices of CCFv3 mapped to two corresponding slices of BARseq cell data classified into 52 cell types. a, b Initial positions of cell markers (circles) from BARseq sections in three hippocampal regions (CA1, CA3, Dentate Gyrus (DG)) and initial positions of equivalent CCFv3 regions in corresponding coronal sections (Z = 877 and Z = 837 out of 1320). c, d Set distance of cell markers to CCFv3 region in initial positions. e, g Set distance of cell markers to CCFv3 region in initial positions. f, h Positions of cell markers (circles) in BARseq sections mapped to CCFv3 coordinates with xIV-LDDMM and positions of equivalent CCFv3 regions in corresponding coronal sections. f, h Set distance of cell markers to CCFv3 region in aligned positions following estimation of diffeomorphism with xIV-LDDMM. Percents of cell markers clipped to 25% in all histograms, with 70%, 50%, and 60% cell markers within 10 µm in CA1, CA3, and DG, respectively in (f) and 95%, 35%, 60% cell markers within 10 µm in CA1, CA3, and DG, respectively in (h).</figcaption></figure>

We also compared our joint image-varifold-based approach in xIV-LDDMM to a classical image-based approach matching, estimating only a geometric transformation and matching based on foreground-background (Fig. 5). We mapped one of the same CCFv3 sections (Z= 837) to an image rendering of the BARseq cell section, with grayscale intensity capturing density of cells over space. Both CCFv3 and BARseq images were discretized at 50 µm resolution to estimate the mapping, with the latter computed with Gaussian smoothing, and image-based LDDMM, as described in ref. 17 was used to estimate optimal diffeomorphic transformation of the atlas to the target. To compare these methods more directly, a 10 µm BARseq image was

transformed with the estimated diffeomorphism to the CCFv3 space (Fig. 5a) as closer to the particle-based representation at full resolution used for estimation of the diffeomorphism in the image-varifold based mapping (Fig. 5c). At this resolution, ~1.5 million pixels are needed to model the BARseq data as an image whereas only ~90k particles are needed in the image-varifold representation. In addition to the computational and memory expense of image handling, we compromise our ability to delineate one cell from another, as evidenced in the blur in the image rendering (Fig. 5a) versus image-varifold representation (Fig. 5c), with 75% of cells having a nearest neighbor less than 23 µm away and 25% having a nearest neighbor less than 3.5 µm away in the
<figure><img src="image_10.png" /><figcaption>**Fig. 4** | Quantification of accuracy in two coronal slices of CCFv3 mapped to two corresponding slices of BARseq cell data classified into 52 cell types. a, b Initial positions of cell markers (circles) from BARseq sections in three hippocampal regions (CA1, CA3, Dentate Gyrus (DG)) and initial positions of equivalent CCFv3 regions in corresponding coronal sections (Z = 877 and Z = 837 out of 1320). c, d Set distance of cell markers to CCFv3 region in initial positions. e, g Set distance of cell markers to CCFv3 region in initial positions. c, d Positions of cell markers (circles) in BARseq sections mapped to CCFv3 coordinates with xIV-</figcaption></figure>

LDDMM and positions of equivalent CCFv3 regions in corresponding coronal sections. f, h Set distance of cell markers to CCFv3 region in aligned positions following estimation of diffeomorphism with xIV-LDDMM. Percents of cell markers clipped to 25% in all histograms, with 70%, 50%, and 60% cell markers within 10 µm in CA1, CA3, and DG, respectively in (f) and 95%, 35%, 60% cell markers within 10 µm in CA1, CA3, and DG, respectively in (h).

We also compared our joint image-varifold-based approach in xIV-LDDMM to a classical image-based approach matching, estimating only a geometric transformation and matching based on foreground-background (Fig. 5). We mapped one of the same CCFv3 sections (Z = 837) to an image rendering of the BARseq cell section, with grayscale intensity capturing density of cells over space. Both CCFv3 and BARseq images were discretized at 50 µm resolution to estimate the mapping, with the latter computed with Gaussian smoothing, and image-based LDDMM, as described in ref. 17 was used to estimate optimal diffeomorphic transformation of the atlas to the target. To compare these methods more directly, a 10 µm BARseq image was

transformed with the estimated diffeomorphism to the CCFv3 space (Fig. 5a) as closer to the particle-based representation at full resolution used for estimation of the diffeomorphism in the image-varifold based mapping (Fig. 5c). At this resolution, ~1.5 million pixels are needed to model the BARseq data as an image whereas only ~90k particles are needed in the image-varifold representation. In addition to the computational and memory expense of image handling, we compromise our ability to delineate one cell from another, as evidenced in the blur in the image rendering (Fig. 5a) versus image-varifold representation (Fig. 5c), with 75% of cells having a nearest neighbor less than 23 µm away and 25% having a nearest neighbor less than 3.5 µm away in the
<figure><img src="image_9.png" /><figcaption>**Fig. 5** | Comparison of mapping CCFv3 section Z = 837 out of 1320 to corresponding BARseq cell data using image-based LDDMM¹⁷ versus image-varifold-based mapping method. **a** Estimated diffeomorphic transformation applied to 10 µm smoothed cell density image to bring into overall alignment with CCFv3 section in CCFv3 coordinates. Source and target images rendered at 50 µm for image-based LDDMM, with matching based on foreground/background in CCFv3 section to smoothed cell density in BARseq. **b** Estimated diffeomorphic transformation applied to 5 µm cell type image highlighting noticeable misalignment in (a) of high cell density areas to layer 1 and white matter structures (corpus callosum and dorsal hippocampal commissure) in CCFv3 section with image-based LDDMM. Layer-by-layer mismatch (-150 µm) observed between characteristic BARseq prescribed layer-specific cell type and CCFv3 delineations. **c** Estimated diffeomorphic transformation with xIV-LDDMM applied to BARseq section to bring into overall alignment with CCFv3 section in CCFv3 coordinates without misalignment in white matter structures. **d** Layer-by-layer matching in CCFv3 partition of cortex to cell types in BARseq following image-varifold based mapping.</figcaption></figure>

BARseq capture. Nevertheless, with a kernel of bandwidth 10 µm used to generate the 10 µm smoothed image and most regions at least on the order of 100 µm in width at their narrowest (e.g. corpus callosum), the majority of cell mass stayed primarily within the initial region it was found.

Regarding alignment, an image-based approach can achieve global alignment of the template to the target (Fig. 5a), with initial differences in scale and shape (Fig. 4b) resolved with the estimated transformation. However, we see clear differences in alignment on the order of 50–150 µm by examination of the overlap in the different layers of the cortex and adjacent corpus callosum (Fig. 5b, d). Importantly, as described in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, in the joint image-varifold-based approach to estimating geometric transformations and latent distributions, we model each atlas region as homogeneous and stationary with respect to space. This gives an optimal alignment between the atlas and target that maximizes similarity in distribution over features across each site in a single atlas region while minimizing the energy of the geometric deformation (diffeomorphism). We would consequently expect this to skew emphasis away from the foreground-background boundaries that typically govern image alignment and instead highlight the underlying assumptions in the architecture of the cartoon atlas, whose boundaries were initially constructed so as to maximize the homogeneity of the region. This becomes clear in comparing the alignment of BARseq to CCFv3 section in areas of low versus high cell density where these areas of low density correctly align with layer 1 and corpus callosum (CC) as a result of the image-varifold based approach (Fig. 5d) but not in the image-based approach (Fig. 5b). Notably, the corpus callosum and dorsal hippocampal commissure appear with equivalently

low cell density on both right and left hemispheres in the image-varifold-based approach (yellow arrow) whereas we see these areas partially covered by CA1 cells, particularly in the right hemisphere in the image-based approach. Layer 1 (white arrow) is covered by high cell density around the entire circumference of the coronal section in the image-based approach, whereas it is visible, with relatively few cells mapping to it, as expected, in the image-varifold-based approach. Supplementary Fig. 4 shows the relative cell density in each of the cortical layers and white matter structures within the neighborhood of the primary visual cortex, with xIV-LDDMM showing more accurate alignment with higher levels (4–6×) in layers 2–6 compared with layer 1 and white matter. Image-based LDDMM, in contrast, yields similar levels of cell density across all of these structures, with the band of high cell density in BARseq layers 2–6 covering the entirety of the cortex, thus generating overestimates (4–6×) of cell density in layer 1. Cell density within CA1 also covers the areas of the corpus callosum, leading to overestimates of density there. Finally, as the image-varifold-based approach jointly considers both variations in the total density of cells and relative distribution of cell type, we see within each layer of cortex delineated in the CCFv3, a majority of corresponding cell types indicating correct layer-by-layer alignment (Fig. 5d). Classification of cells by nearest atlas region yields 88% of layer 2/3 cells (total 443), 86% of layer 4/5 cells (total 436), 87% of layer 5 cells (total 278), and 70% of layer 6 cells (total 608) classified correctly according to cortical layer, within a 1.5 mm² section through the primary visual cortex shown in Fig. 5. In contrast, for a similar section in the result of image-based mapping (Fig. 5b), 38% of layer 2/3 pixels (total 308), 36% of layer 4 pixels (total 269), 64% of layer 5 pixels (total 225), and 78% of layer 6 pixels (total 326) are classified correctly according to nearest cortical layer.
<figure><img src="image_11.png" /><figcaption>**Fig. 5** | Comparison of mapping CCFv3 section Z = 837 out of 1320 to corresponding BARseq cell data using image-based LDDMM¹⁷ versus image-varifold-based mapping method. **a** Estimated diffeomorphic transformation applied to 10 µm smoothed cell density image to bring into overall alignment with CCFv3 section in CCFv3 coordinates. Source and target images rendered at 50 µm for image-based LDDMM, with matching based on foreground/background in CCFv3 section to smoothed cell density in BARseq. **b** Estimated diffeomorphic transformation applied to 5 µm cell type image highlighting noticeable misalignment in (a) of high cell density areas to layer 1 and white matter structures (corpus callosum and dorsal hippocampal commissure) in CCFv3 section with image-based LDDMM. Layer-by-layer mismatch (-150 µm) observed between characteristic BARseq prescribed layer-specific cell type and CCFv3 delineations. **c** Estimated diffeomorphic transformation with xIV-LDDMM applied to BARseq section to bring into overall alignment with CCFv3 section in CCFv3 coordinates without misalignment in white matter structures. **d** Layer-by-layer matching in CCFv3 partition of cortex to cell types in BARseq following image-varifold based mapping.</figcaption></figure>

BARseq capture. Nevertheless, with a kernel of bandwidth 10 µm used to generate the 10 µm smoothed image and most regions at least on the order of 100 µm in width at their narrowest (e.g. corpus callosum), the majority of cell mass stayed primarily within the initial region it was found.

Regarding alignment, an image-based approach can achieve global alignment of the template to the target (Fig. 5a), with initial differences in scale and shape (Fig. 4b) resolved with the estimated transformation. However, we see clear differences in alignment on the order of 50–150 µm by examination of the overlap in the different layers of the cortex and adjacent corpus callosum (Fig. 5b, d). Importantly, as described in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, in the joint image-varifold-based approach to estimating geometric transformations and latent distributions, we model each atlas region as homogeneous and stationary with respect to space. This gives an optimal alignment between the atlas and target that maximizes similarity in distribution over features across each site in a single atlas region while minimizing the energy of the geometric deformation (diffeomorphism). We would consequently expect this to skew emphasis away from the foreground-background boundaries that typically govern image alignment and instead highlight the underlying assumptions in the architecture of the cartoon atlas, whose boundaries were initially constructed so as to maximize the homogeneity of the region. This becomes clear in comparing the alignment of BARseq to CCFv3 section in areas of low versus high cell density where these areas of low density correctly align with layer 1 and corpus callosum (CC) as a result of the image-varifold based approach (Fig. 5d) but not in the image-based approach (Fig. 5b). Notably, the corpus callosum and dorsal hippocampal commissure appear with equivalently

low cell density on both right and left hemispheres in the image-varifold-based approach (yellow arrow) whereas we see these areas partially covered by CA1 cells, particularly in the right hemisphere in the image-based approach. Layer 1 (white arrow) is covered by high cell density around the entire circumference of the coronal section in the image-based approach, whereas it is visible, with relatively few cells mapping to it, as expected, in the image-varifold-based approach. Supplementary Fig. 4 shows the relative cell density in each of the cortical layers and white matter structures within the neighborhood of the primary visual cortex, with xIV-LDDMM showing more accurate alignment with higher levels (4–6×) in layers 2–6 compared with layer 1 and white matter. Image-based LDDMM, in contrast, yields similar levels of cell density across all of these structures, with the band of high cell density in BARseq layers 2–6 covering the entirety of the cortex, thus generating overestimates (4–6×) of cell density in layer 1. Cell density within CA1 also covers the areas of the corpus callosum, leading to overestimates of density there. Finally, as the image-varifold-based approach jointly considers both variations in the total density of cells and relative distribution of cell type, we see within each layer of cortex delineated in the CCFv3, a majority of corresponding cell types indicating correct layer-by-layer alignment (Fig. 5d). Classification of cells by nearest atlas region yields 88% of layer 2/3 cells (total 443), 86% of layer 4/5 cells (total 436), 87% of layer 5 cells (total 278), and 70% of layer 6 cells (total 608) classified correctly according to cortical layer, within a 1.5 mm² section through the primary visual cortex shown in Fig. 5. In contrast, for a similar section in the result of image-based mapping (Fig. 5b), 38% of layer 2/3 pixels (total 308), 36% of layer 4 pixels (total 269), 64% of layer 5 pixels (total 225), and 78% of layer 6 pixels (total 326) are classified correctly according to nearest cortical layer.
Cross-replicate comparison and atlas construction

As described in the Introduction, two key applications of xIV-LDDMM are in cross-sample comparison and atlas construction. We illustrate both of these applications in considering cell-typed MERFISH sections taken from three separate mice at approximately the same coronal level (Fig. 6a-c). CCFv3 section Z=890 out of 1320 was mapped separately to each of these targets with both geometric transformations and latent feature distributions estimated in each case. As in the

classical imaging setting with LDDMM, comparison of the estimated
diffeomorphisms taking the CCFv3 section to each respective target
section offers one metric of similarity between both atlas and target
and across targets³⁷,³⁸ (Supplementary Fig. 5A-C). Here, the
CCFv3 section contracts to a similar extent (-0.5) in all three mice in
areas of the midbrain around the periacqueductal gray matter. In
contrast, levels of contraction/expansion across the cortical layers in
the primary visual area range from 0.7 to 1.4 across the three mice,

<figure><img src="image_10.png" /><figcaption>**Fig. 6** | Cross-replicate comparison and atlas construction from mapping CCFv3 section Z = 890 out of 1320 to MERFISH cell-typed coronal sections from three mice at approximately the same location. **a–c** Position and cell type for all cells measured in each mouse section. **d–f** CCFv3 section geometrically transformed to each target space with the most probable cell type in estimated latent distribution shown for each atlas region. **g–i** Initial CCFv3 section with summative statistics of empirical probability distribution over cell types estimated from the inverse transformation of all three replicates to CCFv3 coordinates. **g** Variance in estimated empirical probability distributions across replicates. White dotted arrow</figcaption></figure>

highlights area around medial geniculate nucleus with variance in cell types
between replicates shown in (a–c). **h** Most likely cell type in empirical distribution
estimated from all three replicates. **i** Spatial variance in empirical probability dis-
tribution per each atlas region with white arrow highlighting region of outer cortex
varying in ependymal cell and excitatory pyramidal neuron distributions across the
space of the region. **j** Spatial variance in empirical cell type distribution following
rigid+scale (left) versus diffeomorphic (right) transformations in each of the three
mice. **k** Reduction in spatial variance in two adjacent areas of cortex (red/blue box
in **j**) following rigid+scale or diffeomorphic transformation in mice.
Cross-replicate comparison and atlas construction

As described in the Introduction, two key applications of xIV-LDDMM are in cross-sample comparison and atlas construction. We illustrate both of these applications in considering cell-typed MERFISH sections taken from three separate mice at approximately the same coronal level (Fig. 6a-c). CCFv3 section Z=890 out of 1320 was mapped separately to each of these targets with both geometric transformations and latent feature distributions estimated in each case. As in the

classical imaging setting with LDDMM, comparison of the estimated
diffeomorphisms taking the CCFv3 section to each respective target
section offers one metric of similarity between both atlas and target
and across targets³⁷,³⁸ (Supplementary Fig. 5A-C). Here, the
CCFv3 section contracts to a similar extent (-0.5) in all three mice in
areas of the midbrain around the periacqueductal gray matter. In
contrast, levels of contraction/expansion across the cortical layers in
the primary visual area range from 0.7 to 1.4 across the three mice,

<figure><img src="image_12.png" /><figcaption>**Fig. 6** | Cross-replicate comparison and atlas construction from mapping CCFv3 section Z = 890 out of 1320 to MERFISH cell-typed coronal sections from three mice at approximately the same location. **a–c** Position and cell type for all cells measured in each mouse section. **d–f** CCFv3 section geometrically transformed to each target space with the most probable cell type in estimated latent distribution shown for each atlas region. **g–i** Initial CCFv3 section with summative statistics of empirical probability distribution over cell types estimated from the inverse transformation of all three replicates to CCFv3 coordinates. **g** Variance in estimated empirical probability distributions across replicates. White dotted arrow</figcaption></figure>

highlights area around medial geniculate nucleus with variance in cell types
between replicates shown in (a–c). **h** Most likely cell type in empirical distribution
estimated from all three replicates. **i** Spatial variance in empirical probability dis-
tribution per each atlas region with white arrow highlighting region of outer cortex
varying in ependymal cell and excitatory pyramidal neuron distributions across the
space of the region. **j** Spatial variance in empirical cell type distribution following
rigid+scale (left) versus diffeomorphic (right) transformations in each of the three
mice. **k** Reduction in spatial variance in two adjacent areas of cortex (red/blue box
in **j**) following rigid+scale or diffeomorphic transformation in mice.
indicating differences in geometric shape. In xIV-LDDMM, the jointly estimated latent feature distributions offer a second metric of similarity (Supplementary Fig. 5D-G). Differences in the predicted cell type with highest probability occur, for instance, in the perirhinal and ectorhinal areas with excitatory granule cells, excitatory pyramidal neurons, or cortical excitatory neurons predominant in the different mice (Fig. 6d-f). The sample variance for each atlas region across the three estimated cell type probability distributions, modeled as vectors, $p \in \mathbb{R}^C$, for $C=33$ cell types and $|p|=1$, measures difference not just in the single cell type with highest probability, but amongst the entire estimated distributions over cell type for each mouse (see "Empirical distribution estimation and cross-replicate statistical comparison" section). Sample variance (Supplementary Fig. 5F, G) is computed as $\frac{1}{N-1} \sum_{i=1}^{N} \|p_i - \bar{p}\|_2^2$ in $\mathbb{R}^C$, with $N=3$ and $\bar{p}$, the sample mean (Supplementary Fig. 5D, E). Notably, the area of the medial geniculate nucleus, while exhibiting excitatory neurons as the predominant cell type across mice, exhibits high variance in total distribution. This region is magnified in Supplementary Fig. 6A, illustrating this variance in relative distribution of excitatory neurons, astrocytes, endothelial cells, and ependymal cells.

Regarding atlas construction, the sample mean cell type distribution across those latent distributions estimated jointly with geometric transformations (Supplementary Fig. 5D, E) gives one potential construction of a cell type atlas over space. A finer-grained atlas can alternatively be achieved in this setting by pulling back each target MERFISH section onto the same CCFv3 section via the inverse estimated diffeomorphism (Fig. 1). This eliminates some of the variance in geometry between targets (Supplementary Fig. 5A-C), facilitating more direct comparison of cell type distributions. Particles (analogous to foreground pixels) in the CCFv3 section serve as a scaffold for resampling mapped cells to with either a nearest neighbor assignment or dispersion of cell mass according to kernel choice (see "Empirical distribution estimation and cross-replicate statistical comparison" section). Atlases can be constructed from the resampling of each target individually, with variance across the replicates' empirical cell type probability distributions (Fig. 6g) similar to those estimated jointly in xIV-LDDMM (Supplementary Fig. 5D-G), but also together for generation of a population average (Fig. 6h). Importantly, each particle, here, is treated independently without assuming homogeneity in the cell type distribution across particles belonging to the same atlas region. This is in contrast to the distributions estimated jointly with geometric transformations (Fig. 6d-f), and consequently results in more locally varying distributions. Overall, ventral areas exhibit higher spatial variance than dorsal areas (Fig. 6i), exhibiting that in such regions, an assumption of homogeneity within these regions may not always be appropriate. The white arrow, in particular, points to an area of cortex where this spatial variance is amongst the highest, with the outermost portion of cortex predominantly composed of ependymal cells versus the innermost portion predominantly excitatory pyramidal cells (see Supplementary Fig. 6B).

The spatial variance in cell type distribution within each atlas region and the variance in distribution between replicates of a given population both influence the statistical power needed to detect biologically relevant differences in these distributions between populations. Classically, they are accounted for in mixed effects models⁴⁷ aimed at detecting differences in group features. For instance, we have looked at differences in atrophy rate of medial temporal lobe structures between control and diseased cohorts in Alzheimer's disease¹⁹,⁴⁸. Here, we observe greater spatial variance across ventral regions than dorsal regions (Fig. 6i), with particularly minimal variance in layers of cortex. Importantly, the use of non-rigid geometric transformations (diffeomorphisms) to pull back targets into CCFv3 coordinates reduces this spatial variance in distribution across cortical layers (Fig. 6j, k) compared with rigid and scaling transformations only in each of the three mice. Consequently, in

settings of comparing groups of replicates under different experimental conditions, we would expect this reduction in variance to facilitate detection of significant differences in cell type distribution per CCFv3 region between them.

Here, the small sample size ($n=3$) of replicates, which are all produced under control conditions, precludes any definitive statistical statement about variations in cell type distributions observed between them and within CCFv3 regions. Nevertheless, summation across all CCFv3 regions of the variance per individual cell type probability between replicates and within each region (Fig. 7) elucidates the relative contribution of each cell (sub)type probability to the variances observed in Fig. 6g, i. We observe, for instance, with both types of variance, a large contribution from the variance in ependymal cell type probabilities (Fig. 7a, b). Specifically between replicates, excitatory neurons exhibit the second largest variance in cell type probabilities, as evidenced in the area of the dentate gyrus (Fig. 7c-e). With regard to spatial variance, we see specific astrocyte subtype probabilities with large variance, exhibited particularly in the areas of CA1 (yellow arrow) and the pons (orange arrow) both medially to laterally within the regions and when considering left versus right hemispheres (Fig. 7f, g). In contrast, other astrocyte subtype probabilities (Fig. 7h) are seen to be consistent over space, in line with our assumption of homogeneity within each CCFv3 region.

## Extension to additional tissue modalities for within and across species comparison

Constructed across institutions with varying combinations of molecular, chemical, genetic, and electrophysiological signals, multiple atlases per species now exist with different levels of granularity and intended applications¹,²,⁴⁹⁻⁵². Given this plethora of atlases, each of which might define a different partitioning scheme over the same area of tissue, questions of comparison and relevance of each atlas to emerging molecular and cellular signatures naturally arise⁵³. While some atlases have been defined in the same coordinate framework—often achieved through existing methods of image registration or manual alignment⁵³—many exist in different coordinate frameworks. Together with mismatches in the number, type, and positioning of partitions, this poses a challenge not only to the evaluation of each atlas ontology's fit to a molecular target but also to the ready comparison of atlas to atlas and the establishment of a clear metric of similarity between them. Hence, a second application of our method rests in its use to map not just atlas to the molecular dataset, but one atlas to another, with the geometric and functional correspondences $(\varphi, \pi)$ yielded by our method serving as an anchor for cross-examination of existing ontologies both within and across species.

To map tissue-scale atlas to tissue-scale atlas, we model both objects as image varifolds, as described in "Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms" section, but where both template and target are constrained to be of constant density ($w_i=1$ for all $x_i$ in foreground tissue). Both geometric transformations and latent feature distributions, $(\varphi, \pi)$ are estimated, with the latter giving distributions over the target atlas features (labels) for each template atlas label, following the assumption of spatial homogeneity with each template atlas label presumed to map onto a single distribution of target atlas labels.

The joint estimation of geometric transformation, $\varphi$ and conditional feature laws, $(\pi_t)_{t \in \mathcal{L}}$ in xIV-LDDMM offers two modes of quantitative comparison between atlas ontologies. First, as in the classical image setting of LDDMM, the determinant of the Jacobian, $|D\varphi|$, of the estimated diffeomorphism, can be used as a metric of how similar the atlas ontologies are, reflective of how much boundaries of partitions move to maximize the overlap between homogeneous regions. However, unlike in classical image settings, the estimation of the additional family of feature laws here, $\pi$, affords a second metric of similarity with
indicating differences in geometric shape. In xIV-LDDMM, the jointly estimated latent feature distributions offer a second metric of similarity (Supplementary Fig. 5D-G). Differences in the predicted cell type with highest probability occur, for instance, in the perirhinal and ectorhinal areas with excitatory granule cells, excitatory pyramidal neurons, or cortical excitatory neurons predominant in the different mice (Fig. 6d-f). The sample variance for each atlas region across the three estimated cell type probability distributions, modeled as vectors, $p \in \mathbb{R}^C$, for $C=33$ cell types and $|p|=1$, measures difference not just in the single cell type with highest probability, but amongst the entire estimated distributions over cell type for each mouse (see "Empirical distribution estimation and cross-replicate statistical comparison" section). Sample variance (Supplementary Fig. 5F, G) is computed as $\frac{1}{N-1} \sum_{i=1}^{N} \| p_i - \bar{p} \|_2^2$ in $\mathbb{R}^C$, with $N=3$ and $\bar{p}$, the sample mean (Supplementary Fig. 5D, E). Notably, the area of the medial geniculate nucleus, while exhibiting excitatory neurons as the predominant cell type across mice, exhibits high variance in total distribution. This region is magnified in Supplementary Fig. 6A, illustrating this variance in relative distribution of excitatory neurons, astrocytes, endothelial cells, and ependymal cells.

Regarding atlas construction, the sample mean cell type distribution across those latent distributions estimated jointly with geometric transformations (Supplementary Fig. 5D, E) gives one potential construction of a cell type atlas over space. A finer-grained atlas can alternatively be achieved in this setting by pulling back each target MERFISH section onto the same CCFv3 section via the inverse estimated diffeomorphism (Fig. 1). This eliminates some of the variance in geometry between targets (Supplementary Fig. 5A-C), facilitating more direct comparison of cell type distributions. Particles (analogous to foreground pixels) in the CCFv3 section serve as a scaffold for resampling mapped cells to with either a nearest neighbor assignment or dispersion of cell mass according to kernel choice (see "Empirical distribution estimation and cross-replicate statistical comparison" section). Atlases can be constructed from the resampling of each target individually, with variance across the replicates' empirical cell type probability distributions (Fig. 6g) similar to those estimated jointly in xIV-LDDMM (Supplementary Fig. 5D-G), but also together for generation of a population average (Fig. 6h). Importantly, each particle, here, is treated independently without assuming homogeneity in the cell type distribution across particles belonging to the same atlas region. This is in contrast to the distributions estimated jointly with geometric transformations (Fig. 6d-f), and consequently results in more locally varying distributions. Overall, ventral areas exhibit higher spatial variance than dorsal areas (Fig. 6i), exhibiting that in such regions, an assumption of homogeneity within these regions may not always be appropriate. The white arrow, in particular, points to an area of cortex where this spatial variance is amongst the highest, with the outermost portion of cortex predominantly composed of ependymal cells versus the innermost portion predominantly excitatory pyramidal cells (see Supplementary Fig. 6B).

The spatial variance in cell type distribution within each atlas region and the variance in distribution between replicates of a given population both influence the statistical power needed to detect biologically relevant differences in these distributions between populations. Classically, they are accounted for in mixed effects models⁴⁷ aimed at detecting differences in group features. For instance, we have looked at differences in atrophy rate of medial temporal lobe structures between control and diseased cohorts in Alzheimer's disease¹⁹,⁴⁸. Here, we observe greater spatial variance across ventral regions than dorsal regions (Fig. 6i), with particularly minimal variance in layers of cortex. Importantly, the use of non-rigid geometric transformations (diffeomorphisms) to pull back targets into CCFv3 coordinates reduces this spatial variance in distribution across cortical layers (Fig. 6j, k) compared with rigid and scaling transformations only in each of the three mice. Consequently, in

settings of comparing groups of replicates under different experimental conditions, we would expect this reduction in variance to facilitate detection of significant differences in cell type distribution per CCFv3 region between them.

Here, the small sample size ($n=3$) of replicates, which are all produced under control conditions, precludes any definitive statistical statement about variations in cell type distributions observed between them and within CCFv3 regions. Nevertheless, summation across all CCFv3 regions of the variance per individual cell type probability between replicates and within each region (Fig. 7) elucidates the relative contribution of each cell (sub)type probability to the variances observed in Fig. 6g, i. We observe, for instance, with both types of variance, a large contribution from the variance in ependymal cell type probabilities (Fig. 7a, b). Specifically between replicates, excitatory neurons exhibit the second largest variance in cell type probabilities, as evidenced in the area of the dentate gyrus (Fig. 7c-e). With regard to spatial variance, we see specific astrocyte subtype probabilities with large variance, exhibited particularly in the areas of CA1 (yellow arrow) and the pons (orange arrow) both medially to laterally within the regions and when considering left versus right hemispheres (Fig. 7f, g). In contrast, other astrocyte subtype probabilities (Fig. 7h) are seen to be consistent over space, in line with our assumption of homogeneity within each CCFv3 region.

## Extension to additional tissue modalities for within and across species comparison

Constructed across institutions with varying combinations of molecular, chemical, genetic, and electrophysiological signals, multiple atlases per species now exist with different levels of granularity and intended applications¹,²,⁴⁹⁻⁵². Given this plethora of atlases, each of which might define a different partitioning scheme over the same area of tissue, questions of comparison and relevance of each atlas to emerging molecular and cellular signatures naturally arise⁵³. While some atlases have been defined in the same coordinate framework—often achieved through existing methods of image registration or manual alignment⁵³—many exist in different coordinate frameworks. Together with mismatches in the number, type, and positioning of partitions, this poses a challenge not only to the evaluation of each atlas ontology's fit to a molecular target but also to the ready comparison of atlas to atlas and the establishment of a clear metric of similarity between them. Hence, a second application of our method rests in its use to map not just atlas to the molecular dataset, but one atlas to another, with the geometric and functional correspondences $(\varphi, \pi)$ yielded by our method serving as an anchor for cross-examination of existing ontologies both within and across species.

To map tissue-scale atlas to tissue-scale atlas, we model both objects as image varifolds, as described in "Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms" section, but where both template and target are constrained to be of constant density ($w_i=1$ for all $x_i$ in foreground tissue). Both geometric transformations and latent feature distributions, $(\varphi, \pi)$ are estimated, with the latter giving distributions over the target atlas features (labels) for each template atlas label, following the assumption of spatial homogeneity with each template atlas label presumed to map onto a single distribution of target atlas labels.

The joint estimation of geometric transformation, $\varphi$ and conditional feature laws, $(\pi_t)_{t \in \mathcal{L}}$ in xIV-LDDMM offers two modes of quantitative comparison between atlas ontologies. First, as in the classical image setting of LDDMM, the determinant of the Jacobian, $|D\varphi|$, of the estimated diffeomorphism, can be used as a metric of how similar the atlas ontologies are, reflective of how much boundaries of partitions move to maximize the overlap between homogeneous regions. However, unlike in classical image settings, the estimation of the additional family of feature laws here, $\pi$, affords a second metric of similarity with
<figure><img src="image_11.png" /><figcaption>**Fig. 7** | Variance in individual cell subtype probabilities across replicates and across space within individual CCFv3 regions. **a** Variance in estimated cell subtype probabilities per CCFv3 region across three replicates summed overall CCFv3 regions. **b** Spatial variance of cell subtype probabilities, estimated empirically from three pulled-back MERFISH sections, per CCFv3 region and summed across all regions. **c–e** Probability of excitatory neuron subtype 2 (star in **a**) for each of the three mice in CCFv3 coordinates. Yellow arrow highlights the area of the dentate gyrus with differences in excitatory neuron subtype 2 probabilities. **f–h** Probability</figcaption></figure>

for astrocyte subtypes 1, 2, and 3 (stars in **b**) in empirical distribution computed from all three mice in CCFv3 coordinates (most likely cell type shown in Fig. 6h). Yellow arrow highlights area of CA1 with differences in astrocyte probability medially to laterally in subtypes 1 and 2 but not 3. Orange arrow highlights differences medially to laterally and left and right in areas of the pons in astrocyte probability for subtypes 1 and 2 but not 3. Black lines indicate boundaries between CCFv3 regions.

the computation of the entropy of the estimated conditional feature distributions.

The first panel in Fig. 8 shows the results of mapping one mouse atlas ontology to another with the Z section 680 in the CCFv3 mapped to the corresponding section in the Kim Lab Developmental atlas, termed Developmental Common Coordinate Framework (DevCCF)⁵¹, (Fig. 8a–c and vice versa Fig. 8d–f). Allen CCFv3 and DevCCF anatomical delineations utilize two distinct ontologies based on cytoarchitecture and genoarchitecture, respectively¹,⁵¹,⁵⁴. Fig. 8a, d depict the geometry of the section under each ontology, with the CCFv3 section hosting ≈140 independent regions and the DevCCF section ≈80. In this setting, both atlases are published in the same coordinate framework, giving $\phi = \text{Id}$ and thus, highlighting, instead, the estimated distributions over the other ontologies. Figure 8b, e depict the estimated conditional probability distributions, $p'_{i}$, $i \in I$, for each atlas section

over the other atlas section's ontology. The label with the highest probability in these distributions is plotted for each simplex in the mesh and which is consistent across each partition of each original atlas, given the homogeneity assumption in our model (i.e. a single $\pi_{\ell}$ for each $\ell \in \mathcal{L}$). The comparatively larger set of labels in the CCFv3 ontology results in labels being omitted from the corresponding estimated set of labels on the DevCCF ontology section Fig. 8e in the isocortex while multiple regions in the CCFv3 ontology carry the same most probable region in the DevCCF ontology. The differences in granularity with which each atlas segments the various areas of tissue are captured by the entropy of the estimated conditional feature distributions, for each simplex of the mesh (Fig. 8c, f). The entropy of the distributions estimated for the DevCCF ontology over the CCFv3 ontology Fig. 8f is on average, higher, than that of the distributions estimated for the CCFv3 ontology Fig. 8c, with probability mass
<figure><img src="image_13.png" /><figcaption>**Fig. 7** | Variance in individual cell subtype probabilities across replicates and across space within individual CCFv3 regions. a Variance in estimated cell subtype probabilities per CCFv3 region across three replicates summed overall CCFv3 regions. b Spatial variance of cell subtype probabilities, estimated empirically from three pulled-back MERFISH sections, per CCFv3 region and summed across all regions. c–e Probability of excitatory neuron subtype 2 (star in a) for each of the three mice in CCFv3 coordinates. Yellow arrow highlights the area of the dentate gyrus with differences in excitatory neuron subtype 2 probabilities. f–h Probability</figcaption></figure>

for astrocyte subtypes 1, 2, and 3 (stars in **b**) in empirical distribution computed from all three mice in CCFv3 coordinates (most likely cell type shown in Fig. 6h). Yellow arrow highlights area of CA1 with differences in astrocyte probability medially to laterally in subtypes 1 and 2 but not 3. Orange arrow highlights differences medially to laterally and left and right in areas of the pons in astrocyte probability for subtypes 1 and 2 but not 3. Black lines indicate boundaries between CCFv3 regions.

the computation of the entropy of the estimated conditional feature distributions.

The first panel in Fig. 8 shows the results of mapping one mouse atlas ontology to another with the Z section 680 in the CCFv3 mapped to the corresponding section in the Kim Lab Developmental atlas, termed Developmental Common Coordinate Framework (DevCCF)⁵¹, (Fig. 8a–c and vice versa Fig. 8d–f). Allen CCFv3 and DevCCF anatomical delineations utilize two distinct ontologies based on cytoarchitecture and genoarchitecture, respectively¹,⁵¹,⁵⁴. Fig. 8a, d depict the geometry of the section under each ontology, with the CCFv3 section hosting ≈140 independent regions and the DevCCF section ≈80. In this setting, both atlases are published in the same coordinate framework, giving $\varphi = \text{Id}$ and thus, highlighting, instead, the estimated distributions over the other ontologies. Figure 8b, e depict the estimated conditional probability distributions, $p'_{i}$, $i \in I$, for each atlas section

over the other atlas section's ontology. The label with the highest probability in these distributions is plotted for each simplex in the mesh and which is consistent across each partition of each original atlas, given the homogeneity assumption in our model (i.e. a single $\pi_{\ell}$ for each $\ell \in \mathcal{L}$). The comparatively larger set of labels in the CCFv3 ontology results in labels being omitted from the corresponding estimated set of labels on the DevCCF ontology section Fig. 8e in the isocortex while multiple regions in the CCFv3 ontology carry the same most probable region in the DevCCF ontology. The differences in granularity with which each atlas segments the various areas of tissue are captured by the entropy of the estimated conditional feature distributions, for each simplex of the mesh (Fig. 8c, f). The entropy of the distributions estimated for the DevCCF ontology over the CCFv3 ontology Fig. 8f is on average, higher, than that of the distributions estimated for the CCFv3 ontology Fig. 8c, with probability mass
<figure><img src="image_12.png" /><figcaption>**Fig. 8 |** Comparison of corresponding sections in CCFv3 and DevCCF atlases both independently and in the context of being mapped to a cell-based MERFISH target. **a, d** Original CCFv3 and DevCCF ontologies at location *Z* = 680 out of 1320. **b** CCFv3 geometry with predicted DevCCF atlas ontology. Delineations of original CCFv3 partitions are outlined in gray. **e** DevCCF atlas geometry with predicted CCFv3 ontology. Delineations of original DevCCF partitions are outlined in gray. **c, f** Entropy of predicted ontologies, with higher entropy values (light indicating less 1:1 correspondence between ontologies. **g** Predicted cell type with highest probability per simplex in CCFv3 atlas following mapping to MERFISH</figcaption></figure>

target (shown in Fig. 3) with XIV-LDDMM. **j** Predicted cell type with a highest probability per simplex in DevCCF atlas following mapping to same MERFISH target with XIV-LDDMM. **h, k** Estimated geometric transformation, $\phi_1$, in each setting applied to each atlas, with areas of expansion (red) and contraction (blue) as measured by the determinant of the Jacobian. White arrow highlights differences in ontologies in amygdala and striatum designation leading to different geometric transformations. **i, l** Entropy of estimated cell type distribution per simplex in atlas. Circled area of the hippocampus highlights differences in atlas ontologies leading to differences in the estimated entropy of cell type distributions.

distributed across ~5–7 different CCFv3 regions for each DevCCF region of the isocortex. Nevertheless, we see close to 1:1 correspondence between CCFv3 and DevCCF labels in the center section of the slice, where the entropy of the estimated distributions is near 0.

Atlas comparison can be conducted both independently and within the context of particular molecular targets. For example, here, we compare corresponding sections of the CCFv3 and DevCCF atlases⁵¹ via mapping each section to the same cell-segmented

MERFISH section and comparing both estimated diffeomorphisms and cell type distributions (Fig. 8g–l). Cell type distributions are compared via visualization of the cell type with the highest probability within each region (Fig. 8g, j) and the overall entropy of the estimated distribution per region (Fig. 8i, l). The areas of the hippocampus (dashed circle) and striatum and amygdala (arrow) are partitioned with different levels of granularity. This leads to different optimal geometric transformations, as characterized by the determinant of the Jacobian
<figure><img src="image_14.png" /><figcaption>**Fig. 8 |** Comparison of corresponding sections in CCFv3 and DevCCF atlases both independently and in the context of being mapped to a cell-based MERFISH target. **a, d** Original CCFv3 and DevCCF ontologies at location *Z* = 680 out of 1320. **b** CCFv3 geometry with predicted DevCCF atlas ontology. Delineations of original CCFv3 partitions are outlined in gray. **c, f** DevCCF atlas geometry with predicted CCFv3 ontology. Delineations of original DevCCF partitions are outlined in gray. **g** Entropy of predicted ontologies, with higher entropy values (light indicating less 1:1 correspondence between ontologies). **h** Predicted cell type with highest probability per simplex in CCFv3 atlas following mapping to MERFISH</figcaption></figure>

target (shown in Fig. 3) with XIV-LDDMM. **j** Predicted cell type with a highest probability per simplex in DevCCF atlas following mapping to same MERFISH target with XIV-LDDMM. **h, k** Estimated geometric transformation, $\phi_1$, in each setting applied to each atlas, with areas of expansion (red) and contraction (blue) as measured by the determinant of the Jacobian. White arrow highlights differences in ontologies in amygdala and striatum designation leading to different geometric transformations. **i, l** Entropy of estimated cell type distribution per simplex in atlas. Circled area of the hippocampus highlights differences in atlas ontologies leading to differences in the estimated entropy of cell type distributions.

distributed across ~5–7 different CCFv3 regions for each DevCCF region of the isocortex. Nevertheless, we see close to 1:1 correspondence between CCFv3 and DevCCF labels in the center section of the slice, where the entropy of the estimated distributions is near 0.

Atlas comparison can be conducted both independently and within the context of particular molecular targets. For example, here, we compare corresponding sections of the CCFv3 and DevCCF atlases⁵¹ via mapping each section to the same cell-segmented

MERFISH section and comparing both estimated diffeomorphisms and cell type distributions (Fig. 8g–l). Cell type distributions are compared via visualization of the cell type with the highest probability within each region (Fig. 8g, j) and the overall entropy of the estimated distribution per region (Fig. 8i, l). The areas of the hippocampus (dashed circle) and striatum and amygdala (arrow) are partitioned with different levels of granularity. This leads to different optimal geometric transformations, as characterized by the determinant of the Jacobian
(Fig. 8b, e), and different predicted cell type distributions (Fig. 8c, f). As discussed above, though both atlases are published as geometrically aligned⁵¹, the diffeomorphism solving the variational problem, here, transforms geometrically the homogeneous regions between the atlas and target. Hence, regions of the amygdala and striatum undergo significant contraction in the optimal mapping of DevCCF but not CCFv3 to MERFISH given the partitioning of this region into fewer and thus larger presumed homogeneous regions in the DevCCF atlas. Just as entropy in estimated distribution can be used to compare atlas to atlas, directly, we examine it here comparatively between atlases as an indication of which regions in which atlas achieve more or less homogeneous cell type distributions (Fig. 8i, l). Here, the hippocampus is more finely partitioned in the Kim atlas, which yields lower entropy distributions over cell types than in those estimated for the CCFv3.

Atlas ontologies can be mapped not just within species but also across them, where both geometric transformations and estimated ontology distributions, together reflect metrics of comparison between the two. As an example, we map a coronal section, *Z* = 537, in the CCFv3 to a coronal section, *Z* = 628 in the Waxholm Rat Brain Atlas⁵², with both sections chosen to correspond as sections through the anterior commissure (Fig. 9). Here, the CCFv3 ontology is comprised of ≈120 regions (Fig. 9a) with ≈30 for the Waxholm atlas section (Fig. 9d). Initial differences in size and shape exist between the two tissue sections (Fig. 9b). After scaling the volume of the mouse brain by 1.5, additional deformation, with magnitude given by the determinant of the Jacobian, [*Dφ*], distorts both internal and external tissue boundaries to align homogeneous regions in each atlas, such as cingulate area to cingulate area (white arrow). Estimated distributions over the Waxholm ontology labels for each region in the CCFv3 are shown in Fig. 9c, f, summarized by the maximum probability label (Fig. 9f) and measures of entropy (Fig. 9c), which highlight in gray,

CCFv3 regions mapping to ≈3-4 Waxholm regions versus those in black achieving 1:1 correspondence.

## Extension to additional cellular-scale modalities for image segmentation

With the computational constraints imposed by modeling spatial transcriptomics datasets as classical images, this work has emphasized the mapping of tissue-scale atlases to molecular datasets generated by emerging image-based rather than the often regularized spot-based spatial transcriptomics technologies. However, as exhibited through the modeling of atlas sections as image varifolds, themselves, both as template and target (see "Extension to additional tissue modalities for within and across species comparison" section), the xIV-LDDMM is equally capable of aligning tissue-scale atlases to alternative molecular and cellular-scale modalities. As an example, we take a DAPI-stained tissue section (Fig. 10), digitized at 2.5 µm resolution, corresponding to the gene-based MERFISH section shown in Fig. 1. The DAPI image (-13 million pixels) is converted to an image-varifold particle representation by discretization of its image values into -35 discrete bins, and selection of foreground pixels with corresponding image values in the later 30 of these bins. Approximately 250k particles are used to represent these pixels, with each particle pertaining to a set of 25 neighboring pixels (5 × 5 square). Particles each carry the distribution of bins into which the image values of these 25 pixels fall, retaining the individual values of each foreground pixel at the highest resolution in contrast to typical image downsampling, in which a pixel only captures the single mean image value of its neighbors.

Figure 10 depicts this particle representation of the DAPI image and its positioning before and after alignment with the corresponding CCFv3 section. Similar to the BARseq mapping illustrated in Fig. 5, we see areas of lower cell density (fewer foreground pixels) versus higher cell density (greater and higher intensity foreground pixels) aligning to

<figure><img src="image_13.png" /><figcaption>**Fig. 9** | **Results of mapping coronal section Z = 537 of CCFv3 mouse atlas to corresponding coronal section of Waxholm Rat Brain Atlas at Z = 628, both chosen to be through the anterior commissure.** a CCFv3 section in CCFv3 coordinates. **b** Waxholm section in Waxholm atlas coordinates. **c** Initial tissue overlap between mouse and rat atlas sections shown in Waxholm atlas coordinates. **d** Resulting overlap between rat section (blue) and mouse section following action of estimated diffeomorphism on mouse section. The determinant of the Jacobian</figcaption></figure>

highlights areas of expansion (red) and contraction (blue) in the mouse section deforming to match the rat section, with white arrow highlighting expansion in the cingulate area needed to match the region in the mouse to the corresponding region in the rat. **e** Entropy for each mouse region's predicted distribution of rat labels. **f** Predicted rat label with the highest probability for each region in CCFv3 mouse ontology.
(Fig. 8b, e), and different predicted cell type distributions (Fig. 8c, f). As discussed above, though both atlases are published as geometrically aligned⁵¹, the diffeomorphism solving the variational problem, here, transforms geometrically the homogeneous regions between the atlas and target. Hence, regions of the amygdala and striatum undergo significant contraction in the optimal mapping of DevCCF but not CCFv3 to MERFISH given the partitioning of this region into fewer and thus larger presumed homogeneous regions in the DevCCF atlas. Just as entropy in estimated distribution can be used to compare atlas to atlas, directly, we examine it here comparatively between atlases as an indication of which regions in which atlas achieve more or less homogeneous cell type distributions (Fig. 8i, l). Here, the hippo-campus is more finely partitioned in the Kim atlas, which yields lower entropy distributions over cell types than in those estimated for the CCFv3.

Atlas ontologies can be mapped not just within species but also across them, where both geometric transformations and estimated ontology distributions, together reflect metrics of comparison between the two. As an example, we map a coronal section, *Z* = 537, in the CCFv3 to a coronal section, *Z* = 628 in the Waxholm Rat Brain Atlas⁵², with both sections chosen to correspond as sections through the anterior commissure (Fig. 9). Here, the CCFv3 ontology is comprised of ≈120 regions (Fig. 9a) with ≈30 for the Waxholm atlas section (Fig. 9d). Initial differences in size and shape exist between the two tissue sections (Fig. 9b). After scaling the volume of the mouse brain by 1.5, additional deformation, with magnitude given by the determinant of the Jacobian, [*Dφ*], distorts both internal and external tissue boundaries to align homogeneous regions in each atlas, such as cingulate area to cingulate area (white arrow). Estimated distributions over the Waxholm ontology labels for each region in the CCFv3 are shown in Fig. 9c, f, summarized by the maximum probability label (Fig. 9f) and measures of entropy (Fig. 9c), which highlight in gray,

CCFv3 regions mapping to ≈3-4 Waxholm regions versus those in black achieving 1:1 correspondence.

## Extension to additional cellular-scale modalities for image segmentation

With the computational constraints imposed by modeling spatial transcriptomics datasets as classical images, this work has emphasized the mapping of tissue-scale atlases to molecular datasets generated by emerging image-based rather than the often regularized spot-based spatial transcriptomics technologies. However, as exhibited through the modeling of atlas sections as image varifolds, themselves, both as template and target (see "Extension to additional tissue modalities for within and across species comparison" section), the xIV-LDDMM is equally capable of aligning tissue-scale atlases to alternative molecular and cellular-scale modalities. As an example, we take a DAPI-stained tissue section (Fig. 10), digitized at 2.5 µm resolution, corresponding to the gene-based MERFISH section shown in Fig. 1. The DAPI image (-13 million pixels) is converted to an image-varifold particle representation by discretization of its image values into -35 discrete bins, and selection of foreground pixels with corresponding image values in the later 30 of these bins. Approximately 250k particles are used to represent these pixels, with each particle pertaining to a set of 25 neighboring pixels (5 × 5 square). Particles each carry the distribution of bins into which the image values of these 25 pixels fall, retaining the individual values of each foreground pixel at the highest resolution in contrast to typical image downsampling, in which a pixel only captures the single mean image value of its neighbors.

Figure 10 depicts this particle representation of the DAPI image and its positioning before and after alignment with the corresponding CCFv3 section. Similar to the BARseq mapping illustrated in Fig. 5, we see areas of lower cell density (fewer foreground pixels) versus higher cell density (greater and higher intensity foreground pixels) aligning to

<figure><img src="image_15.png" /><figcaption>**Fig. 9** | **Results of mapping coronal section Z = 537 of CCFv3 mouse atlas to corresponding coronal section of Waxholm Rat Brain Atlas at Z = 628, both chosen to be through the anterior commissure.** a CCFv3 section in CCFv3 coordinates. **b** Waxholm section in Waxholm atlas coordinates. **c** Initial tissue overlap between mouse and rat atlas sections shown in Waxholm atlas coordinates. **d** Resulting overlap between rat section (blue) and mouse section following action of estimated diffeomorphism on mouse section. The determinant of the Jacobian</figcaption></figure>

highlights areas of expansion (red) and contraction (blue) in the mouse section deforming to match the rat section, with white arrow highlighting expansion in the cingulate area needed to match the region in the mouse to the corresponding region in the rat. **e** Entropy for each mouse region's predicted distribution of rat labels. **f** Predicted rat label with the highest probability for each region in CCFv3 mouse ontology.
<figure><img src="image_14.png" /><figcaption>**Fig. 10** | **Application of image-varifold-based method for mapping CCFv3 section to DAPI-stained image of tissue corresponding to MERFISH section.** a Original DAPI-stained image, digitized at 2.5 µm resolution for tissue section measured with MERFISH technology (Fig. 1). c Image-varifold particle representation (black points) of DAPI-stained image overlaying the corresponding CCFv3 section in their respective initial coordinate spaces. Thresholded foreground pixels from (a) converted to particle image-varifold representation over a feature space of -30 binned grayscale values. d Alignment of CCFv3 section to DAPI</figcaption></figure>

particles following diffeomorphic transformation to the DAPI coordinate space.
b Alignment of CCFv3 section and DAPI particles in image format, with the deformed CCFv3 section image generated by resampling the deformed CCFv3 particles onto a regular 2.5 µm grid. White arrows highlight areas of alignment in the area of the substantia inominata (SI) and layer 1 of the cortex whereas red arrows highlight areas of questionable alignment in the area of the olfactory tubercle (OT).

expected CCFv3 regions (e.g. layer 1 (L1) especially between the hemispheres, and substantia inominata (SI), respectively). Interestingly, the area of the olfactory tubercle (OT) in the CCFv3 section which is initially quite larger than that in the DAPI image remains large following geometric transformation to the target space, with optimal cost minimization favoring this alignment to one in which this area would drastically contract.

Finally, while the estimated geometric transforms here could equivalently pull back the DAPI image into the space of the CCFv3 section, as done in the setting of atlas construction and cross-replicate comparison (see "Cross-replicate comparison and atlas construction" section), we show, instead, the CCFv3 section transformed to the space of the DAPI target as an illustration of how our methodology can be used in settings of image segmentation. This facilitates the extraction of measures in a particular subset of regions, that can be directly compared to other types of measures (e.g. gene expression, cell type distributions) as detected by other types of technologies and equivalently localized to these same CCFv3 regions through mappings computed with xIV-LDDMM, as demonstrated here in "Generalizability of xIV-LDDMM to settings of complete and incomplete gene-based and cell-based data," "Quantitative and comparative evaluation of xIV-LDDMM," and "Cross-replicate comparison and atlas construction" sections. For instance, a comparison of the mean DAPI signal (Supplementary Fig. 7) to the maximally expressed gene per CCFv3 region, as shown in Fig. 1, illustrates a positive correlation between DAPI intensity and *Whrn* expression, particularly in areas bordering the corpus callosum versus negative correlation with *Mdga1* expression in outer cortical layers. Finally, as described in the setting of resampling image varifolds across scales³⁴, we demonstrate the feasibility of translating between particle and image representations in resampling

the deformed CCFv3 particles with Gaussian smoothing onto the same 2.5 µm grid of the original DAPI image (Fig. 10b). Note that this smoothing extends the borders of the atlas section slightly beyond their respective positioning as particles, as seen in the difference in overlay laterally within L1 between image and particle representations (Fig. 10b, d, respectively).

## Discussion

We have introduced, here, a universal method, xIV-LDDMM, for mapping tissue-scale, 'cartoon' atlases to molecular and cellular datasets arising in the context of emerging transcriptomics technologies. We root our method in the modeling of each object as an image varifold, as previously described³⁵, and map across scale and modality by simultaneously estimating a geometric transformation with classical deformation tools of LDDMM¹⁷ and latent distribution over target features for each atlas region, under an assumption of homogeneity. Offering multiple means of implementation and optimization schemes (see "Alternating LDDMM and quadratic program algorithm for joint optimization" section), we have used both point cloud and mesh-based implementations of image-varifold objects, showcasing successful estimation of geometric mappings and latent feature distributions in both cases.

As presented, xIV-LDDMM fills a current need, as highlighted previously⁵⁵, for universal tools that can integrate the diverse types and large quantities of data emerging from the evolution of both transcriptomics and imaging technologies over the last decade. With each technology generating a slightly different perspective and a different set of animal or human samples to compare, a method that can stably handle the format of past, current, and future datasets will be paramount to integrate both new findings with the vast number of
<figure><img src="image_16.png" /><figcaption>**Fig. 10** | **Application of image-varifold-based method for mapping CCFv3 section to DAPI-stained image of tissue corresponding to MERFISH section.** a Original DAPI-stained image, digitized at 2.5 µm resolution for tissue section measured with MERFISH technology (Fig. 1). c Image-varifold particle representation (black points) of DAPI-stained image overlaying the corresponding CCFv3 section in their respective initial coordinate spaces. Thresholded foreground pixels from (a) converted to particle image-varifold representation over a feature space of -30 binned grayscale values. d Alignment of CCFv3 section to DAPI</figcaption></figure>

particles following diffeomorphic transformation to the DAPI coordinate space.
b Alignment of CCFv3 section and DAPI particles in image format, with the deformed CCFv3 section image generated by resampling the deformed CCFv3 particles onto a regular 2.5 µm grid. White arrows highlight areas of alignment in the area of the substantia inominata (SI) and layer 1 of the cortex whereas red arrows highlight areas of questionable alignment in the area of the olfactory tubercle (OT).

expected CCFv3 regions (e.g. layer 1 (L1) especially between the hemispheres, and substantia inominata (SI), respectively). Interestingly, the area of the olfactory tubercle (OT) in the CCFv3 section which is initially quite larger than that in the DAPI image remains large following geometric transformation to the target space, with optimal cost minimization favoring this alignment to one in which this area would drastically contract.

Finally, while the estimated geometric transforms here could equivalently pull back the DAPI image into the space of the CCFv3 section, as done in the setting of atlas construction and cross-replicate comparison (see "Cross-replicate comparison and atlas construction" section), we show, instead, the CCFv3 section transformed to the space of the DAPI target as an illustration of how our methodology can be used in settings of image segmentation. This facilitates the extraction of measures in a particular subset of regions, that can be directly compared to other types of measures (e.g. gene expression, cell type distributions) as detected by other types of technologies and equivalently localized to these same CCFv3 regions through mappings computed with xIV-LDDMM, as demonstrated here in "Generalizability of xIV-LDDMM to settings of complete and incomplete gene-based and cell-based data," "Quantitative and comparative evaluation of xIV-LDDMM," and "Cross-replicate comparison and atlas construction" sections. For instance, a comparison of the mean DAPI signal (Supplementary Fig. 7) to the maximally expressed gene per CCFv3 region, as shown in Fig. 1, illustrates a positive correlation between DAPI intensity and *Whrn* expression, particularly in areas bordering the corpus callosum versus negative correlation with *Mdga1* expression in outer cortical layers. Finally, as described in the setting of resampling image varifolds across scales³⁴, we demonstrate the feasibility of translating between particle and image representations in resampling

the deformed CCFv3 particles with Gaussian smoothing onto the same 2.5 µm grid of the original DAPI image (Fig. 10b). Note that this smoothing extends the borders of the atlas section slightly beyond their respective positioning as particles, as seen in the difference in overlay laterally within L1 between image and particle representations (Fig. 10b, d, respectively).

Discussion

We have introduced, here, a universal method, xIV-LDDMM, for mapping tissue-scale, 'cartoon' atlases to molecular and cellular datasets arising in the context of emerging transcriptomics technologies. We root our method in the modeling of each object as an image varifold, as previously described³⁵, and map across scale and modality by simultaneously estimating a geometric transformation with classical deformation tools of LDDMM¹⁷ and latent distribution over target features for each atlas region, under an assumption of homogeneity. Offering multiple means of implementation and optimization schemes (see "Alternating LDDMM and quadratic program algorithm for joint optimization" section), we have used both point cloud and mesh-based implementations of image-varifold objects, showcasing successful estimation of geometric mappings and latent feature distributions in both cases.

As presented, xIV-LDDMM fills a current need, as highlighted previously⁵⁵, for universal tools that can integrate the diverse types and large quantities of data emerging from the evolution of both transcriptomics and imaging technologies over the last decade. With each technology generating a slightly different perspective and a different set of animal or human samples to compare, a method that can stably handle the format of past, current, and future datasets will be paramount to integrate both new findings with the vast number of
datastores currently available across institutes. The image-varifold framework used here is general enough to model classical tissue-scale imaging data (as demonstrated in our atlas-to-atlas mappings), digital pathology (as demonstrated in our example with the DAPI-stained image) ("Extension to additional cellular-scale modalities for image segmentation" section), and emerging transcriptional data from both image-based and spot-resolution technologies that might generate data complete or partial tissue captures, irregularly sampled across space. Therefore, it provides a gateway for integrating data historically curated through MRI, immunohistochemistry, and other staining and imaging procedures in addition to the emerging transcriptomics methods. Importantly, this spatial integration will enable the correlation of disease signatures at different scales, as measured through different technologies. For instance, we aim to correlate tau tangle and amyloid-beta pathology in Alzheimer's disease to gene expression and imaging biomarkers, as vital for furthering our mechanistic understanding of the disease and developing early diagnostic strategies. Additionally, the specific estimation of diffeomorphic mappings compared with rigid+scale transformations decreases the spatial variance in feature distributions across regions, increasing the statistical power for detecting differences in feature distributions (e.g. gene expression, cell composition) not just across replicates within a single group, but between groups of replicates under different experimental conditions. For instance, we are currently looking at differences in cell type and gene distribution following neonatal binocular enucleation as measured with BARseq in four cases versus control hemispheres⁵⁶. As presented, xIV-LDDMM fills a current need, as highlighted previously⁵⁵, for universal tools that can integrate the diverse types and large quantities of data emerging from the evolution of both transcriptomics and imaging technologies over the last decade. With each technology generating a slightly different perspective and a different set of animal or human samples to compare, a method that can stably handle the format of past, current, and future datasets will be para-

set of animal or human samples to compare, a method that can stably
handle the format of past, current, and future datasets will be para-
mount to integrate both new findings with the vast number of datas-
tores currently available across institutes. The image-varifold
framework used here is general enough to model classical tissue-scale
imaging data (as demonstrated in our atlas-to-atlas mappings), digital
pathology (as demonstrated in our example with the DAPI-stained
image) ("Extension to additional cellular-scale modalities for image
segmentation" section), and emerging transcriptional data from both
image-based and spot-resolution technologies that might generate
data complete or partial tissue captures, irregularly sampled across
space. Therefore, it provides a gateway for integrating data historically
curated through MRI, immunohistochemistry, and other staining and
imaging procedures in addition to the emerging transcriptomics
methods. Importantly, this spatial integration will enable the correla-
tion of disease signatures at different scales, as measured through
different technologies. For instance, we aim to correlate tau tangle and
amyloid-beta pathology in Alzheimer's disease to gene expression and
imaging biomarkers, as vital for furthering our mechanistic under-
standing of the disease and developing early diagnostic strategies.
Additionally, the specific estimation of diffeomorphic mappings
compared with rigid+scale transformations decreases the spatial var-
iance in feature distributions across regions, increasing the statistical
power for detecting differences in feature distributions (e.g. gene
expression, cell composition) not just across replicates within a single
group, but between groups of replicates under different experimental
conditions. For instance, we are currently looking at differences in cell
type and gene distribution following neonatal binocular enucleation as
measured with BARseq in four cases versus control hemispheres⁵⁶. As
presented, xIV-LDDMM fills a current need, as highlighted previously⁵⁵,
for universal tools that can integrate the diverse types and large
quantities of data emerging from the evolution of both tran-
scriptomics and imaging technologies over the last decade. With each
technology generating a slightly different perspective and a different
set of animal or human samples to compare, a method that can stably
handle the format of past, current, and future datasets will be para-
mount to integrate both new findings with the vast number of datas-
tores currently available across institutes. The image-varifold
framework used here is general enough to model classical tissue-scale
imaging data (as demonstrated in our atlas-to-atlas mappings), digital
pathology (as demonstrated in our example with the DAPI-stained
image) ("Extension to additional cellular-scale modalities for image
segmentation" section), and emerging transcriptional data from both
image-based and spot-resolution technologies that might generate
data complete or partial tissue captures, irregularly sampled across
space. Therefore, it provides a gateway for integrating data historically
curated through MRI, immunohistochemistry, and other staining and
imaging procedures in addition to the emerging transcriptomics
methods. Importantly, this spatial integration will enable the correla-
tion of disease signatures at different scales, as measured through
different technologies. For instance, we aim to correlate tau tangle and
amyloid-beta pathology in Alzheimer's disease to gene expression and
imaging biomarkers, as vital for furthering our mechanistic under-
standing of the disease and developing early diagnostic strategies.
Additionally, the specific estimation of diffeomorphic mappings
compared with rigid+scale transformations decreases the spatial var-
iance in feature distributions across regions, increasing the statistical
power for detecting differences in feature distributions (e.g. gene
expression, cell composition) not just across replicates within a single
group, but between groups of replicates under different experimental
conditions. For instance, we are currently looking at differences in cell
type and gene distribution following neonatal binocular enucleation as
measured with BARseq in four cases versus control hemispheres⁵⁶.

In parallel to the development and dispersion of diverse molecular datasets, there has been continued development on the side of
datastores currently available across institutes. The image-varifold framework used here is general enough to model classical tissue-scale imaging data (as demonstrated in our atlas-to-atlas mappings), digital pathology (as demonstrated in our example with the DAPI-stained image) ("Extension to additional cellular-scale modalities for image segmentation" section), and emerging transcriptional data from both image-based and spot-resolution technologies that might generate data complete or partial tissue captures, irregularly sampled across space. Therefore, it provides a gateway for integrating data historically curated through MRI, immunohistochemistry, and other staining and imaging procedures in addition to the emerging transcriptomics methods. Importantly, this spatial integration will enable the correlation of disease signatures at different scales, as measured through different technologies. For instance, we aim to correlate tau tangle and amyloid-beta pathology in Alzheimer's disease to gene expression and imaging biomarkers, as vital for furthering our mechanistic understanding of the disease and developing early diagnostic strategies. Additionally, the specific estimation of diffeomorphic mappings compared with rigid+scale transformations decreases the spatial variance in feature distributions across regions, increasing the statistical power for detecting differences in feature distributions (e.g. gene expression, cell composition) not just across replicates within a single group, but between groups of replicates under different experimental conditions. For instance, we are currently looking at differences in cell type and gene distribution following neonatal binocular enucleation as measured with BARseq in four cases versus control hemispheres⁵⁶. As presented, xIV-LDDMM fills a current need, as highlighted previously⁵⁵, for universal tools that can integrate the diverse types and large quantities of data emerging from the evolution of both transcriptomics and imaging technologies over the last decade. With each technology generating a slightly different perspective and a different set of animal or human samples to compare, a method that can stably handle the format of past, current, and future datasets will be para-

set of animal or human samples to compare, a method that can stably
handle the format of past, current, and future datasets will be para-
mount to integrate both new findings with the vast number of datas-
tores currently available across institutes. The image-varifold
framework used here is general enough to model classical tissue-scale
imaging data (as demonstrated in our atlas-to-atlas mappings), digital
pathology (as demonstrated in our example with the DAPI-stained
image) ("Extension to additional cellular-scale modalities for image
segmentation" section), and emerging transcriptional data from both
image-based and spot-resolution technologies that might generate
data complete or partial tissue captures, irregularly sampled across
space. Therefore, it provides a gateway for integrating data historically
curated through MRI, immunohistochemistry, and other staining and
imaging procedures in addition to the emerging transcriptomics
methods. Importantly, this spatial integration will enable the correla-
tion of disease signatures at different scales, as measured through
different technologies. For instance, we aim to correlate tau tangle and
amyloid-beta pathology in Alzheimer's disease to gene expression and
imaging biomarkers, as vital for furthering our mechanistic under-
standing of the disease and developing early diagnostic strategies.
Additionally, the specific estimation of diffeomorphic mappings
compared with rigid+scale transformations decreases the spatial var-
iance in feature distributions across regions, increasing the statistical
power for detecting differences in feature distributions (e.g. gene
expression, cell composition) not just across replicates within a single
group, but between groups of replicates under different experimental
conditions. For instance, we are currently looking at differences in cell
type and gene distribution following neonatal binocular enucleation as
measured with BARseq in four cases versus control hemispheres⁵⁶. As
presented, xIV-LDDMM fills a current need, as highlighted previously⁵⁵,
for universal tools that can integrate the diverse types and large
quantities of data emerging from the evolution of both tran-
scriptomics and imaging technologies over the last decade. With each
technology generating a slightly different perspective and a different
set of animal or human samples to compare, a method that can stably
handle the format of past, current, and future datasets will be para-
mount to integrate both new findings with the vast number of datas-
tores currently available across institutes. The image-varifold
framework used here is general enough to model classical tissue-scale
imaging data (as demonstrated in our atlas-to-atlas mappings), digital
pathology (as demonstrated in our example with the DAPI-stained
image) ("Extension to additional cellular-scale modalities for image
segmentation" section), and emerging transcriptional data from both
image-based and spot-resolution technologies that might generate
data complete or partial tissue captures, irregularly sampled across
space. Therefore, it provides a gateway for integrating data historically
curated through MRI, immunohistochemistry, and other staining and
imaging procedures in addition to the emerging transcriptomics
methods. Importantly, this spatial integration will enable the correla-
tion of disease signatures at different scales, as measured through
different technologies. For instance, we aim to correlate tau tangle and
amyloid-beta pathology in Alzheimer's disease to gene expression and
imaging biomarkers, as vital for furthering our mechanistic under-
standing of the disease and developing early diagnostic strategies.
Additionally, the specific estimation of diffeomorphic mappings
compared with rigid+scale transformations decreases the spatial var-
iance in feature distributions across regions, increasing the statistical
power for detecting differences in feature distributions (e.g. gene
expression, cell composition) not just across replicates within a single
group, but between groups of replicates under different experimental
conditions. For instance, we are currently looking at differences in cell
type and gene distribution following neonatal binocular enucleation as
measured with BARseq in four cases versus control hemispheres⁵⁶.

In parallel to the development and dispersion of diverse molecular datasets, there has been continued development on the side of
reference atlases to reflect trends in these new measures and integrate these trends across even more samples of particular species. With regard to atlas refinement and creation, as shown in “Cross-replicate comparison and atlas construction” section and the second panel of Fig. 1, the invertibility of the estimated diffeomorphism in the setting of mapping atlas to molecular target, enables the carrying of each target into the same coordinate space of the atlas. Though we illustrate an example of atlas creation from samples of a single modality, each target to which the same atlas section is mapped can be pulled back into the same coordinate space. This effectively allows different molecular modalities to be *indirectly* mapped to one another, where the challenge of matching different molecular features as needed to map one modality directly to another is facilitated through the scaffold of the presumed homogeneous regions in the atlas that anchor the estimated feature distributions for alignment to each respective molecular modality. Furthermore, this enables the construction of a composite set of molecular and cellular-scale data across technologies and modalities (e.g. digital pathology and spatial transcriptomics, as seen in the comparative mappings of the same CCFv3 section to MERFISH (Fig. 1) and DAPI images (Fig. 10)), that can be integrated into a multi-factorial atlas with defined segmentation schemes according to homogeneous regions across these features. Our method also offers a tool for re-examining and comparing existing atlas ontologies in the context of new data⁵⁵, and serves as a means for developing new atlases in the future. As described in “Extension to additional tissue modalities for within and across species comparison” section, examination of the mappings achieved between different atlases and the same molecular target offers an indirect comparison between atlases in the context of a particular molecular setting. However, this comparison can also be made directly in a context-independent setting by harnessing our method to map atlas to atlas. In the field of evolutionary biology, for instance, our method could aid in the mapping and comparison of atlases across species⁵⁷,⁵⁸ and in the field of developmental neurobiology, the available atlases of the brain at different stages of development⁴⁹⁻⁵¹,⁵⁴,⁵⁹. Finally, as shown in “Extension to additional cellular-scale modalities for image segmentation” section, xIV-LDDMM can be applied across molecular and cellular modalities, not just for incorporating them into atlas creation, but also as a means of segmentation in the classical image sense.

While the results presented here survey a wide variety of potential applications of xIV-LDDMM to mapping atlas modalities to diverse targets, there remain uncertainties and potential modes of improvement that are the subject of current and future work. First, we have presented results mapping 2D sections of 3D atlases to corresponding 2D sections of MERFISH data. Both the Allen MERFISH and BARseq data showcased here are part of an entire set of serial sections that span the whole brain. While we have approximated these sections as strictly coronal in pairing them to respective coronal slices of the CCFv3 here, typically such datasets are generated by sectioning on an angled plane, which would be best approximated by the mapping of the entire set of sections as a 3D object to the 3D CCFv3. Consequently, we are optimizing our method to compute mappings of atlas to molecular targets in 3D, where both added dimensions and added magnitudes of data contribute to the theoretical and computational complexity of the problem. Indeed, with ≈6 billion individual transcripts measured across the span of the brain, treatment of this data as a regular lattice image would require on the order of 1000 billion voxels at 1 µm resolution, which is coarser than that needed even to resolve two molecules of mRNA. Hence, it becomes even more vital to treat such data in the particle setting, as presented here, where we capture the sparsity and irregularity of the data in modeling it effectively in its lowest dimension, as 6 billion individual particle measures. Second, though we have highlighted a range of gene-based and cell-based, whole and hemi-brain datasets achieved with image-based technologies, here, the future investigation includes the use of our

method to map data from additional types of technologies, such as the spot-resolution SlideSeq¹⁴, and in additional biological settings, where we might see further disruptions of tissue architecture (e.g. with tumors) and tissue types with varying levels of the organization and therefore varying levels to which homogeneous distributions over molecular features can be seen over physical space (e.g. in heart, breast, and lung).

Indeed, we finally emphasize that central to the model posed here is the underlying assumption that each compartment has a homogeneous distribution over molecular features that is stationary with respect to space. We make this assumption to govern the estimation of latent feature distributions for each atlas compartment to take the atlas into both the physical and feature space of the target. This is in contrast to learning-based approaches that estimate features independent of such assumptions to map atlas to target across feature (e.g. imaging) modalities but often require extensive training datasets and computational power to learn such features¹⁸,²⁰. This assumption stems from the inherent construction of atlases often to delineate regions of particular cell types or with a particular function, and thus, where we see a set of predominant cell types or gene types consistently across the region in the molecular scale data, as in Fig. 3. The successful alignment of tissue-scale atlas to varying molecular scale data demonstrated through the examples shown here supports the benefit of such an assumption in jointly estimating geometric transformations and latent feature distributions. However, as manifest in Fig. 2, where the expression of *Trp53i1* appears to be distributed along a decreasing gradient medial to lateral within the striatum, or explicitly through the spatial variance in empirical distribution shown for each CCFv3 region in Fig. 6, this assumption of spatial homogeneity does not always hold. The subsequent construction of empirical distributions from the pulling back of different molecular targets, as shown in Fig. 6 enables the estimation of feature distributions at a finer scale without the assumption of homogeneity and simultaneously, an evaluation of where and to what extent it holds, with notably, the use of diffeomorphic mappings achieving a further reduction in this spatial variance than rigid+scale transformations alone. Furthermore, the results presented here reflect a particular balance between expected deformation and this homogeneity assumption, imposed by the relative weighting of the separate terms in the cost function. Current work at controlling this balance further includes the addition of a term controlling the divergence of the vector field to the energy defined in the variational problem (10), which leads to solutions more robust to deformation within the interior of the tissue. Future work will also include a more rigorous evaluation of how well this homogeneity assumption holds across different tissue contexts and the effect the given balance between the two terms might have in different settings.

## Methods

### Construction of image-varifold representation for different modalities

As introduced in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, we can represent each image-varifold object as a point cloud (particles) or a triangulated mesh. In the first case, data is modeled as a collection of particles, each with a center $x_i \in \mathbb{R}^2$, and a measure over its feature space, $w_i p_i \in \mathcal{M}(\mathcal{F})$. Particles may be individual mRNA detections, cell centers, or image pixels at the highest resolution, with thresholding or exclusion of particular image values done in the latter case to extract foreground pixels only. Sets of particles can be downsampled to a smaller set of discrete particles by capturing neighboring particle distributions together into one particle, as used in the implementation of the DAPI image (“Extension to additional cellular-scale modalities for image segmentation” section). We refer to the resolution of the rendering to indicate the span of the neighborhood each individual particle captures with the coarser the
reference atlases to reflect trends in these new measures and integrate these trends across even more samples of particular species. With regard to atlas refinement and creation, as shown in “Cross-replicate comparison and atlas construction” section and the second panel of Fig. 1, the invertibility of the estimated diffeomorphism in the setting of mapping atlas to molecular target, enables the carrying of each target into the same coordinate space of the atlas. Though we illustrate an example of atlas creation from samples of a single modality, each target to which the same atlas section is mapped can be pulled back into the same coordinate space. This effectively allows different molecular modalities to be *indirectly* mapped to one another, where the challenge of matching different molecular features as needed to map one modality directly to another is facilitated through the scaffold of the presumed homogeneous regions in the atlas that anchor the estimated feature distributions for alignment to each respective molecular modality. Furthermore, this enables the construction of a composite set of molecular and cellular-scale data across technologies and modalities (e.g. digital pathology and spatial transcriptomics, as seen in the comparative mappings of the same CCFv3 section to MERFISH (Fig. 1) and DAPI images (Fig. 10)), that can be integrated into a multi-factorial atlas with defined segmentation schemes according to homogeneous regions across these features. Our method also offers a tool for re-examining and comparing existing atlas ontologies in the context of new data⁵⁵, and serves as a means for developing new atlases in the future. As described in “Extension to additional tissue modalities for within and across species comparison” section, examination of the mappings achieved between different atlases and the same molecular target offers an indirect comparison between atlases in the context of a particular molecular setting. However, this comparison can also be made directly in a context-independent setting by harnessing our method to map atlas to atlas. In the field of evolutionary biology, for instance, our method could aid in the mapping and comparison of atlases across species⁵⁷,⁵⁸ and in the field of developmental neurobiology, the available atlases of the brain at different stages of development⁴⁹⁻⁵¹,⁵⁴,⁵⁹. Finally, as shown in “Extension to additional cellular-scale modalities for image segmentation” section, xIV-LDDMM can be applied across molecular and cellular modalities, not just for incorporating them into atlas creation, but also as a means of segmentation in the classical image sense.

While the results presented here survey a wide variety of potential applications of xIV-LDDMM to mapping atlas modalities to diverse targets, there remain uncertainties and potential modes of improvement that are the subject of current and future work. First, we have presented results mapping 2D sections of 3D atlases to corresponding 2D sections of MERFISH data. Both the Allen MERFISH and BARseq data showcased here are part of an entire set of serial sections that span the whole brain. While we have approximated these sections as strictly coronal in pairing them to respective coronal slices of the CCFv3 here, typically such datasets are generated by sectioning on an angled plane, which would be best approximated by the mapping of the entire set of sections as a 3D object to the 3D CCFv3. Consequently, we are optimizing our method to compute mappings of atlas to molecular targets in 3D, where both added dimensions and added magnitudes of data contribute to the theoretical and computational complexity of the problem. Indeed, with ≈6 billion individual transcripts measured across the span of the brain, treatment of this data as a regular lattice image would require on the order of 1000 billion voxels at 1 µm resolution, which is coarser than that needed even to resolve two molecules of mRNA. Hence, it becomes even more vital to treat such data in the particle setting, as presented here, where we capture the sparsity and irregularity of the data in modeling it effectively in its lowest dimension, as 6 billion individual particle measures. Second, though we have highlighted a range of gene-based and cell-based, whole and hemi-brain datasets achieved with image-based technologies, here, the future investigation includes the use of our

method to map data from additional types of technologies, such as the spot-resolution SlideSeq¹⁴, and in additional biological settings, where we might see further disruptions of tissue architecture (e.g. with tumors) and tissue types with varying levels of the organization and therefore varying levels to which homogeneous distributions over molecular features can be seen over physical space (e.g. in heart, breast, and lung).

Indeed, we finally emphasize that central to the model posed here is the underlying assumption that each compartment has a homogeneous distribution over molecular features that is stationary with respect to space. We make this assumption to govern the estimation of latent feature distributions for each atlas compartment to take the atlas into both the physical and feature space of the target. This is in contrast to learning-based approaches that estimate features independent of such assumptions to map atlas to target across feature (e.g. imaging) modalities but often require extensive training datasets and computational power to learn such features¹⁸,²⁰. This assumption stems from the inherent construction of atlases often to delineate regions of particular cell types or with a particular function, and thus, where we see a set of predominant cell types or gene types consistently across the region in the molecular scale data, as in Fig. 3. The successful alignment of tissue-scale atlas to varying molecular scale data demonstrated through the examples shown here supports the benefit of such an assumption in jointly estimating geometric transformations and latent feature distributions. However, as manifest in Fig. 2, where the expression of *Trp53i1* appears to be distributed along a decreasing gradient medial to lateral within the striatum, or explicitly through the spatial variance in empirical distribution shown for each CCFv3 region in Fig. 6, this assumption of spatial homogeneity does not always hold. The subsequent construction of empirical distributions from the pulling back of different molecular targets, as shown in Fig. 6 enables the estimation of feature distributions at a finer scale without the assumption of homogeneity and simultaneously, an evaluation of where and to what extent it holds, with notably, the use of diffeomorphic mappings achieving a further reduction in this spatial variance than rigid+scale transformations alone. Furthermore, the results presented here reflect a particular balance between expected deformation and this homogeneity assumption, imposed by the relative weighting of the separate terms in the cost function. Current work at controlling this balance further includes the addition of a term controlling the divergence of the vector field to the energy defined in the variational problem (10), which leads to solutions more robust to deformation within the interior of the tissue. Future work will also include a more rigorous evaluation of how well this homogeneity assumption holds across different tissue contexts and the effect the given balance between the two terms might have in different settings.

## Methods

### Construction of image-varifold representation for different modalities

As introduced in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, we can represent each image-varifold object as a point cloud (particles) or a triangulated mesh. In the first case, data is modeled as a collection of particles, each with a center $x_i \in \mathbb{R}^2$, and a measure over its feature space, $w_i p_i \in \mathcal{M}(\mathcal{F})$. Particles may be individual mRNA detections, cell centers, or image pixels at the highest resolution, with thresholding or exclusion of particular image values done in the latter case to extract foreground pixels only. Sets of particles can be downsampled to a smaller set of discrete particles by capturing neighboring particle distributions together into one particle, as used in the implementation of the DAPI image (“Extension to additional cellular-scale modalities for image segmentation” section). We refer to the resolution of the rendering to indicate the span of the neighborhood each individual particle captures with the coarser the
resolution corresponding to larger and larger neighborhoods of measurements encapsulated in a single particle.

In the case of meshes, each mesh is built from a collection of vertices indexed by the set $I$, $\mathbf{x}=(x_i)_{i \in I}$ with each $x_i \in \mathbb{R}^2$. Each simplex in the mesh is defined from the vertices denoted as $\gamma(\mathbf{x})$ and is paired with a 3-tuple with components that index the vertices of the simplex, $(\gamma(\mathbf{x}), c = (c^1, c^2, c^3) \in \mathcal{I}')$ and determine the center $m(\mathbf{x}) = \frac{1}{3}(x_{c^1} + x_{c^2} + x_{c^3})$. Each triangle simplex is defined by

$$ \gamma(\mathbf{x}) = \left\{ y \in \mathbb{R}^2 : y = \sum_{k=1}^{3} a_k x_{c^k}, a_k \ge 0, \sum_{k=1}^{3} a_k = 1 \right\}, \quad (5) $$

with positive orientation and volume $|\gamma_c(\mathbf{x})| := \frac{1}{2} \| (x_{c^2} - x_{c^1}) \times (x_{c^3} - x_{c^1}) \| > 0$.

The total mesh $\tau$ is the collection of vertices $\mathbf{x}$, and simplices and centers $(\gamma_j(\mathbf{x}), c_j = (c_j^1, c_j^2, c_j^3), m_j(\mathbf{x}))_{j \in J}$, with the simplices indexed by the set $J$, and with the resolution determining the complexity as total numbers of vertices $|I|$ and the number of simplices $|I|$ in the mesh. To complete the image varifold we append to the mesh the spatial density defined over the area of each simplex: $\alpha = (\alpha_j)_{j \in J}$ and the field of probability laws $\mathbf{p} = (p_j)_{j \in J}$ on $\mathcal{F}$. Hence, we denote an image varifold implemented as a mesh similar to the normalized definition of (1) as:

$$ \mu_{\tau} = \sum_{j \in J} w_j(\delta_{m_j} \otimes p_j) \quad (6) $$

with $w_j = a_j|\gamma_j|$ giving the assumed constant density over the area of the simplex. Notably, diffeomorphisms act directly on the vertices, with the center $m_j$ of each simplex shifted according to the movement of each of its vertices. The determinant of the Jacobian is introduced to retain the same spatial density of the object before and after transformation, as in (2) but with it evaluated at the center of each simplex, and again assumed constant over the area of each simplex, approximated as the ratio in the simplex area after and before transformation:

$$ \varphi \cdot \mu_{\tau} = \sum_{j \in J} \alpha_j |D\varphi|_{m_j} |\gamma_j| (\delta_{\varphi(m_j)} \otimes p_j). \quad (7) $$

Additionally, in the setting of meshes, it is convenient to approximate the product $|D\varphi|_{m_j}|\gamma_j|$ by the new area of each simplex following transformation, $|\gamma_j|\varphi_1(\mathbf{x})|$.

Meshes mapped in this work were constructed using Delauney triangulation⁶⁰ on a grid defined over the support of the starting dataset with the size of each square dictated by the input resolution. Spatial density and conditional feature measures, $\alpha, \mathbf{p}$, were associated with the simplices of the mesh following the assignment of each individual data point (e.g. mRNA or cell read) into its single nearest simplex. Meshes were pruned of simplices that both contained fewer than 1 data point and existed outside the largest connected component of simplices containing at least one data point. In this manner, both for atlas images and transcriptomics datasets, resulting simplex meshes spanned the entire tissue foreground.

### Molecular scale varifold norm

We define the space of image varifolds $\mu \in W^*$ to have a norm $\| \cdot \|_{W^*}^2$, and transform the atlas coordinates onto the targets to minimize the norm. The space of varifold norms is associated to a reproducing kernel Hilbert space$^{34,61}$ (see (8) below) defined by the inner product of the space as $\langle \mu, \nu \rangle_{W^*}, \| \mu \|_{W^*}^2 = \langle \mu, \mu \rangle_{W^*}$. To specify the image-varifold norm for $\mu \in W^*$, $\| \cdot \|_{W^*}^2$, it suffices to provide the inner product between Diracs $(\delta_x \otimes \delta_f, \delta_x^* \otimes \delta_f^*)_{W^*} = K((x_f, x^*f))$. For any $\mu$ in (1)

then

$$ \| \mu \|_{W^*}^2 = \sum_{i,j \in J} \sum_{f, f' \in \mathcal{F}} w_i p_i(f) w_j p_j(f') K((x_i f), (x_j f')). \quad (8) $$

Throughout we use the kernel product $K((x_i f), (x^*f')) = K_1(x, x^*)K_2(f, f')$ chosen as a Gaussian over physical space $K_1(x, y) = \exp(-\frac{|x-x'|^2}{2\sigma^2})$ with $K_2(f, f') = 1$ if $f = f'$, 0 otherwise giving:

$$ \| \mu \|_{W^*}^2 = \sum_{j,k \in J} K_1(x_j, x_k) \sum_{f \in \mathcal{F}} w_j p_j(f) w_k p_k(f) \quad (9) $$

for the image-varifold object, as defined in (1) in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section.

### Variational problems

As described in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, the mapping variational problem constructs a diffeomorphism $\varphi: \mathbb{R}^2 \to \mathbb{R}^2$ and feature laws $(\pi_l)_{l \in \mathcal{L}}$ on $\mathcal{F}$ to carry the atlas image varifold onto the target, minimizing the varifold normed difference between them, with norm defined as in “Molecular scale varifold norm” section. We follow LDDMM as described in the image case$^{17}$, parameterizing the diffeomorphism with the smooth time-varying velocity field, $v_t, t \in [0, 1]$, as $\dot{\varphi}_t = v_t \circ \varphi_t$. This gives the variational problem between the atlas image varifold, $\mu_A$, with indexed locations, $I$, and the target image varifold, $\mu_T$, with indexed locations, $J$:

$$ \inf_{v \in L^2([0,1], \mathcal{V}), \quad \pi_l, \ell \in \mathcal{L}} \frac{1}{2} \int_0^1 \| v_t \|_{\mathcal{V}}^2 dt + \| \varphi_1 \cdot \mu_A^\pi - \mu_T \|_{W^*}^2 \quad (10) $$

$$ \text{with } \dot{\varphi}_t = v_t \circ \varphi_t, \quad \varphi_0 = 1d, $$

with $\mu_A^\pi$ depicting the feature-transformed atlas image varifold, as described in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, as $\sum_{i \in J} \delta_{x_i} \otimes w_i p_i^t$ with $w_i p_i^t = \sum_{\ell \in \mathcal{L}} w_\ell p_\ell(\ell) \pi_\ell$, the estimated distribution over target features for indexed location $i$ in the atlas image varifold. The space of smooth time-varying velocity fields giving the flow, $\varphi$, is defined as a reproducing kernel Hilbert space, equipped with norm, $\| \cdot \|_{\mathcal{V}}^2$, ensuring smoothness and invertibility of $\varphi$ as a diffeomorphism$^{62}$.

We solve (10) for optimal $\varphi, \pi$ through either single or alternating algorithms as described in “Alternating LDDMM and quadratic program algorithm for joint optimization” section, with the opportunity to impose priors on the estimated distributions, $\pi$, appropriate to the specific setting. For instance, we typically impose positivity on all values in estimated distributions: $\pi_l(f) \ge 0$ for all $f \in \mathcal{F}$ and $\ell \in \mathcal{L}$. Given prior knowledge of the spatial density of cells or mRNA of a target, we may also specify a range of values for the resulting spatial density to take through constraints on the total measures of features estimated:

$$ w^{\min} \le w_i^* \le w^{\max}, \quad i \in I. \quad (11) $$

This prior can be easily incorporated as a constraint in the setting of the quadratic program used in the alternating algorithmic approach for estimating optimal distributions $\pi$ (see “Alternating LDDMM and quadratic program algorithm for joint optimization” section). In the examples shown here, we typically take $w^{\min}$ to be the 5th percentile of values $w_j, j \in J$ of the target image varifold. Additionally, in the setting of mapping tissue-scale atlases to each other, as discussed in “Extension to additional tissue modalities for within and across species comparison” section, we impose the greater constraint of ensuring
resolution corresponding to larger and larger neighborhoods of measurements encapsulated in a single particle.

In the case of meshes, each mesh is built from a collection of vertices indexed by the set $I$, $\mathbf{x} = (x_i)_{i \in I}$ with each $x_i \in \mathbb{R}^2$. Each simplex in the mesh is defined from the vertices denoted as $\gamma(\mathbf{x})$ and is paired with a 3-tuple with components that index the vertices of the simplex, $(\gamma(\mathbf{x}), c = (c^1, c^2, c^3) \in \mathcal{I}')$ and determine the center $m(\mathbf{x}) = \frac{1}{3}(x_{c1} + x_{c2} + x_{c3})$. Each triangle simplex is defined by

$$ \gamma(\mathbf{x}) = \left\{ y \in \mathbb{R}^2 : y = \sum_{k=1}^{3} a_k x_{c_k}, a_k \ge 0, \sum_{k=1}^{3} a_k = 1 \right\}, \quad (5) $$

with positive orientation and volume $|\gamma_c(\mathbf{x})| := \frac{1}{2} \| (x_{c2} - x_{c1}) \times (x_{c3} - x_{c1}) \| > 0$.

The total mesh $\tau$ is the collection of vertices $\mathbf{x}$, and simplices and centers $(\gamma_j(\mathbf{x}), c_j = (c_j^1, c_j^2, c_j^3), m_j(\mathbf{x}))_{j \in J}$, with the simplices indexed by the set $J$, and with the resolution determining the complexity as total numbers of vertices $|I|$ and the number of simplices $|I|$ in the mesh. To complete the image varifold we append to the mesh the spatial density defined over the area of each simplex: $\alpha = (\alpha_j)_{j \in J}$ and the field of probability laws $\mathbf{p} = (p_j)_{j \in J}$ on $\mathcal{F}$. Hence, we denote an image varifold implemented as a mesh similar to the normalized definition of (1) as:

$$ \mu_{\tau} = \sum_{j \in J} w_j(\delta_{m_j} \otimes p_j) \quad (6) $$

with $w_j = a_j|\gamma_j|$ giving the assumed constant density over the area of the simplex. Notably, diffeomorphisms act directly on the vertices, with the center $m_j$ of each simplex shifted according to the movement of each of its vertices. The determinant of the Jacobian is introduced to retain the same spatial density of the object before and after transformation, as in (2) but with it evaluated at the center of each simplex, and again assumed constant over the area of each simplex, approximated as the ratio in the simplex area after and before transformation:

$$ \varphi \cdot \mu_{\tau} = \sum_{j \in J} \alpha_j |D\varphi|_{m_j} |\gamma_j| (\delta_{\varphi(m_j)} \otimes p_j). \quad (7) $$

Additionally, in the setting of meshes, it is convenient to approximate the product $|D\varphi|_{m_j}|\gamma_j|$ by the new area of each simplex following transformation, $|\gamma_j|\varphi_1(\mathbf{x})|$.

Meshes mapped in this work were constructed using Delauney triangulation⁶⁰ on a grid defined over the support of the starting dataset with the size of each square dictated by the input resolution. Spatial density and conditional feature measures, $\alpha, \mathbf{p}$, were associated with the simplices of the mesh following the assignment of each individual data point (e.g. mRNA or cell read) into its single nearest simplex. Meshes were pruned of simplices that both contained fewer than 1 data point and existed outside the largest connected component of simplices containing at least one data point. In this manner, both for atlas images and transcriptomics datasets, resulting simplex meshes spanned the entire tissue foreground.

### Molecular scale varifold norm

We define the space of image varifolds $\mu \in W^*$ to have a norm $\| \cdot \|_{W^*}^2$, and transform the atlas coordinates onto the targets to minimize the norm. The space of varifold norms is associated to a reproducing kernel Hilbert space$^{34,61}$ (see (8) below) defined by the inner product of the space as $\langle \mu, \nu \rangle_{W^*}, \| \mu \|_{W^*}^2 = \langle \mu, \mu \rangle_{W^*}$. To specify the image-varifold norm for $\mu \in W^*$, $\| \cdot \|_{W^*}^2$, it suffices to provide the inner product between Diracs $(\delta_x \otimes \delta_f, \delta_x^* \otimes \delta_f^*)_{W^*} = K((x_f, (x^*f^*))$. For any $\mu$ in (1)

then

$$ \| \mu \|_{W^*}^2 = \sum_{i,j \in J} \sum_{f, f' \in \mathcal{F}} w_i p_i(f) w_j p_j(f') K((x_i f), (x_j f')). \quad (8) $$

Throughout we use the kernel product $K((x_i f), (x^*f')) = K_1(x, x^*)K_2(f, f')$ chosen as a Gaussian over physical space $K_1(x, y) = \exp(-\frac{|x-x'|^2}{2\sigma^2})$ with $K_2(f, f') = 1$ if $f = f^*$, 0 otherwise giving:

$$ \| \mu \|_{W^*}^2 = \sum_{j, k \in J} K_1(x_j, x_k) \sum_{f \in \mathcal{F}} w_j p_j(f) w_k p_k(f) \quad (9) $$

for the image-varifold object, as defined in (1) in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section.

### Variational problems

As described in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, the mapping variational problem constructs a diffeomorphism $\varphi: \mathbb{R}^2 \to \mathbb{R}^2$ and feature laws $(\pi_l)_{l \in \mathcal{L}}$ on $\mathcal{F}$ to carry the atlas image varifold onto the target, minimizing the varifold normed difference between them, with norm defined as in “Molecular scale varifold norm” section. We follow LDDMM as described in the image case$^{17}$, parameterizing the diffeomorphism with the smooth time-varying velocity field, $v_t, t \in [0, 1]$, as $\dot{\varphi}_t = v_t \circ \varphi_t$. This gives the variational problem between the atlas image varifold, $\mu_A$, with indexed locations, $I$, and the target image varifold, $\mu_T$, with indexed locations, $J$:

$$ \inf_{v \in L^2([0,1], \mathcal{V}), \quad \pi_l, \ell \in \mathcal{L}} \frac{1}{2} \int_0^1 \| v_t \|_{\mathcal{V}}^2 dt + \| \varphi_1 \cdot \mu_A^\pi - \mu_T \|_{W^*}^2 \quad (10) $$

$$ \text{with } \dot{\varphi}_t = v_t \circ \varphi_t, \quad \varphi_0 = 1d, $$

with $\mu_A^\pi$ depicting the feature-transformed atlas image varifold, as described in “Data model and optimization problem: image varifolds and transformations for molecular scales based on varifold norms” section, as $\sum_{i \in J} \delta_{x_i} \otimes w_i p_i^*$ with $w_i p_i^* = \sum_{\ell \in \mathcal{L}} w_\ell p_\ell(\ell) \pi_\ell$, the estimated distribution over target features for indexed location $i$ in the atlas image varifold. The space of smooth time-varying velocity fields giving the flow, $\varphi$, is defined as a reproducing kernel Hilbert space, equipped with norm, $\| \cdot \|_{\mathcal{V}}^2$, ensuring smoothness and invertibility of $\varphi$ as a diffeomorphism$^{62}$.

We solve (10) for optimal $\varphi, \pi$ through either single or alternating algorithms as described in “Alternating LDDMM and quadratic program algorithm for joint optimization” section, with the opportunity to impose priors on the estimated distributions, $\pi$, appropriate to the specific setting. For instance, we typically impose positivity on all values in estimated distributions: $\pi_l(f) \ge 0$ for all $f \in \mathcal{F}$ and $\ell \in \mathcal{L}$. Given prior knowledge of the spatial density of cells or mRNA of a target, we may also specify a range of values for the resulting spatial density to take through constraints on the total measures of features estimated:

$$ w^{\min} \le w_i^* \le w^{\max}, \quad i \in I. \quad (11) $$

This prior can be easily incorporated as a constraint in the setting of the quadratic program used in the alternating algorithmic approach for estimating optimal distributions $\pi$ (see “Alternating LDDMM and quadratic program algorithm for joint optimization” section). In the examples shown here, we typically take $w^{\min}$ to be the 5th percentile of values $w_j, j \in J$ of the target image varifold. Additionally, in the setting of mapping tissue-scale atlases to each other, as discussed in “Extension to additional tissue modalities for within and across species comparison” section, we impose the greater constraint of ensuring