B Borderie

S Piantelli

M F Rivet

Ad R Raduta

E Bonnet

R Bougault

A Chbihi

E Galichet

D Guinet

Ph Lautesse

N Le Neindre

O Lopez

M Parlog

E Rosato

R Roy

G Spadaccini

E Vient

M Vigilante

The behaviour of constrained caloric curves as ultimate signature of a phase transition for hot nuclei

The behaviour of constrained caloric curves as ultimate

signature of a phase transition for hot nuclei

Introduction

One of the most important challenges of heavy-ion collisions at intermediate energies is the identification and characterization of the nuclear liquid-gas phase transition for hot nuclei, which was earlier theoretically predicted for nuclear matter [1]. In the last fifteen years a big effort to accumulate experimental indications of the phase transition has been made. Statistical mechanics for finite systems appeared as a key issue to progress by proposing new first-order phase transition signatures related to thermodynamic anomalies like negative microcanonical heat capacity and bimodality of an order parameter [START_REF]Dynamics and Thermodynamics with Nuclear Degrees of Freedom[END_REF]3,4]. Correlated temperature and excitation energy measurements, commonly termed caloric curves, were the first studied possible signatures of phase transition. However in spite of the observation of a plateau in some caloric curves, no decisive conclusion could be extracted [1,5,6,7]. The reason is that experimentally it is not possible to explore the caloric curves at constant pressure or constant average volume, which is required for an unambiguous phase transition signature. Indeed, theoretical studies show that if many different caloric curves can be generated depending on the path followed in the thermodynamical landscape, constrained caloric curves exhibit different behaviours in presence of a first order phase transition: a monotonous evolution at constant average volume and a back bending of curves at constrained pressures [8,9]. With the help of a simulation able to correctly reproduce most of the experimental observables measured for hot nuclei formed in central collisions (quasi-fused systems, QF, from 129 Xe+ nat Sn, 32-50 AMeV), event by event properties at freeze-out were restored and used to build constrained caloric curves. The definition of pressure in the microcanonical ensemble is presented in section 2. Then, in section 3, simulations to recover freeze-out properties of multifragmentation events [10,11] are briefly described and constrained caloric curves deduced are discussed. Section 4 is dedicated to the description of the thermometer finally chosen and to comparisons of results relative to a caloric curve without any constraints. Constrained caloric curves with temperatures including quantum fluctuations are presented in section 5. Section 6 is devoted to a discussion of the various results. Conclusions are given in section 7.

Pressure in the microcanonical ensemble

Let us consider a gas of weakly interacting fragments (i.e. they interact only by Coulomb and excluded volume), which corresponds to the freeze-out configuration. Within a microcanonical ensemble, the statistical width of a configuration C, defined by the mass, charge and internal excitation energy of each of the constituting M C fragments, writes

W C (A, Z, E, V ) = 1 M C ! χV M C M C n=1 ρ n (ǫ n ) h 3 (mA n ) 3/2 2π Γ(3/2(M C -2)) 1 (detI) (2πK) 3/2M C -4 (mA) 3/2 , ( 1
)
where I is the moment of inertia, K is the thermal kinetic energy, V is the freeze-out volume and χV M C stands for the free volume or, equivalently, accounts for inter-fragment interaction in the hard-core idealization.

The microcanonical equations of state are

T = ∂S ∂E -1 | V,A , P/T = ∂S ∂V | E,A , -µ/T = ∂S ∂A | E,V . (2)
Taking now into account that

S = ln Z = ln C W C and that ∂W C /∂V = (M C /V ) W C , it comes out that P/T = ∂S ∂V = 1 C W C C ∂W C ∂V = 1 V C M C W C C W C = < M C > V . (3)
The microcanonical temperature is also easily deduced from its statistical definition [12]:

T = ∂S ∂E -1 = ( 1 C W C C W C (3/2M C -5/2)/K) -1 = < (3/2M C -5/2)/K > -1 . ( 4
)
where K is the total thermal kinetic energy of the system at freeze-out. As M C , the total multiplicity at freeze-out, is large, the pressure P can be well approximated by Figure 1. Comparison between the experimental velocity spectra (full points) of fragments of a given charge and the simulated ones (histograms). Each row refers to a different fragment charge: starting from the top Z=6, Z=11, Z=18 and Z=27. Each column refers to a different beam energy: starting from the left 32, 39, 45 and 50 AMeV. From [11].

P = 2 3 < K > V .

Event by event freeze-out properties

Starting from all the available asymptotic experimental information (charged particle energy spectra, average and standard deviation of fragment velocity spectra and calorimetry) of selected QF sources produced in central 129 Xe+ nat Sn collisions which undergo multifragmentation, a simulation was performed to reconstruct freeze-out properties event by event [10,11]. The method requires data with a very high degree of completeness (total detected charge ≥93% of the total charge of the system), which is crucial for a good estimate of Coulomb energy. QF sources are reconstructed, event by event, in the reaction centre of mass, from all the fragments and twice the charged particles emitted in the range 60-120 • in order to exclude a major part of pre-equilibrium emission. Dressed excited fragments and particles at freeze-out are described by spheres at normal density. Then the excited fragments subsequently deexcite while flying away. Four free parameters are used to recover the data at each incident energy: the percentage of measured particles which were evaporated from primary fragments, the collective radial energy, a minimum distance between the surface of products at freeze-out and a limiting temperature for fragments. The limiting temperature, related to the vanishing of level density for fragments [13], was mandatory to reproduce the observed widths of fragment velocities. Indeed, Coulomb repulsion plus collective energy plus thermal kinetic energy (directed at random) plus spreading due to fragment decays are responsible for about 60-70% of the observed widths. By introducing a limiting temperature for fragments, thermal kinetic increases, due to energy conservation, which produces the missing percentage for the widths of final velocity distributions. The agreement between experimental and simulated velocity/energy spectra for fragments and for the different beam energies is quite remarkable (see figure 1). Relative velocities between fragment pairs were also compared through reduced relative velocity correlation functions [14,15] (see figure 2). Again a good agreement is obtained between experimental data and simulations, which indicates that the retained method (freeze-out topology built up at random) and parameters are sufficiently relevant to correctly describe the freeze-out configurations, including volumes. Finally one can also note that the agreement between experimental and simulated energy spectra for protons and alpha-particles (see an example in figure 3) is not very good at high energy. This comes from the fact that we have chosen (to limit the number of parameters of the simulation) a single value for the percentage of all measured particles which were evaporated from primary fragments. We shall see in the following section how to correct for temperature measurements derived from protons.

From the simulations it is then possible to recover, event by event, the different quantities needed to build constrained caloric curves, namely the thermal excitation energy of QF hot . Centre of mass energy spectra for protons, deuterons, tritons, 3 He, 4 He and 7 Li for the reaction at 39 AMeV. Open points are the simulated data, while full points are the experimental ones. From [11].

nuclei, E * (total excitation minus collective energy) with an estimated systematic error of around 1 AMeV, the kinetic temperature T kin at freeze-out, the freeze-out volume V (see envelopes of figure 8 from [11]) and the total thermal kinetic energy at freeze-out K. In simulations, Maxwell-Boltzmann distributions are used to reproduce the thermal kinetic properties at freeze-out and consequently the deduced temperatures, T kin , are classical. Constrained caloric curves, which correspond to correlated values of E * and T kin have been derived for QF hot nuclei with Z restricted to the range 80-100, which corresponds to the A domain 194-238, in order to reduce any possible effect of mass variation on caloric curves [7] (see figure 4). Curves for internal fragment temperatures are also shown in the figure. For different average freeze-out volumes expressed as a function of V 0 , the volume of the QF nuclei at normal density, a monotonous behaviour of caloric curves is observed as theoretically expected. The caloric curves when pressure has been constrained exhibit a backbending and moreover the qualitative evolution of curves with increasing pressure exactly corresponds to what is theoretically predicted with a microcanonical lattice gas model [8]. However by extrapolating at higher pressures, one could infer a critical temperature around 20 MeV. Such a value is within the range calculated for infinite nuclear matter whereas a lower value is expected for nuclei in relation with surface and Coulomb effects [1]. At that point, one may wonder if T kin is a relevant thermometer. We shall see, in what follows, that the final choice was to use, for protons thermally emitted at freeze-out, a new thermometer recently proposed for which quantum effects can be included.

E* (A MeV)

Temperatures from proton transverse momentum fluctuations

Very recently a new method for measuring the temperature of hot nuclei was proposed [START_REF] Wuenschel | [END_REF]18]. It is based on momentum fluctuations of emitted particles like protons in the centre of mass frame of the fragmenting nuclei. On the classical side, assuming a Maxwell-Boltzmann distribution of the momentum yields, the temperature T is deduced from the quadrupole momentum fluctuations defined in a direction transverse to the beam axis:

σ 2 = < Q 2 xy > -< Q xy > 2 = 4m 2 T 2 with Q xy = P 2
x -P 2 y ; m and P are the mass and linear momentum of emitted particles. Taking into account the quantum nature of particles, a correction F QC related to a Fermi-Dirac distribution was also proposed [18,19]. In that case σ 2 = 4m 2 T 2 F QC where F QC = 0.2(T /ǫ f ) -1.71 + 1; ǫ f = 36 (ρ/ρ 0 ) 2/3 is the Fermi energy of nuclear matter at density ρ and ρ 0 corresponds to normal density.

Before using the thermometer with protons to build constrained caloric curves, it was important to verify several things. With the classical simulation (freeze-out and asymptotic proton momenta) , it is possible to test the agreement with the proposed classical thermometer. Moreover the effects of secondary decays on temperature measurements can be estimated. Figure 5 shows different caloric curves without constraints. Note that the selection in Z and A of hot nuclei is the same as in the previous section. It was verified that, within statistical error bars, at a given thermal excitation energy transverse momentum fluctation values are the same for our selection or by selecting only a single (A and Z) hot nucleus. Open diamonds refer to classical temperatures calculated from momentum fluctuations for protons thermally emitted at freeze-out. Within the statistical error bars they perfectly superimposed on T kin values (see figure 2 of [START_REF] Borderie | Proc. Int. Workshop on Multifragmentaton and related topics[END_REF]), which just verifies that Maxwell-Boltzmann distribution was correctly implemented in the simulation. Full squares correspond to classical temperatures calculated from momentum fluctuations for protons after the secondary decay stage. We note that the caloric curve is distorted, which means that it is not possible to use experimental data from protons to measure temperatures. Moreover, in this case, quantum corrections for temperatures can not be made because protons are emitted at different stages of deexcitation with different Fermi energy values. In figure 5 classical temperatures calculated from experimental proton data are also shown (full points). As for temperatures calculated from asymptotic proton data of simulations, a monotonous behaviour of the caloric curve is observed. One also note the differences between the two sets of temperature values, which are related to the fact that, as indicated in the previous section, simulations do not describe experimental proton energy spectra very well. Those temperature differences will be used to correct classical temperatures derived from simulated protons at freeze-out. It finally appears that the only way to extract temperatures from proton transverse momentum fluctuations taking into account quantum effects is to use protons thermally emitted at freeze-out. In that case classical temperature values from simulations must be extracted and corrected and then, quantum corrections applied, which needs Fermi energy values. Those values can be estimated from semi-classical calculations (Xn+Sn at 32 AMeV and Sn+Sn at 50 AMeV) [20,21]: protons thermally emitted at freeze-out at time around 100-120 fm/c after the beginning of collisions come from a low density uniform source. For the two incident energies low densities around ρ ∼ 0.4ρ 0 are calculated which corresponds to ǫ f ∼ 20 MeV. We have introduced a systematic error of ± 0.1ρ 0 for the calculation of ǫ f and consequently a systematic error for "quantum" temperatures of ± 0.6-0.5 MeV on the considered thermal excitation energy range. Figure 6 shows the final caloric curve with temperatures from quantum fluctuations (full squares). It exhibits a plateau around a temperature of 10-11 MeV on the E * range 5-10 AMeV. For comparison the caloric curve with classical temperatures derived from the simulation is also shown (open diamonds).

E* (A MeV)

Constrained caloric curves with quantum temperatures

Constrained caloric curves, which correspond to correlated values of E * and quantum corrected temperatures have been determined. Accordingly E * values, which were initially derived from experimental calorimetry including estimated correction for neutrons (see [10]), have been corrected a posteriori using quantum temperatures at freeze-out. Pressure values were also corrected using quantum temperatures. In figure 7 (left) we have constructed caloric curves for two different average freeze-out volumes corresponding to the ranges 3.0-4.0V 0 and 5.0-6.0V 0 where V 0 correspond to the volume of the QF nuclei at normal density. Again as theoretically expected a monotonous behaviour of caloric curves is observed. Figure 7 (right) shows the caloric curves when pressure has been constrained within two domains: 1.3-4.5 and 4.5-7.9 10 -2 MeVfm -3 . Backbending are clearly seen especially for the lower pressure range. For higher pressures the backbending of the caloric curve is largely reduced and one can estimate that the critical temperature is around 12-13 MeV for the selected hot nuclei.

E* (A MeV)

Discussion

First of all, as far as internal temperature for fragments are concerned (see figure 4), one can observe that they perfectly agree, for our A range, with those calculated with the well known "He/Li thermometer" used in [5,7] keeping the proposed prefactor 16. Note that similar temperature values are also obtained from our data (see figure 2 from [START_REF] Borderie | Proc. Int. Workshop on Multifragmentaton and related topics[END_REF]). This is a strong indication that the He/Li thermometer mainly reflects the internal temperature of fragments in relation with the strong secondary emission of alpha-particles and Li isotopes. As indicated in [7] those plateau temperatures can be interpreted as representing the limiting temperatures resulting from Coulomb instabilities of heated nuclei which has long been predicted [22]. Indeed, for thermally equilibrated QF hot nuclei one expects internal temperature for fragments equal to the limiting temperature of the fragmenting system. As a direct consequence internal fragment temperatures must reflect the evolution of limiting temperature with A of hot nuclei, which is indeed experimentally observed [7].

For a finite piece of nuclear matter like a hot nucleus the microcanonical ensemble is the most relevant ensemble and the kinetic/microcanonical temperature is the relevant parameter to build caloric curves and deduce information on a possible phase transition. On the other hand, one may also recall that dispersions in fragment transverse momentum spectra generated in projectile fragmentation were succesfully explained by adding the Fermi momenta of individual nucleons in fragments [23], which is a consequence of the quantum nature of nucleons. For multifragmentation reactions a similar approach was first proposed in [24] to qualitatively explain the different temperatures obtained from various observables (kinetic energies, level populations, isotope ratios). Thus, the use of the thermometer recently proposed based on quantum fluctuations [18] was a good opportunity to better investigate the behaviour of caloric curves especially under constraints in volume and pressure. The behaviour predicted by theoretical works for a first order phase transition is observed and confirms results obtained for other signals: negative microcanonical heat capacity and bimodality of an order parameter.

The last point that we want to discuss concerns information which can be directly deduced from experimental data like proton transverse momentum fluctuations. In that case as previously indicated Fermi energies can not be estimated and only classical temperatures can be calculated. From simulations we have seen previously that secondary decays distorted caloric curves. Results displayed in figure 8 fully confirm this observation: constrained caloric curves built with classical temperatures derived from experimental proton momenta also exhibit a monotonous behaviour and do not show any dependence upon average volumes or pressures, which prevent any direct information.

Conclusions

Several caloric curves have been derived for hot nuclei from quasi-fused systems using a new thermometer based on proton transverse momentum fluctuations including quantum effects. The unconstrained caloric curve exhibits a plateau at a temperature of around 10-11 MeV on the thermal excitation energy range 5-10 AMeV. For constrained caloric curves (volume and pressure) we observe what is expected for a first order phase transition for finite systems in the microcanonical ensemble, namely a monotonous behaviour at constant average volumes and backbending for constrained pressures. After the observation of negative microcanonical heat capacity and bimodality the behaviour of caloric curves is the ultimate signature of a first order phase transition for hot nuclei.

The exit, at high excitation, of the spinodal region and of the coexistence region around respectively 8 and 10 AMeV are in good agreement (within error bars) with spinodal [25] and bimodality [4] signals. Note also that 10 AMeV corresponds to first experimental indications for the onset of vaporization [26,27,28,29].

Finally, one can infer from the present studies that the temperatures obtained with the He/Li thermometer, which was largely used in the past, seem mainly reflect the internal temperatures of fragments in the excitation energy range 5-10 AMeV.

Figure 2 .

 2 Figure 2. Comparison betwwen the experimental (full points) and simulated (histograms) reduced relative velocity correlation functions for fragments. Each panel refers to a different beam energy. From[11].

Figure 3

 3 Figure 3. Centre of mass energy spectra for protons, deuterons, tritons,3 He,4 He and7 Li for the reaction at 39 AMeV. Open points are the simulated data, while full points are the experimental ones. From[11].

10 Figure 4 .

 104 Figure 4. Caloric curves (kinetic temperature versus thermal excitation energy) constrained at average volumes (left) and for selected ranges of pressure (right) and the corresponding internal temperatures for fragments. Error bars correspond to statistical errors. From[START_REF] Borderie | Proc. Int. Workshop on Multifragmentaton and related topics[END_REF].

Figure 5 .Figure 6 .

 56 Figure 5. Caloric curves (classical temperature from proton transverse momentum fluctuations versus thermal excitation energy) for protons (simulation) thermally emitted at freeze-out (open diamonds), for protons (simulation) after the secondary decay stage (full squares), and from protons experimentally measured (full points). Error bars correspond to statistical errors

10 Figure 7 . 10 Figure 8 .

 107108 Figure 7. Caloric curves (quantum corrected temperature versus thermal excitation energy) constrained at average volumes (left) and for selected ranges of pressure (right). Error bars include statistical and systematic errors.

Acknowledgments

The authors are indebted to A. Bonasera for private information and discussions.


Aurélie Courcoul
email: aurelie.courcoul@oniris-nantes.fr

Hervé Monod
email: herve.monod@jouy.inra.fr

Mirjam Nielen
email: m.nielen@uu.nl

Don Klinkenberg
email: d.klinkenberg@uu.nl

Lenny Hogerwerf
email: l.hogerwerf@uu.nl

François Beaudeau
email: francois.beaudeau@oniris-nantes.fr

Elisabeta Vergu
email: elisabeta.vergu@jouy.inra.fr

Modelling the effect of heterogeneity of shedding on the within herd Coxiella burnetii spread and identification of key parameters by sensitivity analysis

Keywords: Q fever, stochastic model, uncertainty and variability, PCA, ANOVA

Coxiella burnetii is the bacterium responsible for Q fever, a worldwide zoonosis. Ruminants, especially cattle, are recognized as the most important source of human infections. Although a great heterogeneity between shedder cows has been described, no previous studies have determined which features such as shedding route and duration or the quantity of bacteria shed have the strongest impact on the environmental contamination and thus on the zoonotic risk. Our objective was to identify key parameters whose variation highly influences C. burnetii spread within a dairy cattle herd, especially those related to the heterogeneity of shedding.

To compare the impact of epidemiological parameters on different dynamical aspects of C. burnetii infection, we performed a sensitivity analysis on an original stochastic model describing the bacterium spread and representing the individual variability of the shedding duration, routes and intensity as well as herd demography. This sensitivity analysis consisted of a Principal Component Analysis followed by an ANOVA. Our findings show that the most influential parameters are the probability distribution governing the levels of shedding, especially in vaginal mucus and faeces, the characteristics of the bacterium in the environment (i.e. its survival and the fraction of bacteria shed reaching the environment), and some physiological parameters related to the intermittency of shedding (transition probability from a non shedding infected state to a shedding state) or to the transition from one type of shedder to another one (transition probability from a seronegative shedding state to a seropositive shedding state).

Our study is crucial for the understanding of the dynamics of C. burnetii infection and optimization of control measures. Indeed, as control measures should impact the parameters influencing the bacterium spread most, our model can now be used to assess the effectiveness of different control strategies of Q fever within dairy cattle herds. the infection (**) * N denotes the herd size 8 ** the infection is assumed extinct when there is no I and C + left until the end of the simulation time 9

Introduction

Q fever is a worldwide zoonosis caused by Coxiella burnetii. This intracellular bacterium infects a wide range of animals and is associated with reproductive disorders in domestic ruminants [START_REF] Plommet | Renoux Fièvre Q expérimentale des bovins[END_REF][START_REF] Aitken | Clinical aspects and prevention of Q fever in animals[END_REF][START_REF] Bildfell | Coxiella burnetii infection is associated with placentitis in cases of bovine abortion[END_REF][START_REF] Berri | Relationships between the shedding of Coxiella burnetii, clinical signs and serological responses of 34 sheep[END_REF]. Goats, sheep and cattle are recognized as the main source of human infection [START_REF] Lyytikainen | An outbreak of sheep-associated Q fever in a rural community in Germany[END_REF][START_REF] Mcquiston | Q fever in humans and animals in the United States[END_REF][START_REF] Schimmer | The use of a geographic information system to identify a dairy goat farm as the most likely source of an urban Q-fever outbreak[END_REF][START_REF] Wallensten | Q fever outbreak in Cheltenham, United Kingdom, in 2007 and the use of dispersion modelling to investigate the possibility of airborne spread[END_REF]. Infected animals shed bacteria through various routes (parturition products, faeces, urine, vaginal mucus, milk) [START_REF] Berri | Relationships between the shedding of Coxiella burnetii, clinical signs and serological responses of 34 sheep[END_REF]9,[START_REF] Arricau-Bouvery | Excretion of Coxiella burnetii during an experimental infection of pregnant goats with an abortive goat strain CbC1[END_REF]. As the bacterium survives very well in the environment, humans can get infected by inhaling contaminated dusts or aerosols [START_REF] Tigertt | Airborne Q fever[END_REF][START_REF] Arricau-Bouvery | Is Q fever an emerging or re-emerging zoonosis?[END_REF][START_REF] Marrie | Q fever -a review[END_REF]. This was recently experienced in the Netherlands where more than 3,500 human cases were reported since 2007 [START_REF] Van Der Hoek | Q fever in the Netherlands: an update on the epidemiology and control measures[END_REF]. Although Q fever is asymptomatic in humans in more than 60% of cases, it can lead to acute or chronic infections and cause flu-like syndrome, hepatitis, pneumonia, endocarditis or abortions [START_REF] Raoult | Clinical and epidemiologic features of 1,383 infections[END_REF][START_REF] Gikas | Q fever: clinical manifestations and treatment[END_REF]. Hence, for public health and economic and animal health concerns, it is important to control C. burnetii infections in livestock herds.

In C. burnetii infections, a great heterogeneity between shedders has been described [START_REF] Durand | Lacteal and placental excretion of Coxiella burnetti, agent of Q fever, in the cow. Importance and prevention[END_REF][START_REF] Bouvery | Experimental Coxiella burnetii infection in pregnant goats: excretion routes[END_REF][START_REF] Rodolakis | Comparison of Coxiella burnetii shedding in milk of dairy bovine, caprine, and ovine herds[END_REF]: the shedding duration and routes, as well as the level of shedding (i.e. the quantities of bacteria shed) are variable between cows. According to Guatteo et al. [START_REF] Guatteo | Coxiella burnetii shedding by dairy cows[END_REF], cows can shed sporadically or persistently, the shedding routes are rarely concomitant and the concentrations of bacteria shed in vaginal mucus or milk can vary from less than 100 Bacteria/g to more than 1,000,000 B/g. Heterogeneity of shedding is known to affect infection dynamics in many diseases [START_REF] Matthews | New approaches to quantifying the spread of infection[END_REF] but it is generally difficult to determine which of its aspects are the most influential. The length of shedding, its route, the quantity of bacteria shed, or other features may all have a strong impact on the environment contamination by C. burnetii and thus on the zoonotic risk in the case of Q fever infection.

The identification of key parameters characterizing the heterogeneity of shedding and the inclusion of this heterogeneity in representations of disease spread within a herd are thus critical for the understanding of the infection dynamics. In addition, for some diseases, the effectiveness of control measures was shown to be dramatically improved by targeting the individuals transmitting the pathogen most (e.g. in the case of Escherichia coli O157 infection [START_REF] Matthews | Heterogeneous shedding of Escherichia coli O157 in cattle and its implications for control[END_REF], of measles epidemics [START_REF] Garske | The effect of superspreading on epidemic outbreak size distributions[END_REF] or of Salmonella transmission [START_REF] Lanzas | The effect of heterogeneous infectious period and contagiousness on the dynamics of Salmonella transmission in dairy cattle[END_REF]). However, understanding and predicting the spread of C. burnetii in a herd or identifying such key parameters cannot be assessed by field experiments alone. In this context, mathematical models and statistical methods are useful tools for understanding how the infection spreads within the herd and how various inputs (such as epidemiological characteristics of infected animals) affect the dynamics [START_REF] Keeling | Modeling ingectious diseases in Humans and animals[END_REF]. Techniques such as sensitivity analysis allow assessing the impact of the uncertainty and variability in the parameters on models outputs and hence determining key factors [START_REF] Saltelli | Sensitivity Anlysis[END_REF]. It consists in studying how the variation in the outputs of the model can be apportioned to different sources of variation, and how the model depends upon the information fed into it.

The aim of our study is to determine, using a sensitivity analysis, the key parameters related to the heterogeneity of shedding whose variation highly influences C. burnetii infection dynamics within a dairy cattle herd. First, the characteristics of the original model we developed for representing the spread of this infection will be summarized. Then, the sensitivity analysis performed, accounting for factors governed by probability distribution, will be detailed. This will be followed by the presentation and the discussion of the results.

Model

The original epidemic model that we developed is presented in details in the Supplementary material.

General description

In summary, the model describes the spread of C. burnetii within a dairy cattle herd and takes into account the heterogeneity of shedding: both the shedding routes and shedding levels (quantities of bacteria shed) are represented. Herd demography is also accounted for. This model is, to our knowledge, the first one proposed in the literature for Q fever spread coupling epidemiological aspects with herd population dynamics. Six different health statuses are defined function of the shedding status (i.e. excretion of bacteria or not), the shedding route if the animal is a shedder and the immunity status (presence or absence of antibodies). See Fig. 1 and Table 1 for parameter description. Each cow is then in one of the six mutually exclusive health states at a given time (Fig. 1): S (susceptible, non-shedder without antibodies), I -(shedder without antibodies), I + (shedder with antibodies), I + milk pers (shedder with antibodies, shedding in milk at higher levels and for a longer period of time than I + , as described by Guatteo et al. [START_REF] Guatteo | Coxiella burnetii shedding by dairy cows[END_REF]), C + (animal with antibodies, non-shedder but still infected), C -(non-shedder without antibodies which was infected and had antibodies in the past). Besides, all shedding cows I ( I -, I + and I + milk pers ) are subdivided according to their shedding routes: (1) I 1 , milk only, (2) I 2 , vaginal mucus and/or faeces, (3) I 3 , both. Susceptible animals can become contaminated through bacteria shed in the environment by infected animals. Once infected, animals progress through the different health states according to specific probability distributions (Table 1). The environment has its own dynamics (as described in the Supplementary material). The necessity to consider herd demography is generated by the interaction between lactation and gestation statuses and shedding routes and levels: cows which recently calved (during the 4 weeks post-calving) have different probability distributions governing shedding routes and levels than other cows. See Table 2 for description of herd demography related parameters.

The model is stochastic, individual-based and in discrete time with a time step of one week, which is appropriate for both epidemiological and herd management processes. The stochasticity has two main sources: for each individual, all the transitions between health states are supposed stochastic and the quantities of bacteria shed in the environment follow discrete distributions, different according to the shedding route.

Initial conditions and parameter values of the standard scenario

At t=0, the herd consists of 50 cows. To initiate the infection cycle, a primiparous I + cow which has just calved is introduced into a wholly susceptible herd.

The epidemiological parameters are put at their standard values (Table 1): parameters m (transition probability I -=> S), q (transition probability I -=> I + ), r 1 (transition probability I + => C + ), s (transition probability C + => I + ) and μ (elimination rate of C. burnetii) come from a study where they were estimated through Bayesian inference using data from five French chronically infected dairy cattle herds [START_REF] Courcoul | Spread of Q fever within dairy cattle herds: key parameters inferred using a Bayesian approach[END_REF]; pIp (proportion of cows going from I -to I + and becoming I +milk pers ) and probav (probability of abortion after a transition S => I -, C + => I + and C -=> I + ) were qualitatively calibrated to match respectively the average number of I +milk pers cows and the distribution of abortions observed in the field (R. Guatteo and A.F. Taurel, 2010, personal communication); for r 2 (transition probability I +milk pers => C + ), it was assumed that the shedding duration of I +milk pers cows (i.e. 1/ r 2 ) was on average 10 times longer than the shedding duration of I + cows (i.e. 1/ r 1 ); ρ mf (proportion of bacteria shed through mucus/faeces arriving into the environment compartment) and ratio

milk mf ρ ρ ⎛ ⎞ ⎜ ⎟ ⎝ ⎠
) were calibrated from expert opinion to match the environmental bacterial load inferred in Courcoul et al. [START_REF] Courcoul | Spread of Q fever within dairy cattle herds: key parameters inferred using a Bayesian approach[END_REF]; for τ (transition probability C + => C -), it was assumed, based on Fournier et al. [28] and Plommet et al.

[29], that the mean life duration of antibodies in cattle was 2 years; probability distributions of shedding related parameters, α, β, β calv , γ, γ calv , (all five representing probability distributions of the shedding routes for the different types of shedders) and Q1, Q2, Q3, Q4 and Q5 (all five representing probability distributions of the shedding levels for the different types of shedders) were qualitatively calibrated from field data (R. Guatteo 2009, personal communication).;The parameters governing the demography and herd management (Table 2) were chosen to represent a standard French dairy cattle herd.

To account for the variability in Q fever infections, 200 repetitions of the standard scenario were run over a 5-year simulation period.

Sensitivity analysis

We conducted a sensitivity analysis to identify the parameters that mostly contributed to the output variability. This analysis, similar to that performed in Lurette et al. [START_REF] Lurette | Sensitivity analysis to identify key parameters influencing Salmonella infection dynamics in a pig batch[END_REF] following the methodology developed in Lamboni et al. [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] involves an original aspect (not present in [START_REF] Lurette | Sensitivity analysis to identify key parameters influencing Salmonella infection dynamics in a pig batch[END_REF]): some factors are governed by probability distributions (such as heterogeneity related parameters) and are not point values, as usual in this kind of methods. Various scenarios were run, each of them being characterized by a specific combination of parameter values, in order to relate the variability obtained for the outputs to that induced by the input parameters.

Outputs and factors

Eight outputs were considered (Table 3): (i) the environmental bacterial load of the main buildings and close pastures for lactating cows E building , (ii) the environmental bacterial load of the specific pastures for dry cows E dry , (iii) the prevalence of milk shedders, (iv) the prevalence of mucus/faeces shedders, (v) the prevalence of shedders in milk in a persistent way, (vi) the seroprevalence, (vii) the number of abortions per herd per year, and (viii) the extinction rate. All these outputs except the number of abortions and the extinction rate were computed weekly over a 5-year period.

Parameters related to the herd demography were fixed at their nominal values of Table 2 since demography and herd management processes are considered as well known. The sensitivity of the model outputs was evaluated with respect to the epidemiological parameters, which are those given in Table 1, except for τ. These 19 parameters presented in Table 1 are thus the inputs of the sensitivity analysis and they will be called factors in the rest of the paper. They belong to two categories: parameters concerning the transitions between health states and the environment (m, q, pIp, r 1 , r 2 , s, μ)

and parameters directly related to the heterogeneity in shedding (i.e. α, β, β calv , γ, γ calv , Q1, Q2, Q3, Q4, Q5, ρ mf and ratio

milk mf ρ ρ ⎛ ⎞ ⎜ ⎟ ⎝ ⎠
). The parameters in the first category were estimated from field data previously [START_REF] Courcoul | Spread of Q fever within dairy cattle herds: key parameters inferred using a Bayesian approach[END_REF], but some uncertainty still remains (due for instance to the limitation of the data). They were included in the sensitivity analysis but with relatively limited ranges of variation. In this study we focused on the latter category of parameters because they directly describe heterogeneity related aspects, which represented our main objective.

Design of experiments

All the designs were generated using R 2.10.1 [START_REF]R: A Language and Environment for Statistical Computing[END_REF] and PLANOR R package [START_REF] Bouvier | PLANOR : an R library for the automatic generation of regular fractional factorial designs[END_REF].

First experiment

We used a fractional factorial experiment design, with four values (called levels) per factor related to the shedding and two levels for the other factors (values in Table 1). As our model is stochastic, we ran the model for each combination of factor levels 30 times. Since the complete factorial design would lead to too many combinations (exactly 30 x 4 12 x 2 7 simulations), a fractional factorial design of resolution V was chosen. Such a design allows to estimate the main effects and two-factor interactions, provided higher order interactions are assumed to be negligible [START_REF] Bailey | Design of comparative experiments[END_REF][START_REF] Kobilinsky | Les plans factoriels[END_REF]. In the present case, a design was obtained with 4,096 scenarios. Thus, we ran 122,880 realizations of the model (i.e. 30 repetitions for each of the 4,096 scenarios).

Second experiment

A complete factorial design for the eight most influential factors according to the first experiment was performed. This enabled us to more accurately quantify the impact of the interactions between these eight factors and also to disentangle potential confounded main effects and interactions. Besides, we determined in this experiment the factors that mostly contributed to the variability of the extinction rate between repetitions. The remaining 11 parameters were put to their standard value (Table 1). For this second study, we ran 2,048 scenarios with 30 repetitions each.

Third experiment

In a third analysis, the influence of the probability distributions of the shedding levels (factors Q) was specifically explored. Since these probability distributions depend on the type of shedder ( I -, I + or I + milk pers ) and route of shedding (milk versus mucus/faeces), this analysis enabled us to explore which type of A fractional factorial design for the 10 new factors Q* with four levels each was generated. These four levels were (0.85, 0.15, 0), (0.6, 0.4, 0), (0.25, 0.25, 0.5) and (0.15, 0.6, 0.25) for the probability to be in (low, mid, high) shedding level respectively. The other parameters were put to their standard values given in Table 1. For this third experiment, we ran 1,024 scenarios with 30 repetitions each.

Analysis of the temporal outputs (of the first, second and third experiments)

In order to compare the influence of factors on the seven outputs which exhibit temporal dynamics (all outputs except the extinction rate), we applied a method developed by Lamboni et al. [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF]. The results are recorded as tables with one row for each scenario and one column for each output time point (260 weekly time points for outputs (i) to (vi) and five annual time points for the abortion number, output

(vii)).
This method allows to simultaneously analyze potentially correlated variables (here the successive time points of a given output). It consists in two main steps. First, a Principal Component Analysis (PCA) is operated in order to provide linear combinations (or components) of the initial variables (here the columns of our tables) explaining the maximum of inertia (i.e. variability) between scenarios. Only the first three principal components (PC) were kept since they are sufficient to cover most variability amongst simulations. The PCA provides to each line of the tables a score on each component. The second step involves an ANOVA, including the main effects and the two-factor interactions for all factors and carried out on the scores of each of the components considered. Sensitivity indices (SI), corresponding to the main effect or to interactions, and total sensitivities (TS), corresponding to the sum of the main effect and the interactions, were calculated for each factor and for each component. This analysis was performed with R 2.10.1 [START_REF]R: A Language and Environment for Statistical Computing[END_REF] and multisensi R package [START_REF] Lamboni | Multisensi: Multivariate Sensitivity Analysis. R package version 1[END_REF].

The analyses were performed on both the mean and standard deviation of the 30 repetitions of each scenario, in order to assess the two sources of variability influencing the outputs: the model intrinsic stochasticity and the parameter variability generated by the factorial designs.

Analysis of the extinction rate (of the second experiment)

An ANOVA was performed to assess the influence of the eight most influential factors on the extinction rate. It was calculated for each scenario defined by the complete factorial design of the second experiment.

Analysis of the outputs at the time point 260 (of the first, second and third experiments)

In order to determine the factors with the highest influence on the output variability as a whole, we performed a joint analysis on the values of the six dynamic outputs (first six lines of Table 3) at the last simulation time step (week 260). This time point was chosen to illustrate the long-term steady-state of the system. Thus, the two-step analysis (PCA followed by ANOVA) was performed twice on six output variables. The first analysis was done on the mean and the second on the standard deviation of the 6 dynamic outputs at time 260.

Results

Infection dynamics of the standard scenario

Over the 200 repetitions of the standard scenario, 34 led to the extinction of infection (defined as the after the initial infection. The ratio between the mean prevalence of shedders and the mean prevalence of milk shedders was around 2.6 in the first weeks of simulation and then it decreased to reach 1.9 at the end of the simulation time. The mean environmental bacterial load E building increased with time corresponding to a mean transition probability from S to I -, p building , equal to 0.45 at the end of the simulation time. On the contrary, the mean environmental bacterial load E dry was close to 0 for the 5 years of simulation (results not shown). The median number of abortions varies from 2 per herd in the first year to 4 per herd in the 5 th year, and a large variability surrounded these values [0-10 for the percentiles 2.5% and 97.5% respectively]. In addition, as shown in Figure 3, the route and the level of shedding of a shedder cow had a great impact on the contamination of the environment. This result is an unsurprising consequence of the model parameterization. As expected, the 'low level' shedding category, although the most common, did not contribute much to the increase of the environmental bacterial load. On the contrary, shedders in mucus/faeces of the 'mid level' category and shedders in milk of the 'high level' category filled the environment in a non negligible way. Above all, shedders in mucus/faeces of the 'high level' category (both non aborting and aborting cows) had the greatest impact.

absence

Influence of the epidemiological factors on the model outputs

The results obtained with 30 runs for each parameter set were robust: the mean and the percentiles 2.5, 50 and 97.5 of our outputs were similar to those obtained with 200 runs (results not shown). The outputs variability due to the model stochasticity (within scenario variability) is illustrated Figure S5 of the Supplementary material.

The variability between scenarios was important. As an example, for the first experiment, the percentiles 2.5%, 50% and 97.5% of the seroprevalence (mean over the 30 repetitions of a given scenario) at the end of the 5-year simulation period were of 0.0%, 45.4% and 94.5% respectively. Those of the environmental bacterial load (mean over the 30 repetitions of a given scenario) were of 0.0, 2.1 and 51.0 respectively (Figure S4 of the Supplementary material).

First experiment

As shown in Table 4, since the inertia obtained for the first principal component (PC) was very high for each of the model outputs (except the joined analysis at time 260), only the results on the first PC are presented.

For the two mean environmental bacterial loads, the factors Q1 (probability distribution of the shedding levels for all the I -and for the I + shedding in mucus/faeces after 4 weeks post-calving), μ (elimination rate of C. burnetii) and ρ mf (proportion of bacteria shed through mucus/faeces arriving into the environment) were the most influential ones.

For the mean prevalences of mucus/faeces and milk shedders, the most sensitive factors were q (transition probability from I -to I + ), s (transition probability from C + to I + representing the intermittency of shedding) and Q1 (probability distribution of the shedding levels for all the I -and for some of the I + in mucus/faeces), whereas the mean prevalence of milk shedders in a persistent way was mostly impacted by pIp (the proportion of cows going from I -to I + and becoming I +milk pers ), q (transition probability from I -to I + ) and r 2 (transition probability from I + milk pers to C + ). Concerning the mean seroprevalence, the factor q (transition probability from I -to I + ) had a TS higher than 60%.

Lastly, the most influential factors of the mean abortion number were q (transition probability from I - to I + ), Q1 (probability distribution of the shedding levels for all the I -and for some of the I + in mucus/faeces), s (transition probability from C + to I + ), μ (elimination rate of C. burnetii) and ρ mf (proportion of bacteria shed through mucus/faeces arriving into the environment).

Globally, the most sensitive two-factor interactions (with a SI higher than 5%) were Q1:q on the variability of the number of abortions, Q1:μ, ρ mf :Q1, ρ mf :μ on the variability of the environmental bacterial loads and q:pIp, q:r 2 , pIp:r 2 on the variability of the prevalence of milk shedders in a persistent way.

Concerning the variability of the standard deviations of the outputs, the same factors as above were identified as the most influential ones for the environmental bacterial loads, the prevalence of milk shedders in a persistent way and the abortion number. The main effect of the factors was always very low (no SI higher than 5%) on the prevalences of mucus/faeces shedders, whereas the part of two-factor interactions was much more important. The most sensitive factors were Q1 (probability distribution of the shedding levels for all the I -and for some of the I + in mucus/faeces), ρ mf (proportion of bacteria shed through mucus/faeces arriving into the environment), q (transition probability from I -to I + ) and Q 3 (probability distribution of the shedding levels for all the I + during the 4 first weeks post-calving).

and for some of the I + in mucus/faeces) and α (probability distribution of the shedding routes for the I - cows).

For the joined analysis on six of the dynamic outputs at time 260, the inertia obtained for the first PC was much lower and the second PC had to be taken into account. For the means, the most influential factors were Q1 (probability distribution of the shedding levels for all the I -and for some of the I + in mucus/faeces), q (transition probability from I -to I + ), s (transition probability from C + to I + ), μ

(elimination rate of C. burnetii) and ρ mf (proportion of bacteria shed through mucus/faeces arriving into the environment) on the first PC and Q1, μ, ρ mf and q on the second one, by order of importance. For the standard deviation, the most sensitive factors were q, Q1, Q 3 (probability distribution of the shedding levels for all the I + during the 4 first weeks post-calving), r 2 (transition probability from I + milk pers to C + ) and pIp (the proportion of cows going from I -to I + and becoming I +milk pers ) on the first PC and Q1, ρ mf and μ on the second.

Second experiment

The eight factors chosen in the second experiment (appearing as the most influential according to the findings of the first experiment) were two of the five shedding level probability distributions, Q1 and Q3, three transition probabilities q, s, and r 2 , the elimination rate of the bacterium in the environment μ, the proportion of bacteria shed through mucus/faeces arriving into the environment, ρ mf and the proportion of cows going from I -to I + and becoming I +milk pers , pIp. The results obtained were globally similar to the results of the first experiment described above (same most influential factors, sometimes in a slightly different order), suggesting that no important interactions were confounded with the main effects in the first analysis.

The most sensitive factors on the extinction rate (frequency distribution amongst scenarios in Figure S3 of Supplementary material) were firstly Q1, μ and ρ mf (with a SI higher than 10%) then q, s, and Q3

(with a lesser SI, but higher than 5%). The most sensitive two-factor interactions were Q1:q and Q1:μ (SI higher than 5%).

Third experiment

Amongst the ten probability distributions of the shedding levels Q*, the most influential factor on the variability of the outputs was Q2*, the probability distribution of the shedding levels for the I - shedding in mucus/faeces (Table 5 and Figure 5). Q4* and Q6*, the probability distributions of the shedding levels for the I + shedding in mucus/faeces before and after the 4 th week post-calving, respectively, also had a significant impact and, to a lesser extent, Q1* , the probability distribution of the shedding levels for the I -shedding in milk, and Q8*, the one for the I + milk pers shedding in mucus/faeces.

Moreover, the only interaction among the five most sensitive terms was Q2*:Q4*. Overall, the factors with the greatest impact were probability distributions of the shedding levels in mucus/faeces.

Discussion

The main objective of this work was the identification, using a sensitivity analysis, of the key parameters related to the heterogeneity of shedding highly impacting C. burnetii infection dynamics within a dairy cattle herd. The underlying step of the achievement of this goal was the elaboration of a model describing the bacterium spread.

We then proposed, in the first part of our study, the first model of C. burnetii spread within a dairy cattle herd taking into account the individual variability of shedding, defined through duration, routes and intensity. Simulated infection dynamics are consistent with field data: at the last time point of the simulated time series (five years after the introduction of the initial infectious case), the median seroprevalence is around 31% [11.0% -41.2% for the 25 th and 75 th percentiles respectively], which is globally consistent with the mean observed seroprevalence in cows (mean: 40%, 25 th and 75 th percentiles: 25% and 51% respectively) of 56 naturally infected French dairy herds [START_REF] Taurel | Q fever within-herd seroprevalence in infected dairy herds: assessment using an ELISA applied to bulk tank milk[END_REF]. At the same simulated time point, the median prevalence of shedders is around 32.3% [0 -56.9% for the percentiles 2.5 th and 97.5 th respectively] whereas in the field, the apparent proportion of shedder cows is 45.5% in Guatteo et al. [9] and 38.9% in Rodolakis et al. [START_REF] Rodolakis | Comparison of Coxiella burnetii shedding in milk of dairy bovine, caprine, and ovine herds[END_REF].

The second and main part of our study consisted in the sensitivity analysis. This approach aims at improving the understanding of complex systems such as stochastic epidemiologic models and at suggesting possible targeted control strategies in livestock infections [START_REF] Ezanno | Sensitivity analysis to identify key-parameters in modelling the spread of bovine viral diarrhoea virus in a dairy herd[END_REF][START_REF] Turner | A semi-stochastic model of the transmission of Escherichia coli O157 in a typical UK dairy herd: dynamics, sensitivity analysis and intervention/prevention strategies[END_REF][START_REF] Lu | Stochastic simulations of a multi-group compartmental model for Johne's disease on US dairy herds with test-based culling intervention[END_REF] or assessing the effect of varying the input parameters on the economic impact associated with the disease [START_REF] Weber | Milk quality assurance for paratuberculosis: simulation of within-herd infection dynamics and economics[END_REF][START_REF] Bruijnis | Assessing economic consequences of foot disorders in dairy cattle using a dynamic stochastic simulation model[END_REF].

To perform sensitivity analyses, we used complete and fractional factorial designs. Alternative approaches are available (see for example [START_REF] Saltelli | Sensitivity Analysis[END_REF] or [START_REF] Monod | Uncertainty and sensitivity analysis for crop models[END_REF]), but factorial designs are very convenient to control which main effects and interactions can be estimated and to manage a mixture of qualitative and quantitative factors. In the present study, some factors (the probability distributions of the shedding routes and levels for the different types of shedders α, β, β calv , γ, γ calv , Q1, Q2, Q3, Q4 and Q5) were not scalars but probability distributions. Performing a sensitivity analysis with such types of factors is to our knowledge fairly new. We used multinomial distributions with three classes and defined the modalities of each such factor as four alternative sets of multinomial probabilities. This choice allowed flexibility as well as a fine control in the probability distributions that were simulated.

To cope with the dynamic and multivariate outputs of the model, the PCA+ANOVA approach [START_REF] Lamboni | Multivariate global sensitivity analysis for dynamic crop models[END_REF] offered additional insight compared to single sensitivity analyses. Multivariate decomposition methods other than PCA could be used, but the key idea is that sensitivity analyses are now performed on synthetic and meaningful output variables. As the infection dynamics is composed of two phases (see on

Figure 2 an initial phase of rapid progression followed by a kind of steady state), the parameters influencing the dynamics may not be the same between the phases. We additionally performed a preliminary sensitivity analysis on the very initial phase of infection (first 26 weeks; results not shown).

The most sensitive factors during the first 15 weeks were highly influenced by the initial conditions but very fast afterwards, the same factors as in experiments one or two were identified as the most influential ones. We then chose to conduct our study on the whole simulation period.

Another aspect concerns the stochasticity of our model which generates complex dynamics. Hence, attention has to be paid to this point in order to appropriately apply sensitivity analysis approaches that are mostly developed for deterministic models. Here, to be able to study with confidence how the variation in the outputs of our model can be apportioned to the uncertainty of epidemiological factors, we checked that the variability due to the model stochasticity (i.e. the within scenario variability) was negligible compared to the variability due to the input parameters variability (i.e. the between scenarios variability). More specifically, to provide reliable analysis on trends, means were calculated on 30 repetitions. Standard deviations were also considered as they can provide complementary information on the most influential factors.

In summary, the method used has the major advantage of allowing to consider temporal varying outputs and thus to identify parameters influencing the dynamics over the 5 year simulation period. Moreover, it allows properly dealing with the model stochasticity.

According to our findings (first and second experiments), we can classify the eight parameters that have most influence on the C. burnetii infection dynamics in three categories, depending on the dynamics aspects they are involved in. These parameters should be precisely assessed in further studies to improve the model accuracy and the understanding of the infection spread. Control strategies impacting them would also strongly influence the infection dynamics.

The first group of key parameters, comprising the parameters related to the transition between health states, slightly influences the different prevalences and the abortion number. The transition probability from I -to I + (q) is a physiological parameter associated to seroconversion. It was estimated based on data from a follow-up of five chronically infected herds [START_REF] Courcoul | Spread of Q fever within dairy cattle herds: key parameters inferred using a Bayesian approach[END_REF]. However, we can assume that those five herds do not cover the whole potential range of variability of this factor, especially at the beginning of the infection when this parameter probably takes a different value depending on how recently the infection occurred in the herd. Further experimental or survey studies focusing on the start of the infection dynamics are needed to improve the knowledge of this parameter, and therefore the accuracy of the model. However, it is unlikely that a control measure could impact this seroconversion parameter.

The transition probability from C + to I + (s) represents the intermittency of shedding. It was inferred from data in the same study [START_REF] Courcoul | Spread of Q fever within dairy cattle herds: key parameters inferred using a Bayesian approach[END_REF] and the estimated values were herd-dependent. It is biologically plausible to assume that a control measure such as vaccination could decrease this parameter and then have an impact on the prevalence of shedders. However, since, to our knowledge, no data on the impact of vaccination on the intermittency of shedding is currently available, further studies are needed. The transition probability from I + milk pers to C + (r2), and the proportion of cows going from I -to I + and becoming I + milk pers (pIp) have an impact on the variability in the prevalence of persistent milk shedders (I + milk pers ). These parameters were not estimated from data, but calibrated so that simulated trajectories of prevalence of persistent milk shedders are consistent with field observations. Indeed, following Guatteo assumed to be five weeks in the standard scenario and two or 13 weeks (two extreme values) in the sensitivity analysis. However, this assumption is difficult to verify. Based on Dutch studies, it seems that within a month, more than 75% of the manure does not contain viable C. burnetii anymore but that the bacterium survives only a few days in the outer layer of the manure [START_REF] Vellema | Q fever in the Netherlands: One Health in relation to Q fever, in humans and animals[END_REF]. It is then difficult to calibrate the elimination rate of C. burnetii in the environment which represents both the natural mortality of the bacterium and its removal due to the periodic cleaning of the cattle housing carried out by the farmer.

However, it seems possible to influence it by implementing environmental control measures such as increased cleaning of the farm.

Lastly, as suggested by the results of the first experiment (Table 4) and detailed in the third experiment (Table 5), the parameters which have the greatest impact on the infection dynamics are those governing the shedding levels (through their associated probability distributions Q*), especially in mucus/faeces.

The specific effects of mean, variance or form of these distributions have not been separately evaluated when adjusting for constant overall transmissibility; the influence of these distributions can therefore be considered as an effect of the mean quantities of bacteria shed. As shedding in mucus/faeces contaminates the environment much more than shedding in milk, it could seem surprising that the parameters governing the probability distribution of the shedding routes (α, β and γ) do not appear to influence the model outputs. This could partly be explained by the parameter values used, especially for the quantity of bacteria released by shedders in low, mid and high levels respectively (Qty). The quantity of bacteria released by a high level shedding cow is so high (compared to mid or low levels), that the probability distributions governing the levels (such as Q in experiments one and two) are by construction more important than those related to the shedding routes. However, although exact values of the quantities of bacteria shed by shedders in low, mid and high levels are unknown, the standard values used in this study were calibrated according to field data, which tend to support our findings.

This third experiment was of high interest as it allows highlighting the importance of a sub-category of animals. Indeed, the distribution of the levels of bacteria shed in mucus/faeces by seronegative shedders (factor Q2*) has a strong impact on all the outputs, including the shedder prevalences and the environmental bacterial loads. The distributions of the levels of bacteria shed in mucus/faeces by seropositive shedders (factors Q4* and Q6*) also have a significant impact on the variability of the shedder prevalences. The predominance of Q2* over Q4* and Q6* can be explained by the larger simulated number of seronegative shedders than seropositive ones. In fact, due to the low standard value of the transition probability from I -to I + (q) and the high standard value of the transition probability from I -to S (m), the number of seropositive cows is very low during the three first years of simulation compared to the number of seronegative cows. The standard values of these two transition probabilities, which were estimated in chronically infected herds, are perhaps not perfectly appropriate to describe the initial phase of an infection and then could lead to overestimation of the influence of seronegative shedders.

As suggested by Matthews et al. for Escherichia coli O157 [START_REF] Matthews | Heterogeneous shedding of Escherichia coli O157 in cattle and its implications for control[END_REF], identifying factors such as age, genetics, reproductive status or other management factors that might predispose an animal to high levels of shedding would be of undisputable interest. Moreover, control measures should aim at reducing the proportion of high shedders in mucus/faeces, such as phase I vaccines seem to do. According to Arricau-Bouvery et al. [START_REF] Arricau-Bouvery | Effect of vaccination with phase I and phase II Coxiella burnetii vaccines in pregnant goats[END_REF], vaccination dramatically reduced both abortions and excretion of bacteria in the milk, vaginal mucus and faeces. In Rousset et al. [START_REF] Rousset | Efficiency of a phase 1 vaccine for the reduction of vaginal Coxiella burnetii shedding in a clinically affected goat herd[END_REF], the vaccine was effective in reducing massive bacterial shedding from a heavily infected goat herd. In the same way, Hogerwerf et al. [START_REF] Hogerwerf | Reduction of Coxiella burnetii prevalence by vaccination of goats and sheep, the Netherlands[END_REF] found that in uterine fluid, vaginal mucus and milk of vaccinated dairy goats, both the prevalence of shedders as well as the concentration of bacteria were reduced. are the sum of all quantities Qty of bacteria shed by all the shedders through all the shedding routes times ρ, the fraction of bacteria shed reaching the environment of the herd. Table 1. Parameters of the epidemiological model: description, standard values and values tested in the first experiment of the sensitivity analysis.

Conclusion

Table 2. Description of the parameters of the herd demography model and their standard values.

Table 3. Description of the outputs of the sensitivity analysis.

Table 4. Results of the first experiment of the sensitivity analysis: for the mean and standard deviation of each output, the proportion of total inertia represented by the first principal component and the corresponding most influential factors are given. Factors with total sensitivities higher than 10%: one cross, factors with total sensitivities ≥ 25%: two crosses, factors with total sensitivities ≥ 50%: three crosses.

Table 5. Results of the third experiment of the sensitivity analysis: for the mean and standard deviation of each output, the proportion of total inertia represented by the first principal component and the corresponding most influential factors are given. Factors with total sensitivities higher than 10%: one cross, factors with total sensitivities ≥ 25%: two crosses, factors with total sensitivities ≥ 50%: three crosses.

Table 1

  shedders and which type of shedding route played a major role in the variability of the outputs. Ten new factors (called Q*) were defined: each factor Q* corresponds to a given shedding level probability distribution for a type of shedder and a shedding route. These 10 factors Q* replaced the former 5 factors Q used in first and second experiments (where a same shedding level probability distribution Q could be used for different types of shedders or different shedding routes). Factors Q* were numbered as follows: for milk and mucus/faeces respectively, Q1* and Q2* refer to the probability distributions of the shedding levels for the I -, Q3* and Q4* refer to those for the I + after 4 weeks post-calving, Q5* and Q6* refer to those for the I + during the 4 first weeks post-calving, Q7* and Q8* refer to those for the milk pers I + after 4 weeks post-calving, and finally Q9* and Q10* refer to those for the milk pers I + during the 4 first weeks post-calving. Thus, former factor Q1 of first and second experiments corresponds to new factors Q1*, Q2* and Q4*, former factor Q2 to Q3*, former factor Q3 to Q5* and Q6*, former factor Q4 to Q8* and former factor Q5 to Q7*, Q9* and Q10*.

  This work led to the identification of key parameters in C. burnetii infection dynamics based on a methodology combining novel aspects of sensitivity analysis and an original model describing the bacterium spread within a dairy cattle herd composed of animals with heterogeneous shedding characteristics. The most influential parameters are associated to the probabilities governing the levels of shedding, especially for mucus/faeces shedders and to the characteristics of the bacterium in the environment. Some physiological parameters related to the intermittency of shedding or to the transition from one type of shedder to another one also play a non negligible role. Our study also highlights parameters that should be precisely assessed and then further investigated to improve the model accuracy and the understanding of the infection spread. Besides, interventions impacting those key parameters would be of great interest. The model developed here can be further used to assess over a longer time scale the effectiveness of different control strategies for C. burnetii infection, such as vaccination.

Figure 2 .

 2 Figure 2. Temporal dynamics of the seroprevalence, prevalence of shedders, prevalence of milk

Figure 3 .

 3 Figure 3. Number of shedders for each shedding route and each shedding level and their contributions in

Figure 4 .

 4 Figure 4. Sensitivity analysis on the mean prevalence in mucus/faeces shedders over time: results of the

Figure 5 .

 5 Figure 5. Sensitivity analysis of the means of six of the outputs (all except the abortion number) for the

  the shedding levels for all the I -and for the I + shedding in 0the shedding levels for all the I + in the 4 first weeks postthe shedding levels for the I +milk pers shedding in mucus/faeces after 4 weeks postthe shedding levels for all the I +milk pers shedding in milk and for the I +milk pers in the 4 first weeks postet al.[START_REF] Courcoul | Spread of Q fever within dairy cattle herds: key parameters inferred using a Bayesian approach[END_REF] b calibrated to match field data (R. Guatteo 2009, personal communication) * includes both the natural mortality of the bacterium and its removal in relation to the periodic cleaning of the cattle housing 2

  of the herd shedding in vaginal mucus and/or faeces, i.e. in milk in a persistent way Proportion of animals I + milk pers , i.e. with antibodies, i.e.

  of runs of a particular scenario leading to an extinction of

  of shedders (animals in I states) or chronically infected animals (in C + state) in the herd)

	occurring in week 60 (median; see Figure S2 in the Supplementary material) after the introduction of the
	initial infected cow (min: week 2, max: week 255). The mean seroprevalence and the mean prevalences
	of shedders increased with time (Fig 2) to reach respectively 34.4% [0 -57% for the percentiles 2.5%
	and 97.5% respectively] and 36.3% [0 -60% for the percentiles 2.5% and 97.5% respectively] five years

Table 2

 2

	4
	Description		Standard value
	Replacement rate (year -1 )		0.355
		lactation 1	0.0057
		lactation 2	0.0052
	Culling rate (week -1 )	lactation 3	0.0065
		lactation 4	0.0067
		lactations 5&6 0.0161
		lactation 1	0.337
		lactation 2	0.252
	Probability distribution at	lactation 3	0.173
	time 0 for the lactation
	numbers of the cows	lactation 4	0.11
		lactation 5	0.088
		lactation 6	0.04

Table 3 7

 3 Output name DescriptionE buildingEnvironmental bacterial load of the main buildings and close pasturesE dryEnvironmental bacterial load of the specific pastures for the dry cowsPrevalence of milk sheddersProportion of animals of the herd shedding in milk, i.e.

Table 4

 4

	10
										Prevalence of
		Environment E building	Environment E dry	Prevalence of milk shedders	Prevalence of m/f shedders	milk shedders in a persistent	Seroprevalence	Abortion number	Joined analysis time point 260
										way
		mean	sd	mean	sd	mean	sd	mean	sd	mean	sd	mean	sd	mean	sd	mean	sd
	inertia PC1	99.1% 96.3% 98.4% 96.8% 94.4% 87.5% 94.1% 80.7% 93.1% 87.7% 97.0% 86.8% 94.0% 88.0% 49.0% 43.8%
	m
	q					++	+	+	+	++	++	+++	++	++	++	+	+
	pIp									++	++						+
	r 1
	r 2									++	++						+
	s					+		+					+	+	+	+
	μ	++	++	++	++									+	+	+
	ρ mf	++	++	++	++				+					+	+	+
	ratio
	ρ milk /ρ mf
	α						+

Acknowledgements

The authors would like to thank Alain Joly and Raphaël Guatteo from the Oniris-INRA group 'Bioaggression, Epidemiology and Risk Analysis' for fruitful discussions on the epidemiological model and on the way to take into account the heterogeneity of shedding. The collaboration between the Faculty of Veterinary Medicine of Utrecht and the French National Institute for Agricultural Research (INRA) was financially supported by the Netherlands Organisation for Scientific Research (NWO) and the French Ministry of Foreign and European Affairs through the Van Gogh Programme.

et al. [START_REF] Guatteo | Coxiella burnetii shedding by dairy cows[END_REF], almost one milk shedder out of three was detected as persistent shedder over three months.

In our case, the mean prevalence of milk shedders at time point 260 weeks is 3.3 times the mean prevalence of persistent milk shedders. While the proportion of cows becoming I + milk pers seems difficult to be decreased, the transition probability from I + milk pers to C + (r 2 ), which rules the shedding duration of persistent milk shedders, could probably be modulated by control strategies. Although according to Astobiza et al. [START_REF] Astobiza | Kinetics of Coxiella burnetii excretion in a commercial dairy sheep flock after treatment with oxytetracycline[END_REF], an oxytetracycline treatment would not limit the duration of bacterial excretion in a dairy sheep flock, further studies are needed in order to determine if vaccination could decrease the duration of the shedding period.

The second group of parameters that influenced the infection dynamics most is related to the characteristics of C. burnetii in the environment. In fact, the proportion of bacteria shed through mucus/faeces and arriving into the environment (ρ mf ) and the elimination rate of C. burnetii in the environment of the farm (μ) have a strong impact on the environmental bacterial load and to a lesser extent on the abortion number. It is very difficult to quantify in practice which proportion of bacteria contained in milk, vaginal mucus or faeces contaminates the environment. Thus, we calibrated the proportion of bacteria shed reaching the environment to match the environmental bacterial load estimated by Courcoul et al. [START_REF] Courcoul | Spread of Q fever within dairy cattle herds: key parameters inferred using a Bayesian approach[END_REF]: the posterior medians for the environmental bacterial loads of each of the five chronically infected herds were comprised between 0.044 [0.005-0.143 for the 95% credible interval (CI)] and 0.558 units of environment [0.201-1.278 for the 95% CI]. Since those herds did not exhibit any clinical signs, we assumed that the simulated environmental bacterial load in herds with abortions should be slightly higher. The median environmental bacterial load at time step 260 weeks is then 0.57 (with the percentiles 2.5 th and 97.5 th equal to 0.00 and 1.49 respectively). Concerning the elimination rate of the bacterium in the environment, given that C. burnetii withstands harsh environmental conditions [START_REF] Marrie | Q fever -a review[END_REF], its life expectancy (1/μ) within the farm in an infectious form was