Rémi Servien
email: remi.servien@inrae.fr

An alternative ranking system for biathlon pursuit races

Keywords: (between 3 and 6): Biathlon, Pursuit race, Ranking system

Biathlon is an Olympic sport combining cross-country skiing with rifle shooting, giving a penalty for each target miss. The biathletes ran different race formats, including the pursuit race. During this race, the biathletes chase the leader with a start time identical to the result of the sprint race previously achieved. So, pursuit involves different skills (such as tactics or management of emotional pressure) that are not present during races with an interval-start procedure like sprint. Nevertheless, final pursuit rankings are strongly correlated to sprint ones, which prevents a spectacular comeback after a disappointing sprint race. We present here an alternative pursuit ranking system that is nearly decorrelated to sprint rankings. This simple ranking system is based on comparisons with previous pursuit results. The current and the alternative rankings were then compared on different pursuit rankings, using a database of 148 results from men pursuit world cups. The alternative ranking was shown to strongly modify a single pursuit ranking but these modifications were smoothed on a whole world cup season.

Advantages and limitations of the alternative ranking system are discussed, paving the way to a fairer modification of the current pursuit ranking to increase surprise and suspense in biathlon pursuit races.

Introduction

Biathlon is an Olympic sport combining 3 to 5 laps of cross-country skiing with rifle shooting.

Between each skiing lap, biathletes complete a shooting session in which they attempt to hit five targets placed at a distance of 50 m, alternating between the prone and standing shooting positions between laps. A penalty (time or skiing distance) is given for each target miss. The biathlete with the shortest overall time wins the race (International Biathlon Union, 2020).

Several different biathlon events exist, in which the individual distance was included as an official Olympic event in 1960, followed by sprint (1980), pursuit (2002), and mass start (2006) (International Biathlon Union, 2020). More precisely, in pursuit races, the 60 best biathletes from the sprint race chase the leader with a start time identical to the result of the sprint race achieved a few days before (i.e. if the second biathlete arrives 12s after the winner of the sprint race, he will start 12s after the first for the pursuit race and so on). So, two of the four current individual Olympic biathlon races involve direct confrontation (mass-start and pursuit), where biathletes are fighting against each other, not versus time. In these events, tactics play a major role and the final ranking is often decided during the last shooting and/or the final skiing sprint.

Furthermore, tight duels during the shootings and the subsequently increased emotional pressure [START_REF] Vickers | Performing under pressure: the effects of physiological arousal, cognitive anxiety, and gaze control in biathlon[END_REF] influence shooting times and accuracies differently than for races with an interval-start procedure. During pursuit or mass-start races, drafting behind other skiers, locating oneself optimally in the crowd also helps maximize the utilization of individual skills (Laaksonen et al., 2018b). Finally, in pursuit races, the skiing speed exerts less impact on the overall performance than in sprint, since the pursuit event involves four bouts of shooting with shorter skiing loops between (Laaksonen et al., 2018b). The pursuit race is therefore expected to reward different skills than sprint or individual races.

Nevertheless, and despite its increasing public audience (EBU, 2019), the biathlon has been sparsely studied, as highlighted by the fact that a search in PubMed with "biathlon" as a keyword currently results in 107 hits, whereas a similar search with "cross-country skiing" (resp. "sport shooting") yields almost 8 (resp. 7) times as many hits. Among these references, the impact of different parameters on shooting accuracy [START_REF] Gallicchio | The influence of physical exercise on the relation between the phase of cardiac cycle and shooting accuracy in biathlon[END_REF][START_REF] Josefsson | Mindfulness and shooting performance in biathlon. A prospective study[END_REF] or the influence of the different biathlon phases on sprint or individual results (Laaksonen et al., 2018a;[START_REF] Luchsinger | The Contribution From Cross-Country Skiing and Shooting Variables on Performance-Level and Sex Differences in Biathlon World Cup Individual Races[END_REF] were extensively examined. Despite their specific aspects, the pursuit and the mass-start races are almost unexplored. Recently, [START_REF] Luchsinger | Contribution from crosscountry skiing, start time and shooting components to the overall and isolated biathlon pursuit race performance[END_REF] investigated the contribution from cross-country skiing, sprint race performance, and shooting components to the pursuit race performance. Sprint race performance was found to be the most influential factor, explaining more than 50% of the final pursuit performance. This result and the fact that the sprint races are the most numerous events (approximatively 40% of the events, 30% being pursuits, 20% mass-starts, and 10% individual races) during a world cup season involve that more than 55% of the final overall world cup results are due to sprint races, which seems very high. Also, the specific skills needed for the pursuit races (tactics, management of emotional pressure …) are not rewarded by the current pursuit ranking, mostly hidden by the importance of the sprint performances. By consequence, decreasing the correlation between sprint and pursuit results is desirable to decrease the influence of the sprint results in the world cup ones, to reward the pursuit specific skills, and to increase suspense and surprise by allowing biathletes that poorly performed during the sprint race to obtain a good pursuit result. An alternative pursuit ranking which minors the influence of the sprint results would therefore be of high interest for biathletes and organizers of international biathlon events. Different rankings than official ones have been recently developed in numerous sports, for example for football teams (Gásquez and Rovuela, 2016), for football players [START_REF] Wolf | A football player rating system[END_REF], for tennis [START_REF] Kovalchik | Extension of the Elo rating system to margin of victory[END_REF], for basketball [START_REF] Barrow | Ranking rankings: an empirical comparison of the predictive power of sports ranking methods[END_REF] … We refer the interested reader to the review of [START_REF] Wunderlich | Forecasting the outcomes of sports events: A review[END_REF] for more details. But, to our knowledge, none of the previous works could be easily adapted to our specific biathlon pursuit problem. Therefore, the current paper aimed to propose an alternative, simple, and fairer ranking for the biathlon pursuit and to investigate its impact on pursuit races and world cup pursuit rankings.

Materials and methods

Data collection

The final results of all sprint and pursuit races are publicly available on the datacenter webpage

Alternative pursuit ranking

As explained previously, during pursuit races, biathletes are racing each other in real-time for a better rank and not racing against time. Therefore, we chose to work using final ranks, not final times. All the pursuit results were gathered to compute final pursuit ranks according to the starting pursuit rank. This global information is given in Figure 1 for all sprint ranking positions and the distributions of the final pursuit ranks are provided in Figure 2 for some starting ranks.

[***Figure 1 We propose an alternative approach to define a fairer pursuit final ranking that will decrease this correlation. For each starting biathlete at a pursuit race 𝑘, a quantity 𝑞 𝑘𝑖 is calculated according to the position of his final results 𝑓 𝑘𝑖 in the final result distribution of all previous starters with the same rank 𝑖. Some of these distributions are plotted in Figure 2. This quantity is given by the following formula

𝑞 𝑘𝑖 = 1 148 ∑ 𝟙 (𝑓 𝑗𝑖 <𝑓 𝑘𝑖 ) 148 𝑗=1
where 𝑓 𝑗𝑖 denotes the final pursuit rank of the biathlete with the starting pursuit rank 𝑖 at the race 𝑗 and 𝟙 (𝑓 𝑗𝑖 <𝑓 𝑘𝑖 ) is the usual indicator function that is equal to one when 𝑓 𝑗𝑖 < 𝑓 𝑘𝑖 , and zero otherwise. Each quantity 𝑞 𝑘𝑖 can be viewed as a quantile of the distribution of the (𝑓 𝑗𝑖 ) 𝑗=1,…,148 .

Then, the quantities (𝑞 𝑘𝑖 ) 𝑖=1,…,60 are ordered, which provides the final ranking of the pursuit race 𝑘. To break the ties, the best current pursuit rank is selected. This rule ensures that the first finisher of the pursuit race will be ranked first at our alternative final pursuit ranking, the best pursuit biathlete will still win. Estimating a quantity (here a performance during a pursuit race) using comparisons to existing distributions (such as the previous pursuit results) is a common approach in non-parametric statistics [START_REF] Wasserman | All of nonparametric statistics[END_REF]. It prevents making false assumptions on the underlying distributions but requires an important number of datasets. In this application, with 148 datasets of previous biathlon pursuit results, this framework seems adapted. All the different pursuit scenarios (including extreme ones) will by consequence be taken into account according to their previous occurrences. Note that, if we consider each of the 148 pursuit races as a ballot, a connection can be made with the well-known Borda count ranking [START_REF] Borda | Mémoire sur les élections au scrutin[END_REF] for the use of all the pursuit races with strong differences for the other part of the alternative ranking computations. We can also notice that there exists a continuum of possible ranking functions between the non-parametric alternative one defined above and the current one. So, if needed, the alternative method could be modified to bring fewer modifications to the pursuit rankings.

But this will complicate the formula provided above and, thus, is beyond the scope of this paper.

Data analyses

All the data analyses were performed using the R freeware, version 3.5.1 (R Core Team, 2019).

The correlations were calculated using Spearman's rank correlation coefficient. For the world cup rankings, we remind that only the first forty biathletes of each race score points, according to the current rules of IBU (International Biathlon Union, 2020). Note that all the races are taken into account to compute all the alternative rankings, even those that took place in the future. This will ensure a good estimation of the different distributions of the final rankings according to the starting ranks.

Results

Study of a specific pursuit race

We first choose to study a specific pursuit race to illustrate the modifications induced by our alternative ranking. We choose the pursuit race that took place at Annecy -Le Grand Bornand (21 December 2019). The results are given in Table 1.

[***Table 1 near here***]

The correlation between the starting rank and the current pursuit rank (resp. the alternative pursuit rank) is equal to 0.82 (resp. 0.20) which highlights the decreased influence of the sprint results on the alternative ranking. If we look at the main modifications we could see that T.

Boe, B. Doll, E. Bjoentegaard, or J. Dale are losing more than 15 ranks with the alternative pursuit ranking. This is due to the fact that they had lost ranks during the pursuit and, therefore, their current good pursuit ranks are mainly due to their good performances in the sprint race.

So, it seems logical that they lose ranks with the alternative ranking. On the contrary, E.

Jacquelin, S. Schempp, and T. Bormolini performed very well during the pursuit race (resp. 14, 22, 28 ranks won during the pursuit race) and deserve their better pursuit rank using the alternative ranking. For example, T. Bormolini will be ranked 6 th with the alternative pursuit ranking whereas it never happened on all the past 148 pursuit races with the current ranking system for a biathlete with the 60 th starting rank, as it could be seen on the last plot of Figure 2.

The computer code used to obtain the results of Table 1 is provided as Supplementary Material with the corresponding dataset. This code could be reused with any pursuit result to compute the alternative rankings in less than a second on an ordinary laptop.

Study of the 2019/2020 pursuit world cup ranking

As explained above, the alternative pursuit ranking can lead to major modifications on a specific pursuit race. Then, we chose to study the 2019/2020 pursuit world cup to analyze the modifications at the scale of a whole season. The first ten biathletes using the two pursuit rankings on all the 2019/2020 races are given in Table 2.

[***Table 2 near here***] First, we could see that there is only a small modification on the podium, J. Boe who was 4 th with the current ranking is now 3 rd whereas Q. Fillon-Maillet who was 3 rd is now 5 th . There is no modification for the first two ranks and eight biathletes are in the two top 10. The strong modifications of the rankings of each pursuit race (as seen in the previous subsection) lead to non-negligible but with less impact world cup ranking modifications. Nevertheless, we can note some important individual modifications for example for E. Garanichev (26 th with the current ranking and 6 th with the alternative) or M. Krcmar (resp. 21 st and 10 th ) who benefit from the alternative ranking unlike T. Boe (resp. 6 th and 11 th ) or S. Desthieux (resp. 7 th and 18 th ).

The number of points with the alternative ranking seems lower than the current one. Indeed, that is an important property of the alternative ranking: the points are awarded to most biathletes as they are less linked to the sprint results (71 biathletes with the current ranking and 81 with the alternative one). But there is a strong correlation of 0.83 between the number of points of each biathlete with the current or the alternative ranking which could explain the relatively small modifications between the two rankings, as mentioned above. strongly modified, for example, a biathlete who was 3 rd with the current ranking is 16 th with the alternative one highlighting the importance of his sprint results in his good current pursuit ranking.

Discussion

Advantages of the alternative ranking

First, the main advantage of this alternative ranking is obviously that the correlation with the starting rank is very low. Therefore, even the 60 th ranked at the end of the sprint had a chance to be on the podium which is not the case with the current ranking. In other words, a biathlete with high pursuit skills who had a bad sprint result has an increased chance of obtaining a good pursuit result with the alternative ranking. But we have to keep in mind that the best pursuit biathlete is still the one that cut the finish line first. It will result in more surprising and contested pursuit races, at each stage of the races, which is desirable for gaining audience [START_REF] Bizzozero | The importance of suspense and surprise in entertainment demand: Evidence from Wimbledon[END_REF]. More generally, it will also decrease the importance of the sprint races on the overall world cup rankings.

Second, even if the alternative ranking deeply modifies each pursuit ranking, each season pursuit world cup ranking is less modified than each single pursuit race. It sounds natural as, even if tactics and head-to-head are of major importance in pursuit races, it remains biathlon with cross-country skiing and shootings. So, the best biathletes are globally the same, the alternative pursuit ranking allows to define the pursuit as a whole discipline with real specialists, not just as a relatively small perturbation of the sprint ranking (as proven in [START_REF] Luchsinger | Contribution from crosscountry skiing, start time and shooting components to the overall and isolated biathlon pursuit race performance[END_REF].

Limitations of the alternative ranking

The first criticism that could be made to the alternative ranking is that it is more complicated than the current one. Nowadays, when you cross the finish line of the pursuit race in 3 rd place, you are ranked 3 rd , whereas with the alternative ranking you need to wait for all the biathletes to finish the race. Even if the alternative ranking is calculated in less than one second at the end of the race, it could be seen as a limiting factor. Nevertheless, this argument needs to be mitigated. First, the winner of the pursuit race is necessarily the winner of the alternative pursuit ranking and is therefore known immediately as he crosses the finish line. Then, for the sprint or individual biathlon races or for other sports such as the decathlon (where you need to refer to a complex points system to see how many points you score, see Cox et al. (2002) for further details) the final ranks are unknown until the last athlete crosses the finish line. This could induce important cliffhangers when biathletes are waiting in the finish area to wait and see if they are or not on the podium. Finally, at each split time, a ranking based on the alternative pursuit ranking could be quickly calculated to inform the biathletes of their rankings.

Another limitation is that, when you have several biathletes that did not start or did not finish the pursuit race despite their presence on the first 60 biathletes of the sprint, it artificially increases the alternative rankings of biathletes that are at the end of the ranking. That could induce unmerited good alternative pursuit rankings for biathletes that have not performed well during the pursuit race but who improved their final rankings thanks to those who gave up. It could be solved by integrating the number of finishers of each pursuit race in the formula to calculate the quantity 𝑞 𝑘𝑖 for example by dividing 𝑓 𝑗𝑖 and 𝑓 𝑘𝑖 by the number of finishers of each pursuit race. But, to keep a very simple formula and as it is uncommon and does not impact the more important highest ranks, it was not taken into account in this paper.

Conclusion

The alternative pursuit ranking presented in this paper is less correlated to the starting ranking than the current one. Some minor limitations remain but, if considered as important, could be easily corrected. This paper paves the way to a fairer modification of the current pursuit ranking that will also increase surprise and suspense in the pursuit races.

  of the IBU: https://biathlonresults.com/. The results were collected on the 15 th December 2020 starting from the 2001/2002 season. All the results taken into account for the men's pursuit world cup were gathered, including world championships and Olympic games before 2014. It provides us 148 different pursuit results. All these results provide us an important database of pursuit results in different conditions that are considered to give an appropriate representation of the different possible pursuit scenarios.

  near here***] [***Figure 2 near here***] These figures emphasize the results of Luchsinger et al. (2020), highlighting the importance of the starting pursuit rank in the final pursuit result.

Figures

  Figures

  Figure 2.

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. Barplots of the final pursuit ranks according to six different starting pursuit ranks.

  This formula is somewhat natural and explainable: indeed, if 𝑞 𝑘𝑖 is equal to zero, it means that, during the previous 148 pursuit races, no biathlete with starting rank 𝑖 achieved a final better (i.e. smaller) rank 𝑓 𝑘𝑖 and so, for all 𝑗, 𝑓 𝑗𝑖 ≥ 𝑓 𝑘𝑖 . So, this biathlete deserves a good alternative final pursuit rank, whatever his starting rank. On the contrary, if 𝑞 𝑘𝑖 is equal to one (i.e. for all 𝑗, 𝑓 𝑗𝑖 < 𝑓 𝑘𝑖 ), then 𝑓 𝑘𝑖 is the worst final pursuit rank achieved by any of the 148 biathletes with this starting rank 𝑖 and it must lead to a poor alternative final pursuit rank.

Tables 327 Table 1 .

 3271 Sprint ranks, current and alternative final pursuit rank for the pursuit race in 328Annecy-Le Grand Bornand in 2019. The gain is the difference between the alternative and the

	329
	330	current pursuit rank.
		Current pursuit rank Name	Sprint rank	Alternative pursuit rank	Gain
		1	BOE Johannes Thingnes	4	1	0
		2	FILLON MAILLET Quentin	3	5	-3
		3	CHRISTIANSEN Vetle	13	2	1
		4	BOE Tarjei	2	25	-21
		5	Sjaastad DOLL Benedikt	1	42	-37
		6	JACQUELIN Emilien	20	4	2
		7	FOURCADE Martin	12	12	-5
		8	BJOENTEGAARD Erlend	5	24	-16
		9	PEIFFER Arnd	21	10	-1
			SCHEMPP Simon	32	3	7
			HORN Philipp	25	7	4
			DALE Johannes	6	38	-26
			LOGINOV Alexander	11	23	-10
			KRCMAR Michal	17	15	-1
			DESTHIEUX Simon	8	35	-20
			EBERHARD Julian	10	29	-13
			WINDISCH Dominik	7	48	-31
			ILIEV Vladimir	22	14	4
			PONSILUOMA Martin	15	30	-11
			PIDRUCHNYI Dmytro	18	27	-7
			LAPSHIN Timofei	19	26	-5
			CLAUDE Florent	23	20	2
			KUEHN Johannes	14	41	-18
			HOFER Lukas	9	49	-25
			TRSAN Rok	47	9	16
			EDER Simon	26	21	5
			PRYMA Artem	34	17	10
			LABASTAU Mikita	46	11	17
			DUDCHENKO Anton	31	19	10
			ELISEEV Matvey	28	33	-3
			PORSHNEV Nikita	24	37	-6
			BORMOLINI Thomas	60	6	26
			RASTORGUJEVS Andrejs	36	22	11
			CLAUDE Fabien	35	31	3
			FEMLING Peppe	56	8	27
			GARANICHEV Evgeniy	30	40	-4
			SEPPALA Tero	33	34	3

Acknowledgments

The author thanks Benoît Rey-Robert for fruitful comments and English corrections and the anonymous reviewer for his/her fruitful comments that helped us to improve the paper.

Note that these small modifications could have a major impact on the overall world cup ranking. Indeed, J. Boe won the overall world cup with 2 points ahead of M. Fourcade. With the alternative ranking, M. Fourcade would have won the overall world cup with the same margin.

Obviously, this is science fiction as the application of the alternative ranking would probably modify the pursuit races. Nevertheless, it could highlight the importance of the sprint results in the overall world cup ranking (J. Boe won 4 of them this season) and the potential impact of the alternative ranking on the overall world cup rankings, mainly when there are few points of difference.

Study of the last ten pursuit world cup seasons

We then studied the pursuit world cup seasons of the ten last years to analyze if the previous remarks could be extended. First, on the pursuit races, the correlations between the starting ranks and the pursuit ranks decreased as seen in the first results subsection: the correlation mean is equal to 0.74 with the current ranking and to 0.06 with the alternative one. Then, we analyzed the last ten pursuit world cup rankings. For all the rankings, there are more biathletes with points with the alternative ranking than with the current one with a mean increase of 11 biathletes. The mean of the difference of points between the first rank and the ranks from 2 to 10 are also all smaller for the alternative pursuit ranking. This would have led to, in most of the cases, closer rankings and more suspense in the last races of the season.

As seen in the previous subsection, the modifications on the podiums of the pursuit world cup rankings are small but not negligible. For 7 seasons we have the same winner, two times the first and the second invert their rankings and for the last one, the 4 th becomes 1 st with the alternative ranking. There are only two identical podiums but, if we compare the name of the first three biathletes, 23 above 30 are shared by the two different rankings. It highlights some important common traits between the two rankings even if some individual rankings could be

Declaration of interest statement

No potential conflict of interest was reported by the author.

Supplementary Material

The computer code and the dataset of the pursuit of Annecy-Le Grand Bornand 2019 are provided as Supplementary Material.


Benoît Sauty

Stanley Durrleman

Progression models for imaging data with Longitudinal Variational Auto Encoders

Keywords: Variational Auto Encoders, Mixed-effects models, Disease progression models, Alzheimer's Disease

Disease progression models are crucial to understanding degenerative diseases. Mixed-effects models have been consistently used to model clinical assessments or biomarkers extracted from medical images, allowing missing data imputation and prediction at any timepoint. However, such progression models have seldom been used for entire medical images. In this work, a Variational Auto Encoder is coupled with a temporal linear mixed-effect model to learn a latent representation of the data such that individual trajectories follow straight lines over time and are characterised by a few interpretable parameters. A Monte Carlo estimator is devised to iteratively optimize the networks and the statistical model. We apply this method on a synthetic data set to illustrate the disentanglement between time dependant changes and inter-subjects variability, as well as the predictive capabilities of the method. We then apply it to 3D MRI and FDG-PET data from the Alzheimer's Disease Neuroimaging Initiative (ADNI) to recover well documented patterns of structural and metabolic alterations of the brain.

Introduction

Estimating progression models from the analysis of time dependent data is a challenging task that helps to uncover latent dynamics. For the study of neurodegenerative diseases, longitudinal databases have been assembled where a set of biomarkers (medical images, cognitive scores and covariates) are gathered for individuals across time. Understanding their temporal evolution is of crucial importance for early diagnosis and drug trials design, especially the imaging biomarkers that can reveal a silent prodromal phase.

In this context, several approaches have been proposed for the progression of scalar measurements such as clinical scores or volumes of brain structures [START_REF] Fonteijn | An event-based model for disease progression and its application in familial Alzheimer's disease and Huntington'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] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF] or series of measurements across brain regions forming a network [START_REF] Bauer | Multiscale modeling for image analysis of brain tumor studies[END_REF][START_REF] Koval | AD Course Map charts Alzheimer's disease progression[END_REF]. These approaches require the prior segmentation and extraction of the measurements from the images. Providing progression models for high dimensional structured data without prior processing is still a challenging task. The difficulty is to provide a low dimensional representation of the data, where each patient's trajectory admits a continuous parametrization over time. It should allow sampling at any time point, be resilient to irregularly spaced instances and disentangle temporal alterations from the inter-patients variability.

Related work

When dealing with high dimensional data, it is often assumed that the data can be encoded into a low dimensional manifold where the distribution of the data is simple. Deep Generative Neural Networks such as Variational Auto Encoders (VAE) [START_REF] Kingma | Auto-encoding variational bayes[END_REF] allow finding such embeddings. Several approaches have explored longitudinal modeling for images within the context of dimensionality reduction.

Recurrent Neural Networks (RNN) provide a straightforward way to extract information from sequential data. Convolutional networks with a recurrent structure have been used for diagnosis prediction using MRI [START_REF] Cui | Rnn-based longitudinal analysis for diagnosis of alzheimer's disease[END_REF] or PET [START_REF] Liu | Classification of alzheimer's disease by combination of convolutional and recurrent neural networks using fdgpet images[END_REF] scans in Alzheimer's Disease (AD). The main caveat of these approaches is that the recurrent structure is highly sensible to the temporal spacing between instances which is troublesome in the context of disease modeling, where visits are often missing of irregularly spaced.

Mixed-effects models provide an explicit description of the progression of each patient, allowing to sample at any timepoint. Through a time reparametrization, all patients are aligned on a common pathological timeline, and individual trajectories are parametrised as small variations (random effects) around a reference trajectory (fixed effects) that can be seen as the average scenario. Now considered a standard tool in longitudinal modeling [START_REF] Koval | AD Course Map charts Alzheimer's disease progression[END_REF][START_REF] Sauty | Riemannian metric learning for progression modeling of longitudinal datasets[END_REF][START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF], mixed-effects models have yet been scarcely used for images within the context of dimension reduction. In [START_REF] Louis | Riemannian geometry learning for disease progression modelling[END_REF], a RNN outputs the parameters of a mixed-effect model that describes patients' trajectories as straight lines in the latent space of a VAE across time.

Self supervised methods have proposed to alleviate the need for labels, in our case the age of the patients at each visit. In [START_REF] Couronné | Longitudinal self-supervision to disentangle inter-patient variability from disease progression[END_REF], the encoder of a VAE learns a latent time variable and a latent spatial variable to disentangle the temporal progression from the patient's intrinsic characteristics. Similarly, in [START_REF] Zhao | Longitudinal self-supervised learning[END_REF], the encoder is penalized with a cosine loss that imposes one direction in the latent space that corresponds to an equivalent of time. Both these methods allow the model to learn a temporal progression that does not rely on the clinical age of the patients, offering potential for unlabeled data, at the cost of interpretability of the abstract timeline and the ability to sample at any given timepoint.

Longitudinal VAEs architectures have been proposed in order to endow the latent space with a temporal structure. Namely, Gaussian Process VAEs (GP-VAE) [START_REF] Casale | Gaussian process prior variational autoencoders[END_REF] introduced a more general prior for the posterior distribution in the latent space, in the form of a Gaussian Process (GP) that depends on the age of the patients [START_REF] Fortuin | Gp-vae: Deep probabilistic time series imputation[END_REF] as well as a series of covariates [START_REF] Ashman | Sparse gaussian process variational autoencoders[END_REF][START_REF] Ramchandran | Longitudinal variational autoencoder[END_REF]. This approach poses challenges as to the choice of parametrization for the Gaussian Process, and does not provide an expected trajectory for each patient.

Diffeomorphic methods provide progression models for images. The main approaches are based off of the geodesic regression framework [START_REF] Banerjee | A nonlinear regression technique for manifold valued data with applications to medical image analysis[END_REF][START_REF] Niethammer | Geodesic regression for image timeseries[END_REF] and allow learning a deformation map that models the effect of time on the images for a given subject. While providing high resolution predictions, these methods show limited predictive abilities further in time when compared to mixed-effects models, that aggregate information from all the subjects at different stages of the disease [START_REF] Bône | Prediction of the progression of subcortical brain structures in alzheimer's disease from baseline[END_REF].

Contributions

In this context, we propose to endow the latent space of a VAE with a linear mixed-effect longitudinal model. While in [START_REF] Louis | Riemannian geometry learning for disease progression modelling[END_REF], the networks predict the random effects from visits grouped by patients, we propose to enrich a regular VAE that maps each individual visits to a latent representation, with an additional longitudinal latent model that describes the progression of said representations over time. A novel Monte Carlo Markov Chain (MCMC) procedure to jointly estimate the VAE and the structure of its representation manifold is proposed. To sum up the contributions, we :

1. use the entire 3D scan without segmentation or parcellation to study relations across brain regions in an unsupervised manner, 2. proceed to dimension reduction using a convolutional VAE with the added constraint that latent representations must comply with the structure of a generative statistical model of the trajectories, 3. provide a progression model that disentangles temporal changes from changes due to inter-patients variability, and allows sampling patients' trajectories at any timepoint, to infer missing data or predict future progression, 4. demonstrate this method on a synthetic data set and on both MRI and PET scans from the Alzheimer's Disease Neuroimaging Initiative (ADNI), recovering known patterns in normal or pathological brain aging.

Methodology

Representation learning with VAEs

Auto Encoders are a standard tool for non-linear dimensionality reduction, comprised of an encoder network q ϕ that maps high dimensional data x P X to z P Z, in a smaller space refered to as the latent space, and a decoder network p θ : z P Z Þ Ñ x P X . VAEs [START_REF] Kingma | Auto-encoding variational bayes[END_REF] offer a more regularized way to populate the latent space. Both encoder and decoder networks output variational distributions q ϕ pz|xq and p θ px|zq, chosen to be multivariate Gaussian distributions. Adding a prior qpzq, usually the unit Gaussian distribution N p0, Iq, on Z allows to derive a tractable Evidence Lower BOund for the log-likelihood ELBO " L recon `βL KL where L recon is the ℓ 2 reconstruction error, L KL is the Kullback-Leibler (KL) divergence between the approximate posterior and the prior on the latent space and β balances reconstruction error and latent space regularity [START_REF] Higgins | beta-vae: Learning basic visual concepts with a constrained variational framework[END_REF].

Longitudinal statistical model

In this section, we propose a temporal latent variables model that encodes disease progression in the low-dimensional space Z. Given a family of observations from N patients tx i,j u 1ďiďN , each observed at ages t i,j for 1 ď j ď n i visits, and their latent representations tz i,j u, we define a statistical generative model with z i,j " p 0 `"e ξi pt i,j ´τi q ‰ v 0 `wi `εi,j where e ξi and τ i , respectively the acceleration factor and the onset age of patient i, allow an affine time warp aligning all patients on a common pathological timeline, and w i P Z is the space shift that encodes inter-subjects variability, such as morphological variations across regions that are independent from the progression. These parameters position the individual trajectory with respect to the typical progression that is estimated at the population level. These three parameters form the random effects of the model ψ r . Vectors w i and v 0 need to be orthogonal in order to uniquely identify temporal and spatial variability. We choose the Gaussian priors for the noise ε i,j " N p0, σ 2 ε q and random effects τ i " N pt 0 , σ 2 τ q, ξ i " N p0, σ 2 ξ q and w i " N p0, Iq. The parameters p 0 P Z, v 0 P Z, t 0 P R are respectively a reference position, velocity and time and describe the average trajectory. Together with the variances σ ε , σ τ , σ ξ , they form the fixed-effects of the model ψ f . We note ψ " pψ r , ψ f q.

Longitudinal VAE

We combine dimension reduction using a regular β-VAE and the aforementioned statistical model to add a temporal structure to the latent space. To do so, we consider a composite loss that accounts for both the VAE loss and the goodnessof-fit of the mixed-effect model:

L " L recon `βL KL `γL align where $ & % L recon " ř i,j ||x i,j ´x i,j || 2 L KL " ř i,j KLpq ϕ pz|x i,j q||N p0, Iqq L align " ř i,j ||z i,j ´ηi ψ pt i,j q|| 2
where z i,j and xi,j are the modes of q ϕ px i,j q and p θ pz i,j q, and η i ψ pt i,j q " p 0 `"e ξi pt i,j ´τi q ‰ v 0 `wi is the expected position of the latent representation according to the longitudinal model and γ balances the penalty for not aligning latent representations with the linear model. Since the loss is invariant to rotation in Z, we set p 0 " 0 and v 0 " p1, 0, ¨¨¨, 0q.

Since L align is a ℓ 2 loss in the latent space, it can be seen as the log-likelihood of a Gaussian prior z i,j " N pη i ψ pt i;j q, Iq in the latent space, which defines an elementary Gaussian Process, and which supports the addition of GP priors in the latent space of VAEs to model longitudinal data [START_REF] Fortuin | Gp-vae: Deep probabilistic time series imputation[END_REF][START_REF] Ramchandran | Longitudinal variational autoencoder[END_REF]. Besides, X can be seen as a Riemannian manifold, the metric of which is given by the pushforward of the Euclidean metric of Z through the decoder, such that trajectories in X are geodesics, in accordance with the Riemannian modeling of longitudinal data [START_REF] Bône | Learning distributions of shape trajectories from longitudinal datasets: a hierarchical model on a manifold of diffeomorphisms[END_REF][START_REF] Gruffaz | Learning riemannian metric for disease progression modeling[END_REF][START_REF] Louis | Riemannian geometry learning for disease progression modelling[END_REF][START_REF] Sauty | Riemannian metric learning for progression modeling of longitudinal datasets[END_REF][START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF]. The metric on X thus allows to recover non linear dynamics, as is often the case for biomarkers. Our approach thus bridges the gap between the deep learning approach to longitudinal data and a natural generalization of well studied disease progression models to images.

Network implementation and estimation. Both the encoder and decoder are chosen to be vanilla convolutional Neural Networks (4 layers of Convolution with stride/BatchNorm/ReLU and transposition for decoder) with a dense layer towards the latent space, as described in Fig. 1. The implementation is in PyTorch and available at https://github.com/bsauty/longitudinal-VAEs.

Algorithm 1: Monte Carlo estimation of the Longitudinal VAE

Input : Longitudinal visits tpxi,j, ti,jqu and hyperparameters β and γ. Output: Estimation of pϕ, θq for the VAE and ψ for the temporal model.

Init

: Initialize pϕ, θq as a regular β-VAE.

Set k " 0 and z 0 " q ϕ pxq; while not converged do Simulation Draw candidates ψ c r " pp.|ψ k f q ; // Sampling with prior pp.|ψ k f q @i, j compute η i ψ k pti,jq ; // Expected latent trajectories

Compute likelihood ratio ω " min ˆ1, The difficulty lies in the joint estimation of pθ, ϕq and ψ, which are codependant since L align depends on η i ψ and ψ depends on the encoded representation z " q ϕ pxq. The longitudinal statistical model is part of a family of geometric models that have been studied in [START_REF] Koval | Statistical learning of spatiotemporal patterns from longitudinal manifold-valued networks[END_REF][START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF]. Given tz i,j u, we can proceed to a Maximum a Posteriori estimation of ψ with the MCMC-SAEM procedure in which the estimation step of an EM algorithm is replaced by a stochastic approximation. See [START_REF] Allassonnière | Construction of bayesian deformable models via a stochastic approximation algorithm: a convergence study[END_REF][START_REF] Kuhn | Coupling a stochastic approximation version of em with an mcmc procedure[END_REF] for details. Given the target trajectories tη i ψ u, the weights from both networks of the VAE is optimized through backpropagation using L with an optimizer with randomized batches. Both estimation schemes are iterative so we designed a Monte Carlo estimator for pϕ, θ, ψq, presented in Algorithm 1, alternating between both schemes.

qpψ c r |z k ,ψ k f q qpψ k r |z k ,ψ k f q if u " Up0, 1q ą ω then ψ k`1
Once calibrated with a training set, we freeze the VAE parameters pϕ, θq and fixed effects ψ f , and learn the individual parameters ψ r , via gradient descent of the likelihood, to personalize the model for new subjects.

Hyperparameters are dim(Z), which should be small enough to allow the mixed-effect model to be interpretable but big enough to reach good reconstruction; β, which should minimize overfitting while not impairing reconstruction quality; and γ, which should also not be too big to avoid loosing contextual information in Z. These parameters were set using grid search. Lastly, the MCMC-SAEM is computationaly inexpensive compared to backpropagation so memory footprint and runtime are similar to training a regular VAE. All training was performed with a Quadro RTX4000 8Go GPU.

Experiments and results

Results on synthetic experiments

We first validated our approach on a synthetic data set of images of silhouettes of dimension 64x64 [START_REF] Couronné | Starmen longitudinal[END_REF] . Over time, the silhouette raises its left arm. Different silhouettes are generated by varying the relative position of the three other limbs, all of them raising their left arm in time. The motion is modulated by varying the time stamp at which the motion starts and the pace of motion. This is done using an affine reparametrization of the time stamp t i,j of the silhouette with Gaussian log-acceleration factor ξ i and onset age τ i . This data set contains N " 1, 000 subjects with n " 10 visits each, sampled at random time-points.

We choose dimpZq " 4 to evaluate the ability of our model to isolate temporal changes (motion of the left arm) from the independent spatial generative factors (the position of the other 3 limbs). Results are displayed in Fig 2 . The 5-fold reconstruction mean squared errors (MSE) (times 10 ´3) for train/test images are 7.88˘.22{7.93˘.29, showing little over-fitting. Prediction error for missing data, when trained on half-pruned data set, is 8.1˘.78 which shows great extrapolation capabilities. A thorough benchmark of six former approaches on this data set was provided in [START_REF] Couronné | Longitudinal self-supervision to disentangle inter-patient variability from disease progression[END_REF] displaying similar MSE to ours. Although a couple of approaches [START_REF] Couronné | Longitudinal self-supervision to disentangle inter-patient variability from disease progression[END_REF][START_REF] Zhao | Longitudinal self-supervised learning[END_REF] also disentangle time from space, ours is the only one to yield the true generative factors, with the direction of progression encoding the motion of the left arm and the 3 spatial directions orthogonal to it encoding legs spreading, legs translation, and right arm position respectively.

Results on 3D MRI and PET scans

We then applied the method to 3D T1w MRI and FDG-PET scans from the public ADNI database (http://adni.loni.usc.edu). For MRI, we selected two cohorts: patients with a confirmed AD diagnosis at one visit at least (N " 783 patients for a total of N tot "3,685 images) and Cognitively Normal (CN) patients at all visits for modeling normal aging (N " 886 and N tot " 3, 205). We considered PET data for AD patients only (N " 570 and N tot " 1, 463). Images are registered using the T1-linear and PET-linear pipelines of the Clinica software [START_REF] Routier | Clinica: an open-source software platform for reproducible clinical neuroscience studies[END_REF] and resampled to 80x96x80 resolution.

We set dimpZq " 16 for both modalities, as it is the smallest dimension that captured the reported dynamics with satisfying resolution. For MRI, error (10 ´3) for train/test reconstruction and imputation on half-pruned data set for the AD model are 14.15 ˘.12{15.33 ˘.23 and 18.65 ˘.76, again showing little over-fit and good prediction abilities. In Fig. 3, the reference trajectory for AD patients reveals the structural alterations that are typical of AD progression. The control trajectory displays alterations more in line with normal aging.

We tested differences in the mean of individual parameters between subgroups using Mann-Whitney U test within a 5-fold cross-validation. AD average onset age occured earlier for women than for men: 72.2˘.4 vs. 73.7˘.6 years, Fig. 2. Synthetic experiment. (a) the average trajectory over time (left to right) on first row, followed by its translation in the directions w1 " p0, 1, 0, 0q and w2 " p0, 0, 1, 0q in the latent space (second and third row). (b) the gradient of the image at p0 " 0 in the 4 directions of the latent space v0, w1, w2, w3. (c) data of a test subject (first row), its reconstructed trajectory (second row) and the ground truth (third row) p ă 3.10 ´7˘5.10 ´8 ). APOE-ε4 mutation carriers experience also earlier onset than non-carriers (71.8˘.2 vs 73.1˘.4, p ă 3.3.10 ´2˘6.10 ´3) and greater pace of progression (.1˘3.10 ´2 vs -.08˘2.10 ´2, p ă 1.4.10 ´4˘6.10 ´3). The normal aging model shows an earlier onset for men (71.2˘.4 vs 73.7˘.6, p ă 3.10 ´10 ˘6.10 ´11 ). These results are in line with the current knowledge in AD progression [START_REF] Gurvich | Sex differences and the influence of sex hormones on cognition through adulthood and the aging process[END_REF][START_REF] Jack | Age, sex, and apoe ε4 effects on memory, brain structure, and β-amyloid across the adult life span[END_REF] and normal aging [START_REF] Coffey | Sex differences in brain aging: a quantitative magnetic resonance imaging study[END_REF]. For PET scans, the 5-fold train/test reconstruction MSE (10 ´2) are 4.71 ˘.32{5.10 ˘.23.

Conclusion

We proposed a generative Variational Autoencoder architecture that maps longitudinal data to a low dimensional Euclidean space, in which a linear spatiotemporal structure is learned to accurately disentangle the effects of time and inter-patient variability, while providing interpretable individual parameters (onset age and acceleration factor). This is the first approach to provide a progression model for 3D MRI or PET scans and it relies on vanilla deep learning architectures that only require the tuning of the loss balance. We showed that it bridges the gap between former approaches to handle longitudinal images, namely GP-VAEs, and Riemannian disease progression models. The method applied to MRI and PET data retrieves known patterns of normal and pathological brain aging but without the need to extract specific biomarkers. It does not only save time but also makes the approach independent of prior choice of biomarkers.

Current work focuses on linking this progression model of brain alterations with cognitive decline, and exploring disease sub-types in the latent space.

ApproximationFig. 1 .

 1 Fig. 1. Images txi,ju are encoded into Z such that the tzi,ju are close to the estimated latent trajectories. Individual trajectories (straight lines) are parametrized with wi,τi and e ξ i as variations around the reference trajectory (orange arrow).

Fig. 3 .

 3 Fig. 3. (a) Sagittal, coronal and axial views of the population trajectory over pathological time (left to right) for the AD cohort. Enlargment of the ventricles and atrophy of the cortex and the hippocampus are visible. Red squares around the hippocampus are positioned at the 5 stages of the Schelten's scale used in the clinics to evaluate AD progression. (b) Coronal view of the estimated normal aging scenario, with matched reparameterized age distribution.Atrophy is also visible but to a smaller extent. As is common in atlasing methods, these images average anatomical details from different subjects to provide a population trajectory, and are thus not as sharp as true images.

Fig. 4 .

 4 Fig. 4. Sagittal, coronal and axial views of the average trajectory for FDG-PET scans, showing decreased level of metabolism across brain regions.

Acknowledgments

This work was funded in part by grant number 678304 (ERC) and 826421 (TVB-Cloud) from H2020 programme, and ANR-10-IAIHU-06 (IHU ICM), ANR-19-P3IA-0001 (PRAIRIE) and ANR-19-JPW2-000 (E-DADS) from ANR.


Louis-Philippe Dormegnie

Adrien Henry

Jean-Hugues Salmon

come

Internal Medicine Flashcard

No more than meets the eye [START_REF] Annamalai | Bilateral choroidal tuberculoma in miliary tuberculosis -report of a case[END_REF]

Main section (case description)

A 33-year old Ivorian woman was admitted for fever and thrombocytopenia. She had been previously treated in France by tocilizumab then anti TNF-α for rheumatoid arthritis since 7 and 2 years respectively. She lost weight three months ago, whereas her last trip to Ivory Coast was three years ago. Physical examination was unremarkable except a hectic fever. Complete Blood count revealed hypochromic microcytic anaemia (6.1 g/dL) and profound thrombocytopenia (11 000/mm3), Serum C reactive protein and ferritin values were 81 mg/L and 636 ng/mL respectively. Haemoglobin level was corrected whereas platelets level was not, once red blood cells and platelet transfusion had been performed.

Bone marrow cytology was normal with presence of megakaryocytes. Platelet associated immunoglobulins were detected by flow cytometry and ELISA (anti GPIIb/IIIa and anti GPIa/IIa). Computerized tomography scan yielded multiple micronodules of lungs, spleen and liver associated with mesenteric lymphadenopathies (Fig. 1A). Funduscopy revealed a single elevated yellowish-white choroidal lesion of the left eye with fluorescein angiographic signs of exudation (Fig. 1B).

Q1

Louis Philippe Dormegnie a, 1 , Adrien Henry a,1 , Jean Hugues Salmon b,2 , Yohan N'Guyen i The corrections made in this section will be reviewed and approved by a journal production editor.

alt-text: Fig. The multiple micronodules in both lungs were consistent with the millet-seed like appearance observed during miliary Tuberculosis, which is a subset of disseminated Tuberculosis. Although less evocative than Choroidal Tubercles (usually multiples, grayish, surrounded by inflammation, with a size inferior to 0.5 disc diameter), the presence here of a choroidal Tuberculoma (always single, yellowish, with exudative detachment and larger size) suggested also strongly

Tuberculosis. It has been described during miliary and disseminated Tuberculosis [START_REF] Annamalai | Bilateral choroidal tuberculoma in miliary tuberculosis -report of a case[END_REF], which were overrepresented (in 5 to 40% of cases) amongst patients treated by anti TNF-α [START_REF] Baronnet | Incidence of tuberculosis in patients with rheumatoid arthritis. A systematic literature review[END_REF]. Systematic funduscopy evidencing Choroidal Tuberculoma allowed here to make (before culture) the diagnosis of Tuberculosis, which is infrequently associated with Immune Thrombocytopenia. Indeed, Immune Thrombocytopenia has only been described in more than 50 cases of all kinds of Tuberculosis so far [START_REF] Weber | Immune thrombocytopenia secondary to tuberculosis: a case and review of literature[END_REF].

Our patient received 12 months of antimycobacterial therapy associated with Corticosteroids because of immune thrombocytopenia. Corticosteroids treatment lasted 4 months because of subsequent Rheumatoid Arthritis flare-ups, which then led to the reintroduction of Etanercept. After the end of antimycobacterial therapy, all micronodules and choroidal tuberculoma had vanished as well as the patient gained all the weight she had lost and did not experience fever or thrombocytopenia recurrence.

Figure Fig. 1 Question: What is the most plausible diagnosis? 2

 12 Figure Fig. 1

  Service de Rhumatologie. Hôpital Maison Blanche. Centre Hospitalier Universitaire Reims. France Service de Médecine Interne, Maladies Infectieuses et Immunologie Clinique. Hôpital Robert Debré. Centre

		c,2,⁎
	yohan.nguyen@wanadoo.fr
	Hospitalier Universitaire Reims. France
	⁎	Corresponding author: Avenue du général Koenig, 51100 Reims.

a Service d'Ophtalmologie. Hôpital Robert Debré. Centre Hospitalier Universitaire Reims. France b c 1 Equally contributed to this work 2 Equally contributed to this work

Acknowledgement

The authors report no source of funding.

Declaration of Competing Interest