Andrea Polzien
email: andrea.polzien@upb.de

Iris Güldenpenning
email: iris.gueldenpenning@upb.de

Matthias Weigelt
email: matthias.weigelt@upb.de

Perception during Action Observation in Combat Sports

Keywords: Action anticipation, visual perception, athletic expertise, gaze behavior, eye tracking ACM Classification, General and reference~Empirical studies, Applied computing~Psychology

Visual perception is of particular relevance for action anticipation within human social interactions. Competitive sports offer a great possibility to study visual perception, as it is necessary to anticipate action intentions of team members and/or opponents within a highly restricted timeframe. We conducted an eye tracking experiment and examined differences between athletes and novices during action observation in combat sports. Participants were instructed to respond to images of punches and kicks by pressing a left or right button. The images depicted different points-in-time during movement execution. Reaction times, error rates, fixation durations, and fixation rates were analyzed. Results showed that participants reacted faster, produced less errors, and used fewer fixations of longer duration at later points-in-time compared to earlier ones. A group effect regarding reaction times indicated expertise differences between combat-sports athletes and novices. However, no expertise-dependent differences in fixation duration and fixation rates were found.

INTRODUCTION

Interactive sports place high demands on visual perception, since visual stimuli are not only used for the control of own actions, but also for the anticipation of the actions of team members and opponents. Complex situations have to be perceived and actions have to be initiated within limited timeframes, which often last only milliseconds. Therefore, it is not surprising that there has been an increased interest in studying visual perception within various kinds of sports during the last years [START_REF] Mann | Perceptual-cognitive expertise in sport. A meta-analysis[END_REF][START_REF] Mcguckian | A systematic review of the technology-based assessment of visual perception and exploration behaviour in association football[END_REF]. The methodological approach most often taken is the expert-novice paradigm, in which the performance of expert athletes is compared to novices. In general, this kind of approach aims to understand what differentiates expert athletes from non-experts, in order to derive recommendations for training and competition, which are assumed to be helpful to develop expertise in a specific field of sports [START_REF] Müller | Expert Anticipatory Skill in Striking Sports. A Review and a Model[END_REF]. With regard to perception and anticipation, the following picture emerges from studies using the expert-novice paradigm: [START_REF] Abernethy | Expertise, visual search, and information pick-up in squash[END_REF] Experienced athletes are able to precisely meet the spatio-temporal conditions of their sports. [START_REF] Abernethy | The relationship between expertise and visual search strategy in a racquet sport[END_REF] In comparison to less experienced athletes, experts are better capable of encoding and retrieving sportspecific information, so that they can recognize action patterns more efficiently. (3) Experienced athletes compared to novices are able to detect relevant objects within their visual field faster and make better use of action-relevant (movement) cues prior to an event [START_REF] Williams | Visual perception and action in sport[END_REF].

For studying visual perception, the eye-tracking paradigm is often used [START_REF] Williams | Perceptual-cognitive expertise in sport: some considerations when applying the expert performance approach[END_REF]. Different quantitative measures, like the number of fixation and fixation duration serve as an indicator for the amount of information an observer acquires. Furthermore, qualitative parameters, like the scan path and fixation location, can be used to analyze areas of interest.

Several eye-tracking studies have examined visual perception in the context of sports expertise. However, there are divergent results regarding differences between experts and novices. Some studies demonstrated that experts use a greater number of fixations of shorter duration, as compared to novices [START_REF] Vaeyens | Mechanisms underpinning successful decision making in skilled youth soccer players: an analysis of visual search behaviors[END_REF][START_REF] Williams | Visual Search Strategy, Selective Attention, and Expertise in Soccer[END_REF], while others found the direct opposite, i.e. fewer fixations of longer duration for experts, as compared to novices [START_REF] Piras | Visual Search Strategy in Judo Fighters During the Execution of the First Grip[END_REF][START_REF] Ripoll | Analysis of information processing, decision making, and visual strategies in complex problem solving sport situations[END_REF][START_REF] Savelsbergh | Visual search, anticipation and expertise in soccer goalkeepers[END_REF], whereas other studies did not find any differences between both groups [START_REF] Abernethy | Expertise, visual search, and information pick-up in squash[END_REF][START_REF] Abernethy | The relationship between expertise and visual search strategy in a racquet sport[END_REF][START_REF] Williams | Anxiety, Expertise, and Visual Search Strategy in Karate[END_REF]. The reasons for these discrepancies seem to be diverse. They reach from typically small sample sizes in eye tracking studies [START_REF] Gegenfurtner | Expertise Differences in the Comprehension of Visualizations. A Meta-Analysis of Eye-Tracking Research in Professional Domains[END_REF], over the particular sport examined, to the complexity of the task and stimulus design [START_REF] Mann | Perceptual-cognitive expertise in sport. A meta-analysis[END_REF][START_REF] Ripoll | Analysis of information processing, decision making, and visual strategies in complex problem solving sport situations[END_REF].

In many sports, experts' use of fewer fixations of longer duration seems to be plausible from a theoretical perspective. The psychological Information-Reduction-Theory (IRT) by Haider and Frensch [START_REF] Haider | The Role of Information Reduction in Skill Acquisition[END_REF][START_REF] Haider | Eye Movement During Skill Acquisition: More Evidence for theInformation-Reduction Hypothesis[END_REF] is based on the assumption that not only information processing, but also information perception changes during skill acquisition. According to IRT, the distinction between relevant and irrelevant information is a crucial factor of learning and, as a consequence, expertise. This allows for a limitation of information processing to task-relevant information, so that the total amount of to-be-processed information is reduced. The restriction to task-relevant information is important, because only a small amount of available information on the retina can be processed at the time [START_REF] Desimone | Neural mechanisms of selective visual attention[END_REF]. For skill acquisition in sports, it is assumed that learning leads to an improving or reorganization of mental action representations [START_REF] Frank | The Representation of Motor (Inter)action, States of Action, and Learning. Three Perspectives on Motor Learning by Way of Imagery and Execution[END_REF][START_REF] Weigelt | The cognitive representation of a throwing technique in judo experts -Technological ways for individual skill diagnostics in high-performance sports. Evidence from eye movements[END_REF]. This knowledge about movements should allow experts to better distinguish between relevant and irrelevant stimulus features and should therefore lead to fewer fixations of longer duration. This strategy would be accompanied by fewer saccades, i.e. eye movements between two fixations. Since it is assumed that information processing is suppressed during saccades [START_REF] Mann | Perceptual-cognitive expertise in sport. A meta-analysis[END_REF], a gaze strategy involving fewer saccades seems to be reasonable.

In the present study, we examined gaze behavior in combat sports. A characteristic of all combat sports is the typical oneto-one situation, which requires a focused attention to a single opponent. In this way, combat sports are different from team and racket sports, which include perception of a larger sports field. To this day, however, only a few eyetracking studies have investigated visual perception in combat sports [START_REF] Milazzo | Effect of Implicit Perceptual-Motor Training on Decision-Making Skills and Underpinning Gaze Behavior in Combat Athletes[END_REF][START_REF] Milazzo | Do karate fighters use situational probability information to improve decision-making performance during on-mat tasks[END_REF][START_REF] Piras | Visual Search Strategy in Judo Fighters During the Execution of the First Grip[END_REF][START_REF] Ripoll | Analysis of information processing, decision making, and visual strategies in complex problem solving sport situations[END_REF][START_REF] Williams | Anxiety, Expertise, and Visual Search Strategy in Karate[END_REF]. These studies came to different results regarding expert-novice differences. For example, Ripoll and colleagues [START_REF] Ripoll | Analysis of information processing, decision making, and visual strategies in complex problem solving sport situations[END_REF] asked participants of different level of expertise to react to filmed sequences of a French Boxer, who was performing different boxing techniques, by manipulating a joystick. Results showed expertise-differences in gaze behavior. Experts used fewer fixations compared to novices and intermediates. However, Williams and Eliott [START_REF] Williams | Anxiety, Expertise, and Visual Search Strategy in Karate[END_REF] conducted a study with expert karate athletes and novices. Participants were instructed to react to attacks by a karate athlete presented on a screen. In this study, no expertise differences regarding number of fixations and fixation durations were found. We therefore asked combat sports athletes and novices to respond to photo-realistic stimuli, depicting four different attacking techniques (composed of straight kicks and punches with either the left or the right leg or arm), by pressing one out of four different buttons. Every attacking technique was presented at one out of four different points-in-time during movement execution. That is, for the initialization of the movement, for two significant changes of the body posture during movement execution, and for the end point of the kick (leg completely stretched out) or punch (arm completely stretched out).

The assumptions made regarding the point-in-time of movement execution were based on earlier studies by Güldenpenning and colleagues, who conducted a subliminal priming experiment with combat sports athletes and novices [START_REF] Güldenpenning | Masked priming of complex movements. Perceptual and motor processes in unconscious action perception[END_REF]. The prime (picture) depicted a combat-sport athlete performing one of two different kicks (front kick vs. roundhouse kick) at one of six different points-in-time during movement execution. Participants were asked to react to a target (picture), which depicted the end position of the front kick or the roundhouse kick. Results showed an influence of the point-in-time of movement execution presented in the prime on reaction times of both, athletes and novices, even though the prime was not consciously perceivable.

Furthermore, another study by Güldenpenning and colleagues revealed expertise-differences between athletes and novices in volleyball [START_REF] Güldenpenning | Athletes and novices are differently capable to recognize feint and non-feint actions[END_REF]. They conducted a priming experiment and asked their participants to respond to two beach-volleyball attack-hits. Results showed that beachvolleyball athletes were able to identify the technique at an earlier point-in-time during movement execution, as compared to novices. However, novices also showed a congruency effect for primes depicting later movement phases. That is, for pictures taken late during movement execution, novices also reacted faster for congruent prime and target pictures. Accordingly, we expected participants to be more familiar with later points-in-time during movement execution, as compared to earlier ones. Therefore, we expected them to react faster, to make fewer mistakes and to use fewer fixations of longer duration during observation of later compared to earlier points-in-time during movement execution. Furthermore, athletes were expected to react faster and to use fewer fixations of longer duration at every point-in-time of movement execution (i.e. general expertise effect).

METHOD Participants

Thirty participants took voluntarily part in the experiment. They received neither course credit nor monetary payment. All participants were naïve to the purpose of the study, had normal or corrected-to-normal vision, and gave written consent in accordance with the 1964 declaration of Helsinki.

Eighteen participants were assigned to the group of novices (6 females, 2 left-handed, mean age 22.6, SD = 1.9) with no experience in combat sports. Another twelve participants (1 female, mean age 27.0, SD = 6.1) were athletes of different kinds of combat sports with an average of 12.7 years (SD = 6.3) of training experience. Importantly, only athletes, who were familiar with the presented attacking techniques, were assigned to the group of combat sports athletes.

Stimuli and Set-Up

The stimulus set consisted of sixteen static images (kicks and punches with the left and right leg/arm at four different points-in-time during the movement execution). The images were taken from videos of a combat sports athlete performing the two attacking techniques. These particular techniques were chosen on the basis of their popularity in different kinds of combat sports. Videos were filmed from a front perspective, so that the perspective matched the common perspective of an athlete viewing the opponent during sparring or competition. Besides the end point of movement execution, three different points-in-time during movement execution for each technique were used as targets, with t1 being the earliest point and t4 being the end point of the movement (see Figure 1). For the purpose of inducing an apparent motion [START_REF] Alaboud | Täuschungshandlungen im Sport[END_REF], as well as directing participants' gaze to the middle of the screen, another image of the combat athlete standing in an initial position was presented prior to the target stimulus (see Figure 2). All stimuli were presented on a 24-inch monitor.

Eye-tracking data were recorded with a remote eye tracker (EyeFollower, Interactive Minds), which was located below the monitor. The EyeFollower possesses a sample rate up to 120 Hz and a gaze position accuracy smaller than 0.45°. Four different buttons on a keyboard ("<" = punch right arm/leg; "y" = punch left arm/leg; "n" = kick right leg/arm; "m" = punch left arm/leg) were used as response buttons. Button assignments to punches and kicks, respectively, were counterbalanced across participants.

Design and Procedure

A 4 x 2 mixed factorial design with the within-subjects factor "point-in-time" (t1, t2, t3, t4) and the between-subjects factor "expertise" (combat sport athletes vs. novices) was used in the present study.

Participants sat approximately 60 cm in front of a computer screen and were instructed to respond as soon and as accurately as possible to the presented kicks and punches by pressing one of four buttons. Each button was assigned to one of the four attacking techniques. A video of each technique was presented once to the participants prior to the experiment, so that participants were able to get familiar with these techniques. Afterwards, each participant completed a block of 32 randomized practice trials, before the 9-point calibration procedure for the eye tracking was conducted and the experiment started. The experiment was separated into two blocks with 120 trials each. Participants could take a rest between the blocks, if needed. Each trial consisted of a blank screen (1000 ms), a picture of the initial position (1000 ms), and the target picture, which remained on the screen until a response was given (see Figure 2). During the practice block, the subsequent trial only started if the participant had pressed the correct response button. During the experimental blocks, the subsequent trial also started after a wrong response and no error feedback was given.

RESULTS

Reaction times, fixation durations, and fixation rates were examined for correct responses, as well as error rates for incorrect responses. Each of the dependent variables was analyzed by mixed analyses of variance for the withinsubjects factor "point-in-time" and the between-subjects factor "expertise". A violation of the sphericity-assumption resulted in a correction of the p-values and degrees of freedom according to Greenhouse-Geisser. All post-hoc ttests were Bonferroni-corrected.

Reaction Times

The reaction times are displayed in Figure 3. A mixed ANOVA for the dependent variable reaction time showed a main effect for the factor point-in-time [F(1.43, 40.02) = 38.75, p < .001, ε = .48, η p 2 = .58]. Reaction times decreased with shorter temporal distances to the end point: They were highest with 1365 ms at t1, were 1121 ms at t2, were 949 ms at t3, and were 757 ms at t4. Post-hoc t-tests showed significant differences between every point in time (all t's > 3.00, all p's ≤ .001). The main effect of expertise was also significant [F(1,28) = 4.42, p = .045, η p 2 = .14]. There was a difference of 218 ms between athletes (917 ms) and novices (1135 ms). The interaction between point-in-time and expertise was not significant [(F(1.43, 28) = .91, p = .38, ε = .48, η p 2 = .03)].

Error Rates

Error rates are presented in Figure 4. A mixed ANOVA for the dependent variable error rate showed a main effect for the factor point-in-time [F(3, 84) = 120.73, p < .001, η p 2 = .81]. Error rates decreased with shorter temporal distances to the end point: They were highest with 57% at t1, were 32% at t2, were 20% at t3, and were 12% at t4. Post-hoc t-tests showed significant differences between every point-in-time (all t's > 3.00, all p's ≤ .001). The main effect of expertise [(F(1, 28) = 2.26, p = .14, η p 2 = .08)] and the interaction between point-in-time and expertise [(F(3, 84) = 1.30, p = .28, η p 2 = .04)] was not significant.

Fixation Durations

The fixation durations can be seen in Figure 5 A mixed ANOVA for the dependent variable fixation duration showed a main effect for the factor point-in-time [F(1.87, 52.37) = 22.02, p < .001, ε = .62, η p 2 = .44], as the fixation duration especially increased for t4. The mean fixation duration was longer at t4 (529 ms) than at t3 (381 ms), t2 (367 ms), and t1 (354 ms). Post-hoc t-tests confirmed these differences between the three earlier points-in-time during movement execution and t4 to be significant (all t's < 5.00, all p's < .001). There was neither a significant main effect for expertise [(F(1, 28) = .36, p = .55, η p 2 = .01)], nor a significant interaction between point-in-time and expertise [(F(1.87, 28) = .39, p = .67, ε = .62, η p 2 = .01)].

Fixation Rates

The fixation rates are presented in Figure 6. A mixed ANOVA for the dependent variable fixation rate showed a main effect for the factor point-in-time [(F(3, 84) = 34.45, p < .001, η p 2 = .55)], as fixation rates decreased from earlier to later points-in-time. The mean fixation rate was smaller at t4 (1.98 fixations/second) than at t3 (2.23 fixations/second), t2 (2.52 fixations/second), and t1 (2.65 fixations/second). Posthoc t-tests confirmed significant differences between every point-in-time (all t's > 2.00, all p's < .001), except for the difference between t1 and t2 (p = .05). There was neither a significant main effect for expertise [(F(1, 28) = .85, p = .37, η p 2 = .03)], nor a significant interaction effect between point-in-time and expertise [(F(3, 84) = 2.27, p = .08, η p 2 = .08)].

DISCUSSION

In this study, we investigated differences in reaction times and gaze behavior between combat-sports athletes and novices, who responded to photo-realistic images of two attacking techniques commonly used in combat sports (straight kicks and punches). Importantly, the images displayed different points-in-time during movement execution. We examined if combat-sports athletes are better able to identify these techniques earlier during movement execution and if they use different perceptual processes (as revealed from eye-tracking measures) than novices.

In accordance with our hypothesis, effects of point-in-time were observed for reaction time, error rates, fixation durations, and fixation rates. Participants reacted slower, produced more errors, and used more fixations of shorter duration at early points-in-time compared to later ones. Higher error rates for early points-in-time of movement execution indicate participants' difficulties to identify the techniques during the early phase of movement execution. Furthermore, the pattern of response errors reassures us that the reaction time results are not due to a speed-accuracy tradeoff. The pattern for reaction times, as well as for gaze behavior, reflects higher effort of visual search for relevant cues during early phases of movement execution compared to later ones, when the action is still unfolding. Additionally, regarding the factor points-in-time during movement execution, the reaction time results are in line with findings from priming and temporal-occlusion studies. Priming studies demonstrated stronger congruency effects with primes from later compared to earlier movement phases [START_REF] Güldenpenning | Masked priming of complex movements. Perceptual and motor processes in unconscious action perception[END_REF][START_REF] Güldenpenning | Athletes and novices are differently capable to recognize feint and non-feint actions[END_REF] and temporal-occlusion studies revealed an improvement of performance in later occlusion conditions compared to earlier ones [START_REF] Loffing | Skill differences in visual anticipation of type of throw in teamhandball penalties[END_REF][START_REF] Rosalie | Timing of in situ visual information pick-up that differentiates expert and near-expert anticipation in a complex motor skill[END_REF][START_REF] Schorer | Foveal and peripheral fields of vision influences perceptual skill in anticipating opponents' attacking position in volleyball[END_REF].

The differences between combat-sports athletes and novices regarding reaction times were also significant, and numerically large, with the difference being 218 ms in favor of the combat sports athletes. In contrast to our expectations, there were no expertise differences regarding fixation durations and fixation rates. This is in contrast to some previous studies reporting differences in gaze behavior between expert athletes and novices [START_REF] Milazzo | Do karate fighters use situational probability information to improve decision-making performance during on-mat tasks[END_REF][START_REF] Ripoll | Analysis of information processing, decision making, and visual strategies in complex problem solving sport situations[END_REF]. The absence of such differences in the present study might have several reasons: First of all, we only used an apparent motion paradigm with two static images, instead of using dynamic video stimuli. Differences between experts and novices have been reported for static images before, but fixation durations seem to be numerically smaller with static images compared to video presentations [START_REF] Mann | Perceptual-cognitive expertise in sport. A meta-analysis[END_REF]. Considering the numerical results of the data for gaze behavior in our study, it becomes clear that there is not so much room for expertise differences. This is especially true for the fixation rates, which were about 2.3 fixations per second, with an average target duration of approximately one second. Second, the absence of expertise differences might relate to the size of the stimuli. Even though previous studies found expertise differences in gaze behavior with smaller stimulus sizes [START_REF] Charness | The perceptual aspect of skilled performance in chess. Evidence from eye movements[END_REF][START_REF] Singer | Visual search, anticipation, and reactive comparisons between highly-skilled and beginning tennis players[END_REF], the majority of studies used large screens with life size images for the stimulus presentation [START_REF] North | Perceiving patterns in dynamic action sequences. Investigating the processes underpinning stimulus recognition and anticipation skill[END_REF][START_REF] Vaeyens | Mechanisms underpinning successful decision making in skilled youth soccer players: an analysis of visual search behaviors[END_REF][START_REF] Williams | Anxiety, Expertise, and Visual Search Strategy in Karate[END_REF]. In this study, stimuli were presented on a 24-inch screen with a visual angle of about 20°. The central field is regarded to have a size of about 10° [START_REF] Wandell | Foundations of vision[END_REF]. This enabled our participants to pick-up a relatively large part of the picture with only one fixation. Furthermore, it is possible that participants used peripheral vision to gather additional information.

Thus, expertise differences regarding gaze behavior should be investigated in future studies with dynamic stimuli on a large screen. Nevertheless, this study provides further insight into the cognitive processing of attacking techniques in combat sports with regard to the information pick-up.

Figure 1 .Figure 2 .

 12 Figure 1. Depicted is a punch with the left arm from early to later points-in-time during movement execution. Figure 2. Sequence of a trial, exemplarily shown for the kick with the left leg.

Figure 3 .

 3 Figure 3. Reaction times (RT) in milliseconds (± SE) for combat sports athletes and novices at all points-in-time during movement execution.

Figure 4 .

 4 Figure 4. Error rates (ER) in percentage (± SE) for combat sports athletes and novices at all points-in-time during movement execution.

Figure 5 .

 5 Figure 5. Fixation durations in milliseconds (± SE) for combat sports athletes and novices at all points-in-time during movement execution.

Figure 6 .

 6 Figure 6. Fixation rates in fixations per second (± SE) for combat sports athletes and novices at all points-in-time during movement execution.

ACKNOWLEDGMENTS

This work was supported by a grant of the German Research Foundation (DFG) awarded to the second author (GU 1683/1-1) of the manuscript.


Marco Lorenzi
email: marco.lorenzi@inria.fr

Maurizio Filippone
email: maurizio.filippone@eurecom.fr

Giovanni B Frisoni
email: giovanni.frisoni@unige.ch

Daniel C Alexander
email: d.alexander@ucl.ac.uk

Sébastien Ourselin
email: s.ourselin@ucl.ac.uk

Sebastien Ourselin

Probabilistic disease progression modeling to characterize diagnostic uncertainty: application to staging and prediction in Alzheimer's disease

Keywords: Alzheimer's disease, Diagnosis, Disease progression modeling, Gaussian process, Clinical trials

Introduction

Neurodegenerative disorders (NDDs), such as Alzheimer's disease (AD), are characterised by the progressive pathological alteration of the brain's biochemical processes and morphology, and ultimately lead to the irreversible impairment of cognitive functions [START_REF] Brookmeyer | Forecasting the global burden of Alzheimer's disease[END_REF]. The correct understanding of the relationship between the different pathological features is of paramount importance for improving the identification of pathological changes in patients, and for better treatment [START_REF] Jack | Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade[END_REF].

To this end, ongoing research efforts aim at developing precise models allowing optimal sets of measurements (and combinations of them) to uniquely identify pathological traits in patients. This problem requires the definition of optimal ways to integrate and jointly analyze the heterogeneous multi-modal information available to clinicians [START_REF] Young | Accurate multimodal probabilistic prediction of conversion to Alzheimer's disease in patients with mild cognitive impairment[END_REF][START_REF] Mwangi | A review of feature reduction techniques in neuroimaging[END_REF][START_REF] Lorenzi | Multimodal image analysis in Alzheimer's disease via statistical modelling of non-local intensity correlations[END_REF]. By consistently analyzing multiple biomarkers that to date have mostly been considered separately, we aim at providing a richer description of the pathological mechanisms and a better understanding of individual disease progressions.

Disease progression modeling (DPM) is a relatively new research direction

for the study of NDD data [START_REF] Fonteijn | An event-based disease progression model and its application to familial Alzheimer's disease[END_REF][START_REF] Jedynak | A computational neurodegenerative disease progression score: Method and results with the Alzheimer's disease neuroimaging initiative cohort[END_REF][START_REF] Donohue | Estimating longterm multivariate progression from short-term data[END_REF][START_REF] Younes | Inferring changepoint times of medial temporal lobe morphometric change in preclinical Alzheimer's disease[END_REF][START_REF] Bilgel | Temporal trajectory and progression score estimation from voxelwise longitudinal imaging measures: Application to amyloid imaging[END_REF][START_REF] Schiratti | A mixed-effects model with time reparametrization for longitudinal univariate manifold-valued data[END_REF][START_REF] Guerrero | Instantiated mixed effects modeling of Alzheimer's disease markers[END_REF][START_REF] Marinescu | A vertex clustering model for disease progression: Application to cortical thickness images[END_REF]. The main goal of DPM consists in revealing the natural history of a disorder from collections of imaging and clinical data by: 1) quantifying the dynamics of NDDs along with the related temporal relationship between different biomarkers, and 2) staging patients based on individual observations for diagnostic and interventional purposes. Therefore, this research domain is closely related to the exploitation of advanced statistical/machine-learning approaches for the joint modelling of the heterogeneous and information available to clinicians: imaging, biochemical, and clinical biomarkers. Differently from the several predictive machine-learning approaches proposed in the past in NDD research, disease progression models aim at explicitly estimating the temporal progression of the biomarkers from normal to pathological stages, to provide a better interpretation and understanding of the natural evolution of the pathology. For this reason it represents a very appealing modeling approach in clinical settings.

The main challenge addressed by DPM consists in the general lack a welldefined temporal reference in longitudinal clinical dataset of NDDs. Indeed, age or visit date information are biased time references for the individual longitudinal measurements, since the onset of the pathology may vary across individuals according to genetic and environmental factors [START_REF] Yang | Quantifying the pathophysiological timeline of Alzheimer's disease[END_REF]. This is a very specific methodological issue requiring the extension and generalization of the analysis approaches classically used in time-series analysis.

To tackle this problem, it is usually assumed that individual biomarkers are measured relatively to an underlying disease trajectory defined with respect to an absolute time axis describing the natural history of the pathology [START_REF] Jedynak | A computational neurodegenerative disease progression score: Method and results with the Alzheimer's disease neuroimaging initiative cohort[END_REF]. Each individual is thus characterized by a specific observation time that needs to be estimated in order to assess the individual pathological stage. According to this statistical setting, we therefore aim at estimating a group-wise disease model defined with respect to an absolute time scale, along with individual time re-parameterisation relative to the group-wise progression. This modeling paradigm has been implemented in a number of approaches proposed in the recent years, either by assuming continuous temporal trajectories of the biomarkers [START_REF] Jedynak | A computational neurodegenerative disease progression score: Method and results with the Alzheimer's disease neuroimaging initiative cohort[END_REF][START_REF] Donohue | Estimating longterm multivariate progression from short-term data[END_REF][START_REF] Younes | Inferring changepoint times of medial temporal lobe morphometric change in preclinical Alzheimer's disease[END_REF][START_REF] Bilgel | Temporal trajectory and progression score estimation from voxelwise longitudinal imaging measures: Application to amyloid imaging[END_REF][START_REF] Schiratti | A mixed-effects model with time reparametrization for longitudinal univariate manifold-valued data[END_REF][START_REF] Guerrero | Instantiated mixed effects modeling of Alzheimer's disease markers[END_REF][START_REF] Marinescu | A vertex clustering model for disease progression: Application to cortical thickness images[END_REF], or by modeling the disease progression as a sequence of discrete events [START_REF] Fonteijn | An event-based disease progression model and its application to familial Alzheimer's disease[END_REF][START_REF] Young | A data-driven model of biomarker changes in sporadic Alzheimer's disease[END_REF]. For example, in [START_REF] Donohue | Estimating longterm multivariate progression from short-term data[END_REF] the authors proposed to model the temporal biomarker trajectories through random effect regression, building on the theory of selfmodeling regression [START_REF] Kneip | Convergence and consistency results for self-modeling nonlinear regression[END_REF], while the authors of [START_REF] Schiratti | A mixed-effects model with time reparametrization for longitudinal univariate manifold-valued data[END_REF] re-frame the random effect regression model in a geometrical setting, based on the assumption of a logistic curve shape for the average biomarker trajectories.

Continuous progression models have been recently extended to the modelling of brain images based on the time-reparameterization of voxel/mesh-based measures [START_REF] Younes | Inferring changepoint times of medial temporal lobe morphometric change in preclinical Alzheimer's disease[END_REF][START_REF] Bilgel | Temporal trajectory and progression score estimation from voxelwise longitudinal imaging measures: Application to amyloid imaging[END_REF][START_REF] Marinescu | A vertex clustering model for disease progression: Application to cortical thickness images[END_REF].

The use of disease progression models for diagnostic purposes is instead less investigated. Predictive models of patient staging were proposed within the setting of the Event Based Model [START_REF] Fonteijn | An event-based disease progression model and its application to familial Alzheimer's disease[END_REF], or still through random effect modeling [START_REF] Guerrero | Instantiated mixed effects modeling of Alzheimer's disease markers[END_REF]. However, the Event Based Model relies on the coarse binary discretization of the biomarker changes, and does not account for longitudinal observations, while the predictive models proposed in [START_REF] Guerrero | Instantiated mixed effects modeling of Alzheimer's disease markers[END_REF] and [START_REF] Schmidt-Richberg | Multi-stage biomarker models for progression estimation in Alzheimer's disease[END_REF] require cohorts with known disease onset, and therefore lack flexibility while being prone to bias due to misdiagnosis and uncertainty of the conversion time.

Furthermore, these methods are generally not formulated in a probabilistic setting, which makes it difficult to account for uncertainties in biomarker progressions and diagnostic predictions. Indeed, the quantification of the variability associated with the biomarkers trajectories, as well as the assessment of the diagnostic uncertainty in de-novo patients, are crucial requirements for decision making in clinical practice [START_REF] Shinkins | Diagnostic uncertainty: dichotomies are not the answer[END_REF].

Nonetheless, the ensemble of this research offers a sight of the potential of these approaches in representing a novel and powerful diagnostic instrument: in this study we thus aim at assessing the ability of DPM in providing a statistical reference for the transition from normal to pathological stages, for probabilistic diagnosis in the clinical scenario. To this end, we reformulate classical DPM within a Bayesian setting in order to allow the probabilistic estimate of the biomarker trajectories and the quantification of the uncertainty of predictions of the individual pathological stage. The resulting probabilistic framework is exploited in an hypothetical clinical scenario, for the estimation of the pathological stage in a de-novo cohort of testing individuals, by assessing the influence of missing observations, biomarkers, and follow-up information.

The manuscript is structured as follows. Section 2.1 formulates DPM based on Bayesian Gaussian Process regression [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], while Section 2.2 illustrates the validation of our model on clinical and multivariate imaging measurements from a cohort of 782 amyloid positive individuals extracted from the ADNI database.

Methods

Statistical setting

This section highlights the statistical framework employed in this study, based on the reformulation of self-modeling regression withing a Bayesian setting. This achieved by 1) defining a random effect Gaussian process regression model to account for individual correlated time series (section 2.1); 2) modeling individual time transformations encoding the information on the latent pathological stage (section 2.1.2); and 3) introducing a monotonicity information in order to impose a regular behaviour on the biomarkers trajectories (section 2.1.3). We finally illustrate in section 2.1.4 how the proposed framework leads to a probabilistic model of disease staging in de-novo individuals, naturally accounting for missing information. Further details on model specification and inference are provided in the Supplementary Section AppendixA.1, while the experimental validation on synthetic data is reported in Supplementary Section AppendixA.2. We represent the longitudinal biomarker's measures associated with each individual j as a multidimensional array (y j (t 1 ), y j (t 2 ), . . . , y j (t k j )) sampled at k j multiple time points t = {t 1 , t 2 , . . . , t k j }. Although different biomarkers may be in reality sampled at different time-points, for the sake of notation simplicity in what follows we will assume, without loss of generality, that the sampling time is common among them. The observations for individual j at a single time point t are thus a random sample from the following generative model:

y j (t) = y j b1 (t), y j b2 (t), . . . , y j b N b (t) (1)
= f (t) + ν j (t) + , (2)
where

f (t) = (f b1 (t), f b2 (t) . . . , f b N b (t)) is the fixed effect function modelling the biomarker's longitudinal evolution, ν j (t) = (ν j b1 (t), ν j b2 (t), . . . , ν j b N b (t))
is the individual random effect, and = ( b1 , b2 , . . . , b N b ) is time-independent observational noise. The group-wise evolution is modelled as a GP, f ∼ GP(0, Σ G ), the individual random effects are assumed to be correlated perturbations ν j ∼ N (0, Σ S ), while the observational noise is assumed to be a Gaussian heteroskedas-

tic term ∼ N (0, Σ ), where Σ is a diagonal matrix diag[σ 2 b1 , σ 2 b2 , . . . , σ 2 b N b
].

Fixed Effect Process

The covariance function Σ G describes the biomarkers temporal variability, and is represented as a block-diagonal matrix

Σ G (f , f ) = diag[Σ b1 (f b1 , f b1 ), Σ b2 (f b2 , f b2 ), . . . , Σ b N b (f b N b , f b N b )],
where each block represents the withinbiomarker temporal covariance expressed as a negative squared exponential

function Σ b (f b (t 1 ), f b (t 2 )) = η b exp -(t1-t2) 2 2 l 2 b
, and where the parameters η b and l b are the marginal variance and length-scale of the biomarker's temporal evolution, respectively.

Individual Random Effects

The random covariance function Σ S models the individual deviation from the fixed effect, and is represented as a block-diagonal matrix Σ S (ν j ,

ν j ) = diag[ Σ j b1 (ν j b1 , ν j b1 ), Σ j b2 (ν j b2 , ν j b2 ), . . . , Σ b j N b (ν j b N b , ν j b N b
)], where each block Σ j b corresponds to the covariance function associated with the individual process ν j b (t). Thanks to the flexibility of the proposed generative model, any form of the random effect covariance Σ S can be easily specified in order to model the subject-specific biomarkers' progression. In what follows we will use a linear

covariance form Σ j b (ν j b (t 1 ), ν j b (t 2 )) = (σ j b ) 2 (t 1 -t)(t 2 -t)
, where t is the average observational time for individual j, when more than 4 measurements are available, and i.i.d. Gaussian covariance form Σ j b (ν j b (t 1 ), ν j b (t 2 )) = (σ j b ) 2 when 2 or 3 measurements are available, while assigning it to 0 otherwise (thus by accounting only for the observational noise σ 2 b ). This choice is motivated by stability concerns, in order to keep the model complexity compatible with the generally limited number of measurements available for each individual.

Individual time transformation

The generative model ( 1) is based on the key assumption that the longitudinal observations across different individuals are defined with respect to the same temporal reference. This assumption may be invalid when the temporal alignment of the individual observations with respect to the common group-wise model is unknown, for instance in the typical scenario of a clinical trial in AD where the patients' observational time is relative to the common baseline, and

where the disease onset is a latent event (past or future) which is not directly measurable. This modeling aspect is integrated by assuming that each individual measurement is made with respect to an absolute time-frame τ through a time-warping function t = φ j (τ ) that models the time-reparameterization with respect to the common group-wise evolution. Model (1) can thus be reparame-terized as

y j (φ j (τ )) = f (φ j (τ )) + ν j (φ j (τ )) + . (3)
The present formulation allows the specification of any kind of time transformation, and in what follows we shall focus on the modelling of a linear reparameterization of the observational time φ j (τ ) = τ + d j . This modeling assumption is mostly motivated by the choice of working with a reasonably limited number of parameters, compatibly with the generally short follow-up time available per individual (cfr. Table 2). Within this setting, the time-shift d j encodes the disease stage associated with the individual relatively to the group-wise model.

Monotonic constraint in random-effect multimodal GP regression

Due to the non-parametric nature of Gaussian process regression, we need an additional constraint on model (3) in order to identify a unique solution for the time transformation. By assuming a steady temporal evolution of biomarkers from normal to pathological values, we shall assume that the biomarker trajectories described by (3) follow a (quasi) monotonic behaviour. This requirement can be implemented by imposing a prior positivity constraint on the derivatives of the GP function. Inspired by [START_REF] Riihimäki | Gaussian processes with monotonicity information[END_REF], we impose a monotonicity constraint by assuming a probit-likelihood for the derivative measurements m(t) associated with the derivative process ḟ (t) = df (t) dt at time t:

p(m(t)| ḟ (t)) = Φ 1 λ ḟ (t) , (4)
with Φ(z) = k=1 to the value of 1e-6. The position of the derivative points was updated at each iteration, according to the changes of the GP domain.

By following a similar construction, we could equally enforce a monotonic behavior to the random effects associated with the individual trajectories. This additional constraint would however come with a cumbersome increase of the model complexity, since it would introduce an additional layer of virtual derivative parameters (with associated location) per individual. Moreover, while we are interested in modeling a globally monotonic biomarker trajectory on the fixed parameters, we relax this constraint at the individual level, since some subjects may be characterised by non strictly monotonic time-series due to specific clinical conditions.

Model likelihood and parameters

Given the sets of individual biomarker measurements y = {(y j (t i ))

k j i=1 } N j=1 ,
and of D control derivatives m = {m b k (t l )} D l=1 at points t = {t l } D l=1 for the progression of each biomarker b k , the random effect GP model posterior is:

p f , ḟ , ν j |y, m = 1 Z p(f , ḟ |t, t )p(ν|t)p(y|f , ν)p(m| ḟ ) = p(f , ḟ |t, t )p(ν|t)p(y|f , ν) k l Φ 1 λ ḟb k (t l ) , (5)
where ν = {ν j } N j=1 . Due to the non-Gaussianity of the derivative term Φ, the direct inference on the posterior is not possible due to its analytically intractable form. For this reason, we employ an approximate inference scheme based on classical approaches to Gaussian process with binary activation functions [START_REF] Nickisch | Approximations for binary Gaussian process classification[END_REF] (AppendixA.1).

Overall, model ( 3) is identified by (N j + 3)N b + N j parameters, represented by the fixed effects and noise

θ G = {η b k , l b k , b k } N b k=1 , by the individual random effects parameters θ j G = {σ j b k } N b k=1
and by the time-shifts d j . In what follows, the optimal parameters are obtained by maximising the approximated log-marginal likelihood derived from the posterior (5) through conjugate gradient descent, via alternate optimization between the hyper-parameters θ G and θ j G , and the individuals' time-shifts d j . Regularization is also enforced by introducing Gaussian priors for the parameters θ G and θ j G .

Prediction of observations and individual staging

Gaussian processes naturally allow for probabilistic predictions given the observed data. At any given time point t * , the posterior biomarker distribution has the Gaussian form p(f * |t * , y, t, m, t ) ∼ N (f * |µ * , Σ * ):

µ * = Σ G (f (t * ), f (t))(Σ joint + Σjoint ) -1 μjoint (6)
Σ * = Σ G (f (t * ), f (t * )) -Σ G (f (t * ), f (t)) (Σ joint + Σjoint ) -1 Σ G (f (t), f (t * )), (7)
where the matrix Σ joint + Σjoint is the joint covariance resulting from the inference scheme detailed in Supplementary Section AppendixA.1 [START_REF] Riihimäki | Gaussian processes with monotonicity information[END_REF].

We also derive a probabilistic model for the individual temporal staging given a set of biomarker observations y * , thanks to the Bayes formula:

p(t * |y * , y, t, m, t ) = p(y * |t * , y, t, m, t ) p(t * )/p(y * |y, t, m, t ), (8)
which we compute by assuming an uniform distribution on t * , and by noting

that p(y * |t * , y, t, m, t ) ∼ N (µ * , Σ * + Σ ).
In particular, the covariance form

Σ G (f (t * ), f (t * )
) can be specified in order to account for incomplete data, and thus generalizes the GP model for predictions in presence of missing biomarker observations. The posterior distribution (8) quantifies the confidence of the model about the individual disease staging, and thus is a valuable information about the precision of the diagnosis. We will also compute the expectation of the distribution p(t * |y * , y, t, m, t ), which provides a scalar value that can be used in subsequent classification methods.

Meterials and Methods

Study Participants

Data used in the preparation of this article were obtained from the ADNI database (http://adni.loni.usc.edu). The ADNI was launched in 2003 as a public-private partnership, led by Principal Investigator Michael W. Weiner, MD. The primary goal of ADNI has been to test whether serial magnetic resonance imaging, positron emission tomography, other biological markers, and clinical and neuropsychological assessment can be combined to measure the progression of MCI and early AD. For up-to-date information, see www.adni-info. org.

Data Processing

We collected longitudinal measurements for the ADNI individuals with baseline values of cerebrospinal fluid (CSF) Aβ amyloid lower than the nominal val- 3.4 (11) 3.9 (0) 3.9 (0) 3.9 (0) 3.8 (43) 3 [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] Table 2: Average follow-up years and percentage of individuals with missing data (in parenthesis).

Model testing on de-novo individuals

The estimated probabilistic disease progression model provides a valuable clinical reference, as it can be used to predict an individual pathological stage, as well as to quantify the biomarkers predictive value, or the influence of missing data. To this end, we estimated the predictive performance of the model in assessing the individual pathological stage with respect to follow-up assessments and missing biomarkers. This was done by estimating the predictive accuracy of the group-wise separation obtained via increasing thresholds of the estimated temporal progression.

Results

Model plausibility

The estimated biomarker progression (Figure 1-A) shows a biologically plausible description of the pathological evolution, compatible with previous findings in longitudinal studies in familial AD [START_REF] Bateman | Clinical and biomarker changes in dominantly inherited Alzheimer's disease[END_REF], and with the hypothetical models of AD progression [START_REF] Jack | Hypothetical model of dynamic biomarkers of the Alzheimer's pathological cascade[END_REF][START_REF] Frisoni | The clinical use of structural mri in Alzheimer disease[END_REF]. The progression is defined on a time scale spanning Finally, when applying [START_REF] Donohue | Estimating longterm multivariate progression from short-term data[END_REF] to the training data we measured a strong agreement between the resulting progression and the one obtained with our method, resulting in a correlation between the corresponding individual time-shifts of 0.94 (p <1e-6) (Supplementary Material AppendixA.2.2).

Assessing diagnostic uncertainty in testing data: an illustrative example.

This section illustrates the use of the model represented in Figure 1 for the quantification of diagnostic uncertainty in testing individuals. We consider the hypothetical scenario where the baseline values for different biomarkers are measured for a given patient, namely FAQ, hippocampal and ventricle volumes.

We assume that the biomarkers values correspond to the 20th percentiles with

350

We now suppose that for the same patient we acquire a follow-up measure-ment for each biomarker at year 1, with values corresponding to the 50th percentiles of the distribution of the training set (i.e. FAQ = 5, normalized Hippo = 4.3e-3, normalized Ventr = 2.7e-2). The right hand side of Figure 2 shows the new prediction based on the joint baseline+follow-up information. For each biomarker the posterior distributions indicate an increase of the predicted disease stage with respect to the baseline scenario, while the prediction uncertainty is generally lower. Although the expectation for the 3 biomarkers is very similar (t F AQ =-2.5, t hippo =-3.5, and t V entr =-3.2), we notice that FAQ leads to the highest diagnostic confidence (sd F AQ =2.6), followed by hippocampal volume (sd Hippo =3.8), and finally by ventricles volume (sd V entr =5.7). Further assessment of the relationship between biomarker variability and model prediction is provided in supplementary Section AppendixA.2.3.

This illustrative example shows that the proposed probabilistic framework represents a valuable instrument for the assessment of the diagnostic value and uncertainty associated with different biomarkers, and can faithfully track the pathological progression of testing individuals along the modeled trajectories, from normal to pathological levels.

DPM for probabilistic diagnosis in ADNI.

We now assess the predictive results of the model when applied to the testing ADNI cohort. Figure 3 shows the individual posterior predictive distributions associated with the testing individuals, and the boxplot of the expected timeshift when using the model as statistical reference through formula (Figure 8).

The figure reports the two different modeling scenarios based on baseline information only (Figure 3-1), and on the complete set of baseline and longitudinal measurements (Figure 3-2). We first note that the group-wise differences between the expected time-shifts are compatible for both scenarios, as shown by the similar boxplot distributions across groups (Figure 3-1b vs 3-2b). The con- As with the training case, for both scenarios the group-wise distribution of the expected time-shift shows a significant separation between the clinical groups according to the increase of the pathological stage (ANOVA, p <1e-4).

Interestingly, the temporal positioning of the non converting MCI lies between controls and MCI converters, and is on average lower than the one of healthy individuals subsequently converted to cognitive impairment. We further tested the model in presence of missing information, by computing the predictions when only one baseline biomarker is available (Figure 5). The predictive outcomes show important variations depending on the con- sidered biomarker, while the confidence bounds for the predictions are usually large, to denote increased uncertainty. We also note that FAQ, ADAS13, and hippocampal volume are the biomarkers leading to the largest group-wise separation, along with the lowest prediction uncertainty. This aspect is quantified in Table 3, reporting the discrimination results with respect to the nominal cut-off For the sake of comparison we finally benchmarked the predictive results

point
provided by the disease progression model with respect those obtained by the classification analysis performed with standard statistical tools, such as a random forest classifier. We note that the comparison of the classification perfor-mance obtained on the heterogeneous data considered in this work is generally challenging, since the proposed DPM 1) accounts for missing data and nonfixed number of time points per individuals, and 2) is formulated in order to consistently handle both longitudinal and cross sectional measurements, either

for training and prediction. To date there is not a consensus on the optimal approach to adopt to tackle these important modelling constraints, while the comparison between the classification performance obtained with complex machine learning methods is currently matter of scientific debate and investigation [START_REF] Mendelson | Selection bias in the reported performances of AD classification pipelines[END_REF].

For this reason we restricted the random forest classification task to a standard statistical setting, in order to essentially provide a reliable benchmark for the classification performance of the proposed disease progression. To this end we trained the random forest on the classification between healthy individuals and AD patients based on the baseline measurements of the training group, while the missing entries in the testing data were imputed via nearest neighbour search, based on the available biomarkers. The classification results are reported in Supplementary Table A.5.

The performance of the random forest classifier is generally inferior to the one obtained with the proposed approach, as witnessed by the consistently lower AUC obtained for all the comparisons. The difference becomes more evident for the more challenging classification problems, such as the identification of conversion in MCI and healthy individuals. This result is indicative of the reliability of the classification results obtained by the proposed disease progression model, especially when considering that the random forest classifier is explicitly optimized to maximize the separation between groups, while the accuracy results reported in Table 3 are based on the empirical choice of a reference threshold in the training population.

DPM staging and chronological age.

We finally compare the relationship between the predicted disease staging in training and testing set and the individual chronological age. We first note that both training and testing clinical groups were matched by age, with the exception of the 5 training healthy subjects converted to MCI (or AD) that were slightly older with respect to the reference training healthy population (p=0.02).

Nevertheless, when comparing the estimated time shift with respect to the chronological age of each individual we didn't report any significant correlation between these measures. Interestingly, the same lack of association is also quantifiable in the testing group (Figure 6). This result, in association with the strong relationship between time shift and clinical condition reported in Section 3.3, let us conclude that the model is describing the biomarker's variation essentially related to the pathological progression, which is orthogonal to the effect of healthy aging quantified by the chronological age. This result points to the effectiveness of the proposed approach in capturing significant effects related to the specific temporal progression of the disease.

Discussion

This progression. We also note that the time shift is a relative measure of disease progression accounting for the biomarker variability observed in the training population. Thus, the point 0 is generally not associated with the conversion to AD, as it is relative to the data initialization (in this case the study baseline).

We illustrated the use of DPM as benchmarking tool for the statistical comparison of biomarkers. The model allows the quantification of the variability associated with the single biomarkers, by identifying the related uncertainty in characterizing the progression from normal to pathological levels. The proposed model can be thus used as a reference for screening and enrichment purposes in clinical trials [START_REF] Lorenzi | Enrichment through biomarkers in clinical trials of Alzheimer's drugs in patients with mild cognitive impairment[END_REF][START_REF] Yu | Operationalizing hippocampal volume as an enrichment biomarker for amnestic mild cognitive impairment trials: effect of algorithm, test-retest variability, and cut point on trial cost, duration, and sample size[END_REF][START_REF] Hill | Coalition against major diseases/european medicines agency biomarker qualification of hippocampal volume for enrichment of clinical trials in predementia stages of Alzheimer's disease[END_REF].

The modelled progression showed that neuropsychological tests generally lead to lower uncertainty for identifying the individual clinical stage, and to the higher separation power between healthy and AD groups. This finding is compatible with the results reported by previous disease progression models applied to ADNI, such as [START_REF] Jedynak | A computational neurodegenerative disease progression score: Method and results with the Alzheimer's disease neuroimaging initiative cohort[END_REF] and [START_REF] Young | A data-driven model of biomarker changes in sporadic Alzheimer's disease[END_REF]. In this latter study ADAS13 consistently appeared among the first events distinguishing the normal disease stages from the pathological ones. Furthermore, our analysis further showed that volumetric measures such as hippocampal and entorhinal volumes provide equivalent if not superior diagnostic performances when tested on the more challenging problem of detecting conversion to dementia from healthy and MCI stages, especially in terms of improved AUC. Nevertheless, some care should be taken in drawing conclusions from the present analysis. Our model was based on the standard volumetric measures provided in the ADNI database, and we cannot exclude that a more precise quantification of morphological brain changes would lead to even better performance of volumetric biomarkers [START_REF] Wolz | Measurement of hippocampal atrophy using 4D graphcut segmentation: application to ADNI[END_REF][START_REF] Cash | Assessing atrophy measurement techniques in dementia: Results from the MIRIAD atrophy challenge[END_REF]. Furthermore, the proposed model was not optimized in order to maximize the classification accuracy between clinical groups. For example, the results reported in Table 3 are based on the choice of the temporal threshold corresponding to the reference value of the 15th percentile of the AD distribution. This cut-off was not optimized to maximize the predictive outcome of the biomarkers, but was rather chosen based on heuristics aimed at illustrating the use of the model for predictive purposes. We thus cannot exclude that the optimization of the temporal threshold would lead to different figures for the classification task. The reported results are therefore indicative of the effectiveness of the model in faithfully representing the clinical spectrum of the disease. We note also that the reported figures are in line with those provided by state-of-art methods in AD classification, without requiring complex parameter optimization procedures, which would introduce additional levels of cross-validation and expose the results to selection bias and generalization issues [START_REF] Mendelson | Selection bias in the reported performances of AD classification pipelines[END_REF].

Thanks to the probabilistic formulation we showed that the use of longitudinal information is important for reducing the uncertainty of the prediction, and thus allowing one to better identify the disease status associated to an individual. This important aspect is in agreement with the generally higher statistical power reported in previous Alzheimer's studies comparing longitudinal measurements to baselines ones [START_REF] Henneman | Hippocampal atrophy rates in Alzheimer's disease. added value over whole brain volume measures[END_REF][START_REF] Frisoni | The clinical use of structural mri in Alzheimer disease[END_REF][START_REF] Xu | Longitudinal analysis is more powerful than crosssectional analysis in detecting genetic association with neuroimaging phenotypes[END_REF].

In this work we focused on the modelling of the progression of amyloid positive individuals. This choice was motivated by the interest in assessing the model performance on an homogeneous clinical population likely to be representative of the Alzheimer's evolution. While the absence of pathological amylod levels seems indicative of non-Alzheimer's pathophysiology [START_REF] Gordon | Longitudinal β-amyloid deposition and hippocampal volume in preclinical Alzheimer's disease and suspected non-Alzheimer disease pathophysiology[END_REF][START_REF] Mormino | Heterogeneity in suspected non-Alzheimer disease pathophysiology among clinically normal older individuals[END_REF], there is currently an active debate on the mechanisms of neurodegeneration not related to brain amyloidosis [START_REF] Jack | Suspected non-Alzheimer disease pathophysiology -concept and controversy[END_REF]. The investigation of these aspects goes beyond the scope of the present work, and future extensions of disease progression modeling will aim at identifying differential progressions underlying sub-pathologies, for example by reformulating the proposed random effect regression within the realm of Gaussian process mixture models [START_REF] Lázaro-Gredilla | Overlapping mixtures of Gaussian processes for the data association problem[END_REF][START_REF] Ross | Nonparametric mixture of Gaussian processes with constraints[END_REF]

Methodological considerations

From the methodological perspective, we proposed a novel probabilistic approach based on Gaussian process regression for disease progression modeling from time-series of biomarker measurements enabling novel applications beyond the state-of-art, such as the probabilistic prediction of disease staging in testing individuals. Furthermore, the model naturally accounts for missing data, and provides uncertainty quantification of the biomarker evolutions. Similarly to [START_REF] Donohue | Estimating longterm multivariate progression from short-term data[END_REF], in this work we focused on the modeling of disease staging represented by a time shift, although the proposed framework can naturally account for more complex time transformations, provided that a sufficient number of time points is available for each individual.

From the methodological point of view, the proposed model extends current approaches to GP-regression by consistently integrating time-reparameterization and monotonic constraints within a random effect regression framework. Monotonic GPs were introduced in [START_REF] Riihimäki | Gaussian processes with monotonicity information[END_REF] as a principled regularization solution to im-prove the plausibility of modeling results. For example, the strength of such a regularization approach in biomedical applications has been illustrated in survival analysis [START_REF] Joensuu | Risk of recurrence of gastrointestinal stromal tumour after surgery: an analysis of pooled population-based cohorts[END_REF]. Our approach extends this framework by consistently integrating a latent time variable parameter within a random effect model formulation.

The idea of estimating a time transformation in a GP regression framework has been previously used by [START_REF] Liu | Estimating replicate time shifts using Gaussian process regression[END_REF] to account for uncertain measurement times to a microarray dataset of mRNA. However, in that work the estimation of the time uncertainty was subject to a strong prior constraint based on the assumption that the unknown biological time must be similar to the measured one. In the application proposed in our work such an assumption is no longer valid and would ultimately lead to implausible estimations. On the contrary, the proposed GP regression is able to recover the underlying time transformation thanks to the proposed monotonicity regularization.

Finally, thanks to the flexibility of the proposed Gaussian process framework, further extensions of the model will enable to consistently integrate a spatio-temporal covariance model, such as the efficient Kronecker form of [START_REF] Lorenzi | Efficient Gaussian process-based modelling and prediction of image time series[END_REF],

to provide a unified framework for jointly modelling time series of images and scalar biomarkers data in a coherent fully Bayesian setting.

Conclusions

This work proposes an extension of DPM for the accurate quantification of the diagnostic uncertainty in Alzheimer's disease. The proposed application shows that DPM provides at the same time a plausible description of the transition from normal to pathological stages along the natural history of the disease, as well as remarkable diagnostic performances when tested on de-novo individuals. The model used in this study can account for any missing data patterns (longitudinal or across biomarkers), and allows to directly quantify the uncertainty related to the missing information. It thus represents a novel and promising tool for the analysis of clinical trials data.

Further Information

The open-source code as well as the proposed predictive model trained on ADNI data will be available at the author's web-page: https://team.inria.

fr/asclepios/marco-lorenzi/. The realization of this study required about 1.5kWh of computing power.

with f joint =    f ḟ    distributed as N       0 0    ,    Σ G (f (t), f (t)) ∂Σ G (f (t),f (t )) ∂t dΣ G (f (t ),f (t)) dt d 2 Σ G (f (t ),f (t )) dt 2       .

AppendixA.1.1. Approximated inference

Due to the non-Gaussianity of the derivative likelihood term, the direct inference on the posterior (A.1) is not possible due to its analytically intractable form. For this reason, we employ an approximate inference scheme based on clas-As shown in [START_REF] Riihimäki | Gaussian processes with monotonicity information[END_REF] for univariate monotonic regression, moments and updates of the approximation parameters can be computed in an analogous manner as in the classical GP classification problem [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF].

AppendixA.1.3. Marginal Likelihood and hyper-parameter estimation

The model's log-marginal likelihood under the EP approximation is:

log L = - 1 2 log |Σ joint + Σjoint | - 1 2 μT joint (Σ joint + Σjoint ) -1 μjoint + k l (µ -kl -μkl ) 2 2(σ 2 -kl ) + σ2 kl ) + k l log Φ( µ -kl λ 2 k + σ 2 -kl ) ) + 1 2 k l log(σ 2 -kl + σ2 kl ). (A.4)
In what follows, the optimal parameters are obtained by maximising log L through conjugate gradient descent, via alternate optimization between the hyper-parameters θ G and θ j G , and the individuals' time-shifts d j . The position of the derivative points was updated at each iteration, according to the changes of the GP domain. Regularisation was also enforced by introducing Gaussian priors for the parameters θ G and θ j G . We note that the block structure of the GP covariance allows the computation of the gradients with respect to the biomarkers' and individual parameters by working on matrices of much smaller dimension than the one of the whole GP, thus considerably improving the numerical stability and the computational efficiency of the optimization procedure.

2. 1 . 1 .

 11 Gaussian process-based random effect modeling of longitudinal progressions In what follows, longitudinal measurements of N b biomarkers {b 1 , . . . , b N b } over time are given for N individuals.

z-∞

  N (x|0, 1) dx. The quantity λ > 0 is an additional model parameter controlling the degree of positivity enforced on the derivative process, with values approaching zero for stronger monotonicity constraint. In what follows, the monotonicity of each biomarker is controlled by placing 10 derivative points equally spaced on the observation domain, and by fixing the N b derivative parameters {λ b k } N b

Figure 1 :

 1 Figure 1: A) Modelled biomarker progression in the training set of 200 Aβ amyloid positive individuals (solid/dashed lines: mean ± sd). B) Posterior prediction for the individual time shift in training data (top: individual time-shift distribution; bottom: group-wise boxplot of the expected time-shift). Healthy individuals are generally displaced at the early stages of the pathology, while the predictions for MCI and AD patients are associated with respectively intermediate and late progression stages. NL: normal individuals, MCI: mild cognitive impairment, AD: Alzheimer's patients.

  roughly 20 years, and is characterized at the initial stages by high-levels of AV45, followed by the abnormality of ventricles volume, of FDG uptake, and of the whole brain volume. These latter measures are however heterogeneously distributed across clinical groups, and with rather large variability. The evolution is further characterized by increasing abnormality of the volumetric measures (especially hippocampal volume), and by the steady worsening of neuropsychological scores such as FAQ. The model thus shows that the transition from normal to pathological levels is essentially characterized by increase of hypometabolism, followed by the pronounced temporal brain atrophy. Moreover, the worsening of the neuropsychological and functional scores closely (almost linearly) follows the progression in the advanced clinical stages. The joint visualization of the temporal progression of the biomarkers with temporal derivative of the modelled average trajectories is shown in Supplementary figure A.10.The illustration confirms that ADAS13 and FAQ are characterised by very sim-ilar longitudinal profiles, and show the largest changes in the latest stages of the pathology (peak of the derivative at t>0). On the contrary, the change in hippocampal volume is more strongly associated with the earlier stages of the pathology. AV45 and ventricles volumes are the least informative and are associated with the lowest changes.

Figure 1 -

 1 Figure 1-B (top) shows the posterior time-shift distributions associated with the individuals. The distributions denote the confidence of the model in associating to each individual a temporal staging with respect to the global pathological progression. The boxplot of Figure 1-B (bottom) reports the group-wise expectation of the individual time-shifts. Healthy individuals (blue) are associated with the early stages of the pathology in both training and testing data, while MCI (purple) and AD patients (red) are characterized by respectively intermediate and late predicted progression stages. The group-wise comparison between the expected time-shifts was statistically significant between each group pairs (ANOVA, p <1e-6). Moreover, the time to conversion to AD in the MCI group was significantly correlated with the disease staging quantified by the expectation of the individual time distributions (R 2 = -0.4, p = 3.8e -4).

Figure 2 :

 2 Figure 2: Illustrative example: posterior prediction of disease staging for a testing individual based on baseline (left), and baseline + follow-up (right) information for three biomarkers: FAQ, hippocampal and ventricles volume. The biomarkes values correspond to the 20th and 50th percentiles of the training-group distribution for respectively baseline and follow-up measures. Adding the follow-up information leads to increased estimates of the disease staging and to generally lower prediction uncertainty. Although the distributions associated with different biomarkers generally lead to similar expectations, FAQ and hippocampal volume lead to the lowest diagnostic uncertainty. Vertical lines: expectation for each posterior distribution.

Figure 3 :

 3 Figure 3: Posterior prediction for the individual time shift in testing data by using i) only the baseline information (1a-b), and ii) the baseline + follow-up information available for each testing subject (2a-b). Healthy individuals are generally displaced at the early stages of the pathology, while the predictions for MCI and AD patients are associated with respectively intermediate and late progression stages. The results are similar for both scenarios, although by adding the follow-up information we largely reduce the uncertainty in the prediction of the individual's pathological stage (subfigure 1a vs 2a). NL: normal individuals, MCI: mild cognitive impairment, AD: Alzheimer's patients.

Figure 4

 4 Figure 4 reports the classification results based on the baseline information only, and on increasing thresholds of the progression time course. Although the model is not optimized to explicitly classify the clinical groups, the simple thresholding based on the model predictions generally shows high face validity with respect to the clinical diagnosis. For all the considered scenarios, the highest accuracy is reached in a time window around the point t = 0, while the area under the ROC curve is .99, .88 and .83 for NL vs AD, MCI converters vs MCI stable, and NL converters vs NL stable, respectively.

Figure 4 :

 4 Figure 4: Predictive accuracy of the model when considering the joint set of available biomarkers measurements. The vertical bar indicates the reference threshold value of t = 1.65, corresponding to the 15th percentile of the time distribution of the training AD group. MCI: individuals with mild cognitive impairment, AD: Alzheimer's patients.

  of t = 1.65, corresponding to the 15th percentile of the distribution of the expected time-shift in the training AD group, as well as the area under the receiving operating characteristic curve (AUC). Although the highest discriminative results are consistently obtained when the biomarkers are used jointly, we note that the neuropsychological tests generally lead to the best predictive performance in identifying AD patients with respect to healthy individuals, followed by brain hypo-metabolism (FDG-PET), and temporal atrophy (Entorhinal and Hippocampal volume). This is related to the lower uncertainty of the modelled progressions, which leads to a more accurate identification of the individual staging along the pathological trajectory. The scenario sensibly changes in the other comparison scenarios (MCI conv vs stable and NL conv vs stable), where the sensitivity of the neuropsychological scores shows an important drop, while the other biomarkers (especially hippocampal and entorhinal volumes) provide comparable or even better discriminative performances. These figures were similar when considering the single biomarkers within the longitudinal setting, where the neuropsychological tests still outperformed the other biomarkers in discriminating the clinical groups (Supplementary Figure A.11).

Figure 5 :

 5 Figure 5: Posterior prediction on testing data by using a single biomarker and the baseline information only.

Figure 6 :

 6 Figure 6: Chronological age (y-axis) vs model staging (x-axis). The estimated time-shift is decorrelated from the chronological age in both training and testing data (p>0.05).

  study explores the use of DPM for probabilistic diagnosis and uncertainty quantification in an hypothetical clinical scenario. The proposed approach is based on the reformulation of DPM through a novel probabilistic approach aimed at leveraging on the longitudinal modeling of disease progression for prediction and quantification of the diagnostic uncertainty in neurodegeneration, by optimally combining the information provided by the several biomarkers into a biologically plausible and intelligible score quantified by the time shift. This work thus extends the previous contributions by proposing DPM as a probabilistic tool for diagnostic purposes, which can be used to quantify staging and predictive uncertainty of de-novo individuals in clinical trials. The disease progression model itself thus can be seen as a novel biomarker of pathological

Figure A. 8 :

 8 Figure A.8: Predictive uncertainty with respect to individual's biomarkers variability

Figure A. 10 :

 10 Figure A.10: Joint temporal progression of the biomarkers (top) and derivative of the modelled average trajectory (bottom). ADAS13 and FAQ are characterised by very similar longitudinal profiles, and show the largest changes in the latest stages of the pathology (peak of the derivative at t>0). On the contrary, the change in hippocampal volume is more strongly associated with the earlier stages of the pathology. AV45 and ventricles volumes are the least informative and are associated with the lowest changes (lowest derivative values).

Table 1 :

 1

	Group	NL	NL converted MCI stable MCI converted	AD
			Training data
	N	67	5	0	53	75
	Age	73 (6)	81.4 (5.2)	/	72 (7.7)	73 (8.5)
	Sex (% females)	61	0	/	43	45
	Education (yrs)	16.2 (2)	17.2 (3)	/	15.8 (2.6)	16 (2.4)
	ADAS13	8.8 (4.5)	13.8 (2.4)	/	22.6 (6.7)	31.3 (8.5)
	FAQ	0.2 (0.6)	0.4 (0.5)	/	5.2 (4.5)	12.9 (7)
	RAVLT learning	5.6 (2.6)	2.2 (1.9)	/	3.2 (2.5)	1.8 (1.7)
	Entorhinal (cm 3 )	3.9 (0.6)	3.7 (0.5)	/	3.2 (0.7)	2.9 (0.6)
	Hippocampus (cm 3 )	7.5 (0.9)	6.7 (0.7)	/	6.2 (0.9)	6 (9.3)
	Ventricles (cm 3 )	36 (20)	57 (26)	/	42 (21)	47 (22)
	Whole brain (cm 3 )	1057 (105)	1106 (116)	/	1040 (107)	1013 (113)
	FDG	6.6 (0.5)	6.1 (0.65)	/	5.7 (0.6)	5.2 (0.64)
	AV45	1.2 (0.2)	1.3 (0.09)	/	1.4 (0.2)	1.4 (0.2)
			Testing data
	N	74	17	243	106	145
	Age	75.3 (5.9)	76.5 (4)	73.3 (7)	73.6 (7.3)	75 (7.9)
	Sex (% females)	55	41	39	40	39
	Education (yrs)	16 (2.9)	16.2 (2.6)	16 (2.8)	16 (3)	15.3 (3.1)
	ADAS13	9.8 (4)	11.7 (3.4)	15.7 (6)	21 (6.1)	29.4 (8.2)
	FAQ	0.5 (1.3)	0.6 (1.6)	2.7 (3.5)	5.1 (4.7)	12.9 (6.8)
	RAVLT learning	5.6 (2.2)	5.6 (2.7)	4.3 (2.5)	2.8 (2.2)	1.8 (1.9)
	Entorhinal (cm 3 )	3.8 (0.4)	3.6 (0.7)	3.6 (0.7)	3.1 (0.7)	2.7 (0.7)
	Hippocampus (cm 3 )	7.2 (0.7)	7.2 (0.8)	6.9 (1)	6 (0.8)	5.7 (0.1)
	Ventricles (cm 3 )	33 (15)	44 (21)	39 (23)	41 (23)	49 (24)
	Whole brain (cm 3 )	1019 (102)	1055 (93)	1056 (100)	992 (110)	972 (124)
	FDG	6.5 (0.62)	6.4 (0.7)	6.3 (0.7)	5.9 (0.6)	5.4 (0.7)
	AV45	1.21 (0.19)	1.4 (0.2)	1.3 (0.19)	1.4 (0.2)	1.4 (0.2)
	The model was trained on a group of 200 randomly selected individuals in-
	cluding healthy volunteers, mild cognitive impairment subjects converted to AD
	(MCI conv), and AD patients having at least one measurement for each of the
	1 adni.bitbucket.io/adnimerge.html

Baseline sociodemographic and clinical information for training and testing study cohort. NL: normal individuals, MCI: mild cognitive impairment, AD: Alzheimer's patients. ADAS13: Alzheimer's Disease Assessment Scale-cognitive subscale, 13 items; FAQ: Functional Assessment Questionnaire; RAVLT learning: Rey Auditory Verbal Learning Test, learning item; FDG: (18)F-fluorodeoxyglucose positron emission tomography (PET) imaging; AV45: (18)F-florbetapir Amyvid PET imaging.

ues of 192 pg/ml. The information was extracted from the ADNIMERGE 1 R package

[START_REF] Core | R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing[END_REF]

(MEDIAN field of the upennbiomk_master table). This preliminary selection is aimed to validate the model on a clinical population likely to represent the whole disease time-span. including a subgroup of MCI non converted to AD during the observational time (MCI stable). The image-derived measures used in the study (volumetric MRI and average uptake values for AV45-and FDG-PET) were the scalar estimates reported in the ADNIMERGE package (adnimerge table). The volumetric measures were scaled by the individual total intracranial volume, and all the biomarkers measurements were converted into quantile scores (0 to 1 for normal to abnormal values), with respect to the biomarkers distribution of the training set. This latter modeling precaution is aimed to avoid spurious correlation between training and testing data due to the combined normalization of the values. The modeling results were evaluated with respect to the baseline diagnostic information reported in the ADNI database, assessed according to the WMS and NINCDS/ADRDA AD criteria [23]. Conversion to MCI or AD was established according to the last follow-up information. Moreover, the MCI group was composed by 138 individuals with baseline diagnosis of early MCI, assessed through the Wechsler Memory Scale Logical Memory II. Among these subjects, 14 of them were in the training group (26% of the total MCI training set size), 2.3. Longitudinal modelling of Alzheimer's disease progression 2.3.1. Model training The model was applied in order to estimate the temporal biomarker evolution and the disease stage associated with each individual in the training set. The plausibility of the model was assessed by group-wise comparison of the predicted time-shift, and by correlation with respect to the time to AD diagnosis for the MCI individuals subsequently converted to AD.

For sake of comparison we also correlated the progression modelled with our approach with respect to the one estimated with the method proposed in

[START_REF] Donohue | Estimating longterm multivariate progression from short-term data[END_REF]

. The method was applied to the training data by using the standard parameters defined in the R package grace 3 (see Supplementary Material AppendixA.2.2 for further details).

3 https://mdonohue.bitbucket.io/grace/

Table 3 :

 3 Classification results by using the reference time threshold of t = 1.65, corresponding to the 15th percentile of the training AD time distribution .

						Biomarker
		all	Hippo Ventr WholeBr Entor FDG AV45 RAVLT FAQ ADAS13
					NL vs AD (145 vs 74)
	Accuracy	.89	.81	.62	.76	.83	.80	.63	.82	.88	.83
	Sensitivity	.83	.84	.52	.9	.82	.74	.82	.76	.84	.75
	Specificity	.98	.76	.80	.46	.83	.89	.46	.94	.97	.98
	AUC	.99	.87	.69	.7	.89	.87	.73	.91	.98	.98
				MCI conv vs MCI stable (106 vs 243)
	Accuracy	.82	.67	.62	.69	.7	.71	.69	.67	.79	.79
	Sensitivity	.65	.85	.5	.89	.74	.65	.37	.56	.63	.54
	Specificity	.90	.59	.68	.60	.68	.73	.75	.71	.86	.9
	AUC	.88	.79	.61	.78	.76	.74	.61	.66	.81	.82
				NL conv vs NL stable (17 vs 74)
	Accuracy	.83	.70	.71	.54	.77	.76	.73	.83	.82	.83
	Sensitivity	.18	.47	.41	.82	.52	.29	.27	.35	.17	.17
	Specificity	.98	.77	.80	.47	.83	.89	.86	.94	.97	.98
	AUC	.83	.71	.65	.63	.74	.65	.65	.7	.63	.68

  Data collection and sharing for this project 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.; 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; 625 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 organisation is the Northern California Institute for Research and Education, 630 and the study is coordinated by the Alzheimer's Disease Cooperative Study at the University of California, San Diego. ADNI data are disseminated by the Laboratory for NeuroImaging at the University of Southern California.

	7. Acknowledgments
	EPSRC grants EP/J020990/01 and EP/M020533/1 support DCA and SO's
	work on this topic. DCA and SO also received support from the European
	Union's Horizon 2020 research and innovation programme under grant agree-
	ment No 666992 (EuroPOND) for this work. MF gratefully acknowledges sup-
	port from the AXA Research Fund.

p f , ḟ , ν j |t, t ∼ GP (f joint |0, Σ joint ) ,

ADAS13: Alzheimer's Disease Assessment Scale-cognitive subscale, 13 items; FAQ: Functional Assessment Questionnaire; RAVLT learning: Rey Auditory Verbal Learning Test, learning item; FDG: (18)F-fluorodeoxyglucose positron emission tomography (PET) imaging; AV45: (18)F-florbetapir Amyvid PET imaging.

AppendixA. Supplementary Information

AppendixA.1. Joint Model: marginal likelihood and inference Given the sets of individual biomarker measurements y = {(y j (t i )) k j i=1 } N j=1 , and of D control derivatives m = {m b k (t l )} D l=1 at points t = {t l } D l=1 for the progression of each biomarker b k , the random effect GP model posterior is:

where ν = {ν j } N j=1 . Thanks to the linearity of GPs under derivation, we have that Cov f (t), ḟ (t ) = dCov(f (t),f (t )) dt , and that the joint distribution p f , ḟ |t, t is again a GP sical approaches to Gaussian process with binary activation functions [START_REF] Nickisch | Approximations for binary Gaussian process classification[END_REF]. Following [START_REF] Riihimäki | Gaussian processes with monotonicity information[END_REF], we compute an approximated posterior distribution q f , ḟ , ν j |y j , m by replacing the derivative likelihood terms with local un-normalized Gaussian approximations:

where

Zkl ,

and Σ is a diagonal matrix with elements σ2 kl . It follows that

785

the marginal posterior has a Gaussian form, q f , ḟ , ν j |y j , m ∼ N (µ, Σ), with µ = Σ Σ-1 μjoint , and Σ = (Σ -1 joint + Σ-1 joint ) -1 , where

AppendixA.1.2. Estimating the EP parameters.

The EP update of the local Gaussian approximation parameters is classically done by iterative moment matching with respect to the product between the cav-790 ity distributions q -k l ḟb k (t l ) and the target likelihood term Φ 1 λ ḟb k (t l ) . In the GP case the cavity distribution has a straightforward Gaussian form:

AppendixA.2. Model benchmarking on synthetic multivariate progressions

We benchmarked the model with respect to synthetic multivariate biomarker progressions. We generated random multivariate sigmoid functions for N b bio-

τ ∈ [0, 15] and α k ∼ N (0, .06), and we sampled N individual noisy trajectories at time points τ j k :

For each individual we used the same initial sampling time point for every biomarker, while the number of samples per biomarker was allowed to independently vary between 1 and 4. The individual time points were subsequently centered by their mean µ j k to obtain shifted time-points

AppendixA.2.1. Results.

Table A.4 reports summary correlations between time-shift estimation and the ground truth individual sampling time. The correlation values are generally high, and increase with lower noise levels. Interestingly, the increase in number of modelled biomarkers is associated with a better performance in recovering the underlying disease staging. We also observe that larger sample sizes are associated with higher correlation values, especially with increasing noise levels.

We note however an exception for the case σ = 0.1 where, although the overall performance is still high, the correlation slightly decreases with N = 100.

AppendixA.2.2. Model benchmarking with respect to grace

The R package grace (v 1.0) was used to estimate the multivariate biomarker progression curves from the training set used in this study, by using default parameters and syntax: g r a c e . s i m u l a t i o n . f i t s <-with ( output_table , . . . s e l e c t=c ( "RID" , " Month " , "DX" , " Hippocampus " , " V e n t r i c l e s " , " WholeBrain " , " E n t o r h i n a l " , "FDG" , "AV45" , "RAVLT. l e a r n i n g " , "FAQ" , "ADAS13" , "ICV . b l " ) ) S e t = s u b s e t ( Set , s e l e c t=c ( "RID" , " Month " , "DX" , " Hippocampus " , " V e n t r i c l e s " , " WholeBrain " , " E n t o r h i n a l " , "FDG" , "AV45" , "RAVLT. l e a r n i n g " , "FAQ" , "ADAS13" ) )

#I d e n t i f y i n g i n d i v i d u a l s with a t l e a s t one measurements f o r each biomarker