Mircea T Sofonea
email: mircea.sofonea@normalesup.org

Lafi Aldakak

Luis Fernando Boullosa

Samuel Alizon

Can Ebola Virus evolve to be less virulent in humans?

Keywords: case fatality ratio, virulence, EBOV, transmission, evolutionary epidemiology, public health, trade-offs, outbreak, adaptation Reference

Understanding Ebola Virus (EBOV) virulence evolution is not only timely but also raises specific questions because it causes one of the most virulent human infections and it is capable of transmission after the death of its host. Using a compartmental epidemiological model that captures three transmission routes (by regular contact, via dead bodies and by sexual contact), we infer the evolutionary dynamics of case fatality ratio (CFR) on the scale of an outbreak and on the long term. Our major finding is that the virus's specific life cycle imposes selection for high levels of virulence and that this pattern is robust to parameter variations in biological ranges. In addition to shedding a new light on the ultimate causes of EBOV's high virulence, these results generate testable predictions and contribute to informing public health policies. In particular, burial management stands out as the most appropriate intervention since it decreases the R 0 of the epidemics, while imposing selection for less virulent strains.

Impact Summary

The severe haemorrhagic fever caused by Ebola Virus (EBOV) usually kills more than one infected individual out of two in the absence of treatment, which makes this pathogen one of the most virulent known to humans. The recent outbreak in West Africa (2013-2016) revealed that the virus is able to spread and persist for months across countries. It is often thought that virulence could be due to the fact that the virus is adapted to its reservoir host. Given that microbes evolve rapidly, it is important to determine whether EBOV virulence is likely to decrease as the virus adapts to its human host. To address this problem, we developed a novel mathematical model tailored to EBOV's life cycle, notably by capturing its three main transmission routes (by regular contact, sexual contact and via dead bodies). We investigated the evolutionary trends of EBOV's virulence on different time scales (outbreak initiation, short term and long term). Our results reveal that the virulence of EBOV might not be due to the maladaptation of the virus, but could rather originate from its unique life cycle. These results are robust to the parameter values chosen. From a public health perspective, burial management stands out as the main leverage to fight the virulence of EBOV, both in the short and long terms.

Introduction

Ebola Virus (EBOV) has been a major source of fear since its discovery in 1976. Until 2013, all outbreaks could be traced to spillover from reservoir hosts [START_REF] Leroy | Fruit bats as reservoirs of Ebola virus[END_REF] and were limited in size. This was attributed to EBOV's extremely high case fatality ratio (CFR), that is the ratio of infected hosts who die from the infection, which we use here as a measure of infection virulence. The dramatic 2013-2016 epidemic in West Africa, which caused more than 28,000 infections and claimed at least 12,000 lives, showed that the virus can persist in the human population for months, therefore raising the question: 'How will the virulence of Ebola Virus evolve in humans?' [START_REF] Kupferschmidt | Imagining Ebola's next move[END_REF].

Being an RNA virus, Ebola is prone to rapid evolution [START_REF] De | Ebolavirus evolution: past and present[END_REF] and in vitro analyses suggest that the virus has evolved during the outbreak towards an increased tropism for human cells [START_REF] Urbanowicz | Human adaptation of Ebola virus during the West African outbreak[END_REF]. It was first thought that host-parasite interactions should always evolve towards benignity because mild strains seem to transmit over a larger period of time than strains that kill their hosts rapidly [START_REF] Méthot | Why do parasites harm their host? On the origin and legacy of Theobald Smith's 'law of declining virulence' -1900-1980[END_REF]. Since the 1980s, evolutionary biologists have shown that parasite virulence can be adaptive because it may correlate with transmissibility or within-host competitiveness (for a review, see [START_REF] Alizon | Adaptive virulence evolution: the good old fitness-based approach[END_REF]. The avirulence theory does remain prevalent in many fields. For instance, some envisage a decrease in EBOV virulence due to host adaptation [START_REF] Kupferschmidt | Imagining Ebola's next move[END_REF], even though studies suggest that the virulence of some human infectious diseases such as HIV or tuberculosis has followed an increasing trend since their emergence [START_REF] Gagneux | Host-pathogen coevolution in human tuberculosis[END_REF][START_REF] Herbeck | Is the virulence of HIV changing? A meta-analysis of trends in prognostic markers of HIV disease progression and transmission[END_REF].

Studying virulence as a potentially adaptive trait for the parasite requires encompassing the whole epidemiological life cycle of the infectious agent [START_REF] Alizon | Adaptive virulence evolution: the good old fitness-based approach[END_REF]. In the case of Ebola Virus, most individuals acquire the infection after direct contact with blood, bodily secretions or tissues of other infected humans whether alive or dead [START_REF] Bausch | Assessment of the risk of Ebola virus transmission from bodily fluids and fomites[END_REF]. This post-mortem transmission route cannot be neglected in Ebola Virus epidemics [START_REF] Chan | Ebola Virus Disease in West Africa -No Early End to the Outbreak[END_REF], although its magnitude is likely to be low compared to direct transmission [START_REF] Weitz | Modeling post-death transmission of Ebola: challenges for inference and opportunities for control[END_REF]. From an evolutionary standpoint, this route might be crucial as well since the timing of life-history events can dramatically affect virulence evolution [START_REF] Day | Virulence evolution and the timing of disease life-history events[END_REF]. Intuitively, if the virus is still able to transmit after host death, virulence will have a smaller effect on the parasite's transmission potential. Moreover, there is increasing evidence that EBOV might also be transmitted through sexual contact even long after the clinical 'clearance' of the infection since the virus can persist in the semen for months [START_REF] Eggo | Duration of Ebola virus RNA persistence in semen of survivors: population-level estimates and projections[END_REF][START_REF] Thorson | Systematic review of the literature on viral persistence and sexual transmission from recovered Ebola survivors: evidence and recommendations[END_REF][START_REF] Uyeki | Clinical management of Ebola virus disease in the United States and Europe[END_REF].

Will EBOV become more virulent by adapting to humans? To address this question, we use mathematical modelling to determine how case fatality ratio affects the risk of emergence, how it evolves on the long and in the short term. To this end, we introduce an original epidemiological model that captures all three transmission routes of the virus in human populations. We parametrize our model with data from the well-documented 2013-2016 epidemics. We also perform sensitivity analyses on conservative biological ranges of parameter values to verify the robustness of our conclusions.

We find that EBOV undergoes selection for higher case fatality ratios due to its life cycle that includes transmission from hosts after death. This result is robust to most parameter values within biological ranges. We also show that short-term evolutionary dynamics of virulence are more variable but consistently depend on the duration of the incubation period. Finally, we investigate how public health interventions may affect EBOV virulence evolution. We find a direct, but limited, effect of safe burials that may decrease the spread of the virus, while favouring less virulent strains over more virulent ones.

Methods

For clarity, most of the technical details are shown in Supplementary materials and this section contains verbal description of the model, Figures illustrating the life cycle and a list of parameter values.

Epidemiological model

Our original compartmental model is based on the classical Susceptible-Exposed-Infected-Recovered (SEIR) model, which we enhanced by adding a convalescent class (C) that allows for sexual transmission [START_REF] Abbate | Potential Impact of Sexual Transmission on Ebola Virus Epidemiology: Sierra Leone as a Case Study[END_REF] and an infected dead body class (D) that allows for post-mortem transmission [START_REF] Legrand | Understanding the dynamics of Ebola epidemics[END_REF][START_REF] Weitz | Modeling post-death transmission of Ebola: challenges for inference and opportunities for control[END_REF]. The model is deterministic and does not include additional host heterogeneity, spatial structure or public health interventions.

We incorporated demography through a host inflow term λ and a baseline mortality rate µ. Susceptible individuals (S) can become infected by regular (i.e. not sexual) contact with symptomatic infected individuals (I) (World Health Organization Ebola Response Team, 2014), by sexual contact with convalescent individuals (C) [START_REF] Mate | Molecular evidence of sexual transmission of Ebola virus[END_REF] and by contact with the dead body of Ebola virus disease (EVD) victims, mostly during ritual practices (D) [START_REF] Chan | Ebola Virus Disease in West Africa -No Early End to the Outbreak[END_REF]. The rates at which these events occur are proportional to β I , β C and β D respectively. As in most models [START_REF] Keeling | Modeling infectious diseases in humans and animals[END_REF], we assumed sexual transmission to be frequency-dependent. For non-sexual transmission, we assumed density-dependent transmission following an analysis of the 2013-2016 epidemic [START_REF] Leander | Modelling Ebola within a community[END_REF], although performed at a smaller scale than ours. The total population size of live hosts is denoted N .

Upon infection, susceptibles move to the so-called 'exposed' class (E), meaning they are infected but not yet infectious. For Ebola Virus infections, this latency period is also the incubation period, i.e. the time from the contamination to the onset of the symptoms. These symptoms arise at a rate ω.

At the end of this latency/incubation period, individuals move to the symptomatic infected compartment (I). They leave this compartment at a rate γ, which we calibrated using the average time elapsed from the onset of the symptoms to death or recovery. We hereafter refer to this as the 'symptomatic period'. The probability that an infected individual survives the infection is 1 -α. The case fatality ratio (CFR), α, is our measure of virulence.

We assumed that infected individuals who survive the infection clear the virus from their bloodstream but not from other fluids such as semen and may therefore still transmit EBOV through sexual contacts [START_REF] Deen | Ebola RNA persistence in semen of Ebola virus disease survivors -preliminary report[END_REF]. These convalescent individuals (C) completely eliminate EBOV at a rate σ. Notice that given the severity of the symptoms [START_REF] Feldmann | Ebola haemorrhagic fever[END_REF] and the fact that the convalescence period is one order of magnitude greater than the symptomatic period, we neglected the sexual transmission from symptomatic infected individuals (I).

Based on the current immunological data [START_REF] Sobarzo | Persistent immune responses after Ebola virus infection[END_REF], we assumed that full elimination of EBOV from convalescent hosts confers lifelong immunity to recovered individuals (R), who do not contribute to the epidemiology.

On the contrary, infected individuals who die from the disease may continue to transmit EBOV if their burial is performed in unsafe conditions, which occurs at a probability θ. There is little data from which to estimate this parameter. However, the proportion of EBOV-positive dead bodies has been estimated to decline from 35% to 5% by the end of 2014 in the most populous county of Liberia [START_REF] Nyenswah | Ebola and its control in Liberia, 2014-2015[END_REF]. We therefore set the default value θ = 0.25. In the analysis, we strongly vary this parameter since it represents an important leverage public health policies have on the epidemics. In the absence of burial team intervention, body fluids from dead bodies remain infectious for a period ε -1 which is known to be less than 7 days in macaques [START_REF] Prescott | Postmortem stability of Ebola virus[END_REF].

Our model, pictured in Figure 1, consists in a set of Ordinary Differential Equations (ODEs) shown in Supplementary Material A.

Table 1 lists the 11 model parameters. Their values were calibrated using data from the 2013-2016 epidemic in West Africa. We worked at a country level and preferentially chose estimates from the Liberia outbreak, because with approximately 10,000 cumulative cases (World Health Organization, 2016), its magnitude lies in between that of Sierra Leone and Guinea. Demography data from Liberia were obtained from publicly available data from the Central Intelligence Agency (Central Intelligence Agency, 2016). The newborn inflow was set such that the disease free equilibrium matches the country's population size.

In Supplementary Material C, we calculate the basic reproduction number of EBOV (denoted R 0 ), which is the average number of secondary infections caused by a single infected individual in a fully susceptible population [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations[END_REF]. By studying the local stability of the system (S1) at the disease free equilibrium, we found that

R 0 = β I γ S 0 + αθβ D ε S 0 + (1 -α) β C σ , (1)
where S 0 = λ/µ is the total population size before the epidemic. The three terms in R 0 correspond to each transmission route: from symptomatic individuals (R 0,I := β I S 0 /γ), infectious bodies (R 0,D := αθβ D S 0 /ε), and convalescent individuals (R 0,C := (1 -α) β C /σ). Notice that the incubation period does not affect R 0 .

Transmission-virulence trade-off and estimated values

Trade-offs are a central component of evolutionary model and, without them, predictions tend to be trivial (e.g. viruses should evolve to maximise their transmission rates and minimise their virulence). Although the EBOV life cycle generates constraints that may lead to non-trivial evolutionary outcomes, we do also allow for an explicit trade-off between CFR and transmission rates. Such a relationship has been shown in several host-parasite systems [START_REF] Alizon | Adaptive virulence evolution: the good old fitness-based approach[END_REF]. The case of HIV is particularly well documented [START_REF] Fraser | Virulence and pathogenesis of HIV-1 infection: an evolutionary perspective[END_REF]: viruses causing infections with higher viral loads have increased probability to be transmitted per sexual contact, while causing more virulent (shorter) infections. As a result, there exists an optimal intermediate viral load that balances the virus benefits of transmission with the costs of virulence, thus maximising the number of secondary infections.

In the case of EBOV, there is indirect evidence that viral load is positively correlated with case fatality ratio (CFR) since survivor hosts tend to have lower viral loads than non-survivors [START_REF] Towner | Rapid diagnosis of Ebola hemorrhagic fever by reverse transcription-PCR in an outbreak setting and assessment of patient viral load as a predictor of outcome[END_REF][START_REF] Crowe | Prognostic indicators for Ebola patient survival[END_REF]. Viral load is thought to correlate with transmission [START_REF] Osterholm | Transmission of Ebola viruses: what we know and what we do not know[END_REF] but demonstrating a clear link is challenging (for HIV, it has required identifying sero-discordant couples in cohorts).

We assumed an increasing relationship between transmission rates and CFR (α) such that:

β H (α) := b H α p , ( 2
)
where b H is a constant factor, p ≥ 0 is a parameter capturing the concavity of the trade-off curve and H stands for any of the compartment I, D or C. The exact value of p results from within-host interactions [START_REF] Alizon | Emergence of a convex trade-off between transmission and virulence[END_REF] but one can identify four kinds of trade-offs: p > 1 corresponds to an amplifying transmission benefit from increasing CFR, p = 1 corresponds to a linear relationship between transmission and CFR, 0 < p < 1 corresponds to a saturating transmission benefit from increasing CFR and p = 0 is the degenerate case without trade-off. From a biological standpoint, we could expect different transmission routes to have different trade-off shapes (p) but, as we show here, our results are largely independent of p. Transmission rates being difficult to estimate [START_REF] Leander | Modelling Ebola within a community[END_REF], we indirectly infer the order of magnitudes of b I , b C and b D by setting the R 0 close to 2, that is its approximate value for the 2014 epidemic (World Health Organization Ebola Response Team, 2014). Since R 0 is the sum of the epidemiological contributions of each transmission route (see equation ( 1)), we added the constraint that, according to previous studies, transmission from symptomatic individuals contributes about ten times more than transmission from dead bodies (World Health Organization Ebola Response Team, 2014; [START_REF] Weitz | Modeling post-death transmission of Ebola: challenges for inference and opportunities for control[END_REF] and one hundred times more than transmission from convalescent individuals [START_REF] Abbate | Potential Impact of Sexual Transmission on Ebola Virus Epidemiology: Sierra Leone as a Case Study[END_REF]. Straightforward calculations (analogous to those done for sensitivity analysis in Supplementary Material E) resulted in attributing the orders of magnitude shown in Table 1.

In the following, the exponent p was left undetermined, but its value was set to 0 for the estimation of b H in the null hypothesis. This leads to R 0 ≈ 1.86, which is very close to the WHO mean estimation for the Liberia epidemic, namely 1.83 (World Health Organization Ebola Response Team, 2014).

Long term evolution

We used the adaptive dynamics framework [START_REF] Geritz | Evolutionarily singular strategies and the adaptative growth and branching of the evolutionary tree[END_REF], which assumes that evolution proceeds by rare and small phenotypic mutations occurring in a monomorphic population that has reached ecological stationarity. Unless the fitness landscape exhibits a branching point (which is not the case here), polymorphism is therefore limited to transient dimorphic phases where the ancestor (hereafter called the 'resident') and its mutant compete. Depending on the outcome of the competition, the system either goes back to the previous endemic equilibrium or reaches a new monormophic equilibrium. The corresponding dynamical system is shown in Supplementary Material A (system (S2) applied to n = 2).

Given the focus of this study, we assumed that the resident and the mutant only differ in their CFR (and in their transmission traits if there was a transmission-virulence trade-off). We denoted the virulence of the mutant and the resident by α and α respectively. α was assumed to be close to α. Since the mutant is rare by definition, its emergence can be assumed not to affect the resident system. We therefore investigated the local stability of the related endemic equilibrium. This only depended on the local stability of the mutant sub-system (E 2 , I 2 , D 2 , C 2 ) to which we applied the next-generation theorem [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations[END_REF][START_REF] Hurford | Next-generation tools for evolutionary invasion analyses[END_REF]. This eventually led (see Supplementary Material C.3) to the mutant relative reproduction number

R α , α = β I (α ) γ + xθβ D (α ) ε + (1-α )β C (α ) σ N (α) β I (α) γ + αθβ D (α) ε + (1-α)β C (α) σ N (α) , ( 3
)
where there total host population size can be approximated (see Supplementary Material B) by

N (α) ≈ (1 -α) S 0 + α β I (α) γ + αθβ D (α) ε . (4)
We then calculated the selection gradient by differentiating the relative reproduction number (equation (3)) with respect to the mutant's trait [START_REF] Otto | A biologist's guide to mathematical modeling in ecology and evolution[END_REF]. Equating the mutant and resident trait value we eventually found

∆ (α) = p α + θb D ε N (α) -b C σ (1 -α) b C σ + b I γ + θb D ε α N (α)
.

(5)

The sign of ∆ indicates the direction in which natural selection acts on the trait depending on the resident's trait.

Short term evolution

Viruses evolve so rapidly that evolutionary and epidemiological dynamics may overlap. The epidemiological Price equation framework is designed to predict short-term evolution based on standing genetic variation [START_REF] Day | A general theory for the evolutionary dynamics of virulence[END_REF][START_REF] Day | Insights from Price's equation into evolutionary epidemiology[END_REF]. Practically, we assumed that the parasite population is initially diverse and consists of n different genotypes, each genotype i being defined by a specific value for several phenotypic traits of interest: the case mortality (α i ), the rate of end of latency (ω i ), the rate of end of the infectious period (γ i ), the rate at which dead bodies cease to be infectious (ε i ), the rate at which convalescent hosts clear the infection (σ i ) and the transmission rates (β D,i , β I,i and β C,i ).

The dynamics of the populations of interest are described by 4n + 1 ODEs that are shown in Supplementary Material A.

After thorough calculations (in Supplementary Material F) we find that, if we neglect mutational bias, mean traits in each compartment vary according to a system of ODEs that involves statistical covariances and variances of traits in the different compartments. The equations involving average CFR are shown in the Results section.

An important assumption of this Price equation approach is that covariance terms are assumed to be constant, which implies that predictions are only valid in the short term. Given that the main limitation of the adaptive dynamics framework relies in its assumption that epidemiological dynamics are fast compared to evolutionary dynamics, combining the two frameworks allows us to get a broader picture.

Results

Virulence and emergence

We first consider the risk for an epidemic to occur as a function of EBOV virulence, trade-off shape and burial management. A disease spreads in a fully susceptible population if its reproduction number (R 0 ) is greater than unity [START_REF] Anderson | The Population Dynamics of Microparasites and Their Invertebrate Hosts[END_REF]. Using our default parameter values (Table 1), we show in Figure 2 that the most virulent EBOV strains are almost always the most likely to emerge, independently of the trade-off exponent p and of proportion of unsafe burials θ. To have R 0 decrease with α, one needs to have neither trade-off nor unsafe burials (θ = p ≈ 0). However, with our default parameter values this decrease is limited.

If we focus on the lowest virulence that may lead to an epidemic (denoted α min ), we find that with our default parameter values burial management can prevent emergence (that is bring R 0 below unity by moving vertically in Figure 2) only if the transmission-virulence trade-off is strong enough (the green, blue and cyan curves).

In the following, we will generally assume that EBOV is adapted enough to persist in the human population (R 0 > 1). Since outbreak originates via spillover from reservoir hosts [START_REF] Leroy | Human asymptomatic Ebola infection and strong inflammatory response[END_REF], it is likely that the virus is maladapted in the first human infections. However, to capture these dynamics, an evolutionary rescue formalism would be more appropriate given the importance of stochastic events and this is outside the scope of this study (for a review, see [START_REF] Gandon | What limits the evolutionary emergence of pathogens?[END_REF].

Long-term virulence evolution

If the selection gradient in equation ( 5) is negative for any CFR (∆(α) < 0), then the virus population evolves towards its lowest virulence that allows persistence (that is α min ). If the gradient is always positive (∆(α) > 0), the CFR evolves towards 1. Intermediate levels of virulence can only be reached if there exists α such that ∆(α) ≥ 0 for α ≤ α , ∆(α) ≤ 0 for α ≥ α and R (α, α ) < 1 for any α = α . We show in Supplementary Material D.4 that this occurs only if the proportion of unsafe burials (θ) and the trade-off parameter (p) are lower than the following boundaries

θ < b I b C ε γσb D and p < b C σ , ( 6
)
Unless these two conditions are met, the selection gradient is always positive and EBOV is expected to always evolve towards higher case fatality ratios (α = 1). Rewriting the first inequality as

θb D S 0 ε < b I S 0 γ × b C
σ highlights that virulence is favoured by natural selection as soon as the post mortem transmission component is greater than the product of the symptomatic and convalescent transmission components.

Figure 3 shows how α is numerically affected by a change in burial management (θ) and trade-off strength (p). Unless the proportion of safe burials is brought below 4%, and unless there is a weak trade-off (p < 0.01), CFR will remain high. Intuitively, this double condition can be understood in the following way. If the trade-off is negligible, the CFR is weakly linked to transmission by regular contact and therefore selection on α only weakly depends on this component of the life cycle. As a consequence, the value of θ governs the relative importance of the two transmission routes that matter. Post mortem transmission always favours higher CFR, whereas the sexual transmission route can be maximised for intermediate levels of virulence (see Supplementary Material H.)

It was not possible to find an explicit expression for the long-term equilibrium virulence (α ), but we found it lies in the following interval (Supplementary Material D.5):

α ∈      p b C σ - θ b D ε b I γ +θ b D ε , b I γ + θ b D ε S 0 + b C σ p (1 + p) b C σ -θ b D ε S 0      . ( 7
)
The lower bound of this interval increases with trade-off strength (p) and intensity of the post mortem transmission route (θb D /ε). If post mortem transmission is strongly reduced, owing to a safer burial management (θ → 0), the lower bound simplifies to pσ/b C . The long-term virulence then appears to be a balance between trade-off strength and sexual transmission intensity, which is consistent with our intuitive explanation of the condition for an intermediate virulence to be selected. In particular, any decrease in the time for convalescent hosts to clear the virus (i.e. increase in σ) will increase the lower bound for CFR.

To further assess the robustness of these results, we performed a sensitivity analysis by varying the relative importance of each transmission route (regular contact, sexual transmission and transmission from dead bodies), while keeping the total value of R 0 constant. As shown in Figure 4, unless the values of p are extremely low, variations in the relative transmission routes is unlikely to be sufficient to move our default value (dotted lines) to the region where low virulences (e.g. α < 50%) are favored (green area). To give a quantitative estimate, the relative importance of transmission via sexual contact (on the vertical axis) would need to be about 40 times greater than the current estimates to bring the current estimate above the blue separatrix, which already assumes a low trade-off and a perfect burial management.

Short term evolutionary dynamics

Reaching an evolutionary equilibrium may take time (especially if strains have similar trait values) and the transient dynamics can be non-trivial because the system is non-linear. The Price equation framework provides a way to qualitatively address the initial trends of average trait values, by considering the initial diversity in the virus population.

If we denote by x H the average value of trait x in compartment H and by cov H (x, y) the statistical covariance between traits x and y in compartment H (which becomes the statistical variance var H (x) if x ≡ y), the dynamics of average virulence in the four infected compartments satisfy the following ODEs (see Supplementary Material F for further details):

dα I dt = -cov I (α, γ) + cov E (α, ω) + α E -α I ω E E • I • , ( 8a
) dα E dt = -cov E (α, ω) + S E • H∈{I,D,C} cov H (β H , α) + α H -α E β H H H • , (8b) dα D dt = -cov D (α, ε) + var I (α) + α I α I -α D θγ I • D • , ( 8c
) dα C dt = -cov C (α, σ) + cov I (α, γ) -γvar I (α) + 1 -α I α I -α C γ I I • C • , ( 8d
)
where

H • := n i=1 H i denotes the total size of compartment H ≡ E, I, D, C.
Focusing on the compartment on which virulence acts, namely the symptomatic individuals, indicates that the short-term evolution of the average virulence α I is mainly governed by the correlations between this trait and the symptomatic and latency periods. More explicitly, equation (8a) states that if the most virulent strains induce the longest symptomatic period and/or the shortest latency periods, the average virulence in I can be expected to increase at the beginning of the epidemic. Intuitively, newly symptomatic individuals are more likely to have been infected with a highly virulent strain.

Equation (8a) contains a third term proportional to α E -α I , which is more difficult to apprehend. Indeed, α E varies as well and follows a complicated ODE that involves not only the correlation with the latency period but also correlations with the transmission rates (equation ( 8b)). This diversity of components make both α E and α I difficult to predict early in the epidemics.

We therefore simulated epidemics numerically according to system (S2). We considered nine scenarios of increasing complexity, three of which are shown in Figure 5 (see Supplementary Material G for more details). During the first six months of an (unmanaged) epidemic, average virulence exhibits wide variations. In most scenarios (panels A and B), it tends to evolve towards the maximum of the range provided by its initial polymorphism. Transient evolution for further virulence in an expanding epidemic have been highlighted by previous models (Day and Proulx, 2004, e.g.) and studies [START_REF] Berngruber | Evolution of virulence in emerging epidemics[END_REF] but this effect was due to a positive correlation between virulence and transmission rates, which is here not required (Figure 5A). Unsurprisingly, the addition of such correlation amplifies the transient increase in virulence observed during the approximately first 300 days after the onset of the epidemic (Figure 5B).

A scenario where average virulence decreases initially is when it is positively correlated with the latency period (Figure 5C). This occurs because less virulent strain have an advantage early in the epidemics by reaching the infectious class earlier. More virulent strains become more frequent again when the value of D begins to take off (Figure 5B).

A secondary result shown in these figures is that dead bodies (in brown) always carry more virulent strains on average.

Discussion

Virulence could be adaptive for EBOV Ebola Virus is one of the deadliest human pathogen [START_REF] Feldmann | Ebola haemorrhagic fever[END_REF]. The recent epidemic in West Africa has shown that it can transmit for months among humans throughout entire countries. As with any microbe (especially RNA viruses), it is likely exposed to fast evolution during the course of the epidemic [START_REF] Bedford | Did a single amino acid change make Ebola virus more virulent?[END_REF]. From a public health standpoint, it is important to predict Ebola virus' next move and the current hope is that the shift from an emerging to an endemic lifestyle could favour less virulent strains [START_REF] Kupferschmidt | Imagining Ebola's next move[END_REF].

Predicting virulence evolution is usually challenging because we tend to lack details about correlations between epidemiological parameters. Furthermore, even when there is data to estimate a trade-off relationship, its exact shape can have huge quantitative and even qualitative effects on the evolutionary dynamics of the trait [START_REF] Alizon | Emergence of a convex trade-off between transmission and virulence[END_REF][START_REF] Svennungsen | Evolutionary branching of virulence in a single infection model[END_REF]. Our results stand out because they are robust both to parameter variation in wide biological ranges (World Health Organization Ebola Response Team, 2014Team, , 2015) ) and also to the type of trade-off assumed. Importantly, the numerical analysis of the model show that our results still hold even if there is no transmission-virulence trade-off at all, as long as the burial management is not perfect (θ > 5%).

In addition to the strong selection on EBOV virulence due to its transmission via dead bodies, another striking result is that decreasing the ratio of unsafe burials is triply effective. First, it decreases the spread of the virus (i.e. its R 0 ). Second, in the short term, it can help limit a transitory increase in virulence. Indeed, a better burial management prevents the spread of the strains hosted by infected corpses, which is the compartment that concentrates most of the virulent strains (as shown by Figure 5). Third, in the long term, decreasing the proportion of unsafe burials is necessary to shift the selective pressure in favour of less virulent strains.

Overall, EBOV is unlikely to evolve to become less virulent because that would require two conditions. First, the proportion of unsafe burials must be brought to a very low value, which we estimate to be lower than 4%. Second, there must be very little or no genetic relationship between EBOV case fatality ratio and transmission rate. This latter condition is particularly frustrating because it cannot directly be addressed by public health policies. Finally, even if these conditions are met, the level of virulence reached in the long term may still be high, especially if sexual transmission is limited. On a more positive note, results from the Price equation approach show that the virus may experience transitory lower levels of virulence before reaching this maximum via a positive genetic correlation between virulence and incubation period. This is somehow unexpected because this latter parameter does not appear in the calculations related to long-term evolutionary or emergence.

Virulence is a shared trait

In evolutionary biology, virulence is usually defined as the decrease in host fitness due to the infection [START_REF] Read | The evolution of virulence[END_REF][START_REF] Alizon | Adaptive virulence evolution: the good old fitness-based approach[END_REF][START_REF] Cressler | The adaptive evolution of virulence: a review of theoretical predictions and empirical tests[END_REF]. Given the speed at which a pathogen kills its host, EBOV's virulence can neither be measured as a decrease in instantaneous fecundity (which would be almost zero) nor as an increase of instantaneous mortality rate (which would tend towards infinity or zero depending on the infection outcome). The case fatality ratio (CFR) therefore appears to be a convenient measurable and epidemiologically relevant proxy for EBOV's virulence.

We focused on the virus side but, like any infection trait, virulence is also determined by the host and its environment. To predict how virulence will change in the future, we should also consider how hosts may change. In the case of EBOV, it was known before the recent epidemics that some people can become immune to the virus without exhibiting any symptoms [START_REF] Leroy | Human asymptomatic Ebola infection and strong inflammatory response[END_REF]. The question remains to know if they can also be infectious. More generally, our assumption of life-long protection could be oversimplifying.

Finally, to make predictions on the long term evolution, we need to factor in how the virus population will evolve in response to the variation in the host population's immune status. Since the immunological status of the host population is determined by that of the virus population, this de facto qualifies as a coevolutionary interaction. Earlier models shows that host heterogeneity in resistance and tolerance can lead to a variety of outcomes [START_REF] Miller | The evolution of parasites in response to tolerance in their hosts: the good, the bad, and apparent commensalism[END_REF][START_REF] Cousineau | Parasite evolution in response to sex-based host heterogeneity in resistance and tolerance[END_REF]. Overall, introducing realistic host heterogeneity, in particular age-dependent or sex-dependent mortality, appears like a relevant extension of this model.

Symptomatic period and lethality

The Price equation framework allows us to study any statistical association between epidemiological traits even if the shape of the dependency is unknown. On the contrary, the derivation of the basic reproduction number as well as the adaptive dynamic framework require the precise expression of epidemiological traits as functions of the focal trait. We therefore consider possible genetic variation on the infectious period 1/γ when studying the short-term evolutionary dynamics on one hand, while assuming that 1/γ is independent from the CFR α when investigating emergence and long-term evolution, on the other hand.

To our knowledge, clinical data allowing to model the variations of rate at which symptomatic period ends (either due to recovery or EVD) γ in function of α is still lacking. Intuitively, the most virulent EBOV strains should proliferate faster than average and therefore be more likely to kill their host, and to kill it more quickly, thus reducing the symptomatic period. This would yield a positive correlation between α and γ. However, it is also possible to envisage the opposite trend: less virulent EBOV strains are more likely to be cleared by the immune system, and to be cleared more quickly, thus reducing the symptomatic period as well. Such a trend would yield a negative correlation between α and γ. Therefore, the assumption that γ is independent from α should not be seen as a pure mathematical simplification but rather as a parsimonious modelling choice driven by the lack of relevant data.

Even in the absence of detailed data, we can further speculate on how allowing the duration of the symptomatic phase to vary could affect our results. Indeed, clinical studies have shown that nonsurvivors from EVD undergo faster disease course, associated with higher RNA copy levels [START_REF] Towner | Rapid diagnosis of Ebola hemorrhagic fever by reverse transcription-PCR in an outbreak setting and assessment of patient viral load as a predictor of outcome[END_REF]. Since the EVD pathogenesis in humans is primarily driven by the impairment of the immune system [START_REF] Hoenen | Ebola virus: unravelling pathogenesis to combat a deadly disease[END_REF][START_REF] Falasca | Molecular mechanisms of Ebola virus pathogenesis: focus on cell death[END_REF], which is consequently unable to stop the replication of the virus, EBOV's lethality can also be seen as a massive reduction in the host's ability to recover. Therefore, if we denote by m the EVD-induced mortality rate, by r the recovery rate and x a withinhost trait that characterizes the EVD course for a given EBOV strain (such as replication rate or host exploitation), the end of symptomatic period rate is γ (x) = m + r (x). Since the CFR is equal to the relative contribution of EVD-induced mortality to the rate at which the symptoms end, it follows that α (x) = m/ (m + r (x)) = m/γ (x). Therefore, the symptomatic period 1/γ can be replaced by α/m in (1) and (3), which in the end increases the contribution of the regular contact transmission route to the reproduction numbers, leading to even greater selection for higher virulence.

Spatial structure

Trait evolution is shaped by contact patterns between hosts [START_REF] Lion | Self-structuring in spatial evolutionary ecology[END_REF]. Regarding the recent Ebola epidemic, the lack of medical personnel and infrastructure in the affected countries played an key role in the spread of the disease as, for example, according to the World Health Organisation, in 2008 Liberia and Sierra Leone had only a density of 0.015 physicians per 1000 inhabitants, when at the same time France had a density of 3.5 and the United States of America 2.4. This was further exacerbated by historical, political and sociological factors [START_REF] Ali | The Social and Political Dimensions of the Ebola Response: Global Inequality, Climate Change, and Infectious Disease[END_REF].

It is difficult to predict how explicitly accounting for spatial structure would affect the results. Indeed, it is generally thought that the more 'viscous' the population, the more virulence is counter-selected [START_REF] Boots | Small worlds' and the evolution of virulence: infection occurs locally and at a distance[END_REF]. However, the life cycle of the parasite and the host demography can create epidemiological feedbacks that alter this prediction by causing virulence to be maximised for intermediate levels of population structures [START_REF] Lion | Are parasites 'prudent' in space?[END_REF]. This is why predicting virulence evolution in a fully spatially structured model constitutes a challenging future work.

Inputs from data

One of the underlying assumptions of our model, which could be tested is that the variation we observe in virulence is at least partially controlled by the virus genetics. This could be done by combining virus sequence data with infection traits (virus load or infection outcome) through a phylogenetic comparative approach [START_REF] Alizon | Phylogenetic approach reveals that virus genotype largely determines HIV set-point viral load[END_REF] or a genome wide association study on the virus genome [START_REF] Power | Microbial genome-wide association studies: lessons from human GWAS[END_REF]. If virus load is confirmed to be partially controlled by the virus genetics and if, as current evidence suggests, it is correlated with virulence [START_REF] Towner | Rapid diagnosis of Ebola hemorrhagic fever by reverse transcription-PCR in an outbreak setting and assessment of patient viral load as a predictor of outcome[END_REF][START_REF] Crowe | Prognostic indicators for Ebola patient survival[END_REF], then studying variations in virus load throughout the 2013-2016 epidemics can help us understand the evolutionary dynamics of virulence. On that note, an experiment consisting in generating pseudovirions based on ancestral or recent EBOV sequences suggests that some of the substitutions observed during the 2013-2016 epidemics may confer increased tropism for human cells [START_REF] Urbanowicz | Human adaptation of Ebola virus during the West African outbreak[END_REF].

Another result of the short-term evolutionary dynamics analysis is that individuals who contract EBOV from dead bodies should have a higher probability of dying than those infected by contact with living infectious individuals. This could be tested by collecting data from individuals where the transmission route is well documented.

Finally, the remote possibility that lower virulence strains will evolve depends on the existence of a transmission-virulence trade-off. Assessing the shape of this trade-off may be, therefore, very valuable. Note that in the case of EBOV, it is not the exact shape that matters but rather the general trend.

Conclusion

This evolutionary epidemiology work shows that EBOV's high virulence, whether it is about emergence, short-term or long-term dynamics, can be explained by its particular life cycle that mixes parasitism and parasitoidism (post mortem transmission). Unfortunately, any long term decrease in virulence is unlikely for West African strains at any time scale, although increasing the safe burial proportion appears to be an optimal response in both the short and long terms. The solid, dashed and dotted lines correspond to α = 1, 0.7 and 0.3 respectively. Parameter vlaues are shown in Table 1. See Supplementary Material D.4 for more details.

Figure 4: Sensitivity analysis of long-term virulence evolution.

The graphic shows the sign of the selection gradient for virulence when varying the relative weight (in orders of magnitude) of the post-mortem transmission component and the sexual transmission component in the overall transmission of EBOV. When the basic reproduction number is set at its upper bound ( R 0 = 2.53, dashed line), the selection gradient at the maximum virulence (α = 1) is positive below the dashed line (red area). When the basic reproduction number is set at its lower bound ( R 0 = 1.26, plain line), the selection gradient is also positive for a range of virulence higher than one half (α ≥ 50%) in the orange area. It is negative for lower virulences (α < 50%) above the dashed line (green area). The unsafe burial proportion and the trade-off exponent are low (θ = 0 and p = 0.1). See Supplementary Material E for more details.

Figure 5: Short-term evolution of CFR with standing genetic variation in three scenarios.

A) Without correlations between traits, B) with a positive correlation between CRF and transmission rate and C) with a negative correlation between CRF and latency period. The CFR averaged over the exposed individuals (α E ) is depicted in cyan, over the symptomatic individuals (α I ) in pink and over the infectious dead bodies (α D ) in brown. Parameter values are shown in Table 1. See Supplementary Material G for details about the simulations.

Hence,

dE dt = 0 ⇐⇒ β I I + β D D + β C C N S -ω E = 0, ⇐⇒ β I ω γ E + β D αθω ε Ẽ + β C (1 -α) ω σ N E S -ω E = 0, ⇐⇒ S = β I γ + αθβ D ε + (1 -α) β C σ N -1 , ( S5
)
and

dS dt = 0 ⇐⇒ λ -µ + β I I + β D D + β C C N S = 0, ⇐⇒ µ + β I ω γ E + β D αθω ε Ẽ + β C (1 -α) ω σ N E = λ S , ⇐⇒ ω E S = λ S -µ, ⇐⇒ E = λ -µ S ω . ( S6
)
It follows that

N := S + E + I + C + R = S + 1 + ω γ + (1 -α) ω σ + (1 -α) ω µ E, = S + 1 ω + 1 γ + 1 σ + 1 µ (1 -α) λ -µ S . (S7)
By combining (S5) and (S7), we can find the exact solution for N . This closed form is excessively large and therefore not shown here. It is however possible to find a approximation of N as a simple function of the model's parameters with some simplifications that are shown hereafter, with a particular treatment of the α = 1 case. the third term in equation (S5), which leads to

S ≈ γε εβ I + αγθβ D . (S9) This then results in N ≈ (1 -α) S 0 + αγε εβ I + αγθβ D . (S10)
Numerical comparisons performed on positive and stable equilibria for realistic parameter sets show that this approximation differs from the exact value by less than 10,000 individuals, which corresponds to a relative error of less than 1%, thus validating the accuracy of this approximation.

B.3 Stationary total population size approximation for α = 1

Here we assume that α = 1 (notice that in this case the trade-off exponent p vanishes).

We then get back to equation ( S5) that becomes such that

S = β I γ + θβ D ε -1 = γε εβ I + γθβ D ,
which shows that the exact value of S coincides with its equation (S9) approximation for α = 1. As for equation (S7), we get

N = 1 - (γ + ω) µ γω S + (γ + ω) λ γω ,
it is straightforward to show numerically (using parameters from Table 1) that, since S and λ µ = S 0 are of the same order of magnitude and that, as already mentioned, life expectancy is several orders of magnitude greater than the latency and symptomatic period, i.e.

1 µ 1 ω + 1 γ , which is equivalent to (γ+ω)µ γω 1, we have N ≈ S = γε εβ I + γθβ D , (S11)
which shows the consistency of equation (S10) even for α = 1. This approximation shows a relative error of about 10 -3 with Table 1 parameters values.

C Reproduction number derivation

The basic reproduction number, R 0 , and the relative reproduction number, R m , are two epidemiological quantifications of the invasion potential of an infectious agent in a fully susceptible population and in an population already infected by an alternative strain, respectively. They emerge from the stability analysis of the disease free equilibrium (DFE) and the endemic equilibrium (EE) respectively. Their threshold value is 1.

The next-generation method [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations[END_REF][START_REF] Hurford | Next-generation tools for evolutionary invasion analyses[END_REF] is the most efficient derivation of these reproduction numbers and proceeds as follows.

C.1 General reproduction number

First, we isolate the ODEs of the infected compartments from the rest of the system (here system (S1)) to obtain

             dE dt = β I I + β D D + β C C N S -ωE, dI dt = ωE -γI, dD dt = αθγI -εD, dC dt = (1 -α) γI -σC. (S12)
Second, we write the Jacobian matrix J that corresponds to this sub-system (S27), by deriving each time-derivative ( dE dt , dI dt , dD dt , dC dt ) with respect to each infected compartment density (E, I, D, C):

J =      -β C CS N 2 -ω β I + β C C N 2 S β D S 1 -C N β C S N ω -γ 0 0 0 αγθ -ε 0 0 (1 -α) γ 0 -σ      , reminding that N := S + E + I + R + C.
Third, we arbitrarily decompose J as a sum of an 'inflow' matrix F and an 'outflow' matrix -V provided that V is non-singular (that is V -1 exists), F and V -1 are non-negative elementwise and the real parts of all eigenvalues of -V are negative. Here, we conveniently choose two matrices that fulfil these requirements:

F =      0 β I + β C C N 2 S β D S β C S N 0 0 0 0 0 0 0 0 0 0 0 0      and V =      -β C CS N 2 -ω 0 0 -β C CS N 2 ω -γ 0 0 0 αγθ -ε 0 0 (1 -α) γ 0 -σ      .
Finally, the general reproductive number R is given by the largest modulus of all eigenvalues of the F.V -1 matrix. Elementary calculations result in the following general result

R = ((1 -α) γεβ C + (εβ I + αγθβ D ) N ) ωSN (γσβ C CS + ((1 -α) γβ C CS (β C CS + N 2 γ) σ) ω) ε . (S13)

C.2 Basic reproduction number

The basic reproduction number R 0 is obtained from R by setting the densities to their values at the disease free equilibrium (DFE), namely S = N = λ µ and C = 0, hence

R 0 = β I γ + αθβ D ε S 0 + (1 -α) β C σ , ( S14
)
Any strain introduced in a fully susceptible host population spreads if and only if R 0 > 1.

C.3 Relative reproduction number

As for the relative reproduction number R m , it is obtained from R by setting the densities to their values at an alternative strain endemic equilibrium (EE), namely S = S, N = N and C = 0 (notice that in such setting N = S + R + E r + I r + C r + E + I + C where the r index denotes compartments of individuals infected by the previously established ('resident'), which may not be empty at EE, making S < N ), hence

R m = β I γ + αθβ D ε + (1 -α) β C σ N S.
It follows that a rare mutant strain of CFR x that appears in a host population endemically infected by a resident strain of CFR y spreads and persists if and only if

R (x, y) := β I (x) γ + xθβ D (x) ε + (1 -x) β C (x) σ N (y) S (y) > 1. (S15)
Moreover, it is possible to eliminate S (y) using equation (S5), leading to

R (x, y) = β I (x) γ + xθβ D (x) ε + (1-x)β C (x) σ N (y) β I (y) γ + yθβ D (y) ε + (1-y)β C (y) σ N (y)
.

(S16)

This formula shows in particular that, because of the two occurrences of N (y), the relative reproduction number is not the ratio between the two basic reproduction numbers, as it is in simpler models [START_REF] Dieckmann | Adaptive dynamics of pathogen-host interactions[END_REF].

We can finally apply approximation (S10) N (y) ≈ (1 -y) S 0 + γεy εβ I (y)+γθyβ D (y) to obtain a closed-form expression for R (x, y) (not shown here).

D Evolutionary analysis of virulence

Investigating the evolutionary trends require to consider trade-offs. From now on, we will then always apply the transmission-virulence trade-off assumed in equation ( 2) and keep in mind that the β H constant case can be retrieved if p = 0.

D.1 Virulence effect on basic reproduction number

It is worth noticing that unless p = 0, we have R 0 = 0 when α = 0. Therefore not all CFR/virulence levels are able to give rise to an epidemic and persist in the population. Indeed, R 0 (α) may not be greater than 1 for all α ∈ [0; 1]. First, let us then study how R 0 is affected by α, by calculating its derivative

dR 0 dα (α) = d dα b I γ + αθb D ε S 0 + (1 -α) b C σ α p , = θb D S 0 ε - b C σ α + p b I γ + αθb D ε S 0 + (1 -α) b C σ α p-1 , which cancels only if α = 0 or α = b C σ + b I S 0 γ b C σ -θb D S 0 ε × p 1 + p =: α • , which lies in ]0, 1[ if and only if θ < εb C σb D S 0 ≈ 2.3% and p < b C σS 0 - θb D ε b I γ + θb D ε ≈ θ=0 5.6 • 10 -3 (numerical values are
given according to Table 1 calibration). If these conditions are not fulfilled, then dR 0 /dα is positive for all CFR. Given these conditions, the value of α • can be approximated by

α°≈ θ→0 σb I S 0 γb C p.
The basic reproduction number at this value is

R 0 (α°) = b I γ + αθb D ε S 0 + (1 -α) b C σ σb I S 0 γb C p p ≈ θ→0 b I S 0 γ , which is a maximum on ]0, 1[ (inequality d 2 R 0 dα 2 (α • ) < 0 has been checked after calculations not shown). Besides, evaluating R 0 for α = 1 leads to R 0 (1) = b I γ + θb D ε S 0 ≥ b I S 0 γ ,
the lower bound being greater than one according to Table 1 estimates, and this holds even with almost half of the b I value and smaller values of γ.

To conclude, for any given values of p and θ, there is always a CFR interval [α min , 1] in which any strain can spread.

Notice also that for α = 0 and p = 0,

R 0 = b I S 0 γ + b C σ ≥ b I S 0 γ ,
and likewise this is greater than one for estimated parameters. Consequently, all CFR values can spread in absence of trade-off.

D.2 Minimum spreadable CFR approximation

α min ∈ [0, 1] is the minimum CFR of EBOV required to spread, i.e. R 0 (α min ) := 1. However, it is not possible to find the exact closed form of α min (as the equation R (x) = 1 involves irreducible terms of both x p and x). It is nonetheless possible to analytically find a lower bound for α min , which we denote by α -(0 ≤ α -≤ α min ≤ 1). First, notice that

R 0 (α) := β I γ + αθβ D ε S 0 + (1 -α) β C σ ≤ β I γ + θβ D ε S 0 + β C σ =: R 0,+ (α) , ( S17
)
where R 0,+ (α) is an over-estimate of R 0 (α). Applying the trade-off from equation (2), we get

R 0,+ (α) = b I γ + θb D ε S 0 + b C σ α p .
Let α -be the CFR such that R 0,+ (α -) = 1, that is

α -= b I γ + θb D ε S 0 + b C σ -1 p . ( S18
)
From equation (S17) and as R 0 and R 0,+ are functions of α, it follows that

R 0 (α -) ≤ R 0,+ (α -) = 1 = R 0 (α min ) ≤ R 0,+ (α min ) ,
thus α -is an analytical under-estimate of α min .

Notice that in absence of trade-off, the closed form α min can be easily found as

α min = 1 -b I S 0 γ -b C σ θb D S 0 ε -b C σ . (S19)

D.3 Virulence effect on stationary densities

From now one, we will assume that α ∈ [α min , 1]. We apply definition from equation (2) to equation (S9), that is

S ≈ γε (εb I + γθb D α) α p . Its derivative with respect to α is d S dα (α) ≈ -(γθb D α + p (εb I + γθb D α)) γεα p-1 (εb I + γθb D α) 2 < 0.
Thus, it comes from equation (S8), that

d N dα (α) = -S 0 + S (α) + α d S dα (α) .
Since S (α) < S 0 for α ≥ α min (any EBOV strain that spreads decreases the number of susceptible individuals), d N dα (α) < 0, that is N (α) is a decreasing function of α, we have

N (α) ≥ N (1) = γε εb I + γθb D .

D.4 Virulence selection gradient

We can finally investigate the virulence selection gradient, ∆, by deriving the relative reproduction number from equation (S16) with respect to the first argument (i.e. the mutant's virulence), which leads to ∂ 1 R, and equalizing the mutant and resident's virulence. After some calculations, we find that

∆ (y) := ∂ 1 R (y, y) = p y + θb D ε N (y) -b C σ (1 -y) b C σ + b I γ + θb D ε y N (y)
.

Since the CFR is bounded by 1, it is expected to evolve towards lower values if and only if ∆ (1) < 0, that is ∆ (1)

≈ p + θb D γ εb I + γθb D - b C σ < 0.
By investigating burial control under the most favourable trade-off, which is no trade-off (p = 0), we find that this condition is equivalent to (for b C σ):

θ < b I b C ε γσb D ≈ 4%.
By investigating trade-off shape under the most favourable burial control (θ = 0), we find that the condition is equivalent to

p < b C σ ≈ 10 -2 .
Moreover, investigating the selection gradient at the lowest spreadable CFR y = α min , we notice that the following lower bound

∆ (α min ) ≥ θb D ε N (1) -b C σ (1 -α min ) b C σ + b I γ + θb D ε α min N (1) is positive if θ > b I b C ε γσb D .
Therefore a necessary condition for having a negative selection gradient on the lowest spreadable CFRs is θ < b I b C ε γσb D .

D.5 Evolutionary attracting virulence estimation

Unless it is equal to the 0 or 1 boundaries (the determination of which only requires the invariant sign of the selection gradient), the evolutionary attracting virulence α is an intermediate CFR value in ]0, 1[ such that ∆ (α ) = 0 (singularity condition), ∂ 1,1 R (α , α ) < 0 (evolutionary stability condition) and d∆ dα (α ) < 0 (convergent stability condition), according to the adaptive dynamics framework [START_REF] Geritz | Evolutionarily singular strategies and the adaptative growth and branching of the evolutionary tree[END_REF]. We therefore investigate the singularity condition that provides an equation α has to satisfy. By writing the selection gradient under a condensed form and defining η :

= b I γ , φ := b D ε , ψ := b C σ , we find that ∆ (α ) = p α + θφ N (α ) -ψ (1 -α ) ψ + (η + θφα ) N (α ) = 0, ⇐⇒ (p -(1 + p) α ) ψ + (pη + (p + α ) θφ) N (α ) = 0. (S20)
As noticed in the previous subsection, this equation has a solution only if p is small enough (for the first term to be negative) and θ is small enough (for the second term not to be too positive). Hereafter, and because we are seeking for an intermediate evolutionary attracting virulence, we assume that p and θ are small enough such that this equation has a solution in ]0, 1[. The complexity of N (α ) prevents us from having a closed-form expression of α . However, α can be bounded by an underestimate α -on the one hand, and an overestimate α + on the other hand.

First, notice that the left hand side of equation (S20) has the following upper bound, replacing N (α ) by its maximum S 0 = λ µ , (p -(1 + p) α ) ψ + (pη + (p + α ) θφ) S 0 , which cancels out only if α is replaced by a greater value we denote by α + (since the first term is a decreasing function of the CFR). This leads to

α + = ((η + θφ) S 0 + ψ) p (1 + p) ψ -θφS 0 . ( S21
)
Notice that α + > 0 requires that ψ (1 + p) > θφS 0 that is θ small. Second and analogously, the left hand side of equation (S20) has the following lower bound, replacing N (α ) by its minimum

N (1) = 1 η+θφ , (p -(1 + p) α ) ψ + pη + (p + α ) θφ η + θφ ,
which cancels out only if α is replaced by a smaller value we denote α -(since the first term is a decreasing function of the CFR). This leads to

α -= (1 + ψ) p (1 + p) ψ -θφ η+θφ . (S22)
Notice that α -< 1 requires that p < ψ -θφ η+θφ . Since b C σ = 10 -2 , the approximations 1 + ψ ≈ 1 and 1 + p ≈ 1 holds, which makes this expression even simpler:

α -≈ p ψ -θφ η+θφ ,
which in turn is positive only if θ < ψη (1+ψ)φ ≈ ψη φ . Therefore, the two conditions on p and θ related to the cancellation of ∆ are retrieved.

As one can see by comparing Figure S1 with Figure 3, α + but moreover α -are accurate estimates of α .

Figure S1: Boundaries of the evolutionary attracting values of CFR. Underestimate (α -, left) and the overestimate (α + , right) values of the evolutionary attractor α as a function of unsafe burial ratio θ and trade-off exponent p. The solid, the dashed and the dotted lines correspond to α = 1, 0.7 and 0.3 respectively. Other parameter values are default (Table 1).

E Sensitivity analysis

For the sake of both generality and graphical readability, we reduced the dimensionality of the parameter space by defining:

η := b I γ , δ := b D ε • γ b I , κ := b C σ • γ b I , ( S23
)
where η is the average number of infectious contacts between one susceptible and one symptomatic individual over the symptomatic individual's symptomatic period, δ is the ratio between the equivalent of η for dead bodies and η itself, and κ is the ratio between the equivalent of η for convalescent individuals and η itself. Quantity η is equal to the symptomatic individual relative (that is normalised by S 0 ) contribution to the basic reproduction number, and is used for defining both δ and κ. Therefore, η is a primary scaling factor that can be eliminated through any estimated value of the R 0 . δ and κ are secondary scaling factors that can be studied independently. First, the basic reproduction number can be rewritten using definitions (S23),

R 0 = ((1 + αθδ) S 0 + (1 -α) κ) ηα p .
Therefore, any given a set of estimated epidemiological data ( α, θ, S 0 , R 0 ) can be used to scale η, while keeping undetermined the trade-off exponent p, that is to say

η = R 0 α -p 1 + α θδ S 0 + (1 -α) κ . (S24)
Rewriting the selection gradient at α = 1 with definitions (S23),

∆ (1) ≈ p + θδ 1 + θδ -κη,
and imputing estimated data with equation (S24),

∆ (1) ≈ p + θδ 1 + θδ - R 0 α -p κ 1 + α θδ S 0 + (1 -α) κ , we find that ∆ (1) < 0 ⇐⇒ ((1 + θδ) p + θδ) 1 + α θδ S 0 + (1 -α) κ < (1 + θδ) R 0 α -p κ, ⇐⇒ ((1 + θδ) p + θδ) 1 + α θδ S 0 < (1 + θδ) R 0 α -p -((1 + θδ) p + θδ) (1 -α) κ.
Further investigation requires to study the inequality

R 0 ≷ p + θδ 1 + θδ (1 -α) α p .
Elementary calculus then shows that for any couple (p, α) ∈ R + × [0, 1], we have the following upper bound

p + θδ 1 + θδ (1 -α) α p ≤ (p + 1) (1 -α) α p ≤ 1,
By definition, epidemiological data originate from outbreaks for which R 0 > 1, therefore we get the inequality used to plot our figure

∆ (1) < 0 ⇐⇒ κ S 0 > ((1 + θδ) p + θδ) 1 + α θδ (1 + θδ) R 0 α -p -((1 + θδ) p + θδ) (1 -α) , ( S25

F Application of the Price equation

Introducing a diversity of n ∈ N non-coinfecting strains of EBOV, system (S1) becomes a set of 4n + 1 ordinary differential equations, where for all i ∈ {1, . . . , n},

                 dS dt = λ -µS -n i=1 (β I,i I i + β D,i D i + β C,i C i ) S, dE i dt = (β I,i I i + β D,i D i + β C,i C i ) S -ω i E i , dI i dt = ω i E i -γ i I i , dD i dt = α i θγ i I i -ε i D i , dC i dt = (1 -α i ) γ i I i -σ i D i . (S26)
The total population size of each class, denoted by H

• := n i=1 H i , therefore satisfies                    dS dt = λ -µS -β I I I • + β D D D • + β C C C • S, dE• dt = β I I I • + β D D D • + β C C C • S -ω E E • , dI• dt = ω E E • -γ I I • , dD• dt = θαγ I I • -ε D D • , dC• dt = γ I I • -αγ I I • -σ D D • . (S27)
By definition, an average value of a trait x in a compartment H is given by x

H = n i=1 x i H i H• .
If we assume that the trait value of a strain is constant and neglect mutational variance (i.e. dx i dt = 0), the dynamics of any trait x in the I compartment are thus given by dx

I dt = n i=1 1 I • dI i dt - I i I 2 • dI • dt x i , = n i=1 (ω i E i -γ i I i ) 1 I • -ω E E • -γ I I • I i I 2 • x i , = n i=1 ω i E i E • -ω I i I • x i E • I • - n i=1 γ I i I • -γ I i I • x i , = E • I • n i=1 ω i E i E • -ω E i E • + ω E i E • -ω I i I • x i - n i=1 (γ -γ) x i I i I • , dx I dt = cov E (x, ω) + x E -x I ω E E • I • -cov I (x, γ) , ( S28
)
where cov indicates a genetic covariance between two traits, x H is the average value of trait x in host compartment X and xy H is the average value of the product xy in the same compartment. This illustrates that it might be difficult to disentangle a trait of interest x with the duration of the latent period (1/ω) if this latter trait varies for different virus genotypes.

Similarly, in the E compartment we have

dx E dt = n i=1 1 E • dE i dt - E i E 2 • dE • dt x i , = n i=1 β I i I i + β D i D i + β C i C i S -ω i E i 1 E • - β I I I • + β D D D • + β C C C • S -ωE • E i E 2 • x i , = S E • n i=1 β I i I i + β D i D i + β C i C i -β I I I • + β D D D • + β C C C • E i E • x i - n i=1 (ω i -ω) x i E i E • , = S E •   H∈{I,D,C} H • n i=1 β H i H i H • -β H H H i H • + β H H H i H • -β H H E i E • x i   -cov E (x, ω) , dx E dt = S E •   H∈{I,D,C} cov H (x, β H ) + x H -x E β H H H •   -cov E (x, ω) . ( S29
)
If we now focus on the dead hosts, we have

dx D dt = n i=1 1 D • dD i dt - D i D 2 • dD • dt x i , = n i=1 (α i θγ i I i -ε i D i ) 1 D • -θαγ I I • -εD • D i D 2 • x i , = θ I • D • n i=1 α i γ i I i I • -αγ I D i D • x i - n i=1 (ε i -ε) x i D i D • , = θ I • D • n i=1 α i γ i I i I • -αγ I I i I • + αγ I I i I • -αγ I D i D • x i - n i=1 (ε i -ε) x i D i D • , dx D dt = cov I (x, αγ) + x I -x D αγ I θ I • D • -cov D (x, ε) . (S30)
Finally, in the convalescent hosts, we have

dx C dt = n i=1 dC i dt 1 C • - dC • dt C i C 2 • x i , = n i=1 ((1 -α i ) γ i I i -σ i C i ) 1 C • -γ I -αγ I I • -σ C C • D i D 2 • x i , = I • C • n i=1 γ i I i I • -α i γ i I i I • -γ I -αγ I C i C • x i - n i=1 σ i -σ C x i C i C • , = I • C • n i=1 I i I • γ i -γ I -α i γ i + αγ I + γ I I i I • -γ I C i C • + αγ I C i C • -αγ I I i I • x i - n i=1 σ i -σ C x i C i C • , dx C dt = I • C • cov I (x, γ) -cov I (x, αγ) + x I -x C γ I -x I -x C αγ I -cov C (x, σ) . (S31)

G Numerical simulations

We explored 8 scenarios. For each, we assume that we have n = 100 EBOV strains. The standing genetic variation for the CFR α i (i ∈ {1, . . . n}) is drawn from a Gaussian distribution with mean α = 0.7 and standard deviation ς = 0.1 We only explored positive correlations between CFR and transmissions rates, consistently with our trade-off hypothesis. We also did not investigate correlations between CFR and post mortem elimination rate because we assume that the period over which an unsafe buried body is still a suitable environment for virion survival is independent from the initial number of virions. Finally, we ignored convalescentrelated variables received since this component of EBOV transmission is much smaller than the others two.

The description of the 9 scenarios is as follows:

1. No genetic correlation between the CFR α and other traits ("all constant" panel, also shown in the main text).

2. Addition of a negative correlation between α and the rate of end of latency period ω.

3. Addition of a positive correlation between α and the transmission rates β H ("+bH" panel) that will be kept for the next six scenarios.

4. Addition of a positive correlation between α and the inverse of the latency period ω ("+bH+O" panel).

5. Reversing the correlation between α and ω ("+bH-O" panel).

6. positive correlations between α and β H , α and ω and α and the inverse of the symptomatic period γ ("+bH+O+G" panel) 7. positive correlations between α and β H and between α and γ, negative correlation between α and ω ("+bH-O+G" panel) 8. positive correlations between α and β H and between α and ω, negative correlation between α and γ ("+bH+O-G" panel) 9. positive correlation between α and β H , negative correlations between α and ω and between α and γ ("+bH-O-G" panel)

Positively and negatively correlated traits were drawn according to the formulas x i = α i α + (1 -) ξ i x and x i = -(max (α) -α) α i α + max (α) -min (α) α-min(α) + (1 -) ξ i x respectively, where x the estimated value of x ≡ β H , ω, γ according to Table 1, The CFR averaged over the exposed individuals (α E ) is depicted in cyan, over the symptomatic individuals (α I ) in pink and over the infectious dead bodies (α D ) in brown.

H Virulence and transmission routes

We have showed that EBOV' reproduction numbers can be split into three additive components that correspond to each of the three transmission routes namely symptomatic (through regular contact with symptomatic individuals), post mortem (through contact with unsafe buried infectious dead bodies) and sexual (through sexual contact with convalescent individuals). According to our trade-off assumption, the intensity of each of these components is modulated by virulence: both symptomatic and post mortem components always increase with virulence while sexual component is maximum for an intermediate virulence level (unless there is no trade-off, in which case the symptomatic component is constant and the sexual component decreases with virulence), as in Figure S3. This come from the fact that virulence increases all transmission rates and infectious bodies inflows, while decreasing the convalescent individuals inflow. Virulence then also acts as an investment cursor between the exclusive post mortem and sexual transmission routes. It thus appears that the cost of virulence is strictly limited to loss in sexual transmission. Therefore, it is only if the sexual component is the dominant route of the virus' life cycle that this cost can really balances with the overall transmission Post mortem (left), symptomatic (middle), and sexual (right) transmission components as a function of virulence and trade-off shape (p = 0 in black, 0.05 in red, 0.5 in orange, 1 in blue and 2 in green).

Figure 2 :

 2 Figure 2: Basic reproduction number as a function case fatality ratio (α), unsafe burial ratio and trade-off shape.Colors indicate five trade-off scenarios: absence (p = 0, in grey), very weak (p = 0.1, in red), concave (p = 0.5, in green), linear (p = 1, in dark blue), and convex (p = 2, in light blue). The width of the coloured regions corresponds to variations in the unsafe burial ratio from completely unsafe (θ = 1, dashed upper bound) to completely safe (θ = 0, solid lower bound). The intersection between the horizontal line and the colored areas indicates the range of α min for each scenario. The dotted gridlines show the α and R 0 estimates from the literature. Other parameter values are in Table1. See Supplementary Material D.4 for more details.

Figure 3 :

 3 Figure 3: Evolutionary stable virulence (α ) as a function of unsafe burial ratio (θ) and tradeoff exponent (p). The solid, dashed and dotted lines correspond to α = 1, 0.7 and 0.3 respectively. Parameter vlaues are shown in Table1. See Supplementary Material D.4 for more details.

  )with values α = 0.7, θ = 0.25, S 0 = 4.44444 • 10 6 and R 0 ∈ [1.26, 2.53], according to reference (World Health Organization Ebola Response Team, 2014;[START_REF] Nyenswah | Ebola and its control in Liberia, 2014-2015[END_REF][START_REF] Abbate | Potential Impact of Sexual Transmission on Ebola Virus Epidemiology: Sierra Leone as a Case Study[END_REF] respectively.

  = 0.5 denotes the strength of the correlation and ξ i a Gaussian random variable with mean 1 and standard deviation ς = 0.1. Initial conditions are given by S (0) = λ µ ≈ 4.4 • 10 6 ind and H i (0) = 1 ind for all i ∈ {1, . . . , n} and all H ≡ E, I, D, C. Results for 8 of the scenarios are shown in FigureS2(scenario 2 is only shown in the main text for space constraint reasons).

Figure S2 :

 S2 Figure S2: Short-term evolution of CFR with standing genetic variation.The CFR averaged over the exposed individuals (α E ) is depicted in cyan, over the symptomatic individuals (α I ) in pink and over the infectious dead bodies (α D ) in brown.

Figure S3 :

 S3 Figure S3: Relative variation of transmission route component intensity as a function of virulence and trade-off shape.Post mortem (left), symptomatic (middle), and sexual (right) transmission components as a function of virulence and trade-off shape (p = 0 in black, 0.05 in red, 0.5 in orange, 1 in blue and 2 in green).

Table 1 : Parameter list, description and default values.

 1 See the main text for further details about the calibration of the transmission constants. Note that the sexual transmission constant is higher because it involves frequency-dependent transmission. ind stands for individuals. When used in the main text or in the appendix, these estimated values are denoted by a hat. Epidemiology of Ebola Virus in humans. A) Epidemiological life cycle of Ebola Virus in humans and B) Population dynamics for default parameters. S, E, I, C, R and D are host densities that correspond to the following compartments: susceptible, exposed (infected but not yet infectious), symptomatic infected, convalescent, recovered (immunised) and dead bodies respectively. N is the total living population size. Lower-case letters are rate and flow parameters, the description of which is given in Table1.

	Notation Description	Default value	Reference or inference
	λ	Host inflow	2 • 10 2 ind.day -1	such that λ/µ ≈ 4 • 10 6 , (Central Intel-
				ligence Agency, 2016)
	µ	Host baseline mortality	4.5 • 10 -5 day -1	(Central Intelligence Agency, 2016)
		rate
	b I	Regular contact transmis-sion (from infectious hosts)	10 -7 ind -1 .day -1
		factor
				with the constrain R 0,I ≈ 10 R 0,D ,
				R 0,I ≈ 10 2 R 0,C and R 0 ≈ 1.8
	b C	Sexual transmission (from convalescent hosts) factor	10 -4 day -1	(World Health Organization Ebola Re-sponse Team, 2014; Weitz and Dushoff,
				2015; Abbate et al., 2016)
	b D	Post mortem transmission (from dead hosts) factor	10 -8 ind -1 .day -1
	ω	Inverse of latency period	10 -1 day -1	(World Health Organization Ebola Re-
				sponse Team, 2014)
	α	Case fatality ratio	0.7	(World Health Organization Ebola Re-
				sponse Team, 2014)
	γ	Inverse of symptomatic in-	0.25 day -1	(Abbate et al., 2016)
		fectious period
	θ	Unsafe burial proportion	0.25	(Nyenswah et al., 2016)
	ε	Inverse of post mortem in-	10 -1 day -1	(Prescott et al., 2015)
		fectious period
	σ	Elimination rate of conva-	10 -2 day -1	(Uyeki et al., 2016; Abbate et al., 2016)
		lescent hosts

Acknowledgements

The authors thank Virginie Ravigné, François Blanquart and Hans Metz for their helpful comments and constructive discussions.

Funding Statement

M. T. Sofonea and S. Alizon acknowledge support from the CNRS, the UM and the IRD.

Competing Interests

We have no competing interests.

Authors' Contributions

H Virulence and transmission routes 17

A Equation systems

A.1 Epidemiology

The epidemiological dynamics are governed by the following set of ODEs:

where N := S + E + I + C + R is the total living population size, which varies with time. Notice that since life expectancy is several orders of magnitudes greater than the latency, the symptomatic and the convalescent periods, mortality rate µ can be neglected when summed with ω, γ or σ.

A.2 Price equation

The dynamics of the populations of interest are described by 4n + 1 ODEs, for all i ∈ {1, . . . , n}:

The total density of each compartment is denoted by a bullet index (•) and its dynamics satisfy

where the bars indicate average values and the superscripts indicate the compartment in which the trait is averaged. We can already notice that the CFR and the rate at which the infectious period ends are difficult to disentangle in this system because we have second order terms (i.e. an average of the product α i γ i ).

B Stationary dynamics B.1 Endemic equilibrium

At equilibrium, all time derivatives of system (S1) cancel out. If we denote by H the corresponding value of density H at this equilibrium, we get S = S 0 = λ µ and E = I = D = C = R = 0 for the disease free equilibrium (DFE).

The endemic equilibrium (EE), on the other hand, is found by assuming non zero values for all H. We thus first get that,

(S4)

B.2 Stationary total population size approximation for α = 1

In this subsection, we assume that α < 1 (the case where α = 1 is treated in the next subsection).

Given that life expectancy is several orders of magnitude greater than the convalescent period, i.e.

Furthermore, since life expectancy is also several orders of magnitude greater than the latency and symptomatic period, i.e. The virulence of EBOV is usually high and its sexual transmission low compared to the two other transmission route [START_REF] Abbate | Potential Impact of Sexual Transmission on Ebola Virus Epidemiology: Sierra Leone as a Case Study[END_REF], which is why we can approximate S by its value by neglecting