Projective LDDMM: Spatially Reconstructing a Story of Rostrally-Dominant Tau in Alzheimer’s Disease

Since Braak’s initial histological observations, it has been recognized that Alzheimer’s disease (AD) neurofibrillary tangles (NFTs) appear in the medial temporal lobe (MTL) of the brain very early in the disease course. MRI-based shape diffeomorphometry markers have demonstrated pre-clinical AD changes in the MTL but it has not been possible to confirm that these MRI changes correspond to the presence of NFTs. Here, we present a method termed Projective LDDMM for aligning sparse measurement profiles of AD pathology (i.e., 2D digital histology images) with 3D MRI. We reconstruct measures of 2D NFT density in the dense metric of 3D MRI, using the Mai Paxinos Atlas coordinates for two cases of advanced AD. Analyses reveal the highest levels of NFT density in the rostral third (10-15 mm) of the hippocampus and the adjoining regions of the entorhinal cortex and amygdala. These findings emphasize the selective vulnerability of MTL subregions in AD, and suggest that high resolution MRI methods might benefit from focusing on the rostral MTL to more closely link these MRI images to AD neuropathology.

with 3D MRI. We reconstruct measures of 2D NFT density in the dense metric of 3D MRI, using the Mai Paxinos Atlas coordinates for two cases of advanced AD. Analyses reveal the highest levels of NFT density in the rostral third (10-15 mm) of the hippocampus and the adjoining regions of the entorhinal cortex and amygdala. These findings emphasize the selective vulnerability of MTL subregions in AD, and suggest that high resolution MRI methods might benefit from focusing on the rostral MTL to more closely link these MRI images to AD neuropathology.
Keywords: Projective LDDMM, Alzheimer's disease, Tau Tangles, Multimodal and Multiscale Image Registration Alzheimer's disease (AD) is the leading cause of dementia worldwide [1]. Diagnosis and characterization of AD in its early stages remain key challenges, as existing technologies limit the identification of the neuropathological patterns thought to emerge years before symptom onset [2,3,4]. In clinical practice, AD is typically first characterized by progressive clinical changes in memory and behavior, and subsequently through imaging changes that indirectly reflect AD neuropathology (i.e. misfolded proteins, tau and amyloid-Beta (Aβ)) [5,6,7]. Efforts to identify and understand the spatiotemporal profile of AD in its early stages have centered on these biomarkers [8]-measures that indirectly reflect the underlying pathology, which are obtainable over the course of disease. Of the methods used, neuroimaging has emerged as a prominent player with the ability to localize pathology non-invasively (e.g. tau/amyloid PET) [9,10], and with proposed surrogates such as shape diffeomorphometric markers (e.g. MRI) [11,12,13]. While these imaging measures have shown consistency with Braak staging [5,6], an accurate rendering of the 3D spatiotemporal profile of tau and Aβ at the micron scale has not been achieved [9,10]. The principle challenge has been integrating the 2D sparse measurements of histology, which are direct measures of disease, to the MRI 3D markers which are at much lower in-plane resolution.
This paper focuses on a new class of image-based diffeomorphometry methods which we term Projective LDDMM for aligning sparse 2D histological profiles to 3D coordinate systems across micron and millimeter scales. The 3D MRI to 2D digital histology mapping is representative of a class of multiscale, multi-modality mapping in biomedical research including traditional light microscopy mapping to dense reference atlases [14,15,16,17], light sheet methods [18,19], and spatial transcriptomics [20,21,22,23]. We formulate the dense mapping of atlases to sparse images problem using the random orbit model of computational anatomy [24,25,26,27] in which the space of dense anatomies I ∈ I is modelled as an orbit of a 3D template under the group of diffeomorphisms. Projective LDDMM models the sparse 2D histological observables not as an element of the orbit I but rather a random deformation in dense 3D coordinates composed with a measurement projection to sparse coordinates. While LDDMM provides the geodesic metric [28,29] on the orbit of 3D anatomies, there is no symmetry between the observable and the template, in general, and there should not be. This departs significantly from the symmetric methods [30,31].
Alignment specifically of the modes of histology to MRI warrants two extensions of the basic model of Projective LDDMM. First, cross-modality similarity modelling is essential. Several strategies for representing image similarity have emerged including cross-correlation [32], mutual information [33], and local textural characteristics [34]. Our approach is to extend previous work [35,36] by modelling a photometric transformation of histology to MRI with Mallat's Scattering Transform [37,38] to represent local radiomic textures at histological scales. Second, histological images carry large numbers of imperfections with tears, image stitching, and lighting variations. Extending previous work [15,36], we introduce Gaussian mixtures models in the image plane of each histological slice to interpret image locations as matching tissue, background, or artifact. We proceed by way of the Expectation-Maximization (EM) algorithm [39] in estimating deformations that prioritize image matching at locations that are, in turn, estimated more likely to be matching tissue.
Here, we use Projective LDDMM to reconstruct the 3D geometries of two sets of 2D histological sections taken from the medial temporal lobe (MTL) of advanced cases of AD. As tau has exhibited stronger predisposition over Aβ for segregating to particular brain regions (ERC, CA1, subiculum) and layers (superficial) of cortex in AD [5], we use machine-learning based methods to detect and quantify neurofibrillary tangles (NFTs) from histological images. Modeling this data in a measure theoretic framework amenable to quantifying trends at different scales [40], we transport these detections to the 3D space via the correspondences yielded by Projective LDDMM.

Projective LDDMM
In the random orbit model of computational anatomy [24], the hidden space of human anatomical images is modelled I : R 3 → R r as an orbit under diffeomorphisms of a template The observables J : R 3 → R q are modelled as a random field with mean due to the randomness of diffeomorphic deformation and measurement process. For different problems of interest, the atlas image is R r -valued with, for instance, r = 1 corresponding to single contrast MRI or r = 6 for diffusion tensor images (DTI) [41]. Likewise, observables are R q -valued with q = 3 for traditional histological stains or q >> 6 for alternative representations such as the scattering transform [37] encoding the meso-scale radiomic textures in histology images (see Section 2.2). In general, the range space of 3D templates versus targets do not have the same dimension, so q ̸ = r.
Projective LDDMM is characterized by the fact that the observable is not dense in the 3D metric of the brain. Rather, the observable(s) result from either optical or physical sectioning, as in histological slice preparation, taking LDDMM into the projective setting akin to classical tomography [42,43]. The sample measured observables J n (·), n = 1, 2, . . . are a series of projections P n of I(·) on the source space X ⊂ R 3 to measurement space Y ⊂ R 2 (or Z ⊂ R 1 ) defined through the class of point-spread functions associating source to target: with J n = P n I • φ −1 + noise .
We adopt measure theoretic notation, p n (y, dx) for describing point-spreads to accommodate those taking the form of generalized functions, such as the delta dirac. Density notation δ(x − x 0 )dx corresponds to the measure notation δ x0 (dx), each evaluated against a test function f (x) ∈ C 0 , yielding f (x 0 ).
The diffeomorphism φ is generated as the solution to the floẇ with velocity field v t , t ∈ [0, 1] controlling the flow constrained to be an element of a smooth reproducing kernel Hilbert space (RKHS) (V, ∥·∥ 2 V ) with the entire path square integrable 1 0 ∥v t ∥ 2 V dt < ∞ ensuring smoothness and existence of the inverse [44].
This gives us the first variational problem of Projective LDDMM.

Histological Sectioning: Crossing Modality and Projective Plane Distortions
Because of tissue section preparation, we attach two transformations on the target planes. Deformation of the tissue section geometry associates diffeomorphisms to the histology, ϕ n ∈ Φ : R 2 → R 2 both rigid and high-dimensional, expanding the dimensions first used in [14] for block sectioning. For crossing from the range space of histology contrast to MRI we introduce a complete expansion, reparameterizing the histology with 48 dimensions via Mallat's Scattering Transform [37,38], and solve for the optimal linear predictor f ωn , matching the scattering transform output to the scalar MRI. Via a subsampled Scattering Transform, color histology images J n : Y → R 3 , Y ⊂ R 2 are resampled at the resolution of MRI to 48-valued vector fields through alternating wavelet convolutions and nonlinear modulus operators across scales (see Appendix C): S : J n (·) →Ĵ n (·) = (Ĵ n,1 (·), . . . ,Ĵ n,48 (·)) T .
S : Jn →Ĵn(·), n = 1, . . . , N fω n :Ĵn → fω n (Ĵn)(·), n = 1, . . . , N For non-rigid modeling of ϕ n , a penalty term 1 0 ∥u n,t ∥ 2 U dt is added to (6), withφ n,t = u n,t • ϕ n,t , ϕ n,0 = Id. For estimating the cross-modality dimensions ω n , we have not mapped a large enough empirical sample to build an informative prior π ω on those high dimensions. Therefore, we treat them in a maximum-likelihood setting, optimizing them with initial conditions of deformation and image plane dimensions fixed, then solve the variational problem over all of the other dimensions with the ω n estimates fixed (see Section 2.2). This avoids collapse of the variational problem in these high dimensional settings.

Optical Sectioning, PET, and Parallel Beam Tomography: Ideal and Non-Ideal Planar and Linear Projections
Confocal optical sectioning reconstruct volumes X ⊂ R 3 with models that are fundamentally 3-D point-spreads p n (x, dx),X ⊂ R 3 [45,46], with imaging focused to n = 1, . . . , N measurement planes with significant blur out of plane. The mean field of the measurement volume J n (x),x ∈ R 3 are given by the projections P n I(x),x ∈ R 3 of (4). Two-dimensional (2D) positron emission tomography (PET) introduces point-spreads p n (y, dx) for reconstruction which are less idealized supported over planes Y ⊂ R 2 with uncertainty perpendicular to the line of flight but as well a second measurement the time-of-flight of the annihilating protons to the detectors [47]. Generally the point-spreads are modelled as cigar shaped two-dimensional Gaussians in the plane oriented by N angles θ n , n = 1, . . . , N , with high fidelity systems having N > 96, with standard deviation of uncertainty significantly larger along the lines of flight than perpendicular to them.

Biological Results: Multi-scale Maps of Tau Tangles in 3D
We present, in this section, the geometric and neuropathological results for corresponding pairs of histology images and MRI, taken postmortem, of MTL tissue from two cases of advanced AD. To acquire these data, several steps are required (see Section 2). Solutions to Variational Problem 2 (Equation 6) yield geometric reconstruction of histologically stained tissue in 3D. Figure  1 illustrates the two sets of 35 individual digitized sections on which NFTs were detected and geometric mappings to 3D were estimated. Subsequent coordination of both MRI and mapped histology with the Mai Paxinos Atlas is demonstrated in Figure 2, with coronal Mai views shown for an example intersecting histological slice.    Accuracy of tau tangle detections was evaluated both at intermediate and final steps of our detection algorithm. Estimates of accuracy in per pixel tau probabilities were computed using 10-fold cross validation on the entire training dataset for each brain sample. Table F3 (see Appendix F) shows accuracy metrics for one brain sample, with mean AUC and of 0.9860 and accuracy of 0.9729.
Accuracy of individual counts of NFTs as output following segmentation by the watershed algorithm (see Section 2.6) were estimated by comparison to manually annotated patches of ERC tissue reserved for validation. Between 5 and 20 mm 2 patches in the region of the ERC on 10 roughly consecutive slices of one brain were selected for annotation. Per pixel annotations were completed by a single individual over the course of 2-3 weeks. Each patch totaled approximately 2, 500, 000 pixels, yielding a total of 25 million for the 10 patches-on the order of the number of voxels for 20 whole brain MRIs at 1 mm resolution. Densities of NFTs (counts per cross-sectional tissue area) for each patch were computed using counts of manually annotated tangles vs. algorithmic output and rescaled to the range, [0, 1], for accurate comparison. Slice-by-slice relative densities are illustrated in Figure 4, exhibiting high levels of similarity between human and machine generated NFT densities. Counts of detected NFTs, cross-sectional tissue area, and MTL subregion (from MRI deformed to 2D) were computed in the space of histology slices. Average NFT densities per MTL subregion, tallied from all histology slices per brain sample, showed highest amounts of NFTs in amygdala, ERC, CA1, and subiculum for both advanced AD samples (see Figure 5,7).
Histological data was modeled as discrete particle measures and subsequently transported to the 3D space of the Mai Paxinos atlas via estimated transformations, ϕ n , φ (see Section 2.7). Diverse modes of resampling yielded distributions of NFT density within, over, and between MTL subregions (see Appendix G). Resampling via Gaussian kernels yielded smoothed NFT densities computed within the dense metric of the brain at approximate resolutions of MRI. Resampling via nearest neighbor kernels projected particle measures to the surface of segmented MTL subregions, such as amygdala, ERC, CA1, and subiculum for visualization in 3D. These are both shown in Figure 5. Average NFT densites per region are summarized in Figure 7. Tau pathology localized not just to particular regions (e.g. amygdala, ERC, CA1, subiculum), but within them. As illustrated in Figure 6, high densities of NFTs in amygdala and ERC concentrated particularly at the border between the two structures, with tau migrating to the inferior, medial boundary of the amygdala. NFT density also appeared to segregate within the hippocampus with highest densities achieved in the anterior third of the hippocampus in both brain samples (see Figure 7). Fig. 6 Posterior view of amygdala-ERC boundary in brain sample 1 (top) and brain sample 2 (bottom). NFT densities projected to and smoothed over surface of each structure independently. Outline of CA1 and subiculum surfaces shown in black mesh. First two rows illustrate NFT densities within 2 mm sliding window along Mai Z axis for amygdala and whole hippocampus (top) and amygdala and ERC (middle). Bottom row illustrates global average NFT density within each MTL subregion. All densities calibrated against simultaneously stained slices from both brains (see Section 2.6).

Algorithm for Solving Projective LDDMM with In Plane Transformation
To solve the Variational Problem 2 for 3D atlas, I, and set of 2D targets (J n , for n = 1, ..., N ) we formulate an algorithm that alternately optimizes for the deformation in 3D space and the geometric transformation in 2D space, while holding the other fixed. The algorithm can be implemented to incorporate increasing complexity as needed first for crossing modalities and second for crossing resolutions, as needed, for instance, to map 3D MRI to 2D histological slices, as is presented here. In its simplest form, I and J are of the same modality, yielding f ωn fixed as identity, and ϕ n is modeled as a rigid motion in plane, as in Lee et. al mapping histological sections to an atlas of the mouse brain [14]. Transformations φ, ϕ n for n = 1, · · · , N are estimated following Algorithm 1.

Return to A
When ϕ n 's are broadened to non-rigid diffeomorphisms, as in [15], each ϕ n is estimated in step B via a separate iteration of LDDMM for each target J n . Here, ϕ n 's encapsulate both rigid and non-rigid components. Separate gradient based methods are used to update each component in step B with velocity fields updated using Hilbert gradient descent as in [50] and linear transform parameters updated by Gauss-Newton [51].

Optimization Algorithm for Solving Projective LDDMM Crossing Modalities and Resolutions via Scattering Transform
Crossing modalities at similar resolution (e.g. 3D MRI and downsampled 2D histology slices) requires re-introduction of f ωn as a polynomial mapping between range spaces of template and target, giving a similar formulation to that used in Tward et. al [15]. In this work, we introduce the Scattering Transform [37] for crossing modalities at differing resolution. We model contrast variations between histology and MRI by (i) representing local radiomic textures via 48 dimensions of histological scales using Mallat's Scattering Transform [37,38], (ii) dimension reduction projecting the 48 dimensional scattering transform onto a 6-dimensional PCA basis, and (iii) solving for the optimal linear predictor f ωn matching the scattering transform 6-dimensions to the scalar MRI. Histology images J n : Y → R 3 , Y ⊂ R 2 are color images resampled to 48-valued vector fields via the Scattering Transform (see Appendix C): S : J n (·) → J s n (·) = (J s n,1 (·), . . . , J s n,48 (·)) T , defined at the resolution of MRI. A basis, b 1 , · · · , b 48 , is computed offline for the space of Scattering coefficients using PCA. A 6-valued vector image is generated from J s n as the projection onto the 6 largest eigenvalue basis elements, yielding the predictor, f ωn as affine including a constant offset: The mean field histology f ω parameterized by linear weights ω ∈ R 7 becomes: f ωn (Ĵ • ϕ n (·)) = ⟨Ĵ • ϕ n (·), ω n ⟩ R 7 , n = 1, . . . , N .
(7b) Figure 8 shows a mean field section using the scattering transform. These linear weights ω n are estimated from initialized ϕ n 's and φ following Algorithm 2.
Initializations of ϕ n and φ are estimated following the approach in Tward et. al [36] in which cubic polynomials are used to match MRI range space to histology range space. Solutions for ω n are then obtained using the pseudo-inverse: Algorithm 2 A: Solve for ω n 's: 1. Initialize φ 1 , ϕ n , n = 1, . . . , N .

Algorithm for Solving Projective LDDMM With Distortions and Missing data via EM Algorithm
The last element of complexity we introduce aims to account for the many distortions and subsampled missing data sections not wholly accounted for by geometric transformations, ϕ n . We formulate a weighted, least-squares problem by introducing weights representing different model interpretations of the pixels in each histology tissue section: measurements (e.g. tissue foreground), artifacts (e.g. tears and distortion), or background, following the example in Tward et. al [36]. Each model is characterized as a Gaussian with standard deviation, σ M , σ A , σ B and mean: with µ A and µ B representing artifacts and background. Weighted least-squares interprets the images weighing each model π j n,k (·), with 3 k=1 π k = 1, giving The weights arise from the E-step of an Expectation-Maximization (EM) algorithm [39], that we use for estimation of parameters θ n = ϕ n in step C of our Algorithm 2. Weights are a function of the previous parameters θ j n = ϕ j n hence giving the iteration: selecting at each point in the image the appropriate model for giving the spatial field of weights. This iteration corresponds to a Generalized EM (GEM) algorithm [39] (see Appendix D for proof). The results highlighted in Section 1.4 were generated following the approach of this section. For select slices in each brain sample, in plane ϕ n 's were estimated in step C of Algorithm 2 as spline transformations via manual landmark placement.

Specimen Preparation and Imaging
Brain tissue samples were prepared by the Johns Hopkins Brain Resource Center. The demographics and pathological staging of each sample analyzed are summarized in Table E, Appendix E. From each formalin immersion fixed brain, a portion of the MTL including entorhinal cortex, amygdala, and hippocampus, was excised in 3-4 contiguous blocks of tissue, sized 20-30 mm in height and width, and 15 mm rostral-caudal (see reconstructed MRI of tissue blocks, Figure 1).
Each block was imaged with an 11T MR scanner at 0.125 mm isotropic resolution and then cut into two or three sets of 10 micron thick sections, spaced 1 mm apart. Each block yielded between 7 and 15 sections per set. Sets of sections were stained with PHF-1 for tau tangle detection, 6E10 for Aβ plaque detection, or Nissl, and digitized at 2 micron resolution.

Segmentations of MTL Subregions
In all brains, MTL subregions were manually delineated using Seg3D [52]. Individual block MRIs were rigidly aligned using an in-house manual alignment tool, and per voxel labels were saved for the composite MRI for each brain. Delineations were deduced from patterns of intensity differences, combined with previously published MR segmentations [53,54,55] and expert knowledge on the anatomy of the MTL. The established borders were applied in three other brains, showing consistent results (in preparation, EX, CC, SM, DT, JT, Alesha Seifert, Tilak Ratnanather, MA, MW, and MM). In one brain sample, corresponding regional delineations were drawn on all histology sections stained with PHF-1 (see Figure 3). Delineations were based on visible anatomical markers and were afterwards confirmed with a corresponding Nissl-stained set of sections. In each of these sections, cytoarchitectonic borders between areas of interest were indicated, independently from the other datasets, using an previously published cytoarchitectonic accounts of the MTL [56,57,58,59,60,61]. Labels were assigned per pixel to 4x-downsampled histology images at a resolution of 32 microns and used to evaluate accuracy of registration (see Section 1.4). Regions of interest include amygdala, entorhinal cortex (ERC), cornu ammonis fields (CA1, CA2, CA3), and subiculum (see Figure 1).

Density Determination of NFTs
Patterns of tau pathology are summarized as total counts of tau tangles (NFTs) per mm 2 of cross-sectioned tissue. NFT counts were computed using a 2-step algorithm: (1) prediction of per pixel probabilities of tau and (2) segmentation of these probability maps into discrete NFTs.
As described previously [62,36], we used a convolutional neural network to model and predict probabilities of being part of a tau tangle for each pixel in a digital histology image. To capture larger contextual features as well as local information for producing per pixel probabilities at high resolutions, we trained UNETs [63] with the architecture described in Table F2 (see Appendix F). Given differences in staining intensity between brain samples, we trained separate UNETs, each with the same architecture, for each brain sample. Training data per brain was generated on every third slice of histology. Between 8 and 24 sample zones, sized 200-by-200 pixels were selected at random until 8 zones covered tissue (not background). Every pixel in each zone was manually annotated, 1 or 0, as part of a tau tangle or not.
Counts of NFTs in each histology slice were generated by segmenting the probability maps output from the trained UNET. Segmentations were computed using an opensource implementation of the watershed algorithm [64] to extract connected components with "high probability" of tau. Each component was defined as an individual NFT, with center, area, and roundness computed as features.
In the results above (Figures 4,5,6), NFT densities are reported on different scales for each brain sample as differences in staining intensity, timing, and handling of tissue samples for each brain resulted in different scales of absolute counts of NFTs. To compare NFT densities across brain samples, subsets of 4-5 additional sections from each brain were selected in the rostral hippocampus at approximate locations of original sections. These two subsets were stained simultaneously to achieve consistency between them. A UNET was trained on the original training data from both brains and NFTs were detected and summed across each of original and new sets of histology slices. Ratios of tau tangles detected on the original vs. new version of each section were computed, yielding average ratios of 4.4 and 1.7 for brain samples from subjects 1 and 2, respectively. Finally, assuming consistent overestimation in tau tangle density measures across sets of slices, measures along the rostral-caudal axis and in MTL subregions were rescaled in each brain according to the average ratio factor to produce comparable NFT densities (see Figure 7).

Particle Representation of Histological Data
We model histology data at the microscopic scale following the generalized measure approach in [40] where each particle of tissue carries a weighted Dirac measure over histology image space and a Dirac measure over the feature space w i δ yi ⊗δ fi , y i ∈ Y ⊂ R 2 and F = R 2ℓ . Weights reflect sampled tissue area captured in each particle measure, defined at the finest scale (µ 0 ) as cross-sectional area in the histology plane w i ∈ {2µm 2 , 0}, computed with thresholding using Otsu's method [65]. The first ℓ dimensions of f i ∈ F denote the number of tau tangles in each of ℓ MTL subregions. The second ℓ dimensions denote the fraction of sampled tissue area (w i ) within each of ℓ MTL subregions. At the finest scale, We transfer measures via diffeomorphisms ϕ n and φ and rigid transformation to the space of template I and the Mai Paxinos Atlas. Discrete weights w i adjust according to in plane expansion/contraction of cross-sectional tissue area, with adjustment at the fine scale by ϕ n given by the varifold action: To cross scales we use the fundamental decomposition of the particle measures with ρ being the density of the model and µ x the field of conditional probabilities on the features. Our tranformation across scales non-linearly rescales space and smooths the empirical feature distributions on the features. Spatial resampling is determined by π(x, x ′ ), the fraction that the particle at x assigns to the particle at x ′ , with R d π(x, x ′ )dx ′ = 1. The smoothing on the field of conditional probabilities occurs according to choice of kernel k((y, α y ), ·): with w 1 (y) = i∈I w i π(x i , y) α y = 1 Spatial resamplings at MRI resolutions (0.125 mm) and over surface boundaries of MTL subregions were achieved through isotropic Gaussian resampling and nearest neighbor resampling, respectively, through choice of π (see Appendix G). Feature reduction occurs via maps β → ψ(β) ∈ F ′ , for probability measures β with k((y, α y ), ·) = δ y ⊗ δ ψ(αy) (·) giving Here, ψ(·) reduces feature dimension by taking empirical distributions α y over each of 2ℓ dimensions to expected first moments for each corresponding dimension, giving F ′ = R 2ℓ with: Total NFT density is computed from the sum of the first ℓ features while NFT density per region is computed from the ratio of feature value j to ℓ + j for any of j = 1, · · · , ℓ MTL subregions.

Surface Smoothing using Laplace Beltrami Operator Basis
Spatial variations in NFT density within MTL subregions are visualized as smooth functions over the surface of each corresponding region. Particle mass belonging to a given subregion volume is "projected" to the surface boundary using a nearest neighbor kernel for π, as defined in Equation 12b (see Section 2.7): We construct functions, g τ (x) and g a (x), to represent the total number of NFTs and cross-sectional area of tissue from discrete particle measures of particles projected to the surface vertices x i ∈ V . We generate smooth representations of NFT density (ĝ τ (·) ga(·) ) using Laplace-Beltrami on each surface [66]: for both g τ (·) and g a (·), where B := {β 1 , · · · , β N } is a basis for the Laplace-Beltrami operator, and k the smoothing constant (see Appendix H).

Discussion
These findings demonstrate that 3D maps of tau tangle density at mm and micron resolution show a strong spatial predominance of tau in the rostral third of the MTL, with especially high densities measured in the amygdala. They can be considered a 'proof of concept' of an approach comprised of three key components. (1) The definition of Projective LDDMM is the key contribution of this work, as a class of image-based diffeomorphometry methods for reconstructing sparse and irregularly sampled data in the dense metric of the 3D brain. (2) The coupling of geometric and contrast transformations generalizes these methods further to images of different modalities. (3) The use of a measure theoretic framework to model data following reconstruction (i) preserves data mass and values and (ii) allows for resampling the dense metric to generate statistical distributions at arbitrary resolution. Evidence of sparsity in original data measurements is maintained in the preservation of data mass while reconstruction and resampling within the dense metric offers smooth interpolation between the data values. We demonstrate the success of this approach at reconstructing spatial profiles of tau tangle density in two MTL samples from individuals with advanced AD. By computing NFT "density" as counts of tau tangles per cross-sectional area, we maintain data values reflective of the original 2D histological space in which they were measured. Reconstruction within the dense 3D metric of the brain allows for resampling of these measures at micron and mm resolution along both regular and irregular sampling "grids", such as over surface boundaries of MTL structures. Reconstruction within the spatial coordinate system of a well-known atlas (e.g. Mai atlas) allows for consistent resampling across brain samples, such as along the rostral-caudal axis, for comparing profiles of pathology.
Both reconstructions illuminate AD as a spatially-oriented disease, which further motivates the need for generating such 3D spatial reconstructions to uncover these significant patterns. Specifically, two main patterns in tau pathology emerged for both brain samples: (1) high levels of NFT density in amygdala and ERC, and (2) a high-to-low gradient of NFT density along the rostral-caudal axis of the hippocampus (see Figure 7), localizing NFT density to the regions where it abuts the Amygdala and ERC. Though both samples reflect end-stage pathology, the rostral-caudal gradient is consistent with the spread of tau outlined in Braak staging [5,6], and supports recent work by Yushkevich et. al that found similar segregation of tau pathology in the anterior vs. posterior portions of hippocampal structures in cases of earlier stage disease [35].
The role of the amygdala in AD has historically been under less investigation than the MTL subregions initially highlighted by Braak (e.g. TEC, ERC, CA fields, and Subiculum) [5]. Nevertheless, a number of studies have suggested its involvement not just with AD but other neurodegenerative diseases involving misfolded proteins [67,68,69]. Its emergence here as a key reservoir of tau pathology complements these studies as well as Yushkevich et. al, who found high burdens of tau pathology in the amygdala [35]. Furthermore, previous work from our group in MRI diffeomorphometry suggested a similar effect of AD pathology in the amygdala as in the ERC based on shape markers reflective of atrophy [70,71]. Both basolateral and basomedial regions of the amygdala demonstrated significant atrophy, consistent with the appearance of higher tau densities, here, in the amygdalar regions adjacent to the ERC ( Figure 6).
Though this work has been limited to the analysis of two brain samples, the methods we've described and the distributions of tau pathology that have emerged suggest a number of avenues for expanding on these efforts. Projective LDDMM, as we've formulated it here, can easily see value in applications to light sheet microscopy, spatial transcriptomics, and other emerging imaging modalities in which sparsely sampled data sets must be reconstructed in the dense space of some atlas for full appreciation of the spatiotemporal resolution of these modalities. Within the realm of AD research, we foresee its use in drawing additional correspondances between microscopic neuropathology and imaging modalities, such as MRI, PET and DTI, in the validation and development of biomarkers.
The segregation of tau pathology to the rostral most third of the hippocampus, the amygdala, and the ERC also suggests direction of future study in AD. First, MRI images that routinely capture the entirety of the hippocampus might be refined to target this smaller section of the hippocampus, amygdala and ERC, offering greater specificity to AD and sensitivity to neuropathological changes, as higher resolutions could be achieved with a narrower region of capture. Second, the spatiotemporal distribution of neuropathology at the earliest stages of AD remains evasive and yet most significant in efforts to achieve earlier diagnosis. Emerging evidence of the neuropsychiatric syndromic complex known as Mild Behavioral Impairment (MBI) amongst individuals prior to the onset of AD [72,73] suggests a role for the amygdala, as an emotional and behavioral organ of the brain, particularly in these early stages. Hence, we are currently using our methods to reconstruct the distributions of tau pathology in earlier cases of AD to elucidate whether the amygdala shows similarly high levels of tau in these stages.
In summary, AD research, like so many other fields of neurodegeneration and neurobiology, harbors a disparity between knowledge, communication, and scientists working at the microscopic level and those at the macroscopic level.
This disparity manifests prominently in incomplete understanding of the 3D spatiotemporal progression of pathology (tau and Aβ) in AD. This, in turn, has prohibited the progression of MRI biomarkers for earlier diagnosis of AD, as they have not been adequately linked to corresponding 3D patterns of neuropathology. We have developed the method we call Projective LDDMM to address this gap.
Supplementary information. This manuscript is accompanied by a set of Supplementary Notes (Appendices) (referenced A-H in the main text).

Availability of Data and Materials
Due to size limitations, imaging data is available upon request.

Code Availability
Code used for training and applying UNET in tau tangle detection can be found here: https://github.com/twardlab/ADproject. Code for solving Projective-LDDMM, building measure theoretic data representations, and resampling across scales can be found here: https://github.com/kstouff4/ projective-lddmm.

Authors' Contributions
MA, SM, and JT recruited brain samples and oversaw preparation of histology and imaging. MW, CC, and EX manually segmented histology images and MRI. KS, MM, and DT developed methods. KS and DT implemented algorithm. KS and MM drafted manuscript. All authors contributed to editing final manuscript.

Appendix A Classical Tomography
The Radon Transform is used in classical tomography [48] to describe the generation of sinograms as projected images at different angles. The transform is typically written in functional notation as an integral along one dimension, with Lebesgue measure as: R θ I(z) = R I(z sin(θ) + z cos(θ), −z cos(θ) + z sin(θ))dz.
In the notation introduced in Section 1.1, we extend this to an integral over all of R 2 , using Dirac delta measures that assign nonzero measure only to lines in this same set: The line integral indexed by (θ, z) is modeled as a single projection P n I(z), z ∈ R indexed by a set of n = 1, ..., N and is given as in, equation 4, with point spread, p n (z, dy) = L θn (z) δȳ(dy)dȳ, as defined in section 1.2. We show below its equivalence to the line integral as defined in above (A1).

Appendix B Registration Accuracy
Similar to other groups [35], we evaluated accuracy of alignment between 2D histology and 3D MRI by looking at sets of discrete points (pixels) labeled in 2D versus corresponding voxels labeled in 3D and subsequently deformed to 2D (see Figure 3). Here, our sets of points were the sets of pixels (voxels) within a particular MTL subregion as delineated on 2D histology images and 3D MRI (see Section 2.5). We restricted our attention to particular subregions of interest (amygdala, ERC, CA1, subiculum), and measured accuracy by Dice Score and 95th Percentile Hausdorff distance for each region on each slice of one brain. Average overlap scores were 0.85, 0.82, 0.74, 0.65 for amygdala, ERC, CA1, and subiculum, respectively, while average 95th percentile Hausdorff distance was 1.886 mm, 1.039 mm, 1.972 mm, and 1.746 mm for amygdala, ERC, CA1, and subiculum, respectively.

Appendix C Scattering Transform
The Scattering Transform, by Mallat and Bruna [37,38], defines a cascade of alternating non-linear and non-commuting operators that map functions J n (·) ∈ L 2 (R d ) to representations, J ′ n (·). A scattering propogator takes signals down a path of alternating convolutions with wavelets (localized waveforms) and modulus operations. The path, p, is defined by a set of parameters λ ∈ 2 Z that scale a mother wavelet, ψ, so as to capture lower and lower frequencies.
We use 16 paths, p i ∈ P of length 1 or 2 and a high pass Gaussian filter, with width dilated according to λ i , in place of a traditional wavelet to achieve a representation both translation and rotation invariant in addition to Lipschitz continuous to small deformations. Each of the R,G,B channels of histology images are propagated independently along the same paths. Histology images are downsampled by a factor of 32 to reach the approximate resolution of MRI. Together, the subsampling and scattering of each channel yield a total of 48 scattering coefficients for each pixel in the downsampled histology image.

Appendix D EM Algorithm for Optimization
The iterative algorithm in step C (Algorithm 2) is based on the EM algorithm implying it is monotonic in the cost. The complete-data likelihood for each histology plane n = 1, 2, . . . N as a function of parameters, θ = ϕ n is where µ A and µ B represent artifact and background. The E-step takes the conditional expectation of the complete data log-likelihood with respect to the incomplete data, and the previous parameters θ old . The M -step generates our sequence of parameters: with π n,k (·) = E 1 k (L(·))|I • φ −1 1 (·, z n ), θ old ; .
The spatial field of weights π n,k (·) is the conditional expectation of the indicator 1 k (L(·)). We implement the Generalized EM (GEM) algorithm (see [39]) solving the maximization step: which generate a sequence with increasing log-likelihood:

Appendix F UNET Details
To identify individual tau tangles in histology images, we trained a UNET to predict per pixel probabilities of tau and subsequently segmented probability maps into discrete "tangles" using opencv's implementation of the watershed algorithm [64]. The architecture is summarized in table E. 10-fold cross validation was used to estimate the accuracy of UNET probability predictions prior to segmentation with the watershed algorithm. Table F3 demonstrates per-fold and average performance for the trained UNET for one of our brain samples.

Appendix G Resampling in Mai Atlas Space
To compute and compare NFT density measures across brain samples, we rigidly aligned all samples to the reference brain in the Mai Paxinos Atlas [77]. We used a manual alignment tool, created in-house, to select optimal alignments between surface renderings of the Hippocampus, Amygdala, and ERC of our brain samples and that of the Mai brain. Distributions of NFT density were computed in the coordinates of the Mai atlas, according to choice of π(x, x ′ ) governing physical spatial spread, and k((y, α y ), ·) governing smoothing over conditional feature distributions (Equation 12b). In all cases, total mass (2D cross-sectional tissue area of histological images) and total number of NFTs were conserved. To achieve this, initial NFT feature values (counts per MTL subregion) were reformulated following physical transformation (Equation 10b) as counts of NFTs per MTL subregion per weight of particle (i.e. total cross-sectional area): We highlight three different modes of resampling. Volumetric resampling (e.g. at mm resolution) was computed with a 3D isotropic Gaussian kernel with width, σ, and with new particles in a regular lattice, x ′ ∈ X ′ .
Resampling over 2D manifolds (e.g. the surface of ERC, Amygdala, CA1, or Subiculum) was computed using a nearest neighbor kernel, assigning all weight (tissue area) and NFTs from a particle at the fine scale to a single particle on the 2D manifold (e.g. vertex of a triangular mesh).
π(x, x ′ ) := 1 if x ′ = arg min X ′ ∥x − x ′ ∥ 2 0 otherwise (G7) Finally, resampling to a regular 1D lattice (e.g. the rostral-caudal axis of the human brain) was computed using an anisotropic Gaussian kernel to spread particle mass widely in two dimensions and narrowly in the third, with dimensions treated independently.
π(x, x ′ ) : In each case, feature reduction occurred via computation of expected first moments, as described in section 2.7.