A multivariate nonlinear mixed effects model for longitudinal image analysis: Application to amyloid imaging
Murat Bilgel, Jerry L. Prince, Dean F. Wong, Susan M. Resnick, Bruno M. Jedynak
AA multivariate nonlinear mixed effects model forlongitudinal image analysis: Application to amyloidimaging (cid:73) , (cid:73)(cid:73) Murat Bilgel a,b,c, ∗ , Jerry L. Prince a,b,d,e , Dean F. Wong e , Susan M. Resnick c ,Bruno M. Jedynak f a Image Analysis and Communications Laboratory, Johns Hopkins University School ofEngineering, Baltimore, MD, USA b Dept. of Biomedical Engineering, Johns Hopkins University School of Medicine, Baltimore,MD, USA c Laboratory of Behavioral Neuroscience, National Institute on Aging, National Institutes ofHealth, Baltimore, MD, USA d Dept. of Electrical and Computer Engineering, Johns Hopkins University School ofEngineering, Baltimore, MD, USA e Dept. of Radiology, Johns Hopkins University School of Medicine, Baltimore, MD, USA f Dept. of Mathematics and Statistics, Portland State University, Portland, OR, USA
Abstract
It is important to characterize the temporal trajectories of disease-relatedbiomarkers in order to monitor progression and identify potential points ofintervention. These are especially important for neurodegenerative diseases, astherapeutic intervention is most likely to be effective in the preclinical diseasestages prior to significant neuronal damage. Neuroimaging allows for the measure-ment of structural, functional, and metabolic integrity of the brain at the level ofvoxels, whose volumes are on the order of mm . These voxelwise measurementsprovide a rich collection of disease indicators. Longitudinal neuroimaging studiesenable the analysis of changes in these voxelwise measures. However, commonlyused longitudinal analysis approaches, such as linear mixed effects models, donot account for the fact that individuals enter a study at various disease stages (cid:73) DOI: 10.1016/j.neuroimage.2016.04.001 (cid:73)(cid:73)
This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/ . ∗ Corresponding author at: National Institute on Aging, Laboratory of Behavioral Neu-roscience, 251 Bayview Blvd., Suite 100, Rm 04B316, Baltimore, MD 21224, USA. Phone:+1-410-558-8151. Fax: +1-410-558-8674.
Email address: [email protected] (Murat Bilgel)
Preprint submitted to NeuroImage October 15, 2018 a r X i v : . [ s t a t . M E ] A p r nd progress at different rates, and generally consider each voxelwise measureindependently. We propose a multivariate nonlinear mixed effects model for esti-mating the trajectories of voxelwise neuroimaging biomarkers from longitudinaldata that accounts for such differences across individuals. The method involvesthe prediction of a progression score for each visit based on a collective analysisof voxelwise biomarker data within an expectation-maximization framework thatefficiently handles large amounts of measurements and variable number of visitsper individual, and accounts for spatial correlations among voxels. This scoreallows individuals with similar progressions to be aligned and analyzed together,which enables the construction of a trajectory of brain changes as a function ofan underlying progression or disease stage. We apply our method to studyingcortical β -amyloid deposition, a hallmark of preclinical Alzheimer’s disease, asmeasured using positron emission tomography. Results on 104 individuals witha total of 300 visits suggest that precuneus is the earliest cortical region toaccumulate amyloid, closely followed by the cingulate and frontal cortices, thenby the lateral parietal cortex. The extracted progression scores reveal a patternsimilar to mean cortical distribution volume ratio (DVR), an index of globalbrain amyloid levels. The proposed method can be applied to other types oflongitudinal imaging data, including metabolism, blood flow, tau, and structuralimaging-derived measures, to extract individualized summary scores indicatingdisease progression and to provide voxelwise trajectories that can be comparedbetween brain regions. Keywords:
Longitudinal image analysis, progression score, amyloid imaging
1. Introduction
It is important to characterize the temporal trajectories of disease-relatedbiomarkers in order to monitor progression and to identify potential points of in-tervention. Such a characterization is especially important for neurodegenerativediseases, as therapeutic intervention is most likely to be effective in the preclinical disease stages prior to significant neuronal damage. For example, in Alzheimer’s2isease, brain changes evident in structural, functional, and metabolic imagingmay occur more than a decade before the onset of cognitive symptoms (Batemanet al., 2012), with cortical amyloid- β (A β ) accumulation being one of the earliestchanges (Jack et al., 2013; Sperling et al., 2014a; Villemagne et al., 2013). Such brain changes can be measured using neuroimaging techniques and can be trackedover time at the individual level via longitudinal studies.Given the focus on preventing and delaying the onset of incurable neurodegen-erative diseases, the emphasis of clinical trials has shifted to studying clinicallynormal individuals with positive biomarkers, for example those exhibiting brain amyloid in the case of AD, in order to identify early intervention opportunitiesin the preclinical stages of disease (Sperling et al., 2014b). It is important to de-termine the temporal trajectories of hypothesized biomarkers in the early diseasestages in order to better understand their associations with disease progression.Current neuroimaging methods allow for the characterization of the brain at the mm level, generating hundreds of thousands of measurements that can beused as potential biomarkers of neurodegenerative diseases. Understanding thetemporal trajectories of these voxelwise measurements can provide clues intodisease mechanisms by identifying the earliest and fastest changing brain regions.Changes in voxelwise neuroimaging measurements over time are commonly studied using linear mixed effects models (Bernal-Rusiel et al., 2012, 2013;Ziegler et al., 2015). Univariate linear mixed effects models use time or age tocharacterize changes in a single imaging measure. However, time or age may notbe the appropriate metric for measuring disease progression due to variabilityacross individuals. While covariates can be included in linear mixed effects models to account for this variability, choosing the correct set of covariatesis difficult and covariates generally have a more complicated association withdisease progression than the assumed linear relationship of linear mixed effectsmodels. Instead, this variability can be accounted for by aligning individualsin time based on their longitudinal biomarker profiles within a multivariate framework. This is the premise of the Disease Progression Score method, whichhas been applied to studying changes in cognitive and biological markers related3o Alzheimer’s disease (Jedynak et al., 2012, 2014; Bilgel et al., 2014). It isassumed that there is an underlying progression score (PS) for each subject visitthat is an affine transform of the subject’s age, and given this PS, it is possible to place biomarker measurements across a group of subjects onto a commontimeline. The affine transformation of age removes across-subject variability inbaseline biomarker measures as well as in their rates of longitudinal progression.Each biomarker is associated with a parametric trajectory as a function of PS,whose parameters are estimated along with the PS for each subject. This allows one to “stitch” data across subjects to obtain temporal biomarker trajectoriesthat fit an underlying model (Fig. 1).Previous approaches have used certain cognitive measures, such as ADAS-Cog (Caroli and Frisoni, 2010; Yang et al., 2011), MMSE (Doody et al., 2010) orCDR-SB (Delor et al., 2013) as a surrogate for disease progression to delineate the trajectories of other AD-related cognitive measurements. These methodsoperate with the assumption that disease progression is reflected by a singlecognitive measurement rather than a profile of multiple measurements, andtherefore are inherently limited in their characterization of disease evolution.Younes et al. (2014) fitted a piecewise linear model to longitudinal data assuming that each biomarker becomes abnormal a certain number of years before clinicaldiagnosis, and this duration was estimated for each biomarker to yield longitudi-nal trajectories as a function of time to diagnosis. A quantile regression approachwas employed by Schmidt-Richberg et al. (2015) to align a sample of cognitivelynormals and mild cognitively impaired (MCI) with a sample of MCI and AD, and then to estimate biomarker trajectories. These approaches assume that allindividuals are on a path to disease and require knowledge of clinical diagnosis.Therefore, they are not suitable for studying the earliest changes in individualswho have not converted to a clinical diagnosis. Donohue et al. (2014) applied aself-modeling regression model within a multivariate framework to characterize the longitudinal trajectories of a set of cognitive, CSF, and neuroimaging-basedbiomarkers. This approach allows for across-subject variability only in the ageof onset, not in progression speed. Models incorporating fixed effects as well4s individual-level random effects have been proposed to study ADAS-Cog (Itoet al., 2011; Schiratti et al., 2015b) and regional cortical atrophy (Schiratti et al., changes in longitudinal biomarker measures as well as the appropriate thresholdsfor separating normal from abnormal measures (Fonteijn et al., 2012; Younget al., 2014). These methods characterize longitudinal biomarker trajectoriesin a discrete framework rather than a continuous one. Schiratti et al. (2015a)proposed an extension to their earlier approach to model multiple measures together. Biomarker trajectories are assumed to be identical except for a shiftalong the disease timeline, and this assumption prevents hypothesis testing re-garding rate of change across biomarkers. Furthermore, biomarkers are assumedto be conditionally independent given the subject-level random effects, butthis assumption is not realistic when biomarkers are voxel-based neuroimaging measurements.Here, we adapt the disease progression score principle to studying longitudinalneuroimaging data by making substantial innovations to the progression scoremodel and parameter estimation procedure. First, voxelwise imaging measuresconstitute the biomarkers in the model, and are analyzed together in a multi- variate framework. Studying progression at the voxel level rather than usingregion of interest (ROI)-based measures allows for the discovery of patterns thatmay not be confined within any given ROI. Second, since voxelwise imagingmeasures have an underlying spatial correlation, we incorporate the modelingof the spatial correlations among the biomarker error terms. Modeling spatial correlations makes the inference of the subject-specific progression scores lesssusceptible to the inherent correlations among the voxels. Third, we incorporatea bivariate normal prior on the subject-specific variables that define the relation-ship between age and PS. The prior allows a better modeling of the variance5ithin individuals and enables the incorporation of individuals with a single visit into the model fitting procedure. Fourth, instead of using an alternatingleast-squares approach for parameter estimation as presented by Jedynak et al.(2012), we formulate the model fitting as an expectation-maximization (EM)algorithm, which guarantees convergence to a local maximum and allows for anefficient model fitting framework for a large number of biomarkers. Finally, we present a statistical framework for comparing the onset and rate of progressionacross different regions. This paper extends our previous approach for analyzinglongitudinal voxelwise imaging measures using the progression score frameworkby incorporating a prior on the subject-specific variables, presenting a hypothe-sis testing framework for determining biomarker ordering, and performing an extensive validation of the method (Bilgel et al., 2015c).We first show using simulated data that the model parameters are estimatedaccurately and that modeling spatial correlations improves parameter estimation.We then apply the method to distribution volume ratio (DVR) images derivedfrom Pittsburgh compound B (PiB) PET imaging, which show the distribution of cerebral fibrillar amyloid. Models fitted using data for 104 participants with atotal of 300 PiB-PET visits reveal that the precuneus and frontal cortex show thegreatest longitudinal increases in fibrillar amyloid, with smaller increases in lateraltemporal and temporoparietal regions, and minimal increases in the occipitalcortex and the sensorimotor strip. Our results suggest that the precuneus is the earliest cortical region to accumulate amyloid. The results are consistentacross the two hemispheres, and the estimated PS agrees with a widely usedPET-based global index of brain amyloid known as mean cortical DVR. Thepresented method can be applied to other types of longitudinal imaging data tounderstand voxelwise trajectories and to quantify each individual scan against the estimated progression pattern. 6 . Method Our goal is to characterize the progression of disease or an underlying processas measured using a collection of relevant biomarkers. Disease or process stage, as indicated by a progression score (PS), s , is intrinsically related to time t ,measured as the age of a subject. Since individuals differ in their onset andrate of progression, the relationship between s and t varies across individuals.We model the progression s as a linear function of time t for each individualand allow for the prediction of separate slopes and intercepts to account for this variability across individuals.Generally, there is a particular presentation of symptoms and biomarkermeasurements at a given progression stage. Furthermore, as the disease orprocess progresses, there is a particular temporal progression of the biomarkers.In this work, we consider voxelwise PET measures as biomarkers and model the temporal trajectory of each biomarker, or voxel, as a linear function of theprogression s . Considering voxelwise neuroimaging measures as biomarkers mayappear unusual; however, these measures fit the NIH definition of a biomarker:“a characteristic that is objectively measured and evaluated as an indicator ofnormal biological processes, pathogenic processes, or pharmacologic responses to a therapeutic intervention” (Biomarkers Definitions Working Group, 2001). Asillustrated in Figure 1, PS aligns longitudinal biomarker measures better thanage since it accounts for differences across individuals in rates as well as baselinelevels of progression. After this alignment in time, the estimated biomarkertrajectories can be compared on the common PS scale. In the following subsections, we describe the progression score model in detail.7 igure 1: Illustration of the biomarker alignment concept in the progression score model. Thebiomarkers we consider in this work are PET measures of cerebral amyloid across a total of K ≈ ,
000 voxels.
Top:
Progression score (PS) aligns longitudinal measures better thanage, and allows for the estimation of a trajectory for each biomarker/voxel (in gray).
Bottom:
Estimated biomarker trajectories can be compared on the common PS scale.
The progression score s ij for subject i at visit j is assumed to be an affinetransformation of the subject’s age t ij : s ij = α i t ij + β i = q Tij u i , (1)where q ij = t ij , and u i = α i β i . The subject-specific variables, α i and β i contained in the vector u i , are assumed to follow a bivariate normal distribution,i.e., u i ∼ N ( m , V ), which are independent and identically distributed acrosssubjects. This model accounts for differences between subjects in the rate ofprogression via α , and in the baseline levels of disease progression via β . The prior covariance V is modeled as a 2 × V , given by ν , ensures that V ≡ V ( ν ) is positive8efinite (Pinheiro and Bates, 1996). Let U = U U U be an upper triangularmatrix such that V = U T U . If the diagonal elements of U are constrained to bepositive, then this Cholesky decomposition is unique. To ensure that the diagonal elements of U are positive in an unconstrained optimization framework, we usethe natural logarithm of the diagonal elements of U as parameters. We thenvectorize the upper triangular elements (including the diagonal) of U to obtainthe parameter vector ν = [log U U log U ] T , which uniquely parameterizes V and ensures its positive definiteness. The collection of K biomarker measurements form the K × y ij forsubject i at visit j . Longitudinal trajectories associated with these biomarkersare assumed to be linear and parameterized by K × a and b : y ij = a s ij + b + (cid:15) ij . (2)Here, a = [ a , a , . . . , a K ] T , b = [ b , b , . . . , b K ] T , and (cid:15) ij ∼ N K (0 , R ) is theobservation noise. (cid:15) ij are assumed to be independent and identically distributedacross subjects and visits. The matrix R is assumed to have the form R = Λ C Λ, where Λ is a diagonalmatrix with positive diagonal elements λ and C is a correlation matrix param-eterized by ρ . This parameterization guarantees that R is a positive definitematrix (Galecki and Burzykowski, 2013). For ease of notation, we let Λ ≡ Λ( λ ), C ≡ C ( ρ ), and R ≡ R ( λ , ρ ). When the biomarkers under consideration have a spatial organization, i.e.,if they are voxelwise measurements from medical images, then the correlationmatrix C can be described as a function of the spatial distance d ≡ d ( k, k (cid:48) )between pairs of voxels indexed by k and k (cid:48) as well as the spatial correlationparameter ρ . Possible univariate parameterizations (i.e., ρ = ρ ∈ R ) of C are C isa valid correlation matrix (Galecki and Burzykowski, 2013). Table 1: Spatial correlation functions.
Exponential C kk (cid:48) = e − d/ρ Gaussian C kk (cid:48) = e − ( d/ρ ) Rational quadratic C kk (cid:48) = d/ρ ) Spherical C kk (cid:48) = (cid:18) − dρ + (cid:16) dρ (cid:17) (cid:19) ( d < ρ ) d ≡ d ( k, k (cid:48) ) is the spatial distance between voxels indexed by k and k (cid:48) . The overall model, diagrammatically summarized using plate notation inFig. 2, is described by the following equations: s ij = q Tij u i (3) y ij = a s ij + b + (cid:15) ij (4) u i ∼ N ( m , V ( ν )) (5) (cid:15) ij ∼ N K (0 , R ( λ , ρ )) . (6)While this model is a mixed effects model since it incorporates the fixed effects a , b as well as the individual-level random effects u i , and is nonlinear in the pa-rameters, it departs from the form of the nonlinear mixed effects model describedby Lindstrom and Bates (1990). Therefore, instead of pursuing a restricted max-imum likelihood approach, we use an expectation-maximization (EM) approach, as described below.Let θ be the collection of model parameters m , ν , a , b , λ , ρ . The complete10 igure 2: Probabilistic model of the progression score using plate notation. The progressionscore s ij establishes the link between age t ij and the voxelwise observations y ij . Circlesindicate variables, and known measures are shaded. Rectangles indicate that the same modelapplies across visits (inner rectangle) and individuals (outer rectangle). Arrows indicatedependencies between variables and parameters. log-likelihood for the model is (cid:96) ( y , u ; θ ) = (cid:88) i (cid:96) ( y i , u i ; θ ) (7)= (cid:88) i log f ( y i | u i ; θ ) + log f ( u i ; θ ) (8)= − (cid:88) i,j log | πR | − (cid:88) i,j ( y ij − Z ij u i − b ) T R − ( y ij − Z ij u i − b ) − (cid:88) i log | πV | − (cid:88) i ( u i − m ) T V − ( u i − m ) , (9)where Z ij = aq Tij . y i is a vector of all biomarker measures stacked across allvisits of individual i , and y ij is a vector of all biomarker measures for individual i at visit j . Given n individuals, we use the summation notation (cid:80) i,j as ashorthand for (cid:80) ni =1 (cid:80) v i j =1 , where v i is the number of visits for individual i .Since the marginal likelihood f ( y ; θ ) involves an integral over all possiblevalues of u i , maximizing it directly is difficult. Therefore, we use the expectation-maximization (EM) approach, where we consider the subject-specific variables u i } as hidden. EM allows us to reformulate the maximization problem in termsof the complete log-likelihood (cid:96) ( y , u ; θ ). The observations y include biomarkermeasurements { y ij } at each visit. The unknown parameters are the subject-specific variable distribution parameters m and ν , the trajectory parameters a and b , and the noise covariance parameters λ and ρ . ( y , u ) are the complete data. Let θ (cid:48) = { m (cid:48) , ν (cid:48) , a (cid:48) , b (cid:48) , λ (cid:48) , ρ (cid:48) } be the previousparameter estimates. By Proposition A.1, the E-step integral is proportional to (cid:80) i (cid:82) Φ(˜ u i ; ˆ u (cid:48) i , Σ (cid:48) i ) (cid:96) ( y i , ˜ u i ; θ ) d ˜ u i , where Φ is the multivariate normal probabilitydensity function with meanˆ u (cid:48) i = (cid:88) j Z (cid:48) Tij R (cid:48)− Z (cid:48) ij + V (cid:48)− − (cid:88) j Z (cid:48) Tij R (cid:48)− ( y ij − b (cid:48) ) + V (cid:48)− m (cid:48) , (10)and covariance Σ (cid:48) i = (cid:16)(cid:80) j Z (cid:48) Tij R (cid:48)− Z (cid:48) ij + V (cid:48)− (cid:17) − . Note that for individuals witha single visit, (cid:80) j Z (cid:48) Tij R (cid:48)− Z (cid:48) ij is a singular matrix. Considering u i as parameters,as in Jedynak et al. (2012) or Bilgel et al. (2015c), is equivalent to assumingthat V has infinitely large diagonal elements (i.e., an uninformative uniform prior) such that its inverse disappears. Therefore, it is not possible to computeˆ u (cid:48) i for individuals having only one visit using this approach. On the other hand,incorporation of a bivariate normal prior on the subject-specific variables u i allows Eq. 10 to be computed for individuals with a single visit.Evaluation of the E-step integral involves second moments of a Gaussian random variable. We ignore the terms that do not depend on θ as they will not12e relevant in the maximization step and obtain: Q (cid:0) θ , θ (cid:48) (cid:1) = − (cid:88) i,j log | R |− (cid:88) i,j ( y ij − Z ij ˆ u (cid:48) i − b ) T R − ( y ij − Z ij ˆ u (cid:48) i − b ) − (cid:88) i,j Tr (cid:0) Z Tij R − Z ij Σ (cid:48) i (cid:1) − (cid:88) i log | V |− (cid:88) i (ˆ u (cid:48) i − m ) T V − (ˆ u (cid:48) i − m ) − (cid:88) i Tr (cid:0) V − Σ (cid:48) i (cid:1) . (11) Here, we provide the update equations for the EM algorithm obtained bymaximizing Q ( θ , θ (cid:48) ) with respect to each parameter, and provide derivations for these update equations in the Appendix. The update equations depend onprevious parameter estimates θ (cid:48) = { m (cid:48) , ν (cid:48) , a (cid:48) , b (cid:48) , λ (cid:48) , ρ (cid:48) } as well as the progressionscore estimates s (cid:48) ij = q Tij ˆ u (cid:48) i , where ˆ u (cid:48) i is as given in Eq. 10: a = ( (cid:80) i v i ) (cid:16)(cid:80) i,j y ij s (cid:48) ij (cid:17) − (cid:16)(cid:80) i,j y ij (cid:17) (cid:16)(cid:80) i,j s (cid:48) ij (cid:17) ( (cid:80) i v i ) (cid:16)(cid:80) i,j q Tij Σ (cid:48) i q ij + s (cid:48) ij (cid:17) − (cid:16)(cid:80) i,j s (cid:48) ij (cid:17) , (12) b = (cid:16)(cid:80) i,j y ij (cid:17) (cid:16)(cid:80) i,j q Tij Σ (cid:48) i q ij + s (cid:48) ij (cid:17) − (cid:16)(cid:80) i,j y ij s (cid:48) ij (cid:17) (cid:16)(cid:80) i,j s (cid:48) ij (cid:17) ( (cid:80) i v i ) (cid:16)(cid:80) i,j q Tij Σ (cid:48) i q ij + s (cid:48) ij (cid:17) − (cid:16)(cid:80) i,j s (cid:48) ij (cid:17) , (13) m = 1 n (cid:88) i ˆ u (cid:48) i , (14) ν = arg max ν Q ( θ , θ (cid:48) ) , (15) λ , ρ = arg max λ , ρ Q ( θ , θ (cid:48) ) . (16)Note that if C is fixed to be the identity matrix, a closed form solution for λ exists, as given in Equation A.15. Once the optimal parameters are found, the subject-specific variables are predicted using Eq. (10). As described in Proposition A.2, there are certain reparameterizations thatyield identical models. For example, one can multiply the trajectory slope13arameters by 2 and divide all progression scores by 2 (which is achieved by dividing all α and β values by 2) without altering the model. This is the scalingdegree of freedom. There is also a translation degree of freedom. We account forthese degrees of freedom and anchor the model by calibrating the progression scorescale. We calibrate such that baseline PS has a mean of 0 and a standard deviationof 1. This involves replacing the model parameters { m , V ( ν ) , a , b , R ( λ , ρ ) } with w m + z , w V ( ν ) , w a , b − zw a , R ( λ , ρ ) , where − zw = n (cid:80) i s i is themean progression score at baseline, and w = (cid:113) n (cid:80) i (cid:0) s i − n (cid:80) i s i (cid:1) is thestandard deviation of baseline progression scores. We standardize the subject-specific estimates α i , β i , and s ij accordingly: α ∗ i = wα i , β ∗ i = wβ i + z , and s ∗ ij = ws ij + z . This reparametrization yields a PS scale where 0 corresponds to the sample average at baseline and the variance of PS at baseline is 1. We first fit the model assuming that C = I K × K and Λ is a diagonal matrixwith positive elements λ . We denote the estimated Λ in this model as ˆΛ. Theestimates obtained from this model for the parameters a , b , m , V are used as initializations in the model where C = C ( ρ ). In this model where correlations aretaken into account, we assume that Λ = λ ˆΛ, where λ is an unknown parameterto be estimated. The spatial correlation function among those presented inTable 1 that results in the highest log-likelihood value is chosen for the finalmodel. The EM algorithm is implemented in MATLAB 8.1 and Statistics Toolbox 8.2(The MathWorks Inc., Natick, MA). Our code is freely available online. We use bootstrapping via Monte Carlo resampling to estimate confidenceintervals for each model parameter. We sample with replacement from the ρ at its estimate on the whole sample to enable faster computation. We compared our model to a linear mixed effects (LME) model that includedrandom intercepts and slopes at each voxel. The LME model for subject i , visit j and voxel k is given by y ijk = ( η k + η ik ) t ij + ( γ k + γ ik ) + ε ijk , (17)where η k γ k are the fixed effects, η ik γ ik ∼ N (0 , Ξ k ) are the random effects and ε ijk ∼ N (0 , σ k ) is the observation noise. We used the LME implementation inMATLAB Statistics Toolbox 8.2. We simulated visits such that the sample was similar to our PET data interms of number of visits per subjects and age range. We generated a data setwith 100 individuals, each with up to 7 visits with 5 × × θ = { m , ν , a , b , λ , ρ } at values close to those we observed in exploratory modelsfitted to DVR data. We generated u i from a bivariate normal distributionwith mean m and variance V ( ν ). The progression score for each visit was thencomputed as s ij = q Tij u , and observations were generated using the PS model.We performed 1000 bootstrap iterations to obtain 95% confidence intervals foreach parameter and subject-specific variable. We computed the cosine similarityfor each variable for each bootstrap experiment using φ T (cid:98) φ (cid:107) φ (cid:107) (cid:107) (cid:98) φ (cid:107) , (18)15here (cid:98) φ is the estimate of the variable of interest ( a , b , α , β , or s ) and φ is thecorresponding ground truth value. A value of 1 indicates a perfect estimate. We used longitudinal positron emission tomography (PET) data for partic- ipants from the Baltimore Longitudinal Study of Aging (Shock et al., 1984)neuroimaging substudy (Resnick et al., 2000). Participant demographics arepresented in Table 2. Starting in 2005, amyloid-beta (A β ) PET scans were ac-quired on a GE Advance scanner over 70 minutes following an intravenous bolusinjection of Pittsburgh compound B (PiB), yielding dynamic PET scans with 33 time frames. PET images were reconstructed using filtered backprojection witha ramp filter, yielding a spatial resolution of approximately 4.5 mm full width athalf max at the center of the field of view (image matrix = 128 × × the first two minutes to remove motion (Jenkinson et al., 2002). For registra-tion purposes, we obtained static images by averaging the first 20 minutes ofthe dynamic PiB-PET scan. Follow-up PiB-PET scans were rigidly registeredonto the baseline PiB-PET within each participant using the 20-minute averageimages. Baseline magnetic resonance images (MRIs) were rigidly registered onto their corresponding 20-minute PiB-PET average, and their FreeSurfer segmen-tations (Dale et al., 1999; Desikan et al., 2006) were transformed accordingly.Distribution volume ratio (DVR) images were calculated in the native spaceof each PiB-PET image using the simplified reference tissue model with thecerebellar gray matter as reference tissue (Zhou et al., 2003). The MRIs coreg- istered with the PET were deformably registered (Avants et al., 2008) onto astudy-specific template (Avants et al., 2010; Bilgel et al., 2015b) and transformedto 4 mm isotropic MNI space using a pre-calculated affine transformation. Theresulting mappings were applied to the DVR images that have been registeredto baseline to bring them into the MNI space. We used all voxels within the brain mask in the MNI space to fit the PS model.16 able 2: Participant demographics. MCI = mild cognitive impairment, SD = standarddeviation. Characteristic N = 104Baseline age in years, mean (SD) 77.0 (7.9)Range 55.7–93.4Female, n (%) 48 (46%)PiB-PET scans, n 300PiB-PET per subject, n 2.9 (1.8)Range 1–7Years between first and last scan, mean (SD) 3.3 (2.9)Range 0.0–9.0Mean cortical DVR at last visit (SD) 1.12 (0.19)PiB+ at last visit, n (%) 42 (40%)MCI diagnosis at last visit, n (%) 3 (3%)Dementia diagnosis at last visit, n (%) 4 (4%)Mean cortical DVR is a widely used measure obtained from PiB-PET imagesfor quantifying the level of brain amyloid. We computed mean cortical DVR foreach PiB-PET image by averaging the voxelwise DVR values across cingulate,frontal, parietal, lateral temporal, and lateral occipital cortices, excluding the sensorimotor strip. We used a mean cortical DVR threshold of 1.06, which wascomputed from a two-class Gaussian mixture model fitted on baseline measures,to separate individuals into PiB- and PiB+ groups, as described in Bilgel et al.(2015a). The confidence intervals obtained via bootstrapping allow for hypothesistesting. For the amyloid images, we focus on studying the precuneus, sinceprevious cross-sectional amyloid imaging studies have suggested that precuneushas the highest deposition levels (Mintun et al., 2006) and provided preliminary17vidence that it may be the most rapid accumulator (Rodrigue et al., 2012) among cortical regions. Our specific hypotheses are as follows:1. The precuneus has the highest amyloid load along stages of amyloidaccumulation.2. The precuneus accumulates amyloid faster than other cortical regions.We refer to the progression scores calculated using the DVR images as A β -PS.To test the first hypothesis, we compare the amyloid levels in the precuneuswith other regions at various A β -PS values spanning the range observed in ourdata set. Our purpose in performing these comparisons is to shed light onto thetemporal ordering of changes in different cortical regions. For each bootstrapexperiment, we average the predicted DVR levels ˆ y k ( s ) = a k s + b k within eachROI to obtain ˆ y r ( s ) = 1 | ROI r | (cid:88) k ∈ ROI r ˆ y k ( s ) , (19)where k is the voxel index and r is the ROI index. We then compute the followingstatistic for each bootstrap: T = ˆ y precuneus ( s ) − max r (cid:54) =precuneus ˆ y r ( s ) . (20)We reject the null hypothesis that the amyloid load in the precuneus at A β - PS = s is not different from other regions at significance level γ if the two-sided100(1 − γ ) confidence interval of the test statistic T , computed using the bootstrapestimates of T , does not contain 0. The smallest value of γ such that the two-sided 100(1 − γ ) confidence interval of the test statistic T contains 0 is the p -valueof the test. To test the second hypothesis, we pursue a similar approach using thetrajectory slope parameter a . For each bootstrap experiment, we obtain theaverage a k per cortical ROI to obtain a r , where r the ROI index. We thencompute the following statistic for each bootstrap: T = a precuneus − max r (cid:54) =precuneus a r . (21)18he condition for the rejection of the null hypothesis that the rate of amyloidaccumulation in the precuneus is not different from that of other regions is asdescribed previously.
3. Results
The Akaike information criterion (AIC) was 39 . × for the LME model,32 . × for the PS model where C = I K × K , and − . × for the PS modelwhere C = C ( ρ ), indicating that the PS model where spatial correlations aremodeled fits the data the best. We computed the percentage of variables thatare “correct” by counting the variables whose ground truth values fell within the
95% confidence interval computed via bootstrapping. For example, for a and b , a value of 10% indicates that 10% of the biomarkers had a correct estimatefor a k and b k . For α and β , 10% indicates that 10% of the individuals had acorrect estimate for α i and β i . For s , 10% indicates that 10% of the visits hada correct estimate for s ij . Simulation results are presented in Table 3. When there are correlations in the data, not modeling them yields inaccurate estimatesfor the trajectory slope parameter a , the subject-specific variables α i , β i , andthe progression scores s ij . Using the correlation model improves these estimatessignificantly. In our preliminary analysis of the noise spatial correlation structure usingthe semivariogram (Cressie and Hawkins, 1980; Bilgel et al., 2015c), we observedthat the rational quadratic had the best fit to the empirical semivariogram(Inline Supplementary Fig. 1). Rational quadratic also yielded the fit withlargest log-likelihood, and thus was selected as the correlation function for the model. AIC was − . × for the LME model, − . × for the model where C = I K × K , and − . × for the model where C = C ( ρ ), indicating that thePS model where spatial correlations are modeled fits the data the best. The19 able 3: Simulation results. Mean cosine similarity values across 1000 bootstrap experiments(with standard deviation) and percentage of variable elements that are correct based on 95%confidence intervals are presented. C = I K × K C = C ( ρ )Variable Cosine similarity % correct Cosine similarity % correct a . ± . . ± . b . ± . . ± . α . ± . . ± . β . ± . . ± . s . ± . . ± . Figure 3: Estimated trajectory slope parameters a k vs ground truth. Dashed line indicates x = y . Gray band corresponds to the 95% confidence intervals obtained from bootstrapping.Estimates are shown in blue if their 95% confidence interval intersects the x = y line, and inred otherwise. Results from the model where (a) C = I K × K and (b) C = C ( ρ ). spatial correlation parameter ρ for the rational quadratic correlation model wasestimated to be 4.5 mm, and λ was estimated to be 0.929. A comparison of the PS values obtained from the different correlation models are presented in theSupplementary Material (Inline Supplementary Fig. 2). Model fitting using 104subjects with a total of 300 visits each having about 30 ,
000 brain voxels tookapproximately 5 seconds on an Intel Xeon 8-core 3.1 GHz machine with 128 GBRAM for the case where voxels are assumed to be independent (i.e., C = I K × K ). igure 4: Predicted progression scores s ij vs ground truth. Dashed line indicates x = y . Grayband corresponds to the 95% confidence intervals obtained from bootstrapping. Estimates areshown in blue if their 95% confidence interval intersects the x = y line, and in red otherwise.Results from the model where (a) C = I K × K and (b) C = C ( ρ ). The following model fitting with C = C ( ρ ) took approximately 30 minutes perEM iteration, with convergence achieved in 2 to 6 iterations depending on thecorrelation structure. β -PS to mean cortical DVR We compared the subject-specific variables obtained from model fitting to empirical values obtained from mean cortical DVR, which is the widely usedmeasure for quantifying levels of brain amyloid (Fig. 5). For each individualwith at least two visits, we fit a line to the mean cortical DVR data to estimatethe slope of amyloid accumulation as well as the intercept. The subject-specificvariable α i , which represents the rate of amyloid progression, explained 62% of the variability in the empirical slope of amyloid accumulation computed fromlongitudinal mean cortical DVR. We observed a much higher correlation betweenmean cortical DVR and A β -PS ( R = 0 .
95 at baseline and 0 .
96 at last visit).When plotted against age, mean cortical DVR and A β -PS revealed similarpatterns (Fig. 6). The progression as measured using voxelwise DVR data is not linearly associated with age, and the A β -PS was able to capture this.21 igure 5: Correlation of estimated subject-specific variables with mean cortical DVR measures.The line of best fit is shown in blue, and its 95% confidence band in gray. (a) Rate of annualchange in mean cortical DVR vs. α , the predicted rate of change in amyloid progressionscore ( R = 0 . β , the progression score intercept( R = 0 . β -PS at baseline ( R = 0 . β -PS at last visit ( R = 0 . The trajectory slope parameters a k obtained from the PS model (Fig. 7)revealed symmetric patterns across the cerebral hemispheres. The precuneus andfrontal lobe showed the greatest increases in DVR with A β -PS, smaller increases in lateral temporal and temporoparietal regions, and minimal increases in theoccipital lobe and the sensorimotor strip. Voxelwise trajectories are furtherillustrated in Fig. 8, which shows the predicted DVR values on the cortical22 igure 6: (a) Mean cortical DVR and (b) A β -PS plotted against age. Longitudinal data pointsare connected by lines within each subject. Different colors indicate different subjects.Figure 7: Slope parameters a k obtained from voxelwise PS model projected onto the corticalsurface. For each unit increase in A β -PS, the DVR value at voxel k increases by a k . surface at three A β -PS levels.In order to better investigate regional trends, we averaged the PS model results within each cortical ROI and plotted these ROI averages as a function ofA β -PS (Fig. 9). We used bootstrap results to compute 95% confidence bandsfor these ROI averages. Based on the 95% confidence band of the precuneusillustrated in Fig. 9, precuneus appears to be the earliest accumulator and has thehighest amyloid levels through late stages of amyloid accumulation (A β -PS ≥ We present confidence bands for other cortical regions in Fig. 10.Based on our hypothesis testing procedure, we found that precuneus had thehighest amyloid levels at A β -PS values of -0.5, 0, 1, 2, and 3 (all p < . p = 0 . of the levels of amyloid at A β -PS= 0 and rates of accumulation across ROIs is23 igure 8: Predicted DVR levels at A β -PS = − . top ), 0.4 ( middle row ), and 1.5 ( bottom ).Figure 9: Regional trajectories as function of A β -PS. The PS model was used to make voxelwisepredictions at a range of A β -PS values, and these predictions were averaged within each ROIto obtain regional trajectories. The dashed lines indicate the 95% confidence band obtainedusing bootstrap results for the precuneus. presented in Fig. 11. 24 igure 10: Regional trajectories as function of A β -PS. The PS model was used to makevoxelwise predictions at a range of A β -PS values, and these predictions were averaged withineach ROI to obtain regional trajectories. The dashed lines indicate the 95% confidence bands forthe cortical regions. Estimated trajectories with their 95% confidence bands are superimposedon observed longitudinal data (in gray). The trajectories we obtained using the PS model were consistent with theresults of the LME model (Inline Supplementary Figs. 3, 4). However, the fixed effects obtained from the LME model did not describe the observed data aswell as the trajectory slopes obtained from the PS model (Inline SupplementaryFig. 5). Using the LME model estimates and a bootstrapping procedure similarto the one we applied for the PS model, we found that the precuneus had thehighest regional amyloid levels compared to other cortical regions at ages 70, 80, and 90 (all p < . p = 0 . igure 11: Comparison of levels of amyloid at A β -PS= 0 and rates of amyloid accumulationacross cortical regions. The intercept parameter b obtained from the PS model was averagedwithin each ROI to obtain the regional amyloid levels at A β -PS= 0, and the trajectory slopeparameter a was averaged within each ROI to obtain regional rates. but this was not statistically significant ( p = 0 .
4. Discussion
We presented a statistical model for estimating longitudinal trajectories ofvoxelwise neuroimaging data. Our model is based on the concept that age isnot a good metric for disease progression since each individual has his or herown onset and rate of progression of disease. We accounted for these inter- individual differences to temporally align individuals based on their collection ofvoxelwise measurements. As a result of this temporal alignment, we obtaineda metric of disease or underlying biological process stage as reflected in theneuroimaging data, and we call this metric “progression score”. By analyzingvoxelwise neuroimaging data as a function of progression score rather than time, we constructed trajectories that better represent changes occurring with diseasestage.Simulation results showed that the progression score model fits the databetter than the linear mixed effects model, as evidenced by AIC. Furthermore,26he results showed that modeling spatial correlations is important for extractingaccurate summary scores from data with underlying correlations. In the absenceof correlation modeling when the underlying data are correlated, estimates ofthe trajectory slopes a as well as subject-specific variables α, β , and s wereadversely affected. The difference in the subject-specific variables was mainlydue to inaccurate estimates of the prior covariance V . Note thatCov( y ij ) = aq Tij V q ij a T + R, (22)which explains why both a and V are affected when correlations are not modeledthrough the covariance matrix R . This highlights the importance of modelingspatial correlations in data where such effects are evident, especially in order to estimate correct trajectory slopes and individualized progression scores. Tofurther understand this phenomenon, we conducted simulations (data not pre-sented) where we duplicated biomarkers with lower signal-to-noise ratios andrepeated model fitting including these duplicate biomarkers. We observed thatthe parameter estimates as well as the PS values were biased if we assumed independence across biomarkers. This might be because the model interpretseach of the duplicated biomarkers as a separate piece of evidence towards thecomputation of individual progression scores. However, when we included aproper noise correlation model, we were able to recover our original results onthe collection of unique biomarkers. The proposed EM framework simplifies the estimation of spatial correlationparameters, and accurately estimates the trajectory parameters a and b as wellas predicting the progression scores s . The performance of the fitting procedurein predicting α and β is worse than for s . This may be due to the modeling of α and β as random variables that do not directly contribute to the observations y . We observed a similar pattern in the amyloid data: the agreement betweenestimated A β -PS and mean cortical DVR was greater than the agreement betweenestimated α and the rate of change of mean cortical DVR per individual. Thissuggests that the progression scores are more reliable than the subject-specificvariables. A β -PS was highly correlated with mean cortical DVR, a widely used β -PS is a meaningful score extracted from voxelwise imagingdata. Unlike mean cortical DVR, there are no a priori assumptions regardingwhich regions or voxels should be included in the computation of A β -PS. Thisproperty of PS allows for the discovery of new patterns in longitudinal data. Trajectories obtained from the PS model captured a wider dynamic range ofDVR values and had a better fit to the observed data compared to the trajectoriesobtained using the LME model fixed effects. This difference between the PSand LME models underscores the fact that age is not an appropriate metric forstaging amyloid accumulation. By aligning individuals based on their amyloid scans, the PS model enables a better temporal metric of amyloid staging.A β -PS allows for the exploration of longitudinal voxelwise trajectories withina hypothesis testing framework due to its underlying statistical model. Ourresults suggest that the precuneus exhibits the earliest cortical amyloid changes,but that its rate of amyloid accumulation does not differ significantly from other cortical regions. Based on a qualitative evaluation of our estimated trajectories,we found that amyloid accumulation in the precuneus is followed by cingulateand frontal cortices, then by lateral parietal cortex, followed by insula and lateraltemporal cortex. We observed minimal amyloid accumulation in visual cortex,hippocampal formation and the sensorimotor strip, which agrees with previous reports that these regions accumulate amyloid in later stages of AD (Braak andBraak, 1991). Previous reports have highlighted precuneus as an early amyloidaccumulator (Mintun et al., 2006; Rodrigue et al., 2012) as well as frontal,cingulate, and parietal regions (Jack et al., 2008). Contrary to our findinghighlighting precuneus as the earliest cortical accumulator, a study of cognitively normal adults found that the right medial frontal cortex accumulates amyloid theearliest, closely followed by bilateral precuneus based on cross-sectional voxelwiseanalyses (Villeneuve et al., 2015). We presented regional trajectories averagedbilaterally in this work; however, hemisphere-specific trajectories were consistentwith our bilaterally-averaged trajectories and we did not observe that medial frontal cortex precedes precuneus. Further studies with larger samples will be28nstrumental in elucidating the regional progression of amyloid accumulation.One of the strengths of our model for the progression scores s is that thereare no a priori assumptions on the global form of s ( t ), which can be thoughtof as an underlying function that describes the evolution of progression score over time. Instead, the model makes linear local approximations to s ( t ) perindividual. The linear form we use makes the EM approach analytically tractableand enables the discovery of the global form of s ( t ). On the other hand, thelinear relationship we assume between s and t is also a limitation since it may bean oversimplification over long follow-up durations. In order to capture dynamics over longer periods more accurately, it is necessary to investigate the relationshipbetween the progression scores s and time t , and select an appropriate functionto link these variables.Another limitation of our model is the assumption of linear biomarker trajec-tories. This assumption is not reflective of the fact that voxelwise DVR has a theoretical minimum of 1 (while in practice, there are DVR values lower than 1due to noise). Furthermore, several studies of longitudinal amyloid depositionhave found evidence for a ceiling effect in amyloid deposition in late AD, sug-gesting a sigmoid trajectory for amyloid levels (Villemagne et al., 2013; Villainet al., 2012). Our use of a linear rather than sigmoid trajectory inevitably results in inaccuracies in PS calculation for individuals whose voxelwise data lie alongthe plateaus of the sigmoid in reality. The linear trajectory assumption mayalso prevent us from characterizing subtler differences in the temporal orderingof the onset of amyloid accumulation across different regions. However, lineartrajectories, unlike sigmoid trajectories, yield closed-form update equations for trajectory parameters, greatly facilitating the model fitting procedure.The spatial covariance functions we used in our analyses are covariancefunctions of strictly stationary isotropic processes. The underlying spatial noiseprocess is affected by the PET scanner, image reconstruction algorithm, ra-diotracer delivery and binding in different brain regions, kinetic parameter estimation algorithm, and registration methods used to bring all scans in align-ment. To accurately model these influences on noise, complex noise spatial29ovariance models are needed. In this work, we made simplifying assumptions onthe noise properties to facilitate the study of longitudinal trajectories of amyloidimages. As a result of these simplifications, our model does not accurately capture noise properties; however, by eliminating the inaccurate assumption ofindependence across voxels, our approach improves the accuracy of the estimatedtrajectories and progression scores. Refined spatial noise models incorporatingnon-stationarity and anisotropy may further improve the estimation accuracy.Our model can be applied to studying higher resolution images. To reduce computational memory burden, it may be necessary to impose sparsity on thespatial correlation matrix, which can be done by imposing a correlation valueof 0 at large distances. The sparsity property of the correlation matrix can betaken advantage of to yield Cholesky decompositions using less memory andtime. In conclusion, the progression score model allows for extracting summaryscores from longitudinal data for each scan. In this work, we have extended theprogression score model proposed by Jedynak et al. (2012) to voxelwise imagingdata by accounting for spatial correlations and enabling efficient handling of thelarge number of voxels through the EM framework. The incorporation of a prior on subject-specific variables allowed for the inclusion of individuals with a singlevisit in the model. Our method can be extended to the analysis of other types ofimaging data, and in cases where a summary score such as mean cortical DVRis not available, the progression score estimates can be highly informative forprogression staging.
5. Acknowledgment
We thank Kalyani Kansal for her thorough review of the mathematicalderivations and assistance with code testing. This research was supported inpart by the Intramural Research Program of the National Institutes of Healthand by the Michael J. Fox Foundation for Parkinson’s Research, MJFF Research
Grant ID: 9310. 30 ppendix A.Proposition A.1. f (˜ u i | y i ; θ (cid:48) ) ∝ Φ(˜ u i ; ˆ u (cid:48) i , Σ (cid:48) i ) , (A.1) where Φ( · ; µ , Σ) denotes the probability density function of a multivariate normalwith mean µ and covariance matrix Σ , ˆ u (cid:48) i = (cid:88) j Z (cid:48) Tij R (cid:48)− Z (cid:48) ij + V (cid:48)− − (cid:88) j Z (cid:48) Tij R (cid:48)− ( y ij − b (cid:48) ) + V (cid:48)− m (cid:48) , (A.2) and Σ (cid:48) i = (cid:16)(cid:80) j Z (cid:48) Tij R (cid:48)− Z (cid:48) ij + V (cid:48)− (cid:17) − is a covariance matrix.Proof. Independence across visits allows us to write f ( y i | ˜ u i ; θ (cid:48) ) = (cid:89) j f ( y ij | ˜ u i ; θ (cid:48) )= (cid:89) j Φ( y ij ; Z (cid:48) ij ˜ u i + b (cid:48) , R (cid:48) ) ∝ exp − ˜ u Ti (cid:88) j Z (cid:48) Tij R (cid:48)− Z (cid:48) ij ˜ u i − (cid:88) j ( y ij − b (cid:48) ) T R (cid:48)− Z (cid:48) ij ˜ u i , (A.3)where we have ignored the terms that do not depend on ˜ u i . By Bayes’ rule, f (˜ u i | y i ; θ (cid:48) ) ∝ f ( y i | ˜ u i ; θ (cid:48) ) f (˜ u i ; θ (cid:48) ) (A.4) ∝ exp − ˜ u Ti (cid:88) j Z (cid:48) Tij R (cid:48)− Z (cid:48) ij + V (cid:48)− ˜ u i − (cid:88) j ( y ij − b (cid:48) ) T R (cid:48)− Z (cid:48) ij + m (cid:48) T V (cid:48)− ˜ u i (A.5) ∝ Φ(˜ u i ; ˆ u i , Σ (cid:48) i ) (A.6)Below we derive the EM algorithm update equations given in Section 2.3:31 olving for the intercept parameter b The value of b that solves ∂Q∂ b = R − (cid:80) i,j ( y ij − Z ij ˆ u (cid:48) i − b ) = 0 is given by b = 1 (cid:80) i v i (cid:88) i,j ( y ij − Z ij ˆ u (cid:48) i )= 1 (cid:80) i v i (cid:88) i,j (cid:0) y ij − a s (cid:48) ij (cid:1) . (A.7)Plugging in the expression for a from Equation A.9 yields b = (cid:16)(cid:80) i,j y ij (cid:17) (cid:16)(cid:80) i,j q Tij Σ (cid:48) i q ij + s (cid:48) ij (cid:17) − (cid:16)(cid:80) i,j y ij s (cid:48) ij (cid:17) (cid:16)(cid:80) i,j s (cid:48) ij (cid:17) ( (cid:80) i v i ) (cid:16)(cid:80) i,j q Tij Σ (cid:48) i q ij + s (cid:48) ij (cid:17) − (cid:16)(cid:80) i,j s (cid:48) ij (cid:17) . (A.8) Solving for the slope parameter a The value of a that solves ∂Q∂ a = R − (cid:104)(cid:80) i,j ( y ij − b ) s (cid:48) ij − a (cid:16)(cid:80) i,j q Tij Σ (cid:48) i q ij + s (cid:48) ij (cid:17)(cid:105) = 0is given by a = 1 (cid:80) i,j q Tij Σ (cid:48) i q ij + s (cid:48) ij (cid:88) i,j ( y ij − b ) s (cid:48) ij . (A.9)Plugging in the expression for b from Equation A.7 yields a = ( (cid:80) i v i ) (cid:16)(cid:80) i,j y ij s (cid:48) ij (cid:17) − (cid:16)(cid:80) i,j y ij (cid:17) (cid:16)(cid:80) i,j s (cid:48) ij (cid:17) ( (cid:80) i v i ) (cid:16)(cid:80) i,j q Tij Σ (cid:48) i q ij + s (cid:48) ij (cid:17) − (cid:16)(cid:80) i,j s (cid:48) ij (cid:17) . (A.10) Solving for the subject-specific variable mean parameter m The value of m that solves ∂Q∂ m = − (cid:80) i V − m + (cid:80) i V − ˆ u (cid:48) i = 0 is given by m = 1 n (cid:88) i ˆ u (cid:48) i . (A.11) Solving for the subject-specific variable variance parameter ν We use numerical optimization to estimate ν : ν = arg max ν Q ( θ , θ (cid:48) ) (A.12)= arg min ν (cid:88) i (cid:0) log | V | + (ˆ u (cid:48) i − m ) T V − (ˆ u (cid:48) i − m ) + Tr (cid:0) V − Σ (cid:48) i (cid:1)(cid:1) . (A.13)We can reconstruct V by undoing the log-Cholesky parametrization steps usingthe estimated parameter vector ν = [ ν , ν , ν ] T as V = e ν ν e ν e ν ν e ν . (A.14)32 olving for the noise covariance parameters λ and ρ If C = I K × K (i.e. ρ = 0 and fixed), then the diagonal elements λ k of Λ can be estimated as: λ k = (cid:115) (cid:80) i v i (cid:88) i,j (cid:2) ( y ijk − a k s (cid:48) ij − b k ) + a k q Tij Σ (cid:48) i q ij (cid:3) . (A.15)If C is not fixed to be the identity matrix, then in general it is not possible toobtain closed-form solutions for λ and ρ and we must use numerical optimizationover Q . When λ is a high-dimensional vector, this is not feasible. Therefore,we simplify the parametrization of the noise covariance matrix R as R ( λ, ρ ) = λ ˆΛ C ( ρ )ˆΛ, where we now consider ˆΛ as a fixed diagonal matrix with positivediagonal entries. We fix ˆΛ at the estimate of Λ from the model with C = I K × K . ρ is the correlation matrix parameter as defined previously, and λ > Q to estimateonly two parameters: λ, ρ = arg max λ,ρ Q ( θ , θ (cid:48) ) . (A.16)Note that since the EM algorithm does not require the maximization of Q ( θ , θ (cid:48) )but simply requires an increase in this function with each iteration, performingnumerical optimization until convergence is not necessary. Proposition A.2.
Consider the model given in Section 2.1.5 with θ = { m , V, a , b , R } .The reparametrization θ ∗ = w m + z , w V, w a , b − zw a , R , where w, z ∈ R , w (cid:54) = 0 yields an equivalent model.Proof. The incomplete log-likelihood islog f ( y ; θ ) = 12 (cid:88) i (log | π Σ i | − log | πV | − v i log | πR | ) − (cid:88) i,j ( y ij − b ) T R − ( y ij − b ) − (cid:88) i m T V − m + 12 (cid:88) i ˆ u Ti Σ − i ˆ u i . (A.17)33ere, we consider the case where z = 0 for algebraic simplicity. For thereparameterized model, we obtain ˆ u ∗ i = w ˆ u i and Σ ∗ i = w Σ i . Therefore, eachof the terms in log f ( y ; θ ∗ ) remain the same as those in log f ( y ; θ ), yielding log f ( y ; θ ∗ ) = log f ( y ; θ ). 34 nline Supplementary Figure 1: Preliminary analysis of the noise spatial correlation structureusing the semivariogram. The empirical semivariogram computed using the residuals fromthe model where C = I K × K is shown in blue. We fitted the semivariograms correspondingto the exponential, Gaussian, rational quadratic, and spherical correlation structures listedin Table 1 to the empirical semivariogram. The rational quadratic function had the best fitto the empirical curve based on the sum of squared error of the fitted semivariogram over a100 mm distance calculated using 30 equidistant points. nline Supplementary Figure 2: Bland-Altman plots comparing A β -PS values obtained usingdifferent spatial correlation models. We fitted the model using the “no correlation” structure( C = I K × K ) and four spatial correlation structures (exponential, Gaussian, rational quadratic,and spherical). We then used Bland-Altman plots to assess the differences in PS computedusing these correlation structures. In the Bland-Altman plots, the x -axis is the average of thePS values computed using the two methods being compared, and the y -axis is the difference.The labels above and to the left of the figures describe the spatial correlation structures beingcompared. Each blue dot corresponds to a visit. The solid black line indicates the meandifference between PS values computed using the methods, and the dashed red lines indicateits 95% confidence band. The Bland-Altman plots indicate that the “no correlation”, Gaussian,and spherical structures yield similar PS values. Furthermore, the exponential and rationalquadratic structures also yield similar PS values, but these differ from those obtained usingthe other three structures. nline Supplementary Figure 3: Slope parameters η k (fixed effects) obtained from voxelwiseLME model projected onto the cortical surface. DVR value at voxel k increases by η k per yearon average.Inline Supplementary Figure 4: Regional trajectories obtained from the LME model as functionof age. The LME model was used to make voxelwise predictions at a range of age values basedon the estimated fixed effects, and these predictions were averaged within each ROI to obtainregional trajectories. The dashed lines indicate the 95% confidence bands for the corticalregions. nline Supplementary Figure 5: Regional trajectories obtained from the LME model as functionof age. The LME model was used to make voxelwise predictions at a range of age values basedon the estimated fixed effects, and these predictions were averaged within each ROI to obtainregional trajectories. The dashed lines indicate the 95% confidence bands for the corticalregions. Estimated trajectories with their 95% confidence bands are superimposed on observedlongitudinal data (in gray). nline Supplementary Figure 6: Comparison of levels of amyloid at age 77 (the mean baselineage of the sample) and rates of amyloid accumulation across cortical regions using results ofthe LME model. We used the fixed effect estimates to obtain voxelwise amyloid levels at age77 as ˆ y k = η k ×
77 + γ k , and averaged ˆ y k within each ROI to obtain regional amyloid levels atage 77. The fixed effects η k were averaged within each ROI to obtain regional rates. eferencesReferences Avants, B.B., Epstein, C.L., Grossman, M., Gee, J.C., 2008. Symmetric dif-feomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Medical Image Analysis 12,26–41.Avants, B.B., Yushkevich, P., Pluta, J., Minkoff, D., Korczykowski, M., Detre,J., Gee, J.C., 2010. The optimal template effect in hippocampus studies ofdiseased populations. NeuroImage 49, 2457–2466.
Bateman, R.J., Xiong, C., Benzinger, T.L.S., Fagan, A.M., Goate, A., Fox, N.C.,Marcus, D.S., Cairns, N.J., Xie, X., Blazey, T.M., Holtzman, D.M., Santacruz,A., Buckles, V., Oliver, A., Moulder, K., Aisen, P.S., Ghetti, B., Klunk,W.E., McDade, E., Martins, R.N., Masters, C.L., Mayeux, R., Ringman, J.M.,Rossor, M.N., Schofield, P.R., Sperling, R.A., Salloway, S., Morris, J.C., 2012.
Clinical and biomarker changes in dominantly inherited Alzheimer’s disease.New England Journal of Medicine 367, 795–804.Bernal-Rusiel, J.L., Greve, D.N., Reuter, M., Fischl, B., Sabuncu, M.R., 2012.Statistical analysis of longitudinal neuroimage data with linear mixed effectsmodels. NeuroImage 66, 249–260.
Bernal-Rusiel, J.L., Reuter, M., Greve, D.N., Fischl, B., Sabuncu, M.R., 2013.Spatiotemporal linear mixed effects modeling for the mass-univariate analysisof longitudinal neuroimage data. NeuroImage 81, 358–370.Bilgel, M., An, Y., Lang, A., Prince, J., Ferrucci, L., Jedynak, B., Resnick,S.M., 2014. Trajectories of Alzheimer disease-related cognitive measures in a longitudinal sample. Alzheimer’s & Dementia 10, 735–742.Bilgel, M., An, Y., Zhou, Y., Wong, D.F., Prince, J.L., Ferrucci, L., Resnick,S.M., 2015a. Individual estimates of age at detectable amyloid onset for riskfactor assessment. Alzheimer’s & Dementia [In press].40ilgel, M., Carass, A., Resnick, S.M., Wong, D.F., Prince, J.L., 2015b. Defor- mation field correction for spatial normalization of PET images. NeuroImage119, 152–163.Bilgel, M., Jedynak, B., Wong, D.F., Resnick, S.M., Prince, J.L., 2015c. Temporaltrajectory and progression score estimation from voxelwise longitudinal imagingmeasures: Application to amyloid imaging, in: Ourselin, S., Alexander, D.C.,
Westin, C.F., Cardoso, M.J. (Eds.), Lecture Notes in Computer Science (9123),Information Processing in Medical Imaging, pp. 424–436.Biomarkers Definitions Working Group, 2001. Biomarkers and surrogate end-points: Preferred definitions and conceptual framework. Clinical Pharmacology& Therapeutics 69, 89–95.
Braak, H., Braak, E., 1991. Neuropathological staging of Alzheimer-relatedchanges. Acta Neuropathologica 82, 239–259.Caroli, A., Frisoni, G.B., 2010. The dynamics of Alzheimer’s disease biomarkersin the Alzheimer’s Disease Neuroimaging Initiative cohort. Neurobiology ofAging 31, 1263–1274.
Cressie, N., Hawkins, D.M., 1980. Robust estimation of the variogram. Journalof the International Association of Mathematical Geology 12, 115–125.Dale, A., Fischl, B., Sereno, M., 1999. Cortical surface-based analysis: I.Segmentation and surface reconstruction. NeuroImage 194, 179–194.Delor, I., Charoin, J.E., Gieschke, R., Retout, S., Jacqmin, P., 2013. Modeling
Alzheimer’s disease progression using disease onset time and disease trajectoryconcepts applied to CDR-SOB scores from ADNI. CPT: Pharmacometrics &Systems Pharmacology 2, e78.Desikan, R.S., S´egonne, F., Fischl, B., Quinn, B.T., Dickerson, B.C., Blacker,D., Buckner, R.L., Dale, A.M., Maguire, R.P., Hyman, B.T., Albert, M.S.,
Killiany, R.J., 2006. An automated labeling system for subdividing the human41erebral cortex on MRI scans into gyral based regions of interest. NeuroImage31, 968–980.Donohue, M.C., Jacqmin-Gadda, H., Le Goff, M., Thomas, R.G., Raman, R.,Gamst, A.C., Beckett, L.A., Jack, C.R., Weiner, M.W., Dartigues, J.F., Aisen,
P.S., 2014. Estimating long-term multivariate progression from short-termdata. Alzheimer’s and Dementia 10, S400–S410.Doody, R.S., Pavlik, V., Massman, P., Rountree, S., Darby, E., Chan, W., 2010.Predicting progression of Alzheimer’s disease. Alzheimer’s Research & Therapy2.
Fonteijn, H.M., Modat, M., Clarkson, M.J., Barnes, J., Lehmann, M., Hobbs,N.Z., Scahill, R.I., Tabrizi, S.J., Ourselin, S., Fox, N.C., Alexander, D.C., 2012.An event-based model for disease progression and its application in familialAlzheimer’s disease and Huntington’s disease. NeuroImage 60, 1880–1889.Galecki, A., Burzykowski, T., 2013. Linear model with fixed effects and correlated errors, in: Linear Mixed-Effects Models Using R. Springer New York, NewYork, NY. Springer Texts in Statistics. chapter 10, pp. 177–196.Ito, K., Corrigan, B., Zhao, Q., French, J., Miller, R., Soares, H., Katz, E.,Nicholas, T., Billing, B., Anziano, R., Fullerton, T., 2011. Disease progressionmodel for cognitive deterioration from Alzheimer’s Disease Neuroimaging
Initiative database. Alzheimer’s and Dementia 7, 151–160.Jack, C.R., Knopman, D.S., Jagust, W.J., Petersen, R.C., Weiner, M.W., Aisen,P.S., Shaw, L.M., Vemuri, P., Wiste, H.J., Weigand, S.D., Lesnick, T.G.,Pankratz, V.S., Donohue, M.C., Trojanowski, J.Q., 2013. Tracking pathophys-iological processes in Alzheimer’s disease: an updated hypothetical model of dynamic biomarkers. Lancet Neurology 12, 207–216.Jack, C.R., Lowe, V.J., Senjem, M.L., Weigand, S.D., Kemp, B.J., Shiung, M.M.,Knopman, D.S., Boeve, B.F., Klunk, W.E., Mathis, C.a., Petersen, R.C., 2008. C PiB and structural MRI provide complementary information in imaging42f Alzheimer’s disease and amnestic mild cognitive impairment. Brain 131,
Jedynak, B.M., Liu, B., Lang, A., Gel, Y., Prince, J.L., 2014. A computationalmethod for computing an Alzheimer’s disease progression score; experimentsand validation with the ADNI data set. Neurobiology of Aging 36 Supplement,S178–S184.Jenkinson, M., Bannister, P., Brady, M., Smith, S., 2002. Improved optimization for the robust and accurate linear registration and motion correction of brainimages. NeuroImage 17, 825–841.Lindstrom, M., Bates, D., 1990. Nonlinear mixed effects models for repeatedmeasures data. Biometrics 46, 673–687.Mintun, M.A., Larossa, G.N., Sheline, Y.I., Dence, C.S., Lee, S.Y., Mach, R.H.,
Klunk, W.E., Mathis, C.A., DeKosky, S.T., Morris, J.C., 2006. [C]PIB in anondemented population: potential antecedent marker of Alzheimer disease.Neurology 67, 446–452.Pinheiro, J.C., Bates, D.M., 1996. Unconstrained parametrizations for variance-covariance matrices. Statistics and Computing 6, 289–296.
Resnick, S.M., Goldszal, A.F., Davatzikos, C., Golski, S., Kraut, M.A., Metter,E.J., Bryan, R.N., Zonderman, A.B., 2000. One-year age changes in MRIbrain volumes in older adults. Cerebral Cortex 10, 464–472.Rodrigue, K.M., Kennedy, K.M., Devous Sr., M.D., Rieck, J.R., Hebrank, A.C.,Diaz-Arrastia, R., Mathews, D., Park, D.C., 2012. β -amyloid burden in healthy aging. Neurology 78, 387–395. 43chiratti, J.B., Allassonniere, S., Colliot, O., Durrleman, S., 2015a. Learningspatiotemporal trajectories from manifold-valued longitudinal data, in: Cortes,C., Lawrence, N., Lee, D., Sugiyama, M., Garnett, R. (Eds.), Advances inNeural Information Processing Systems 28, pp. 2395–2403. Schiratti, J.B., Allassonniere, S., Routier, A., Colliot, O., Durrleman, S., 2015b.A mixed-effects model with time reparametrization for longitudinal univariatemanifold-valued data, in: Ourselin, S., Alexander, D.C., Westin, C.F., Cardoso,M.J. (Eds.), Lecture Notes in Computer Science (9123), Information Processingin Medical Imaging, pp. 564–575.
Schmidt-Richberg, A., Guerrero, R., Ledig, C., Molina-Abril, H., Frangi, A.F.,Rueckert, D., 2015. Multi-stage biomarker models for progression estimationin Alzheimer’s disease. Lecture Notes in Computer Science (9123), InformationProcessing in Medical Imaging 9123, 387–398.Schulam, P., Wigley, F., Saria, S., 2015. Clustering longitudinal clinical marker trajectories from electronic health data: Applications to phenotyping andendotype discovery, in: AAAI Conference on Artificial Intelligence.Shock, N.W., Greulich, R.C., Andres, R., Arenberg, D., Costa Jr., P.T., Lakatta,E.G., Tobin, J.D., 1984. Normal human aging: The Baltimore Longitudi-nal Study of Aging. Technical Report. U.S. Government Printing Office.
Washington, DC.Sperling, R., Mormino, E., Johnson, K., 2014a. The evolution of preclinicalAlzheimer’s disease: implications for prevention trials. Neuron 84, 608–622.Sperling, R.A., Rentz, D.M., Johnson, K.A., Karlawish, J., Donohue, M., Salmon,D.P., Aisen, P., 2014b. The A4 Study: Stopping AD Before Symptoms Begin?
Science Translational Medicine 6, 228fs13.Villain, N., Ch´etelat, G., Grassiot, B., Bourgeat, P., Jones, G., Ellis, K.A.,Ames, D., Martins, R.N., Eustache, F., Salvado, O., Masters, C.L., Rowe,C.C., Villemagne, V.L., 2012. Regional dynamics of amyloid- β deposition44n healthy elderly, mild cognitive impairment and Alzheimer’s disease: a voxelwise PiB-PET longitudinal study. Brain 135, 2126–2139.Villemagne, V.L., Burnham, S., Bourgeat, P., Brown, B., Ellis, K.A., Salvado,O., Szoeke, C., Macaulay, S.L., Martins, R., Maruff, P., Ames, D., Rowe, C.C.,Masters, C.L., 2013. Amyloid β deposition, neurodegeneration, and cognitivedecline in sporadic Alzheimer’s disease: a prospective cohort study. Lancet Neurology 12, 357–367.Villeneuve, S., Rabinovici, G.D., Cohn-Sheehy, B.I., Madison, C., Ayakta, N.,Ghosh, P.M., La Joie, R., Arthur-Bentil, S.K., Vogel, J.W., Marks, S.M.,Lehmann, M., Rosen, H.J., Reed, B., Olichney, J., Boxer, A.L., Miller, B.L.,Borys, E., Jin, L.W., Huang, E.J., Grinberg, L.T., DeCarli, C., Seeley, W.W.,
Jagust, W., 2015. Existing Pittsburgh Compound-B positron emission tomog-raphy thresholds are too high: statistical and pathological evaluation. Brain138, 2020–2033.Yang, E., Farnum, M., Lobanov, V., Schultz, T., Raghavan, N., Samtani, M.N.,Novak, G., Narayan, V., DiBernardo, A., 2011. Quantifying the pathophysi- ological timeline of Alzheimer’s disease. Journal of Alzheimer’s Disease 26,745–753.Younes, L., Albert, M., Miller, M.I., 2014. Inferring changepoint times ofmedial temporal lobe morphometric change in preclinical Alzheimer’s disease.NeuroImage: Clinical 5, 178–187.
Young, A.L., Oxtoby, N.P., Daga, P., Cash, D.M., Fox, N.C., Ourselin, S., Schott,J.M., Alexander, D.C., 2014. A data-driven model of biomarker changes insporadic Alzheimer’s disease. Brain 137, 2564–2577.Zhou, Y., Endres, C.J., Braˇsi´c, J.R., Huang, S.C., Wong, D.F., 2003. Linearregression with spatial constraint to generate parametric images of ligand- receptor dynamic PET studies with a simplified reference tissue model. Neu-roImage 18, 975–989. 45iegler, G., Penny, W.D., Ridgway, G.R., Ourselin, S., Friston, K.J., 2015.Estimating anatomical trajectories with Bayesian mixed-effects modeling.NeuroImage 121, 51–68.745