Alexandre Bône

Olivier Colliot

Stanley Durrleman

Learning the spatiotemporal variability in longitudinal shape data sets

Keywords: Longitudinal data, Statistical shape analysis, Large deformation diffeomorphic metric mapping, Medical imaging, Disease progression modeling

In this paper, we propose a generative statistical model to learn the spatiotemporal variability in longitudinal shape data sets, which contain repeated observations of a set of objects or individuals over time. From all the short-term sequences of individual data, the method estimates a longterm normative scenario of shape changes and a tubular coordinate system around this trajectory. Each individual data sequence is therefore (i) mapped onto a specific portion of the trajectory accounting for differences in pace of progression across individuals, and (ii) shifted in the shape space to account for intrinsic shape differences across individuals that are independent of the progression of the observed process. The parameters of the model are estimated using a stochastic approximation of the expectation-maximization algorithm. The proposed approach is validated on a simulated data set, illustrated on the analysis of facial expression in video sequences, and applied to the modeling of the progressive atrophy of the hippocampus in Alzheimer's disease patients. These experiments show that one can use the method to reconstruct data at the precision of the noise, to highlight significant factors that may modulate the progression, and to simulate entirely synthetic longitudinal data sets reproducing the variability of the observed process.

Introduction

Motivation

Video sequences of smiling faces, repeated measurements of growing plants or developing cells, medical images collected at multiple visits from a population of patients affected by a chronic disease: all these examples can be understood as data collections where individual instances of a common underlying process are observed at multiple timepoints. Such collections are called longitudinal data sets.

The individual processes are thought to result from random variations of a common underlying process (or few of them). Because of the dynamic nature of the observed processes, one might decompose the variability in two components: the dynamic or temporal variability on the one hand, and the time-independent or spatial variability on the other hand. In our examples, the variability in the pace of growth or in age at disease onset is understood as temporal variability. By contrast, there are also intrinsic inter-individual differences in height, weight or shape that are independent of the pace at which the plant grows or the disease progresses, which we call spatial variability. The main difficulty here is that growth or disease progression affect also height, weight or shape, so that the differences between two observations of two different samples are due to (i) the fact that the two individuals are observed at different stages of the process, and (ii) that they have different intrinsic characteristics. Disentangling these two sources of variability would not be possible if one had only one observation per individual. Having repeated observations of the individuals over time, as in longitudinal data sets, implies that one could decompose the changes due to the progression of the process from those due to intrinsic differences that are independent of the progression.

The goal of this paper is to propose a statistical learning method which can describe the spatiotemporal variability in a longitudinal data set. We focus here on shape data, where the shape may be encoded by an image, or by geometrical objects extracted from images such as curves, surface meshes or segmented volumes.

One of the main difficulties is that the experimental design often provides little control on the temporal sampling of the observations. We are interested here by processes for which there is no clear marker of progression, such as the progression of neurodegenerative diseases for which the age at disease onset is hard to determine. Therefore there is no easy way to re-align in time the individual data sequences to analyze the inter-individual variability at each stage of the process. By contrast, the method needs to learn how the individual data sequences position themselves in relation to each other. Furthermore, the follow-up period of the observations rarely covers the whole process, but often just a small part of it. In clinical studies for examples, patients may be followed for few years whereas the disease may progress over decades. Eventually, dealing with shape data raises the need for generic representations of such data that can be included in computational approaches.

Related work

Structured data like shapes can be advantageously represented as elements of curved spaces, such as Riemmanian manifolds, in order to account for the prior on their structure. Either defined by invariance [START_REF] Kendall | Shape manifolds, procrustean metrics, and complex projective spaces[END_REF][START_REF] Su | Statistical analysis of trajectories on riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF][START_REF] Su | Rate-invariant analysis of trajectories on riemannian manifolds with application in visual speech recognition[END_REF] or topologypreserving properties [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF][START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data[END_REF][START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF][START_REF] Pennec | A riemannian framework for tensor computing[END_REF], shape spaces define distance metrics adapted to the geometry of a well-identified class of objects, such as brain magnetic resonance images or segmented organs. These data representations allows the generalization of the mean-variance analysis [START_REF] Allassonnière | Bayesian mixed effect atlas estimation with a diffeomorphic deformation model[END_REF][START_REF] Gori | A Bayesian Framework for Joint Morphometry of Surface and Curve meshes in Multi-Object Complexes[END_REF][START_REF] Pennec | Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements[END_REF][START_REF] Zhang | Bayesian estimation of regularization and atlas building in diffeomorphic image registration[END_REF], which learns the geometrical distribution of a cross-sectional data set in terms of an average shape, and variability-encoding parameters. Typical healthy or pathological configurations can be summarized in this manner, thus opening the way to automatic diagnosis at the individual level. Time-series data sets, consisting in the repeated observation of the same object at successive time-points, can be described by generalized regression approaches on the same shape spaces [START_REF] Banerjee | A nonlinear regression technique for manifold valued data with applications to medical image analysis[END_REF][START_REF] Fishbaugh | Geodesic regression of image and shape data for improved modeling of 4D trajectories[END_REF][START_REF] Fletcher | Geodesic regression and the theory of least squares on riemannian manifolds[END_REF][START_REF] Hinkle | Polynomial regression on riemannian manifolds[END_REF][START_REF] Lorenzi | 4D registration of serial brain's MR images: a robust measure of changes applied to Alzheimer's disease[END_REF][START_REF] Niethammer | Geodesic regression for image time-series[END_REF]. A time-continuous scenario of geometrical transformation is then estimated, offering in turn individualized interpolation and extrapolation methods. The statistical analysis of longitudinal data sets requires to extend the concept of generalized mean-variance for such time series. In other words, it requires the definition of a statistical distribution of curves drawn on a shape space.

Shape spaces are usually equipped with a differential structure of infinite dimension. In particular, the large deformation diffeomorphic metric mapping (LDDMM) approach defines shape spaces as orbits of template shapes under the action of an infinite-dimensional parametric group of diffeomorphisms of the 2D/3D ambient space [START_REF] Younes | Shapes and Diffeomorphisms[END_REF]. With this approach, the geometrical differences between two objects are captured by estimating the diffeomorphic transformation that warps one into the other. More recent works propose finite-dimensional approaches built on the same principles: [START_REF] Zhang | Finite-dimensional lie algebras for fast diffeomorphic image registration[END_REF] uses truncated Fourier transforms to build a finitedimensional Lie algebra, and [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF] constructs a finite-dimensional Riemannian manifold based on a set of self-interacting particles.

Such structures are favorable to the analysis of longitudinal data sets because they naturally offer the parallel transport operator [START_REF] Lorenzi | Schild's ladder for the parallel transport of deformations in time series of images[END_REF], which allows to compare tangent-space vectors at distant points in a relevant manner. This operator is key to compare trajectories on the manifold, and therefore to analyze longitudinal data. In [START_REF] Su | Statistical analysis of trajectories on riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF][START_REF] Su | Rate-invariant analysis of trajectories on riemannian manifolds with application in visual speech recognition[END_REF] for instance, trajectories on manifolds are compared by parallel-transporting their initial velocity vectors back to some privileged point of the manifold, thereby handling the spatial variability if a reference configuration and reference time-point is known. A similar approach is followed in [START_REF] Kim | Riemannian nonlinear mixed effects models: Analyzing longitudinal deformations in neuroimaging[END_REF] where medical images are analyzed in a voxel-wise fashion, or also in [START_REF] Singh | Hierarchical geodesic models in diffeomorphisms[END_REF] with the co-adjoint transport instead of the parallel transport. In [START_REF] Cury | Spatiotemporal analysis for detection of presymptomatic shape changes in neurodegenerative diseases: Initial application to the genfi cohort[END_REF], the variability of a large longitudinal data set of thalamus shapes is analyzed by transporting individual residual deformations along a common and pre-computed trajectory, back to a baseline point. In [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF][START_REF] Schiratti | A bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations[END_REF] the authors define the expparallelization operator which extends the notion of parallel lines to Riemannian manifolds. The works [START_REF] Bône | Learning distributions of shape trajectories from longitudinal datasets: a hierarchical model on a manifold of diffeomorphisms[END_REF][START_REF] Koval | Statistical learning of spatiotemporal patterns from longitudinal manifoldvalued networks[END_REF] build on this operator to analyze dynamic networks and shape objects respectively. Other approaches propose to work on a space of trajectories, such as in [START_REF] Muralidharan | Sasaki metrics for analysis of longitudinal data on manifolds[END_REF] where the Sasaki metric is used to define distances between geodesic curves on a manifold, or in [START_REF] Chakraborty | Statistics on the space of trajectories for longitudinal data analysis[END_REF] which requires the same number of observations per subject.

If parallel transport allows to spatially align manifoldvalued trajectories, a temporal alignment mechanism is also needed for data sets with variability in the individual progression dynamics. For instance, two patients developing the same neurological disease have no reason to reach the same disease stage at the same age, nor to have synchronous progressions. A solution is to use time-warp functions, which define a mapping between an abstract common reference time frame and the individual time lines [START_REF] Bilgel | A multivariate nonlinear mixed effects model for longitudinal image analysis: Application to amyloid imaging[END_REF][START_REF] Bône | Learning distributions of shape trajectories from longitudinal datasets: a hierarchical model on a manifold of diffeomorphisms[END_REF][START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data[END_REF][START_REF] Kim | Riemannian nonlinear mixed effects models: Analyzing longitudinal deformations in neuroimaging[END_REF][START_REF] Koval | Statistical learning of spatiotemporal patterns from longitudinal manifoldvalued networks[END_REF][START_REF] Marinescu | A vertex clustering model for disease progression: application to cortical thickness images[END_REF][START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF][START_REF] Schiratti | A bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations[END_REF]. In [START_REF] Su | Statistical analysis of trajectories on riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF][START_REF] Su | Rate-invariant analysis of trajectories on riemannian manifolds with application in visual speech recognition[END_REF], the authors build on the square-root velocity fields framework to quotient the space of spatiotemporal paths by diffeomorphic time-warps. In [START_REF] Nader | Monotonic gaussian process for spatio-temporal trajectory separation in brain imaging data[END_REF], a monotonic Gaussian process is built from a set of temporal sources.

Contributions

In this paper, we propose a method that learns an average progression and its spatiotemporal variability from a longitudinal shape data set. The average progression takes the form of a geodesic curve in the finite-dimensional Riemannian approximation of the LDDMM framework of [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. The concept of exp-parallelization introduced in [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF][START_REF] Schiratti | A bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations[END_REF] is then applied in this context to define a tubular coordinate system, also called Fermi coordinates, around the average geodesic. The average trajectory and its coordinate system are automatically learned by the method, so that every individual data sequence is mapped to a specific portion of the average trajectory to account for the temporal variability, and shifted in the shape space to account for the spatial variability. The calibration of the resulting generative statistical model is done by adapting a stochastic approximation EM method. This paper extends the conference paper [START_REF] Bône | Learning distributions of shape trajectories from longitudinal datasets: a hierarchical model on a manifold of diffeomorphisms[END_REF], with finer modeling of the variability in the individual paces of progression, and an original optimization method for accelerated model calibration. The proposed approach is validated on a simulated data set, illustrated on a facial expression recognition task, and applied to hippocampus shape progression modeling in Alzheimer's disease.

Section 2 defines the concept of shape spatiotemporal coordinate systems, which allows the introduction of the generative statistical model in Section 3. Section 4 details the calibration, personalization and simulation algorithms, which are evaluated and illustrated in Section 5. These experiments will evaluate the goodness-of-fit of the model, the relevance of the representation of the spatiotemporal variability for the identification of factors explaining this variability, and the ability of the model to generate synthetic data sets that reproduce the observed variability in the training data set.

Shape spatiotemporal reference frame

Within LDDMM frameworks, shape are positioned with respect to a reference shape, often called atlas or template. A coordinate system is defined in the tangent-space at the template shape. We propose here to replace the template (which is a single shape) by a curve (i.e. a shape trajectory), and the coordinate system by a tubular spatiotemporal coordinate system centered around the template trajectory. We first review the usual construction of a static template shape before extending it to the spatiotemporal case.

2.1 Positioning a shape with respect to a static atlas Positioning a target shape y with respect to a static reference y 0 is called the registration problem. Deformation-based mor-phometry solves it by estimating a diffeomorphism φ 1 of the ambient space R d (d = 2 or 3) that transforms y 0 into y, which we note φ 1 y 0 = y. In the context of LDDMM, diffeomorphisms are constructed by following the streamlines of dynamic vector fields t

→ v t ∈ C ∞ 0 (R d , R d ) over [0, 1]: ∂ t φ t = v t • φ t with φ 0 = Id. ( 1
)
Following the approach in finite-dimension of [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF], we further assume that any v t writes as the Gaussian convolution of p momentum vectors m t = m

(1) t , ..., m

(p) t ∈ R d over a corresponding set of control points c t = c (1) t , ..., c (p) t ∈ R d : v t : x ∈ R d → p k=1 g c (k) t , x • m (k) t ∈ R d (2) with g : x, x ∈ R d → exp x -x 2
2 /σ 2 the Gaussian kernel function of kernel width σ > 0. Many other diffeomorphisms constructed in this manner might actually transform y 0 into φ 1 y 0 : we call solution of the registration problem the most regular transformation i.e. that minimizes its "kinetic" energy:

1 2 1 t=0 v t 2 Gc t = 1 2 1 t=0 m t • G ct • m t (3)
where ∀t ∈ [0, 1], G ct is the p × p "kernel" symmetric positive-definite matrix of general term g[c

(k) t , c (l)
t ], and (.) is the matrix transposition. Such energy-minimizing curves, also called geodesics, are such that the control points and the momentum vectors trajectories are fully determined by their initial values and the following Hamiltonian equations [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF]:

ċt = G ct • m t ; ṁt = - 1 2 ∇ ct m t • G ct • m t (4)
where ∇ x (.) is the gradient operator with respect to x. Assuming that there exist a diffeomorphism φ 1 constructed according to equations (1, 2) such that φ 1 y 0 = y, this last result allows to compactly represent the positioning of y with respect to y 0 with a set of p control points c 0 and attached momenta m 0 . In other words, m 0 is the coordinate of y in the coordinate system defined by (c 0 , y 0 ). In practice, the perfect registration constraint φ 1 y 0 = y is relaxed, and we call solution to the registration problem the extremal-path diffeomorphism φ 1 that warps y 0 as close as possible to y, for some extrinsic error measure d E (y 0 , y). In this paper, the following choices are considered for d E , depending on the nature of y 0 and y:

the 2 metric for meshes with point-to-point correspondence (i.e. the sum of squared differences between point positions), the current metric [START_REF] Charon | Fidelity metrics between curves and surfaces: currents, varifolds, and normal cycles[END_REF][START_REF] Vaillant | Surface matching via currents[END_REF] for oriented surface meshes without point-to-point correspondence (details are given in appendix for the reader's convenience).

Noting y 0 as a collection y

(1) 0 , ..., y (K) 0

of K points of R d , φ 1 acts independently and directly on each point y

(k) 0 ac- cording to φ 1 y (k) 0 = φ 1 • y (k)
0 . Note that the methodology introduced in this section can be adapted straightforwardly to image data, by defining the action of the diffeomorphisms φ on the image I as I • φ -1 and using the sum of squared differences between image intensities as the error measure.

Riemannian structure

Let c 0 be a set of p control points. We define:

D c0 = φ 1 | ∂ t φ t = v t • φ t , φ 0 = Id, v t = Conv(c t , m t ) ( ċt , ṁt ) = Ham(c t , m t ), m 0 ∈ R p×d (5)
where Conv(., .) and Ham(., .) are compact notations for the convolution operator defined by equation ( 2) and the Hamiltonian equations ( 4 . Shape spatiotemporal reference frame y 0 , c 0 , m 0 , t 0 with respect to which a shape y admits coordinates t ∈ R, v ∈ v ⊥ 0 ⊂ T γ(t0) D c0 . Three spaces are involved: the manifold of control points R p×d (top), the manifold of diffeomorphisms D c0 (middle), and the shape submanifold S y0,c0 of the extrinsic shape space E (bottom). The momenta m 0 , m and the velocity fields v 0 , v are in one-to-one correspondence. The velocity field v, also called space-shift, is parallel-transported along the geodesic γ by the operator t → P v γ (t). Figure 2 illustrates the effect of parallel transport on D c0 . [Right]. Illustrations of the manifolds abstractly depicted on the left side of the figure. Each row displays two example elements of the corresponding geodesic (solid black lines on the left panel). The two columns correspond respectively to the times t 0 and t.

! "

Fig. 2: [Bottom]. Illustration of a shape geodesic t → γ(t) y 0 : the man-like shape (solid black contour) raises his left arm. This geodesic is parametrized by a single set of control points c 0 (black dots) and attached momentum vectors m 0 (bold blue arrows), to which corresponds the velocity field v 0 (light blue arrows). A second set of momentum vectors m (bold red arrows) attached to the same control points c 0 parametrizes the exp-parallelization of this shape geodesic.

[Top]. Exp-parallel shape curve t → η(t) y 0 to the shape geodesic γ y 0 : the exp-parallelization transfers the arm-raising motion from one man-like shape to another.

with the local metric G -1 φ(c0) , D c0 has the structure of a Riemannian manifold of dimension p × d. The tangent-space at φ is the set of velocity fields obtained by convolving any set of momentum on φ(c 0 ):

T φ D c0 = Conv(φ(c 0 ), m) | m ∈ R p×d . ( 6
)
The geodesics of D c0 are the curves t → φ t of constant kinetic energy (see equation (3)) i.e. such that the corresponding control points and momenta trajectories t → c t , m t satisfy the Hamiltonian equations (4). We define the exponential operator on D c0 :

Exp t0,t φ : v 0 ∈ T φ D c0 → φ t ∈ D c0 (7)
where φ t is the diffeomorphism reached at time t by the geodesic path obtained by integration from some reference time t 0 ∈ R with initial conditions φ(c 0 ), m 0 such that v 0 = Conv(φ(c 0 ), m 0 ), and φ 0 = φ. The momentum vector m 0 is the dual of the velocity field v 0 . The particular case Exp 0,1 φ corresponds to the usual Riemannian exponential map and will be noted Exp φ . Diffeomorphisms φ ∈ D c0 act on shapes of the ambient space y through the action previously defined. Let y 0 be a reference shape. We define its orbit under the action :

S y0,c0 = D c0 y 0 = {φ y 0 | φ ∈ D c0 }. (8)
S y0,c0 is a submanifold of the extrinsic shape space E in which is defined the distance d E .

2.3 Positioning a shape with respect to a dynamic atlas Instead of positioning shapes with respect to a static atlas y 0 , we aim now to position shapes with respect to a shape geodesic t → γ(t) y 0 , where γ is a geodesic of D c0 of the form γ(t) = Exp t0,t Id (v 0 ) with v 0 = Conv(c 0 , m 0 ). Similarly to cylindrical coordinates in Euclidian spaces, under some conditions (see [START_REF] Hirsch | Differential topology[END_REF][START_REF] Manasse | Fermi normal coordinates and some basic concepts in differential geometry[END_REF]) a shape y ∈ S y0,c0 admits a unique spatiotemporal coordinate, also known as Fermi coordinates, t ∈ R and v ∈ T γ(t0) D c0 such that v⊥ γ(t 0 ):

y = ExpP v γ (t) y 0 with ExpP v γ (t) = Exp γ(t) P v γ (t) (9)
where P v γ (t) denotes the parallel transport of v along γ from t 0 to t. The curves γ and η : t → ExpP v γ (t) are said expparallel, and the mapping γ → η is called exp-parallelization along v [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF][START_REF] Schiratti | A bayesian mixed-effects model to learn trajectories of changes from repeated manifold-valued observations[END_REF]. In other words, a choice of y 0 , c 0 , m 0 , t 0 defines a spatiotemporal reference frame, with respect to which a shape y can be unambiguously positioned in terms of a time t and a velocity field v orthogonal to v 0 = γ(t 0 ). The time t is the temporal component of the coordinate which positions the shape along the reference trajectory given by the direction v 0 . The velocity v is the spatial component of the coordinate, which positions the shape in the hyperplane that is orthogonal to v 0 . This decomposition can also be understood as the orthogonal projection of y onto the onedimensional submanifold γ y 0 , hence the condition v⊥v 0 .

Figures 1 and2 illustrate this concept of spatiotemporal reference frame, in which any shape y admits the coordinates t, v. Note that the time-point t 0 does not play any particular role, in the sense that y can be described in the same manner for any other choice t 0 ; a one-to-one transformation of the spatiotemporal reference frame can actually be derived as:

t = t + t 0 -t 0 and v = P v γ (t 0 ). (10)
In general, the target shape y might not exactly belong to S y0,c0 . Similarly to the static atlas case, equation ( 9) is relaxed and we call solution to the longitudinal registration problem the pair t, v such that y 0 is warped as close as possible to y, in the sense of the extrinsic metric d E .

3 Statistical model for longitudinal data sets of shapes

Hierarchical generative model

Let {y i,j , t i,j } i,j be a longitudinal data set of shapes, which is the collection of repeated individual measurements y i,1 , ..., y i,ni for i = 1, ..., n, where each shape y i,j corresponds to a time t i,j ∈ R. Measurements are considered as sample points along individual trajectories, which are in turn considered exp-parallel to a reference geodesic curve, therefore having a constant spatial coordinate in the spatiotemporal reference frame centered around this reference geodesic. Noting y 0 , c 0 , m 0 , t 0 the parameters of the spatiotemporal coordinate system, v 0 = Conv(c 0 , m 0 ) and γ : t → Exp t0,t c0 (v 0 ) the reference geodesic, the statistical model writes:

ExpP vi γ ψ i (t i,j ) y 0 iid ∼ N E y i,j , σ 2 ),
where

ψ i : t → α i • (t -τ i ) + t 0 , v i = Conv(c 0 , m i ), m i = A 0,m ⊥ 0 • s i , and
α i iid ∼ N [0,+∞[ (1, σ 2 α ), τ i iid ∼ N (t 0 , σ 2 τ ), s i iid ∼ N (0, 1) (11)
where the noise distribution N E µ, σ 2 ) is defined such that the likelihood is proportional to p(y) ∝ exp(-d E (y-µ) 2 /2σ 2 ). Model [START_REF] Bône | Learning low-dimensional representations of shape data sets with diffeomorphic autoencoders[END_REF] is hierarchical in the sense that individual trajectories t → ExpP vi γ • ψ i (t) are independently defined as spatiotemporal transformations of a common, population-level geodesic t → γ(t).

The time-warp functions ψ i encode the temporal variability of the observed individual trajectories in terms of pace of progression α i and onset time τ i . They map the index t i,j of the j-th shape of the i-th individual (e.g. the age of the subject at a given visit), to a time-point ψ i (t i,j ) on the reference geodesic (e.g. the disease stage).

The spatial variability is encoded by the space-shifts v i ∈ v ⊥ 0 ⊂ T γ(t0) D c0 along which γ is exp-parallelized. Those space-shifts admit dual representations under the form of the momenta m i , which are assumed to derive from q source parameters s i = s

(1) i , ..., s (q) i , in the spirit of independent component analysis (ICA) [START_REF] Hyvärinen | Independent component analysis[END_REF]. The orthogonality v i ⊥v 0 , necessary for the identifiability of the model, is ensured by the projection of each column of the (p • d) × q mixing matrix A 0 onto the hyperplane m ⊥ 0 of R p×d for the cometric G c0 , noted A 0,m ⊥ 0 . The individual parameters are modeled as independent samples from normal distributions:

a truncated normal distribution with fixed mean for the acceleration factor α i , allowing the identifiability of m 0 ; a normal distribution for the onset time τ i ; a normal distribution with fixed mean and variance for the sources s i , allowing the identifiability of y 0 and A 0,m ⊥ 0 respectively.

These parameters define individual trajectories as random spatiotemporal transformations of a common reference trajectory. The spatial and temporal transformations commute, in the sense that ∀t ∈ R, ExpP vi γ•ψi = ExpP vi γ • ψ i . The population trajectory is fully parameterized by the template shape y 0 , the control points c 0 , the momenta m 0 and the reference time t 0 . The individual variability is unambiguously represented by two reduced sets of scalar parameters: the acceleration α i and the onset time τ i for the temporal part, and the sources s

(1) i , ..., s (q)
i for the spatial part. In practice, it is possible to choose a number of sources q p × d much lower than the dimension of the tangent-space T γ(t0) D c0 while still capturing most of the geometrical variability in the data.

Mixed-effects and Bayesian modeling

We further specify the formulation of the model [START_REF] Bône | Learning low-dimensional representations of shape data sets with diffeomorphic autoencoders[END_REF] to fit the framework of mixed-effects models. We distinguish:

-the fixed-effects θ = (θ 1 , θ 2 ) with θ 1 = (t 0 , σ τ , σ α , σ )
and θ 2 = (y 0 , c 0 , m 0 , A 0 ), also called the model parameters, the random effects z = (z i ) i where z i = (α i , τ i , s i ).

We choose to work in a Bayesian framework, in order to theoretically ensure the existence of the maximum a posteriori (MAP) estimate of the parameters θ m . Such priors also regularize and guide the estimation procedure thanks to reasonable and mild prior assumptions on the optimal fixed effects values. The following standard conjugate distributions are selected as Bayesian priors on the model parameters:

t 0 ∼ N (t 0 , ς 2 t ), y 0 ∼ N (y 0 , ς 2 y ), σ 2 τ ∼ IG(m τ , ς 2 τ ), c 0 ∼ N (c 0 , ς 2 c ), σ 2 α ∼ IG(m α , ς 2 α ), m 0 ∼ N (m 0 , ς 2 m ), σ 2 ∼ IG(m , ς 2 ), A 0 ∼ N (A 0 , ς 2 A ),
where IG(., .) denotes the inverse-gamma distribution.

4 Algorithms: calibration, personalization, simulation

Objectives

Given a longitudinal data set of shapes {y i,j , t i,j } i,j that we may note more compactly {y, t}, we formulate three algorithmic objectives:

-Calibration, which consists in computing the MAP parameters θ m , unconditionally to any random effect z:

θ m = argmax θ p {y}, z, θ ; {t} • dz. ( 12
)
-Personalization, which consist in computing the MAP random effects z m that best represent some longitudinal shape data set {y, t} (which may or may not be the one used for calibration), given the calibrated model θ m :

z m = argmax z p {y}, z, θ m ; {t} . (13)
-Simulation, which consist in generating a new data set {y s } that resembles the original data set {y}.

We give now the details of the algorithms to solve these optimization problems. Their implementation is freely available in the software Deformetrica (find the install instructions and the documentation at www.deformetrica.org).

Computation of the complete log-likelihood

Evaluating the joint log-likelihood log p({y}, z, θ ; {t}) = n i=1 ni j=1 log p(y i,j , z i , θ ; t i,j ) for some set of parameters θ and random effects z = (z i ) i is central for both calibration and personalization algorithms. The computationally most intensive part is the computation of the conditional log-likelihood log p(y i,j |z i , θ ; t i,j ), which amounts to synthesize the candidate data for the current values of the fixed and random-effects (θ, z i ) and measure its discrepancy with the true observation y i,j . The synthesis of the data follows the generative model introduced in Section 2 and essentially requires the integration of ordinary differential equations.

Algorithm 1 details the procedure, where |E| denotes the dimension of the extrinsic shape shape E, and F (.) the cumulative distribution function of the standard Gaussian. The "source index" refers to the ICA components.

Calibration

Initialization procedure for model calibration

A good choice of initial parameters θ [0] and latent variables z [0] improves the convergence speed of the calibration algorithm. We propose in this section an initialization procedure that combines several elementary shape analysis tools. Given a longitudinal data set of shapes {y i,j , t i,j } i,j : Algorithm 1: Compute the complete log-likelihood. input : Longitudinal data set of shapes {y, t} = {y i,j , t i,j } i,j .

Population parameters θ = y 0 , c 0 , m 0 , a 0 , t 0 , σ τ , σ α , σ .

Individual parameters z = (z i ) i with z i = α i , τ i , s i . output: The complete log-likelihood Q = log p({y}, z, θ ; {t}). Set Q = 0.
// initialization / * compute the squared residuals 2 i,j for each visit * /

Compute the initial velocity field v 0 = Conv(c 0 , m 0 ). Compute the geodesic γ : t → Exp t 0 ,t Id (v 0 ). // see [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF] for the source index l = 1 to q Compute the l-th

column of A 0,m ⊥ 0 , projecting Col l (A 0 ) on m ⊥ 0 . Compute the initial velocity field w l = Conv c 0 , Col l (A 0,m ⊥ 0 ) . Compute the parallel transport w l : t → P w l γ (t). // see [41] end for the individual index i = 1 to n for the visit index j = 1 to n j Compute the time-warped age ψ i,j = α i • (t i,j -τ i ) + t 0 . Compute the initial velocity field v i,j = q l=1 s (l) i • w l (ψ i,j ). Compute φ i,j = Exp γ(ψ i,j ) (v i,j ) • γ(ψ i,j ). // see [22] Compute the squared residual 2 i,j = d E (y i,j , φ i,j y 0 ) 2 .
/ * add the model log-likelihood log p(yi,j |zi, θ ; ti,j

) * / Update Q ← Q -1 2 |E| • log σ 2 + 2 i,j /σ 2 . end / * add the random effects log-likelihood log p(zi|θ) * / Update Q ← Q-1 2 log σ 2 τ +(τ i -t 0 ) 2 /σ 2 τ + s i 2 2 +log σ 2 α + log(1 -F (-1/σ α )) 2 + (α i -1) 2 /σ 2 α . end / * add the Bayesian prior log-likelihood log p(θ) * / Update Q ← Q -1 2 (t 0 -t 0 ) 2 /ς 2 t + m τ (log σ 2 τ + ς 2 τ /σ 2 τ ) + m α (log σ 2 α + ς 2 α /σ 2 α ) + m (log σ 2 + ς 2 /σ 2 ) . // log p(θ1) Update Q ← Q -1 2 y 0 -y 0 2 2 /ς 2 y + c 0 -c 0 2 2 /ς 2 c + m 0 -m 0 2 2 /ς 2 m + A 0 -A 0 2 2 /ς 2 A . // log p(θ2)
1. estimate a Bayesian atlas model (see [START_REF] Gori | A Bayesian Framework for Joint Morphometry of Surface and Curve meshes in Multi-Object Complexes[END_REF]) from the baseline shapes {y i,1 } i , to get an approximate populationlevel average geometry y 0 , c 0 as well as n space-shift momenta m i mapping this geometry to the baseline observations, and an estimate of the noise level σ ; 2. for i = 1, ..., n, estimate a geodesic regression model (see [START_REF] Fishbaugh | Geodesic regression of image and shape data for improved modeling of 4D trajectories[END_REF]) from the individual time-series {y i,j } j , then parallel transport (see [START_REF] Louis | Parallel transport in shape analysis: a scalable numerical scheme[END_REF]) the computed individual initial momenta back to the mean geometry y 0 , c 0 along the corresponding space-shift m i , and finally compute the Euclidean average of those w i to get an approximate population-level mean momenta m 0 = w i i ; 3. for i = 1, ..., n, initialize the individual temporal param-

eters with τ i = t i,j j , α 2 i = w i •G c 0 •m 0 m 0 •G c 0 •m 0 if this value is
positive and α i = 1 otherwise, then compute σ τ and σ α according to equations ( 16) and ( 17) respectively; 4. solve a standard ICA problem with q components from the collection of space-shift momenta w i,m ⊥ 0 preliminarily projected on the orthogonal space to m 0 , and set A 0 as the estimated mixing matrix; 5. shoot forward the mean geometry y 0 , c 0 in the direction m 0 with length t 0t 0 where t 0 = t i,1 i and t 0 = t i,j i,j to get longitudinally centered estimates c 0 , y 0 , m 0 , and parallel-transport the q columns of A 0 along the same geodesic to obtain A 0 ; 6. personalize the model given by the initial parameters

θ [0] = (y 0 , c 0 , m 0 , A 0 , t 0 , σ τ , σ α , σ ) to obtain z [0] .

The MCMC-SAEM-GD algorithm

Calibration is a computationally-intensive task for mostly two reasons. First, the optimized variable θ is of high di-

mension |θ| = 4 + |y 0 | + d • p • (2 + q)
where d is the dimension of the ambient space, p the number of control points, q the number of sources, and |y 0 | the number of vertices necessary to describe the template mesh. Second, the optimized function requires the computation of the integral over the latent variables. The term p({y}, z, θ ; {t}) can only be evaluated for some given random-effect values z, by solving sets of ordinary differential equations (see Algorithm 1). In this paper, we propose to address this computational challenge by combining the MCMC-SAEM algorithm with gradient descent (GD). The backbone of this algorithm is the SAEM algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of the em algorithm[END_REF], which is a stochastic approximation (SA) of the classical expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]: are alternated a stochastic simulation step z [k] ∼ p(z|{y}, θ [k-1] ; {t}) of the latent variables followed by a deterministic update of the model parameters θ [k] ← θ (z [k] ).

Algorithm 2: Calibration with MCMC-SAEM-GD.

input : Dataset y. Initial parameters θ [0] and z [0] . Sequence of step-sizes (ρ

[k] ) k . Sampling variances (σ (b) ) b . output: Estimation of θ m ≈ θ [k] . Set k = 0 and S [0] 1 = S 1 (z [0] ). // initialization, eq. (14) repeat Set k ← k + 1.
/ * block Gibbs symmetric random walk sampling * / foreach random variable z (b) in (z (1) , z (2) , z

(3) ) = (α, τ, s) do Draw a candidate z (b) ∼ N (z [k-1](b) , σ 2 b ). Let z = (z [k](1) , ..., z [k](b-1) , z (b) , z [k-1](b+1) , ...). Compute the ratio ω = log p({y},z,θ [k-1] ;{t}) p({y},z [k-1] ,θ [k-1] ;{t}) . // alg. 1 Draw u according to the uniform distribution u ∼ U (0, 1). if log u < ω then z [k](b) ← z (b) else z [k](b) ← z [k-1](b) . end Adapt the proposal variances (σ (b) ) b .
// see [START_REF] Atchade | An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift[END_REF] / * analytical update rule for θ1

(classical SAEM) * / Set S [k] 1 ← S [k-1] 1 + ρ [k] • S 1 (z [k] ) -S [k-1] 1
. // eq. ( 14)

Set θ

[k] 1 ← θ 1 (S [k] ) // eqs. (15) - (18)
/ * gradient-descent-based update heuristic for θ2

* / Solve θ 2 = argmax θ 2 p({y}, z [k] , θ [k]
1 , θ 2 ; {t}) by GD. // alg. 1

Set

θ [k] 2 ← θ [k-1] 2 + ρ [k] • θ 2 -θ [k-1] 2
.

// heuristic until convergence In [START_REF] Kuhn | Coupling a stochastic approximation version of em with an mcmc procedure[END_REF], the authors introduce the MCMC-SAEM algorithm, where the simulation step is replaced by a Markov chain Monte-Carlo (MCMC) step while still preserving the theoretical convergence properties. In this paper, an analytical update rule θ cannot be found for all the parameters θ: we use a gradient descent approach to overcome this difficulty, and we name MCMC-SAEM-GD the global resulting algorithm. Algorithm 2 gives a high-level pseudo-code of the proposed procedure. The sufficient statistics write:

S t = 1 n n i=1 τ i , S α = 1 n n i=1 (α i -1) 2 , (14)
S τ = 1 n n i=1 τ 2 i , S = 1 |E| • n • n i i n i=1 ni j=1 2 i,j ,
where 2 i,j = d E {y i,j , ExpP vi γ ψ i (t i,j ) y 0 } 2 and n i i is the average number of longitudinal observations per subject. The update rules write:

t 0 = ς 2 t S t + σ τ 2 n t 0 • ς 2 t + σ τ 2 n -1 (15)
σ τ = S τ -2 t 0 S t + t 0 2 + m τ n ς 2 τ 1 2 • 1 + m τ n -1 2 (16)
σ α = S α + m α n ς 2 α 1 2 • 1 - f (-1/σ α )/σ α 1 -F (-1/σ α ) + m α n -1 2 (17)
σ = S + m |E| n n i i ς 2 1 2 • 1 + m |E| n n i i -1 2 (18)
where f (.) is the probability density function of the standard normal distribution. Both the coupled set of equations ( 15)-( 16) and the implicit equation [START_REF] Debavelaere | Clustering of longitudinal shape data sets using mixture of separate or branching trajectories[END_REF] can easily be solved by iterative update. Equation ( 18) is closed-form.

Implementation details

The sequence of ρ [k] required by Algorithm 2 is chosen to be constantly equal to 1 in a preliminary "burn-in" phase of the calibration procedure, and then decreases with the iterations with an exponential decay. The fanning numerical scheme is used to compute the parallel transport along geodesics in a scalable manner [START_REF] Louis | Parallel transport in shape analysis: a scalable numerical scheme[END_REF][START_REF] Louis | A fanning scheme for the parallel transport along geodesics on riemannian manifolds[END_REF][START_REF] Younes | Jacobi fields in groups of diffeomorphisms and applications[END_REF]. A block Metropolis-Hastingwithin-Gibbs approach is used for the MCMC sampling step, where each variable α i , τ i and s i are successively sampled. Several transition kernels can be chained in order to decrease the correlation between z [k-1] and z [k] . Proposal variances are dynamically adapted during the iterations to ensure that the acceptation rates remain close to 30% [START_REF] Atchade | An adaptive version for the metropolis adjusted langevin algorithm with a truncated drift[END_REF]. The optimization problem for the update of θ 2 is solved by steepest gradient descent. The gradients of the complete log-likelihood are obtained by autodifferentiation using the PyTorch library. The PyKeops library [START_REF] Charlier | An efficient kernel product for automatic differentiation libraries, with applications to measure transport[END_REF] implements smart autodifferentiation methods for convolution intensive computations to avoid memory overflows.

Personalization

Once the model is calibrated using a training data set, any individual data sequence {y i,j , t i,j } j can be reconstructed by fitting the model (whether it was part of the training set or not, i.e. i ≤ n or i > n respectively). This procedure, called here personalization, consists in solving the optimization problem defined by equation [START_REF] Charlier | An efficient kernel product for automatic differentiation libraries, with applications to measure transport[END_REF]. Note that all individuals can be treated independently i.e. equation ( 13) is equivalent to solving several sub-problems of the form z m i = argmax zi log p({y i,j } j , z i , θ m ; {t i,j } j ). The computed optimal latent variables z m i give in turn the spatiotemporal coordinates of the individual trajectory in the reference frame of the calibrated model θ m . We use the L-BFGS optimization method [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF], where gradients are automatically computed using the PyTorch autodifferentiation library.

Simulation

The purpose of the simulation is to take advantage of the generative nature of the model to generate an entirely synthetic data set that reproduces the characteristics of the original training data set.

Given a longitudinal data set, the calibration followed by the personalization to the training data yields a normative model of progression, a spatiotemporal coordinate system (both being encoded by the parameters θ m ) and the coordinates of each individual in this reference frame (i.e. z m = (z m i ) i ). We denote p(z m , {t}) the empirical joint distribution of those individual parameters and of the corresponding time-indices {t}. We simulate synthetic data {y s } by sampling random variables from this empirical distribution, and generate data following the generative model [START_REF] Bône | Learning low-dimensional representations of shape data sets with diffeomorphic autoencoders[END_REF].

We use statistic functions ζ (most often not sufficient, similarly to [START_REF] Marin | Approximate bayesian computational methods[END_REF]) to evaluate to which extend the simulated data resemble the training data:

ζ({y s }) ≈ ζ({y}), with {y s } iid ∼ p({y}|z s , θ m ; {t s }) (z s , {t s }) iid ∼ p(z m , {t}) . ( 19
)
For visualization purposes, we may choose to ignore the calibrated variance of noise σ m and replace it with the degenerated value σ = 0. This choice will generate smoother shapes, and we will call such simulations "without noise".

Experiments

Validation on synthetic shape data

In this section, the calibration, personalization and simulation algorithms are validated on a synthetic data set in 2D.

The simulation algorithm is first used in Section 5.1.1 to generate a synthetic shape data set from a chosen ground truth model. We use then the calibration method to infer the model parameters from the synthetic data set. The performance and the stability of the calibration algorithm is evaluated in various settings in Section 5.1.2. The calibrated model is personalized in Section 5.1.3, and the learned individual parameters are compared to the true values. Eventually, we re-simulate a synthetic data set from the calibrated model, and assess in Section 5.1.4 to which extend this new synthetic data set has similar statistics as the original data set.

Simulating synthetic data from a ground truth model

We choose values of fixed effects θ = (y 0 , c 0 , m 0 , A 0 , t 0 , σ τ , σ α , σ ), which specifies a normal distribution of shape trajectories. The chosen geometrical parameters y 0 , c 0 , m 0 , A 0 are shown in Figure 3. In addition, we choose t 0 = 0, σ τ = 2, σ α = 0.2 and σ ∈ {0.00, 0.01, 0.02, 0.03, 0.05}.

We use the generative model to simulate a total of n ∈ {50, 100, 200} individual trajectories and to sample them at several time-points {t i,j } ni j=1 . We draw the number of observations for each individual n i ≥ 2 according to a shifted Poisson distribution with parameter E(n i ) ∈ {3, 5, 7, 9}. We

-2 -1 0 1 2 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 -2 0 2 -2 -1 0 1 2 -2 0 2 -2 -1 0 1 2 -2 0 2 -2 -1 0 1 2 -2 0 2 -2 -1 0 1 2
Fig. 3: Visualization of the parameters θ 1 = (y 0 , c 0 , m 0 , A 0 ). The template shape y 0 is in solid black, the control points c 0 are the five dot points in either blue or red, the momenta m 0 is the bold blue arrow, and the four columns of A 0 are the bold red arrows. The velocity fields corresponding to the momenta m 0 or the geometrical components of A 0 are respectively represented with light blue or light red arrows. The green dots and arrows on the top figure mark the four landmark positions that will be considered for the statistic ζ, in order to validate the simulation algorithm in Section 5.1.4. finally impose that the individual time-points {t i,j } j are uniformly distributed in the observation interval [t i,1 , t i,ni ] = [t i,0 -∆t i /2, t i,0 + ∆t i /2], where both the observation time window ∆t i = t i,nit i,1 and the mid-point t i,0 are drawn according to normal distributions:

∆t i iid ∼ N E(n i )-2, σ 2 τ and t i,0 iid ∼ N (t 0 , σ 2 τ )
. Figure 4 displays some generated data in the reference case where σ = 0.02 and E(n i ) = 7.

Model calibration

The model calibration outputs are the estimated population parameters θ = (θ 1 , θ 2 ) with θ 1 = (t 0 , σ τ , σ α , σ ) and θ 2 = (y 0 , c 0 , m 0 , A 0 ). They are expected to be close to the MAP θ m defined by equation [START_REF] Chakraborty | Statistics on the space of trajectories for longitudinal data analysis[END_REF].

Computing the MAP. Because only a finite number of data points are available and a Bayesian prior p(θ) is assumed on the parameters θ, the MAP θ m does not correspond exactly to the ground truth parameters θ t . The calibrated parameters θ are expected to converge towards the corresponding MAP parameters θ m when the number of iterations goes to infinity -and this section experimentally verifies it, when the MAP parameters θ m are known to converge towards the ground truth parameters θ t when the number of observations goes to infinity. For each configuration of ground truth parameters θ t and particular random sampling z t of the generative model they define, θ m 1 can be analytically computed with equations (15-18) and θ m 2 approximated by a steepest gradient descent approach initialized with θ t and z t (see Algorithm 2). The calibration error between θ and θ m is analyzed and discussed in details in the rest of this section, when the statistical error between θ m and θ t is only computed in the reference configuration. Using the performance metrics introduced in the following paragraph, the second line of Table 1 (in italicized text) gives the corresponding quantitative normalized distances, which remain below 6% in all cases.

Normalized error metrics. The error for the scalar parameters θ 1 = (t 0 , σ τ , σ α , σ ) is measured by the absolute difference between estimated and MAP values. The error for t 0 is normalized by the characteristic population observation window, that we define as 2

• (1 + σ α ) • [E(∆t i )/2 + σ τ ].
The remaining errors for σ τ , σ α , σ are respectively normalized by the true standard deviations σ τ = 2, σ α = 0.2, and the estimated noise level by the Bayesian atlas model [START_REF] Gori | A Bayesian Framework for Joint Morphometry of Surface and Curve meshes in Multi-Object Complexes[END_REF] computed during the initialization pipeline described in Section 4.3.1. The error on the template shape y 0 is assessed as the maximum point-to-point residual distance, and normalized by the conservative value of 3 spatial units for the characteristic size of the considered shape (see Figure 3). The control points c 0 and the momenta m 0 are jointly evaluated through the 2 distance of the estimated velocity field v 0 = Conv(c 0 , m 0 ) to the MAP value, normalized by the 2 norm of this MAP velocity field. Finally, the convergence The reference configuration corresponds to a noise level σ = 0.02, an average number of visits per subject E(n i ) = 7, q = 4 allowed components of geometrical variability, and n = 100 input subjects. The second line gives the discrepancy between the ground truth (used for generating the data) and the MAP (used for evaluating the calibration performance), in the reference case. We call this discrepancy the statistical error, by opposition to the calibration error. of modulation matrix A 0 is assessed by measuring the mismatch between the sets of space-shifts that can be generated from the pair of parameters (c 0 , A 0 ): (i) the estimated modulation matrix is first re-projected on the MAP control points, (ii) the linear subspaces generated by the columns of the MAP and estimated modulation matrices are defined, (iii) the matrix representations of the projectors over those subspaces are computed, (iv) the average of the four greatest eigenvalues of their difference captures the mismatch be-tween those projectors, (v) the result is normalized by the largest eigenvalue of the MAP projector.

y 0 (%) c 0 , m 0 (%) c 0 , A 0 (%) t 0 (%) σ τ (%) σ α (%) σ (
Evaluation setups and results. In addition to the previously introduced setups, configurations with varying allowed number of sources q ∈ {2, 4, 6} are also evaluated. We call the configuration with σ = 0.02, q = 4, n = 100, and E(n i ) = 7 the reference one. Augmented with the 11 configurations differing from this reference by a single parameter, a total of 12 calibration problems are defined. Each is solved 10 times by running the stochastic MCMC-SAEM-GD algorithm.

Figure 5 plots the evolution of the error metrics across the allowed 200 iterations for the reference configuration: the black lines correspond to the 10 different runs, and in green is represented their mean and standard deviation. The algorithm is stable i.e. converges to similar results at each run: the final standard deviation of the error is smaller than 10% of the maximal (initial) error, for all parameters. The two regimes of the algorithm can be identified: the burnin phase for the first half of the iterations where the stepsizes ρ [k] remained fixed to 1, followed by the concentration phase where the step-sizes decrease geometrically. We can finally notice that σ α is estimated with more variance than other parameters, and than σ τ in particular. This suggests that adding higher-order components to the time-warp functions ψ i would come with estimation challenges. Table 1 gives the average final error metrics and the associated standard deviations. Those standard deviations remain in all but one case below 3%, underlining the stability of the estimation algorithm. In most cases, the parameters are estimated with less than 10% of error: exceptions only appear in the configurations with very low noise levels σ ≤ 0.01 or underestimated number of geometrical sources q = 2. Interestingly, a higher noise level σ does not necessarily correlate with a degraded estimation of the parameters. The presence of noise can actually help the algorithm to better explore the space of parameters: it seemed that local minima may be harder to escape for very low levels of noise. In particular, the estimation performance of the noise variance σ 2 improves when the true value increases. Table 1 studies also the impact of the length of the observation period E(n i ). Long periods generally favour more accurate estimation of the parameters v 0 = Conv(c 0 , m 0 ) which encodes the direction of the progression, and (σ τ , σ α ) which capture its dynamical variability. However, because of compensation mechanisms that may take place at the individual level between α i and τ i , it is rather the joint quantity σ τ +σ α that is clearly better estimated when E(n i ) increases than σ τ and σ α independently.

The same table compares also the estimation quality when the true number of sources is underestimated (q = 2), perfectly chosen (q = 4) or overestimated (q = 6). The reconstruction ability of our model, measured by σ , increases with q, and seems to saturate once the optimal number of sources has been reached. The large estimation error made on σ in the case q = 2 comes from the fact that data was simulated from exactly four geometrical components of comparable importance (see Figure 3), thus creating a strong reconstruction performance thresholding effect when choosing q < 4. One can expect smoother variations of the reconstructive performance on real data sets, which do not result from the exact simulation of the generative model. Parameters are in majority less well estimated in the q = 2 configuration, and at comparable distance to the MAP in the two remaining ones. Finally, Table 1 shows that the number of training subjects n has a major influence over the quality of the estimation. Almost all metrics are improved in the configuration with n = 200 subjects.

In conclusion, the proposed MCMC-SAEM-GD algorithm successfully solves our model calibration problem in varied configurations. The stochastic procedure is stable across independent repetitions. The presence of noise in the training data is well-handled, and actually seems to act as a good regularizer for the estimation procedure. An underestimated number of sources does not harm the convergence of the procedure, but mostly impairs the reconstruction ability of the learned model. This number should therefore be gradually increased to meet the reconstruction goals of the experimenter, keeping in mind that an intrinsic optimal performance will be reached when q is large enough. Finally, increasing the number of subjects or the number of visits are both beneficial for model calibration.

Personalization after calibration

Once calibrated, the longitudinal shape models are personalized to the training data. The estimated individual parameters α i , τ i , s i are compared to their true value. In order to be comparable with the true sources, the estimated sources are first brought back to the cotangent space defined by the true control points c t 0 by solving Conv(c t 0 , m t i ) = Conv(c 0 , m i ). Figure 6 plots the estimated z i against the true corresponding values. The acceleration factors are well aligned on the bisector. The onset ages and sources are also estimated with a low variance, but with a non-negligible bias.

i ) i R 2 = 0.965 -2 -1 0 1 2 true (s i ) i -3 -2 -1 0 1 2 estimated (s i ) i R 2 = 0.999 geometrical component 1 geometrical component 2 geometrical component 3 geometrical component 4
Fig. 6: Comparison of the estimated individual parameters z i = (α i , τ i , s i ) after personalization of the mean calibrated model to the simulated observations, in the reference scenario. In each scatter plot, the identity is represented by the solid black line. The R 2 value for the sources is an average over the four geometrical components.

This effect is due to the fact that the estimation of individual parameters during personalization may compensate for some error made during the estimation of population parameters during calibration. Time-shifts τ i may compensate an error on the reference time t 0 . Acceleration factors α i may compensate for an error in the norm of v 0 . Sources s i may compensate for an error in the norm of the columns of A 0 .

These effects do not question the identifiability of the model, but rather suggest that, for a finite number of observations, the likelihood may have a rather flat maximum, for which a range of parameter values may reconstruct data almost equally well. Finally, two outliers can be noticed in Figure 6 for the pace of progression α i , as well as for the onset age τ i . These outliers correspond to extremely reduced windows of observation t i,ni -t i,1 , respectively equal to 0.07 and 0.20, when the theoretical mean is equal to 5.

Table 2 summarizes the results in all configurations, giving for each of the twelve considered setups the median error and associated median absolute deviation on z i = (α i , τ i , s i ) when personalizing the average calibrated model. The median is reported instead of the mean because it is more robust to outliers. Focusing on the estimation variability, it appears that the sources s i are the best estimated parameters, followed by the pace of progression α i and the onset ages τ i . The results are given for the reference scenario plus eleven perturbed scenarii, where either the noise level σ , the average number of visits per subject E(n i ), the allowed number of geometrical components q or the number of subjects n is varied. The estimation of the pace of progression α i quickly deteriorates with increasing levels of noise σ , reaching almost 25% of the true standard deviation σ α = 0.2 in the most noisy configuration. The estimation of the onset ages τ i and sources seems more robust, with no clear tendency. The estimation of the pace α i improves when the number of visits per subject E(n i ) increases. The same trend can be noticed for the onset age τ i , although with a reduced amplitude. The sources s i remain well-estimated in all scenarii. No clear difference can be noticed between the reference q = 4 and the over-estimated number of geometrical components case q = 6, suggesting that adding components does not hamper the personalization of a calibrated model. However, underestimating this number of components with q = 2 deteriorates the estimation of the sources s i , and the dynamical parameters α i and τ i to a lesser extend. As in the previous section, we interpret this large performance drop due to the fact that data was simulated according to exactly four geometrical sources of similar magnitude (see Figure 3): in real data sets, one may expect the estimation performance to change more smoothly with q. Finally, an increased number of subjects n allows a better performance of the personalization algorithm, especially for the onset age τ i and source s i parameters.

Simulation after calibration and personalization

After calibration and personalization, the learned model and empirical distribution of the random effects can be used to simulate entirely synthetic shape trajectories. Figure 4 gives some randomly selected samples from such simulated trajectories for the reference scenario, where (see equation ( 19)):

the fixed effects θ m are averages over the 10 calibrations; the random effects z s are drawn according to independent normal distributions with mean and standard deviations equal to the values given by Table 2; the visit ages t s are drawn according to the true procedure based on the average calibrated values for t 0 , σ τ , and the empirical average n i i for E(n i ) (see Section 5.1.1).

Dynamic facial expression

Data and preprocessing

The Birmingham University 3D dynamic facial expression database [START_REF] Yin | A high-resolution 3d dynamic facial expression database[END_REF] gathers short video sequences from 101 subjects (of which 58 female, 43 male). Each subject mimics in 6 distinct sequences basic emotions which are Anger, Disgust, Fear, Happiness, Sadness and Surprise. For each of those 606 sequences we uniformly extract 8 frames spanning from the first to the 36-th one, which correspond to a subsampling of the first 1.4 seconds of each video. We do not work directly with the images, but with a set of 75 semi-automatically extracted landmarks, which come with this data set. Every set of 3D landmarks is registered to a reference one by similarity-based Procrustes alignment.

Model calibration: learned emotion models

We learn 6 distinct longitudinal atlas models: one per emotion, calibrated on the n = 101 sequences of n i = 8 frames for all subjects i. We choose q = 10 sources. Figure 8 shows in green the estimated average scenario for each emotion. Qualitatively, those average scenarii show a typical pattern of facial expression. The Disgust, Fear, Happiness and Surprise models feature large displacements in the area of the mouth in particular. The Sadness expression is more mute, with a subtle displacement of the eyebrows. The Anger model shows a combined displacement of both eyes and eyebrows.

Gender-specific emotion patterns

The estimated models are personalized to the corresponding training data sets, giving for each sequence an optimal z i = (α i , τ i , s i ). We only focus on exploiting the individual source parameters s i ∈ R q = R 10 in this section. For each model, we fit a 1D partial least square regression model for predicting the gender from a linear combination of the sources variables s i [START_REF] Abdi | Partial least square regression (pls regression)[END_REF]. We then test whether the linear combination of the sources are significantly different between men and women using a Student t-test. All p-values are smaller than 10 -5 , thus showing significant differences in the geometry of the face between genders that are independent of the pattern of expression. Figure 8 shows the typical scenario for men and women, which are built by translating the mean scenario in the direction of the average of the sources for each gender (in black). For all emotion models, male subjects tend to have wider faces than females, as it can very clearly be seen in the area of the cheeks or of the nose for the Anger and Surprise models. Table 3: Average confusion matrix across 5-fold linear discriminant classification. The sequence features consist in a 12-scalar vector that stacks the 6 pairs of dynamical parameters α i , τ i obtained by personalizing the 6 emotion models. The average accuracy is 67.08 %.

Application to classification

We propose to automatically recognize the emotion from a sequence based on the personalization of each facial expression model to the sequence. We propose here to use the dynamic variables α i , τ i for classification. More precisely, we perform a 5-fold cross-validation ensuring that each group is gender-balanced. For each split:

six longitudinal shape models are learned on the training sequences for each emotion; these models are personalized to all the 606 sequences:

for each sequence a total of 6 z i vectors are therefore estimated; for each sequence, the estimated temporal parameters α i , τ i are stacked into vectors of 6 × 2 = 12 scalars; these feature vectors are used to train and test a simple linear discriminant classifier on the corresponding train and test sequences.

Table 3 gives the confusion matrix obtained with this procedure, averaged over the 5 folds. The average classification accuracy is 67.08 %, above the chance level which amounts to 16.67 %. For comparison, [START_REF] Amor | 4-d facial expression recognition by learning geometric deformations[END_REF] reported an average accuracy of almost 100 %, [START_REF] Sun | Facial expression recognition based on 3d dynamic range model sequences[END_REF] of 90.44 %, and [START_REF] Fang | 4d facial expression recognition[END_REF] of 74.63 %.

We emphasize however that our performance is achieved:

using the default linear discriminant analysis from the sklearn library, without any hyperparameter tuning as in [START_REF] Amor | 4-d facial expression recognition by learning geometric deformations[END_REF] with random forest, in [START_REF] Sun | Facial expression recognition based on 3d dynamic range model sequences[END_REF] with hidden markov model or in [START_REF] Fang | 4d facial expression recognition[END_REF] with radial support vector machine; on all the 606 available sequences, without any manual selection of a subset of 60/101 subjects as it is done in [START_REF] Amor | 4-d facial expression recognition by learning geometric deformations[END_REF][START_REF] Sun | Facial expression recognition based on 3d dynamic range model sequences[END_REF] or of 507/606 sequences as done in [START_REF] Fang | 4d facial expression recognition[END_REF]; based only on 12 intuitive scalar features per sequence, that encode how an individual emotional pattern dynamically compares to population models of basic emotions.

From this experiment that has not been particularly tuned to achieve best classification performance, we conclude that our model captured shape characteristics that are specific to We select all the T1-weighted MRIs of subjects that were diagnosed as presenting mild cognitive impairements at some visit, and diagnosed as converted to Alzheimer's disease at some later visit. See Table 4 for summary statistics. This data set amounts to a total of 1993 visits from n = 322 subjects. Second-take "re-test" MR images are available for 1838 of those visits and will be used to estimate the noise the volumetric segmentations of the hippocampus computed with FreeSurfer are transformed into meshes using the aseg2srf script of July 20094 , the resulting meshes are decimated by a 88% factor using Paraview 5.4.1 5 [2], they are aligned using the previously-computed global affine transformation estimated with the FSL software, residual pose differences among subjects are removed by rigidly aligning the meshes from the baseline image of each subject to the corresponding hippocampus mesh in the Colin27 atlas image, this transformation with 6 degrees of freedom being computed with the GMMReg script of June 20086 [START_REF] Jian | Robust point set registration using gaussian mixture models[END_REF], the same transformation is finally used to align the meshes from the follow-up images of the same subject.

Models of atrophy of the hippocampus

We calibrate two longitudinal shape models on all the 1993 meshes of the left and right hippocampus respectively, choosing in both cases q = 8 sources. The deformation kernel width is set to σ = 10 mm. The current distance is used to compute distances between meshes without point correspondence, with a kernel width of σ E = 5 mm [START_REF] Charon | Fidelity metrics between curves and surfaces: currents, varifolds, and normal cycles[END_REF][START_REF] Vaillant | Surface matching via currents[END_REF]. Figure 9 shows the estimated average progression, which consists in an overall atrophy of both the left and right hippocampus with a specific deformation of their shape. It is worth noting that we reconstruct here the progressive atrophy of the hippocampus over more than 30 years of disease progression although patients have never been observed for more than few years. This can be achieved because the method automatically re-aligns in time the data of patients that are at different, but unknown, disease stage.

Personalization to unseen data

We assess the reconstruction performance of the calibrated models using a 5-fold cross-validation. The n = 322 subjects are split into 5 groups; 2 × 5 distinct shape models are calibrated on the training sets for the left and right hippocampus. Those models are then personalized to the unseen test subjects. To assess the goodness of fit, we measure the residual errors and compared the distribution of such errors with the noise distribution. This noise distribution is determined by measuring the distance between the two meshes extracted from the "test" and "re-test" images acquired from the same patient the same day, thus capturing all the variability due to varying image quality and its consequence in the processing. Figure 10 shows the superimposition of the distribution of the residual errors with the distribution of the differences between the meshes of the test and re-test images. The reconstruction errors are on average smaller than the intrinsic uncertainty on the data, and with a lower variance as well. The model allows therefore to reconstruct individual data at the precision of the noise. It is worth noting that this could be achieved using a reduced set of 2 × 10 scalars, which are for each hippocampus the pace of progression α i , the onset age τ i , and the eight sources s i .

Association with co-factors

We calibrate and personalize the models on whole data set, and aim to study how some genetic, biological and environmental co-factors may modulate the progression of Alzheimer's disease in patients. We therefore aim to find correlations between individual variables z i = (α i , τ i , s i ) and the following factors: gender, APOE-4 carriership, presence of amyloid plaques, education level and marital status.

To this end, the parameters α i and τ i are regressed against the five considered cofactors, and two-tailed t-tests are performed on the coefficients. A 2-block partial least square regression model [START_REF] Abdi | Partial least square regression (pls regression)[END_REF] is used to regress the eight sources s i against the five cofactors in a one-dimensional projection space. A two-tailed t-test is then performed on the weights of the multivariate regression of the linear combination of sources against the cofactors. For each case, the obtained five p-values are corrected with the Benjamini-Hochberg false discovery rate procedure [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF].

The obtained correlations for both left and right hippocampus are summarized in Table 5. The two first rows indicate that the atrophy of the hippocampus develops faster and starts earlier in female subjects. Male and female subjects present significantly different shape of their hippocampus regardless of its atrophy due to aging or disease progression. Figure 11 presents the corresponding mode of geometrical variability. Hippocampal atrophy also starts earlier in carriers of at least one 4 allele of the APOE gene, with an effect size of almost three years. The atrophy occurs at an accelerated pace in amyloid-positive subject, as well as for APOE-4 allele carriers and married subjects but only in a significant manner in the left hemisphere of the brain. Finally, the atrophy occurs earlier in married subjects, as well as in educated subjects.

The results obtained by correlating the estimated individual parameters z i with the genetic and biological factors are in line with current knowledge. The results obtained with respect to the marital status are more surprising, and should probably be taken with care as the non-married group, which represents less than 20% of the considered 299 subjects (see Table 5) is very heterogeneous. It gathers widowed, divorced, or never married subjects. Finally, we show that the atrophy starts earlier also in subjects with higher level of education. This fact is not as counter-intuitive as it appears, and actually is in line with the cognitive reserve theory [START_REF] Stern | Cognitive reserve and alzheimer disease[END_REF], which supports the idea that education can help to compensate damaged brain anatomy at the clinical level, maintaining unaltered cognitive capacities for a period of time. In other words, cognitive decline would be delayed with respect to the onset of brain atrophy in educated subjects. Since, in addition, the age at diagnosis is not correlated with the number of years of education in our dataset (r = -0.02 and p = 0.70 according to a two-tailed test based on Pearson's correlation coefficient), this explains why the subjects present an increased atrophy of their hippocampi for an increased education: they enrolled with a more advanced stage of anatomical pathology, after some years of compensation.

Simulation of hippocampus atrophy due to AD

The calibrated models and the empirical distribution of random effects z i estimated by their personalization to the training data are used to simulate synthetic progressions of the hippocampus. In order to validate such a simulation method, the simulated trajectories are sampled at several ages, and the empirical distribution of the volumes of the simulated hippocampus are compared to the distribution of the original hippocampus. The volume is commonly used as a biomarker in clinical studies, and we aim to assess if the simulated cohort could be used instead of the original one.

To do so, we simulate the same number of subjects as in the training cohort (n = 322) with the same number of timepoints and same time interval between visits. Note that we do not use the age at baseline, so that the sequence of observation time-points in the synthetic subjects may be shifted in time compared to the real ones. We simulate according to the empirical distribution of the individual parameters z i and the age at baseline. There exists indeed a correlation between the estimated time-shift τ i and the baseline age of the enrolled subjects t i,1 , as they tend to be included in the study at similar disease stage. To be more precise:

the empirical joint distribution of the time-related parameters α i and τ i augmented with the age at baseline t i,1 is computed using a kernel density estimation method; the empirical joint distribution of the time-related parameters augmented with the sources s i is captured by fitting a multivariate Gaussian distribution.

A simulated data set is then created by applying 322 times the following procedure:

draw the acceleration factor α i , the onset age τ i and the baseline age t i,1 from the corresponding kernel density; draw the sources s i from the multivariate Gaussian conditional distribution with respect to its already-drawn timerelated parameters; draw without replacement the sequence of visits of one subject i.e. the number of visits and the time intervals between them; sample the individual hippocampus trajectory defined by z i = (α i , τ i , s i ) at the baseline age t i,1 and the followup visits.

This protocol is repeated for both the left and right hippocampus, and for men and women (meaning that the estimation of the empirical distributions is done for both genders separately).

Figure 12 shows the volume distributions of the raw, reconstructed and simulated data. The cumulative distribution functions associated to the simulated and reconstructed distributions of hippocampal volumes are superimposed. This result suggests that for this volume statistic, the simulated and true data set could be used interchangeably. Raw and reconstructed distributions does not superimpose so well, because the model reconstructs smooth shapes whereas raw meshes often have small protrusion pointing outward of the surface which tend to bias volume computation (see Figure 13). This volume difference between the raw and reconstructed meshes amounts on average to 84.5 mm 3 for the left hippocampus and 67.3 mm 3 for the right hippocampus. Now validated, the simulation algorithm could be used to synthesize a data set of left and right hippocampus of any number of subjects, with any desired visit sampling. The proposed gender-wise split further allows to achieve any desired male-female balance.

Conclusion

We proposed a statistical modeling approach that represents individual data sequences as samples along continuous trajectories, these trajectories being considered as spatiotemporal perturbations of a population-average progression. The spatial warp is defined thanks to the exp-parallelization operator on manifolds. The time warps are affine time-reparameterizing functions. The spatial and temporal individual parameters position the progression of each subject in a spatiotemporal reference frame centered around the average trajectory of the population.

We proposed calibration, personalization and simulation algorithms to address different statistical questions. The calibration algorithm combines the MCMC-SAEM stochastic approach with gradient descent to estimate the underlying common process and its spatiotemporal variability from a longitudinal data set of shapes. It does not require a common time reference to be available across individual processes, which furthermore may be observed each for only short periods of time. Personalizing such calibrated models to a new individual data yields quantitative, low-dimensional and interpretable measures of how the progression of an individual deviates from a normative scenario. These parameters include an acceleration factor and a time-shift on the one hand, and geometrical sources of variability on the other hand. Such individual parameters offer relevant features for classification or correlation tasks, in a post-processing step. The generative nature of the proposed model naturally offers a simulation algorithm, which can generate entirely synthetic data sets. Such data set may be sampled at any desired temporal frequency, for any number of subjects and with a full control over the population characteristics, for instance in terms of gender balance.

We emphasize that the proposed modeling approach is able to deal with meshes without any assumption on their topology, in particular without assuming point-to-point correspondence. It may be extended easily to deal with images or other geometric primitives, provided that one can define a metric between such objects.

The three proposed algorithms were validated in varied simulated configurations, demonstrating their ability to retrieve the true parameters or reproduce the original data distribution. They were illustrated on a data set of facial expressions, showing the relevance of the learned normative scenarios and the potential of the spatiotemporal parameters for classification. We apply the method also to large medical data set of patients that develop Alzheimer's disease. The average scenarii of atrophy for the hippocampus subcortical structures are in line with current medical knowledge. Individual sequences are successfully parametrized by 10 scalar spatiotemporal coordinates in the calibrated reference frames. Correlating these coordinates with genetic, biological and environmental factors gives valuable insights into protective factors influencing age at onset or pace of progression. We also evidence typical shape differences across sub-groups, which are independent of the shape changes due to ageing or disease progression.

The calibration algorithm is computationally intensive: estimating a model of hippocampus progression took around a day. Our code is already parallel, combines CPU and GPU together, and offers a fine-grained initialization pipeline. Fur-ther pure optimization of our code (among which multi-GPU support, fast Fourier transforms for convolutions) is planned, as well as evaluating the performance of variational methods for calibration -which are not trivial to implement in a longitudinal context without a fixed number of observations per individual.

As for any modeling approach, our model relies on some assumptions. For instance, subjects are considered to follow trajectories that are parallel to the population average. This hypothesis may be alleviated by introducing drift parameters to model a progressive deviation from the average scenario. Such a development would add to the complexity of the model, which may require to have even more data to be calibrated. Further extensions would consider also to estimate not only one representative trajectory at the population level but several of them, for instance by estimating a mixture model along the lines of [START_REF] Debavelaere | Clustering of longitudinal shape data sets using mixture of separate or branching trajectories[END_REF]. Nonetheless, it is worth noting that in its current form the model is able to reconstruct data at the precision of the noise.

The model also builds on the LDDMM framework for modeling shape variability. This framework relies also on some assumptions on the geometry of the shape space. Future work will consider to learn such geometry from the data instead of relying on prior assumptions, along the lines of [START_REF] Bône | Learning low-dimensional representations of shape data sets with diffeomorphic autoencoders[END_REF] for instance. Learning other parameters such as the number of sources, using automatic model selection methods for instance, would also add to the usability of the method. where (.) is the transposition operator. This inner product defines in turn a distance metric on currents:

d E (C, C ) = C, C E + C , C E -2 • C, C E .

A.2 Practical discrete case

In practice, y is described by a finite set of T triangles in R 3 of centers c 1 , ..., c T and corresponding surface normal vectors n 1 , ..., n T . We further assume that g E is a Gaussian kernel of radius σ E . The current transform equation then writes:

C(y)(x) = T k=1 exp -x -c k 2 2 σ 2 E • n k
for any x ∈ R 3 . Similarly, the inner product formula becomes:

C(y), C(y ) E = T k=1 T l=1 exp -c l -c k 2 2 σ 2 E • n k n l
which fully specifies the distance metric d E that can be implemented in practice to measure the discrepancy between any pair of currents.

Fig. 1 :

 1 Fig.1:[Left]. Shape spatiotemporal reference frame y 0 , c 0 , m 0 , t 0 with respect to which a shape y admits coordinates t ∈ R, v ∈ v ⊥ 0 ⊂ T γ(t0) D c0 . Three spaces are involved: the manifold of control points R p×d (top), the manifold of diffeomorphisms D c0 (middle), and the shape submanifold S y0,c0 of the extrinsic shape space E (bottom). The momenta m 0 , m and the velocity fields v 0 , v are in one-to-one correspondence. The velocity field v, also called space-shift, is parallel-transported along the geodesic γ by the operator t → P v γ (t). Figure2illustrates the effect of parallel transport on D c0 . [Right]. Illustrations of the manifolds abstractly depicted on the left side of the figure. Each row displays two example elements of the corresponding geodesic (solid black lines on the left panel). The two columns correspond respectively to the times t 0 and t.

Fig. 4 :

 4 Fig. 4: Illustration of the evaluation procedure for model calibration, and subsequent personalization or simulation from the learned model. The ground truth population geodesic is plotted in black on the central line. From this model are simulated n = 100 individual spatiotemporal trajectories: three randomly-picked samples are plotted in black on the top lines. The population geodesic of the calibrated model is plotted in green on the central line, superimposed with the ground truth geodesic. This calibrated model can then be personalized to the training observations as plotted in red, or leveraged to simulate new spatiotemporal trajectories that resemble the original data set as plotted in blue.

Fig. 5 :

 5 Fig.5: Evolution of the error metrics across the 200 allowed iterations of the MCMC-SAEM algorithm, for the reference configuration: noise standard deviation σ = 0.02, q = 4 estimated components of geometrical variability, learning on a data set composed of n = 100 with on average E(n i ) = 7 visits per subject, spanning 5 time units. The 10 solid black curves correspond to 10 independent runs of the same -stochastic -MCMC-SAEM algorithm; the bold green curve is their average and the light green region indicates the associated standard deviation. The algorithm consistently converges towards similar parameters at each run, and those estimated parameters are satisfyingly close to the MAP estimate.

Fig. 7 :

 7 Fig.7: Distribution of the position of landmarks of interest in the raw (i.e. original), reconstructed (by personalization of the calibrated model) and simulated data sets, for the reference configuration. Those landmarks of interest are indicated by green dots on Figure3. The simulated distributions are similar to the corresponding raw ones, suggesting that the spatiotemporal variability of the original data set has been successfully captured.

Figure 7

 7 Figure 7 compares the distribution of vertical or horizontal positions of the tips of the original (see Section 5.1.1), reconstructed (see Section 5.1.3) and simulated observations. Those landmarks of interest are indicated by green dots and arrows on Figure 3, and form the statistic ζ introduced in equation (19). A total of 1,000 subjects are simulated, when only 100 were available for model calibration. The three distributions largely overlap, indicating that the learned distribution of shape trajectories reproduces the true distribution.

Fig. 8 :

 8 Fig.8: Learned emotion spatiotemporal models. The population geodesic is plotted in green, and the shifted progressions along the gender mode of geometrical variability are plotted in black.

  (a) Left hippocampus mean progression. (b) Right hippocampus mean progression.

Fig. 9 :

 9 Fig. 9: Typical model of hippocampus atrophy from MCI to Alzheimer's disease stage. Physiological ages (from left to right, in years): 58.6, 63.0, 67.4, 71.8, 76.2, 80.7, 85.1, 89.5, 93.9.

  in the data. All those 1993 + 1838 = 3831 images are pre-processed exactly in the same manner, starting with the longitudinal pipeline of FreeSurfer 1 (version 5.3.0)[START_REF] Fischl | Measuring the thickness of the human cerebral cortex from magnetic resonance images[END_REF][START_REF] Fischl | Whole brain segmentation: automated labeling of neuroanatomical structures in the human brain[END_REF]. The skull-stripped brains are then aligned with an affine 12degrees-of-freedom transformation onto the Colin27 average brain 2 with FSL 5.0 3[START_REF] Woolrich | Bayesian analysis of neuroimaging data in fsl[END_REF]. Meshes of the left and right hippocampus are obtained from the original images as follow:

  Left hippocampus. The mean error is 68.5 ± 15.9 mm 2 for the shape model, and 83.2 ± 36.0 mm 2 for the re-test measurement.

  Right hippocampus. The mean error is 69.8 ± 15.0 mm 2 for the shape model, and 85.2 ± 40.1 mm 2 for the re-test measurement.

Fig. 10 :

 10 Fig.10: Comparison of the generalization error to unseen data of the learned shape models and the intrinsic measurement error. The discrepancies between meshes are computed with the current metric with σ E = 5 mm, without assuming any point-to-point correspondence.

Fig. 11 :

 11 Fig. 11: Superposition of the male-like (in blue) and the female-like (in pink) hippocampus geometries, in two standard views. The letters L, R, A, P respectively indicate the left, right, anterior and posterior directions.

  Left hippocampus: the mean volume is 2958 ± 779 mm 3 for raw data, 2863 ± 693 mm 3 for reconstructed data, and 2865 ± 746 mm 3 for simulated data.

  Right hippocampus: the mean volume is 3081 ± 862 mm 3 for raw data, 3014 ± 754 mm 3 for reconstructed data, and 3063 ± 763 mm 3 for simulated data.

Fig. 12 :

 12 Fig. 12: Distribution of the left and right hippocampal volume in the raw, reconstructed and simulated data set. The simulated volume distribution is very close to the volume distribution of the reconstructed data set. The remaining bias between those two distributions and the one corresponding to the raw data comes from the smoothing behavior of the current noise model, leveraged to deal with noisy meshes without point correspondence. See Figure 13.

Fig. 13 :

 13 Fig. 13: Several views of a single example of the reconstruction of a right hippocampus structure by the longitudinal shape model. The reconstruction is the smooth white structure, and the raw data point is plotted in red.

  Institute for Research and Education, and the study is coordinated by the Alzheimer's Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.where dσ(x) denotes an infinitesimal surface element of y. The inner product of E on currents therefore writes:C(y), C(y ) E = y y n (x ) •g E (x, x )•n(x)•σ(x)•dσ (x )

Table 1 :

 1 Final average normalized performance metrics and associated standard deviations, obtained after 10 independent runs of the MCMC-SAEM algorithm in varied configurations. The algorithm is run for 200 iterations in all configurations.

	%)

Table 2 :

 2 Median of the residual errors and associated median absolute deviation times 1.4826 for the estimated individual parameters, expressed in percentage of the corresponding ground-truth standard deviations σ α = 0.2, σ τ = 2 and σ s = 1.

		∆α i (%)	∆τ i (%)	∆s i (%)
	reference	4.4 ± 11.9	44.2 ± 15.3	-11.7 ± 3.3
	σ = 0.00	-3.8 ± 4.3	42.4 ± 24.8	-5.4 ± 6.8
	σ = 0.01	-2.1 ± 6.4	11.8 ± 9.4	-11.1 ± 12.0
	σ = 0.02	4.4 ± 11.9	44.2 ± 15.3	-11.7 ± 3.3
	σ = 0.03	-7.6 ± 15.7	25.5 ± 15.0	-6.7 ± 2.5
	σ = 0.05	-21.8 ± 23.4	4.4 ± 14.7	-9.1 ± 3.0
	E(n i ) = 3	9.3 ± 32.7	-12.1 ± 19.3	-8.8 ± 3.2
	E(n i ) = 5	2.9 ± 13.0	3.2 ± 20.8	-3.6 ± 7.5
	E(n i ) = 7	4.4 ± 11.9	44.2 ± 15.3	-11.7 ± 3.3
	E(n i ) = 9	-8.2 ± 7.3	46.0 ± 16.1	-7.6 ± 2.3
	q = 2	-1.6 ± 18.0	48.3 ± 50.7	-2.4 ± 69.4
	q = 4	4.4 ± 11.9	44.2 ± 15.3	-11.7 ± 3.3
	q = 6	9.8 ± 12.3	45.3 ± 15.8	-11.7 ± 3.3
	n = 50	-12.9 ± 8.1	24.8 ± 18.3	-6.9 ± 6.1
	n = 100	4.4 ± 11.9	44.2 ± 15.3	-11.7 ± 3.3
	n = 200	9.8 ± 10.5	13.2 ± 10.9	-13.1 ± 2.4

Table 4 :

 4 Summary statistics of the medical data set of Alzheimer's disease patients.

	Number of subjects	322
	Number of visits	1993
	Average number of visits per subject (± std)	5.8 (± 2.4)
	Average age (± std)	74.0 (± 6.7)
	Sex ratio (F/M in %)	41.2 / 58.8
	Amyloid status (+/-/unknown in %)	73.2 / 7.1 / 19.7
	APOE carriership (%)	65.2
	Education (mean ± std, in years)	15.9 (± 2.8)
	Marital status (married/not married in %)	80.9 / 19.1

Table 5 :

 5 Significant associations of individual parameters with genetic, biological and environmental factors: effect sizes and significance levels of the adjusted p-values (thresholds 5%, 1%, 0.1%, 0.01%). Time-shifts τ i are in months, others have no units. Directions of space-shifts are not signed. The 23 subjects (out of n = 322) without amyloid information have been discarded.

					left hippocampus	right hippocampus
	genetic	gender APOE-4	female vs. carrier vs.	male non-carrier	α i ×1.23 τ i -12.4 months [**] -8.7 months [**] ×1.21 s i ± 0.54 [***] ± 0.57 α i ×1.22 [*] τ i -35.8 months [***] -32.5 months [**] [**] [*] [****] s i
	biological	amyloid	positive vs.	negative	α i ×1.52 τ i s i	[**] ×1.67	[*]
	environmental	marital education	married vs. nb. of years	non-married of education	α i ×1.14 τ i -42.5 months [***] -36.3 months [**] [*] s i α i τ i -3.7 months/y [**] -5.1 months/y [***] s i

available at: https://surfer.nmr.mgh.harvard.edu

available at: http://www.bic.mni.mcgill.ca/ServicesAtlases/Colin27

available at: https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/

available at: https://brainder.org

available at: www.paraview.org

available at: https://github.com/bing-jian/gmmreg

Acknowledgements The research leading to this publication has been funded in part by the European Research Council (ERC) under grant agreement No 678304 (LEASP), European Union's Horizon 2020 research and innovation programme under grant agreement No 666992 (EuroPOND) and No 826421 (TVB-Cloud), and the program "Investissements d'avenir" ANR-10-IAIHU-06 (IHU ICM) and ANR-19-P3IA-0001 (PRAIRIE 3IA Institute) The facial expression data set at the basis of Section 5.2 was built and shared by the Binghamton University. The authors warmly thank Pr. Lijun Yin for granting data access, and Peng Liu for his help in downloading the data set. Regarding Section 5.3, data collection and sharing was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). 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.; Cogstate; 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 & Development, LLC.; Johnson & Johnson Pharmaceutical Research & Development LLC.; Lumosity; Lundbeck; Merck & 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 funds 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

The theory of currents has been introduced in [START_REF] Vaillant | Surface matching via currents[END_REF], and is used in this paper to define a distance metric between pairs of meshes without any assumption on their topology, and in particular without assuming point-to-point correspondence. See also [START_REF] Charon | Fidelity metrics between curves and surfaces: currents, varifolds, and normal cycles[END_REF] for more details.

A.1 Continuous theory

Let y be a surface mesh, that we represent as an infinite set of tuples (x, n(x)) where x is a point of R 3 , and n(x) the normal vector of y at this point. Let g E : R 3 × R 3 → R be a positive-definite kernel operator, and E the associated reproducing kernel Hilbert space. We define the current transform C(y) : R 3 → R 3 ∈ E of y as: