Wei Liu

Xiao-Gang Lu
email: xglu@t.shu.edu.cn

Qing-Miao Hu
email: qmhu@imr.ac.cn

Hao Wang
email: haowang@i.shu.edu.cn

Yi Liu
email: yiliu@t.shu.edu.cn

Pascal Boulet
email: pascal.boulet@univ-amu.fr

L U Xiao-Gang

H U Qing-Miao

Boulet Pascal

Atomic bonding and electronic stability of the binary sigma phase

Keywords: A. sigma phase, C. atomic bonding, C. electronic stability, C. DOS, C. ELF, D. firstprinciples calculations

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

Introduction

The sigma phase is a non-stoichiometric intermetallic compound, that crystallizes in a tetragonal structure with 30 atoms distributed on five nonequivalent sites denoted as 2a, 4f, 8i1, 8i2 and 8j [START_REF] Hall | The sigma phase[END_REF][START_REF] Joubert | Crystal chemistry and Calphad modeling of the σ phase[END_REF][START_REF] Liu | Influence of nearest neighbor atoms and coordination polyhedron on atomic volume of sigma phases[END_REF] (see Fig. 1). It mostly forms in transition element systems, but it also exists in two systems containing a main-group element, i.e., A-Al (A=Nb, Ta). The sigma phase has attracted great experimental and theoretical research interests. The thermodynamic [START_REF] Dupin | A thermodynamic database for Ni-base superalloys[END_REF][START_REF] Mathieu | CALPHAD description of the Mo-Re system focused on the sigma phase modeling[END_REF] and molar volume [START_REF] Liu | Modeling of molar volume of the sigma phase involving transition elements[END_REF] databases of the sigma phase were built successively to predict the precipitation and molar volume of the sigma phase in multicomponent systems by using the CALPHAD method [START_REF] Kaufman | Computer Calculations of Phase Diagrams[END_REF][START_REF] Saunders | CALPHAD (calculation of phase diagrams): a comprehensive guide[END_REF][START_REF] Lukas | Computational thermodynamics: the CALPHAD method[END_REF]. Physical properties of the sigma phase, such as site occupancy [START_REF] Joubert | Crystal chemistry and Calphad modeling of the σ phase[END_REF][START_REF] Sluiter | Prediction of Site Preference and Phase Stability of Transition Metal Based Frank-Kasper Phases[END_REF][START_REF] Liu | Influencing factors of atomic order in the binary sigma phase[END_REF], elastic properties [START_REF] Bonny | Elastic properties of the sigma W-Re phase: A first principles investigation[END_REF], molar volume [START_REF] Liu | Influence of atomic mixing and atomic order on molar volume of the binary sigma phase[END_REF], enthalpy of formation [START_REF] Crivello | χ and σ Phases in Binary Rhenium-Transition Metal Systems: a Systematic First-Principles Investigation[END_REF][START_REF] Liu | Influence of atomic order on the enthalpy of formation and bulk modulus of the sigma phase[END_REF] and magnetism [START_REF] Cieślak | Magnetic moments and Curie temperatures of σ-FeCr alloys[END_REF][START_REF] Cieślak | Magnetic ordering above room temperature in the sigma-phase of Fe66V34[END_REF] have drawn significant concern as well.

Electronic structure analyses are effective approaches to explore the physical and crystallographic properties as well as the stability [START_REF] Crivello | χ and σ Phases in Binary Rhenium-Transition Metal Systems: a Systematic First-Principles Investigation[END_REF] of the sigma phase. Berne et al. [START_REF] Berne | Ab initio study of metastability in refractory metal based systems[END_REF] reported that coordination number 12 (CN12) sites are preferable for atoms with filled, nearly filled or empty d shells. The reason is that CN12 sites have approximately icosahedral symmetry, causing high degeneracy of d-like electronic levels [START_REF] Berne | Ab initio study of metastability in refractory metal based systems[END_REF][START_REF] Sluiter | Site occupation in the Cr-Ru and Cr-Os σ phases[END_REF]. Grånäs et al. [START_REF] Grånäs | Theoretical study of the Mo-Ru sigma phase[END_REF] proposed that, for Mo-Ru sigma compounds, Mo atoms prefer to occupy large CN sites, favoring the stronger covalent bonding between atoms.

In the present work, the first principles plane wave pseudo-potential method and exact muffin-tin orbitals method in combination with coherent potential approximation were used to reveal the bonding between atoms of the sigma phase. To facilitate explanation, all the studied binary sigma compounds were designated as A-B, where atom A has a larger atomic size than atom B. The ordered and disordered state is designated with atomic occupancy as presented in Table 1.

Methodology and calculation details

First-principles calculations were performed by using the exact muffin-tin orbitals method in combination with coherent potential approximation (EMTO-CPA) and plane wave pseudo-potential method implemented in the Vienna ab initio simulation package (VASP). The EMTO-CPA method was adopted mainly to deal with chemical disordered sigma compounds. By using this method, we can conveniently obtain the electronic density of states (DOS) for the sigma compounds with different atomic occupancies on lattice sites. We used VASP to calculate the electron localization function (ELF)

for the completely ordered sigma compounds. The calculation details are as follows.

EMTO-CPA calculations

The Green's function technique is used to solve the one electron Kohn-Sham equation within the exact muffin-tin orbitals (EMTO) method [START_REF] Vitos | Computational quantum mechanics for materials engineers: the EMTO method and applications[END_REF][START_REF] Vitos | Total-energy method based on the exact muffin-tin orbitals theory[END_REF]. The optimized overlapping muffin-tin approximation is conducted when dealing with the effective potential in the one-electron equation.

Additionally, the full charge density (FCD) method [START_REF] Vitos | Computational quantum mechanics for materials engineers: the EMTO method and applications[END_REF] is used to correct the total energy. The wave function is expanded by using the basis sets of the exact muffin-tin orbitals [START_REF] Vitos | Computational quantum mechanics for materials engineers: the EMTO method and applications[END_REF][START_REF] Vitos | Total-energy method based on the exact muffin-tin orbitals theory[END_REF]. The coherent potential approximation (CPA) [START_REF] Soven | Coherent-potential model of substitutional disordered alloys[END_REF][START_REF] Gyorffy | Coherent-potential approximation for a nonoverlapping-muffin-tinpotential model of random substitutional alloys[END_REF][START_REF] Vitos | Anisotropic Lattice Distortions in Random Alloys from First-Principles Theory[END_REF] is incorporated within the EMTO code, which facilitates the calculations dealing with chemical disordered alloys.

In the present work, Green's function was calculated for 16 complex energy points distributed exponentially on a semicircular contour. We adopted the scalar-relativistic and soft-core approximations. The electronic exchange-correlation functional within the generalized-gradient approximation (GGA) as parameterized by Perdew et al. [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] was used. The Brillouin zone was sampled by a uniform k-point mesh (3×3×6) without any smearing technique. The DOS were obtained at the equilibrium volume.

VASP calculations

The plane wave method with projector augmented wave (PAW) pseudo-potentials [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] based on density-functional theory (DFT) [START_REF] Dreizler | Density Functional Theory[END_REF] is implemented in VASP [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF]. The exchange-correlation functional within the generalized gradient approximation (GGA) as parameterized by Perdew and Wang namely GGA-PW91 [START_REF] Perdew | Accurate and simple analytic representation of the electron-gas correlation energy[END_REF] was employed. The plane-wave cutoff energy was set as 400 eV. The k-point meshes (8×8×15) for Brillouin zone sampling were constructed using the Monkhorst-Pack scheme [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF]. To facilitate the calculations, ZenGen script-tool [START_REF] Crivello | a tool to generate ordered configurations for systematic first-principles calculations: The Cr-Mo-Ni-Re system as a case study[END_REF] was used to automatically generate input files. The ELF was obtained using the relaxed structure at the equilibrium volume.

Results and discussion

In order to gain better understanding of the chemical bonding in the sigma phase, both electronic density of states (DOS) and electron localization function (ELF) calculations were conducted for several binary sigma phase systems. The calculation results indicate that the atomic bonds of the sigma phase exhibit both metallic and covalent characters for both A-Al (A=Nb, Ta) and transition metal systems (TM-TM). The covalent bonding for A-Al (A=Nb, Ta) is stronger than that for TM-TM.

Electronic stability of A-Al (A=Nb, Ta) sigma systems

To clarify the electronic stability of A-Al (A=Nb, Ta) sigma systems, the total electronic density of states (TDOS) of the Nb65.9Al34.1, Nb66.7Al33.3, Ta65.5Al34.5, Ta66.7Al33.3 and Ta75.6Al24.4 sigma compounds were computed with the site occupancies taken from literature [START_REF] Joubert | Crystal chemistry and Calphad modeling of the σ phase[END_REF][START_REF] Brown | The structure of the σ-phase Nb2Al[END_REF][START_REF] Boulineau | Structural characterization of the Ta-rich part of the Ta-Al system[END_REF][START_REF] Edshammar | The σ-phase Ta2Al[END_REF] as presented in Fig. 2. We can observe that for all the compounds, the TDOS at the Fermi level are finite, which indicates metallic bonding between atoms. For Nb65.9Al34.1, Nb66.7Al33.3, Ta65.5Al34.5 and Ta66.7Al33.3, the Fermi level falls within the pseudogap, which results from hybridization between the electronic states, and represents a covalent bonding between atoms [START_REF] Hu | Theoretical investigation of the ω -related phases in TiAl -Nb /Mo alloys[END_REF][START_REF] Ravindran | Electronic structure, chemical bonding, phase stability, and ground-state properties of YNi 2 -x ( C o / C u ) x B 2 C[END_REF]. Normally, covalent bonding is stronger than metallic bonding [START_REF] Hu | Theoretical investigation of the ω -related phases in TiAl -Nb /Mo alloys[END_REF].

For the completely ordered A66.7Al33.3 (A=Nb, Ta) compounds, the pseudogaps at the Fermi level are deeper, and the densities of states at the Fermi level are smaller than that for Nb65.9Al34.1 or Ta65.5Al34.5/Ta75.6Al24.4 with atomic disorder. It indicates a stronger atomic bonding and thus a larger electronic stability for A66.7Al33.3 (A=Nb or Ta) than for other Nb-Al or Ta-Al configurations.

Additionally, the completely ordered A66.7B33.3 (i.e., B2aA4fA8i1B8i2A8j) compound holds a small mismatch between atom and lattice site, namely large atom A occupy large CN sites (i.e., 4f, 8i1 and 8j); small atom B occupy small CN sites (i.e., 2a, 8i2) [START_REF] Liu | Influencing factors of atomic order in the binary sigma phase[END_REF]. Therefore both the electronic stability and size mismatch favor the stability of the completely ordered A66.7Al33.3 (A=Nb, Ta) compounds. Hence, the sigma compounds of A-Al (A=Nb, Ta) systems are observed experimentally to be highly ordered (see eg. Ref. [START_REF] Joubert | Crystal chemistry and Calphad modeling of the σ phase[END_REF][START_REF] Liu | Influencing factors of atomic order in the binary sigma phase[END_REF] and references therein).

To further understand the electronic stability of the A66.7Al33.3 (A=Nb, Ta) compounds, we analyzed in detail their TDOS and partial DOS (PDOS) as presented, respectively, in Figs. 2 and3, respectively.

We observe that the TDOS of Ta66.7Al33.3 shift to lower energy as compared to that of Nb66.7Al33.3 due to a stronger directional d bonding between Ta atoms [START_REF] Hu | Geometric and electronic structure of Ti2AlX (X=V, Nb, or Ta)[END_REF].

Besides, there is a peak located at around -0.22 Ry for Nb66.7Al33.3 and -0.26 Ry for Ta66.7Al33.3, which comes from the bonding state made up of TM-d (TM=Nb, Ta) and Al-p orbitals [START_REF] Zou | Structural, electronic, and magnetic properties of 3 d transition-metal aluminides with equiatomic composition[END_REF] as can be observed from the PDOS in Fig. in the figure, were adopted from the literature, namely Ref. [START_REF] Joubert | Crystal chemistry and Calphad modeling of the σ phase[END_REF] for Nb65.9Al34.1, Ref. [START_REF] Brown | The structure of the σ-phase Nb2Al[END_REF] for Nb66.7Al33.3, Ref. [START_REF] Boulineau | Structural characterization of the Ta-rich part of the Ta-Al system[END_REF] for Ta65.5Al34.5 Ref. [START_REF] Edshammar | The σ-phase Ta2Al[END_REF] for Ta66.7Al33.3 and Ref. [START_REF] Boulineau | Structural characterization of the Ta-rich part of the Ta-Al system[END_REF] for Ta75.6Al24. (Color online).

The ELF of the completely ordered A66.7Al33.3 (A=Nb, Ta) compounds were also calculated as presented in Fig. 4, where the range of ELF values is from 0 to 0.76. Theoretically, when ELF values equal 0.5, the Pauli repulsion at that position has the same value as in a uniform electron gas of the same density; when ELF values are close to 0 and 1, the local Pauli repulsion is very large and very small compared to that in a uniform electron gas [START_REF] Kohout | Electron localization function for transitionmetal compounds[END_REF], respectively. In other words, ELF=0.5 represents the electron-gas like pair probability; ELF=0 and ELF=1 represent perfect delocalization and localization, respectively [START_REF] Shi | Segregation effects of Y, Ti, Cr and Si on the intergranular fracture of niobium[END_REF]. For more details see Ref. [START_REF] Kohout | Electron localization function for transitionmetal compounds[END_REF][START_REF] Savin | ELF: The Electron Localization Function[END_REF]. In general, metallic bonding is supposed to be delocalized; covalent bonding is supposed to be localized [START_REF] Shi | Segregation effects of Y, Ti, Cr and Si on the intergranular fracture of niobium[END_REF][START_REF] Gao | Electronic origin of the anomalous solid solution hardening of Y and Gd in Mg: A first-principles study[END_REF].

For the completely ordered A66.7Al33.3 (A=Nb, Ta), where Nb and Ta occupy 4f, 8i1, 8j sites, and Al occupy 2a, 8i2 sites, the ELF results in Fig. 4

Influence of electronic configuration on atomic bonding for TM-TM sigma systems

To figure out the influence of electronic configuration, we calculated the DOS for Cr-Fe, Cr-Co, Re-Mn, Os-Cr, Re-Cr, Re-V, Ru-Cr, Mo-Mn, Mo-Re and Mo-Os binary sigma systems, but to avoid redundancy only the results of typical Cr-Fe, Mo-Re, Re-Mn and Re-V systems are shown in the present work.

The Cr-Fe system can serve as a prototype when investigating the influence of the valence electron factor on atomic bonding, where Fe (3d 6 4s 2 ) has two more valence electrons than Cr (3d 5 4s 1 ). Fig. 5 (a) presents the site-decomposed DOS of the two Cr49.5Fe50.5 sigma compounds one with site occupancies from the literature and the other in hypothetically disordered state. For clarity, only 2a and 4f sites are presented, as other sites also support the following argument. It indicates that when more Cr atoms occupy a specific crystal site, the atomic bonding on the site becomes stronger, which manifest itself by a deeper pseudogap and a smaller density of states at the Fermi level.

The Re-Mn system can serve as a prototype when investigating the influence of the electron shell factor on atomic bonding, where Re (5d 5 6s 2 ) has two more electron shells than Mn (3d 5 4s 2 ). Fig. 5 (b)

shows the site-decomposed DOS of the two Re51.2Mn48.8 sigma compounds one with site occupancies from the literature and the other in hypothetically disordered state. It indicates that when more Re atoms occupy a specific crystal site, the atomic bonding on the site becomes stronger.

In fact, the above arguments on Cr-Fe and Re-Mn are universal for the binary TM-TM sigma phase.

For a TM-TM sigma compound, the constituent atom, with less valence electrons or more electron shells, holds higher electronic stability. When increasing its atomic occupancy on a specific site, the atomic bonding on the site becomes stronger. To further validate the above statement, we calculated both the site-decomposed DOS and ELF for the completely ordered Cr66.7Fe33.3 and Re66.7Mn33.3 sigma compounds. The site-decomposed DOS in Fig. 6 indicate a stronger bonding between Cr atoms and a weaker bonding between Fe atoms for Cr66.7Fe33.3, and a stronger bonding between Re atoms and a much weaker bonding between Mn atoms for Re66.7Mn33.3. The above results agree well with the ELF calculations for Cr66.7Fe33.3 and Re66.7Mn33.3

shown in Figs. 7 and8, respectively. Note that for the ELF of Cr66.7Fe33.3, the features of the outercore region, spatially separated from the common valence region, is due to the unequal occupation of d orbital [START_REF] Kohout | Electron localization function for transitionmetal compounds[END_REF].

Discussion

Site occupancy preference on CN12 sites of the sigma phase

Berne et al. [START_REF] Berne | Ab initio study of metastability in refractory metal based systems[END_REF] and then Crivello et al. [START_REF] Crivello | χ and σ Phases in Binary Rhenium-Transition Metal Systems: a Systematic First-Principles Investigation[END_REF] clarified that CN12 sites are preferable for atoms with filled, nearly filled or empty d shells, the reason being that CN12 sites have approximate icosahedral symmetry causing high degeneracy of electronic d-like levels [START_REF] Berne | Ab initio study of metastability in refractory metal based systems[END_REF][START_REF] Sluiter | Site occupation in the Cr-Ru and Cr-Os σ phases[END_REF]. Thus atoms with half-filled d shells would have high density of states when occupying CN12 sites and consequently occupy high CN sites [START_REF] Crivello | χ and σ Phases in Binary Rhenium-Transition Metal Systems: a Systematic First-Principles Investigation[END_REF]. However, the present calculations show disagreement with the above argument.

For Re-V sigma system, where Re (5d 5 6s 2 ) bears half-filled d shells, we calculated the site decomposed DOS of the two Re78V22 compounds with different site occupancies as shown in Fig. 9.

When the five crystal sites are equally occupied by 78 at.% Re, as presented in Fig. 9 (a), the site decomposed DOS of each site near the Fermi level are close to one another, with a relatively higher value on 4f site (i.e., CN15 site). Apparently, CN12 sites do not cause high density of states. On the other hand, when CN12 sites are totally occupied by Re, as presented in Fig. 9 (b), it leads to low density of states on CN12 sites. This is against the argument that atoms with half-filled d shells would have high density of states when occupying CN12 sites [START_REF] Crivello | χ and σ Phases in Binary Rhenium-Transition Metal Systems: a Systematic First-Principles Investigation[END_REF].

The present calculations indicate that for the TM-TM binary sigma phase, the constituent atom with more valence electrons or less electron shells would have high density of states near the Fermi level, no matter whether it occupies large or small CN sites. In fact, CN12 sites are preferable for the constituent atom with smaller atomic size, more valence electrons or electron shells [START_REF] Liu | Influencing factors of atomic order in the binary sigma phase[END_REF]. This is related to the tendency of electron loss or gain of the two constituent atoms. For more details see Ref. [START_REF] Liu | Influencing factors of atomic order in the binary sigma phase[END_REF].

Influence of atomic bonding on molar volume of the sigma phase

The molar volume of the sigma phase varies with the site occupancy even at a fix composition. The influence of site occupancy on molar volume is related to the atomic bonding, and can be described by the volume difference between the completely disordered and ordered state of the A66.7B33.3 compounds (∆Vm=Vdisorder-Vorder) [START_REF] Liu | Influence of atomic mixing and atomic order on molar volume of the binary sigma phase[END_REF]. For most sigma phase systems, the volume of the completely ordered state is smaller than that of the disordered one, i.e., ∆Vm is positive. The only exception is for Mo-Re system, where ∆Vm is negative [START_REF] Liu | Influence of atomic mixing and atomic order on molar volume of the binary sigma phase[END_REF].

It needs to be clarified that for A66.7B33.3 in the completely ordered state, large atom A occupies large CN sites (i.e., 4f, 8i1 and 8j) and small atom B occupies small CN sites (i.e., 2a and 8i2). Comparing to the disordered state, there are more atoms A on large CN sites, and more atoms B on small CN sites.

According to the conclusion drawn in Section 3.2, when constituent A bears more electron shells or less valence electrons than B, the atomic bonding on large CN sites is stronger, and that on small CN sites is weaker for the completely ordered state than for the disordered one. The bonding effect on large CN sites should be more prominent, and thus leads to a smaller ordered volume than the disordered one [START_REF] Liu | Influence of atomic mixing and atomic order on molar volume of the binary sigma phase[END_REF].

On the contrary, when constituent A bears less electron shells or more valence electrons than B, it leads to a larger ordered volume than the disordered one. The only different case is the Mo-Re system, where Mo (4d 5 5s 1 ) has one less electron shell and valence electron than Re (6d 5 6s 2 ). It should be the electron shell factor that makes the volume of the ordered state larger than that of the disordered one [START_REF] Liu | Influence of atomic mixing and atomic order on molar volume of the binary sigma phase[END_REF]. We calculated the site-decomposed DOS for the two Mo48Re52 sigma compounds one with site occupancies from the literature (denoted as 'exp') and the other in hypothetically disordered state (denoted as 'dis') as presented in Fig. 10. The results indicate that the atomic bonding on 4f site of the ordered Mo48Re52 (exp) is obviously weaker than that of the disordered Mo48Re52 (dis), which is caused by the electron shell factor. For other crystal sites (i.e., 2a, 8i1, 8i2, 8j), no remarkable difference can be observed due to the competition between the electron shell factor and valence electron one. Thus, a smaller volume for the ordered Mo48Re52 (exp) is mainly induced by weaker atomic bonding on 4f site.

Conclusion

The present work based on the electronic density of states (DOS) and electron localization function (ELF) calculated by using first-principles methods allowed us to clarify the atomic bonding characteristics of the binary sigma phase including A-Al (A=Nb, Ta) and transition metal systems (TM-TM). The main conclusions are summarized as follows:

1. The atomic bonds of the binary sigma phase exhibit both metallic and covalent characters. The covalent bonding for A-Al (A=Nb, Ta) is stronger than that for TM-TM.

2. The completely ordered A66.7Al33.3 (A=Nb, Ta) sigma compounds have higher electronic stability and smaller size mismatch between atom and lattice site than the disordered configurations.

Hence, the A-Al (A=Nb, Ta) sigma compounds are found highly ordered in nature.

3. For a TM-TM sigma compound, the constituent atom, with more electron shells or less valence electrons, holds higher electronic stability. When increasing its atomic occupancy on a specific site, the atomic bonding on the site becomes stronger.

4. For Mo-Re system, the reason why the molar volume of the ordered Mo48Re51 compound is larger than that of the disordered Mo48Re51 is that the atomic bonds on 4f site for the ordered compound are weaker than that of the disordered one. It is dictated by the electron shell factor of the constituents, namely Mo bears less electron shells than Re.

Acknowledgments

Fig. 1 .

 1 Fig. 1. Crystal structure of the sigma phase. Atoms occupying different Wyckoff positions (i.e., 2a, 4f, 8i1, 8i2 and 8j) are indicated by different colors. (Color online).

3 .

 3 The Nb(d)-Al(p) bonding state is more localized than Ta(d)-Al(p), as indicated by a steeper bonding peak for Nb(d)-Al(p). Next to the peak, slightly higher in energy, there is a drop of the density of states, as presented in Fig.2. It is generated by the hybridization of the TM-d and Al-p orbitals[START_REF] Hu | Geometric and electronic structure of Ti2AlX (X=V, Nb, or Ta)[END_REF].Additionally, several small peaks appear on the top of the major bonding states just below the Fermi level from around -0.2 to 0 Ry. The first, energy-lowest small peak (at around -0.144 Ry for Nb66.7Al33.3 and -0.160 Ry for Ta66.7Al33.3) is mainly made up of TM-d and Al-p orbitals; the other small peaks correspond to the directional d bonding.

Fig. 2 .

 2 Fig.2. Total electronic density of states (TDOS) for A-Al (A=Nb, Ta) sigma compounds from EMTO-CPA calculations. The site occupancies of the compounds, shown in the inserted Table in the figure, were adopted from the literature, namely Ref.[START_REF] Joubert | Crystal chemistry and Calphad modeling of the σ phase[END_REF] for Nb65.9Al34.1, Ref.[START_REF] Brown | The structure of the σ-phase Nb2Al[END_REF] for Nb66.7Al33.3, Ref.[START_REF] Boulineau | Structural characterization of the Ta-rich part of the Ta-Al system[END_REF] for Ta65.5Al34.5 Ref.[START_REF] Edshammar | The σ-phase Ta2Al[END_REF] for Ta66.7Al33.3 and Ref.[START_REF] Boulineau | Structural characterization of the Ta-rich part of the Ta-Al system[END_REF] for Ta75.6Al24.4. 'ord' represents the completely ordered state. The vertical line indicates the Fermi level. (Color online).

Fig. 3 .

 3 Fig. 3. Partial electronic density of states (PDOS) for the completely ordered A66.7Al33.3 (A=Nb, Ta) sigma compounds (i.e., Al2aA4fA8i1Al8i2A8j) from EMTO-CPA calculations. The vertical line indicates the Fermi level.(Color online).

Fig. 4 .

 4 Fig. 4. Electron Localization Function (ELF) of the completely ordered Nb66.7Al33.3 (upper half) and Ta66.7Al33.3 (lower half) sigma compounds namely Al2aA4fA8i1Al8i2A8j (A=Nb, Ta). (a,d), (b,e) and (c,f) depict the ELF across (001), (110) and (-110) slice planes with distance from origin of 0 Å, 7.05 Å / 7.03 Å (Nb-Al/Ta-Al) and 0 Å, respectively. The color scale is given on the left. The scaling is between 0 and 0.76 (the maximum ELF for Nb66.7Al33.3 and Ta66.7Al33.3 is about 0.76 and 0.69, respectively). (Color online).

Fig. 5 . 5 * 8 *Fig. 6 .

 5586 Fig. 5. Site-decomposed DOS for (a) Cr-Fe and (b) Re-Mn sigma compounds from EMTO-CPA calculations.'exp' and 'dis' indicate compounds with site occupancies from the literature (Cr49.5Fe50.5[START_REF] Cieślak | Influence of composition and annealing conditions on the site-occupation in the σ-phase of Fe-Cr and Fe-V systems[END_REF] and Re51.2Mn48.8[START_REF] Joubert | Crystal chemistry and Calphad modeling of the σ phase[END_REF]) and in hypothetically disordered state, respectively. The vertical lines indicate the Fermi level. (Color online).

Fig. 7 .

 7 Fig. 7. Electron Localization Function (ELF) of the completely ordered Cr66.7Fe33.3 sigma compound namely Fe2aCr4fCr8i1Fe8i2Cr8j. (a), (b) and (c) depict the ELF across (001), (110) and (-110) slice planes with distance from origin of 0 Å, 6.14 Å and 0 Å, respectively. The color scale is given on the left. The scaling is between 0 and 0.47. Crystal sites with white labels are in the slice planes crystal site with black labels are in upper or lower layers of the slice planes. (Color online).

Fig. 8 .

 8 Fig. 8. Electron Localization Function (ELF) of the completely ordered Re66.7Mn33.3 sigma compound namely Mn2aRe4fRe8i1Mn8i2Re8j. (a), (b) and (c) depict the ELF across (001), (110) and (-110) slice planes with distance from origin of 0 Å, 6.48 Å and 0 Å, respectively. The color scale is given on the left. The scaling is

Fig. 9 . 20 * 20 *

 92020 Fig. 9. Site-decomposed DOS for Re78V22 sigma compounds from EMTO-CPA calculations. The site

  relative to Fermi level (Ry) *Y2a(Re) Y4f(Re) Y8i1(Re) Y8i2(Re) Y8j(Re) in the inserted Table in the figure, are artificially hypothetical. The vertical lines indicate the Fermi level. (Color online).

Fig. 10 .

 10 Fig. 10. Site-decomposed DOS for Mo48Re52 sigma compounds from EMTO-CPA calculations, where 'exp' and 'dis' indicate compounds with site occupancies from the literature [45] and in hypothetically disordered state, respectively. The vertical line indicates the Fermi level. (Color online).

Table 1

 1 Site occupancies of the completely ordered A2B and disordered AxB1-x sigma compounds.

	Wyckoff position	2a 4f 8i1 8i2 8j
	Coordination number (CN)	12 15 14 12 14
	Completely ordered state (A2B) [2]	0	1	1	0	1
	Disordered state (AxB1-x)	x	x	x	x	x

Nb 65.9 Al 34.1 Nb 66.7 Al 33.3 (ord) Ta 65.5 Al 34.5 Ta 66.7 Al 33.3 (ord) Ta 75.6 Al 24.4

  4. 'ord' represents the completely ordered state. The vertical line indicates the Fermi level. (Color online).

			AxAly	*Y2a(A) Y4f(A) Y8i1(A) Y8i2(A) Y8j(A)
		16	Nb65.9Al34.1 0.09	0.9	0.99	0.05 0.96
	(states/Ry)	12	Nb66.7Al33.3 Ta65.5Al34.5 0.04 0.95 0.96 0 1 1 Ta66.7Al33.3 0 1 1	0 0.03 0.98 1 0 1
	of states	8	Ta75.6Al24.4 0.32 0.97 0.96	0.34 0.97
	Total density	4
		0
			-0.6	-0.5		-0.4	-0.3	-0.2	-0.1	0.0
				Energy relative to Fermi level (Ry)

*Ys(A): site occupancy of atom A on crystal site s ( i.e. 2a, 4f, 8i 1 , 8i 2 , 8j).

  This work was supported by National Science Foundation for Young Scientists of China (Grant i.e. 2a, 4f, 8i 1 , 8i 2 , 8j). ). X.-G. Lu acknowledges the High Performance Computing Center of Shanghai University for providing the calculation facility.

	number: 51801119	states (states/Ry)	16 18 20	*Ys(Mo): site occupancy of Mo on crystal site s 0.113 0.751 0.638 0.1 0.715 0.495 0.495 0.495 0.495 0.495 Mo48Re52 *Y2a(Mo) Y4f(Mo) Y8i1(Mo) Y8i2(Mo) Y8j(Mo) exp dis	1 1 Mo 48 Re 52 8i 2 8i 1 dis 1 1 exp 2a 4f	8i 1 8i 2 2a 4f
		of					8j	1	8j
		density	14						4f
		Site-decomposed	10 12	-0.20	-0.15	-0.10	-0.05	0.00	0.05
					Energy relative to Fermi level (Ry)

(


Ramsès Djidjou-Demasse
email: ramses.djidjoudemasse@ird.fr

Arnaud Ducrot

Nicole Mideo

Gaëtan Texier

Understanding dynamics of Plasmodium falciparum gametocytes production: Insights from an age-structured model

Keywords: Within-host model, malaria, gametocytemia, parasitemia, infectiousness

 Biologically, the goal of these models is to understand what governs the severity of infections, the patterns of infectiousness, and the variation thereof across individual hosts. Mathematically, these models are based on dynamical systems, with standard approaches ranging from K-compartments ordinary differential equations (ODEs) to delay differential equations (DDEs), to capture the relatively constant duration of replication and bursting once a parasite infects a host red blood cell. Using malariatherapy data, which offers fine-scale resolution on the dynamics of infection across a number of individual hosts, we compare the fit and robustness of one of these standard approaches (K-compartments ODE) with a partial differential equations (PDEs) model, which explicitly tracks the "age" of an infected cell. While both models perform quite similarly in terms of goodness-of-fit for suitably chosen K, the K-compartments ODE model particularly overestimates parasite densities early on in infections when the number of repeated compartments is not large enough. Finally, the K-compartments ODE model (for suitably chosen K) and the PDE model highlight a strong qualitative connection between the density of transmissible parasite stages (i.e., gametocytes) and the density of host-damaging (and asexually-replicating) parasite stages. This finding provides a simple tool for predicting which hosts are most infectious to mosquitoes -vectors of Plasmodium parasites-which is a crucial component of global efforts to control and eliminate malaria.

Introduction

Malaria has always been a public health problem and, since the discovery of malaria parasites in human blood by Charles Laveran in 1880, remains so despite more than 100 years of research. Malaria continues to have a significant impact on the world with over 400,000 deaths alone each year [56]. It is a vector-borne disease caused by five plasmodial species: Plasmodium falciparum, P. vivax, P. malariae, P. ovale, and P. knowlesi, with P. falciparum being the most pathogenic species infecting humans [START_REF] Khoury | Within-host modeling of blood-stage malaria[END_REF].

The malaria parasite has a complex life cycle involving sexual reproduction occurring in the insect vector [START_REF] Alano | Sexual differentiation in malaria parasites[END_REF] and two stages of infection within a human (or animal host), a liver stage [START_REF] Frevert | Sneaking in through the back entrance: The biology of malaria liver stages[END_REF] and blood stage [START_REF] Bannister | The ins, outs and roundabouts of malaria[END_REF]. Human infection starts by the bite of an infected mosquito, which injects the sporozoite form of Plasmodium during a blood meal. The sporozoites enter the host peripheral circulation, and rapidly transit to the liver where they infect liver cells (hepatocytes) [START_REF] Frevert | Sneaking in through the back entrance: The biology of malaria liver stages[END_REF]. The parasite replicates within the liver cell before rupturing to release extracellular parasite forms (merozoites), into the host circulation, where they may invade red blood cells (RBCs) to initiate blood stage infection [START_REF] Miller | Malaria biology and disease pathogenesis: Insights for new treatments[END_REF]. Then follows a series of cycles of replication, rupture, and re-invasion of the RBC. Some asexual parasite forms commit to an alternative developmental pathway and become sexual forms (gametocytes) [START_REF] Russell | The Encyclopedia of Medical and Veterinary Entomology[END_REF]. Gametocytes can be taken up by mosquitoes during a blood meal where they undergo a cycle of sexual development to produce sporozoites [START_REF] Alano | Sexual differentiation in malaria parasites[END_REF], which completes the parasite life cycle.

The classical model of within-host parasite multiplication in malaria infections was formulated by Anderson el al. [START_REF] Anderson | Non-linear phenomena in host-parasite interactions[END_REF]. This model tracks uninfected red blood cells (RBCs), parasitized RBCs (pRBCs) and merozoites. The pioneer work of Anderson el al. [START_REF] Anderson | Non-linear phenomena in host-parasite interactions[END_REF], has been further developed in several directions including in particular immune response, see for instance [1, 25, 28-30, 36, 44, 45] for human malaria infection. We also mention discrete-time models such as in [START_REF] Dietz | Mathematical model of the first wave of plasmodium falciparum asexual parasitemia in non-immune and vaccinated individuals[END_REF]. Those models use an exponential process to describe the rate of rupture of pRBCs and, as a consequence, then fail to capture realistic lifetimes of the pRBCs on short time scales [START_REF] Saul | Models for the in-host dynamics of malaria revisited: Errors in some basic models lead to large over-estimates of growth rates[END_REF]. One reason for this is that they are essentially Markovian, i.e. 'memoryless', a RBC that has been parasitized for 40 hours has the same probability of producing merozoites as e.g. a RBC parasitized less than a hour ago. Moreover, those models are treating some processes that are likely to be kinda continuous as occurring only in a narrow window (e.g., the development of parasites within RBCs and the rupture of pRBC followed by the merozoites release phenomenon).

To correct this issue, some models of malaria infection include K-compartments ordinary differential equations (ODEs) representing a progression through a parasite's developmental cycle, e.g. [START_REF] Gravenor | Reply to: Models for the in-host dynamics of malaria revisited: Errors in some basic models lead to large over-estimates of growth rates[END_REF][START_REF] Iggidr | Global Analysis of New Malaria Intrahost Models with a Competitive Exclusion Principle[END_REF][START_REF] Saralamba | Intrahost modeling of artemisinin resistance in Plasmodium falciparum[END_REF][START_REF] Zaloumis | Assessing the utility of an anti-malarial pharmacokinetic-pharmacodynamic model for aiding drug clinical development[END_REF], or delay differential equations (DDEs) to capture the time pRBCs take to mature before producing new merozoites, e.g. [START_REF] Cao | Modeling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection[END_REF][START_REF] Hoshen | Mathematical modelling of the within-host dynamics of Plasmodium falciparum[END_REF][START_REF] Kerlin | Preferential Invasion by Plasmodium Merozoites and the Self-Regulation of Parasite Burden[END_REF][START_REF] Mckenzie | An integrated model of Plasmodium falciparum dynamics[END_REF][START_REF] Su | Periodicity and synchronization in blood-stage malaria infection[END_REF]. Other approaches are the use of partial differential equations (PDEs) to track the agestructure of the pRBC population [START_REF] Antia | The dynamics of acute malaria infections. I. Effect of the parasite's red blood cell preference[END_REF][START_REF] Cromer | Low red cell production may protect against severe anemia during a malaria infection-Insights from modeling[END_REF][START_REF] Demasse | An Age-Structured Within-Host Model for Multistrain Malaria Infections[END_REF][START_REF] Khoury | Within-host modeling of blood-stage malaria[END_REF]. It is shown in [START_REF] Fonseca | Comparison of Mathematical Frameworks for Modeling Erythropoiesis in the Context of Malaria Infection[END_REF] that DDEs perform better than the ODEs in representing the dynamics of red blood cells during malaria infection.

The K-compartments ODE model can be interpreted as the application of the method of stages (or the "linear chain trick") to the life cycle of pRBC, e.g. see [START_REF] Feng | Epidemiological models with non-exponentially distributed disease stages and applications to disease control[END_REF][START_REF] Hurtado | Generalizations of the 'Linear Chain Trick': Incorporating more flexible dwell time distributions into mean field ODE models[END_REF][START_REF] Iggidr | Global Analysis of New Malaria Intrahost Models with a Competitive Exclusion Principle[END_REF] and references therein. One problem with the K-compartments ODE model is how to decide upon the number of repeated compartments to capture the realistic dynamics of the pRBCs [START_REF] Gravenor | A Model for Estimating Total Parasite Load in Falciparum Malaria Patients[END_REF]. Determining the distribution of mean waiting times across compartments for the ODE model is also an issue. If the compartments can be considered equivalent to the developmental stages of pRBCs, then parasites might not spend equal time in each stage.

We first introduce both mathematical models (PDE and K-compartments ODE) and define the model's parameters and outputs. Next, using gametocyte production as a proxy variable of infectiousness, we compared the model outputs from a PDE stagestructured formulation to those from classical K-compartments ODE. Furthermore, the output of both mathematical models is used to qualitatively recover the time course of parasitemia, defined as the proportion of all infected RBCs among the total number of RBCs. Finally, the K-compartments ODE model (when K is properly chosen) and the PDE model are used to highlight a strong qualitative connection between gametocyte density and parasitemia.

2 Material and method

Data and methodology

Our analysis is based on data collected from malariatherapy taken in [START_REF] Eichner | Genesis, sequestration and survival of Plasmodium falciparum gametocytes: Parameter estimates from fitting a model to malariatherapy data[END_REF]. Malaria inoculation was a recommended treatment for neurosyphilis between 1940 and 1963. We also refer to [START_REF] Chernin | The Malariatherapy of Neurosyphilis[END_REF] for a review paper on malariatherapy and the knowledge gained in the understanding of malaria infection. The data we shall use consist in daily records of gametocyte density for twelve patients. Although malariatherapy has been dismissed for obvious ethical reasons, the advantages to use such data are multiple. Indeed, patients are naive to malaria infection and the dynamics are not perturbed by anti-malarial treatments. Let us notice that such data have been widely used in the literature and in particular to estimate mathematical model parameters. We refer to [START_REF] Eichner | Genesis, sequestration and survival of Plasmodium falciparum gametocytes: Parameter estimates from fitting a model to malariatherapy data[END_REF] and the references therein.

The method we shall develop consists of devising a mathematical model to describe the intra-host development of the infection and fitting the model to the available data. The output of the mathematical model will allow us to access various quantities related to the time course of the infection, including parasitemia.

Mathematical model

As discussed above, we now present the mathematical model we shall use to recover parasitemia for twelve patients from observed time courses of gametocyte density. We shall describe the within-host malaria infection coupled with red blood cells (RBCs) production as well as immune effectors. Fig. 1 presents the flow diagram of the model considered in this note. Our model is divided into four parts: (i) uninfected RBC (uRBCs) dynamics; (ii) changes in parasite stage or parasite maturity; (iii) Gametocyte production and dynamics and (iv) immune response dynamics.

For uRBCs dynamics, we divide cells into three age classes: reticulocyte (young), mature and senescent. All three ages are vulnerable to P. falciparum infection. This can be different for other species of Plasmodium. Although we focus in this work on the case of P. falciparum, the model described below could be applied to study other species such as P. vivax or P. malariae, which have specific RBC-age preferences [START_REF] Paul | The evolutionary ecology of Plasmodium[END_REF]. Such agestructured dynamics for uRBC are well known in the literature, see for instance [START_REF] Mcqueen | Host Control of Malaria Infections: Constraints on Immune and Erythropoeitic Response Kinetics[END_REF]. For the parasites, we consider stage-structured dynamics for their development within pRBC. [START_REF] Mcqueen | Host immune constraints on malaria transmission: Insights from population biology of within-host parasites[END_REF], Λ 0 is the RBC production rate from the marrow source. In our model, the parameter 1/µ r→m (resp. 1/µ m→s , 1/µ s→d ) is the time spent in RBC reticulocyte (resp. mature, senescent) class. A continuous parameter a denotes the time since the concerned RBC is parasitized: ring stage (0 < a < 26 hours), trophozoite (26 < a < 38 hours) and schizont (38 < a < 48 hours). In the case of P. falciparum infection, one has (γ r = γ m = γ s = 1) while for P. vivax one has (γ r = 1, γ m = γ s = 0) and for P. malariae (γ r = γ m = 0, γ s = 1) [START_REF] Paul | The evolutionary ecology of Plasmodium[END_REF].

Here the stage is a continuous variable representing the time since the concerned RBC is parasitized. Such a continuous stage structure will allow us to track the development of parasites within RBCs, but also to have a refined description of the pRBC rupture and of the merozoites release phenomenon. We also emphasize that such a model easily allows for inclusion of anti-malarial treatments acting on only some parasite developmental stages.

Uninfected RBC dynamics. We denote by R r (t), R m (t) and R s (t) respectively the density of reticulocytes, mature RBCs and senescent RBCs at time t.

In the absence of malaria parasites, the evolution of circulating red blood cells is assumed to follow a discrete age maturation system of ordinary differential equations that take the form

     dRr(t) dt = Λ 0 -µ r→m R r (t), dRm(t) dt = µ r→m R r (t) -µ m→s R m (t), dRs(t) dt = µ m→s R m (t) -µ s→d R s (t).
(

) 1
The parameters 1/µ r→m , 1/µ m→s and 1/µ s→d respectively denote the average duration of RBCs in the reticulocyte, mature and senescent age classes while Λ 0 represents the normal value of the RBC production from marrow source (i.e. the production rate of RBC). System (1) can also be found in [START_REF] Mcqueen | Host Control of Malaria Infections: Constraints on Immune and Erythropoeitic Response Kinetics[END_REF]. The parameters of this system are selected from [START_REF] Hetzel | The within-host cellular dynamics of bloodstage malaria: Theoretical and experimental studies[END_REF][START_REF] Mcqueen | Host Control of Malaria Infections: Constraints on Immune and Erythropoeitic Response Kinetics[END_REF] (see Table 1) so that in the absence of parasites, the equilibrium age distribution is given by (R * r ; R * m ; R * s ) = (62.50; 4853; 83.30) × 10 6 cell/ml.

(

This leads to the homeostatic equilibrium concentration of RBC (R * r + R * m + R * s ) around 4.99 × 10 9 cells/ml which is in the range expected for humans.

Parasite dynamics with stage-structured formulation (PDE model).

Here we consider the interaction between free merozoites together with the circulating RBCs. We, respectively, denote by m(t), p(t, a) and G(t) the density of merozoites, parasitized RBC, and mature gametocytes at time t. The variable a denotes the time since the concerned RBC is parasitized (i.e. a2 a1 p(t, a)da corresponds to the density of pRBC at time t which are infected since the time a 1 < a < a 2 ). The system we shall consider reads as:

               p(t, 0) = βm(t) j=r,m,s γ j R j (t), ∂ t p(t, a) + ∂ a p(t, a) = -(µ(a) + d 0 ) p(t, a), ṁ(t) = (1 -α G ) ∞ 0 rµ(a)p(t, a)da -µ m m(t) -βm(t) j=r,m,s γ j R j (t), Ġ(t) = α G ∞ 0 rµ(a)p(t, a)da -µ G G(t). (3)
We briefly sketch the interpretation of the parameters arising in [START_REF] Alonso | A Research Agenda to Underpin Malaria Eradication[END_REF]. Parameters d 0 , µ m and µ G , respectively, denote the natural death rates for uRBC, for free merozoites and for mature gametocytes. Function µ(a) denotes the additional death rate of pRBC due to the parasites at stage a and leading to the rupture. The rupture of pRBC at stage a results in the release of an average number r of merozoites into the blood stream, so that pRBC then produce, at stage a, merozoites at rate rµ(a). Together with this description, the quantity ∞ 0 rµ(a)p(t, a)da corresponds to the number of merozoites produced by pRBC at time t. The parameter β describes the contact rate between uRBC and free merozoites. Parameters γ k with k = r, m, s describe the age preference of parasites' targets. Here we shall be concerned in P. falciparum infection that does not have any preference for RBC so that γ r = γ m = γ s = 1. However when considering P. vivax infection one has γ r = 1 and γ m = γ s = 0, so that target RBCs mostly consist in reticulocyltes while when P. malariae infection is concerned then target RBCs are mostly senescent cells, that is γ r = γ m = 0 and γ s = 1 [START_REF] Paul | The evolutionary ecology of Plasmodium[END_REF]. The parameter α G represents the proportion of merozoites from a bursting asexual schizonts that will enter the gametocyte compartment, i.e., are "committed" to the gametocyte developmental pathway. For simplicity, we have ignored the age structure of gameteocytes and consider G as capturing mature, measurable gametocytes.

Parasite dynamics with K-compartments ODE formulation (ODE model).

For the ODE model formulation, we consider K stages for the pRBC before rupture and set p = (p 1 , p 2 , p 3 , • • • , p K ), such that p j (t) denotes the concentration of pRBC at time t. Then, setting ż = dz dt the ODE model writes

                         ṗ1 (t) = βm(t)
j=r,m,s

γ j R j (t) -(µ 1 + d 1 ) p 1 (t), ṗ2 (t) = µ 1 p 1 (t) -(µ 2 + d 2 ) p 2 (t), . . . ṗK (t) = µ K-1 p K-1 (t) -(µ K + d K ) p K (t), ṁ(t) = (1 -α G )rµ K p K (t) -µ m m(t) -βm(t)
j=r,m,s

γ j R j (t), Ġ(t) = α G rµ K p K (t) -µ G G(t), (4)
wherein 1/µ i the duration of the i-stage and d i the death rate of pRBC. The number of stages K is variable and other parameters and state variables are the same as for the PDE model.

The immune responses. Following [START_REF] Dietz | Mathematical model of the first wave of plasmodium falciparum asexual parasitemia in non-immune and vaccinated individuals[END_REF], here we consider two immune responses (IRs) controlling the growth of the parasite population: (i) an innate IR S I (t) at time t representing the effect of the pro-inflammatory cytokine cascade and (ii) an adaptive IR S A (t) at time t. The effect of the innate IR is a function of the present parasite (merozoite) density that takes the form

S I (t) = m(t) m(t) + S * I , (5)
where S * I is the critical parasite density at which the current multiplication factor is reduced by 50%.

The adaptive IR is a function of the cumulative parasite density; this function is determined by two host-specific parameters and one constant: (1) S * A is the critical cumulative parasite density at which the current multiplication factor is reduced by 50%;

(2) ∆ 0 = 16 days is the average delay required by adaptive IR to become effective [START_REF] Dietz | Mathematical model of the first wave of plasmodium falciparum asexual parasitemia in non-immune and vaccinated individuals[END_REF], i.e., for time t before ∆ 0 the cumulative density is set to zero (the adaptive IR has no effect and S A (t) = 0 for t ≤ ∆ 0 ) and (3) ∆ 1 = 8 days is the delay that determines the last term in the cumulative density for times t ≥ ∆ 0 , i.e.,

S A (t) =              t ∆0 m(s)ds t ∆0 m(s)ds + S * A , ∆ 0 ≤ t < ∆ 0 + ∆ 1 ; ∆0+∆1 ∆0 m(s)ds ∆0+∆1 ∆0 m(s)ds + S * A , t ≥ ∆ 0 + ∆ 1 . (6)
Thus, including these two IR effects, the dynamics of asexual parasite concentration m(t) should be replaced in Models ( 3) and ( 4) respectively by: [START_REF] Bansal | Immune Responses in Malaria Transmission[END_REF] and

ṁ(t) = (1-α G ) ∞ 0 rµ(a)p(t, a)da-   µ m + β j=r,m,s γ j R j (t) + S A (t)   m(t)-S I (t),
ṁ(t) = (1 -α G )rµ K p K (t) -   µ m + β j=r,m,s γ j R j (t) + S A (t)   m(t) -S I (t). ( 8
)
Initial conditions. For both PDE and ODE models, the initial RBCs are assumed to be at their homeostatic equilibrium distribution in the absence of parasites given by (2

), i.e., R r (0) = R * r ; R m (0) = R * m ; R s (0) = R * s .
The above models are also assumed to be free of pRBCs at the initial time, and the initial density of malaria parasites is such that m(0) = m 0 , with m 0 a positive constant. These initial conditions are summarized in Table 2.

Parasitemia. The output of both mathematical models can be used to recover the time course of parasitemia, defined as the proportion of all infected RBC among the total number of RBC. Using the notation of the model, the parasitemia at time t, denoted by P (t) is calculated as follows

P (t) = ∞ 0 p(t, a)da ∞ 0 p(t, a)da + j=r,m,s R j (t) PDE model or K l=1 p l (t) K l=1 p l (t) + j=r,m,s R j (t) ODE model . ( 9
)
3 Results

Development of parasites within RBCs and rupture of pRBCs

An important characteristic of P. falciparum is the development of parasites within RBCs. The parasite within a RBC then takes an average of 48 hours to mature and release free merozoites. With a sequential progression through K stages of parasite maturity before the rupture of the pRBC, the ODE model quantifies the average parasite's development period by

1 µ 1 + • • • + 1 µ K = 48 hours, (10)
where 1/µ i is the waiting time across the i-stage of maturity. Indeed, the probability of the pRBC of being in the i-stage after a hours of infection is then given by D i (a) = P(τ i > a), where τ i denotes the duration within the i-th compartment. One assumption of the K-compartments ODE model is that variables τ i ,s are independents and exponentially distributed with parameter µ i (i.e. D i (a) = e -µia , without taking into account other mechanisms such as natural mortality), such that (10) is satisfied. Thus,

K i=1 E[τ i ] = K i=1 1/µ i = 48 hours.
With the PDE model, the development of parasites within RBCs is characterized by the rupture function µ(a), which takes the form

µ(a) = 0 if a < 48 hours, µ if a ≥ 48 hours,
where a is the age of the pRBC and µ is a positive parameter. With such formulation, the overall average development period is ≈ 48 hours as for the ODE model. Indeed, let D(a) = exp -a 0 µ(σ)dσ the probability that a pRBC remains parasitized after a hours (without taking into account other mechanisms such as the natural mortality). Then, the average parasite's development period is

∞ 0 D(a)da = 48 + 1 µ .
Here we fix, for e.g. µ = 10, such that

∞ 0 D(a)da = 48 + 1 µ ≈ 48.
The value of µ is therefore not strictly significant as soon as the last approximation holds.

Consequently, the PDE model formulation allows to continuously track the development of parasites within RBCs and then to have a refined description of the pRBC rupture followed by the merozoites release phenomenon. By contrast, besides the issue of determining the maturation probability {D i } i=1,••• ,K , such a continuous process is quite difficult to capture with the ODE model with K repeated stages.

One option for defining the maturation probability {D i } i=1,••• ,K can be obtained through a "linear chain trick" formulation. Indeed, by assuming that the duration of each repeated compartments is the same (i.e., µ i = µ 0 , for all i = 1, • • • , K), by (10), we then have D i (a) = e -aK/48 , for all i = 1, • • • , K. The total duration before rupture becomes T K = K i=1 D i . Here D i ,s are independent and identically distributed with exponential law of parameter µ 0 = K/48. Hence, T K follows a Gamma distribution Γ(K, µ -1 0 ) = Γ K, 48 K . We recover that the mean value of T K is 48h and also that its variance is given by var(T K ) = 48 2 /K. The latter quantity tends to 0 as K → ∞ meaning that T K → 48 as K → ∞. As a consequence of the above computations, when K is very large the probability that a pRBC remains parasitized after a hours, is approximately given by P(T K ≥ a) ≈ 1 if a ≤ 48 else 0, that is close to D(a) when µ is large (Fig. 2). However, the main problem with such formulation is that a very large number of repeated compartments is necessary to capture a quite realistic development process of parasites within a RBC. For instance, the probability that a pRBC is still parasitized after 48 hours of infection is approximately 37%, 45%, 48%, respectively for K = 1, 10, 50 (Fig. 2). 0.0 0.2 0.4 0.6 0.8 1.0 q q q q q q q q q q q µ=1 µ=3 µ=10 k=1 k=50 k=200 0 48

Time since sequestration (h)

Within-host reproduction number

The basic reproduction number, usually denoted as R 0 , is defined as the total number of parasites arising from one newly pRBC introduced into an uninfected host. It can be used to study the spread of the malaria parasite in an uninfected host, and the parasite will spread if R 0 > 1. The R 0 of the K-compartments ODE and PDE models write

R 0 =                β µ m + β j=r,m,s R * j (1 -α G )r K i=1 µ i µ i + d i j=r,m,s R * j , ODE β µ m + β j=r,m,s R * j (1 -α G )r µ µ + d 0 e -48×d0 j=r,m,s R * j , PDE. (11)
We refer to [START_REF] Demasse | An Age-Structured Within-Host Model for Multistrain Malaria Infections[END_REF][START_REF] Iggidr | Global Analysis of New Malaria Intrahost Models with a Competitive Exclusion Principle[END_REF] for details on the derivation of [START_REF] Cao | Modeling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection[END_REF].

While the probability for merozoites production for each infection cycle is always 1 for the ODE model, such probability is µ µ+d0 for the PDE model. One reason for this is that the K-compartments ODE model is essentially Markovian, i.e. 'memoryless'. With the K-compartments ODE model, a RBC that has been parasitized for 40 hours has the same probability of producing merozoites as e.g. a RBC parasitized less than a hour ago. However, parameter µ can then be chosen such that these probabilities are close to unity. For instance, for the PDE model, µ µ+d0 ≈ 1 as soon as µ is sufficiently large compared to d 0 . Therefore, with the value of µ = 10 introduced in the previous section, we have µ µ+d0 ≈ 0.99.

Finally, one of the main differences between the R 0 expressions of both models is the probability at which pRBCs survive the 48 hours of the parasite's development period for each infection cycle. While such probability is quantified by the term e -48×d0 for the PDE model [START_REF] Cao | Modeling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection[END_REF], it is K i=1 µi µi+di for the ODE model [START_REF] Cao | Modeling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection[END_REF]. However, in some configurations of parameters µ i ,s and for K sufficiently large we can have K i=1 µi µi+di ≈ e -48×d0 . For instance, by assuming that the duration of each repeated compartments is the same for the ODE model (i.e., µ i = µ 0 , for all i = 1, • • • , K), equality [START_REF] Bousema | Determinants of Malaria Transmission at the Population Level[END_REF] gives K/µ 0 = 48 such that

K i=1 µ i µ i + d i = (1 + d 0 /µ 0 ) -K = exp -K ln 1 + 48d 0 K ≈ exp(-48d 0 ) if K 1.

Fitting the model parameters with data

The model presented above is solved numerically by using finite volume numerical schemes (implemented with the MatLab Programming Language). The model is then fitted to the data for the time course of gametocytes of the patients. To fit our model, let us observe that most of the parameters are estimated from the literature [START_REF] Anderson | Non-linear phenomena in host-parasite interactions[END_REF][START_REF] Eichner | Genesis, sequestration and survival of Plasmodium falciparum gametocytes: Parameter estimates from fitting a model to malariatherapy data[END_REF][START_REF] Hetzel | The within-host cellular dynamics of bloodstage malaria: Theoretical and experimental studies[END_REF][START_REF] Mcqueen | Host Control of Malaria Infections: Constraints on Immune and Erythropoeitic Response Kinetics[END_REF][START_REF] Mcqueen | Host immune constraints on malaria transmission: Insights from population biology of within-host parasites[END_REF].

Table 1 provides the values we shall use for the fixed parameters. Three parameters need to be estimated from the data, these are: the proportion of parasitized cells that produce asexual merozoites (α G ), the merozoite initial density (m 0 ), and the duration of sexual stage (1/µ G ). Additionally to the above three parameters, the number K of compartments is also estimated for the ODE model. These parameters are adjusted from the data for each patient by using a least square method. Basically, we find the values which minimize the difference between the ODE model prediction gametocyte density and the observed data by using MatLab nonlinear least-squares solver lsqcurvefit. Those optimal parameters for the ODE model are then used to run the PDE model. The superposition of the data and gametocyte density output of the mathematical models are presented in Fig. 3, while the estimated parameter values for each patient are given in Table S1.

Comparison of ODE and PDE model outputs

We have presented two modelling frameworks to properly model the within-host infection of malaria. Within this context, we compare a classical model based on ordinary differential equations (ODE) with a model based on partial differential equations (PDE). Our first observation is on the parameterisation of both models. More precisely, a good description of the rupture of pRBC requires at least one additional parameter K for the repeated compartments, see (3) versus (4). Such parameter K is necessary to capture the delay in the production or quantification of gametocytes imposed by the development of parasites within RBCs for each infection cycle. This delay in gametocytes production is nicely highlighted by the PDE model formulation (Fig. 3). Through a "linear chain trick" formulation, it is then possible to find the number K of repeated compartments such that the ODE model can slow down the production dynamics of gametocytes. Indeed, by assuming that the duration of each repeated compartment is the same, we can suitably find the parameter K such that both models perform quite similarly in terms of goodness-of-fit for gametocyte production. Overall, for our dataset, we find between 47 to 68 compartments for a goodness-of-fit of the ODE model (Fig. 3). While both models perform quite similarly in terms of goodness-of-fit for a suitable value of K, the K-compartments ODE model particularly overestimates parasite densities early on in infections when the number of repeated compartments is not large enough (Fig. 3). Importantly, the number K of compartments for the goodness-of-fit of the ODE model is quite large and can be variable across individuals (Fig. 3). According to the infection dynamics, our comparative results show that the PDE model and the K-compartments ODE model (when K is suitably chosen) reproduced, at least qualitatively, the true dynamics of malaria infection parasitemia. Indeed, although we do not have parasitemia real data for patients considered here, but in a qualitative comparison to some studies (e.g., [START_REF] Childs | Dissecting the determinants of malaria chronicity: Why within-host models struggle to reproduce infection dynamics[END_REF]), the dynamics of both models seem to mimic qualitatively the parasitemia dynamics when the number of compartments for the ODE model is adequately defined (Fig. 4).

Relationship between parasitemia and gametocyte density

Our mathematical model has been fitted to the available data for each patient under consideration, which consists of gametocyte densities over time. We now use the output of this mathematical model to recover the time course of parasitemia, defined as the proportion of all infected RBC among the total number of RBC, see [START_REF] Bousema | Epidemiology and Infectivity of Plasmodium falciparum and Plasmodium vivax Gametocytes in Relation to Malaria Control and Elimination[END_REF]. The time course of parasitemia, P (t), computed from our model are presented in Fig. 5 for each patient together with the fitted gametocyte trajectories. As is observed for each patient, the relationship between these curves exhibits two different regimes. During some period of time [2, T 0 ], the two curves are increasing with rather similar shape up to a time shift (of length 2 days). This means that, in this increasing regime, the gametocyte density at time t depends on the parasitemia at time t -2, a delay which reflects the life cycle of the parasites inside the RBCs. After this period of increasing parasitemia and gametocyte density, namely after time T 0 , both curves are decreasing and the shapes seem to depend upon the specific patient considered. To make these comments more quantitative, we introduce the following formula from an estimation of the gametocyte density G(t) from the parasitemia P (t):

G(t) = k 1 P (t -2) θ1 if 2 ≤ t ≤ T 0 days, k 2 P (t) θ2 if T 0 ≤ t ≤ 30, ( 12
)
where k 1 , k 2 , θ 1 and T 0 are four positive parameters while θ 2 is a negative parameter. Let us mention that this 2 days delay between gametocyte density and parasitemia should not to be confused, for example, with the time to distinguish the mature gametocytes via microscopy. Such time is longer than 2 days and can be well captured by the PDE model or the ODE model when K is adequately chosen (Fig. 3). While there is a precise biological relationship between parasitemia and gametocytes density at some point in the future, here we seek a robust statistical relationship and so the delay need not match with, e.g., gametocyte maturation time.

To determine the unknown parameters k 1 , k 2 , θ 1 , θ 2 and the changing time T 0 for each patient, we perform a least square analysis. More specifically, we adjust these parameters through a logarithmic scale, that is, through the following formula

log 10 G(t) = log 10 k 1 + θ 1 log 10 P (t -2) if 2 ≤ t ≤ T 0 , log 10 k 2 + θ 2 log 10 P (t) if T 0 < t ≤ 30.
To be more precise, our analysis couples the estimate for the time parameter T 0 , at which the above formula is changing formulation, together with two linear regressions on each part of the graph. We find that parameters k 1 , k 2 , θ 1 and T 0 are remarkably robust with respect to individuals while θ 2 depends upon each individual. The results of this analysis as well as the estimated parameters are presented in Fig. S4 and summarized by Table S2 for each patient. The quality of the fit is quantified using the coefficient of determination R 2 (for linear regression). It is computed for each patient and for the two regimes independently. This adjustment metric is computed using the sample of points induced by the time discretization of the partial differential equation model. For our four cases, this coefficient of determination R 2 is approximately 0.99 in the first part of the curve and even closer to one in the second part of the graph.

Coming back to the adjusted parameters described in Table S2, one may observe that the four parameters k 1 , k 2 , θ 1 and T 0 have robust values with respect to patients while the parameter θ 2 depends on the patient. Using the average values of the adjusted parameters on the set of malariatherapy data, we derive the following clinical formula:

G(t) = 3.843 × 10 7 • P (t -2) +1.0304 if 2 < t ≤ T 0 days, 2.981 × 10 9 • P (t) -0.0470 if T 0 ≤ t ≤ 30 days, (13)
with T 0 = 14.5636 ± 0.0064 days. The relative error for ( 13) is such that

∆G G 2 ≤ 2.3884 × 10 -4 + 1.0617 ∆P P 2
, where ∆P P is the relative error on the measurement of parasitemia. Therefore, if ∆P P ≤ 5%, then ∆G G < 5.38%. From practical point of view, notice that formula (13) can be really useful to estimate the gametocyte density from the parasitemia measurement without necessarily using the quite 'complex' mathematical model described in this note.

Discussion

Many models of within-host malaria infection dynamics have been formulated since the pioneer work of Anderson el al. [START_REF] Anderson | Non-linear phenomena in host-parasite interactions[END_REF] in 1989. These models are based on dynamical systems, with standard approaches ranging from ordinary differential equations (ODEs), to delay differential equations (DDEs) or partial differential equations (PDEs). Most ODE model formulations [25, 28-30, 36, 44, 45] assume an exponential process to describe the rate of pRBCs rupture and therefore fail to capture realistic lifetimes of the pRBCs. This issue is somewhat corrected when the development of parasites within RBCs and rupture of pRBCs are modeled either by a set of K-compartments ODE [START_REF] Gravenor | Reply to: Models for the in-host dynamics of malaria revisited: Errors in some basic models lead to large over-estimates of growth rates[END_REF][START_REF] Iggidr | Global Analysis of New Malaria Intrahost Models with a Competitive Exclusion Principle[END_REF][START_REF] Saralamba | Intrahost modeling of artemisinin resistance in Plasmodium falciparum[END_REF][START_REF] Zaloumis | Assessing the utility of an anti-malarial pharmacokinetic-pharmacodynamic model for aiding drug clinical development[END_REF], or DDEs [START_REF] Cao | Modeling the dynamics of Plasmodium falciparum gametocytes in humans during malaria infection[END_REF][START_REF] Fonseca | Comparison of Mathematical Frameworks for Modeling Erythropoiesis in the Context of Malaria Infection[END_REF][START_REF] Hoshen | Mathematical modelling of the within-host dynamics of Plasmodium falciparum[END_REF][START_REF] Kerlin | Preferential Invasion by Plasmodium Merozoites and the Self-Regulation of Parasite Burden[END_REF]. Other approaches are the use of PDEs to track the infection history of a pRBC [START_REF] Antia | The dynamics of acute malaria infections. I. Effect of the parasite's red blood cell preference[END_REF][START_REF] Cromer | Low red cell production may protect against severe anemia during a malaria infection-Insights from modeling[END_REF][START_REF] Demasse | An Age-Structured Within-Host Model for Multistrain Malaria Infections[END_REF][START_REF] Khoury | Within-host modeling of blood-stage malaria[END_REF]. Using gametocyte production, parasitemia (the proportion of all infected RBC among the total number of RBC) as proxy variables and malariatherapy data, we found that the PDE model and the ODE model perform similarly in terms of goodness-offit when a suitable value of K is chosen (Fig. 3). Some disadvantages of the ODE model are that the number K of compartments required to achieve a good fit is quite large and can be variable across individuals, and the ODE model particularly overestimates parasite densities early on in infections when the number of repeated compartments is not large enough (Fig. 3). Similar comparison holds for the parasitemia dynamics (Fig. 4).

Not least, the K-compartments ODE model (for suitably chosen K) and the PDE model highlight a strong qualitative connection between gametocyte density and parasitemia. From a practical point of view, such a relation given by ( 13) can be really useful to estimate the gametocyte density from the parasitemia measurement without necessarily using the quite 'complex' mathematical model described in this note.

Here, immune-mediated parasite killing is only considered for merozoites. This choice for immunity targeting merozoites, rather than parasitized red blood cells is mostly because it is a lot easier with our PDE model formulation, particularly in terms of parameterization. However, in some studies, e.g. [START_REF] Dietz | Mathematical model of the first wave of plasmodium falciparum asexual parasitemia in non-immune and vaccinated individuals[END_REF], parasite levels are not distinguished by merozoite and parasitized red blood cells, such that immunity is acting against merozoite and parasitized red blood cells. Also note that, while there is evidence that the mature live gametocytes evade the immune response clearance pathway, immune-mediated parasite killing of immature gametocytes has been shown in some studies [START_REF] Bansal | Immune Responses in Malaria Transmission[END_REF]. Finally, taking the best fit parameters and then altering immunity parameters between their minimum, median, and maximum values estimated in [START_REF] Dietz | Mathematical model of the first wave of plasmodium falciparum asexual parasitemia in non-immune and vaccinated individuals[END_REF] has very little impact on the model outputs, namely, parasitemia and gametocytes density (figures not shown). However, this can be explained by the fact that (i) merozoites are only short lived and (ii) the relatively short-term validation of the model presented here.

Reducing infections in mosquitoes-vectors of Plasmodium parasites-is a crucial component of global efforts to control and eliminate malaria [START_REF] Alonso | A Research Agenda to Underpin Malaria Eradication[END_REF]. Because a strong correla-tion exits between the gametocyte density within a host and infectivity of mosquitoes [START_REF] Bousema | Epidemiology and Infectivity of Plasmodium falciparum and Plasmodium vivax Gametocytes in Relation to Malaria Control and Elimination[END_REF][START_REF] Churcher | Predicting mosquito infection from Plasmodium falciparum gametocyte density and estimating the reservoir of infection[END_REF][START_REF] Collins | A retrospective examination of mosquito infection on humans infected with Plasmodium falciparum[END_REF][START_REF] Graves | Measurement of malarial infectivity of human populations to mosquitoes in the Madang area, Papua New Guinea[END_REF], progress towards this goal would be bolstered by quantifying gametocytes and identifying highly infectious hosts [START_REF] Bousema | Determinants of Malaria Transmission at the Population Level[END_REF][START_REF] Stone | Assessing the infectious reservoir of falciparum malaria: Past and future[END_REF][START_REF] Transmission | malERA: An updated research agenda for characterising the reservoir and measuring transmission in malaria elimination and eradication[END_REF]. On the other hand, parasitemia is easily quantified by light microscopy and therefore is more technically accessible, particularly in regions where malaria is endemic. Therefore, quantifying the relationship between the gametocyte density and parasitemia is of great interest to define more simple tools for the prediction of mosquito infection. The results presented in this note provide one such tool.

From a public health or population dynamics point of view, the time course of the disease at the between-host level is strongly related to the basic reproduction number (also denoted here by R 0 ). At the between-host level, the R 0 is defined as the number of secondary infections from a single infected individual introduced in a fully susceptible population. This important metric can be estimated from real data but also using mathematical models. The simplest (deterministic) mathematical model reads as the Ross system of equations from which one can compute this threshold number R 0 as follows:

Between-host R 0 = R V H 0 × R HV 0 with R V H 0 = abd M , R HV 0 = macd H .
Note that the above R 0 is for between-host malaria dynamics and is not for the withinhost models presented here. In the above formula of R 0 , m represents the number of mosquitoes per person, a denotes the mosquito biting rate, b and c denote the per bite transmission probability respectively from mosquito to human and from human to mosquito, while d H and d M correspond respectively to the human recovery and the mosquito death rates. Although more ingredients can been included into the mathematical model, leading to different formulations for R 0 , the above expression contains the main important parameters [START_REF] Macdonald | The Measurement of Malaria Transmission[END_REF][START_REF] Mandal | Mathematical models of malaria -a review[END_REF][START_REF] Ruan | On the Delayed Ross-Macdonald Model for Malaria Transmission[END_REF][START_REF] Smith | Ross, macdonald, and a theory for the dynamics and control of mosquito-transmitted pathogens[END_REF]. Parameters b and c serve as the link between within-and between-host dynamics, since the transmission rates from (to) a host will depend on the dynamics of what is happening within that host (vector). Furthermore, there is a clear relationship between gametocyte density (G) and the transmission probability per bite from human to mosquito (c) [START_REF] Bousema | Mosquito Feeding Assays to Determine the Infectiousness of Naturally Infected Plasmodium falciparum Gametocyte Carriers[END_REF][START_REF] Bousema | Epidemiology and Infectivity of Plasmodium falciparum and Plasmodium vivax Gametocytes in Relation to Malaria Control and Elimination[END_REF][START_REF] Churcher | Predicting mosquito infection from Plasmodium falciparum gametocyte density and estimating the reservoir of infection[END_REF][START_REF] Collins | A retrospective examination of mosquito infection on humans infected with Plasmodium falciparum[END_REF][START_REF] Drakeley | Host haematological factors influencing the transmission of Plasmodium falciparum gametocytes to Anopheles gambiae s.s. mosquitoes[END_REF][START_REF] Johnston | Malaria's Missing Number: Calculating the Human Component of R0 by a Within-Host Mechanistic Model of Plasmodium falciparum Infection and Transmission[END_REF][START_REF] Stepniewska | Plasmodium falciparum gametocyte dynamics in areas of different malaria endemicity[END_REF]. From a practical point of view, the parameter c is difficult to estimate in a relevant way. Indeed, an efficient measurement of c requires a good measure of the gametocyte density which is quite difficult to obtain in practice. Indeed, while we can get either gametocytes or parasitemia from microscopy, but with gametocytes tending to be at lower densities (sometimes orders of magnitude lower), there is a detectability issue. Thus, a simple way to estimate the gametocyte density will help to infer the parameter c. Thanks to formula [START_REF] Chernin | The Malariatherapy of Neurosyphilis[END_REF] proposed here, we then have robust relationships between parasitemia (easier to measure) and gametocyte density, at least during the first days of infection (approximately the first two weeks).

Overall, the proposed model for the dynamics of gametocytes is probably valid for the first asexual wave which last approximately 40 days for each patient in this malariatherapy dataset. This relatively short-term validation is enough for the aim of the current study. However, for the long-term gametocyte dynamics, we need to bring more complexity in the model proposed here. Indeed, the conversion probability of asexual parasites to circulating gametocytes (α G ) should be considered to vary among successive waves of asexual parasitaemia to tackle such issue of the long-term gametocyte dynamics [START_REF] Eichner | Genesis, sequestration and survival of Plasmodium falciparum gametocytes: Parameter estimates from fitting a model to malariatherapy data[END_REF], par-ticularly since smear-positive asymptomatic malaria infections detectable by microscopy are an important gametocytes reservoir and often persist for months [START_REF] Lin | The role of submicroscopic parasitemia in malaria transmission: What is the evidence?[END_REF].

The robustness of the model proposed here, especially the formula linking parasitaemia and gametocyte density, is only guaranteed during the first two weeks after infection. Beyond this time, this estimate is highly variable from one patient to another. This variability is explained, at least in part, by the variability of the duration of sexual stage (1/µ G ) which governs the decrease in gametocyte density (Table S1). One interpretation of the variation in this parameter is that there exists variation in how well individuals clear gametocytes or kill them through immune responses [START_REF] Doolan | Acquired Immunity to Malaria[END_REF][START_REF] Lin | The role of submicroscopic parasitemia in malaria transmission: What is the evidence?[END_REF]. Less variability is expected at the beginning of infection, where the whole system is less constrained by immunity. This is likely to be true for the malariatherapy patient data presented here, since hosts were initially naive. However in high transmission settings, acquired immunity-particularly in older hosts-may obscure the relationship between early gametocyte and asexual parasite densities our work has revealed. Furthermore, fine scale, longitudinal data could assess the applicability of this relationship across settings and age groups, although such data is understandably difficult to obtain.

Finally, while our work reveals a simple tool for linking aspects of the early dynamics of malaria infections, it also offers specific suggestions for how best to mathematically describe those infection dynamics more broadly. Both PDE and K-compartments ODE models have been adopted to capture the subtleties of malaria parasite life cycles in bloodstage infections [START_REF] Khoury | Within-host modeling of blood-stage malaria[END_REF]. Our work provides more understanding among using the PDE model or the ODE model for the within-host human malaria dynamics. Population of senescent RBC 83 × 10 6 RBC.ml -1 p(0, .)

Population of pRBC for the PDE model 0 cells.ml -1 p j Population of pRBC for the ODE màdel 0 cells.ml -1 G(0)

Population of mature gametocyte 0 cells.ml -1 m(0)

Population of malaria parasites variable B Supplementary tables -0.0203 ± 0.0003 -0.0188 ± 0.0002 -0.0170 ± 0.0002 -0.0336 ± 0.0003 T 0 (days) 20.0083 ± 0.0075 14.6917 ± 0.0083 14.2625 ± 0.0083 14.5458 ± 0.0025

Figure 1 :

 1 Figure 1: (S 1 ) The RBC development chain, (S 2 ) the parasite development chain. T D = average duration (± one standard deviation) spent in an RBC age class given in[START_REF] Mcqueen | Host immune constraints on malaria transmission: Insights from population biology of within-host parasites[END_REF], Λ 0 is the RBC production rate from the marrow source. In our model, the parameter 1/µ r→m (resp. 1/µ m→s , 1/µ s→d ) is the time spent in RBC reticulocyte (resp. mature, senescent) class. A continuous parameter a denotes the time since the concerned RBC is parasitized: ring stage (0 < a < 26 hours), trophozoite (26 < a < 38 hours) and schizont (38 < a < 48 hours). In the case of P. falciparum infection, one has (γ r = γ m = γ s = 1) while for P. vivax one has (γ r = 1, γ m = γ s = 0) and for P. malariae (γ r = γ m = 0, γ s = 1)[START_REF] Paul | The evolutionary ecology of Plasmodium[END_REF].

Figure 2 :

 2 Figure 2: The probability that a pRBC is still parasitized after a hours of infection for both PDE and ODE models. Different values of K correspond to the maturation probability P(T K ≥ a) for the ODE model, while values of µ are for the function D(a) of the PDE model.

Figure 3 :

 3 Figure 3: Comparison with data and mathematical model output for gametocyte density for patients G161, G104, S1050 and G299. The ODE model is illustrated for three values of K (K = 1, K = 100, and an intermediate value corresponding to the optimal K for each patient). Comparison for other patients is provided by Figure S1.

Figure 4 :

 4 Figure 4: The time course of parasitemia (in percentage) for patients G161, G104, S1050 and G299. The ODE model is illustrated for three values of K (K = 1, K = 100, and an intermediate value corresponding to the optimal K for each patient). Other patients are provided by Figure S2.

Figure 5 :

 5 Figure 5: The time course of parasitemia (in percentage) and gametocyte density computed from the PDE model for patients G161, G104, S1050 and G299. The ODE model is illustrated only for the optimal K of each patient. Other patients are provided by Figure S3.

Figure S2 :

 S2 Figure S2: The time evolution of parasitemia (in percentage) density for other patients. The ODE model is illustrated for three values of K (K = 1, K = 100, and an intermediate value corresponding to the optimal K for each patient).

Figure S3 :

 S3 Figure S3: The time evolution of parasitemia (in percentage) and gametocyte density curves computed from the PDE model for other patients. The ODE model is illustrated only for the optimal K of each patient.

Figure S4 :

 S4 Figure S4: Comparison between model prediction and the linear regression based on formula[START_REF] Chernin | The Malariatherapy of Neurosyphilis[END_REF]

Table 1 :

 1 Fixed model parameters

	Parameters Description (unit)	Values	References
	Λ 0	Production rate of RBC (RBC/h/ml)	1.73 × 10 6	[4, 42]
	1/µ r→m	Duration of the RBC reticulocyte stage (h)	36	[41]
	1/µ m→s	Duration of the RBC mature stage (day)	116.5	[41]
	1/µ s→d	Duration of the RBC senescent stage (h)	48	[41]
	β	Infection rate of uRBC (RBC/ml/day)	6.27 × 10 -10 [4]
	d 0	Natural death rate of uRBC (RBC.day -1 )	0.00833	[4]
	µ m	Decay rates of malaria parasites (RBC.day -1 )	48	[29]
	r	Merozoites multiplication factor (dimensionless)	16	[4]
	α G	Proportion of sexual merozoites (dimensionless)	0.05	[42]
	S * I S * A	Innate IR density for 50% of parasite killing (cells.µl -1 ) Adaptive IR density for 50% of parasite killing (cells.µl -1 ) 20.4 2,755	[18] [18]
	∆ 0	Delay required by adaptive IR to be effective (day)	16	[18]

Table 2 :

 2 Initial values for the model

	Variables Description	Initial Values
	R r (0)	Population of reticulocytes RBC	62 × 10 6 RBC.ml -1
	R m (0)	Population of mature RBC	4.85 × 10 9 RBC.ml -1
	R s (0)

Table S1 :

 S1 Patient-specific parameters estimated from the data

		G54	G221 S1050 S1204 S1300 G104
	µ G (×10 -3 ) 2.27 0.4632	1.99	1.5	1.07	0.667
	α G (×10 -8 ) 13.02 1.206	1.161 1.289	13.47	13.02
	m 0 (×10 7 )	2.5	1	1	2.97	2.5	6
		G139 G161	G305 G338	G299 S1326
	µ G (×10 -3 ) 1.25	1.25	0.526 0.667	0.80	1.43
	α G (×10 -8 ) 0.134 118.57 12.84 1.321 1.1682 1.382
	m 0 (×10 7 )	2.5	5	0.7	2	0.05	1.1

Table S2 :

 S2 Estimated parameters for the shape (12) -0.0141 ± 0.0004 -0.0477 ± 0.0005 -0.0336 ± 0.0005 -0.0280 ± 0.0003 T 0 (days) 14.7125 ± 0.0075 13.2917 ± 0.0083 14.4958 ± 0.0083 13.4666 ± 0.-0.0336 ± 0.0003 -0.0170 ± 0.0002 -0.0336 ± 0.0003 -0.0214 ± 0.0003 T 0 (days) 14.5375 ± 0.0075 13.2625 ± 0.0083 14.5625 ± 0.0083 13.3875 ± 0.

		G221	G54	S1050	S1204
	log 10 (k 1 )	2.9586 ± 0.0037	2.9326 ± 0.0039	2.9520 ± 0.0037	2.9502 ± 0.0037
	θ 1	1.0321 ± 0.0028	1.0312 ± 0.0030	1.0315 ± 0.0029	1.0315 ± 0.0028
	log 10 (k 2 )	2.9437 ± 0.0002	3.8802 ± 0.0005	2.9101 ± 0.0004	2.9148 ± 0.0003
	θ 2
					0025
		S1300	G104	G139	G161
	log 10 (k 1 )	2.9861 ± 0.0037	2.9821 ± 0.0041	2.8024 ± 0.0037	2.9272 ± 0.0039
	θ 1	1.0315 ± 0.0028	1.0304 ± 0.0031	1.0315 ± 0.0028	1.0312 ± 0.0029
	log 10 (k 2 )	3.9442 ± 0.0003	3.9663 ± 0.0002	1.9360 ± 0.0003	4.8774 ± 0.0002
	θ 2
					0025
		G299	G305	G338	S1326
	log 10 (k 1 )	2.9712 ± 0.0034	2.9854 ± 0.0037	2.9897 ± 0.0041	2.9893 ± 0.0037
	θ 1	1.0301 ± 0.0027	1.0315 ± 0.0028	1.0304 ± 0.0031	1.0315 ± 0.0028
	log 10 (k 2 )	2.9063 ± 0.0003	3.9607 ± 0.0001	2.9739 ± 0.0002	2.9474 ± 0.0003
	θ 2

Acknowledgements

Authors thank Samuel Alizon, Laurence Ayong, Carole Eboumbou and Christophe Rogier for comments and suggestions to improve the manuscript.

Code availability

The code (with the MatLab Programming Language) used to simulate the model can be accessed through the Zenodo platform at https://doi.org/10.5281/zenodo.5526271

A Supplementary figures