G Fubiani

) R Agnello

I Furno

L Garrigues

Ph Guittienne

G Hagelaar

A Howling

R Jacquier

B Reman

A Simonin

F Taccogna

Negative hydrogen ion dynamics inside the plasma volume of a linear device: estimates from Particle-In-Cell calculations

Negative hydrogen or deuterium ions are the precursor particles used to generate a high power beam of neutrals in order to heat the tokamak plasma core of magnetic fusion devices, inject current and to some extent control instabilities. In the case of ITER for instance, the negative ions are produced inside a high power large volume low pressure tandem type magnetised ion source and extracted toward an electrostatic accelerator which accelerates them to 1 MeV before entering a neutraliser converting the ions into a neutral beam. This so-called Neutral Beam Injector (NBI) relies on the production of negative ions on the surface facing the plasma of the ion source extraction electrode. The latter is covered by a cesium layer in order to increase the negative ion yield. The use of cesium is currently an issue as it may diffuse outside of the source and induce secondary particle production or voltage breakdowns inside the accelerator vessel requiring a regular maintenance in a nuclear environment. In this work, we analyse numerically with a 2.5D Particle-In-Cell model the production rate and transport of negative ions in a linear device used as an ion source. The negative ions are generated via a dissociative attachement process with a hydrogen molecule in the volume of a magnetised cesium-free plasma. The linear device in the model has a large aspect ratio with a radius of 5 cm and a length of 100 cm and the magnetic field strength ranges from 100 to 400G. We show that the shape and depth of the plasma potential profile may be controlled by biasing the end-plates which in turn strongly influence the residence time of the electrons and hence the negative ion yield. We observe the formation of large scale rotating structures when the positive ions become magnetised with a rotation velocity in the kHz range. TABLE I. Electron collisions. # Reaction Cross section ref. 1 e + H → e + H (elastic) 20-24 2 e + H → e + H (inelastic, 4 proc.) 25 3 e + H → 2e + H + 25 4 e + H2 → e + H2 (elastic) 26 5 e + H2 → 2e + H + 2 25 6 e + H2 → 2e + H + + H (2 proc.) 25 7 e + H2 → e + H2 (inelastic, 16 proc.) 25,27-32 8 e + H2 → e + 2H (3 proc.) 25,33 9 e + H + 3 → 3H 25 10 e + H + 3 → H + H2 25 11 e + H + 3 → e + H + + 2H 25 12 e + H + 3 → e + H + + H2 25 13 e + H + 2 → 2H 25 14 e + H + 2 → e + H + + H (2 proc.) 25,33 15 e + H + 2 → 2e + 2H + 33 16 e + H -→ 2e + H 25 17 e + H * 2 → H -+ H 33 18 e + H + 2 → e + H + 2 (Coulomb) 34 19 e + H + → e + H + (Coulomb) 34 20 e + H + 3 → e + H +

I. INTRODUCTION

Negative hydrogen isotope ions may be produced either by the dissociative attachment of a vibrationally excited hydrogen molecule or via the tunnelling process of electrons from a metal conduction band toward a neutral atom or positive ion impacting its surface. The former is a two step process where one region of the plasma source has typically a higher electron temperature of about 5-10 eV in order to excite vibrationally the hydrogen molecule by electron impacts, e + H 2 (ν i ) → e + H 2 (ν j ) where j > i and ν is the vibrational quantum state of the molecule, while in a second area the electron temperature must drop to ∼ 1 eV to optimally produce negative ions by a dissociative attachment process 1,2 , e.g., e + H 2 (ν j ) → H -+ H. A low electron temperature region also lengthens the mean-free-path of a given negative ion before being destroyed by collisions with electrons or the background gas. In practice, a plasma source with a strong gradient in electron temperature is generated by implementing a magnetic filter field between the discharge (region where the external power is coupled to the plasma) and the extraction region such as to decrease the mobility of the magnetised electrons and induce a local loss of electron kinetic energy by inelastic a) Electronic mail: gwenael.fubiani@laplace.univ-tlse.fr collisions. The magnetic field is transverse to the discharge axis. Another way to efficiently produce negative ions is through the interaction of atoms or positive ions with metallic walls. In high brightness ion sources, cesium is added in order to decrease the metal surface work function and hence increase the production yield of negative ions. Cesium is one component of negative ion sources used for fusion applications 3 as well as in high energy linear particle accelerators [4][5][6] , devices used for neutron generation 7 , tandem accelerators and accelerator based mass spectrometry. [START_REF] Yoneda | proceedings of the Ninth International Conference on Accelerator Mass Spectrometry[END_REF][START_REF] Alton | [END_REF] In the case of fusiontype plasma sources, the extracted negative ion current produced by volume processes is only a small fraction of the total (< 10 -20%); most ions originating from the surface of the extracting electrode. 10 Cesium is currently an issue in Neutral Beam Injector (NBI) concepts as it typically diffuses outside of the ion source across the extraction apertures. In the accelerator vessel, cesium can induce voltage breakdowns, possibly negative ion beamlet halos or the production of secondary particles which will absorb power from the power supplies and can deposit a power density on the accelerator parts exceeding the damaging threshold for the material. 11 Cesium hence requires a regular maintenance of the injector in a nuclear environment. The estimates for ITER and DEMO provide an evaporation rate of about 10 mg/h per oven in the ion source. 3,[START_REF] Fantz | Proceedings of the 13th International Symposium on Fusion Nuclear Technology (ISFNT-13)[END_REF] As a consequence, there is currently a research effort in the low temperature plasma community to find alternatives to cesium either by studying other low work function materials [START_REF] Cartry | [END_REF] or different ion source configurations providing a plasma confinement optimal for producing negative hydrogen isotope ions through volume processes in a cesium free plasma. One such candidate is a linear device and recent experiments on the Resonant Antenna Ion Device (RAID) 14 at the Swiss Plasma Center (SPC) in Ecole Polytechnique Fédérale de Lausanne (EPFL) has shown the capability to produce a maximum negative ion density of ∼ 2 × 10 16 m -3 measured in both hydrogen and deuterium at a radius of 6 cm for a Radio-Frequency (RF) power (provided by a helicon antenna) of 5 kW. 15,16 The peak electron density is about 2 × 10 18 m -3 and the background gas pressure 0.3 Pa. The RAID testbed is a cylinder of 20 cm radius and 1.8 m length with a magnetic field strength aligned with the axis of up to 800G. The volume of the source is ∼ 0.2 m 3 and 5 kW is hence equivalent to a power density of ∼ 25 kW/m 3 (which is non-uniformly distributed). In comparison, the BATMAN one driver ITER-prototype ion source generates a line-of-sight (LOS) averaged negative ion density of 7 × 10 16 m -3 in the presence of cesium measured by Cavity Ring Down Spectroscopy (CRDS) at 0.45 Pa in hydrogen about 2 cm from the extraction electrode (so-called "plasma grid") for a power density of ∼ 1 MW/m 3 (60 kW of RF power and a plasma volume of ∼ 50 litres). [START_REF] Wimmer | AIP Conference Proceedings[END_REF] Helicons are more efficient in terms of plasma production but the difference in power densities between the two concepts is also due to the increased plasma confinement in the case of a linear device (this in principle can favour the production of negative ions in the plasma volume). Electrons are strongly magnetized in the latter, RF power is optimally absorbed near the discharge axis (which is aligned with the magnetic field vector) and the E × B and diamagnetic drifts are closed (i.e., in the azimuthal direction). Losses consequently occur mainly on the surfaces short-circuiting the magnetic field lines. Note that the distance traveled by electrons along the field lines is typically of meter scale due to the large aspect ratio of these devices. Furthermore, in a tandem-type fusion-prototype device such as BATMAN or ELISE, the negative ions are produced on the cesium covered plasma grid (PG) surface, i.e., in the direct vicinity of the extraction apertures, and hence can be efficiently extracted toward the accelerator vessel. In a linear device, the negative ion density typically peaks near the edges of the region where the helicon antenna power is coupled to the plasma. The question of an optimal extraction of the negative ions remains an open problem at the present time; this likely requires a helicon antenna radius similar to the one of the linear device to shift the negative ion density maximum toward the cylinder side wall.

In this work, we analyse the production yield and transport dynamics of negative ions in a linear device using a 2-dimensional (2D) Particle-In-Cell model including losses in the third dimension calculated analytically (method commonly dubbed "2.5D" for that reason). We show that the residence time of the magnetised electrons can be significantly increased by biasing both end-plate electrodes with the same voltage. The latter are perpendicular to the magnetic field lines (which are straight) and hence a bias voltage may be used to control the electron losses along these lines. The residence time of the electrons is strongly correlated with the yield for producing negative ions in the plasma volume via the two step process consisting of the vibrational excitation followed by the dissociative attachement of the hydrogen isotope molecules. In addition, we observe the formation of large scale rotating structures at frequencies of about tens of kHz when the positive ions become somewhat magnetised. In the next section, we describe in detail the numerical model, the plasma chemistry which we implemented, and the simulation domain. In Sec. III, we analyse the plasma properties when the end-plates are biased negatively so as to reduce electron losses and hence increase their confinement. In Sec. IV, we compare the model predictions to the experimental measurements on the RAID device. In Sec. V, we study the impact of the magnetic field strength on the plasma dynamics and notably the appearance of instabilities. Lastly in Sec. VI, we discuss the possibility of extracting negative ions from a slit aperture embedded on the cylinder wall of the linear device. The numerical resolution does not allow us to model the details of negative ion extraction across the aperture. We analyse instead the effect on the 2D plasma potential profile of biasing positively an electrode on the wall; the aim being to assess whether or not one can channel negative ions from inside the plasma volume toward the aperture.

II. NUMERICAL MODEL

The numerical model is an explicit electrostatic Particle-In-Cell (PIC) algorithm with Monte-Carlo Collisions (MCC). Maxwell's equations are solved in the plasma volume within the electrostatic limit with the plasma potential derived from Poisson's equation. The magnetic field is provided by electromagnets and its profile is prescribed in the model. The algorithm is explicit and as a consequence, the source term for Poisson's equation (the total charge density ρ) is calculated from the particle trajectories estimated at the previous time step. We assume that the magnetic field lines are straight, parallel to the linear device axis (see fig. 1). The field strength B 0 is constant. The simulation domain and hence the numerical mesh are 2-dimensional (2D) with particle losses in the third dimension calculated analytically to account for a 3D volume. [START_REF] Fubiani | [END_REF] In that direction, the rates of positive ion losses are estimated from the Bohm velocity which is the average ion velocity at the entrance of the plasma sheath. The loss frequency is calculated as follows,

ν L = 2κu B /L z , (1)
where u B = eT e (x, y)/m i is the local Bohm velocity, T e the electron temperature, L z the length of the linear device and κ = n s /n the ratio between the plasma density at the sheath edge n s and the average plasma density n along (Oz). Equation (1) was derived by volume averaging the continuity equation in the 1D limit assuming quasi-neutrality. κ 0.5 in a low pressure plasma with ambipolar diffusion (along the magnetic field lines) and nearly collisionless ions. [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] In this work the negative ion to electron density ratio is at maximum around 10% and the negative ion average kinetic energy ∼ 0.3 eV. For this reason, we considered that this ion specie has hence a minor impact on the plasma potential profile and Bohm velocity at the sheath entrance. [START_REF] Lieberman | Principles of Plasma Discharges and Materials Processing[END_REF] The loss rate for the electrons is evaluated by following their trajectories in the third direction. The latter gyrate around the magnetic field lines (electrons are strongly magnetized in our case) along (Oz) and are lost to the wall if their kinetic energies exceed the difference between the wall and plasma poten-

tials 1 2 m e v 2 z ≥ e (φ -V ep )
, where m e and e are the electron mass and charge, v z its velocity in the z-direction, V ep the end-plate bias voltage and φ(x, y) the plasma potential calculated on the 2D mesh, respectively. In 2.5D, the plasma parameters are consequently averaged values over the length L z .

Electrons absorb a given power which is an external parameter to the model. We assume that a fraction N h = N e ν h ∆t of the total number of macro-electrons N e in the power absorption region will couple the external power to the plasma during a time lag ∆t and furthermore that the absorbed energy per electron is identical on average, i.e., T h = 2 3 ∆E h with,

∆E h = P abs ∆t N m N h , (2)
where we defined an artificial heating temperature T h . N m is the macroparticle weight (the number of real electrons per macro-electron) and ν h a heating frequency (which is arbitrary), respectively. We assume that the power absorption profile is Gaussian and the macroelectrons have hence a probability

p c = exp -x 2 i + y 2 i 2σ 2 h ( 3
)
of undergoing a heating collision. x i and y i are the macro-particle coordinates and σ h is the standard deviation, respectively. Their velocity increments ∆v x , ∆v y and ∆v z are sampled from a Maxwellian distribution with a temperature T h such that,

v i = v i + ∆v i , (4)
and i = x, y, z. In practice, one loops over the N e macroelectrons inside the power absorption region. A candidate electron is selected randomly inside the particle array and the actual occurence of a collision is verified via an additional random number r between 0 and 1 such that r ≤ p c (x i , y i ). The process is repeated until N h collisions have materialized. Equation ( 4) preserves the average velocity of the electron distribution function, that is, v i = v i . On average, the electron energy is increased by 1 2 m e ∆v 2 i . Note that one can add a velocity increment along a single direction instead to mimic the power absorbed on average by an electron interacting with the RF wave generated by a time varying current in a coil. This method generates a Maxwellian Electron Energy Distribution Function (EEDF) in hydrogen (see sec. III B).

The complex physical-chemistry resulting from the interaction between electrons, ions and the hydrogen gas are included in the model. The neutral dynamics are not simulated and we implement a uniform gas density and temperature profile instead with a given background density n H2 for the hydrogen molecules and n H for the atoms corresponding to a dissociation rate n H /n H2 ∼ 0.25 (which is typically observed in ITER-like ion sources. [START_REF] Fubiani | [END_REF] ) Collisions are modelled between charged particles and neutrals via the Monte-Carlo (MC) technique where we derive a velocity for the latter from a Maxwellian distribution at a temperature T n while for the interaction between electrons and positive ions (Coulomb collisions), the actual velocities of both the incident and target macroparticles are considered during the scattering process. 35,36 The probability for the occurence of a reaction is calculated from the tabulated cross-sections and particles are scattered in their center of mass after the collision. For more details, please see ref. [START_REF] Fubiani | [END_REF] . The set of re- actions is presented in tables I and II. [START_REF] Fubiani | [END_REF]37,38 Table I corresponds to the collision processes associated with electrons. Reactions #2, 6, 7, 8 and 14 combine inelastic processes. Reaction #2 regroups the excitation of the hydrogen atom from the ground state to the electronic level n = 2-5. [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF] Reaction #7 combines the ground state excitation of the hydrogen molecule H 2 (X 1 Σ + g ; ν = 0) to the vibrational levels ν = 1 -3 [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF]31 , electronic levels (for all ν

) B 1 Σ u , B 1 Σ u , B 1 Σ u , C 1 Π u , D 1 Π u , D 1 Π u , a 3 Σ + g , c 3 Π u , d 3 Π u 25
, Rydberg states 32 and lastly rotational levels J = 2 27,28 and 3. 29,30 Reaction #17 models the generation of negative ions in the ion source volume, which are a byproduct of the dissociative impact between an electron and molecular hydrogen H 2 (ν ≥ 4). [START_REF] Janev | Collision processes in low-temperature hydrogen plasmas[END_REF] We do not calculate self-consistently the concentration of excited species in the model. To estimate the volume production of negative ions, we assume that 1.5% of H 2 molecules (except in sec. IV) are excited in vibrational levels ν ≥ 4. This is in accordance with the H 2 vibrational distribution function calculated with the 0D model of ref. 39 . Changing this parameter modulates the magnitude of the negative ion density. Table II summarizes the collision processes of heavy ions with neutrals. Reaction #9 corresponds to the excitation of the hydrogen molecule from the ground state to vibrationally excited levels ν = 1 -2 40,[START_REF] Krstić | Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms[END_REF] and to the rotational levels J = 2 -3. [START_REF] Phelps | [END_REF] To our knowledge there are no reliable data available for the elastic collision between H + 3 and neutral atoms (reaction #2), we consequently use the same cross-section as in reaction #1.

The dimensions of the grid must be of the order of the Debye length and one needs to resolve the electron plasma frequency in a PIC-MCC model otherwise

! ! ! " " #$ # FIG.
1. Schematic illustration of the simulation domain. The numerical mesh is 2D and in the (Ox), (Oy) directions with particles losses along (Oz) estimated semi-analytically via (1) the Bohm frequency for the positive ions while (2) trajectories of negative charges are followed with losses occurring if their kinetic energy in the z-direction exceeds the difference between the plasma and end-plate (Vep) potentials. In all the calculations except for sec. IV (comparison with the experiments), the radius of the linear device is R = 5 cm, length Lz = 100 cm, an absorbed power of 2.4 kW and a standard deviation for the power profile σ h = 0.92 cm, respectively. The neutral hydrogen gas species densities are fixed (uniform profiles) with nH 2 = 8 × 10 19 m -3 , nH = 2 × 10 19 m -3 together with temperatures of TH 2 = 0.1 eV and TH = 0.3 eV. Negative ions are produced by the dissociative attachment of a vibrationally excited hydrogen molecule. We assume that 1.5% of H2 molecules are in a vibrational state ν ≥ 4. Lastly, the magnetic field is uniform and parallel to the z-axis with a strength ranging from 100 to 400G. Vep is typically negative to increase the residence time of an electron along a given magnetic field line.

the algorithm would be numerically unstable. In this work, the numerical resolution for the simulations ranges between 256 2 and 1024 2 grid nodes which corresponds to plasma densities ∼ 10 14 and 10 15 m -3 , respectively. Plasma densities on axis in the RAID experiments or in the ITER prototype ion sources are typically around ∼ 10 18 m -3 and hence about 3 orders of magnitude larger than in the simulations (or similarly an electron Debye length 30 times smaller). We are interested in studying particle dynamics in quiescent plasmas and as a consequence the plasma behaviour in the quasi-neutral volume will be similar regardless of the plasma density assuming that the particle collision mean-free-paths (mfp) stay constant. The principal difference is the length of the plasma sheath but its influence remains limited as long as its dynamics are preserved and the quasi-neutral volume largely exceeds the area occupied by the sheath. In this work, the sheath is collisionless and non-magnetised in both the model and the experiments for a magnetic field strength of up to ∼ 100G. For larger magnetic fields, the question of the magnetisation of the electrons in the simulated sheath is an open problem (because we model lower plasma densities). Note that it is considered essential not to strongly magnetise the negative ions in order to facilitate their sideways extraction from the plasma through a slit aperture on the linear device cylinder wall and also due to the beam transport properties in the accelerator which will be deteriorated in case of magnetisation. The mfp for the collisions between charged particles and neutrals is independent of the simulated plasma density but not for the case of Coulomb interactions between electrons and ions. One hence needs to multiply the Rutherford scattering cross-sections by a scaling factor to correctly calculate the particle mfp, which is simply the ratio of the modelled plasma density divided by the target one. Another simpler and equivalent approach is to replace the vacuum permittivity constant ε 0 instead by ε * (i.e., a scaling by a factor α = ε * /ε 0 ) in Poisson's equation (in this case the plasma density in the model is equal to the experimental value one wishes to simulate),

-∇ 2 φ = ρ/ε * (5)
where ε * is the permittivity in the model which only affects the length of the plasma sheath (artificially increased by √ α) in a quiescent plasma. Hereafter, we chose the latter approach. The inclusion of the electronelectron scattering process is left for future work. Note that it does not induce any transport across the magnetic field lines even in the presence of a density gradient. [START_REF] Golant | Fundamentals of plasma physics[END_REF][START_REF] Nicholson | Introduction to plasma theory[END_REF] Figure 1 shows the simulation domain in the model. The linear device has a length L z = 1 m and a radius R = 5 cm. The numerical mesh is 2D Cartesian in the (Ox) and (Oy) directions while particle losses along the third dimension -(Oz) -are calculated semianalytically as stated above. The area in light gray corresponds to the plasma location and the one in dark grey to the power absorption region using the profile defined by eq. ( 3). The magnetic field lines are straight, in the z-direction, and are terminated on the end-plates which are biased by a voltage V ep . The magnetic field strength in the calculations analysed in this work ranges from 100 to 400G resulting in a situation where ions are either not or somewhat magnetised. Electrons are in all cases strongly magnetized with, for instance, a Hall parameter h e = Ω ce /ν e ∼ 40 for |B z | = 100G where Ω ce = eB z /m e is the electron cyclotron frequency and ν e the total collision counterpart (which includes elastic, inelastic and ionisation). Electrons are hence well confined by the magnetic field and their residence time can be further increased by biasing negatively the end-plate potential V ep . Ions are magnetised for |B z | = 400G with a corresponding Hall parameter h i ∼ 5.

III. INCIDENCE OF THE END-PLATE POTENTIAL ON THE PLASMA PROPERTIES

In this section, we analyse the impact on the plasma properties of biasing negatively the metallic end-plates. We start by describing in detail a configuration where the end-plates are floating before performing a parameter scan. Note that we considered smaller dimensions than in the RAID experiment as it allows us to implement a higher numerical resolution.

A. Floating end-plates

Figure 2 shows the electron density, temperature, plasma potential and negative ion density 2D profiles for a plasma column with floating end-plates. The side wall of the cylinder is grounded. The numerical resolution of the calculation is 256 2 grid nodes, 40 particles-per-cell (ppc), an average plasma density n p 5.5 × 10 17 m -3 , an absorbed power of 2.4 kW, a standard deviation for the power density profile of σ h = 0.92 cm, a heating frequency ν h = 10 7 s -1 , a magnetic field strength of 100G and a scaling factor α = 3200 (corresponding to an electron Debye length to grid size ratio ∆x i /λ De 0.5), respectively. The background gas specie is hydrogen with a density of n H2 = 8×10 19 m -3 , temperature T H2 = 0.1 eV for the molecules, n H = 2 × 10 19 m -3 and a temperature T H = 0.3 eV for the atoms. These temperatures are estimates based on experimental measurements and PIC-MCC calculations of fusion-type RF powered tandem type magnetised negative ion sources. [START_REF] Fubiani | [END_REF] Hydrogen atoms were found to have typically a larger temperature than the molecules due (1) to the Franck-Condon mech- anism which produces fast atoms as a byproduct of a dissociative collision between H 2 and electrons 33 and (2) in low background gas pressure working conditions as a result of the recombination of fast ions impacting surfaces when the plasma potential is large (typically ∼ 40V inside the discharge of ITER-prototype plasma sources for instance). As stated above, the proportion of hydrogen molecules in a vibrational state ν ≥ 4 (H * 2 ) is assumed to be 1.5% of the background density. Negative ions are produced inside the ion source volume by the dissociative attachment process of H * 2 (reaction #17 of table I). Lastly, the time step is ∆t 0.2 ns and convergence is reached in about 100 µs. The floating potential at steady state is V ep -14V and is calculated self-consistently from the charges accumulating on the end-plates,

V ep = C -1 s q s (z ep , T ) , (6)
where z ep is the physical location of the plates, q s (z ep , T ) the total charge collected on the surfaces during a time interval T for each particle species s and C a numerical capacitance. The capacitance should be chosen large enough to modify the potential value only after an integration over a significant number of time steps (this will reduce fluctuations). The voltage in eq. ( 6) self-adjusts until steady-state conditions are reached. We assumed a value C = 10 -6 F, corresponding for instance to a potential increase of ∆V ep = 1V for I = 1A of net current The normalisation factor φm is equal to 4.1V for a bias voltage of -5V, -0.4V for a floating end-plate, -3.4V and lastly -7.7V for end-plate potentials of -20V and -25V, respectively. The potential is somewhat uniform inside the quasi-neutral plasma region in the last two cases, with variations below 0.1V in the calculation.

impacting the surface over ∼ 5000 time steps ∆t (which is equivalent to about 1 µs) with

∆V ep = I∆t C . (7)
Note that transient effects occurring on time scales shorter than the capacitor response will be smoothed out (which is not an issue in our case as we are investigating steady-state situations ou slower plasma dynamics).

The amplitude of the plasma potential in fig. 2 is found to be close to 0V; its profile has a u-shape radially with a difference between the maximum at r = x 2 + y 2 = 0 and the minimum around 2.8 cm of ∆φ 0.15 eV. This may confine electrostatically some positive ions (these particles will be lost on the end-plates instead). The plasma density is maximum on the linear device axis with n p 1.5 × 10 18 m -3 while the electron temperature T e is 6.5 eV and decreases to about 1 eV at the location where the negative ion density profile peaks. The latter has a donut-like shape with a density of n H -10 16 m -3 on axis and a maximum 6 times larger at a radius r = 3.4 cm. This is the consequence of the cross-section of reaction #17 in table I, which is largest for T e 0.3 eV (σ 10 -15 cm 2 ), is halved for an electron temperature of 1 eV and is 4 orders of magnitude lower for 5 eV. [START_REF] Janev | Elementary Processes in Hydrogen-Helium Plasmas[END_REF]

B. Scan versus the end-plate voltage

The simulations performed in this subsection are for a magnetic field strength of 100G. Electrons are strongly magnetized with a Hall parameter h e ∼ 40 (corresponding to a gyroradius ρ e = v th /Ω e 0.6 mm for T e = 3 eV). Electrons are hence attached to a magnetic field line and slowly diffuse perpendicular to the latter through collisions against ions (Coulomb) and neutrals (the ionisation fraction being less than 1% in the model). The magnetic field lines are short-circuited by the end-plates and furthermore electrons can gain energy axially, i.e., along (Oz), and hence escape across the plasma sheath in that direction. One way to reduce losses and consequently to lengthen the residence time of the electrons is to bias negatively the end-plates such as to increase the height of the sheath potential drop. This is shown in fig. 3 where the average electron residence time τ e is increased by a factor ∼ 10 when the end-plate voltage varies from V ep = 0V to -20V while the average electron density increases from n p 1.7 × 10 17 m -3 to 6 × 10 17 m -3 (2.4 kW absorbed power), respectively. The simulation parameters are provided in fig. 1 together with a numerical resolution of 256 2 grid nodes, 40 ppc, a heating frequency ν h = 10 7 s -1 and a scaling factor α = 3200. τ e is calculated by integrating the electron continuity equation over the plasma volume (assuming steady-state conditions), that is,

I w = S en e u e dS = V en e ν i dV , (8)
where w is the electron current impacting the device walls, S is the corresponding wall surface, u e the electron velocity at the wall, V the plasma volume and ν i the net electron production frequency inside the plasma (which accounts for the balance between ionisation and losses). Further assuming that τ e = ν -1 i is constant over the entire volume, one get,

τ e = eπR 2 L z n e I w (9)
with,

n e = πR 2 L z -1 V n e dV . (10)
Forcing electrons to diffuse across the magnetic field lines and to be lost on the linear device cylinder walls instead of on the end-plates has an impact on the plasma potential profile as shown in fig. 4. The height of the plasma potential decreases when the end-plate bias voltage is increased (in absolute value). The potential profile eventually becomes a well for the ions (it has a u-shape); the transition occuring when the end-plates are floating in the model. Positive ions are then lost on the end-plates because the depth of the potential well largely exceeds the ion temperature. The volume averaged electron density as well as the electron current impacting the device walls are saturating when the end-plates' bias voltage is greater than a couple of times the electron temperature (in electron-Volts) as shown in fig. 3

(a) and (b).

The electron residence time is correlated with the probability of generating a negative ion via the dissociative attachment of the hydrogen molecule by electron impacts. The longer the electrons stay in the plasma, the higher the negative ion production yield. Figure 5 shows the negative ion current lost in volume by collisions (reaction #16 of table I, #10, 11, 15 and 16 of table II, respectively) and on the walls of the device (including the endplates). The production rate is the sum of the two currents. The total negative ion current jumps from ∼ 2A for end-plates at ground potential up to ∼ 10A for -20V, respectively, which is the maximum. For comparison the discharge current is ∼ 12A (summed over the positive ion currents) and the electron current lost on the walls is ∼ 7.5A for the same bias (-20V), indicating that the negative ion yield may be substantial. Note that this does not translate into a large negative ion density because a large fraction of these ions is destroyed through collisions with electrons and neutrals. As a consequence, negative ion losses inside the plasma volume amount to more than 50% of the production rate regardless of the value of the end-plate bias voltage as shown in fig. 5. This is a disadvantage of volume production of negative ions compared with surface production on the extractor electrodes. The ion to electron volume averaged density ratio n H -/ n e is less than 10% in the calculations (see fig. 3). Furthermore, the plasma potential confines negative charges for configurations where the end-plate voltage lies between ground potential and floating conditions (fig. 4) but beyond that point, the potential accelerates negative ions radially. The latter hence increasingly convey the negative charges toward the cylinder wall, which is correlated in turn with a decrease of both the losses in volume and the negative ion density (see fig. 3). Note that there aren't any negative ion losses on the end-plates.

The model is simulating plasma densities of ∼ 5 × 10 17 m -3 on average with a plasma sheath larger than real due to the implementation of a scaling factor α = ε * /ε 0 in order to artificially increase the size of the electron Debye length and hence reduce the number of grid nodes. Figure 6 shows the plasma potential profile versus a numerical resolution of 256 2 and 512 2 nodes corresponding to a factor α = 800 and 3200, respectively. In both cases, the potential well is accurately reproduced with a similar width (defined as the cylinder radius minus the plasma sheath length) and depth. Lastly, an a-posteriori verification of the EEDF shows that the numerical technique simulating the RF power absorption by the macro-electrons in the model generates a local Maxwellian distribution. A detailed inclusion of the RF-plasma coupling dynamics might result in a non-Maxwellian EEDF which will affect the inelastic collision rates, impacting in turn the dissociative attachment and hence the negative ion density but the conclusions derived in this work from the parameter scans will remain identical. Figure 7 shows the EEDF in both the region of power absorption and around the location where the negative ion density is maximum for the case of an end-plate bias voltage of -8V.

IV. COMPARISON WITH EXPERIMENTS ON THE RAID LINEAR DEVICE

Figure 8(a) plots the normalised electron density n e /n em , temperature T e /T em and negative ion density n H -/n m profiles obtained with the 2.5D PIC-MCC model for the conditions of the experiments on the RAID testbed with n em = 1.5 × 10 18 m -3 , T em = 5.2 eV and n m = 1.7 × 10 16 m -3 , respectively. The absorbed power is assumed to be 3.1 kW with a Gaussian profile and a standard deviation of σ h = 1.5 cm (best fit in order to reproduce the experimental electron temperature profile). RAID has a 20 cm radius and a plasma column length L z = 1.5 m. The calculation of fig. 8 corresponds to a magnetic field strength of 200G and end-plates at ground potential. The simulation parameters are a mesh of 512 2 nodes, 100 ppc, a heating frequency ν h = 10 8 s -1 and a scaling factor α = 3200. The characteristics of the hydrogen gas are given in fig. 1. Furthermore, we fixed the proportion of H * 2 which are vibrationally excited (ν ≥ 4) to 0.25% of the total density (which is an external parameter to the model) to approximately fit the peak value of the negative ion density observed in the experiments.

The electron density and temperature in the experi- ment are derived from axially movable Langmuir probe (LP) measurements calibrated with interferometry while the absolute negative ion density is obtained using Cavity Ring Down Spectroscopy (CRDS) [START_REF] Agnello | [END_REF] combined with LP-assisted laser photodetachment. 16 The experimental electron and negative ion density radial profiles together with the electron temperature (dashed-lines) are shown in fig. 8(a) for comparison with the estimates from the 2.5D PIC-MCC calculation (solid lines) while the plasma potential is plotted in fig. 8(b). n em = 2 × 10 18 m -3 , T em = 4.3 eV, n m = 1.1 × 10 16 m -3 and the error bar for the plasma potential is ±1V in the experiment. The model is hence capable of reproducing the general characteristics observed in the experiments, i.e., a similar amplitude and profile for the plasma parameters. In addition, the maximum plasma potential on axis is about 11.5V in the calculation compared to ∼ 12V for a radius R < 6 cm in RAID. The negative ion density profile peaks in both cases in the area where the electron temperature drops to ∼ 1 eV and is significantly lower on the discharge axis. This indicates an annular profile similar to the one of fig. 2. Note that the ions are somewhat magnetized which results in the formation of rotating structures in the model in the area where the plasma pressure gradient is significant (R < 6 cm). A similar behaviour has been observed in the experiments on RAID but in Argon gas and for a magnetic field strength of 800G. 50

V. EFFECT OF THE ION MAGNETISATION ON THE PLASMA PROPERTIES

In this section, we analyse the incidence on the plasma characteristics of magnetising the positive ions. One obtains quiescent plasma conditions in the linear device model when the ions remain unmagnetized as shown in fig. 2. Increasing the strength of the magnetic field lengthens the residence time of the electrons as their transverse mobility is lowered (the Hall parameter evolves from h e ∼ 40 to ∼ 160 between 100G and 400G, respectively). This in turn results in higher plasma densities (for the same absorbed power) and narrower plasma profiles. This is shown by comparing fig. 2 and fig. 9 which corresponds to identical simulations except for the magnetic field strength, which is 400G in the latter case. Note that we have performed a convergence test for this configuration and did not find any noticeable differences when increasing the number of particles per cell to 400 instead of 40 or changing the number of grid nodes from 256 2 to 512 2 , respectively. Similar conclusions were drawn when studying the effect of the artificial heating frequency ν h (in a configuration with a simplified physical-chemistry) where we varied the fraction of macro-electrons in the heating region which absorbs the external power from ν h ∆t = 0.1% to 1% without loss of generality.

A transition in plasma behaviour occurs when the ions become somewhat magnetized. [START_REF] Futtersack | Modélisation fluide du transport magnétisé dans les plasmas froids[END_REF] In fig. 9, one observes the formation of large scale structures and the plasma is rotating clockwise with a period of ∼ 75 µs (13 kHz frequency). The entire plasma volume (except the sheath), including the rotating arm is quasi-neutral and hence the Debye length is not a critical parameter, justifying the use of a scaling factor for the vacuum permittivity in the model without any detrimental effects on the plasma dynamics. The diamagnetic drift velocity,

v Dβ = - ∇P β × B q β n β B 2 , (11)
leads to a counterclockwise rotation for the electrons (β corresponds to any charged particle species and q β its respective charge). The radial electric field E r is positive inside the bulk where the rotating arm is located while the pressure gradient is negative everywhere. The structure hence rotates azimuthally in the E×B direction (taking B z > 0). The arm rotation appears to be approximately rigid, i.e., with the same frequency regardless of its radial location. The corresponding velocity hence varies linearly radially and is about 2 km/s for instance at mid-radius which is the same order of magnitude as the average ion velocity; the latter being defined as,

u i = β (n β u β )/ β n β , (12)
where u β is the β ion specie average velocity and n β the corresponding density. The ions also rotate clockwise with a velocity profile comparable to the E × B drift. Their respective root-mean-square (RMS) velocities, u 2 i and v 2 E calculated inside a disk of radius r ≤ 4 cm and a center aligned with the discharge axis (in order to exclude the plasma sheath) are about 1 km/s in V ep (V) both cases where

v E = E × B B 2 . ( 13
)
For comparison, the RMS proton sound speed u 2 B 13 km/s (u B = eT e /m i ) and the RMS average electron velocity u 2 e ∼ 18 km/s. Figure 10 plots the rotating arm velocity together with the RMS ion and E × B drift velocities versus the polarisation of the end-plate bias plates. One observes that u 2 i and v 2 E are similar and decrease for larger bias voltages (in absolute value) which is caused by a lower plasma potential amplitude as shown in fig. 4. The arm velocity is also decreasing but at a smaller pace with respect to the ion velocity. The detailed assessment of the physics leading to the formation of these large scale rotating structures and the associated modified plasma transport properties (notably the question of anomalous current) is left to future work. Note that we observe similar behaviours in electro-positive gases, which have also been detected experimentally [START_REF] Rebont | [END_REF]53 , and hence the negative ions do not seem to be seeding these instabilities although they likely affect the latter (this is currently an open question). The negative ions have a ring-like density profile (fig. 9) which is also rotating but distorted by the varying electric field strength along their path. The electronegativity is n H -/n e ∼ 0.5 near the peak indicating that the negative ions represent a significant fraction of the negative charges in that area.

VI. EXTRACTION OF NEGATIVE IONS

Figure 5 shows that negative ions are trapped inside the plasma volume and destroyed by collisions with charged particles and neutrals when the plasma potential amplitude is above the linear device cylinder wall voltage (which is at ground potential in the calculation). Biasing the end-plates negatively acts as a control knob chang- ing the shape of the radial plasma potential profile such that negative ions may be directly extracted. The latter approach is not interesting for an ion source as the ion radial flux has a cylindrical symmetry. In a Neutral Beam Injector, one would need to extract the ions through slit or cylindrical apertures embedded on the cylinder surface and connected to an electrostatic accelerator in order to accelerate the particles to high energies (1 MeV in the case of ITER for instance). Note that a similar extraction scheme for negative ions produced in a linear device has been studied experimentally using SF 6 gas. 54 The numerical resolution of the PIC algorithm cannot implement the details of the aperture geometry. We therefore chose to model instead the impact on the negative ion dynamics of a small surface on the cylinder wall which is biased positively with respect to the plasma potential. This mimics the penetration of the accelerator extraction voltage inside the plasma, in the vicinity of the meniscus.

The biased surface has a width of 1 cm in the model and a voltage of 5V. The numerical resolution is 1024 2 grid nodes, 40 ppc, -10V end-plate bias potential, a heating frequency ν h = 10 7 s -1 and a scaling factor α = 50. The remaining parameters are identical to those of fig. 1. The maximum value of the electron Debye length in the cal-culation is about 100 µm and the time step ∆t 50 ps, respectively. Figure 11 plots the 2D plasma potential profile inside the linear device together with the negative ion current density. The plasma potential increases sharply near the biased slit as expected (RHS of the plot) but there is also a minimum around mid-radius with respect to the amplitude on axis (∆φ 0.4V). This traps a sub-population of negative ions in that area which therefore cannot be extracted. The average negative ion current density on the slit aperture is j H -35 A/m 2 (i.e., 7% of the negative ion current produced in the plasma volume), which is an order of magnitude lower than the ITER NBI requirement for instance. j H -is correlated with a maximum negative ion density of 4 × 10 16 m -3 about 1.5 cm from the cylinder's surface. The possibility to extract a larger current remains an open question as the negative ions are produced inside the plasma volume as opposed to the direct vicinity of the apertures where the ions are scattered off a cesiated surface after the impact of a neutral or a positive ion. The negative ion density profile peaks in the region where the electron temperature falls below 1 eV in the model. The location of the peak may hence be controlled and shifted toward the cylinder walls by modifying the extent of the power absorption region.

VII. CONCLUSION

We studied the feasibility of producing and extracting a hydrogen negative ion current generated inside the plasma volume of a linear device (through physicalchemistry processes between charged particles and the background gas) with a magnitude relevant to the requirements of magnetic fusion Neutral Beam Injectors (NBI). Linear devices are interesting configurations as the electrons are strongly magnetised with field lines aligned with the discharge axis. The large aspect ratio of the discharge (defined as the axial length over the cylinder radius) implies that electrons will oscillate between the end-plates located at the extremities of the cylinder and diffuse slowly across the magnetic field. The residence time of the electrons inside the ion source has an incidence on the production yield of the negative ions. Negative ions are produced inside the plasma in a two step process consisting of (i) an inner volume with an electron temperature between 5 and 10 eV where the hydrogen background gas is excited vibrationally and (ii), an outer volume, with a temperature below 1 eV where the negative ions are generated as a byproduct of the dissociative attachment of a hydrogen molecule due to a collision with an electron. We showed that biasing negatively the end-plates increases the electron residence time by an order of magnitude and as a consequence the negative ion density may be increased by a factor ∼ 3 for a given absorbed RF power. Biasing the end-plates negatively causes a reduction of the amplitude of the plasma potential which is a limiting factor for the negative ion density as the potential profile switches to a well when the end-plate voltage is such that one collects a higher ratio of positive ions compared to electrons. In that case, the plasma potential extracts radially the negative ions. One observes that the negative ion density has a ring-like shape with a maximum at the location where the electron temperature falls below 1 eV. This is consistent with experimental measurements on the RAID linear device. The position of the maximum may be moved by adjusting the radius of the power absorption region in the model. We simulated the extraction of negative ions by biasing positively (+5V) a surface on the cylinder wall with a width of 1 cm. The extracted ions amounted to about 7% of the total current produced in the source volume but this remained about an order of magnitude lower that the value fulfilling the requirement for the ITER NBI for instance (which is ∼ 250 A/m 2 ). To which extent increasing the helicon power is correlated also with a larger negative ion density in the experiments (before reaching any saturation mechanisms induced among others by neutral depletion in low pressure conditions, physical-chemistry, etc.) remains an open question. Lastly, large scale rotating structures appear when the ions become somewhat magnetised. These instabilities rotate at a frequency between 13 and 29 kHz in the calculation. The latter increases in correlation with a larger plasma potential amplitude suggesting a relationship with the E × B drift and hence the ion velocity. The detailed understanding of the physical mechanisms seeding these instabilities is left to future work.

16 FIG. 2 .

 162 FIG. 2. 2D electron density, temperature, plasma potential and negative ion density profiles for the conditions of fig. 1. |Bz| = 100G and the end-plates are floating with Vep -14V.

3 )FIG. 3 .

 33 FIG. 3. (a) Electron residence time τe (blue solid line) inside the plasma of the linear device, volume averaged electron density ne (in red) and (b), volume averaged negative ion density nH-(in blue) together with the electron current impacting the cylinder walls (including the end-plates) versus the endplate bias voltage Vep.

FIG. 4 .

 4 FIG.4. Normalized radial plasma potential profiles versus the end-plate bias voltage for the conditions of fig.1. The floating potential corresponds to a voltage of -14V. The normalisation factor φm is equal to 4.1V for a bias voltage of -5V, -0.4V for a floating end-plate, -3.4V and lastly -7.7V for end-plate potentials of -20V and -25V, respectively. The potential is somewhat uniform inside the quasi-neutral plasma region in the last two cases, with variations below 0.1V in the calculation.

FIG. 5 .

 5 FIG. 5. Negative ion current either impacting the discharge walls (Iw) or lost inside the plasma volume (I vol ) versus the end-plate bias voltage Vep.

2 FIG. 6 .

 26 FIG.6. Radial plasma potential profile for an end-plate bias voltage Vep = -20V and a scaling factor either of α = 3200 (corresponding to a 256 2 grid nodes resolution) or 800 (512 2 nodes), respectively.

FIG. 7 .

 7 FIG.7. Electron Energy Distribution Function (EEDF) profile in the power absorption region (r ≤ 2 cm) and near the maximum in negative ion density (2 cm ≤ r ≤ 3 cm). The straight lines are Maxwellian fits. The magnetic field strength is 100G and the end-plate bias voltage -8V, respectively.

FIG. 8 .

 8 FIG. 8. Normalised radial profiles for the electron density ne/nem, temperature Te/Tem and negative ion density n H -/nm in (a). Comparison between the RAID experiments (dashed lines) versus the 2.5D PIC-MCC calculation (solid lines). nem = 1.5 × 10 18 m -3 , Tem = 5.2 eV, nm = 1.7 × 10 16 m -3 in the simulations and nem = 2 × 10 18 m -3 , Tem = 4.3 eV and nm = 1.1 × 10 16 m -3 for the experiments, respectively. The plasma potential is displayed in (b). The error bar is ±1V. Lastly, the absorbed power profile in the model is Gaussian with a standard deviation of 1.5 cm and an integrated power of 3.1kW.

17 FIG. 9 .

 179 FIG. 9. 2D electron density, temperature, negative ion density and plasma potential profiles for the conditions of fig. 1. |Bz| = 400G and the end-plates are floating with Vep -14V. The volume occupied by the plasma sheath is slightly below 10% of the simulation domain due to the use of a scaling factor α = 3200.

FIG. 10 .

 10 FIG. 10. Rotation velocity of the large scale structures measured at mid-radius r = 2.5 cm (grey), root-mean-square (RMS) average ion velocity (black) and E × B drift velocity (red) versus the linear device end-plate bias voltage. The RMS velocities are calculated inside a radius r ≤ 4 cm.

2 )FIG. 11 .

 211 FIG. 11. 2D plasma potential and negative ion current density profiles inside the plasma volume of the linear device. The current streamlines are shown in grey (lower plot). The simulation corresponds to a cylinder radius of 5 cm, an endplate potential of -10V, a slit aperture width of 1 cm and a bias of 5V, respectively. The numerical resolution is 1024 2 grid nodes. j H -35 A/m 2 on the slit aperture.

TABLE II .

 II Heavy particle processes.

	# Reaction	Cross section ref.
	1 H + 3 + H2 → H + 3 + H2 (elastic)	43
	2 H + 3 + H → H + 3 + H (elastic)
	3 H + 2 + H2 → H + 3 + H	42,43
	4 H + 2 + H2 → H2 + H + 2	43
	5 H + 2 + H → H + 2 + H (elastic)

44 6 H + + H → H + H + 45 7 H + + H → H + + H (elastic) 45 8 H + + H2 → H + + H2 (elastic) 43 9 H + + H2 → H + + H2 (inelastic, 4 proc.) 40-43 10 H -+ H → e + 2H 25 11 H -+ H → e + H2 25 12 H -+ H2 → H -+ H2 (elastic) 42 13 H -+ H → H -+ H (elastic) 42 14 H + + H -→ 2H (2 proc.) 25 15 H + + H -→ H + 2 + e 25 16 H -+ H2 → H2 + H + e 25 17 H -+ H → H + H - 46

(Coulomb)[START_REF] Fubiani | [END_REF]

ACKNOWLEDGMENTS

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This work was supported in part also by the Swiss National Science Foundation, Commissariat à l' Énergie Atomique (CEA) and the French Fédération de Recherche sur la Fusion Magnétique (FR-FCM).

AVAILABILITY OF DATA

The data that support the findings of this study are available from the corresponding author upon reasonable request.


Céline Bonnet

Panhong Gou

Simon Girel

Vincent Bansaye

Catherine Lacout

Karine Bailly

Marie-Hélène Schlagetter

Evelyne Lauret

Sylvie Méléard

Stéphane Giraudier
email: stephane.giraudier@aphp.fr

Combined biological and modelling Approach of Hematopoiesis: From Native to stressed Erythropoiesis

We developed here a new mathematical modelling approach which reconciliates native and stress hematopoiesis, focused on hematopoietic stem and progenitor compartments. We first proposed an erythropoiesis model leaning on a minimum of 6 cell-amplification compartments, able to reproduce native erythropoiesis. A phenylhydrazine-induced hemolytic stress was next applied, and in vivo data were used to estimate model parameters through an optimized algorithm and integrated regulatory processes. A reduction of all stem cell compartments was observed, due to a drastic differentiation without proliferation during 7 days, followed by a huge proliferation in all compartments including long-term hematopoietic stem cells, before returning to normal values. The 6-compartment model applied to stress erythropoiesis was complexified by integrating regulatory processes. Thus, regulation of hematopoiesis appears dispensable during native hematopoiesis but mandatory for stress hematopoiesis. In conclusion this multi-step and time-dependent model of immature hematopoiesis opens new opportunities to understand development of normal or pathological hematopoiesis.

INTRODUCTION

Hematopoiesis provides the lifelong supply of mature blood cells derived from a rare population of bone marrow (BM) multipotent hematopoietic stem cells (HSCs) [START_REF] Orkin | Hematopoiesis: an evolving paradigm for stem cell biology[END_REF][START_REF] Till | A direct measurement of the radiation sensitivity of normal mouse bone marrow cells[END_REF]. The lineage relationship between HSCs and mature cells was first proposed as a multistep process in which generations of diverse blood cells are coupled with movement through consecutive Hematopoietic Stem and Progenitor Cell (HSPC) amplification motors, from HSC to mature blood cells. HSPC division and differentiation are controlled by extracellular signals, and by intracellular networks. In the classical model, successive HSPC compartments, Long-Term HSC (LT-HSC), Short-Term HSC (ST-HSC), MultiPotent Progenitor (MPP), etc.. along hematopoietic differentiation process are defined on their immunophenotypic level, typically by the expression of cell surface marker combinations [START_REF] Eich | OMIP-059: Identification of Mouse Hematopoietic Stem and Progenitor Cells with Simultaneous Detection of CD45.1/2 and Controllable Green Fluorescent Protein Expression by a Single Staining Panel[END_REF][START_REF] Notta | Distinct routes of lineage development reshape the human blood hierarchy across ontogeny[END_REF]. In this context, since the 1960s, HSCs are typically functionally defined by their ability to sustain multi-lineage engraftment for an extended period of time upon serial transplantation into irradiated recipient mice [START_REF] Till | A direct measurement of the radiation sensitivity of normal mouse bone marrow cells[END_REF]. Recently, new generation of experimental tools allowing in situ analysis of HSC output has challenged the classical model by pointing out fundamental differences between physiological unperturbed hematopoiesis and stress hematopoiesis [START_REF] Busch | Fundamental properties of unperturbed haematopoiesis from stem cells in vivo[END_REF][START_REF] Carrelha | Hierarchically related lineage-restricted fates of multipotent haematopoietic stem cells[END_REF][START_REF] Chapple | Lineage tracing of murine adult hematopoietic stem cells reveals active contribution to steady-state hematopoiesis[END_REF][START_REF] Pei | Polylox barcoding reveals haematopoietic stem cell fates realized in vivo[END_REF][START_REF] Pei | Resolving Fates and Single-Cell Transcriptomes of Hematopoietic Stem Cell Clones by PolyloxExpress Barcoding[END_REF][START_REF] Rodriguez-Fraticelli | Clonal analysis of lineage fate in native haematopoiesis[END_REF][START_REF] Säwen | Murine HSCs contribute actively to native hematopoiesis but with reduced differentiation capacity upon aging[END_REF][START_REF] Sun | Clonal dynamics of native haematopoiesis[END_REF]. In particular, in native hematopoiesis, progenitors directly downstream of LT-HSCs (ST-HSC, MPP…) serve as a major, nearly self-renewing source of day-to-day hematopoiesis, rendering blood and immune system less dependent on LT-HSC, in agreement with the quiescence of adult LT-HSCs [START_REF] Busch | Unperturbed vs. post-transplantation hematopoiesis: both in vivo but different[END_REF][START_REF] Busch | Fundamental properties of unperturbed haematopoiesis from stem cells in vivo[END_REF][START_REF] Foudi | Defining hematopoietic stem and progenitor cell turnover by analysis of histone 2B-GFP dilution[END_REF][START_REF] Schoedel | The bulk of the hematopoietic stem cell population is dispensable for murine steady-state and stress hematopoiesis[END_REF]. In addition, recent single cell profiling analysis has unveiled the tremendous heterogeneity of each previously defined compartment leading to a continuum in hematopoiesis process [START_REF] Giladi | Single-cell characterization of haematopoietic progenitors and their trajectories in homeostasis and perturbed haematopoiesis[END_REF][START_REF] Laurenti | From haematopoietic stem cells to complex differentiation landscapes[END_REF][START_REF] Paul | Transcriptional Heterogeneity and Lineage Commitment in Myeloid Progenitors[END_REF] Electronic copy available at: https://ssrn.com/abstract=3777468 [START_REF] Velten | Human haematopoietic stem cell lineage commitment is a continuous process[END_REF]. However, punctuated transitions across this continuous gene expression landscape may still exist and represent functionally distinct groups of cells [START_REF] Liggett | Unraveling Hematopoiesis through the Lens of Genomics[END_REF].

Mathematical modelling integrating both unperturbed and stress hematopoiesis should be a powerful approach to address key questions in hematopoiesis, able to provide qualitative and quantitative insights into stem cell dynamics and fate commitment. To date, most of hematopoiesis modelling is based on the existence of few cell compartments (stem cells, progenitors and mature cells) [START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF][START_REF] Crauste | Mathematical study of feedback control roles and relevance in stress erythropoiesis[END_REF]). Rodewald's group has introduced a major breakthrough by modelling hematopoiesis based upon stem compartments from LT-HSC to MPP in unperturbed hematopoiesis [START_REF] Busch | Fundamental properties of unperturbed haematopoiesis from stem cells in vivo[END_REF]. Nevertheless, they notice that their steady state model cannot recapitulate stress hematopoiesis such as 5-fluorouracile (5-FU) treatment, suggesting that unperturbed and stress hematopoiesis do not follow the same biological rules.

Here we want to recapitulate in a single mathematical model stress and steady-state hematopoiesis, based on in vivo experiments. We focused our attention on erythroid response since erythroid cells represent ~95% of blood cells. We chose a "peripheral stress" through hemolytic anemia (phenylhydrazine administration known to destroy mature red blood cells, RBC) to avoid direct alterations of the stem cell compartments [START_REF] Hara | Erthropoietic precursors in mice with phenylhydrazineinduced anemia[END_REF][START_REF] Klinken | Phenylhydrazine stimulates lymphopoiesis and accelerates Abelson murine leukemia virusinduced pre-B cell lymphomas[END_REF]. So far, most of studies focused their attention on committed progenitors, the effects on more immature compartments being poorly addressed.

Here, we demonstrate that stress and unperturbed hematopoiesis can be combined in a single system based on a 6-compartment modeling. We conducted in silico modelling and in vivo experiments to integrate homeostasis and stress hematopoiesis. We observed that transient hematopoietic stress impacts all primitive HSPC compartments, with a large rebound effect inducing a delay in the hematopoietic system. A steady state model (i.e. with constant parameter Electronic copy available at: https://ssrn.com/abstract=3777468 values) cannot explain the fast return of RBC to equilibrium and a model whose regulation depends directly on the different compartment sizes cannot explain the observed rebound effect.

Then, these effects were taken into account in our in silico model: We integrated timedependent regulations depending indirectly on compartment sizes. Self-renewal and differentiation dynamics were modulated through the dynamics of two regulators (u r and u d ).

The model was finely calibrated using a novel stochastic optimization algorithm (CMA-ES) (Hansen and Ostermeier, 2001) based on the experimental observations .

RESULTS

We first developed a mathematical model able to recapitulate steady-state hematopoiesis from quiescent LT-HSC to RBC. Modelling quiescence/kinetics of proliferation/differentiation processes requires a theoretical number of compartments between LT-HSC and RBC. However, deciphering an exact compartment number remains elusive since nature and number of compartments are not standardized from one study to another. Using a mathematical modelling, we deduced that the theoretical number of compartments should rely upon at least 3 criteria for each compartment: 1-size, 2-duration of proliferative steps and 3-self-renewal versus differentiation capacities.

We assumed that an accurate model should integrate stochasticity and amplification in each compartment as well as plasticity in the number of intermediate compartments. Then, a stochastic multi-type branching process was designed to recapitulate native (unperturbed) hematopoiesis, leading to derivation of a macroscopic deterministic system. In vivo steady-state Electronic copy available at: https://ssrn.com/abstract=3777468 parameters applied to this model helped us to determine the minimal number of amplification motors (k cell compartments) necessary to obtain a lifelong stable RBC production.

Modelling of unperturbed early erythropoiesis

Our model describes a hierarchy of k cell compartments evolving over time. Each cell of a given compartment divides with a constant rate independently from the others and is characterized by its division rate, its self-renewal and differentiation capacity, and for the last compartment (corresponding to RBC), its death rate. The entire dynamics of hematopoiesis presented in Figure 1A and 1B is summarized as follows:

1-For any i=1, , k-2, a cell belonging to the i compartment divides at rate τi>0 in two cells of the same compartment with a probability 0<p i r<1 (self-renewal) or in two cells of compartment i+1 with a probability µi p i d (=µi (1-p i r)) (erythroid differentiation). The p i d probability corresponds to differentiation, and µi to commitment towards erythroid rather than nonerythroid lineages. The i=1 compartment corresponds to LT-HSC and i=k to mature RBC.

2-A (k-1)-type cell "divides" at rate τ k -1 >0. Upon division, the cell renews with probability 0<pr k-1 <1, and gives rise to 2 n k-type cells with probability pd k-1 =1-pr k-1 . This dynamic summarizes the n mitosis with negligible self-renewal observed at the end of erythroid differentiation. This amplification was simplified as 2 n RBC issued from 1 cell of the previous compartment since terminal erythroid differentiation has been largely modeled before and is not our main concern.

3 -A k-type cell corresponding to mature RBC dies at rate τk >0. The estimation HSPC death rate is estimated to 2-4% per day and then was neglected in our modeling [START_REF] Domen | The role of apoptosis in the regulation of hematopoietic stem cells: Overexpression of Bcl-2 increases both their number and repopulation potential[END_REF].

We were aware that stem cell division could be asymmetric (a cell of type i giving birth to two daughter cells, one of type i and the other of type i+1). However, from our data, it was not Electronic copy available at: https://ssrn.com/abstract=3777468 possible to distinguish symmetrical from asymmetrical divisions. We therefore neglected the latter for parameter identifiability purpose.

Furthermore, presuming that LT-HSC compartment is at equilibrium, the number of selfrenewing LT-HSC is equal to the number of differentiating LT-HSC, i.e p 1 d=p 1 r. Moreover, for any i ∈ {2,., k-1}, p i d-p i r >0 to ensure the stability of the system. For any type i=2,.,k-1, the icompartment differentiation factor is defined as Di=p i d-p i r. Since p i d+p i r =1, Di =1-2 p i r =2 p i d -1.

According to the literature [START_REF] Ethier | Markov Processes: Characterization and Convergence[END_REF], LT-HSC number was large enough to approximate the stochastic model by a differential equation system (Figure 1C); therefore compartment sizes at equilibrium can be expressed (Figure 1D). In particular, the i compartment size relies upon the one of compartment i-1 multiplied by 1/Di, meaning that the smaller the Di, the bigger the i compartment size.

As a first step, we assumed that the Di are equal to a positive D value and the µi probabilities equal to 1. Then a direct relationship can be computed between the LT-HSC number, their division rate, the D parameter, the n number of final mitoses, the RBC number, and their death rate (Figure 1E). The minimal number of amplification motors / compartments was computed based on this formula. Such method was drafted in Dingly et al. 2007 [START_REF] Dingli | Compartmental architecture and dynamics of hematopoiesis[END_REF] with a simple model which does not include temporal dynamics.

Determination of the minimal compartment number

The determination of the theoretical minimum number of compartments was based on Figure 1E. The RBC number and their death rate are known (ie≈10 10 cells/mouse, and 40-day life expectancy). LT-HSC number was evaluated from the proportion of LT-HSC (Lin-Sca-1 + c-kit+, LSK CD150 + CD48 -) estimated as described in Figure S1. Knowing that a C57Bl6 mouse displays ≈250x10 6 BM cells [START_REF] Boggs | The total marrow mass of the mouse: a simplified method of measurement[END_REF][START_REF] Chervenick | Quantitative studies of blood and bone marrow neutrophils in normal mice[END_REF][START_REF] Mahajan | A quantitative assessment of the content of hematopoietic stem cells in mouse and human endosteal-bone marrow: a simple and rapid method for the isolation of mouse central bone marrow[END_REF], the Electronic copy available at: https://ssrn.com/abstract=3777468

LT-HSC number was around 15248 +/-2200/mouse. The LT-HSC division rate was computed from the proportion of LT-HSC retaining their phenotype after 20 hours of BrdU incorporation at steady state (1.82 in [0-4,7] % of cycling LT-HSC) (Figure 1F and Figure S1). Hence LT-HSC division rate was approximated by 1/100 per day. Assuming that Di must increase as a function of i and since Dk=1 (no self-renewal), Di belongs to [D2;1]. As in Busch's report [START_REF] Busch | Fundamental properties of unperturbed haematopoiesis from stem cells in vivo[END_REF], the second compartment corresponds to the ST-HSC compartment. Then minimal value for D2 (D2≥0.16) was computed (Figure 1D) with bounds for division rate (Figure 1F and Table 1).

The n number of mitosis between the two last compartments was assessed as follows: we assume that the k-1 compartment is upstream of BFU-E (Burst-Forming-Unit Erythroid)enriched cells known to perform 11 cell divisions to in vitro generate mature RBC [START_REF] Li | Rate of Progression through a Continuum of Transit-Amplifying Progenitor Cell States Regulates Blood Cell Production[END_REF]. Therefore "n" should be over 12.

Taking into consideration all these parameters (Figure 1E), we computed that if k=5, then D=0.12; if k=6, then D=0.23; if k=7, then D=0.35. Therefore, the minimum compartment number to obtain D ∈ [0.16;1] was 6. As previously assessed, the last steps (MEP to RBC) was simplified as 2 n RBC from 1 MEP. We focused our attention on HSPC compartments and roughly summarized behavior of committed erythroid progenitors. To match these 6 compartments to cytometric analysis of hematopoiesis, we then in vivo evaluated the cell number and division rates of LT-HSC, ST-HSC, MPP, CMP, MEP, and MEP to RBC compartments (Table 2).

Parameter calibration of the mathematical model for steady-state erythropoiesis.

We first gave bounds to the division rate per day. By considering variations due to BrdU experimental procedure, we used bounds rather than median values to delineate division rates Electronic copy available at: https://ssrn.com/abstract=3777468 (Figure 1F). Bounds offered the opportunity to include standard deviation from one experiment to another and variability of the time interval between BrdU injection and analysis (20h+/-4h).

We then computed appropriate values of all parameters (Table 2). We took into consideration 1-population size (Table 1), 2μi probabilities of differentiation toward the erythroid lineage (based on Busch's data), 3-decreased probability of self-renewal along the differentiation process, and 4-bounds of the division rate of each cell type. The µi erythro-myeloid commitment probabilities were obtained considering Busch et al. data (Busch et al., 2015) and dividing the differentiation rate toward erythro-myeloid lineages by the sum of differentiation rates toward all lineages. Based on decreasing probability of self-renewal in the different compartments along the differentiation process, Di were computed allowing calibration of the model (Figure 1D). Our results (Table 2) showed that this steady state erythropoiesis model was able to reproduce unperturbed native hematopoiesis for at least two years, the normal mouse lifetime.

Steady-state model applied to acute erythroid stress.

In order to validate the 6-amplification compartment model in stress hematopoiesis, we induce an acute hemolytic anemia using phenylhydrazine (PHZ) administration (Figure S2A). After a 40%-hemolysis stress observed at day 3, a rapid recovery of RBC was evidenced after 10-16 days (Figure 2A). Our 6-step model was applied to recapitulate RBC recovery but the recovery time exceeded one month, in total discrepancy with observed biological data (Figure 2A). This indicated that, although able to fulfill steady state hematopoiesis, this model was unable to recapitulate the hematopoietic stress response just like previously observed by Busch et al. using their own model [START_REF] Busch | Fundamental properties of unperturbed haematopoiesis from stem cells in vivo[END_REF]. We then considered that our model requires additional regulations during stress hematopoiesis.

Electronic copy available at: https://ssrn.com/abstract=3777468

Analysis of compartment kinetics during acute hematopoietic stress.

In order to decrypt biological mechanisms responsible for this fast RBC recovery, we analyzed all compartment modifications after PHZ injection in blood, BM, and spleen. Briefly, during the first 3 days after injection, we observed a drastic RBC fall, and an increase in platelets and leukocytes counts (Figure 2A-C). While total number of BM cells was slightly reduced at day 3 (-25%) and enhanced at day 10 (+33%) (Figure 2G), spleen weight strongly increased 24h after PHZ injection (Figure 2D-F), becoming 3 times larger than normal at day 3, as previously observed after erythroid [START_REF] Hara | Erthropoietic precursors in mice with phenylhydrazineinduced anemia[END_REF][START_REF] Sanchez | Differential amplification of murine bipotent megakaryocytic/erythroid progenitor and precursor cells during recovery from acute and chronic erythroid stress[END_REF] and social stress in mice [START_REF] Mckim | Social Stress Mobilizes Hematopoietic Stem Cells to Establish Persistent Splenic Myelopoiesis[END_REF]. The three hematopoietic tissues were analyzed to assess proportion and absolute cell numbers in the different compartments (LT-HSC, ST-HSC, MPP, CMP, and MEP) as well as their percentages into cell cycle (BrdU incorporation), and apoptosis (AnnexinV labelling) during the hemolytic phase. In BM, a decrease in the size of all compartments (from LT-HSC to CMP) was evidenced. BM MEP population remained relatively stable during this phase, while a huge increase in the number of splenic MEP was observed 3 and 5 days after PHZ injection (Figure 2G). No difference in BrdU incorporation was noticed at day 3 in the different BM compartments except for MEP, which displayed a 2fold increased cell proliferation, as compared to D0 (Figure S2F). To further get insight into mechanisms involved in the drastic decrease of compartment size, we analyzed apoptosis. No increase in the proportion of AnnexinV-labelled cells was observed whatever the progenitor compartment examined, allowing us to exclude any excess of cell death in the BM progenitor compartments (Figure S3B). Blood cell analysis performed at day 1 and 3 after PHZ injection showed a 15-fold increase in the number of mononuclear cells (Figure 2B) yet with a negligible proportion of Lin-cells, precluding an egress of LSK from BM to blood (Figure S3A). In spleen, a stress erythropoietic organ in mouse, the number of mature cells increased (3-fold) (Figure 2E) whereas LSK and CMP compartments displayed no significant variations in cell Electronic copy available at: https://ssrn.com/abstract=3777468 number, ruling out any role of splenic immature progenitors at day 3 (Figure 2G). Lastly, to confirm stemness properties of LT-HSC assessed by cytometry, competitive BM transplantation with 50% of CD45.1 BM cells from untreated animals and 50% of total BM cells from day 3-PHZ-treated mice were performed in CD45.2 9.5 Gy-irradiated recipients.

Blood chimerism analysis 3 months post-transplantation showed a decreased proportion of PHZ-treated cells, confirming the decrease number of LT-HSC in BM of day3-PHZ-treated mice (Figure S3C).

The decrease of BM progenitor compartment sizes, without increased proliferation or apoptosis or egress to the spleen or the blood, together with a stable number of MEP displaying a twofold increased proliferation in BM and spleen, was the result of a rapid differentiation process, like a flush, to quickly compensate loss of RBC through a versatile mechanism.

Recovery of progenitor compartment after PHZ stress.

To assess dynamics of progenitor compartments during recovery after PHZ treatment, we analyzed at days 3, 5, 7, 10, 16 and 28 the different compartment sizes and proliferation rates.

Recovery was almost effective in all compartments at day 7-10 (Figure 4A). Surprisingly, at day 10, all compartments displayed a 2-3-fold increased cell numbers before returning to normal values at day 28. In parallel, BrdU analysis showed that all BM compartments exhibited a drastic increased cell proliferation at day-10 (the time point of the hematocrit normalization).

Thereafter, BrdU incorporation returned to normal values at days 16 and 28 (Figure S2C-G).

Modelling steady-state and stress hematopoiesis by adding compartment regulations

The steady state model previously introduced could not recapitulate the fast return of RBC to equilibrium. Therefore, to model all at once unperturbed and stress hematopoiesis, we Electronic copy available at: https://ssrn.com/abstract=3777468 introduced regulations by assuming that previously defined parameters depend indirectly upon compartment size through dynamics of two regulators (u r and u d ) modulating self-renewal and differentiation, bringing positive and negative effects on compartment sizes, respectively (Figure 3A).

Dynamics of the two regulators followed a double system of production and clearance. The production is given by Michaelis-Menten type functions depending on BM cell numbers. The clearance of the regulator after its production is given by a linear degradation term (Figure 3B).

This modelling induced a delay between dynamics of the regulators and changes in compartment sizes, explaining the huge oscillations in cell numbers after PHZ injection.

Assuming that the proliferation rate during late erythropoiesis is regulated in a specific manner (Erythropoietin, Stem Cell Factor, etc…), RBC production dynamics from MEP is controlled with a third type of regulation. MEP differentiation gives birth to 2 n(t) RBC summarizing the n(t) last mitosis with a Hill function (Michaelis-Menten's type) of the RBC number (Figure 3C). In addition, we hypothesized that differentiation towards a non-erythroid lineage was not affected by all these regulators. Parameters of the regulated-model were calibrated using a stochastic optimization algorithm based on a CMA-ES method (Covariance Matrix Adaptation in Evolution Strategy), which minimizes discrepancy between in silico simulated and experimental data (see Suppl. data). Optimal parameters obtained by this method are presented in Table 3 and the corresponding cell dynamics in Figure 4B. The later matched with the experimental data in all compartments in stress condition (the oscillation effect) as well as in steady state (second part of the trajectory). Moreover, division rate values were deduced as function of time for each cell type (Table 4).

Electronic copy available at: https://ssrn.com/abstract=3777468

DISCUSSION

Our present work provides a mathematical model able to reconcile steady state and stress hematopoiesis. Furthermore, it highlights a mechanism by which transient stress exposure leads to a protracted influence on all progenitor compartments including LT-HSC, leading to a fast repopulation of mature compartment.

To decipher mechanisms involved in the normalization of RBC number at steady state and after an acute stress, we first developed a compartmental model of hematopoiesis. The theoretical number of compartments sufficient to fulfill stable erythropoiesis was relied upon the number of LT-HSC and mature RBC. Our modelling proposes that 6-amplification compartments are sufficient to ensure steady state erythropoiesis. We hypothesized that the last compartment was resumed from MEP to RBC as a simple amplification factor according to the n number of mitosis required from 1 MEP to RBC. Based on the division number between BFU-E and erythrocytes, this n number could vary from 11 to 14 [START_REF] Li | Rate of Progression through a Continuum of Transit-Amplifying Progenitor Cell States Regulates Blood Cell Production[END_REF]. We neglected megakaryocytic differentiation that theoretically only corresponds to 20% of the MEP differentiation [START_REF] Mori | Prospective isolation of human erythroid lineage-committed progenitors[END_REF] and other myelo-lymphoid cells representing quantitative minorities. Then, this 6-compartment model parallels the generally assessed compartment hierarchy with LT-HSC, ST-HSC, MPP, CMP, MEP and RBC.

Our "stress model" was initially the destruction of RBCs by chemical hemolysis. Nevertheless, this process also induced a huge inflammatory response leading to changes in all blood cell types and a huge cytokine storm in the next few days after hemolysis (data not shown),

confirming that stress hematopoiesis relied on regulatory mechanisms different from those of unperturbed hematopoiesis, as previously illustrated by the megakaryocytic and the granulocytic differentiation processes [START_REF] Haas | Inflammation-Induced Emergency Megakaryopoiesis Driven by Hematopoietic Stem Cell-like Megakaryocyte Progenitors[END_REF][START_REF] Kang | Deregulated Notch and Wnt signaling activates early-stage myeloid regeneration pathways in leukemia[END_REF]. Based on our data and modelling, we demonstrate here that (u r and u d ) regulations of compartment proliferation

Electronic copy available at: https://ssrn.com/abstract=3777468

and differentiation are required to reproduce stress response while dispensable for steady-state erythropoiesis, these regulators remaining constant in the latter.

Previous mathematical models based on in vivo experiments mostly described hematopoiesis dynamics at steady state. These studies rely upon parameters such as cell division rates or probabilities of self-renewing [START_REF] Abkowitz | In vivo kinetics of murine hemopoietic stem cells[END_REF][START_REF] Busch | Fundamental properties of unperturbed haematopoiesis from stem cells in vivo[END_REF][START_REF] Foudi | Defining hematopoietic stem and progenitor cell turnover by analysis of histone 2B-GFP dilution[END_REF][START_REF] Klose | Hematopoietic Stem Cell Dynamics Are Regulated by Progenitor Demand: Lessons from a Quantitative Modeling Approach[END_REF][START_REF] Mackey | Cell kinetic status of haematopoietic stem cells[END_REF][START_REF] Manesso | Dynamical modelling of haematopoiesis: an integrated view over the system in homeostasis and under perturbation[END_REF][START_REF] Sawai | Hematopoietic Stem Cells Are the Major Source of Multilineage Hematopoiesis in Adult Animals[END_REF]. Models integrating steady state as well as in vivo stress hematopoiesis rarely reproduce oscillations observed in the size of the progenitor compartments apart from the modelling [START_REF] Angulo | Investigating the role of the experimental protocol in phenylhydrazine-induced anemia on mice recovery[END_REF][START_REF] Crauste | Adding self-renewal in committed erythroid progenitors improves the biological relevance of a mathematical model of erythropoiesis[END_REF][START_REF] Crauste | Mathematical study of feedback control roles and relevance in stress erythropoiesis[END_REF][START_REF] Klose | Hematopoietic Stem Cell Dynamics Are Regulated by Progenitor Demand: Lessons from a Quantitative Modeling Approach[END_REF][START_REF] Loeffler | A mathematical model of erythropoiesis in mice and rats. Part 1: Structure of the model[END_REF][START_REF] Manesso | Dynamical modelling of haematopoiesis: an integrated view over the system in homeostasis and under perturbation[END_REF][START_REF] Marciniak-Czochra | Modeling of asymmetric cell division in hematopoietic stem cells--regulation of self-renewal is essential for efficient repopulation[END_REF][START_REF] Roeder | Competitive clonal hematopoiesis in mouse chimeras explained by a stochastic model of stem cell organization[END_REF]. These reports modelized the regulation dynamics by a differential equation with a delay, enabling to predict the oscillating behavior of progenitor compartments characterized by a regulated self-renewal. However they integrated only a theoretical 3-compartment model (stem cells/progenitors/mature erythrocytes) based on RBC recovery data after stress. Other regulation-based models have been developed but did not integrated in vivo data [START_REF] Kirouac | Dynamic interaction networks in a hierarchically organized tissue[END_REF][START_REF] Mahadik | A computational model of feedbackmediated hematopoietic stem cell differentiation in vitro[END_REF].

In contrast, our present 6-compartment model was designed to reproduce the in vivo dynamics observed throughout unperturbed and stress erythropoiesis, using 8 differential equations, including two for regulator dynamics and naturally created a delay in the oscillation phenomenon due to the duration needed for regulator production and clearance and then better match with hematopoiesis process.

Little is known regarding in vivo consequences of a "peripheral" stress on upstream progenitor compartments such as LT-, ST-HSC, MPP, CMP, and MEP. Our experiments demonstrated that after an acute stress destroying RBC, the first compensatory mechanism entails all compartments, from LT-HSC to CMP, and surprisingly, induces as a first step, differentiation

Electronic copy available at: https://ssrn.com/abstract=3777468

without proliferation just like a flush. Thereafter, all depleted compartments enter into proliferation during a recovery phase, the division rate of each progenitor compartment reaching its maximum value around day 11 but still remained lower than 3 divisions per day (in accordance with biological data) (Table 4). This recovery phase not only replenishes the different "flushed" compartments but exceeds steady-state values. Normalization of all compartment sizes occurs in a third phase just like spring oscillation around the standard values.

These oscillations in cell number from each compartment are in accordance with in vivo data previously reported in post-chemotherapy stressed hematopoiesis [START_REF] Busch | Fundamental properties of unperturbed haematopoiesis from stem cells in vivo[END_REF]. In contrast, the MEP compartment was not subjected to such a high regulation process in the bone marrow probably because of the well-known spleen stress erythropoiesis process that follows a different kind of regulation as reported here and in previous reports [START_REF] Perry | Maintenance of the BMP4-dependent stress erythropoiesis pathway in the murine spleen requires hedgehog signaling[END_REF].

In conclusion, we modeled the plasticity of hematopoiesis recapitulating steady-state as well as stress hematopoiesis, no compartment being "protected" from stress consequences. Our mathematical model suggests that regulators of hematopoiesis are dispensable during steadystate hematopoiesis. Thus, our model opens the way to better understand malignant clonal development and invasion: Stem cell disorders just like leukemias and more generally cancers could arise from successive stress if pathological clonal cells react differently from normal cells in response to aggression.

MATERIALS AND METHODS

Animals

We used adult 6-12 weeks old C57BL/6 Ly5.2 and Ly5.1 mice obtained from Janvier (Le Genest, France). Mice were bred and maintained in pathogen-free conditions in our Institute Animal facilities Unit for Laboratory Animal Medicine at the Université de Paris. All procedures performed were approved by local Committee on the Use and Care of Animals and by the Institutional Animal Care and Use Committee.

Hematologic evaluations

Orbital plexus blood was collected in EDTA tubes from anesthetized mice. Blood cell counts were determined using automated blood coulter (MS9, Schloessing Melet, Cergy-Pontoise, France). Blood samples were used for cytometry analysis. Serum were collected, stored at -20°C before assessing cytokine levels.

BM and spleen Analysis

Mice were killed by cervical dislocation. Bones (femurs, tibias, pelvis and humerus) and spleen were harvested, muscle and tendon tissue were removed using a scalpel and kimwipes. BM and spleens were used for cell count, and cytometry.

BM fraction was flushed out using a syringe containing 1xPBS complemented with 2% Fetal Bovine Serum (FBS). The resulting cell suspension were filtered through a 40µM cell strainer (Corning, NY, USA) and pelleted by centrifugation.

For mechanical grinding, spleens were smashed and ground between rough sides of frosted glass slides, and cells were collected in DMEM containing 2%FBS. After incubation in a 24well plate for 30min at 37°C in a humidified incubator, cell suspensions were passed through a 40µm cell strainer, and cells were re-suspended in DMEM containing 10%FBS. Suspensions were subjected to cytometry analysis after RBC lysis.

Electronic copy available at: https://ssrn.com/abstract=3777468

Cytometry analysis

Erythrocytes were lysed using ACK Lysing Buffer (Lonza, Basel, Switzerland) before flow cytometry stainings. Total BM and spleen cells were stained with Zombie Viability kit for 15min at room temperature, and thereafter stained with biotinylated anti-mouse Lineage

Cocktail antibodies for 30 min (Biolegend). After washing, cells were stained for 30 min using the following monoclonal antibodies in the BD Horizon Brilliant staining buffer: anti-CD117/c-Kit (2B8)-BB700 and anti-CD34 (RAM34)-AF647 (BD Pharmingen); anti-Ly6a/Sca-1 (D7)-BV510; anti-CD150 (TC15-12F12.2)-PE-Cy7; anti-CD48 (HM48-1)-BV711 ; Streptavidin APC-Cy7 (Biolegend).

For proliferation analysis, surface stainings were performed as described above, and cells were proceeded for cell proliferation according to manufacturer's recommendations (BD Pharmingen BrdU Flow kits) with BrdU-AF488 (BD Biosciences, Franklin Lakes, NJ, USA).

Data acquisition and/or data analysis were performed on the Cochin Cytometry and Immunobiology Facility Cytometry analyses on a Fortessa cytometer (Becton Dickinson), and analyses were done on Kaluza software.

To assess cell viability, surface labelled cells were resuspended in AnnexinV binding buffer, and AnnexinV (FITC) was added 15min before analysis (FITC-conjugated AnnexinV labeling detection kit, BD Pharmingen).

Reconstitution assays

Donor cells were isolated from 6-12-week-old B6.SJL-Ptprc a Pepc b /BoyCrl Congenic mice (Ly5.1) and Ly5.2 mice after Phenylhydrazine (PHZ) (60 mg/Kg, IP) injection or not. From these mice, hind limbs were extracted and cleaned. Total BM were flushed, then passed through a 70µM cell strainer to obtain a single cell suspension. A mix of cells from PHZ-treated mouse and untreated mouse was prepared in a 1:1 ratio. A mix of 5x10 6 cells were resuspended in a total volume of 0.2 mL and transplanted into 5 irradiated (9Gy) CD45.1 C57/BL6 recipient mice. Chimerism analysis was performed three months post-reconstitution using cytometry.

Statistical analysis

Results are presented as mean+/-SD. Data were analyzed through 2-tailed Student t-test.

TABLES Table 1: Determination of the cell number of each compartment per mouse and the % of

BrdU positive cells in steady state hematopoiesis. Each value is the mean of at least 6 mice.

Table 2: Determination of the division rate/day, the size, and the Differentiation factor for each compartment at steady state. Electronic copy available at: https://ssrn.com/abstract=3777468

Michaelis-Menten's type constant of the ur and ud regulator. mr and md : clearance rate of the ur and ud regulators. (C). Regulation of the terminal erythroid production (from MEP to RBC).

2 n(t) : number of RBC (k-type cells) generated after differentiation of a (k-1)-type cell at time t. 2 (See supplemental information).

Electronic copy available at: https://ssrn.com/abstract=3777468 Table 1 presents the % of BrdU + cells of each compartment at the steady state, which corresponds to the proportion of cells having performed at least one division during 20 hours .

Here, we will relate these informations to the steady-state division rates of each cell type: τi * .

Before BrdU injection, the system is at equilibrium and all cells are unlabeled by BrdU. The numbers of BrdU+ and unlabeled i-type cells at time t of the experiment are denoted by Ni m (t)

and Ni(t). Then Ni m (0)=0 and Ni(0)=xi * .

During the experiment, when an unlabeled cell of i type divides, it becomes BrdU+. Thus the dynamics of unlabeled i-type cell number is expressed as an exponential decrease at τi* rate.

Throughout the experiment, the system is considered to be at steady state. Therefore, for any i = 1,...,k-1, and for any t ∈[0,24h], Ni m (t) + Ni(t) = xi * . The formula in Figure 1F is then deduced.

2-Calibration of the regulated-erythropoiesis model

Parameters of the regulated model were calibrated from data of the 6 population sizes dynamics after PHZ induced hemolysis. In order to model the return to the equilibrium (given in Figure 4B), we considered day 3 as the initial time. Hence initial conditions of the model are given by : x1(day 3) = 0.006 . 10 6 cells/mouse, x2(day 3) = 0.0152 . 10 6 cells/mouse, x3(day 3) = 0.02 . 10 6 cells/mouse, x4(day 3) = 0.093 . 10 6 cells/mouse, Electronic copy available at: https://ssrn.com/abstract=3777468

x5(day 3) = 0.395 . 10 6 cells/mouse, x6(day 3) = 0.94 . 10 10 cells/mouse.

Using a stochastic optimization algorithm based on a CMA-ES method (Covariance Matrix Adaptation Matrix in Evolution Strategy), developed on Python by the INRIA team RandOpt (Hansen and Ostermeier, 2001) and least-squares cost function, we minimized the discrepancy between in silico simulated and experimental data. Results were summarized in Table 3. In such algorithms, identifiability is numerically accepted when identical results for a large number of independent runs are obtained.

Figure 1 :

 1 Figure 1 : Mathematical parameters required to model steady state Erythropoiesis. (A)

Figure 2 :

 2 Figure 2: Analysis of Blood, Spleen and BM parameters after PHZ-mediated chemical

Figure 3 .

 3 Figure 3. Mathematical parameters required to Model stress Erythropoiesis . (A) Effects

Figure 4 :

 4 Figure 4: Evolution of experimental and model parameters after PHZ-mediated chemical

  Figure1, Bonnet et al.Electronic copy available at: https://ssrn.com/abstract=3777468

Figure

  Figure 3. Bonnet et al.

Figure S2 :

 S2 Figure S2: Analysis of BrdU incorporation in HSPCs at different time points after PHZ

Figure S3 :C

 S3 Figure S3: Analysis of peripheral blood, BM and spleen after PHZ administration. (A)

Table 3 : Constants of proliferation and erythroid differentiation rate of compartments ; regulator parameter values and mitosis number from MEP to mature RBC after stress.

 3

Table 4 : Division rates at different time points in the different compartments after PHZ administration.

 4

Electronic copy available at: https://ssrn.com/abstract=3777468

ACKNOWLEDGMENTS

The authors are grateful to C. Chomienne for helpful discussions and improving the English manuscript. The authors greatly acknowledge the Cochin Cytometry and Immunobiology Facility and the IRSL Animal Facility, Université de Paris.

Supported by grants from the INSERM and INCa 2018 and fellowships from The Research ministry of China.

Supplemental Figure 1 Electronic copy available at: https://ssrn.com/abstract=3777468