I Nascimento

D Rodrigues

N Winter

B Gicquel

M Pizza

R Rappuoli

L Leite
email: lccleite@butantan.gov.br

Adjuvanted Recombinant BCG Vaccine increases protection against TB

Vaccination against Mycobacterium tuberculosis (TB) has been only partially successful based on the vaccine M. bovis Bacillus Calmette-Gue ´ rin (BCG), protecting children efficiently against early manifestations of TB. However, its protective immunity wanes in 10-20 years, being limited against adult pulmonary TB. Therefore several strategies are being pursued for the development of improved vaccines against TB. It is considered that to obtain a high protection level against TB, it is necessary to induce a high Th1 cellular immune response against the mycobacteria. Most of the strategies attempt to increase the CD4+ -Th1 and CD8+ characteristics of the immune response. Assays evaluate the efficacy of TB vaccines in experimental animals by pulmonary challenge with MTB in mice and then in Guinea pigs. Universally, BCG provides a 1 log reduction in lung colonization in mice at 1 month, considered to be partially protective. The most promising vaccines being evaluated have been based on different forms of recombinant mycobacteria or the BCG-prime and boost with different antigens. We investigated whether the well-know adjuvant properties of bacterial toxins and toxin derivatives are maintained when expressed in BCG. We have observed that the expression of these derivatives in recombinant BCG could modulate the immune response induced against the co-administered rBCG strains towards a Th2 or Th1 pattern depending on the derivative. Therefore, we investigated the effect of the expression in rBCG of a very potent Th1-driving toxin derivative -the E. coli heat labile enterotoxin derivative, LTK63. The immunization of mice with the rBCG-LTK63 strain induced a very strong Th1-predominant immune response against the mycobacteria, with production of high levels of IFN-γ and TNF-α . This resulted in a 2 log reduction in bacterial colonization in mice following challenge with MTB, or 1 log reduction in comparison to BCG. We have now expressed LTA-K63 in BCG, which resulted in a 1.3 log reduction in relation to BCG immunization, following challenge with MTB. These results confirm that the expression of Th1 driving adjuvants in BCG can modulate the immune response induced against the mycobacteria, inducing higher protection against TB

THEME OF THE CONFERENCE

The EMBO conference series "Tuberculosis 2012" is an international congress on Tuberculosis and its causative organism Mycobacterium tuberculosis to be held from September 11-15, 2012 at the Conference Centre of the Institut Pasteur, Paris.


domain: math info

text:
HAL is

< ∞, where {F n } is an increasing sequence of σ-algebra. We assume that

n (θ 0 , ν 0 ), where g

(1) n and g

(2) n are F n-1 -measurable and may be nonlinear in andg (2)

n (θ, ν) is assumed asymptotically negligible relatively to g

(1)

n (θ). We show the strong consistency and the asymptotic normality of the Weighted Conditional Least Squares Estimators (WCLSE) of θ 0 ,

k , where λ k is F k-1 -measurable, ν is any estimation of the nuisance parameter ν 0 , The proofs are based on a new strong law of large numbers for a class of submartingales (SLLNSM). We illustrate the theoretical results with examples such as a nonstationary GARCH(1,1) model and branching processes modelling some biological phenomena.


Guillaume Vignat
email: guillaume.vignat@centralesupelec.fr

Daniel Durox

Kevin Prieur

Sébastien Candel

An Experimental Study into the Effect of Injector Pressure Loss on Self-Sustained Combustion Instabilities in a Swirled Spray Burner

Keywords: Combustion Instabilities, Swirled flames, Spray flames, Injection Dynamics, Pressure Loss

Combustion instabilities depend on a variety of parameters and operating conditions. It is known, especially in the field of liquid rocket propulsion, that the pressure loss of an injector has an effect on its dynamics and on the coupling between the combustion chamber and the fuel manifold. However, its influence is not well documented in the technical literature dealing with gas turbine combustion dynamics. Effects of changes in this key design parameter are investigated in the present article by testing different swirlers at constant thermal power on a broad range of injection velocities in a well controlled laboratory scale single injector swirled combustor using liquid fuel. The objective is to study the impact of injection pressure losses on the occurrence and level of combustion instabilities by making use of a set of injectors having nearly the same outlet velocity profiles, the same swirl number and that establish flames that are essentially identical in shape. It is found that combustion oscillations appear on a wider range of operating conditions for injectors with the highest pressure loss, but that the pressure fluctuations caused by thermoacoustic oscillations are greatest when the injector head loss is low. Four types of instabilities coupled by two modes may be distinguished: the first group features a lower frequency, arises when the injector pressure loss is low and corresponds to a weakly coupled chamber-plenum mode. The second group appears in the form of a constant amplitude limit cycle, or as bursts at a slightly higher frequency and is coupled by a chamber mode. Spontaneous switching between these two types of instabilities is also observed in a narrow domain.

Introduction

Combustion instabilities induced by an acoustic coupling are a recurring problem in high performance devices like aero-engines, liquid rocket motors and gas turbines for energy production. A substantial research effort has been dedicated to the understanding of the driving mechanisms and to the development of models, predictive tools and design control methods that could be used to avoid or suppress these dynamical phenomena. Much of the early work was motivated by spectacular failures of liquid rockets [START_REF] Sutton | Rocket propulsion elements[END_REF][START_REF] Summerfield | [END_REF][START_REF] Harrje | Liquid propellant rocket combustion instability[END_REF][START_REF] Yang | Liquid Rocket Engine Combustion Instability[END_REF]. Research in the recent period has focused on issues related to gas turbines operating in the premixed mode (see for example the reviews [START_REF] Candel | [END_REF]6]). The effort in control has mainly been aimed at increasing the acoustic damping rate in the chamber using quarter wave tubes, Helmholtz resonators, acoustic liners, bias flow perforated plates, etc. [START_REF] Yang | Liquid Rocket Engine Combustion Instability[END_REF][7][START_REF] Bellucci | ASME Turbo Expo 2001[END_REF][START_REF] Zhao | [END_REF][10][11][START_REF] Lieuwen | Combustion Instabilities In Gas Turbine Engines[END_REF]]. An alternative route has been used that consisted in changing the acoustic mode structure and resonance characteristics of the system by placing baffles in the chamber. This method has been widely used in rocket engines to remove transverse oscillations (see for example [START_REF] Oefelein | [END_REF]). A segmentation of the annular plenum of the combustor has also been used to reduce instabilities coupled by azimuthal modes in an industrially produced gas turbine [START_REF] Tiribuzi | ASME Turbo Expo 2007[END_REF]. Other studies have explored the possibility of changing the flame pattern and flame dynamics to reduce their sensitivity to perturbations and diminish the gain of the process that drives the instability, leading to the suppression of the unstable acoustic-combustion coupling [START_REF] Lieuwen | Combustion Instabilities In Gas Turbine Engines[END_REF][START_REF] Noiray | [END_REF][16][17]. A significant research effort has also been made to develop active control methods for combustion instability reduction [18,19], with successful laboratory scale demonstrations and industrial application in some gas turbines [START_REF] Lieuwen | Combustion Instabilities In Gas Turbine Engines[END_REF] but these methods have generally been considered to be too complex for a generalized implementation.

One method that has been extensively exploited in rocket engines has consisted in implementing a pressure drop in the injection elements to decouple the upstream manifold dynamics from those of the thrust chamber [START_REF] Sutton | Rocket propulsion elements[END_REF][START_REF] Summerfield | [END_REF][START_REF] Harrje | Liquid propellant rocket combustion instability[END_REF]. Relative injector pressure drop values of the order of 10 to 15% are typically adopted in engineering practice. Such head losses cannot be used in gas turbines, because of their impact on the overall engine thermodynamic efficiency. It is however important to examine effects of injector pressure losses on combustion instability and see how this parameter influences the stability of the system, an issue that has not been extensively investigated. In gas turbine combustors, the injection pressure of the liquid fuel is much higher than the chamber pressure to promote atomization and one may assume that the fuel manifold is insensitive to pressure perturbations in the chamber. However, pressure losses in the air stream have relatively low values of the order of a few percent (typically 3 to 4%). When modifying the head loss level in a swirled premixed burner by varying the pressure drop across the swirling injector Polifke et al. [20] observed that for certain operating conditions, an increase in air flow rate in the injector led to the triggering of an instability. They also noted that this corresponded to a change in the slope of the operating line linking the air flow rate to the pressure drop and attributed this effect to a change in the recirculation zone at the injector outlet and thus to a change of flame geometry that in turn led to a change in the flame response to perturbation promoting instability.

The question that can be raised is whether an increase in the pressure drop of the air stream across injectors, without any significant change in flame shape and dynamics, can change the stability characteristics of the system possibly reducing or suppressing instabilities. This question is investigated in the present article. It is of course linked to the possible acoustic coupling between the upstream plenum and the combustion chamber. It is known that these cavities are weakly coupled if there are large changes in cross sections [21][START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Bourgouin | ASME Turbo Expo 2014[END_REF]. The coupling may also be changed by the injector pressure drop that modifies the injector impedance. It is then worth considering injectors with different head losses but similar swirl numbers conserving the flame shape and examining the flame dynamics and corresponding self-sustained oscillations. A single injector configuration is used in the present investigation (designated as "SICCA-Spray") but the analysis is carried out in a framework where tests on this geometry are used to prepare investigations on a more complex multiple injector annular configuration "MICCA-Spray" at EM2C laboratory. Strong azimuthally coupled combustion instabilities have already been observed in MICCA-Spray [START_REF] Prieur | [END_REF] and it was shown that injectors operating near the pressure antinodes were essentially driven by velocity perturbations induced by the injector response to the azimuthal mode. It is then logical to use the iteration loop shown in Fig. 1 to characterize the dynamics of injectors, measure flame describing functions (FDF) using acoustic forcing or examine the dynamics of the system under self-sustained, longitudinally coupled, combustion oscillations. The FDF can then be included in the stability analysis as exemplified in [25] in a configuration featuring longitudinal acoustic modes or implemented in combination with a Helmholtz solver as explored in [26]. The pressure loss across the injector is here considered as a design parameter that can be varied to promote or reduce unstable oscillations. Results obtained in single element tests can be used to guide experimentation on the MICCA-Spray annular system. This article describes experiments carried out only on the single injector spray swirled burner SICCA-Spray. This configuration is briefly described in Sec. 2. Four different injectors are designed (Sec. 3) with the same swirl numbers, in order to have identical flame shapes, but different pressure drops (up to 30% difference). A stability map of this system is constructed in Sec. 4 by keeping the thermal power constant and exploring a broad range of injection velocities. Results are interpreted in Sec. 5.

Experimental Setup

The SICCA-Spray experiment [START_REF] Prieur | Thermoacoustic instabilities in gas turbines and rocket engines: industry meets academia[END_REF] is sketched in Fig. 2. This single injector version of the MICCA-Spray annular combustion chamber [START_REF] Prieur | [END_REF] comprises a combustion chamber, a plenum, fed with compressed air, and an injector fed with air from the plenum and with liquid n-heptane fuel. The chamber is formed by a quartz tube of length l c = 315 mm. The base of the chamber has a 15 mm high section made of brass. The injection system comprises a main body (grey), an interchangeable swirler and an atomizer for the liquid fuel. The atomizer produces a hollow cone shaped spray of liquid fuel droplets. Their Sauter mean diameter measured 2.5 mm above the injector outlet is d 32 = 35 µm [START_REF] Prieur | Thermoacoustic instabilities in gas turbines and rocket engines: industry meets academia[END_REF]. The injector ends in a conical convergent section, with an exit radius of R in j = 4 mm. The swirler is manufactured using a stereolithography 3D printer. It consists in a hollow cylinder with six tangential cylindrical channels. Key dimensions and notations are shown in Fig. 3. Two Monacor SP-6/108 Pro driver units are mounted on the plenum. They are not used in this study, and their electrical circuit is left open. A straight tubular section upstream of the injector is used for measurements in the plenum. The fuel mass flow rate is measured using a Bronkhorst mini CORI-FLOW mass flow meter of relative accuracy ±0.4 %, and similarly, the air mass flow rate is measured using a Bronkhorst El-Flow thermal mass flow meter of relative accuracy ±0.9 %. In this article, the bulk velocity u b = ṁ/(πρR 2 in j ) is defined as the 1D velocity at the injector outlet under cold flow conditions.

Head loss is measured using a Kimo MP111 differential pressure gauge, with a 1% relative precision. Three Bruel & Kjaer type 4938 microphones with type 2670 preamplifier detect acoustic pressure signals. Their relative accuracy is 1% and their cutoff frequency of 70 kHz is far greater than the sampling frequency of 16 384 Hz used for data acquisition. Microphone MC1 is used to measure the pres-

Design of the injectors

Four injection systems were designed by changing the swirler units, while trying to keep the same flame shape and outlet velocity profiles as the reference injector used in [START_REF] Prieur | [END_REF]. To do so, a large parameter sweep was performed on parameters d and R 0 of the swirlers (see Fig. 3). A total of 38 swirlers were manufactured, the velocity profile at their outlet were measured using the LDA system and the flame shape was determined using CH* chemiluminescence. Their experimental swirl number S is calculated using the conventional expression

S = 2R in j 0 U θ U x r 2 dr R in j 2R in j 0 U x 2 rdr (1)
Among the 38 swirlers, four were selected and named S a to S d . Table 1 gathers key geometrical characteris- tics and aerodynamic performances. These four swirlers have very similar swirl numbers, and a nearly identical outlet mean velocity profile, shown in Fig. 4. The RMS velocity profiles (not presented here) show only a very small difference between the different swirlers. The injector head loss ∆P is shown in Fig. 5 for each swirler. Its evolution with the bulk velocity follows the usual scaling law: ∆P = 1/2 ρσu 2 b . The value of σ for each swirler is indicated in Tab. 1.

Figure 6 shows a direct visualization of CH* chemiluminescence from the four flames at (P = 6.2 kW, φ = 0.98, u b = 36 m s -1 ). These four images are obtained by averaging over 900 frames with 0.1 ms exposure time each. It shows that the four flames have similar "M"shapes and this has been verified for wide range of operating points.

Different types of combustion dynamics

The presence of self-sustained combustion instabilities in the SICCA-Spray experiment is investigated for each of these injection systems. The thermal power is kept constant at P = 6.2 kW, while the air mass flow rate is varied, so that u b varies between 29 and 51 m s -1 , corresponding to global equivalence ratios between φ = 1.18 and 0.69. The acoustic pressure near the combustor backplane p c is measured by microphone MC1.

Five distinct types of dynamical regimes may be distinguished depending on the injection system and operating point: closest to 445 Hz, but it does vary on a wide range depending on the operating conditions, reaching as low as 425 Hz.

Figure 8 shows the type of instability that was observed for each swirler in SICCA-Spray when varying the air mass flow rate while keeping the power constant at P = 6.2 kW. Great care was taken to wait for the experiment to reach thermal equilibrium in order to eliminate the thermal environment impact on stability, an effect that is well known (see for example [START_REF] Prieur | [END_REF]). No hysteresis phenomena were observed.

Swirler S a with the highest head loss is unstable over a wider range of operating conditions. Only type (II) and (III) instabilities at f = 445 Hz are observed with this swirler. The burst regime appears for both the leanest and the richest operating points. As the head loss decreases with swirlers S b , S c and S d , the range of unstable operating conditions narrows. For swirler S b , the first oscillations are in the burst regime, before switching to a limit cycle at 418 Hz near the stoichiometry. However, the instability rapidly switches to 445 Hz as the air flow rate is increased. For swirler S c , the first unstable point appears at stoichiometric conditions, and a limit cycle at 418 Hz is observed. This type (I) instability arises over a broad range of operating conditions. Switching between 418 Hz and 445 Hz limit cycles takes place for φ = 0.83. For leaner operating conditions, as for the other swirlers, a 445 Hz limit cycle is established giving rise to bursts at lower air flow rates. Swirler S d , compared to swirler S c , behaves similarly, but has a slightly narrower range of unstable conditions.

Figure 9 shows the amplitude of the acoustic pressure near the combustor backplane for the operating points appearing in Fig. 8. The level is always higher when the thermoacoustic oscillation frequency is around 418 Hz (red squares on Fig. 9). It reaches a peak level between 3400 and 4200 Pa while it is always inferior to 2400 Pa when the frequency is around 445 Hz (blue circles in Fig. 9).

Discussion

The coexistence of two types of combustion instabilities in a system can often be attributed to the coexistence of two flames having a different shapes and dynamical responses [20,[START_REF] Hermeth | [END_REF]29]. To investigate this point, phaseaveraged CH* chemiluminescence imaging of swirler S c at operating condition (P = 6.2 kW, φ = 0.83, u b = 42 m s -1 ) was used. At this operating point, type (IV) instabilities with alternating limit cycles were observed (see Fig. 7).

In Fig. 10, the phase averaging is conditioned with respect to the thermoacoustic oscillation signal. This allows a direct comparison of flame shapes between the two combustion instabilities at a single operating point. The flames appear identical in the two types of oscillations, indicating that the coexistence of the two instabilities is not caused by a change in flame shape.

It is then natural to inquire if this may be linked to the acoustic eigenmodes of the system by making use of the Helmholtz solver of COMSOL Multiphysics. The geometry includes an end correction for the quartz con- finement tube ∆l = 0.4 × d c where d c is the diameter of the tube (see Fig. 2). At that distance from the chamber exit, the acoustic pressure vanishes. The passive loudspeakers in the plenum are not well characterized acoustically. In order to evaluate the influence of the plenum on the acoustic structure of the system, the acoustic impedance of the whole plenum was experimentally measured by making use of the three microphone method described by Chung and Blaser [30]. A third microphone replaces the hot wire in this experiment, and the combustion chamber and injector are replaced by a loudspeaker mounted at the end of a tube used to generate plane waves in the system. The impedance is reconstructed in a plane defined by the hot wire position in Fig. 2. The plenum can then be replaced by an impedance boundary condition in the simulations. All other boundary conditions correspond to rigid walls. The speed of sound and fluid density are assumed to be uniform in the chamber and the plenum respectively.

Acoustic response simulations are compared with experimental data in Tab. 2 for swirler S d . In these experiments, a loudspeaker placed next to the SICCA-Spray setup provides acoustic excitation. In a first case, in the left column of Tab. 2, the experiment is carried out with air at room temperature in the absence of flow. In the second (third column of Tab. 2), the combustion chamber is filled with a mixture of helium and air in a suitable proportion that yields a speed of sound coinciding with that prevailing under hot fire conditions. Experimental data (columns 1 and 3) are in reasonably good agree-Figure 10: Phase averaged CH* chemiluminescence image of the flame in type (V) instability. The PiMax4 camera is triggered using the MP2 microphone signal. The averaging is conditioned by the frequency of the limit cycle oscillations, with the 418 Hz limit cycle at the top, and the 445 Hz limit cycle at the bottom. As the oscillation stays at each frequency for several seconds between switchings, well separated acquisition can be easily performed for each frequency. The images are obtained in the same manner as in Fig. 6. Swirler S c at (P = 6.2 kW, φ = 0.83, u b = 42 m s -1 ).

Table 2: Acoustic simulation for swirler S d : acoustic modes are determined experimentally by exciting the SICCA-Spray experiment with an external loudspeaker. Three operating conditions are considered: in the first column of the table, the plenum and the combustor are filled with air at room temperature, in the third column, the combustor is filled with a mixture of 24.4% air and 75.6% helium by volume such that the sound speed in this mixture is identical to that of the average temperature of the hot gases in the chamber (1000 K). Column 5 corresponds to acoustic simulations performed by considering that the combustor is filled with combustion products at a mean temperature of 1000 K. The experimental data resolution is ∆ f = ±14 Hz.

Air

Air ment with simulation results (respectively columns 2 and 4). The frequency differences are less than 30 Hz in these two cases. One may then calculate the eigenfrequencies under hot fire conditions by considering swirler S d and assuming a mean temperature T = 1000 K in the chamber. This yields two eigenmodes at 461 Hz and at 478 Hz. These frequencies are not quite those observed experimentally but this could be expected since the instability frequencies are displaced with respect to the modal frequencies by the flame response. The two modes have nearly the same spatial structure. However the lowest eigenfrequency (461 Hz) corresponds to a mode that is coupled with the upstream plenum while the higher eigenfrequency (478 Hz) essentially pertains to a chamber mode (see the details Fig. 11). The higher amplitude of the acoustic pressure in the plenum is confirmed in Fig. 12 where the axial evolution of the acoustic pressure amplitude is shown. This acoustic simulation indicates that two closely spaced acoustic modes coexist in the system. The lower frequency mode has a higher relative amplitude in the plenum. When the coupling between the plenum and the combustion chamber is strengthened, that is when the injector pressure loss is low, instabilities can more easily lock on this lower frequency mode, coupled with the plenum, leading to type (I) and type (IV) oscillations. In the experiment this pertains to the lower frequency oscillations at 418 Hz. When the injector pressure loss is higher, coupling between the cavities is diminished, leading to instabilities preferentially locking on the higher frequency, inducing type (II), and (III) combustion instabilities. These are manifested in the experiment at a frequency 445 Hz with in some cases a possible switching between the two eigenfrequencies (type (IV) oscillations).

Conclusion

Effects of injector pressure losses are investigated in this article to see how changes in this parameter modify the combustion dynamics of a generic system. Four injectors are tested in a single element combustor fed with liquid fuel and air. The velocity profiles at the injector outlets are quite similar and the corresponding swirl numbers are nearly the same. Thus the steady flame shapes obtained with these four injectors are close to each other. However, their pressure losses vary by about 30%.

This setup allows to study the influence of the head loss as a design parameter while keeping the other conditions constant. The instability map of the burner is determined by examining the dynamical regimes and recording the various types of oscillations. The pressure loss in the injector has an impact on the occurrence, type and intensity of thermoacoustic instabilities. In this system, the range of global equivalence ratio corresponding to an unstable regime is wider for injectors with lower pressure loss. Four types of instabilities, coupled by two acoustic modes (418 Hz and 445 Hz) are observed. Calculations of the acoustic eigenmodes of the system indicate that two modes, that are close in frequency, coexist in this system. Both appear to have a quarter wave like structure. Oscillations at a frequency around 418 Hz, generally occurs when the injector pressure loss is low, and corresponds to a mode that weakly couples the chamber with the plenum. The second, at a frequency around 445 Hz, arises when the head loss is larger, and corresponds to a chamber mode. In a narrow domain, one also observes a random switching between the two modes. This bistable behavior might have been attributed to a change in flame response. However the flame shape is essentially identical in both types of instabilities and switching is linked to the existence of the two eigenmodes.

Figure 1 :

 1 Figure 1: Test methodology combining experiments on a single injector system and on the annular model scale combustor MICCA-Spray. The objective is to guide experimentation by tailoring the dynamical response.

Figure 2 :Figure 3 :

 23 Figure 2: Schematic representation of the SICCA-Spray setup.

  backplane. It is mounted on a water-cooled wave guide. This microphone is at 290 mm from the chamber backplane. Microphones MP1 and MP2 are mounted flush to the wall of a straight tubular section of the plenum upstream of the injector. They are separated by 50 mm. Combined with a Dantec Dynamics miniCTA constant temperature hot wire anemometer (cutoff frequency 10 kHz), they are used to measure acoustic pressure and velocity upstream of the injector. A Dantec Dynamics FlowExplorer two component Laser Doppler anemometer (LDA) measures velocities 2.5 mm downstream the combustor backplane, under non-reactive conditions and without the quartz tube. The theoretical optical probe volume of this LDA system measures 0.14 mm × 0.14 mm × 0.23 mm. Measurements are spaced every 0.25 mm. Uncertainty on the measured mean velocities, taking statistical uncertainties into account, is less than 0.8 m s -1 . Finally, a Princeton Instruments PiMax4 camera is used for flame visualization. An optical bandpass filter centered on λ = 430 nm is used to restrain visualization to light emitted by CH*.

Figure 4 :

 4 Figure 4: Mean axial (U z , top) and azimuthal (U θ , bottom) velocity profiles measured 2.5 mm downstream of the injector under cold flow conditions outlet using LDA. u b = 43 m s -1 .

Figure 5 :

 5 Figure 5: Head loss ∆P of each swirler for the range of bulk air velocity used in this study. It is the difference between the pressure in the plenum, measured at position MP2, and the ambient pressure, in the absence of combustion.

Figure 6 :

 6 Figure 6: CH* chemiluminescence image of the stable flame obtained by reducing the combustion chamber length to 165 mm. The lineof-sight integrated images obtained by the PiMax4 camera are Abel transformed. Yellow and white correspond to high light intensities while gray corresponds to no light. No information is available near the backplane due to the opaque brass section at the bottom of the combustion chamber. Positions of the quartz confinement tube and conical injector outlet are also indicated. Swirler S a to S d . Operating condition: (P = 6.2 kW, φ = 0.98, u b = 36 m s -1 ).

Figure 7 :

 7 Figure 7: Typical pressure records p c measured by microphone MC1 at the combustor backplane illustrating type (I) to (IV) instabilities, as indicated on the left of each graph. The power is kept constant P = 6.2 kW. Type (I) instability is mainly observed with swirlers S b , S c and S d ; type (II) instability is observed with the four swirlers; type (III) instability is observed with the four swirlers near the borders of the unstable domain; type (IV) instability is only observed with swirler S c .

Figure 8 :

 8 Figure 8: Instability map showing the different types of combustion instabilities observed as a function of equivalence ratio φ and head loss coefficient σ of the swirlers. Red square: Limit cycle at 418 Hz, Open blue circle: Limit cycle at 445 Hz, Solid blue circle: Bursts, Red diamond: Switching between type I and type II, Dot: No oscillation.

Figure 9 :

 9 Figure 9: Acoustic pressure amplitudes measured at the combustor backplane during thermoacoustic oscillations. Symbols are identical to those used in Fig. 8.

Figure 11 :

 11 Figure 11: Acoustic pressure amplitude obtained in the simulations. The simulated domain is shown on the left for a frequency 461 Hz. On the right, the injector is shown in detail for the two eigenfrequencies of interest. The acoustic pressure is normalized by the acoustic pressure at the pressure antinode of each mode. Color version available online.

Figure 12 :

 12 Figure 12: Axial evolution of the the acoustic pressure amplitude for the two modes. The acoustic pressure is normalized by its maximum in order to have the same scale for both curves. 461 Hz mode in red, 478 Hz mode in blue. z = 0 corresponds to the combustor backplane.

Table 1 :

 1 Dimensions and characteristics of the swirlers considered in this study.

  & Helium Hot fire Exp. Sim. Exp.

				Sim.	Sim.
	250	248	260	246	240
	284	273	434	464	461
			490	480	478
	534	541	657	646	647
	744	744	914	892

Acknowledgements

This work benefited from the support of SafranTech, and of project FASMIC ANR-16-CE22-0013 of the French National Research Agency (ANR).


Mathias Christina

Mickael Mézino

Lionel Le Mézo

Sandrine Auzoux

Jean-Christophe Soulié

Matthew Jones

Christophe Poser

Pascal Marnotte

Eward Gerardeaux

Le Mezo

Modeling the influence of extreme climate variability on sugarcane yield gap in Reunion Island

Keywords: extreme climate, Sugarcane, Yield gap, Crop model, Reunion

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The key role of crop models is to help understand and predict the effects and interactions between climate, soil and management on crop development and yield. Crop models can provide insights into the main emerging agricultural challenges such as food security, sustainability and how to cope with the possible negative effects of climate change (Spiertz, 2012). Among these challenge, extreme weather events are expected to increase worldwide and are expected to reduce crop production (Field et al., 2012;Porter et al., 2014). The already increasing incidence and intensity of extreme climate events are seen as the principle immediate threat to global crop production (Meehl et al., 2000;Rosenzweig e tal., 2001;Olesen et al., 2007;Urban et al., 2012;Min et al., 2011;Lobell et al., 2008Lobell et al., , 2013)). Therefore, anticipating and calculating their effect on crop yields is important for topics ranging from food security to the economic viability of biomass products.

Among other crops, sugarcane represents a major source of sugar, ethanol and material for electric energy productions in tropics and sub-tropics. In the energy sector, the demand for ethanol and electricity is likely to increase, as it could replace fossil fuels and reduce greenhouse gas emissions (Goldemberg et al., 2014;Leal et al., 2013). Due to their agroclimatic zones, sugarcane plantations faced various extreme events such as extreme temperature (e.g. India, Mall et al., 2016), drought (e.g. Australie, Inman-Bamber et al., 2005), flood (e.g. Florida, Gilbert 2008) or storm (e.g. Beuzelin et al., 2009), which could be responsible for strong yield loss (e.g. 54-64% under 3-month flood, Gilbert et al., 2008) up to cane failure (e.g. under drought, Inman-Bamber et al, 2005).

Yield Gap (YG) analysis, based on the differences between different crop yield types, provides a basis for identifying the main factors affecting crop yield and providing valuable information to guide the interventions of growers, researchers, public agencies and stakeholders (van Ittersum et al., 2013). The following types of crop yields can be considered: the potential yield (YP), which is the yield of a given cultivar only determined by solar radiation, air temperature, photoperiod, CO2 concentration, genotype and plant populations ; the water-limited yield (YW) which to the same than YP, but additionally influenced by water deficit resulting from the rainfall amount and distribution along the crop cycle, as well as soil water holding capacity and evapotranspiration; and finally the actual or average yield (YA) which is defined as the yield actually obtained by farmers, considering fertilizer limitations as well as soil fertility, pests and diseases damages. According to Dias and Senthelas (2018), the difference between YP and YW results in the YG caused by water deficit (YGW), whereas the difference between YW and YA gives the YG caused by crop management and hazards (YGM).

While the use of crop models as a tool to analyze temporal and spatial variation of yield gaps appears to be straightforward, relatively few examples of this type of application can be found in literature for sugarcane. Most of these studies aimed to better separate the effect of climate and management on sugarcane yield (Garside et al., 1997 ;Jones et al., 2015 ;Cheeroo-Nayamuth et al., 2011;Van den Berg and Singels, 2013) but, to our knowledge, none of them have tried to link extreme events and yield gap.

Reunion island is a particularly adapted agro-climatic zone to assess the impact of extreme weather events variability on sugarcane yield gap. Around 23 thousand ha are cultivated with sugarcane in Reunion Island, representing about 58% of total agricultural areas (DAAF, 2015). Planted area with sugarcane ranged from sea level m to >1000 m above sea level resulting in large differences in temperature, radiation or wind conditions. The Reunion Island also presents a very strong precipitation variability over the island, with very low (<30 mm y -1 ) or high precipitation (>5000 mm y -1 ) depending on location. Moreover, growers face cyclones regularly, as well as fluctuating dry periods during winter.

Considering that yield gap studies are of high importance for guiding strategic actions to improve crop yields, in particular for tropical islands with such diverse environmental conditions and extreme events, the objectives of this study were i) to assess the magnitude of sugarcane yield gaps due to water deficit and crop management in Reunion Island over the last 20 years and, ii) to assess the impact of climate variability and extreme events intensity on sugarcane yield gap.

Material & Method

Cane management data over Reunion Island

Reunion Island is commonly divided into 5 agro-climatic zones (Fig. 1). As presented in Fig. 1 and Fig. S1, two of them are rainfed in majority: Beaufonds (East part of the island) and Bois Rouge (North East) while the others are irrigated for around 50% of the area: Savanna (North West) and Le Gol (South West) and one is an intermediary zone with 20% of the area being irrigated: Grand Bois (South).

A detailed identification of 10,403 sugarcane fields was made over Reunion Island (Fig. 1d), which should represent nearly all the actual sugarcane fields in 2017. The identification was performed by combining a remote-sensing approach and the use of databases from several institutes (Centre Technique Interprofessionnel de la Canne et du Sucre, CTICS) or state agency (Direction de l'Alimentation, de l'Agriculture et de la Forêt,). At the end, all cane fields were associated with a planted area (m²), an irrigation management (drip or spray irrigation when irrigated) and planted cane varieties. The irrigation management was obtained from the CTICS database (2008, resumed in Fig S1). When not identified, the irrigation management was assumed to be a spray irrigation in Savanna, Le Gol and Grands Bois, and was assumed not to be irrigated in the east part of the island (Beaufonds and Bois Rouge, personal communication with the CTICS). The sugarcane varieties in each field were obtained from the DAAF database in 2017 (resumed in Fig. S1). Between 20 and 40% of sugarcane fields were not identified with a given variety due to a lack of information. When we did not know the variety, we assumed that the planted variety was the R570 in Savanna, Le Gol and Grand Bois, and the R579 in Beaufonds and Bois Rouge. We choose these two varieties when no data were available because there are the most planted variety in these agro-climatic zones (CTICS, personal communications).

Climate data

Based on the annual precipitation and altitude range maps (Fig. 1a,b), 33 climatic areas where defined over Reunion Island. These maps were obtained through the AWARE Atlas (https://smartis.re/p/AWARE) developed by the CIRAD institute in Reunion Island. In each climatic zone, the centroid was assessed and daily climatic data (precipitation -PPT, global radiation -RG, potential evapotranspiration -ETP, minimum -TMIN and maximum temperature -TMAX) were obtained using the METEOR software (https://smartis.re/METEOR). The software interpolated the daily weather data at requested localization based on a station network from the CIRAD and Météo-France (Fig. S2). Daily weather data in each of the 33 climatic zones were interpolated from 1998 to 2018. Wind speed data used for the storm indice calculation (see section 2.5) were not interpolated but obtained in 5 meteorological stations (one per agro-climatic zone) at similar distance and altitude from the coast (except for the North-West zone where cane are grown in altitude, Fig. S2). A detailed soil water holding capacity map -SWC (Fig. 1c) was also obtained from the AWARE Atlas. Soil water holding capacity over Reunion Island range from 10 to 200 mm in cane areas. In the current study, PPT, TMEAN (daily average temperature), ETP and RG were averaged per year in each cane field. Then, the average annual climatic data and SWC in a given agro-climatic zone was calculated as the average field-scale data weighted by the field planted areas. Additionally, a climatic variable integrating the sum of annual precipitation and the quantity of irrigated water was defined (PPTIRR).

MOSICAS model presentation

The MOSICAS model (Martiné and Todoroff, 2002) simulates the daily sugarcane growth at a field scale for given management scenarios. MOSICAS, a deterministic thermoradiative type model that accounts for water stress, consists of growth and carbon balance modules linked to a water balance module. The modeling of Leaf Area Index (LAI) is a 'big leaf' type, and is based on the balance between daily increase (growth) and daily decrease (senescence) in leaf surface area. This balance is mainly driven by temperature and is limited by radiation with respect to leaf mass and water stress. MOSICAS converts the daily intercepted photosynthetically active radiation (PAR) into daily gain in total dry mass. This conversion process is affected by temperature and water stress coefficients, as well as a maintenance respiration coefficient and the variety maximum radiation conversion efficiency. The partitioning of total dry mass gain to root and aboveground dry mass gain is regulated by dynamic allometric fractions that depend on thermal time. Then, the partitioning of daily increase in aboveground dry mass to daily increase in millable stalk dry mass is regulated by dynamic allometric fractions that depend on the amount of aboveground biomass. The partitioning of daily increase in millable stalk mass to daily increase in structures is regulated by a dynamic allometric fraction that depends on the amount of existing millable stalk mass, temperature and water stress indexes. The remaining biomass not allocated to structures is stored as stalk sucrose. Above-ground and stalk water content is assumed constant (parameter) at the beginning of the growth and decreased with thermal time. The model code, executable and manual are open-source (General Public Licence) and available on Github repositories (https:// https://github.com/ChristinaMathias/MOSICAS). The main model equations driving cane growth are described in supplementary information (S3).

MOSICAS calibration

Cane varieties were calibrated using the ECOFI open database (Christina et al., 2020b). The ECOFI database includes 58 trials in total (60 varieties) carried out by CIRAD and partners from 1988 to the present, in Reunion Island and a few other countries (e.g. Guadeloupe, South Africa). The link between the ECOFI database and the MOSICAS model was performed using R script to create model input (also available in https://github.com/ChristinaMathias/MOSICAS). The model was calibrated using various measurements (Fig. S4): LAI, fraction of intercepted radiation (faPAR), aboveground dry and fresh mass, stalk dry and fresh mass, as well as aboveground and stalk water content. Calibration of parameters followed an iterative hierarchical process where priorities of decreasing order were the leaf area dynamics and interception, radiation conversion into biomass, biomass partitioning and sensitivity to water stress (Supplementary Material S3 -Calibration Procedure and Results, Fig. S3). The R570, R579, R585 and R586 cane varieties were calibrated using 20, 3, 9 and 9 ratoon crop cycles, respectively. In total, 21 parameters were calibrated. The LAI base temperature of 3 other varieties (R581, R582 and R583) were obtained from literature as used in the web version of the model (https://smartis.re/MOSIWEB). As no experimental data were available to parameterized other planted variety in Reunion Island, we assumed that additional sugarcane varieties were identical to the R570 or R579, as it is historically the most planted varieties in the island. Details on parameter values for the different varieties are given in Table S1.

Sugarcane yield simulations & observations

Actual fresh sugarcane yield (YA) and the average date of half-campaign for cane harvest were obtained from sugarcane industries from 1998 to 2018 in the 5 agro-climatic zones presented previously. To match observations, sugarcane growth was simulated using the MOSICAS model in each field from 1998 to 2018. In each agro-climatic zone, all fields were assumed to have the same harvest date, defined as the date of half-campaign. Considering that the harvest campaign last around 6 months in Reunion Island, we considered that a sugarcane field was harvested at five different dates (half-campaign, +/-1 months and +/-2 months), and the average yield was used. As an approximation, all sugarcane fields over Reunion Island were assumed to be ratoon crop as no information on new planted cane were available on the whole period. When field were irrigated, irrigation management was assumed to be 12,000 m 3 ha -1 year -1 (fixed dose of 25 mm week -1 ), which represents common practices in Reunion, with an efficiency of 1 for drop irrigation and 0.9 for spray irrigation (i.e. means that 10% of irrigated water is lost outside the field). Five soil layers were defined in the MOSICAS model. Each layer had identical welting point, field capacity and saturation, avoiding thus differences in terms of water drainage. To take into account the difference of soil water holding capacity (Fig. 1c), the soil depth was adjusted to match this value.

Two types of simulations were performed at field scale: sugarcane yield at the thermoradiative potential (YP, not taking into account water stress) and sugarcane yield taking into account water stress (YW). The MOSICAS model does not take into account fertilizer stress or other biotic stresses such as weed pressure. All components of cane growth were simulated but only the fresh stalk yield (equivalent to millable stalk) was used in this study and defined as "cane yield". The average YP and YW in a given agro-climatic zone was calculated as the average field-scale YP and YW weighted by the field planted areas. Total Yield Gap (YGT) was defined as the difference between YP and the actual cane yield (YA), averaged at the agro-climatic zone scale. Water-limited yield gap (YGW) was defined as the difference between YP and YW. Finally, YG caused by crop management and hazards (YGM) was defined as the difference between YW and YA.

Extreme events index

Based on daily climatic data, Extreme Climatic Index (ECI) were defined at the agro-climatic scale: extreme precipitation (ECIPPT), extreme temperature (ECITMAX), extreme drought (ECIDRY) and extreme wind speed (ECIWIND) index. For precipitation, ECIPPT was defined at the field scale as the proportion of daily precipitation above a threshold (defining extreme precipitations) between two harvest dates. The precipitation threshold was assessed as the 97.5% quantile of measurements from 1998 to 2018 across space and was equal to 28.2 mm d -1 . Using the same method, ECITMAX was defined as the proportion of extreme daily event with TMAX higher than 31.7°C. Considering ECIDRY, in each field all drought events between two harvest dates were characterized by their duration in days. A threshold of 12 days for drought duration was defined using the 97.5% quantile across space and time. Finally, ECIDRY was defined as the number of drought event higher than 12 days in each field between two harvest days. Average ECIPPT, ECITMAX and ECIDRY per agro-climatic zone were then calculated as the average field scale ECI weighted by the field area. Considering ECIWIND, interpolated field scale wind speed was not available and only one station per zone was used The ECIWIND index was calculated using the daily maximum wind speed with the same method as previous ECI with a threshold of 16 m s -1 .

Statistical and data analyses

The whole dataset is accessible via Dataverse (Christina et al., 2020c). All statistical analyses were made using R 3.6 (R Development Core Team 2019). Weather data (Fig 2 .) and ECI (Fig 4 ) were compared among agro-climatic zones using Tukey HSD test (multicomp R package). Differences were considered statistically significant if p < 0.05. To improve data visualisation (Fig. 2, 3 and 5), smooth regression were performed using local polynomial regression fitting with variability bands (loess.sd function from the msir R package with a constant smoothing parameter of 0.5). Principal components analyses (PCA) were used to visually illustrate separately the effects of climate variables (Fig. 3), on one hand, and extreme climate indices (Fig. 5), on the other hand. PCA was performed with FactoMineR and factoextra R packages. Correlation between the first two axes of the PCA and YGT and YGM were tested using a spearman correlation test (cor.test function). Finally, a hierarchical clustering analysis was performed on the ECIPPT and ECIWIND components to illustrate the effect of cyclones (HCPC function in FactoMineR package). Based on this analysis, 4 groups were isolated and the difference of YGT, YGM and actual yield loss (i.e. the difference between the maximum YA and the YA of each year in a given agro-climatic zone) among these groups were tested using a Dunn multiple test comparison, as the residual normality and homoscedasticity conditions were not satisfied (dunn.test R package). In this Dunn test case, differences were considered statistically significant if p < 0.025.

Results

Sugarcane yield and yield gap variability

Actual cane yield (YA) was similar among agro-climatic zones (74 t ha -1 in average) except for the North-West zone (57 t ha -1 in average, Fig 2a). This North-West zone was characterized by a low average daily temperature (20.4°C, Fig 2b) and low annual precipitation (614 mm y -1 , Fig 2c). The East and North-East zones were the warmest (22.3°C) and wettest areas (1801 mm y -1 ). Nevertheless, actual sugarcane yield showed an annual variability. The coefficient of variation due to annual variability was similar among agro-climatic zones and varied from 8 % to 11 %. Moreover, for a given zone the difference of yield between the best year and the worst one varied between 20 t ha -1 (North West) and 31 t ha -1 (East).

In average, the simulated potential (YP) and water-limited (YW) yields at the agro-climatic zone scale were 170 (+/-19) and 106 (+/-33) t ha -1 , leading to a total yield gap (YGT, i.e. difference between YP and YA) of 96 t ha -1 (Fig 2d). This yield gap resulted from a 64 t ha -1 water limited yield gap (YGW, difference between YP and YW) and from a 32 t ha -1 yield gap due to management and hazard (YGM, i.e. difference between YW and YA), in average. While YP showed little annual variation, with a coefficient of variation ranging from 2 to 3 % depending on the agro-climatic zone, YW presented a higher variability over years, ranging from 13 % (North East) to 26 % (North West). In comparison, YA had annual variations ranging from 8 to 11 %. Total yield gap (YGT) was positively correlated with CI1 value (ρSPEAR = 0.56, p < 10 -6 ), with YGT reaching around 120 t ha -1 in warm and high ETP environments (Fig 3b). In colder environments (i.e. low values of CI1), YGT was relatively constant (around 90 t ha -1 ). YGT was also positively correlated (ρSPEAR = 0.23, p=0.021) with CI2, but it varied less from 90 t ha -1 in low precipitation environments to 110 t ha -1 in high precipitation ones (Fig 3c).

Following a similar pattern, management yield gap (YGM) was positively correlated with CI1 (ρSPEAR = 0.57, p < 10 -6 ). YGM increased from 0 to 80 t ha -1 , with increasing temperature and ETP. YGM was also strongly correlated with CI2 (ρSPEAR = 0.60, p < 10 -6 ). In low precipitation environments, YGM was null but reached 80 t ha - 1 for the years and zones with the highest precipitations.

In warmer and higher ETP environment, YGT increased mainly due the increase in YGM while YGW (difference between YGT and YGM) remained relatively constant. In higher precipitation environments, YGT also increased due to increase in YGM but YGW decreased. Total yield gap (YGT) was negatively correlated with ECI1 value (ρSPEAR = -0.56, p < 10 -6 ), with YGT reaching around 120 t ha - 1 in environments with extreme rainfall events (Fig 5b). In drier environments (i.e. high values of ECI1), YGT was relatively constant (86 t ha -1 ). YGT was not monotony correlated with ECI2 but the lowest YGT values were found in environments with medium wind speed events (Fig 5c).

Following a similar pattern, management yield gap (YGM) was negatively correlated with ECI1 (ρSPEAR = -0.8, p < 10 -6 ). YGM increased from 0 to 100 t ha -1 , with increasing extreme rainfall events. Similarly, to YGT, YGM was not correlated with ECI2 but the lowest values were found to medium wind speed environments.

Fig 5. Principal component analysis (PCA) on extreme climatic index defined in section 2.5 (a): the extreme climatic indices were calculated based on high precipitation (ECIPPT), high wind (ECIWIND), maximum temperature (ECITMAX) and drought length (ECIDRY). The change in total (YGT) and management (YGM) yield gaps with the two first axis of the PCA (ECI2) are shown in b and c. The spearman correlation between the yield gaps and the PCA axis are indicated. A polynomial smooth regression was added (line) with standard deviation (transparent). Arrows were added to help the figure interpretation.

3.4. Yield loss under extreme rainfall and wind speed Based on a hierarchical clustering analysis on ECIPPT and ECIWIND variables, four clusters of yearzone units were defined (Fig 6a). The cluster 1 included environments with low ECIPPT and ECIWIND. The cluster 2 included environments with higher ECIWIND values but with low ECIPPT. The clusters 3 and 4 both included environments with high ECIPPT, but only the cluster 4 also showed extreme high wind speed events. Management yield gap (YGM), total yield gap (YGT) and actual yield loss (i.e. difference between the actual yield of the best year of a given zone with the considered year) all increased from group 1 to 4 (Fig 6). In average, YGT, YGM and actual yield loss were equal to 92, 16 and 12 t ha -1 in cluster 1 with low rainfall and wind extreme events. In environments with high rainfall and wind extreme events, YGT increased by 11 and 22 %, YGM by 265 and 393%, and actual yield loss increased by 58 and 85%, in clusters 3 and 4 respectively. In years with high precipitation and wind extreme events (clusters 3 and 4), the loss in actual cane yield was 19 t ha -1 in average, compared to the year with highest yield, and could reach 28 t ha -1 .

Fig. 6. Clustering analysis on extreme indices for precipitation (ECIPPT) and wind speed (ECIWIND) variables (a)

: four clusters were defined. Mean management yield gap, total yield gap and actual yield loss depending on these clusters are presented in b,c,d. The actual yield loss is defined as the difference between the YA of the best year and the YA of each year in a given agro-climatic zone. Different letters indicated significant difference using a dunn test analysis.

Discussion

Yield gap in Reunion Island

The assessed total yield gap (i.e., the difference between thermoradiative potential cane yield and actual cane yield) over Reunion Island in our study was 96 t ha -1 on average (ranging from 67 to 132.6 t ha -1 depending on the years and zones). In this study, all fields were assumed to be ratoon crop which could slightly increase the yield gap at the regional scale. Nevertheless, this yield gap was weaker than the yield gap assessed in other countries. As, an example, Dias and Sentelhas ( 2018) assessed an average yield gap of 134 t ha -1 over Brazil. Similarly, Zu et al. ( 2018) assessed a regional average yield gap in China of 156.3 t ha -1 , varying from 74.82 t ha -1 to 214.2 t ha -1 , depending on climate conditions and crop management practices. Our results showed that the proportion of the yield gap due to water deficit and crop management was strongly variable. In average, water deficit account for 66% of the total yield gap in Reunion Island, which was similar with has been found in other climates (e.g. 73% in Brazil, Dias and Senthelos). As shown in this study, this proportion differed spatially and annually, with 88% in average in the dry west coast and 50% in the east wet coast (Fig S5). Part of this spatial variability could be explained by irrigation difference between east and west, but YGW was only ~25 t ha -1 higher in rainfed fields that in irrigated ones in average whatever the precipitation (Fig S6). As a comparison, in the neighbour Mauritius Island, previous studies have shown more than 100 t ha -1 yield differences between irrigated and rainfed modeled water-limited yield (Cheeroo-Nayamuth et al. 2011). Another part of this variability could be explained by soil water holding capacity (SWC), in particular because irrigation appeared to be inefficient for SWC lower than 75 mm (Fig. S7). Nevertheless, while at the field scale SWC variability strongly influenced YGW (Fig S7 ), average sugarcane field SWC per agroclimatic zone was only varying from 85 to 104 mm in our study. Consequently, the remaining factor explaining YGW differences in the agro-climatic zones in Reunion Island could be caused by factors not taken into account in Mosicas model such as increasing biotic and abiotic pressure on yields link to climate (weed pressure, management practices…).

Increasing yield gap with source availability

Our study highlighted how yield gap increase with thermoradiative sources (TMEAN, RG, ETP) and water source (precipitation + irrigation, PPTIRR) availability in Reunion Island. Such observation means that sugarcane growth cannot reach its full potential even with increasing radiation use potential. This observation suggested than other limiting factors for sugarcane growth occurred, such as fertilization, weed pressure or water losses.

In Reunion Island, limiting fertilization is not likely to be a main factor explaining increasing YG with source availability, because growers received fertilization advice based on subsidized soil analysis in Reunion Islands (Serdaf program, Versini et al., 2018). Nevertheless, this factor has to be considered as inadequate nutrition has been highlighted as a major impact explaining yield gap in other countries such as China (Zu et al., 2018) or South Africa (Van de Berg and Singels, 2013) with similar fertilization rates. Rather than fertilizers amount, a fertilizer use efficiency decreasing with higher rain (due to nitrate lixiviation, volatilization or surface erosion for example) could explain the increasing YG. As an example, Poultney et al. ( 2019) assessed the fertilizer use efficiency as only 22% during a rainy year in the North-East coast of Reunion.

Increased thermoradiative and water source availability is also likely to increase weeds pressure on sugarcane growth (Patterson et al., 1999;Ibrahim 1984), which is a main constraint in most tropical climate (Oerke 2006(Oerke , 2004)). Temperature and precipitation are the primary variables controlling the distribution of vegetation (Woodward and Williams, 1987). Weeds share the same trophic level with crops and compete with crops for the resources. Additionally, growth of weeds of warm season is generally very responsive to small increases in temperature (Flint et al., 1984;Flint and Patterson, 1983;Patterson, 1993). Consequently, increase in weed pressure could be a factor explaining the increase in yield gap observed in Reunion Islands under elevated temperature and precipitation.

Finally, the water use efficiency under higher precipitation must be considered due to the irregular topography of agricultural basin, particularly in high altitude. The model could underestimate the water loss due to drainage or runoff in these steep reliefs. Nevertheless, we argue that this effect in high altitude on water-limited yield should be limited because the fields in highest altitudes already showed the highest proportion of precipitation loss into drainage in our simulation (Fig S8), probably due to the lower ETP and water consumption by canes in altitude. Consequently, we hypothesize that the increasing yield gap with increasing thermoradiative and water sources would mainly be due to decreasing fertilizer use efficiency and increasing weed pressure.

Managing sugarcane under extreme events

This study suggests that storms and cyclones are main factors affecting the yield gap and inducing yield loss in Reunion Island. Total yield gap increased by 20 t ha -1 when high extreme precipitation and wind speed events occurred, corresponding to a similar increase in actual yield loss. Issues with storms is not limited to this part of the Indian Ocean, other productive countries such as Philippines with extreme winds (Stromberg et al., 2011) or USA with flood (Deren et al., 1991) are sensitive to these extreme events which can cause considerable economic damage. In Reunion island drought and warm event intensity is relatively low, inducing weak effect on yield gap in this study. Nevertheless, other countries have been shown to be much more sensible to extreme temperature (e.g. India, Mall et al., 2016) or drought (e.g. Australia, Inman-Bamber et al., 2005).

The effect of storms on yield is variable and it could directly influence crop physiology or impact crop yield through other biotic and abiotic pressure. As an example, Gilbert et al. 2008 showed that long-term flooding in Florida causes prolonged anoxia in the root zone leading to reductions of 10-78% in sugarcane leaf nutrient concentration, as well as leaf area index (LAI) decrease. Additionally, Beuzelin et al. (2009) highlighted how the increase in pest severity after a storm in Louisiana was a crucial factor explaining low yield values. Extreme winds can also directly impact crop losses through stalk growth termination (Moore et al., 1984).

Whatever the cause (cane physiology or breakage), opening the canopy can have drastic consequences when multiple storms can occur, such as in Reunion Island. In particular, these conditions promote soil erosion and weed development. The pouring of the canes allows light to reach the ground which favours the emergence of weeds and the recovery of plants kept under cover, such as vines (climbing weed or vine weed). These infested areas are therefore difficult to weed because they are not easily spotted and are inaccessible due to the mass of vegetation.

Considering that climate change is likely to increase extreme events occurrence (Cai et al., 2015), there is need to adapt cultural practice to reduce sugarcane sensitivity to storm. A first possibility already highlighted in 1940s by Hebert and Arcenaux (1940) is the use of wind-resistant clones. In their study in Louisiana, field survey showed large difference of broken cane (varying from 0.5 to 21% of total cane) depending of cultivars. Secondly, management practices should be promoted to limit soil erosion and weed pressure after storms. For example, the use of cover crops to limit erosions (Perret et al., 1996, Normand 1998) or improve weed control (Christina et al., pre-print) has been investigated in Reunion island and other sugarcane countries, but it could only be used at the beginning of the sugarcane crop cycle (three months before canopy closure). Another possibility is to promote the conservation of straw mulch in the field. In a review, conservation of straw has been shown to considerably reduce soil erosions and weed pressure, particularly when the amount of straw is superior to 8 t ha -1 (Carvalho et al., 2017). Finally, one could expect a potential increase in fertilizer leaching and runoff after storms. Not surprisingly, nitrogen loss through runoff or leaching showed high variations depending on location, climate and fertiliser rate (Thornburn et al., 2011) and tend to increase with higher rainfall (Reading et al., 2019). In comparison, the relation between extreme rainfall events and fertilizer loss remained poorly investigated, probably because famers naturally avoid fertilizing during periods of climatic risk.

Increasing extreme weather events worldwide are an immediate threat to global crop production (Meehl et al., 2000;Rosenzweig e tal., 2001;Olesen et al., 2007;Urban et al., 2012;Min et al., 2011;Lobell et al., 2008Lobell et al., , 2013)). Promoting crop management to make agrosystems more resilient to climatic hazards is a topical issue in the tropics.

Conclusion

In average, average total yield gap over Reunion Island was predicted as 96 t ha -1 , of which 66% was due to water deficit. This yield gap increased up to 120 t ha -1 with increasing thermoradiative and water source availability, highlighting the influence of others biotic and abiotic limiting factors. The use of extreme climate indices allows us to assess the impact of tropical cyclones and storm on yield gap and yield loss. In average, yield loss was around 19 t ha -1 in the zone and years with the highest precipitation and wind speed events, but this loss could reach 28 t ha -1 . Considering the considerable economic damage implied by yield loss under extreme events, there is a need to adapt cultural practice to frequent storms in tropical areas, particularly in the context of climate change. Where 𝑟𝑢𝑒𝑚𝑎𝑥 𝑝 is the maximum radiation use efficiency, 𝑠𝑤𝑑𝑓 𝑟𝑢𝑒 , the soil water deficit factor for radiation use, and 𝑘 𝑎𝑔𝑒 and 𝑘 𝑇 some parameter defining the effect of maintenance respiration and temperature, respectively.

The effect of temperature on radiation conversion is defined as:

𝑘 𝑇 = 1 -𝑟𝑢𝑒𝑡𝑘 𝑝 * √|𝑇 𝑑𝑖𝑢 -𝑟𝑢𝑒𝑜𝑝𝑡 𝑝 |
Where 𝑟𝑢𝑒𝑜𝑝𝑡 𝑝 and 𝑟𝑢𝑒𝑡𝑘 𝑝 are plant parameters that define the optimal temperature for radiation use efficiency and a parameter used in the response of this efficiency to temperature, respectively. 𝑇 𝑑𝑖𝑢 defines the average temperature during daylight:

𝑇 𝑑𝑖𝑢 = 0.25 𝑇𝑚𝑖𝑛 + 0.75 𝑇𝑚𝑎𝑥

The effect of maintenance respiration is calculating based on the total dry mass (𝑑𝑚𝑡 𝑑 )

𝑘 𝑎𝑔𝑒 = 1 -(𝑑𝑚𝑡 𝑑 * 𝑝01 𝑝 10 4 ) 3

Where 𝑝01 𝑝 is a plant parameter to take into account maintenance respiration

• Dry mass partitioning

The daily total dry mass accumulated (𝑔𝑑𝑚𝑡) is firstly allocated to aerial part vs roots. In case of plant crop, the dry mass allocated to roots (𝑘 𝑟𝑎𝑐 ) decreases with thermal time:

𝑘 𝑟𝑎𝑐 = 𝑝𝑟𝑎𝑐𝑑𝑒𝑏 𝑝 -𝑝𝑟𝑎𝑐𝑑𝑒𝑐 𝑝 * 𝑠𝑑𝑑 And 𝑘 𝑟𝑎𝑐 = min (𝑘 𝑟𝑎𝑐 , 𝑝𝑟𝑎𝑐𝑓𝑖𝑛 𝑝 )

Where 𝑝𝑟𝑎𝑐𝑑𝑒𝑏 𝑝 , 𝑝𝑟𝑎𝑐𝑑𝑒𝑐 𝑝 and𝑝𝑟𝑎𝑐𝑓𝑖𝑛 𝑝 are plant parameters (Table S1), and 𝑠𝑑𝑑 is the sum of degree days since planting (base temperature = 0°C).

In case of ratoon crop, the dry mass allocated to roots is assumed constant and equal to 𝑝𝑟𝑎𝑐𝑓𝑖𝑛 𝑝 .

Consequently, the aerial dry mass (𝑑𝑚𝑎𝑒𝑟 𝑑 ) and daily dry mass allocated to aerial parts are equal to:

𝑑𝑚𝑎𝑒𝑟 𝑑 = 𝑑𝑚𝑎𝑒𝑟 𝑑-1 + 𝑔𝑑𝑚𝑎𝑒𝑟 And 𝑔𝑑𝑚𝑎𝑒𝑟 = 𝑔𝑑𝑚𝑡 * (1 -𝑘 𝑟𝑎𝑐 )

The aboveground dry mass is partitioned into millable stalk dry mass (𝑑𝑚𝑠𝑡𝑚 𝑑 ), as follows:

𝑑𝑚𝑠𝑡𝑚 𝑑 = 𝑑𝑚𝑠𝑡𝑚 𝑑-1 + 𝑔𝑑𝑚𝑠𝑡𝑚 𝑔𝑑𝑚𝑠𝑡𝑚 = 𝑔𝑑𝑚𝑎𝑒𝑟 * 𝑝𝑠𝑡𝑒𝑚𝑓𝑖𝑛 𝑝 * (1 -𝑒 - 𝑝𝑠𝑡𝑒𝑚𝑑𝑒𝑐 𝑝 1000

(𝑑𝑚𝑎𝑒𝑟 𝑑 -𝑝𝑠𝑡𝑒𝑚𝑑𝑒𝑏 𝑝 ) )

Where 𝑝𝑠𝑡𝑒𝑚𝑓𝑖𝑛 𝑝 , 𝑝𝑠𝑡𝑒𝑚𝑑𝑒𝑐 𝑝 and 𝑝𝑠𝑡𝑒𝑚𝑑𝑒𝑏 𝑝 are plant parameters (Table S1).

The model also simulated the sucrose accumulation into the stalk as well as aboveground water content but these processes were not used into this study.

• Cane yield: fresh cane mass Finally, the fresh stalk mass (𝑓𝑚𝑠𝑡𝑚) and the stalk humidity are calculated as follows:

𝑓𝑚𝑠𝑡𝑚 = 𝑑𝑚𝑠𝑡𝑚 100 -ℎ𝑢𝑚𝑠𝑡𝑚 ℎ𝑢𝑚𝑠𝑡𝑚 = ℎ𝑢𝑚𝑠𝑡𝑒𝑚𝑑𝑒𝑏 𝑝 -ℎ𝑢𝑚𝑠𝑡𝑒𝑚𝑑𝑒𝑐 𝑝 * 𝑑𝑎𝑠 Where 𝑑𝑎𝑠 is the number of day since stalk appearance, and ℎ𝑢𝑚𝑠𝑡𝑒𝑚𝑑𝑒𝑏 𝑝 and ℎ𝑢𝑚𝑠𝑡𝑒𝑚𝑑𝑒𝑐 𝑝 are plant parameters (Table S1).

The stalk appearance is defined as the day when the aboveground dry mass reached a value higher than 𝑝𝑠𝑡𝑒𝑚𝑑𝑒𝑏 𝑝 .

• Water stress index Simulation of the soil water balance relies on a CERES type water module, including soil evaporation, soil drainage, plant evaporation and rainfall infiltration. The detailed of these equations can be found in the source code of MOSICAS (https://github.com/ChristinaMathias/MOSICAS). Two soil water deficit factors are calculated for 𝐿𝐴𝐼 𝑑 growth (𝑠𝑤𝑑𝑓 𝐿𝐴𝐼 ) and radiation conversion into dry mass (𝑠𝑤𝑑𝑓 𝐷𝑀 ): Where 𝑝𝑠𝑖𝑙𝑏 is a pression parameter defined as 8 according to Brisson 1992 and Slabber 1982, 𝑡𝑚𝑝 the maximum transpiration of the day according to actual canopy, sthydbio p and sthydcroip are plant parameters defining the sensitivity to water stress.

𝑠𝑤𝑑 is the soil water deficit and calculated as follows:

𝑠𝑤𝑑 = 𝑆𝑊𝐶 𝐹𝐶 -𝑊𝑃

Where 𝑆𝑊𝐶 is the actual soil water content, 𝐹𝐶 the field capacity and 𝑊𝑃 the welting point.

Fig 1 .

 1 Fig 1. Soil, climate and cane field data over Reunion Island. Annual precipitation (a) and altitude range (b) were used to define climatic areas. Soil water holding capacity over Reunion Island (c) range from 10 to 200 mm in agricultural areas. Each cane fields over the island were identified and associated with an agro-climatic zone (d). Sources: a: Météo-France & Cirad ; b: IGN & SHOM ; c: Raunet 1988 ; d: ASP 2016 from IGN ; production: Cirad AIDA, 2019; EPSG:32740 -WGS84 / UTM 40Sud.

Fig 2 .

 2 Fig 2. Average actual cane yield (a), daily temperature (b) and annual precipitation (PPT, c) over the 1998-2018 period depending on the agro-climatic zone define in Fig 1d. The annual evolution of cane yield (YA, black), simulated water-limited yield (YW, red) and potential yield (YP, blue) is presented in (d) and averaged in (e). A polynomial smooth regression over years was added (line) with standard deviation (transparent color).

3. 2 .

 2 Climate variability and yield gap Based on a principal component analysis (PCA) we summarized 5 soil and climatic variables into 2 climatic indexes (CI1 and CI2, Fig 3a). The first component (CI1, 67.9% of the explained variance) was mainly determined via the average daily temperature (TMEAN), global radiation (RG), potential evapotranpiration (ETP) and soil water holding capacity (SWC). The annual water input (precipitation + irrigation, PPTIRR) contributed to both IC1 and the second component (CI2, 18.1% of the explained variance).

Fig 3 .

 3 Fig 3. Principal component analysis (PCA) on average soil and climatic variables per agro-climatic zones (a): Average daily temperature (TMEAN), global radiation (RG), potential evapotranspiration (ETP), annual water input (precipitation + irrigation, PPTIRR) and soil water holding capacity (SWC). The change in total (YGT) and management (YGM) yield gaps with the two first axis of the PCA (CI1 and CI2) are shown in b and c. The spearman correlation between the yield gaps and the PCA axis are indicated. A polynomial smooth regression was added (line) with standard deviation (transparent). Arrows were added to help the figure interpretation.

3. 3 .

 3 Extreme climate events and yield gap variability Extreme events occurrence varied among agro-climatic zones (Fig 4). Extreme rainfall (ECIPPT) occurred more in the East and South coasts of Reunion Island than in the West Coast (Fig 4a). The North-East and South coasts underwent frequent extreme wind speed events (ECIWIND) than the rest of the island (Fig 4b). Extreme drought events (ECIDRY) showed the opposite pattern than ECIPPT, with the East Coast with the lowest values and the West Coast with the highest values (Fig 4c). Finally, considering, extreme daily temperature (ECITMAX), the East Coast was little impacted, contrary to the South and North West zones, and even more to the South-West agro-climatic zone (Fig 4d).

Fig. 4 .

 4 Fig. 4. Bean plot density of extreme events indices depending on the agro-climatic zones: extreme rain (ECIPPT, a), wind speed (ECIWIND, b), drought (ECIDRY, c) and temperature (ECITMAX, d). A Tukey test was performed and significant differences (p<0.05) among zones are indicated by different letters a,b,c.

Fig

  Fig S2. Location of the weather station network used to interpolate daily weather data. The stations with a red circle were used to calculate the extreme wind speed index. The black cross represent the center of each climatic zone, where weather data were interpolated.

Fig. S4 .

 S4 Fig. S4. MOSICAS calibration of 4 cane varieties using the ECOFI database. The model was calibrated using leaf area index, sugar yield, aboveground dry and fresh mass, stalk dry and fresh mass and aboveground and stalk water content.

Fig

  Fig S7. Field scale yield gap due to water deficit (YGW) in irrigated and rainfed fields depending on soil water holding capacity.

Fig

  Fig S8. Change in proportion of precipitation loss (drainage / precipitation) at the field scale depending on altitude (black point). A polynomial smooth regression over altitude was added (line) with standard deviation (transparent blue).

Table S1 .

 S1 Parameter values and description used in the MOSICAS model

				Cane variety
	Parameter	description	R570	R579	R585	R586	Source
	humaerdeb	Initial above ground biomass water content (%)	85.51	85.51	87.17	84.41	Calibrated
	humaerdec	Rate of decrease in above ground biomass water content	0.0058	0.0058	0.0054	0.0044	Calibrated
	humaerfin	begining in decrease of above ground biomass water content	914.39	914.39	823.97	1024.19 Calibrated
	humaertb	Base temperature for above ground biomass water content	13.70	13.70	11.88	11.43	Calibrated
	humstemdeb Initial stem water content (%)	88	88	88	88	fixed
	humstemdec	Rate of decrease in stalk water content	0.0704	0.0704	0.0704	0.0704	fixed
	ke	extinction coefficient	0.482	0.482	0.482	0.482	Calibrated
	laicroi	growth rate of lai	0.0040	0.0042	0.0042	0.0042	Calibrated
	laimax	Maximum LAI	7	7	7	7	fixed
	laitb	base temperature for lai calculation	12	11.2	12	12	Calibrated
	laiwksen	sensitivity of lai to water stress	0.0100	0.0110	0.0110	0.0110	Calibrated
	p01	coefficient for maintenance effect on conversion.	1.0665	0.0968	1.1998	0.5001	Calibrated
	p04	coefficient for diffused radiation effect on conversion.	0	0	0	0	fixed
		Extinction coefficient of daily
	pstrudec	fraction of millable stalk dry mas	1.2162	7.9221	3.8088	1.3913	Calibrated
		allocated to structures
	pstrufin	Final daily fraction of stalk dry mass allocated to structures	0.6914	0.5842	0.5135	0.6209	Calibrated
		Temperature threshold from which
	pstrutb	fraction of stalks dry mass allocated	25.37	18.39	22.52	23.67	Calibrated
		to structures is decreasing.
		Temperature effect on daily fraction
	pstrutcroi	of millable stalks dry mas allocated	0.0081	0.0181	0.0062	0.0104	Calibrated
		to structures
	ptigdeb	Beginning of millable stalk dry mass appearance	499.68	579.01	199.42	292.23	Calibrated
		Extinction coefficient of daily
	ptigdec	fraction of aboveground dry mass	1.8846	2.0213	1.2420	1.6556	Calibrated
		allocated to millable stalk
	ptigfin	Final daily fraction of aboveground dry mass allocated to millable stalk	0.9606	0.8906	1.1930	0.7953	Calibrated
		conversion coefficient of
	ruemax	intercepted photosynthetic radiation	2.6154	3.2381	3.1907	3.1049	Calibrated
		into total dry mass
	ruetopt	Optimum temperature for conversion.	30.4	30.4	30.4	30.4	Calibrated
	ruetk	Effect of temperature on conversion.	0.00037	0.00037	0.00037	0.00037 Calibrated
	sthydbio	Sensitivity to water stress index (swdf1) for mass accumulation.	0.6	0.6	0.6	0.6	Calibrated
	sthydcroi	Sensitivity to water stress index (swdf2) for growth. Plant parameter	0.500	0.514	0.514	0.514	Calibrated
	pracfin	Proportion of dry mass allocated to roots in ratoon	0.1	0.1	0.1	0.1	fixed

Acknowledgements

This study was funded by the FEADER (Fonds européen agricole pour le développement rural) and the FEDER 2014-2020 program (Fonds européen de développement régional). We are thankful to the Direction de l'Alimentation, de l'Agriculture et de la Forêt (DAAF) and the Centre Technique Interprofessionnel de la Canne et du Sucre (CTICS) for sharing the data.

Supplementary Material S1 -Cane management data in Reunion Island

Supplementary Material S3 -Description of cane growth main equations in MOSICAS model

• LAI calculation

In the daily time step mosicas model, the leaf area index (𝐿𝐴𝐼 𝑑 ) for a given day 𝑑 is a variable that grows as a function of temperature, water stress and intra-specific competition. The daily growth in lai (𝑔𝑙𝑎𝑖) is computed as follow:

Where 𝑙𝑎𝑖𝑐𝑟𝑜𝑖 𝑝 is a plant parameter, 𝐿𝐴𝐼 𝑑-1 , the 𝐿𝐴𝐼 of the previous day, 𝑙𝑎𝑖𝑚𝑎𝑥 𝑝 the maximum leaf area and is a plant parameter, 𝑠𝑤𝑑𝑓 𝐿𝐴𝐼 the soil water deficit factor for 𝐿𝐴𝐼 𝑑 . 𝑇 𝐿𝐴𝐼 is the accumulated thermal energy during day and is defined as:

Where 𝑇𝑚𝑎𝑥 𝑑 and 𝑇𝑚𝑖𝑛 𝑑 are the maximum and minimum daily temperature, and 𝑙𝑎𝑖𝑡𝑏 𝑝 is a plant parameter and represents the base temperature for growth in 𝐿𝐴𝐼 𝑑 .

The leaf mortality due to water stress (𝑠𝑙𝑎𝑖) is calculated as follows:

Where laiwksenp is a plant parameter and represents a sensitivity to water stress.

Finally, 𝐿𝐴𝐼 𝑑 is calculated as:

𝐿𝐴𝐼 𝑑 = 𝐿𝐴𝐼 𝑑-1 + 𝑔𝑙𝑎𝑖 + 𝑠𝑙𝑎𝑖

• Intercepted radiation

The daily photosynthetic active radiation intercepted by the canopy (𝑝𝑎𝑟𝑖) is calculated based on the bear-lamber law:

Where 𝑅 𝑔 is the daily global incident radiation and 𝑘𝑒 𝑝 , a plant parameter defining the extinction coefficient, 𝑝04 𝑝 , a plant parameter that takes into account the diffuse radiation effect and 𝑅 𝑒 , the global extraterrestrial radiation.

• Radiation conversion into dry mass

The total dry mass (𝑑𝑚𝑡 𝑑 ) at a given day 𝑑 is calculated as the sum of the total dry mass of the previous day and the daily growth in dry mass (𝑔𝑑𝑚𝑡). The daily 𝑔𝑑𝑚𝑡 is calculated as: 𝑔𝑑𝑚𝑡 = 𝑟𝑢𝑒𝑚𝑎𝑥 𝑝 * 𝑝𝑎𝑟𝑖 * 𝑠𝑤𝑑𝑓 𝑟𝑢𝑒 * 𝑘 𝑎𝑔𝑒 * 𝑘 𝑇

And:

𝑑𝑚𝑡 𝑑 = 𝑑𝑚𝑡 𝑑-1 + 𝑔𝑑𝑚𝑡

And 𝑡𝑚𝑝 the maximum transpiration of the day is calculated as follows:

Where 𝐸𝑃𝑇 is the daily potential evapotranspiration (mm) and Kmax p the maximum crop coefficient (plant parameter).

Supplementary Material S4 -Complementary results

Fig


ON THE ROLE OF DIFFERENTIAL ALGEBRA IN BIOLOGICAL MODELING

Keywords: AMS 2000 Subject Classification: 12H05, 37N25, 62P10 differential algebra, cellular biology

FRANC ¸OIS BOULIER

Differential algebra is an algebraic theory for studying systems of polynomial ordinary differential equations (ODE). Among all the methods developed for system modeling in cellular biology, it is particularly related to the well-established approach based on nonlinear ODE. A subtheory of the differential algebra, the differential elimination, has proved to be useful in the parameters estimation problem. It seems however still more promising in the quasi-steady state approximation theory, recent results show.

Background

Differential algebra [START_REF] Robert | Differential Algebra and Algebraic Groups[END_REF][START_REF] Fels | Differential Algebra[END_REF] is an algebraic theory for studying systems of polynomial ordinary differential equations (ODE). Among all the methods developed for system modeling in cellular biology, it is particularly related to the well-established approach based on nonlinear differential equations [START_REF] Conrad | Modeling Molecular Interaction Networks with Nonlinear Differential Equations[END_REF][START_REF] Leloup | Modeling the molecular regulatory mechanism of circadian rhythms in Drosophila[END_REF]. A very classical way for deriving a system of polynomial ODE from a genetic circuit consists in first, translating the circuit as a system of generalized chemical reactions, second translating the chemical system as a system of polynomial ODE by means of the mass-action law.

The parametric ODE systems derived from these generalized chemical reaction systems are often too complicated for further analysis and need to be reduced. They are often overparameterized, which makes fitting methods to determine the parameters values heavy to carry out and unreliable: the parameters values which reproduce a behaviour of interest are usually far from unique. See the impressive analysis of [START_REF] Dassow | Exploring modularity with dynamical models of gene networks[END_REF] about their own outstanding work on the Segment Polarity Network. Moreover, the large number of parameters and variables makes their qualitative analysis difficult to achieve: many important dynamical properties of polynomial ODE systems do not depend on precise values of the parameters [START_REF] De | Qualitative Approaches to the Analysis of Genetic Regulatory Networks[END_REF]. Discussing the presence of these properties in terms of the parameters values becomes dramatically difficult as the number of parameters and variables increases.

For the important purpose of understanding which parts of the system contribute the most to some property of interest [START_REF] Kell | The Role of Modeling in Systems Biology[END_REF], it is preferable to deal with smaller systems. More conceptually, reducing the size of a model is one way to achieve the important goal: a model should not involve more information than the biological data it comes from.

Applications to the model reduction problem

The model reduction problem is well-known in biological modeling. Though not formulated in these terms, all the approaches listed in [START_REF] De | Qualitative Approaches to the Analysis of Genetic Regulatory Networks[END_REF] address this issue. Even in the particular case of polynomial ODE models derived from generalized chemical reactions systems, many methods [START_REF] Okino | Simplification of Mathematical Models of Chemical Reaction Systems[END_REF] exist (lumpings, sensitivity analysis, multiple time-scales analysis, . . . ).

Among all these ones, the quasi-steady state approximation (QSSA) theory relies on the assumption that some of the chemical reactions are much faster than the other ones. The idea of quasi-steady state approximation is simple: study the dynamics of the slow reactions, assuming that the fast ones are at quasiequilibrium, thereby removing from the ODE system, the differential equations which describe the evolution of the variables at quasi-equilibrium. Quasi-steady state approximation has two advantages:

1. it reduces the number of ODE occuring in the system under study ; 2. it transforms stiff ODE systems as nonstiff ones.

It was extensively studied [START_REF] Okino | Simplification of Mathematical Models of Chemical Reaction Systems[END_REF][START_REF] Van | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF][START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF] and references therein. All these approaches turn out to be equivalent. All amount to a two time-scales analysis. None of them is formulated in algorithmic terms.

The differential algebra [START_REF] Robert | Differential Algebra and Algebraic Groups[END_REF][START_REF] Fels | Differential Algebra[END_REF] theory makes the quasi-steady state approximation algorithmic, in the particular setting of ODE systems derived from generalized chemical reactions systems, as shown for the first time by [START_REF] Franc ¸ois Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF]. Observe that, for general systems of ODE, the QSSA is not algorithmic. The method over which this new algorithm relies is the differential elimination theory [START_REF] Franc ¸ois Boulier | Representation for the radical of a finitely generated differential ideal[END_REF][START_REF] Franc ¸ois Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF][START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF][START_REF] Seidenberg | An elimination theory for differential algebra[END_REF], pioneered by Prof. G. Carrà-Ferro [START_REF] Carra-Ferro | Gröbner bases and differential ideals[END_REF]. After the quasi-steady state approximation step, the reduced model is still overparameterized (it is a raw reduced model). The differential algebra theory then permits to reduce the set of parameters of the raw reduced model, leading to a reduced model. This exact simplification can be algorithmically performed by determining Lie symmetries of the system [START_REF] Sedoglavic | Reduction of Algebraic Parametric Systems by Rectification of their Affine Expanded Lie Symmetries[END_REF].

The sketched above method (model reduction by QSSA followed by reduction of the parameter set) was recently applied over a family of models [START_REF] Franc ¸ois Boulier | Applying a rigorous quasi-steady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF][START_REF] Franc ¸ois Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF], giving conditions for the presence of oscillations.

Applications to the optimisation theory

Differential elimination helps solving the following problem [START_REF] Franc | Differential Elimination and Biological Modelling[END_REF][START_REF] Franc ¸ois Boulier | proceedings of the ICPSS conference[END_REF][START_REF] Noiret | Utilisation du calcul formel pour l'identifiabilité de modèles paramétriques et nouveaux algorithmes en estimation de paramètres[END_REF]: given a parametric ODE system and experimental data for some of the variables, how to estimate the parameters values ?

The idea consists in eliminating the variables for which no experimental data are available (non observed variables) then in applying linear least squares to compute a first estimate of the observed variables. This first estimate can be used as a starting point for a nonlinear least squares methods such as the Levenberg-Marquardt algorithm [START_REF] Gill | The Levenberg-Marquardt Method[END_REF]. The first estimate is usually not very precise because it requires the numerical computation of derivatives of the observed variables from the experimental data.

Put in a nutshell, differential elimination transforms nonlinear least squares problems into linear ones. It is not practical in cellular biology because it requires too accurate experimental data. However, whenever it applies, it improves the classical Newton-like methods by guessing a starting point.

In the particular case of generalized chemical reactions systems, another improvement can be imagined using the algorithm obtained in the quasi-steady state approximation theory. In the nonlinear fitting methods, the parameters are the variables. The sets of large and small parameters have to change during the process. Because of these changes of parameters, the ODE systems to numerically integrate, in order to compute the error, often get stiff at runtime: stiffness is often caused by the presence of different time scales in the ODE systems. Stiffness considerably slows down the whole optimisation problem. The algorithm for performing the quasi-steady state approximation could then be applied in order to replace at runtime the stiff ODE systems by nonstiff ones. The error computation of the nonlinear fitting methods could then be performed over the reduced system, speeding up the overall process.

Conclusion

Differential algebra has indeed a role to play for system modeling in cellular biology. However, this role is local: differential algebra tools must be used in connection with many other tools. This motivates the development of standalone, easy to plug-in software. The BLAD libraries [START_REF] Franc | The blad libraries[END_REF], developed by the author, are an attempt in that direction.