M Kuvshinov
email: maxim.kuvshinov@coria.fr

B Lecordier

S Idlahcen

A Cessou

B Renou

S Hochgreb

Experimental study of the effect of droplets on opposed stretched premixed flames

Keywords: laminar flame, counterflow, spray, strained, droplet, premixed, PIV

The effect of fuel droplets on the burning velocity of strained laminar premixed flames was investigated experimentally using optical diagnostics. The twin counterflow burner configuration was used to create premixed acetone vapour flames with the addition of fuel droplets. Particle image velocimetry (PIV) was used to measure 2D velocities for nominal equivalence ratios in the range (0.8-1.21) as a function of strain rate (250-550 s 1 ). Measurements of reference flame speeds upstream of the flames were made for both reference methane/air flames and acetone vapour/droplet flames, and compared to simulated values of the purely gaseous flames. Spontaneous Raman spectroscopy was used to measure temperature and major species molar fractions across reference methane flames. The results were shown to be in good agreement with the simulations for most of the species, except CO and CO 2 , which had too low a signal-to-noise ratio in the product zone. Preliminary measurements of acetone spray/vapour flames showed the capability of the setup to measure Raman signals in the presence of droplets at lowered laser energy. High laser energies led to prompt ignition of the droplets. Potential means of resolving the problem are suggested.

Introduction

Two-phase spray combustion plays an important role in industrial processes, energy production and storage, as well as aeronautical and automotive transportation. Injection of fuel via sprays is actively used in internal combustion engines, gas turbines and power generators. The interaction of liquid and gaseous-phase fuels is complex; models typically involve the solution of Eulerian gas phase equations overlaid with Lagrangian solutions for droplets. Recent studies [START_REF] Vié | Proceedings of the Combustion Institute[END_REF][START_REF] Xie | [END_REF] suggested that the effect of droplets could be taken into account directly by parameterizing the effect of their presence on the burning rate. Previous experimental studies on the role of droplets on the burning rate of vaporized mixtures were performed in closed combustion chambers on spherically expanding flames [3,[START_REF] Nomura | Proceedings of the Combustion Institute[END_REF]. These suggested that heat release rates were maximum for a particular droplet size, which also varied with the equivalence ratio. A number of computational studies have considered the role of droplets on laminar premixed and diffusion flames, which also revealed a similar behaviour [START_REF] Neophytou | [END_REF]6]. The transient nature of spherically propagating flames does not lend itself to detailed measurements to provide evidence for such behaviour.

The present study provides a platform for steady measurements of velocities, temperatures and species on the structure of strained laminar flames with added droplets, using measured droplet sizes and liquid fraction concentrations.

The present study is based on initial work by Mc-Grath et al. [START_REF] Mcgrath | Proceedings of the Combustion Institute 000[END_REF] using the same configuration of acetone/air mixtures to examine the effect of the presence of droplets on flame propagation, compared to an analogous vapour-only mixture. The reference flame speeds at various equivalence ratios were determined as a function of strain rate for droplet-laden spray and vaporized acetone/air laminar premixed flames using particle image velocimetry (PIV). Droplet size and count statistics were obtained using Phase Doppler Particle Anemometry (PDPA), which also yielded the estimated fuel liquid fraction. The study has shown that under lean conditions, incomplete fuel vaporization resulted in leaner vapor-phase equivalence ratios, resulting in lower reference laminar flame speeds. In the case of richer mixtures, large enough droplets were able to cross the flame fronts leading to a fuel surplus in the opposite vapor flame, leading to a reduction in its adiabatic flame temperature and a richer product or reactant zone. One of the limitations of the study was a lack of temperature control of the inlet mixtures. Due to the co-axial vertical arrangement of the burners, the exhaust gases led to transient heating of the top burner, influencing the rate of droplet evaporation and overall flame speed.

The present study builds on that prior work, and incorporates much finer temperature control for both top and bottom flames, allowing for much longer and wellcharacterised detailed measurements.

This study reports on the results of work in progress for these measurements, including (a) PIV measurements for the modified acetone-vapour-droplet system and (b) preliminary Raman measurements of opposed methane-air flames, in preparation for Raman measurements of acetone flames.

Experimental methods

Experimental setup

Fig. 1 shows the experimental setup for the twin counterflow burner configuration. The burner nozzles converge the flows from the diameter of 71 mm to 14 cm over a distance of 68 mm to provide a top-hat velocity profile at the exit. The distance between outlets was set to 24 mm (L/D n = 1.7). Each nozzle is surrounded by a 10 mm wide circular annulus through which a co-flow of nitrogen is delivered to quench the flames and suppress perturbations from the surrounding environment. The top nozzle is surrounded by a stainless steel watercooling jacket to ensure thermal stability. Thermocouples (RS-pro, type K ± 1.5 C) were placed 30 mm from the nozzle exits to measure the fluid temperatures.

The burners were operated in two modes, A and B, producing either pure gas methane flames or acetone spray/vapour flames respectively. In mode A both burners had equivalent fluid delivery systems. Air flow was supplied by a mass flow controller (MFC) (Bronkhorst, F-202AV-M10) into a fluidised bed seeder for particle image velocimetry. Methane was added into the seeded flow and metered by an MFC (Bronkhorst, F-201CV-20K) in order to provide a premixed air/methane mixture at ambient conditions. Nitrogen was used as a sheath flow, supplied by MFCs (Bronkhorst, F-202AV-M10).

The fluid delivery systems in mode B were different top and bottom. The bottom burner was used to produce a vapour flame. Acetone vapour was generated using a controlled evaporation mixing chamber (CEM)(Bronkhorst, W-303A). Liquid acetone was delivered from a pressurised vessel and monitored by a Coriolis mass flow controller (Bronkhorst, M13) into the CEM, which was set to above acetone evaporation temperature of 60 C. An airflow of 5 SLPM was delivered by an MFC (Bronkhorst, F-201CV-20K) into the air inlet of CEM in order to carry the generated vapor. The vapor-containing flow was then mixed with the seeded flow, as described in mode A. The top burner was used to produce an acetone spray flame. An ultrasonic atomiser (Sono-Tek, 8700-120) was used to create a dilute acetone spray, injected at the inlet of the mixing chamber, 275 mm away from the nozzle exit. A Coriolis mass flow controller (Bronkhorst, M13) with a pump (Tuthill, DGS.19) delivered the liquid acetone to the atomiser. A separate airflow was passed through the second fluidized bed seeder by an MFC (Bronkhorst). Particle-seeded flow was injected into the mixing chamber 25 mm downstream from the atomizer outlet, resulting in an acetone and particle droplet-laden air mixture, with droplet Sauter mean diameter range of 65-75 µm and Stokes number range of 0.500-2.477, increasing with increasing strain rate [START_REF] Mcgrath | Proceedings of the Combustion Institute 000[END_REF].

Operating conditions

The mixture compositions were defined by the nominal equivalence ratio, f , based on the ratio of injected mass flow rates of fuel and air. For mode A, f = [0.8, 1.0, 1.2]. For mode B, f = [0.75, 0.84, 0.93, 1.02, 1.12, 1.21, 1.3]. For each nominal equivalence ratio, measurements were made for fluid exit velocities, v ex = [1.52-2.60 m/s] in increments of 0.21 m/s for mode A and v ex = [1.62-2.77 m/s] in increments of 0.32 m/s for mode B. This resulted in the range of strain rates K = 250-550 s 1 away from flame extinction or nozzle rim attachment.

Gas phase velocity measurements

The velocity fields were determined using PIV. The tracer particle illumination was provided by a doublepulsed Photonics Industries DM60-532 DH laser at the rate of 4 kHz, and 0.7 W of power per laser head. The beam was shaped into a 1⇥24 mm 2 sheet. Images were captured by a Phantom v2012 CCD camera, with a 1280 x 800 sensor array, operating in double frame mode at 4 kHz. The camera's field of view was focused by a model K2 Distamax objective. A narrow 532 nm bandpass filter was attached to the objective to reduce background emission. The pixel resolution was 22.7 µm per pixel as calibrated by the target plate. The resultant field of view was set to 17.9⇥28.9 mm 2 . Zirconium powder with an additive of silica oxide (3% by weight) was used to generate gas phase PIV tracer particles (5 µm) with a Stokes number below 0.18. The time interval between the two laser pulses was varied between 54-100 µs, decreasing with increasing flow rates, in order to obtain optimum spatial resolution. A series of 1000 images were collected for every flow rate and equivalence ratio condition for statistical convergence. The images were processed by Davis 10.2.1 using a multipass crosscorrelation algorithm. Velocity vector fields were calculated by 3 time pass 16⇥16 px 2 interrogation window with a 50% overlap, resulting in vector spacing of 180 µm and interrogation region of 360 µm.

The methods developed in [START_REF] Law | Symposium (International) on Combustion[END_REF][START_REF] Zhu | Symposium (International) on Combustion[END_REF] were employed in this study for the determination of reference flame speed and strain rates. First, a stagnation streamline was determined, as the streamline along which the radial velocity is zero. The two minimum axial velocities observed upstream of the acceleration produced by each flame were set as the reference flame speeds S L,ref for the top and bottom flame respectively. The respective strain rates were then determined from the gradient of radial velocities, a, via relation K = 2a at the location of the S L,ref .

Gas phase temperature and mixture composition measurements

Measurements of temperature and gas composition were performed by spontaneous Raman Scattering spectroscopy (SRS). To overcome the weak efficiency of SRS, the molecular excitation was provided by a Nd:YAG (Agilite Continuum) operating at 10 Hz. The laser output was 1.5 J with a top-hat pulse profile and duration of 1 µs. The beam was focused by a f = 1000 mm focal length lens resulting in a beam waist of 200 µm (1/e 2 ) at focus. The isotropically generated Raman and Rayleigh scattered emission signals were collected at 90 by two telescopes. The first telescope collimated the emission via 100 mm diameter, f = 150 mm achromat lens, which was then focused onto a slit by a f = 300 mm, 100 mm lens, resulting in a probe volume length of 3.2 mm. Rayleigh emission was suppressed by a notch filter (NF03-532E, Semrock, OD = 6, at 532nm, FWHM = 17 nm) placed in the collimation part of the secondary telescope, consisting of two 2-inch, f = 200 and f = 250 lenses respectively. A periscope was then used to rotate the emission from the horizontal beam parallel to the vertical entrance slit of the spectrograph. The singleshot SRS spectra were acquired with the IsoPlane SCT 320 spectrograph (Princeton Instruments) which had a grating of 600 grooves/mm, resulting in spectral resolution of 0.67 nm [START_REF] Guichard | Proceedings of the European Combustion Meeting[END_REF].

The detector camera used was a full-frame backilluminated CCD (Pixis 400b, Princeton Instruments) with a 1340⇥400 sensor array and pixel size of 20 µm. A series of 500 images at 3.3 Hz were collected to reach statistical convergence.

The collection was gated by fast electro-optical shutter made-up of Pockels Cells (PCS) to reduce the contribution of flame emission to the background. The PCS were positioned in the collimated part of the second telescope. The gating time for the PCS was set to 1 µs to overlap with the duration of the laser pulse, minimising the camera's exposure to flame emission. An achromatic half-wave plate (AHWP10M-600, Thorlabs) was placed in front of the spectrograph and oriented to maximise its efficiency.

Results

Gas phase velocity measurements

Fig. 2 shows the reference laminar flame speed as a function of strain rate for f = [0.8, 1.0, 1.2], methane gas flames (mode A). Dark grey and white circles represent the single frame and mean values for the top flame, while light grey and white squares represent the single frame and mean values for the bottom flame. The means are grouped by the v ex conditions. Black lines represent the numerical values obtained with Cantera onedimensional model, using GRI 3.0 mechanism. The values for the top flame are predominantly above the data for the bottom flame, since the mixture temperature at the exit of the top burner is 6±0.81 C higher than at the exit of the bottom burner as a consequence of heat transfer to the burner's body from the product gases. The velocity gradients for all of the conditions are similar, indicating that the temperature difference is the underlying cause of the value discrepancies. The simulations were performed for the measured inlet mixture temperature of the bottom burner. There is a good agreement between the model and the experimental results for lean and stoichiometric conditions, with the rich condition showing a sharper gradient increase with strain rate. The use of high-speed PIV allowed the detection of periodic synchronised anti-phase oscillation of the two flame fronts at 130 Hz, which becomes more apparent and higher strain rate conditions. Its effects are clearly observed in the f = 1.0 condition, where starting from 400 s 1 the spread of data from the mean dramatically increases. The same behaviour is observed in f = 1.2. These oscillations are undesirable, and their origin is currently being investigated.

Fig. 3 shows the reference laminar flame speed as a function of strain rate for f = [0.83, 1.02, 1.21], acetone spray/vapour flames (mode B). The data is represented in the same fashion as in Fig. 2. Black lines represent the numerical values obtained with Cantera one dimensional model at 28 C, using mechanism of Pichon et al. [START_REF] Pichon | [END_REF]. Here, the effect of the oscillation in flame fronts is much more apparent, with S L,ref varying by more than 20 m/s. This is 4 times higher than the variation presented in the original study [START_REF] Mcgrath | Proceedings of the Combustion Institute 000[END_REF]. Moreover, the disagreement between the experimental and numerical data is assumed to be caused by the oscillation of flame fronts. However, different phenomena compared to mode A can be observed. All of the conditions demonstrate that on average the reference flame speed is lower for the top flame than for the bottom flame. This may be due to incomplete vaporization of the droplets, which makes the inlet mixture both leaner and colder, reducing the effective equivalence ratio and resulting in lower reaction rates. The temperature of the vapour mixture at the exit of the bottom nozzle was measured as 28.5±3 C throughout all the measurements. However, the temperature of the spray mixture at the exit of the top nozzle varied between 8.7-14.1 C, with colder temperatures present for the richer conditions. Condition f = 1.21 shows the largest disagreement between the two flames. This condition represents the largest temperature difference between the two mixtures. Incomplete vaporiza-tion of droplets in the top flame under rich conditions leads to a leaner flame, and an increase in the reference velocity. The droplets from the top flame may also be traversing across the two flame fronts [START_REF] Mcgrath | Proceedings of the Combustion Institute 000[END_REF], delivering additional fuel to the bottom mixture, in turn making it richer and colder, thus reducing the reference flame speed. However, the contrary is observed, with the reference velocities for the top flame being lower than for the bottom. In order to determine the dominant effect for this discrepancy, PDPA measurements are planned to calculate the liquid fraction of acetone in the mixture, allowing to obtain the effective equivalence ratio and track the droplet population distribution throughout the flame fronts.

Results of gas phase thermometry and species composition

Preliminary measurements using SRS were used to determine the temperature and species composition in reactant and product zones of the counterflow methane flames (mode A). Measurements were performed at nozzle exit velocities of 1.73, 2.17 and 2.60 m/s for the equivalence ratios of 0.8, 1.0 and 1.2. The tempera-Fig. 4. Axial temperature and species molar fraction concentration measurements for f = 1.0, v ex = 2.17 m/s (K ⇡ 440s 1 ), counterflow methane flames. ture and species composition values were determined by least-square fitting of the convolution of theoretical spectra of each species of interest with the instrument functions to the single-shot experimental spectra. Further details of the exact procedure and databases of fitting parameters are available in Ref. [START_REF] Guichard | Proceedings of the European Combustion Meeting[END_REF].

Fig. 4 shows the experimentally determined profiles of temperature and species molar concentrations for f = 1.0 counterflow methane flames at a nozzle velocity of 2.17 m/s. The measurements were extracted from the average of the values acquired by fitting individual single-shot spectra. The values for the top and bottom flames are shown from left to right respectively. It can be seen that the profiles are symmetric, with equal temperature and species molar concentrations as expected. The measurement fluctuations represented by the error bars in the reactants low-temperature zone, are much lower when compared to those in the product gases, as the signal is much stronger in regions of higher gas density.

Fig. 5 shows a zoomed section of profiles for the top flame shown in Fig. 4, along with results of 1-D counter-Fig. 6. Experimental (symbols) and simulated (dashed lines) molar concentration profiles of CO (red) and CO 2 (blue) and a sum of CO and CO 2 (purple) for counterflow methane flames with v ex = 1.73 m/s (K ⇡ 340 360 s 1 ) for different equivalence ratios. flow premixed flame Cantera simulations performed using a GRI 3.0 mechanism. The values for temperature, CH 4 , N 2 , and H 2 O are in good agreement with the simulations in both reactant and product regions. However, there is a discrepancy in the shapes of the gradients, with experimental data exhibiting sharper gradients in values than experimental values. The horizontal orientation of the two flames is parallel to that of the excitation beam. It was noticed that the axial position of the flame varied vertically within a 0.9 mm range, leading to a large variation within the high gradient zone. Very good agreement is found between simulations and measurements of temperature and H 2 O molar fractions, but not for CO and CO 2 . Fig. 6 shows the experimental and numerical profiles for molar fractions of CO, CO 2 and their sum as a function of equivalence ratio for the inlet fluid velocity of 1.73 m/s (K ⇡ 440 s 1 ). For the lean and stoichiometric conditions, measured CO molar fractions are higher and of CO 2 lower than expected from theory. The signal strengths of CO and CO 2 are very weak in the product zone, as shown in Figs. 78, leading to a high measurement uncertainty. The spectral region for CO 2 Raman signal has a high background, and the CO peak is further complicated by its overlap with the rotational Raman emission of N 2 at high temperatures. Finally, in the rich condition, (f = 1.2), the presence of strong C2 Swan band emission leads to significant overlap with the CO signal, resulting in divergence from the expected equilibrium molar fractions of around 25 %.

The SRS method was then applied for preliminary measurements on an acetone spray/methane gas flame, f = 1.0, v ex = 1.73 m/s, for the top and bottom burner respectively. The top burner was run in mode B while the bottom burner was run in mode A. At full laser energy of 1.5 J, frequent and violent optical breakdown events were observed when a droplet passed through the beam volume, as shown in Fig. 9. Shockwaves from the event were sufficiently strong to extinguish the flame or to cause ignition of the gaseous mixture. It was found that such events did not occur at energies below 0.5 J, which was set as the laser energy from then on. Fig. 10 shows two single instantaneous spectra of acetone spray flame, f = 1.0, v ex = 1.73 m/s (K ⇡ 347 s 1 ), 4 mm from the spray nozzle exit, in the reactant zone, separated in time by 3 seconds. Expected signals of N 2 and O 2 are observed, along with the stretching mode CH 3 of acetone at 2750 cm 1 . For the majority of 500 spectra in the series, the magnitude of CH 3 mode remained lower than that of the N 2 , indicating that it comes from the acetone vapor that was created by the droplet evaporation. Sudden spikes of the CH 3 emission, along with unknown mode at 3596 cm 1 are observed, as shown in Fig. 11 which arise from the droplets crossing the beam. The magnitude of these modes varies significantly with time. The large magnitude of the CH 3 peak may be linked to a sudden increase of acetone molecular density in a droplet, resulting in liquid phase Raman, or laser-induced droplet evaporation, leading to locally high concentrations. Another strong mode is observed at 1711 cm 1 , which also varies in magnitude but independent of the CH 3 mode, as shown in Fig. 11. Its origin is currently being investigated. Fig. 12 shows an average of 500 spectra of the acetone spray/ methane gas flame, f = 1.0, v ex = 1.73 m/s (K ⇡ 347 s 1 ), 12 mm from the spray nozzle exit, prod-

Conclusions

The effect of fuel droplets on the burning velocity of strained laminar premixed flames was investigated experimentally using optical diagnostics. The twin counterflow burner configuration was used to create premixed acetone vapour flames with the addition of fuel droplets.

Particle image velocimetry (PIV) was used to measure 2D velocities for nominal equivalence ratios in the range (0.8-1.21) as a function of strain rate (250-550 s 1 ). Measurements of reference flame speeds upstream of the flames were made for both reference methane/air flames and acetone vapour/droplet flames, and compared to simulated values of the purely gaseous flames. Experimental values of the methane/air flame were in agreement with simulated data for the bottom flame. The values for the top flame had the same gradient but were vertically offset as the top mixture temperature was higher due to heating transfer from the exhaust gases to the burner. The measurements of the acetone spray/vapour flame have shown an inverse relation. The latent heat of droplet evaporation significantly reduced the mixture temperature of the top flame compared to the bottom, resulting in lower reference velocity values, respectively. The discrepancy in the reference velocity values of the two flames grew with an increased nominal equivalence ratio due to a higher fraction of liquid droplets. Spontaneous Raman spectroscopy was used to measure temperature and species molar fractions across reference methane flames. The molar fraction and temperature results were shown to be in good agreement with the simulations for most of the species, except CO and CO 2 , which had too low a signal-to-noise ratio in the product zone. Preliminary measurements of acetone spray flame showed the capability of the setup to measure Raman signals in the presence of droplets at lowered laser energy. High laser energies led to prompt ignition of the droplets. Two types of spectra were detected, reactants with acetone vapour and reactants with an Acetone droplet crossing the probe volume. Expected products were detected past the flame front along along with strong acetone emission, showing the capability of droplets to cross the flame front.

Fig. 1 .

 1 Fig. 1. a) Configuration of the burners and flow delivery for two modes. Mode A (blue) methane/methane gas flames, mode B (red) acetone spray/vapour flames. b) Sample velocity magnitude flow-field with the extraction of the axial velocity profile along the stagnation streamline shown in c).

Fig. 2 .

 2 Fig. 2. Reference laminar flame speed S L,ref as a function of strain rate K for methane gas flames (mode A). Dark grey and light grey points correspond to single measurements of the top and bottom burners respectively. Circles and squares represent average values for top and bottom flames grouped by the v ex . Black lines correspond to Cantera simulations.

Fig. 3 .

 3 Fig.3. Reference laminar flame speed S L,ref as a function of strain rate K for acetone spray/vapour flame (mode B). Dark grey and light grey points correspond to single measurements of the top and bottom burners respectively. Circles and squares represent mean values for top and bottom flames grouped by the v ex . Black lines correspond to Cantera simulations computed using mechanism of Pichon et al.[START_REF] Pichon | [END_REF]

Fig. 5 .

 5 Fig. 5. Axial temperature and species molar fraction concentration measurements (symbols) for f = 1.0, v ex = 2.17 m/s (K ⇡ 440s 1 ) counterflow top methane flame. CH 4 , CO 2 and CO values were magnified 6 times. Dashed lines show results of 1D Cantera simulations.

Fig. 7 .

 7 Fig. 7. Mean and instantaneous Raman spectra for a counterflow methane flame at f = 1.0, v ex = 1.73 m/s (K ⇡ 360 s 1 ) in the product zone.

Fig. 8 .

 8 Fig. 8. Mean and instantaneous Raman spectra for a counterflow methane flame at f = 1.2, v ex = 1.73 m/s (K ⇡ 347 s 1 ) in the product zone.

Fig. 9 .Fig. 10 .

 910 Fig. 9.Optical breakdown of droplets in an acetone spray/methane gas flame, f = 1.0, v ex = 1.73 m/s (K ⇡ 347 s 1 ).

Fig. 11 .

 11 Fig. 11. Variation of CH 3 (top) and unknown modes (588 nm) (bottom) during a series of spectra for acetone spray/methane gas flame, f = 1.0, v ex = 1.73 m/s (K ⇡ 347 s 1 ), 4 mm from the spray nozzle exit, reactant zone.

  Fig. 12. Averaged spectrum of acetone spray/ methane gas flame, f = 1.0, v ex = 1.73 m/s, 12 mm from the spray nozzle exit, product zone.uct zone. The acetone spray and methane flames were generated by the top and bottom burner respectively. The CO 2 , H 2 O and N 2 modes are present, as expected from the product zone. However, a presence of a strong CH 3 mode is still observed, suggesting that droplets continue to pass through the flame front and produce strong Raman signal emission, which contributes to the averaged spectra. The emission at 3600 cm 1 is present only in the presence of CH 3 mode. It should be noted that the emission at 1711 cm 1 , previously observed in the reactants zone, triples in magnitude, in contrast to a nearly 7 times reduction in gas density.

		9000
		8500
		8000
		7500
	Counts	7000	CH 3 (acetone)
		700	CO 2
			H 2 O
		650
		600	1500 2000 2500 3000 3500 4000 4500
			Wavenumber (cm 1 )

Acknowledgements

Financial support from the Région Normandie under the Chaire d'Excellence DOLFIN (DrOplet-Laden Flame INteraction) is gratefully acknowledged.


Nemo Fournier
email: nemo.fournier@icm-institute.org

Stanley Durrleman

A Multimodal Disease Progression Model for Genetic Associations with Disease Dynamics

Keywords: Multimodal Disease Progression Modelling, Alzheimer's Disease, Genetic Associations

We introduce a disease progression model suited for neurodegenerative pathologies that allows to model associations between covariates and dynamic features of the disease course. We establish a statistical framework and implement an algorithm for its estimation. We show that the model is reliable and can provide uncertainty estimates of the discovered associations thanks to its Bayesian formulation. The model's interest is showcased by shining a new light on genetic associations.

Introduction

The clinical courses of neurodegenerative pathologies such as Alzheimer's or Parkinson's Diseases span multiple years and encompass intricate evolution of patients' cognitive abilities, physiological biomarkers and brain structure. Longitudinal studies are an essential tool for clinicians to uncover the diseases' mechanisms. In such studies, biomarkers and cognitive scores of patients are repeatedly measured at different times and need to be analyzed together, usually with a two-sided scope. First, to describe the general process at play across a whole cohort of patients: this is population-level modelling and allows to describe the average course of the disease. A second layer aims at explaining and predicting the variability observed among individuals: this is personalized-level modelling.

Mixed-effect frameworks are widely adopted to address these multi-layered prospects, offering to disentangle fixed effects (population level) from random effects (individual level) to explain the variability of the disease. Linear mixedeffects models are the simplest instances of such models. Generalized linear and non-linear mixed-models are now often prefered to account for the neurodegenerative diseases' peculiarities, and most state of the art disease progression models (e.g. [START_REF] Donohue | Estimating long-term multivariate progression from short-term data[END_REF][START_REF] Schiratti | A Bayesian mixedeffects model to learn trajectories of changes from repeated manifold-valued observations[END_REF][START_REF] Raket | Statistical Disease Progression Modeling in Alzheimer Disease[END_REF]) belong to these categories. They are indeed better suited to describe phenomena whose complex dynamic spans multiple years. They have been used with success to describe the natural history of diseases [START_REF] Koval | AD Course Map charts Alzheimer's disease progression[END_REF] or make individualized predictions, for instance to enrich clinical trials [START_REF] Maheux | Forecasting individual progression trajectories in Alzheimer's disease[END_REF]. A general formulation is as follows, where η is a non linear mapping between timepoints and clinical markers, parametrized by fixed-effects α and random effects β i (methods differ by the chosen non-linearity η and how α and β parametrize the disease course):

y i = η α (t i | β i )
Inter-patient variability is thereby modelled through random perturbations β i around a fixed reference α. However, it is known that some of this clinical variability between patients is explained by external factors (and thus hardly explained entirely by random perturbations). Genetic mutations or external factors such as gender, family history, education or socio-economics levels can influence the course of pathologies. Accumulating evidence suggests that the variability induced by such covariates stems from general mechanisms shared across the population, for instance in Alzheimer's and Parkison's diseases [START_REF] Livingston | Dementia prevention, intervention, and care: 2020 report of the Lancet Commission[END_REF][START_REF] Komarova | High Degree of Heterogeneity in Alzheimer's Disease Progression Patterns[END_REF][START_REF] Greenland | The clinical heterogeneity of Parkinson's disease and its therapeutic implications[END_REF][START_REF] Jutten | for the Alzheimer's Disease Neuroimaging Initiative: Finding Treatment Effects in Alzheimer Trials in the Face of Disease Progression Heterogeneity[END_REF]. In the presented models, observed covariates c i are not taken into account and only the repeated observations y i are modelled as a function of the patient's ages t i . Thus, random effects might be such that E [β i | c i ] = 0. This shows that some signal present in covariates to explain the progression of the disease has not been fully exploited.

Our contribution is to propose a slight change in this mixed-effect paradigm to allow non-linear models to also be influenced by these variables. Instead of estimating a fixed effect α (parametrizing the average disease course) as well as random effects β i , we introduce a link function f ϕ that can predict, given a set of covariates c i (e.g. sex, education level, SNP arrays, genetic risk scores), an expected trajectory of the disease conditionned by these covariates.

The main difference between the standard approach and our method is that the previously introduced fixed-effects α are now estimated for each subject as a deterministic function of their covariates f ϕ (c i ). It also differs from accounting for the heterogeneity through hierarchical progression models [START_REF] Young | Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with Subtype and Stage Inference[END_REF][START_REF] Poulet | Mixture Modeling for Identifying Subtypes in Disease Course Mapping[END_REF]) since covariates are, in our case, supervisingly used during model calibration and used to navigate through a continuum of disease models, instead of having defined clusters.

We demonstrate the value of this approach by adapting a general modelling framework, namely a non-linear Bayesian model: the Disease Course Mapping (DCM) [START_REF] Schiratti | A Bayesian mixedeffects model to learn trajectories of changes from repeated manifold-valued observations[END_REF]. We show that accounting for time-independent covariates in the longitudinal modelling with this approach can be done in a reasonable statistical setting. A stochastic estimation algorithm can be devised and we propose an instantiation and implementation of our model, which we validate first on synthetic. We then use clinical data from the Alzheimer's Disease Neuroimaging (ADNI) cohort and further demonstrate the clinical interest of the method by estimating new associations between genetics and disease dynamics.

Method

We derive here an algorithm that learns to model repeated observations while accounting for the heterogeneity explained by additional covariates.

A Generic Mixed-Effects Geometric Model

In their seminal paper [START_REF] Schiratti | A Bayesian mixedeffects model to learn trajectories of changes from repeated manifold-valued observations[END_REF], Schirrati et al. introduced a generic framework to model a dataset y i,j of multimodal longitudinal measurements. Here y i,j is a vector of N biomarkers measured for the i-th subject at their j-th visiti.e. at age t i,j . Each observation y i,j is assumed to lie on Riemannian submanifold (M, g) of R N . The average course of the disease is posited to be such that individual progressions stem from a geodesic trajectory γ 0 on the manifold surface. Geodesic equations imply that γ 0 is entirely characterized by its initial position p 0 ∈ M and speed v 0 ∈ T p0 M (tangent space of M at p 0 ) at time t 0 . Individual trajectories are obtained from this reference trajectory γ 0 via a temporal reparametrization t → ψ i (t), used to derive what we name a disease age and enables registering patient's chronological ages onto a common disease timeline. Spatial effects w i are applied to the reference trajectory thanks to an exp-parallelization procedure that identifiably deforms geodesics in the manifold space. We denote η wi γ0 the resulting geodesic. We refer to [START_REF] Schiratti | A Bayesian mixedeffects model to learn trajectories of changes from repeated manifold-valued observations[END_REF] for extensive details on the geometric properties of these operations (such as commutativity and identifiability of both temporal and spatial effects).

The choice of the manifold's metric shapes the geodesic trajectories and thus the disease model [START_REF] Schiratti | A Bayesian mixedeffects model to learn trajectories of changes from repeated manifold-valued observations[END_REF][START_REF] Gruffaz | Learning Riemannian metric for disease progression modeling[END_REF][START_REF] Sauty | Riemannian Metric Learning for Progression Modeling of Longitudinal Datasets[END_REF]. Clinical knowledge of Alzheimer's Disease suggests that sigmoid shapes as sound candidates to model biomarkers' evolutions (see [START_REF] Jack | Tracking pathophysiological processes in Alzheimer's disease: An updated hypothetical model of dynamic biomarkers[END_REF] for clinical considerations, or [START_REF] Meisl | In vivo rate-determining steps of tau seed accumulation in Alzheimer's disease[END_REF] for the logistic dynamic of imaging-derived features such as brain-averaged protein loads backed by prion-like diffusion hypothesis). We therefore consider a product-metric g p such that geodesic are sigmoids:

g p (u, v) = u • M (p) • v with M (p) = 1 p 2 (1-p) 2 ,
which gives the trajectories of equation [START_REF] Allassonnière | Construction of Bayesian deformable models via a stochastic approximation algorithm: A convergence study[END_REF]. The resulting trajectories and the geometric interpretation of the (v 0 , p 0 , t 0 ) parameters are also presented in figure 1.

η wi γ0 (ψ i (t i,j )) (k) = 1 + 1 p (k) 0 -1 exp - v (k) 0 ψ i (t i,j ) + w (k) i p (k) 0 (1 -p (k) 0 ) -1 (1)

Covariate Association and Statistical Framework

We provide here a statistical instantiation of the previous geometric model. As described, given a geodesic trajectory γ 0 (fully specified by its position p 0 and speed v 0 at initial time t 0 ), and a set of random effect ψ i and w i , individual trajectories of an individual i observed at times (t i,j ) j are modelled by the curve η wi γ0 (ψ i (t i,j )). We propose that γ 0 (which represents the reference disease course,

Biomarker 2 = Fig. 1: A two feature model, with geodesic trajectories on the manifold (left) and the biomarkers observation space (right). This provides the intuition over the effect of the initial position p 0 and the initial velocity v 0 at time t 0 .

as a fixed-effect of the model) is to be computed for each subject i from the measured covariates c i as:

γ 0,(i) ∼ = p 0,(i) , v 0,(i) , t 0,(i) = f ϕ (c i )
Where f belongs to a parametrized family of functions and ϕ are its parameters treated as the new fixed-effect of the model. The individual effects to register this computed γ 0,(i) onto observations are characterized by two random effects: an acceleration factor ξ i and a time-shift τ i such that ψ i (t) = e ξi t -t 0,(i) -τ i . On top of these are space-shifts w i ∈ R N , computed thanks to an ICA:

w i = As i ,
where A is a latent matrix of independent directions (fixed effect) and s i is the corresponding individual latent source vector (random effect).

Our hierarchical statistical model treats the fixed and random effects as a set of latent variables z which is the reunion of the population and individual variables z pop = {ϕ, A} and z indiv = {(s i ) i , (τ i ) i , (ξ i ) i }. We posit the following priors on these latent parameters, where θ hyper = {σ ϕ , σ A } are fixed hyperparameters and θ model = ϕ, A, σ 2 τ , σ 2 ξ , σ 2 are the parameters of the model to be estimated:

ϕ ∼ N ϕ, σ 2 ϕ A ∼ N A, σ 2 A ξ i ∼ N 0, σ 2 ξ τ i ∼ N 0, σ 2 τ s i ∼ N (0, 1)
A non-informative prior is used over these model parameters due to the lack of a-priori knowledge. We seek to maximize a posteriori the joint-likelihood under the following additive Gaussian noise modelling

y i,j = η Asi fϕ(ci) (ψ i (t i,j )) + ε i,j q(y, z, θ model ) = q(y | z, θ model ) data attachment × q(z | θ model , θ hyper ) regularization of the latent variables × q prior (θ model | θ hyper )
prior on model parameters (taken non informative)

It can be shown that the model's likelihood function lies in the curved exponential family. That is there exist two smooth functions Φ and Ψ functions of θ model and a measurable sufficient statistics function S(y, z) of the data and the latent realizations such that log q(y, z, θ model ) factors as:

log q (y, z, θ model ) = -Φ(θ model ) + S(y, z), Ψ (θ model )
This allows estimating our model with a Monte-Carlo Markov-Chain Stochastic Approximation version of the Expectation Maximization algorithm (MCMC-SAEM) while enjoying theoretical guarantees of convergence [START_REF] Allassonnière | Construction of Bayesian deformable models via a stochastic approximation algorithm: A convergence study[END_REF]. The expectation phase is therefore built upon a sampling scheme to sample from the posterior distribution of the latent parameters (namely Metropolis-Hastings within Gibbs sampler). The maximization phase follows update rules established by finding critical points of θ → -Φ(θ) + Sp , Ψ (θ) ( Sp is the stochastic approximation of the sufficient statistics built at step p), which yields analytic expressions.

We choose to parametrize the link function as a linear mapping between covariates and dynamic parameters f ϕ (c i ) = ϕ slope • c i + ϕ intercept . This will provide an interpretable model to explain the general processes linking covariates to dynamic features such as the base pace of the disease or average onset time. The coefficients of ϕ that correspond to the mapping between the covariates c i and v 0 measure how much a given covariate impact the progression speed of each feature, and can be analyzed easily. Model parameters are initialized by setting their intercept to the models learned by a regular DCM model without considering covariates, while latent parameters are initialized at random.

Evaluation, Clinical Results and Discussion

Simulated Data

We used the generative abilities of the DCM [START_REF] Koval | AD Course Map charts Alzheimer's disease progression[END_REF] to simulate multimodal longitudinal datasets with covariates influencing the dynamic of the progression. To this end, we fixed some reference models corresponding each to a slightly different pure form of a fictional disease. Then, covariates were simulated either:

as binary covariates that directly dictated which hardcoded model is used to simulate the repeated measurements (covariate thought as a mutation-status or sex for instance). -as continuous covariates, influencing the simulated progression by using them as convex coefficients in a combination of the reference models. The covariates are seen as continuous risk factors of following one form or another.

Our simulated datasets typically included 500 subjects with an average of 5 visits and an average follow-up duration of about 5±2 years and a measurement noise of around 5%. Such an experiment is summarized in figure 2, where three continuous covariates were simulated on the [0, 1] range, the first one being a risk to develop a motor form of the disease, the second one a memory-form risk and the third a covariate without any influence on the disease.

In this example, our calibrated model correctly matches each covariate to its simulated effect. For instance: the coefficients of the link function related to the disease initial speed v 0 (we refer to fig. 1 for its interpretation) associated with the memory-risk covariate show an acceleration of the decline in memory (multiplicative factor of 1.32 [1.15, 1.62] -credible interval at 95%) contrasted to the two other biomarkers (factor of 0.79 [0.73, 0.90] for the motor and 0.89 [0.80, 1.01] for the language). These intervals are represented in fig. 2b. These two other features are slowed down relatively to the memory in order for the model to capture the change in the slope-ratio of different features on the fixed effects. If we translate the effects of these coefficients into effects on slope-ratio, we obtain indeed obtain that the ratio between memory-speed and motor-speed goes from 4.22 [3.54, 5.02] (no extra memory-risk) to 7.05 [5.54, 11.26] (maximum extra memory-risk). The ground truth change of slope (from the reference models) was from 4.48 (no particular risk) to 9.19 (full risk) and is therefore covered by our credible intervals.

Similarly, the coefficients associated with the irrelevant covariate did not capture any significant effect (factors of 1.03 [0.94, 1.12], 1.05 [0.98, 1.12] and 1.04 [0.97, 1.11] for the memory, motor and language features), which validates the ability of the model to discard covariates without influence on the disease dynamic.

Multimodal Clinical Data

Data used in the preparation of this article were obtained from the Alzheimer's Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). ADNI was launched in 2003 as a public-private partnership led by Michael W. Weiner, MD. For up-to-date information, see www.adni-info.org. We selected subjects that eventually converted to an MCI or AD stage during their follow-up. This amounted to 1440 patients for a total of 9343 visits. The follow-up duration was 4.069 (± 3.190) years, with a baseline age of 73.683 (± 7.508) years old. We processed and included biomarkers relevant to monitor AD progression:

-Two cognitive scores: the Mini-Mental State Exam (MMSE) and the AD Assessment Scale-Cognitive (ADAS-Cog). We normalized and inverted them so that they both cover the [0, 1] interval, (1 being the highest abnormality). -Hippocampus and Ventricles volumes, measured by structural T1 MRI and normalized by patient's Intracranial Volume (ICV). As for the cognitive scores, these measurements were rescaled to [0, 1] interval. -Contrasted PET imaging derived brain-averaged amyloid β42 and phosphorylated τ proteins loads, also rescaled to [0, 1].

APOE-ε4

. We calibrated our model by including a covariate known to modulate Alzheimer's Disease course, namely the patient's APOE mutation status. The results of this model are showcased in figure 3. It shows that the APOE mutation, which is a known risk factor for AD, has a clear effect on the disease dynamic: in the obtained disease course map, the mutation is associated with earlier and faster abnormalities on most biomarkers. We also investigate the learned linked function f ϕ and its coefficients dictating the interaction between the covariates Fig. 2: Risk factors to develop an acute memory-form or acute motor-form of the disease are sampled continuously from the [0, 1] interval, and used as convex coefficients to combine three reference models (a standard form, a motor-dominant form and an acute memory form), yielding the continuum of possible trajectories presented in (a). We also sample an irrelevant covariate that is never used to modulate the disease course. In (a) is the resulting model, which we can visualize for any combination of covariates (two combinations are presented). We also plot the credible intervals (at 95%) for the coefficients linking covariates to the speed of progression on each feature (multiplicative effect, thus 1.0 stands for no influence while a coefficient of 1.5 stands for an expected progression speed greater by 50%).

and the speed of progression. This is presented in figure 3. This showcases how the contribution of the APOE to the speed of progression (as in the coefficient linking the covariate to the coordinates of v 0 ) is different among biomarkers.

Features can be grouped into a cognitive scores group, more impacted than the other features by the mutation ( Wide Association Study (GWAS) [START_REF] Kunkle | Genetic metaanalysis of diagnosed Alzheimer's disease identifies new risk loci and implicates Aβ, tau, immunity and lipid processing[END_REF]. We included them in our model as covariates. The results suggest that being associated with the diagnosis does not inform a priori on the influence of each SNP on the disease course. In figure 4 we show that, even though all these SNP were selected for being significatively associated with the diagnosis, they can exhibit differences in their association with the disease dynamic. Fig. 4: Analysis of the learned interaction between some of the SNP included in the analysis and the speed of progression of each of the 6 measured features. Some SNPs present no significant interaction with the progression speed of any of the variables (e.g. rs138727474T) to SNPs that are associated with a group of feature (e.g. cognitive domain for rs114812713C) or single features (e.g. rs11932324A or rs286604821A, in either a protective or risk-inducing direction).

M M S E

Conclusion

We proposed a framework to adapt a state of the art Bayesian non-linear mixedeffect disease progression model to capture the effects of external covariates into the disease dynamic. We implemented an estimation algorithm, and show that it reliably provides new interpretable measures of interaction between covariates and the disease course. For instance, we recover the (clinically known) association between the APOE-ε4 mutation and cognitive dysfunction. In particular, its use on genetic data (either single mutation status or SNP arrays) could help to go beyond associations with the sole diagnosis and provide complementary tools to GWAS.

Fig. 3 :

 3 Fig.3: Model of Alzheimer's Disease obtained from ADNI data. Left: estimated parameters of the link function ϕ (mode of the posterior and 95% credible interval). Right: difference in the trajectory conditioned by the mutation status (the represented trajectories are for 0 copy vs 2 copies of the APOE-ε4 allele).

A D A S 1 3 VContribution

 3 e n tr ic le s V o lu m e H ip p o c a m p u s V o lu m e P E T A m y lo id tr ic le s V o lu m e H ip p o c a m p u s V o lu m e P E T A m y lo id tr ic le s V o lu m e H ip p o c a m p u s V o lu m e P E T A m y lo id tr ic le s V o lu m e H ip p o c a m p u s V o lu m e P E T A m y lo id tr ic le s V o lu m e H ip p o c a m p u s V o lu m e P E T A m y lo id lo a d P E T T a u lo a d Effect of rs11932324_A Effect of rs114812713_C Effect of rs138190086_A Effect of rs138727474_T Effect of rs28660482_A

This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the Investissements d'avenir program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute), ANR-19-JPW2-000 (E-DADS), and ANR10-IAIHU-06 (IHU ICM), as well as by the European Research council reference ERC-678304 and the H2020 programme via grant 826421 (TVB-Cloud).