Sara Garbarino
email: sara.garbarino@inria.fr

Marco Lorenzi
email: marco.lorenzi@inria.fr

Investigating hypotheses of neurodegeneration by learning dynamical systems of protein propagation in the brain

Keywords: Neurodegeneration, Causal model, Dynamical systems, Protein propagation, Gaussian process, Brain connectivity

A framework for spatio-temporal modeling of protein dynamics over brain networks from short term imaging data is proposed (GPPM-DS).

• GPPM-DS enables the investigation of bio-mechanical hypotheses governing disease progression via Bayesian model comparison.

• GPPM-DS provides new insights on the mechanisms of amyloid deposition in Alzheimer's disease, indicating the "Accumulation-Clearance-Propagation" model as the best suited dynamical system for interpretation of amyloid dynamics.

• GPPM-DS achieves accurate predictions of individual protein deposition in unseen data and provides plausible simulations of protein propagation.

Introduction

It is common hypothesis that while under normal conditions an efficient clearance process allows the brain to control the accumulation of toxic misfolded proteins, in neurodegenerative disorders (ND) the equilibrium between accumulation and clearance is broken [START_REF] Bateman | Human amyloid-beta synthesis and clearance rates as measured in cerebrospinal fluid in vivo[END_REF]. In this case, misfolded proteins may aggregate in plaques, eventually propagating between regions, leading to cellular dysfunction, disruption of synaptic connections, and neuronal loss [START_REF] Soto | Protein misfolding, aggregation, and conformational strains in neurodegenerative diseases[END_REF]. Despite the involvement of distinct proteins in different ND, the process of protein misfolding is thought to remain similar, although still not completely understood [START_REF] Sweeney | Protein misfolding in neurodegenerative diseases: implications and strategies[END_REF][START_REF] Editorial | Focus on neurodegenerative disease[END_REF]: proteins aggregates could selfpropagate and subsequentially spread the pathology between cells and tissues along pathways largely overlapping with functional or structural brain networks [START_REF] Jucker | Self-propagation of pathogenic protein aggregates in neurodegenerative diseases[END_REF][START_REF] Brettschneider | Spreading of pathology in neurodegenerative diseases: a focus on human studies[END_REF][START_REF] Fornito | Connectomics: a new paradigm for understanding brain disease[END_REF].

While a comprehensive model describing protein dynamics on the whole brain and along the whole span of the diseases is still missing, the current dominant paradigm for the study of proteinopathies in-vivo consists in associations analysis. For example, the relationship between pathological states and protein burden in the brain is usually quantified by means of correlation or regression models [START_REF] Näslund | Correlation between elevated levels of amyloid beta-peptide in the brain and cognitive decline[END_REF][START_REF] Rentz | Cognition, reserve, and amyloid deposition in normal aging[END_REF][START_REF] Palmqvist | Accuracy of brain amyloid detection in clinical practice using cerebrospinal fluid beta-amyloid 42: a cross-validation study against amyloid positron emission tomography[END_REF][START_REF] Hansson | CSF biomarkers of Alzheimer's disease concord with amyloid-beta PET and predict clinical progression: A study of fully automated immunoassays in BioFINDER and ADNI cohorts[END_REF][START_REF] Melzer | Beta amyloid deposition is not associated with cognitive impairment in Parkinson's disease[END_REF][START_REF] Müller | Amyloid-beta PET-Correlation with cerebrospinal fluid biomarkers and prediction of Alzheimer´s disease diagnosis in a memory clinic[END_REF]. Nevertheless, findings issued from these studies shed rather limited insights about the pathophysiological mechanisms underlying the observations. Orthogonal analysis paradigms are based on simulations, and have been proposed to investigate the dynamics of disease progression from a mechanistic point of view. The ambition of these approaches consists in discerning the generative spatio-temporal physical processes underlying the observations. These models thus offer the possibility to investigate the pathology from a causal perspective formalized through bio-mechanical relationships among observations [START_REF] Zhou | Predicting regional neurodegeneration from the healthy brain functional connectome[END_REF][START_REF] Iturria-Medina | Multifactorial causal model of brain (dis)organization and therapeutic intervention: Application to Alzheimer's disease[END_REF]. For example, a variety of generative models based on dynamical systems (DS) has been proposed for describing the kinetics governing the dynamical processes of accumulation, clearance and propagation of proteins [START_REF] Raj | A network diffusion model of disease progression in dementia[END_REF][START_REF] Iturria-Medina | Epidemic spreading model to characterize misfolded protein propagation in aging and associated neurodegenerative disorders[END_REF][START_REF] Raj | Network diffusion model of progression predicts longitudinal patterns of atrophy and metabolism in Alzheimer's disease[END_REF][START_REF] Oxtoby | Data-Driven sequence of changes to anatomical Brain connectivity in sporadic Alzheimer's Disease[END_REF][START_REF] Cauda | Brain structural alterations are distributed following functional, anatomic and genetic connectivity[END_REF][START_REF] Weickenmeier | The multiphysics of prion-like diseases: progression and atrophy[END_REF]. DS define protein kinetics through (typically very large) systems of non-linear differential equations (ODE) which encode the underlying bio-mechanical processes in a set of kinetic parameters governing the system. Many DS models define the propagation dynamics through diffusion equations [START_REF] Raj | A network diffusion model of disease progression in dementia[END_REF][START_REF] Raj | Network diffusion model of progression predicts longitudinal patterns of atrophy and metabolism in Alzheimer's disease[END_REF][START_REF] Cauda | Brain structural alterations are distributed following functional, anatomic and genetic connectivity[END_REF]. This modelling choice allows to reduce the number of parameters to be estimated, but comes at the expenses of an oversimplification of the dynamics governing the dynamical processes of protein propagation. Indeed, while the pathological kinetics may be assimilated to diffusive processes in short term observations, the long term evolution of NDs are unlikely to have diffusive properties. For example, the asymptotically constant behaviour of NDs may not be described by the stationary and constant rate of change specified by diffusion equations. With the aim of overcoming this drawback, more recent works describe the processes of protein misfolding by combining classical population dynamics equations with anisotropic diffusion to describe constant and simultaneous propagation and aggregation [START_REF] Weickenmeier | The multiphysics of prion-like diseases: progression and atrophy[END_REF], or by developing DS with protein concentration-dependent effects [START_REF] Garbarino | Modeling and inference of spatio-temporal protein dynamics across brain networks[END_REF].

Prediction and personalisation stem naturally from the ordinary differential equation structure of DS, which allows extrapolation from known initial conditions and prescribed kinetic parameters. When studying brain protein measurements in-vivo, the possibility of comparing models issued from different DS is of paramount importance for investigating the bio-mechanical hypotheses underlying disease progression. This is a crucial aspect, for observational data can be equally well described by conceptually different mechanistic hypotheses, and is especially important in the clinical setting, when direct intervention on the brain system is ultimately not possible. Unfortunately, since DS are usually identified by a large number of kinetic parameters and different relationship across variables, a consistent theoretical framework for DS comparison is often missing.

Another source of complexity in the modeling of protein propagation in-vivo arises from the lack of a precise definition of the time axis, thus preventing models of disease dynamics to reproduce patterns of protein propagation compatible across the entire history of ND. A well-defined temporal reference is indeed not available in typical longitudinal data sets of ND: the patient's time of onset is generally unknown, while the rate of pathological progression is highly variable across individuals. Moreover, ND are commonly studied through observational studies focusing on specific clinical stages, collecting individual measurements over a limited time-span, resulting in few time-points collected for each subject. To overcome the problem of modeling short term measurements with undefined temporal reference, several data-driven disease progression models (DPM) have been proposed [START_REF] Jedynak | A computational neurodegenerative disease progression score: method and results with the Alzheimer's disease Neuroimaging Initiative cohort[END_REF][START_REF] Villemagne | Amyloid beta deposition, neurodegeneration, and cognitive decline in sporadic Alzheimer's disease: a prospective cohort study[END_REF][START_REF] Donohue | Estimating long-term multivariate progression from short-term data[END_REF][START_REF] Young | A data-driven model of biomarker changes in sporadic Alzheimer's disease[END_REF][START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF][START_REF] Lorenzi | Probabilistic disease progression modeling to characterize diagnostic uncertainty: application to staging and prediction in Alzheimer's disease[END_REF][START_REF] Marinescu | DIVE: A spatiotemporal progression model of brain pathology in neurodegenerative disorders[END_REF][START_REF] Garbarino | Differences in topological progression profile among neurodegenerative diseases from imaging data[END_REF], based on the concept of selfmodeling regression [START_REF] Lawton | Self Modeling Nonlinear Regression[END_REF]. These approaches allow to reconstruct biomarkers trajectories along an ideal long term disease progression by optimally "stitching" together short term individual measurements. Each subject is then characterized by specific time parameters quantifying their pathological stage with respect to the estimated long term group-wise evolution. These models provide a description of temporal progression across biomarkers, usually expressed in terms of regression curves, or event orderings, without however elucidating the kinetics and relationships across them. This aspect indicates a limited ability of DPM in providing an understanding of the pathological mechanisms.

To date, no modeling framework allows inference and comparison of mechanistic hypotheses of protein processes across the long term evolution of ND. The problem is challenging since it requires to simultaneously account for short term observations to reconstruct the long term disease progression, and to estimate the kinetic parameters specified by high dimensional dynamical systems.

In this study we tackle this issue by formulating Gaussian process progression modeling for dynamical systems (GPPM-DS): a framework for specifying hypothetical models of protein dynamics, fitting them to the data and comparing their evidence through Bayesian model comparison. As a result, we can identify the most plausible dynamics underlying the measurements, and provide an interpretable model of the potential causal relationship across variables. This framework enables to formalize and test novel hypotheses on disease dynamics, and to challenge them against alternative ones. This idea introduces a paradigm shift in the analysis of neurodegeneration in neuroimaging data analogous to the one brought by Dynamic Causal Modeling [START_REF] Friston | Dynamic causal modelling[END_REF] for the hypothesis-driven analysis in functional neuroimaging. Compared to this latter framework, here we tackle the additional problem of the lack of a well-defined time axis for the pathology, and present a completely novel theoretical framework for scalable inference.

Overall, in this work we provide a framework for i) testing hypothesis-driven models of protein dynamics against longitudinal data, ii) simulating protein propagation along the whole disease progression time span, and iii) predicting individual protein deposition in unseen data. We test GPPM-DS on a variety of synthetic data and evaluate its performances in estimating protein evolution as well as kinetic and individual time parameters, for each proposed DS, as compared to standard DPM based on monotonic constraints [START_REF] Lorenzi | Probabilistic disease progression modeling to characterize diagnostic uncertainty: application to staging and prediction in Alzheimer's disease[END_REF]. Results on synthetic data can be found in the Supplementary Material.

We then demonstrate GPPM-DS on the modeling of the evolution of cerebral amyloid protein accumulation from AV45-PET data of Alzheimer's disease (AD) subjects from the ADNI data set. We show that GPPM-DS identifies the non-linear process of accumulation, clearance and propagation as the best suited DS for bio-mechanical description of amyloid dynamics.

Methods

GPPM-DS is described as a Bayesian non-parametric constrained regression problem (Section 2.1), where constraints are imposed on the dynamics of protein evolution and are given by a DS for protein dynamics (Section 2.2), while the long term data are reconstructed from short term observation by using time reparametrization techniques typical of the DPM framework. Bayesian inference techniques are used to solve the GPPM-DS problem (Section 2.3). Specifically, we determine a lower bound for the generally intractable model posterior. We resort to optimizing an approximation of the marginal function, and imposing sparsity constraints on the model parameters. Efficient optimisation is performed through stochastic gradient descent via backpropagation and mini-batch optimisation.

Bayesian non-parametric constrained regression framework

Let us assume to have S subjects (s = 1, ..., S), each one with associated measurements of protein concentrations Y s , at different brain regions and at different time-points over a short term time span t s = {t s 1 , ..., t s T s }, where T s denotes the total number of measurements for a given subject s. We then assume that each individual measurement is obtained with respect to an absolute time-frame τ through a function t = g s (τ ) modeling the time-reparameterization with respect to the common group-wise evolution. This means assuming that the short term time-points are obtained through reparametrization of some unknown time-points τ s k ∈ I ⊆ R over a long term span, such that:

g s : I → J s ⊂ I τ s k → t s k = g s (τ s k ). (1)
Finally, we assume that the measurements Y s are realizations of an unknown group-wise process f , describing the temporal evolution of protein concentrations. According to these assumptions, we model the observations Y s via mixed-effect regression with time reparametrization:

Y s (g s (τ )) = f (g s (τ )) + ν s (g s (τ )) + . (2)
A variety of models has been proposed to tackle this kind of regression problem, either via parametric or non-parametric techniques. Here we use GP Progression Model [START_REF] Lorenzi | Probabilistic disease progression modeling to characterize diagnostic uncertainty: application to staging and prediction in Alzheimer's disease[END_REF], which is a non-parametric Bayesian DPM describing f as a Gaussian process (GP), and g s as a translation parameterized by an individual time-shift: g s (τ s k ) = t s k + d s . Finally, ν s are assumed to be Gaussian random effects N (0, φ s ).

We introduce constraints on the dynamics of the model, enforcing the concentrations' evolution to a general DS, governed by a functional H θ , which depends on M kinetics parameters {θ i } M m=1 = θ. This means specifying a family of admissible functions whose derivatives evaluated at points t satisfy the DS:

H = {f (t) : ḟ (t) = H θ (f (t), t)}. (3)
The GPPM-DS model is finally described as equation ( 2) subject to constraints (3).

Dynamical systems for modeling misfolded proteins

We introduce three different DS accounting for protein dynamics: i) diffusion (Diff): a purely diffusive model of constant propagation, based on the work presented in [START_REF] Raj | A network diffusion model of disease progression in dementia[END_REF]; ii) reaction-diffusion (RD): a model where aggregation and propagation of proteins are simultaneous, constant and opposite, and the total process eventually reaches a plateau, to reproduce the DS proposed by [START_REF] Weickenmeier | The multiphysics of prion-like diseases: progression and atrophy[END_REF]; and iii) accumulation clearance and propagation (ACP): a model where propagation is triggered when the aggregation saturates, and then both aggregation and propagation eventually reach a plateau [START_REF] Garbarino | Modeling and inference of spatio-temporal protein dynamics across brain networks[END_REF]. The three chosen models offer a range of varying complexity in terms of understanding of the underlying pathology: from the simplest one (Diff), which allows to estimate purely diffusive effects of an hypothetical constant propagation of proteins, to the most complex (ACP), which incorporates the effects of multivariate propagation, aggregation saturation and clearance. In this context, the ACP model is related to the highest range of parameters describing protein propagation, and so potentially offers the highest capacity in describing the pathological progression. We aim at identifying, amongst these three models, the one that best explain the dynamics of the data, while at the same time accounting for model complexity.

Here we consider the brain as a system of N interconnected regions, where each region i (i = 1, ..., N ) is characterized by its concentration of proteins f i (τ ) along the whole long term time interval I. We describe the network of brain connections as a graph, whose adjacency matrix is denoted by A, its degree matrix by D, and the corresponding Laplacian matrix by H = D -A. In the following, denotes the Hadamard product between either matrices or vectors.

The three DS are defined as follows:

Diff :

ḟ (τ ) = H θ D (f (τ ), τ ) (4)
= (A K θ D ) Hf (τ ).

Here the matrix K θ D encodes region-wise, constant rates of protein propagation across adjacent regions and is described as

(K θ D ) ij = k ij ,
where k ij is the (symmetrical) rate of propagation between regions i and j. The Diff model is linear and produces patterns of protein concentration that accumulate indefinitely over time. The set of kinetic parameters for the Diff model is

θ D = (k ij ).
RD :

ḟ (t) = H θ RD (f (τ ), τ ) (5) = (A K θ RD ) Hf (τ ) + R θ RD f (τ ) (υ -f (τ )).
This model includes both propagation and aggregation mechanisms, which are constant and simultaneous. We assume that no aggregation nor propagation occur in healthy conditions, while protein plaques aggregation develops when the accumulation-clearance equilibrium breaks. The model is comprised of two terms: a standard diffusion term (A K θ RD ) Hf (τ ) for describing constant protein propagation (see above) and a reaction term R θ RD f (t) (υ-f (τ )) for describing protein aggregation. The second term balances constant propagation with a constant aggregation term encoded by the matrix R θ RD = k t I, where k t is the rate of total aggregation and is constant across regions. This term eventually reaches a plateau when protein concentration get to a maximal concentration threshold υ. This equation is known as the Fisher-Kolmogorov equation [START_REF] Adomian | Fisher-Kolmogorov equation[END_REF]. The set of kinetic parameters for the RD model is

θ RD = (k ij , k t , υ). ACP : ḟ (τ ) = H θ ACP (f (τ ), τ ) (6) = (A K θ ACP (f , τ )) Hf (τ ) + R θ ACP (f , τ )f (τ )
The ACP model describes the three processes of accumulation, clearance and propagation of proteins allowing for non-constant effects. We assume again no aggregation nor propagation to occur in healthy conditions, and aggregation to develop when the accumulation-clearance equilibrium breaks. The propagation term K θ ACP is concentration-dependent: the toxic protein concentration in each region saturates when reaching a first critical threshold γ = (γ 1 , ..., γ N ), and subsequently triggers propagation towards the connected regions. Propagation also reaches a plateaus when passing a second critical threshold η = (η 1 , ..., η N ). This can be modeled by setting

(K θ ACP ) ij = k ij ACP (f (τ ), γ, η), where k ij ACP are sigmoid functions of the constant propagation rates k ij k ij ACP (f j (τ ), γ j , η j ) = k ij 1 + e -l 1 (f j (τ )-γ j ) 1 + e l 2 (f j (τ )-η j ) . (7)
For the aggregation term R θ ACP , we hypothesize a similar concentration-dependent threshold above which the total aggregation process reaches a plateau. In this context, (R θ ACP ) ij is re-written as

(R θ ACP ) ij = k t ACP (f (τ ), η) if i = j 0 otherwise; (8)
where

k t ACP (f i (τ ), η i ) = k t 1 + e l 2 (f i (τ )-η i ) . (9)
The key difference with the previous RD model, is that the ACP model does not assume aggregation and propagation to be simultaneous nor constant, but rather hypothesizes the existence of critical values of protein concentrations, different for each region, at which the regional concentrations saturate and subsequently trigger propagation. The set of kinetic parameters for the ACP model is

θ ACP = (k ij , k t , γ, η).
Figure 1 shows examples of Diff, RD and ACP models with varying key parameters on a toy data set of N = 3 synthetic brain regions. In Figure 1A), for the linear Diff model, the highest the rates of propagation k ij , the fastest the proteins travel amongst the simulated brain regions, whose protein concentration tend to balance over time. Figure 1B)-C) show trajectories obtained by varying the k ik (1B) and υ (1C) parameters for the RD model (while the other parameter is kept constant). When k ij increases, the model predicts a behaviour similar to the Diff model in 1A), with fastest propagation amongst region and concentration trajectories that balance over time. In Figure 1C) we note how accumulation reaches a plateau once the maximal concentration threshold parameter υ is hit. Figure 1D)-F) show protein trajectories predicted by the ACP model with varying parameters. At increasing k ij again correspond fastest propagation amongst regions (1D), while when varying the l 1 and l 2 parameters, the trajectories display an increasingly sigmoidal shape (1E). Finally, when the parameters γ and η increase, the concentrations saturate, triggering propagation, and reach a plateau, later in the progression (1F).

Inference

Solving the GPPM-DS system of equations ( 2)-( 3) requires inference of the kinetic parameters θ, the group-wise dynamics f , and, for each subject s, of the parameters of the individual random effects φ s and of the time-shift parameters d s . We first define F s as the realization of f at times t s , and Ḟ s as the set of realizations of f and of its derivatives at (generally different) times u s . To ease the notation, we assume both f and its derivatives to have realizations on the same time-points t s , but computation extends easily to more complex scenarios. We also indicate by F , ν, Ḟ , d and φ the collections of F s , ν s , Ḟ s , d s and φ s for all the subjects (s = 1, ..., S). Following [START_REF] Lorenzi | Constraining the Dynamics of Deep Probabilistic Models[END_REF][START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF], we describe the constrained regression problem in a Bayesian setting and solve the inference problem for θ and F by determining a lower bound for the marginal

p(Y, H |φ, d, t, , ζ) = p(Y|F , φ, d, t, )p(H | Ḟ , θ, d, t, ζ) p(F , Ḟ |φ, d, t)p(θ)dF d Ḟ dθ, (10)
where

p(F , Ḟ |φ, d, t)dF = p( Ḟ |F )p(F |φ, d, t). (11)
Analogously to [START_REF] Lorenzi | Constraining the Dynamics of Deep Probabilistic Models[END_REF], we assume the two likelihood terms (the data fidelity term p(Y|F , φ, , t, d, ) and the constraint term p(H | Ḟ , θ, t, d, ζ)) to be respectively Gaussian and Student-t with respective variance and scale parameters and ζ.

Due to the general intractable form of (10), we approximate the GP F via random features (RF) expansion [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF]. Specifically, a GP with radial basis function covariance can be expressed as F ≈ h(tΩ)W, where Ω is a linear projection of the input t into the RF space specified by trigonometric activation functions h(•) = (cos(•), sin(•)), and W are the regression parameters [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF]. Such approximation extends to the derivatives of the GP thanks to the chain rule [START_REF] Lorenzi | Constraining the Dynamics of Deep Probabilistic Models[END_REF]. As a result, both GP realization and its derivatives can be identified by the same parameters W and Ω. The RF approximation allows to replace inference on F with inference on W and Ω. Following [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF], we optimize [START_REF] Palmqvist | Accuracy of brain amyloid detection in clinical practice using cerebrospinal fluid beta-amyloid 42: a cross-validation study against amyloid positron emission tomography[END_REF] through variational inference of W and θ, and assume Ω to be sampled from the prior with fixed randomness. This leads to the optimization of the following lower bound of the log-marginal (or evidence lower bound, ELBO):

log(p(Y, H |φ, d, t, , ζ)) ≥ E q(W) [log(p(Y|Ω, W, φ, d, t, ))] + E q(W)q(θ) [log(p(H |Ω, W, θ, d, t, ζ))] -DKL(q(W)|p(W)) -DKL(q(θ)|p(θ)). (12)
Here, we perform variational inference on the kinetic parameters θ and group-wise dynamics F (and so W), while finding a maximum likelihood estimate of time shift parameters d and of the individual random effects φ. The DKL(q|p) term represents the Kullback Leibler divergence between the prior p and the variational approximation q, and we assume q(W) and q(θ) to be Gaussians N (µ θ , σ θ ) and N (µ W , σ W ). In this way the last two terms of equation ( 12) have closed forms [START_REF] Cutajar | Random feature expansions for deep Gaussian processes[END_REF]. Thanks to separability and analytical forms of likelihood and DKL terms, the ELBO ( 12) can be optimized through stochastic gradient descent via backpropagation and mini-batch [START_REF] Kingma | Variational dropout and the local reparameterization trick[END_REF]. Full optimization is performed with stochastic gradient descent with adaptive moment estimation (Adam) [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], through the alternate optimization of i) the approximated posterior over W, GP parameters and random effects: q(W), and φ;

ii) the individual time-shifts parameters d;

iii) the approximated posterior over the kinetic parameters and likelihood parameters of the constraints: q(θ) and ζ.

Sparsity on the kinetic parameters

A well-known limitation of fitting non-linear DS concerns the identifiability of the associated kinetic parameters θ. Here, inspired by previous work on dropout and variational dropout [START_REF] Kingma | Variational dropout and the local reparameterization trick[END_REF][START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF][START_REF] Antelmi | Sparse Multi-Channel Variational Autoencoder for the Joint Analysis of Heterogeneous Data[END_REF], we impose parsimonious representations of the model kinetic parameters through sparsity constraints (sparse GPPM-DS). As we will see, sparse GPPM-DS sets to zero the kinetic parameters which are associated with high estimated variance -parameters which can be considered un-identifiable, given the observed data. It has been shown [START_REF] Kingma | Variational dropout and the local reparameterization trick[END_REF] that a way to impose sparsity on Bayesian neural networks (i.e. regularize their weights w) is to use the so-called local reparametrization trick: re-parametrize the Gaussian approximation q(w) ∼ N (µ w , σ w ) of the weights of the approximated posterior p(w) as q(w) ∼ N (µ w , αµ 2 w ), where α is the dropout rate. The authors in [START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF] proposed a form for the prior on w consistent with the optimization of the ELBO associated with new q: the log-scale uniform function p(|w|) ∝ 1 |w| . In our case, we are interested in imposing sparsity constraints on the kinetic parameters θ. For this reason, we apply the local re-parametrization trick to the Gaussian approximation of the posterior distribution of the kinetic parameters q(θ) ∼ N (µ θ , σ θ ), obtaining q(θ) ∼ N (µ θ , αµ 2 θ ). This modification promotes sparsity, as large values of θ i correspond to large uncertainty expressed by the variance term αµ 2 θ i , indicating that the associated kinetic parameter θ i can be set to zero. In what follows, when using sparse GPPM-DS, the associated D KL (q(θ)|p(θ)) is computed according to the numerical approximation as described in [START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF].

Simulation Results

We tested GPPM-DS in both sparse and full versions on a variety of synthetic data sets and compared its performances in recovering kinetic parameters, simulated evolution and time reparametrization parameters, for each of the three proposed DS, as compared to standard DPM. As reference DPM, we implement a GPPM based on monotonic constraints (GPPM-mono) [START_REF] Lorenzi | Probabilistic disease progression modeling to characterize diagnostic uncertainty: application to staging and prediction in Alzheimer's disease[END_REF]: it is, as the model proposed in this work, a Gaussian process-based random effect modelling of longitudinal progressions, and simply assumes a steady temporal evolution of biomarkers from normal to pathological values, enforcing the biomarker trajectories to follow a monotonic behaviour. GPPM-mono provides us with a "null progression" hypothesis, in which no propagation parameters are accounted for, while regions are simply assumed to steadily evolve from normal to pathological values. Specifically (full details can be found in the Supplementary Material, section S1):

1. we generated ground truth data using each of the proposed DS and then attempt reconstruction using every GPPM-DS, demonstrating, in terms of both RMSE and ELBO, that the candidate model is generally optimal when it corresponds to the ground truth DS (Supplementary section S1.1).

2. We tested the variational dropout scheme for identifying zero-valued kinetic parameters and demonstrate that the sparse model has good accuracy in identifying such parameters (Supplementary section S1.2).

3. We analyzed the performance of the sparse models in predicting unseen data as compared to the full models, demonstrating that they perform similarly in terms of RMSE while requiring only a fraction of parameters (Supplementary section S1.3).

Moreover, in order to probe whether the mechanistic inference is truly driven by biological dynamics, we performed a set of simulations in which the ground truth data, generated using one of the propagation models, have been shuffled 100 times, thus obtaining 100 data sets with no biological signal (Supplementary section S1.5). Then, we reconstructed with every GPPM-DS and showed that, in terms of ELBO, the three models best reconstruct the data roughly the same number of times. These results show that under the null hypothesis of no biological signal underlying the spatial data, the models are roughly equivalent from the point of view of Bayesian model comparison.

Results

GPPM-DS relies on the gradient-matching approach presented in [START_REF] Lorenzi | Constraining the Dynamics of Deep Probabilistic Models[END_REF] that overcomes the issue of using costly numerical integration procedures for inferring ODE parameters [START_REF] Macdonald | Gradient Matching Methods for Computational Inference in Mechanistic Models for Systems Biology: A Review and Comparative Analysis[END_REF], thus allowing scalable inference of protein progression jointly with associated kinetic parameters. The proposed framework is formulated as a constrained regression problem in a Bayesian nonparametric setting, where the protein progression is modeled as a Gaussian process (GP), while bio-mechanical processes are defined as constraints on the protein dynamics through DS expressed by ODE. The Bayesian setting allows for uncertainty quantification of the protein dynamics while, to achieve tractability, the inference problem is solved via stochastic variational inference. The framework also provides a principled theory for model comparison via assessment of model evidence: here we investigate use of the evidence lower bound (ELBO) as a surrogate for model comparison. Figure 2 shows a schematic representation of our framework.

We propose GPPM-DS by comparing three different DS of amyloid protein dynamics in AD. We investigate a purely diffusive model of constant propagation (Diff), based on the preliminary work of [START_REF] Raj | A network diffusion model of disease progression in dementia[END_REF]; a reaction-diffusion model (RD), where aggregation and propagation of proteins are simultaneous, constant and opposite, and the total process eventually reaches a plateau, to reproduce the DS proposed by [START_REF] Weickenmeier | The multiphysics of prion-like diseases: progression and atrophy[END_REF]; a non-linear accumulation clearance and propagation model (ACP) [START_REF] Garbarino | Modeling and inference of spatio-temporal protein dynamics across brain networks[END_REF], where propagation is triggered when aggregation reaches saturation, and then both aggregation and propagation reach a plateau. For each DS, GPPM-DS provides a complete description of the protein dynamics, which can be subsequently used for simulating and predicting proteins changes over time through forward integration. To address the ill-posed nature of the identification of the DS parameters [START_REF] Saccomani | Identifiability of Nonlinear ODE Models in Systems Biology: Results from Both Structural and Data-Based Methods[END_REF], we further leverage on variational dropout techniques [START_REF] Kingma | Variational dropout and the local reparameterization trick[END_REF][START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF][START_REF] Antelmi | Sparse Multi-Channel Variational Autoencoder for the Joint Analysis of Heterogeneous Data[END_REF] for imposing parsimonious representations of the model parameters through sparsity constraints (sparse GPPM-DS).

We show first that GPPM-ACP model outperforms the others in terms of ELBO and predictive accuracy of amyloid deposition in unseen follow-up data. This is confirmed by the analysis of the null scenario case where the model is optimised with respect to the sole monotonic constraint of the trajectories (GPPM-mono) [START_REF] Lorenzi | Probabilistic disease progression modeling to characterize diagnostic uncertainty: application to staging and prediction in Alzheimer's disease[END_REF]. Full results on GPPM-mono are shown in the Supplementary Material, Figures S5 andS9. We then show that GPPM-ACP model allows bio-mechanical interpretation of amyloid dynamics in AD, while providing plausible description of the pathological evolution.

Data Acquisition and Preprocessing

We analyzed AV45-PET brain imaging data of 770 subjects from ADNI, with a total of 1477 longitudinal data points. The neuronal connections along which protein propagation occurs is approximated by estimating the average structural connectomes from 24 young and healthy volunteers of the HCP data set. In this case, the substrate for protein propagation is assumed to be unaffected by the pathology.

HCP data: Data used in the preparation of this work were obtained from the MGH-USC Human Connectome Project database. We collected 3D T1w and DTI of 24 age and gendermatched subjects. Details on the selected subjects and the pipeline for structural connectome generation are described in [START_REF] Oxtoby | Data-Driven sequence of changes to anatomical Brain connectivity in sporadic Alzheimer's Disease[END_REF]. We averaged the 24 connectomes and obtained an average young, healthy connectome on 82 FreeSurfer regions [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF]. For computational reasons, and in order to remove weak connections, we averaged left and right hemispheres (remaining with 41 regions) and we set to 0 all the weights below the estimated average weights across nodes, and to 1 the weight above. We collected clinical, demographic and AV45-PET SUVr data from the adnimerge package (adni.bitbucket.io [START_REF]R package version 0.0.1 of AD-NIMERGE: Alzheimer's Disease Neuroimaging[END_REF]) for subjects with diagnosis label "Dementia", Mild Cognitive Impairment" (MCI) or "Cognitively Normal" (CN). As we are interested in modeling amyloid evolution in well phenotyped Alzheimer's trajectories, we selected, among the MCI or Dementia subjects, only the one effectively displaying amyloid accumulation. On the other hand, we selected the whole cohort of CN (both amyloid positive: CN+ and amyloid negative: CN-), in order to maximize the time-span of the data. The threshold for amyloid positivity was established at amyloid level in the cerebrospinal fluid < 192 ng/l [START_REF] Almdahl | Cerebrospinal fluid levels of amyloid beta 1-43 mirror 1-42 in relation to imaging biomarkers of Alzheimer's disease[END_REF]. The information was extracted from the adnimerge R package (MEDIAN field of the upennbiomkmaster table). ADNI AV45-PET SUVr data are already computed on a brain anatomical parcellation estimated with the software Freesurfer [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF], and normalized against cerebellar uptake. For each region we discarded white matter, ventricular and cerebellar regions, remaining with 82 regions, and averaged the SUVr of each region from both hemispheres, remaining with 41 regions. Demographic and clinical details are shown in Table 1. Finally we split the data set in two parts. A training data set D 1 was defined by including all the longitudinal data for each subject up to the second-to-last time points. The remaining time-points were included in a testing data set D 2 . Subjects with one measurement only were included in D 1 . The resulting data set D 1 includes 984 longitudinal measures; D 2 contains 493 cross-sectional measures. We run the models on D 1 , estimating protein dynamics, GP parameters and individual time-shifts, and used D 2 to validate model predictions.

Group

GPPM-DS for modeling amyloid deposition in Alzheimer's disease

The proposed GPPM-DS framework can be effectively used to compare different hypotheses of protein kinetics. We analyze AV45-PET data with three different kinetic models, whose mathematical description can be found in the Method Section 4.2: Diff model of equation ( 4), RD model of equation ( 5) and ACP of equation [START_REF] Brettschneider | Spreading of pathology in neurodegenerative diseases: a focus on human studies[END_REF]. We use the ELBO as surrogate for model comparison (as for synthetic data, see Supp. Mat.): we compared the ELBO while training of D 1 , and the RMSE for prediction on the unseen data (D 2 ). We computed the RMSE for each individual prediction (averaging over the biomarkers), and bootstrapped over the kinetic parameters 100 times, obtaining RMSE distributions. Individual predictions on D 2 are performed by forward integration of the associated dynamical system with the estimated kinetic parameters, from given initial conditions (the individual baseline values) for the sparse model (i.e. having removed the kinetic parameters with dropout threshold of p < 0.65, see synthetic data results, Supp. Mat.).

GPPM-ACP model has best performances according to both ELBO and predictive accuracy on follow-up data (Table 2). The GPPM-Diff model produces the worst RMSE for prediction, which may be due to the unrealistic patterns of indefinite accumulation that the underlying dynamical system allows for. Table 2 also reports RMSE for prediction for each disease status. Results are consistent across diagnosis for GPPM-ACP and GPPM-RD, while for GPPM-Diff predictions deteriorate sensibly for MCI and AD subjects.

The last line of Table 2 reports results for the null model (GPPM-mono), which produces the worst ELBO. Moreover, by not estimating any dynamical system, it does not provide individual predictions on D 2 via forward integration, and so the RMSE on prediction is not provided.

Model

-mono - - - - - - 13 (0.09) 0.13 (0.10) 0.13 (0.10) 0.13 (0.10) -11874 GPPM
Tab. 2: RMSE mean (and std) for individual prediction on unseen data set D 2 and ELBO of the models while training on D 1 . Units for the RMSE is normalised AV45-PET UPTAKE.

Figure 3 shows long term trajectories estimated with the three models on three sample regions of interest, precuneus, supramarginal and lingual, over the estimated long term time axis (roughly 25 years). The three regions of interest were selected to be representative of, respectively, early, intermediate and late disease progression (times in which regions reach abnormality are shown in Figure 5). We note that while the estimated long term trajectories may appear similar across DS, the underlying dynamic properties are different by construction. Figure 3 also shows individual measurements, colored according to the diagnosis and shifted according to the estimated time reparametrization parameters. These parameters are similar across DS, and produce similar diagnostic separation (Supp. Fig. S9). The long term trajectories of the whole set of regions are in Supp. Figures S6, S7 and S8 (S5 for the GPPM-mono model). We note that most regions show an initial plateau at 10 years, and only few are still accumulating amyloid in the time-span estimated by the models. Plots of the RMSE distributions for the individual predictions of the GPPM-DS models across biomarkers are in Supp. Figures S10, S11 and S12.

GPPM-ACP: kinetic parameters and propagation pathways

For each region i, the GPPM-ACP model estimates saturation (γ i ) and plateau (η i ) thresholds with associated variability (equation (7, Method Section 4.2). The associated times to reach saturation and plateau can be computed as t γ i such that f i (t γ i ) = γ i , and

t η i such that f i (t η i ) = η i .
Figure 4A) shows GP functions along with the estimated γ i and η i for three brain regions (precuneus, supramarginal and lingual). The estimated times t γ i and t η i vary across regions and reflect the respective role in the process of amyloid deposition in AD. For instance, we note how the precuneus region saturates quite early, triggering propagation to connected regions. This result is compatible with previous findings in histo-pathological and imaging studies of amyloid deposition in AD [START_REF] Braak | Neuropathological stageing of Alzheimer-related changes[END_REF][START_REF] Thal | Phases of Ab-deposition in the human brain and its relevance for the development of AD[END_REF][START_REF] Murray | Clinicopathologic and 11C-Pittsburgh compound B implications of Thal amyloid phase across the Alzheimer's disease spectrum[END_REF][START_REF] Thal | Estimation of amyloid distribution by [18 F] flutemetamol PET predicts the neuropathological phase of amyloid beta-protein deposition[END_REF][START_REF] Fantoni | A systematic review and aggregated analysis on the impact of amyloid PET brain imaging on the diagnosis, diagnostic confidence, and management of patients being evaluated for Alzheimer's disease[END_REF][START_REF] Grothe | In vivo staging of regional amyloid deposition[END_REF]. The GPPM-ACP model also estimates k ij ACP , the region-wise time-dependent propagation parameters describing amyloid propagation among connected regions -see equation [START_REF] Fornito | Connectomics: a new paradigm for understanding brain disease[END_REF], Method Section 4.2. With reference to the red/purple/light blue vertical bars in Figure 4A), Figure 4B) shows the connectomes where the edges' colors are proportional to the values of the estimated kinetic parameters. The three regions of interest are highlighted, and we can observe that each region's propagation is triggered when the region saturates, gets stronger while propagation is ongoing, and finally reaches a plateau for the regions that reach a plateau. Figure 4C) shows the same results as in Figure 4B), but the connectome is visualised as an adjacency matrix, where gray/black cells correspond to present/missing connections amongst regions, while colors for the lingual, precuneus and supramarginal connections are proportional to the values of the estimated kinetic parameters. An animated version of Figure 4 is available as Supp. data.

GPPM-ACP: saturation thresholds as a proxy of time-to-abnormality

The estimated kinetic parameters can be interpreted in term of disease progression, for example in terms of ordering in which regions become abnormal. Figure 5A) shows the regional ordering induced by the time to reach abnormality, measured as the time at which maximal separation between CN-and AD subjects was measured for that region. Figure 5B) shows the regional ordering induced by the estimated time to reach saturation t γ i . The two orderings share a pattern compatible with previous findings in amyloid deposition, reporting frontal and parietal areas as the first regions involved in amyloid deposition [START_REF] Mawuenyega | Decreased clearance of CNS b-amyloid in Alzheimer's disease[END_REF][START_REF] Rodrigue | b-Amyloid burden in healthy aging: regional distribution and cognitive consequences[END_REF]. This result highlights that the regions that first reaches saturation and start propagation are also the first one to show abnormality. The last involved regions are the subcortical ones, such as the thalamus [START_REF] Irvine | Protein aggregation in the brain: the molecular basis for Alzheimer's and Parkinson's diseases[END_REF].

Fig. 5: Saturation thresholds as a proxy of time-to-abnormality. A) Regional ordering, according to the time to reach abnormality, measured as the time at which maximal separation between CN-and AD subjects was measured for each region, with associated variability. B) Regional ordering, according to the time to reach saturation t γ i and associated variability -see equation [START_REF] Fornito | Connectomics: a new paradigm for understanding brain disease[END_REF], Method Section 4.2. In both panels, the 3 regions of interest -precuneus, supramarginal and lingual, are highlighted in blue, orange and green, respectively.

GPPM-ACP: personalisation of protein dynamics

For each subject, the DS structure of GPPM-ACP allows to integrate the dynamics over time given an initial condition (the individual baseline measurements) thus obtaining a vector field governing forward and backward evolution in time associated with the individuals. Figure 6A) shows the estimated vector field for the dynamics associated with two regions of interest: precuneus and lingual. The field was obtained by integrating the associated ACP dynamical system while setting the other biomarkers constant to their mean values. We can appreciate the non-linear dynamics of the ACP model. Given the baseline values, we follow the predicted dynamics over time associated with the variability on the kinetic parameters. By plotting the subject's follow-up data of the testing data set D 2 , we note that the model can reliably predict the unseen point. Specifics on the selected subject can be found in Supp. Table S10. Figure 6B) shows the predicted cumulative amyloid deposition in the whole brain. An animated version of Figure 6 is available as Supp. data.

Discussion

We presented a framework for spatio-temporal modeling of protein dynamics over brain networks from short term imaging data, which enables the investigation of bio-mechanical hypotheses governing disease progression via Bayesian model comparison. The framework leverages on Bayesian non-parametric regression techniques, and is coupled with a variational inference approach for scalable inference, in this way providing uncertainty quantification of the kinetic parameters governing protein dynamics and of the long term ND progression. Constraints on the protein dynamics are enforced via DS, naturally allowing prediction of protein changes over time, and so ultimately enabling realistic personalized simulation of pathological evolution. When applied to AV45-PET brain imaging data, our framework provides new insights on the mechanisms of amyloid deposition in AD, indicating the ACP model as the most accurate dynamical system for bio-mechanical interpretation of amyloid dynamics. The model provides plausible simulations of protein propagation, and achieves accurate predictions of individual protein deposition in unseen data.

Our results further indicate that fundamentally different propagation mechanisms can be associated with similar progression patterns. This is a natural consequence of the ill-posed nature of the modeling problem tackled in this work. For this reason, the diagnostic separation along the disease progression axis is mildly affected from the choice of the underlying mechanistic process. This is confirmed by the analysis of the limit case where the model is optimized with respect to the sole monotonic constraint of the trajectories (GPPM-mono). Although not associated with any specific mechanistic hypothesis, the diagnostic separation obtained with this simple model is still very similar.

As implemented in our framework, Bayesian model comparison through the inspection of the ELBO accounts for model complexity. Indeed, highly parameterized models are automatically penalised by the Kullback-Leibler divergence term, which scales with the number of parameters. This is a direct consequence of the fact that the lower bound is a surrogate of the model evidence. To further mitigate the ill-posedness nature of the modelling problem, we also introduced a sparsity constraint through variational dropout, which we have shown to enhance model parsimony by considerably reducing the number of effective parameters. For all these reasons, our theoretical setup provided us with formal and practical guarantees against overfit. As a demonstration of this aspect, we also showed that the ELBO is associated with out of sample predictive accuracy.

We observe that, while we have implemented mechanistic models attempting to explain amyloid spread in humans based on common and ackowledged hypotheses [START_REF] Bateman | Human amyloid-beta synthesis and clearance rates as measured in cerebrospinal fluid in vivo[END_REF][START_REF] Soto | Protein misfolding, aggregation, and conformational strains in neurodegenerative diseases[END_REF], alternative hypotheses have been presented [START_REF] Whittington | Alzheimer's Disease Neuroimaging Initiative, 2018. Spatiotemporal distribution of beta-amyloid in Alzheimer disease is the result of heterogeneous regional carrying capacities[END_REF], describing the distribution of amyloid in Alzheimer's disease as the result of heterogeneous regional carrying capacities.

We note that a limitation in the use of DS concerns the identifiability of the kinetic parameters, whose analysis typically requires symbolic computation [START_REF] Chis | Structural identifiability of systems biology models: a critical comparison of methods[END_REF] and can become prohibitive as the system size and non-linearity increase. In particular, while for linear dynamical systems there are well-established frameworks for establishing either structural, local, or data-driven identifiability, relying on either the power series method [START_REF] Pohjanpalo | System identifiability based on the power series expansion of the solution[END_REF] or the Laplace transform [START_REF] Delbary | Compartmental analysis of dynamic nuclear medicine data: models and identifiability[END_REF], proof of identifiability for general non-linear dynamical systems is still missing [START_REF] Saccomani | Identifiability of Nonlinear ODE Models in Systems Biology: Results from Both Structural and Data-Based Methods[END_REF]. A number of methods has been proposed in recent years [START_REF] Villaverde | Structural identifiability of dynamic systems biology models[END_REF][START_REF] Thomaseth | Local identifiability analysis of nonlinear ODE models: how to determine all candidate solutions[END_REF], for analyzing identifiability of sub-classes of non-linear dynamical systems: to our knowledge however, no method can deal with nonlinear system with time-varying parameters, which is, for instance, the case of the ACP model. Further, the proposed framework uniquely combines the DS with a disease progression models for estimating long term protein trajectories from short term data. This adds a second layer of parameters for which identifiability analysis is required: the individual time parameters and the parameters for the protein trajectories fit. Although we do not provide here formal investigation of identifiability of the GPPM-DS framework, the proposed Bayesian framework naturally allows for uncertainty quantification on the model parameters. Also, we leverage on variational dropout techniques [START_REF] Kingma | Variational dropout and the local reparameterization trick[END_REF][START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF][START_REF] Antelmi | Sparse Multi-Channel Variational Autoencoder for the Joint Analysis of Heterogeneous Data[END_REF] for imposing parsimonious representations of the model parameters through sparsity constraints by means of the sparse GPPM-DS. From a biological perspective, imposing sparsity on the kinetic parameters implies retaining only the brain connections along which protein propagation is strong. Since brain networks have "small-world" organization, segregation and integration are balanced in such a way that connections are somehow redundant -in this way, by retaining only connections along which propagation is strong, we may be discarding weak paths while the flow of propagation could still able to reach the other regions.

We assume here propagation to take place along the structural connectome [START_REF] Iturria-Medina | Epidemic spreading model to characterize misfolded protein propagation in aging and associated neurodegenerative disorders[END_REF][START_REF] Raj | Network diffusion model of progression predicts longitudinal patterns of atrophy and metabolism in Alzheimer's disease[END_REF][START_REF] Warren | Molecular nexopathies: a new paradigm of neurodegenerative disease[END_REF][START_REF] Crossley | The hubs of the human connectome are generally implicated in the anatomy of brain disorders[END_REF][START_REF] Prusiner | A ing role for prions in neurodegenerative diseases[END_REF]. Anatomical connectivity networks are a natural choice for propagation models as they estimate physical connections between brain regions, rather than the correlations estimated in functional [START_REF] Ogawa | Brain magnetic resonance imaging with contrast dependent on blood oxygenation[END_REF] and gray matter structural covariance networks [START_REF] Alexander-Bloch | The convergence of maturational change and structural covariance in human cortical networks[END_REF]. To estimate the structural connectomes we rely on tractography, which is prone to false positive and negative connections [START_REF] Maier-Hein | The challenge of mapping the human connectome based on diffusion tractography[END_REF][START_REF] Thomas | Anatomical accuracy of brain connections derived from diffusion MRI tractography is inherently limited[END_REF]. Nevertheless, here we take an average connectome over multiple young and healthy subjects, which we believe can provide an accurate anatomical reference, although not accounting for the effects of the disease on the connectome itself [START_REF] Oxtoby | Data-Driven sequence of changes to anatomical Brain connectivity in sporadic Alzheimer's Disease[END_REF]. Future models could approximate connectivity impairment arising from white matter damage commensurate with amyloid deposition. Further, limitations of current diffusion MRI tractography techniques makes it impossible to consider non symmetrical connectivity matrices, and as a consequence anterograde or retrograde protein propagation processes, which could potentially present different dynamics. Nonetheless, one major advantage of our framework is that it naturally allows the use of different or complementary networks that provide directionality, as fMRI or EEG, by introducing a non-symmetrical graph adjacency matrix A and a non-symmetrical propagation matrix. K θ . The caveat is that this would increase the total number of parameters to be estimated by the model. Finally, we note that the chosen ADNI dataset is suboptimal for longitudinal amyloid data analysis, as longitudinal amyloid measurement is more stable when using a white matter reference region [START_REF] Chen | Improved power for characterizing longitudinal amyloid-beta PET changes and evaluating amyloid-modifying treatments with a cerebral white matter reference region[END_REF][START_REF] Landau | Measurement of longitudinal beta-amyloid change with 18F-florbetapir PET and standardized uptake value ratios[END_REF][START_REF] Brendel | Improved longitudinal [18F]-AV45 amyloid PET by white matter reference and VOI-based partial volume effect correction[END_REF] and using a longitudinal pipeline [START_REF] Reuter | Avoiding asymmetry-induced bias in longitudinal image processing[END_REF]. Similarly, the chosen dataset precludes an investigation as to how resolution (i.e. number of ROIs)/atlas affects model fitting.

The ideas we propose here may apply to a much larger range of neurological diseases with proteinopathies and alternative models of propagation. The proposed framework will be used in the future to investigate a model of the natural history of AD, integrating a comprehensive panel of biomarkers: imaging, cognitive, genetic and demographic data. To this end, we aim at developing DS jointly accounting for biomarkers of protein dynamics (including AV1451-PET for Tau accumulation in AD), hypometabolism and atrophy [START_REF] Abi-Nader | Monotonic Gaussian Process for spatio-temporal disease progression modeling in brain imaging data[END_REF]. This extension would require the definition of appropriate mechanistic constraints within-and across-different protein spreading processes. Of course, the experimental validation of this extension requires the availability of sufficient longitudinal samples. Also, the GPPM framework rests on the assumption of consistent regional propagation across individuals, which is the case for amyloid protein [START_REF] Grothe | In vivo staging of regional amyloid deposition[END_REF], while tau as well as MRI-derived atrophy patterns are far more heterogeneous [START_REF] Tam | A highly predictive signature of cognition and brain atrophy for progression to Alzheimer's dementia[END_REF][START_REF] Murray | Neuropathologically defined subtypes of Alzheimer's disease with distinct clinical characteristics: a retrospective study[END_REF]. Future work could be devoted to study protein kinetics on more homogeneous subgroups or to explore the variability among different subgroups.

Finally, the application of the model in clinical trials, although of certain interest, falls beyond the scope of the current study. A future validation of this framework with clinical partners is currently under study.

GPPM-DS availability

GPPM-DS source code will be available upon acceptance. Data used in the preparation of this work are from ADNI. ADNI is a public-private partnership. All ADNI data are shared without embargo through the LONI Image and Data Archive (https://ida.loni.usc.edu/login.jsp) a secure research data repository. Interested scientists may obtain access to ADNI imaging, clinical, genomic, and biomarker data for the purposes of scientific investigation, teaching, or planning clinical research studies. Access is contingent on adherence to the ADNI Data Use Agreement. For up-to-date information please see http://adni.loni.usc.edu/wp-content/uploads/howto-apply/ADNI-DSP-Policy.pdf.

Fig. 1 :

 1 Fig. 1: Examples of Diff, RD and ACP models with varying key parameters on a toy data set of N = 3 synthetic brain regions. A) Diff model with varying rates of propagation k ij . For each biomarker, k ij is a random realisation from a Gaussian centered in (from left to right): 0.05, 0.5 and 1. B) RD model with varying rates of propagation k ij . For each biomarker, k ij is a random realisation from a Gaussian centered in (from left to right): 0.05, 0.5 and 1. The threshold parameter υ is set at 1. C) RD model with varying threshold parameter υ. For each biomarker, υ is a random realisation from a Gaussian centered in (from left to right): 0.5, 1 and 3. The k ij parameters are set at 0.5. D) ACP model with varying rates of propagation k ij . For each biomarker, k ij is a random realisation from a Gaussian centered in (from left to right): 0.05, 0.5 and 1.The threshold parameters γ and η ares set at 0.6 and 0.9 respectively. The sigmoid parameters l 1 and l 2 are set at 3. E) ACP model with varying sigmoid parameters l 1 and l 2 . For each biomarker, l 1 and l 2 are a random realisation from a Gaussian centered in (from left to right): (0.5, 0.5),[START_REF] Sweeney | Protein misfolding in neurodegenerative diseases: implications and strategies[END_REF][START_REF] Sweeney | Protein misfolding in neurodegenerative diseases: implications and strategies[END_REF] and[START_REF] Jucker | Self-propagation of pathogenic protein aggregates in neurodegenerative diseases[END_REF][START_REF] Jucker | Self-propagation of pathogenic protein aggregates in neurodegenerative diseases[END_REF]. The k ij parameters are set at 0.5. The threshold parameters γ and η are set at (0.6, 0.9). F) ACP model with varying threshold parameters γ and η. For each biomarker, γ and η are a random realisation from a Gaussian centered in (from left to right): (0.1, 0.3), (0.6, 0.9) and (2, 3). The k ij parameters are set at 0.5. The sigmoid parameters l 1 and l 2 are set at 3.

Fig. 2 :

 2 Fig. 2: Schematic representation of the proposed GPPM-DS framework. Regional protein concentrations f i and f j are collected for k subjects over a short term time span t, encoded in a measurement array Y k (t) (A). The dynamics of such concentrations are described in terms of a functional H θ , with unknown parameters θ, encoded in a dynamical system for the vector of concentrations f (B). The proposed framework estimates such parameters encoding the strength of propagation (D) and the long term protein concentrations with respect to the estimated long term time axis (C).

Fig. 3 :

 3 Fig. 3: Modeling amyloid deposition in Alzheimer's disease. Long term trajectories and individual short term measurements estimated with the three GPPM-DS models on three sample regions of interest: precuneus, supramarginal and lingual, selected to be representative of early, intermediate and late progression. Specifically: in A) results for GPPM-Diff of equation (4); B) results for GPPM-RD of equation (5); C) results for GPPM-ACP of equation (6). The black solid line represents the average trajectory, while the red dashed ones ± 3 standard deviations. In blue CN-subjects, in orange CN+ subjects, in green MCI subjects and in red AD subject.

Fig. 4 :

 4 Fig. 4: Kinetic parameters and propagation pathway. A) Long term trajectories estimated by GPPM-ACP; mean saturation threshold parameters (horizontal dashed lines); mean plateau threshold parameters (horizontal dotted lines) and corresponding distributions on the y-axis; for 3 regions of interest: precuneus (in blue), supramarginal (in orange) and lingual (in green). B) Estimated time-dependent propagation parameters to/from the 3 regions of interest along the associated anatomical connections, sampled at 3 times (corresponding to the red/purple/light blue vertical bars shown in A). The colors of the edges of the connectome are proportional to the values of the estimated kinetic parameters k ij ACP , which is the region-wise time-dependent propagation parameters describing amyloid propagation among connected regions. Glass brain images obtained with Nilearn [65] (available at https://nilearn.github.io/). C) Estimated time-dependent propagation parameters to/from the 3 regions of interest overlaid on the binary adjacency matrix of the connectome, sampled at 3 times (corresponding to the red/purple/light blue vertical bars shown in A). The adjacency matrix shows black/gray entries for missing/present connections amongst regions. The colors of the entries (i, j) describing connections to/from the three regions of interest are proportional to the values of k ij ACP .

Fig. 6 :

 6 Fig. 6: Personalisation of protein dynamics. For each subject, the DS structure of GPPM-ACP allows to integrate the dynamics over time given an initial condition (the individual baseline measurements), thus obtaining a vector field governing forward and backward evolution in time associated with the individuals. A) Streamlines for the estimated amyloid deposition dynamics of precuneus and lingual regions for a sample healthy individual converted to MCI at 8 years from baseline, and subsequently to AD at 10 years. Black lines: predicted dynamics. Red dashed lines: associated variability. Colored dots: 4 observed time-points at 4, 6, 8, 10 years. Star shaped point: unseen follow-up at 11 years; B) Predicted cumulative amyloid deposition for the time-points highlighted in A). Brain images obtained with BrainPainter [55] (available at brainpainter.csail.mit.edu/).

Acknowledgements

SG acknowledges financial support from the French government managed by L'Agence Nationale de la Recherche under Investissements d'Avenir UCA JEDI (ANR-15-IDEX-01) through the project "AtroProDem: A data-driven model of mechanistic brain Atrophy Propagation in Dementia". ML was supported by the French government, through the UCA JEDI and 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR-19-P3IA-0002, ANR-19-CE45-0006). This work is supported by the Inria Sophia Antipolis -Méditerranée, "NEF" computation cluster. Data collection and sharing for this project was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer's Association; Alzheimer's Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research and Development, LLC.; Johnson and Johnson, Pharmaceutical Research and Development LLC.; Lumosity; Lundbeck; Merck and Co., Inc.; Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research is providing fuND to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer's Disease Cooperative Study at the University of California, San Diego. ADNI data are disseminated by the Laboratory for NeuroImaging at the University of Southern California. The HCP project (Principal Investigators: Bruce Rosen, M.D., Ph.D., Martinos Center at Massachusetts General Hospital; Arthur W. Toga, Ph.D., University of California, Los Angeles, Van J. Weeden, MD, Martinos Center at Massachusetts General Hospital) is supported by the National Institute of Dental and Cranio-facial Research (NIDCR), the National Institute of Mental Health (NIMH) and the National Institute of Neurological Disorders and Stroke (NIND). Collectively, the HCP is the result of efforts of co-investigators from the University of California, Los Angeles, Martinos Center for Biomedical Imaging at Massachusetts General Hospital (MGH), Washington University, and the University of Minnesota.

Competing interests

The authors declare no competing interests.

Authors contributions

SG and ML conceived the study and designed the experiments. SG performed the analysis that produced the connectomes, collated the data and performed the experiments with inputs from ML. SG and ML wrote the manuscript and provided the interpretation of the results.

Competing interests

The authors declare no competing interests.