key: cord-289325-jhokn5bu
authors: Lachiany, Menachem; Louzoun, Yoram
title: Effects of distribution of infection rate on epidemic models
date: 2016-08-11
journal: Phys Rev E
DOI: 10.1103/physreve.94.022409
sha: 
doc_id: 289325
cord_uid: jhokn5bu

A goal of many epidemic models is to compute the outcome of the epidemics from the observed infected early dynamics. However, often, the total number of infected individuals at the end of the epidemics is much lower than predicted from the early dynamics. This discrepancy is argued to result from human intervention or nonlinear dynamics not incorporated in standard models. We show that when variability in infection rates is included in standard susciptible-infected-susceptible ([Formula: see text]) and susceptible-infected-recovered ([Formula: see text]) models the total number of infected individuals in the late dynamics can be orders lower than predicted from the early dynamics. This discrepancy holds for [Formula: see text] and [Formula: see text] models, where the assumption that all individuals have the same sensitivity is eliminated. In contrast with network models, fixed partnerships are not assumed. We derive a moment closure scheme capturing the distribution of sensitivities. We find that the shape of the sensitivity distribution does not affect [Formula: see text] or the number of infected individuals in the early phases of the epidemics. However, a wide distribution of sensitivities reduces the total number of removed individuals in the [Formula: see text] model and the steady-state infected fraction in the [Formula: see text] model. The difference between the early and late dynamics implies that in order to extrapolate the expected effect of the epidemics from the initial phase of the epidemics, the rate of change in the average infectivity should be computed. These results are supported by a comparison of the theoretical model to the Ebola epidemics and by numerical simulation.

A goal of many epidemic models is to compute the outcome of the epidemics from the observed infected early dynamics. However, often, the total number of infected individuals at the end of the epidemics is much lower than predicted from the early dynamics. This discrepancy is argued to result from human intervention or nonlinear dynamics not incorporated in standard models. We show that when variability in infection rates is included in standard susciptible-infected-susceptible (SIS) and susceptible-infected-recovered (SIR) models the total number of infected individuals in the late dynamics can be orders lower than predicted from the early dynamics. This discrepancy holds for SIS and SIR models, where the assumption that all individuals have the same sensitivity is eliminated. In contrast with network models, fixed partnerships are not assumed. We derive a moment closure scheme capturing the distribution of sensitivities. We find that the shape of the sensitivity distribution does not affect R 0 or the number of infected individuals in the early phases of the epidemics. However, a wide distribution of sensitivities reduces the total number of removed individuals in the SIR model and the steady-state infected fraction in the SIS model. The difference between the early and late dynamics implies that in order to extrapolate the expected effect of the epidemics from the initial phase of the epidemics, the rate of change in the average infectivity should be computed. These results are supported by a comparison of the theoretical model to the Ebola epidemics and by numerical simulation. DOI: 10.1103/PhysRevE.94.022409

An important element in theoretical epidemiology is the epidemic threshold, which specifies the condition for an epidemic to grow. In mean-field epidemiological models, the concept of the basic reproductive number, Ro, has been systematically employed as a predictor for epidemic spread and as an analytical tool to study the threshold conditions [1] [2] [3] [4] . An advantage of Ro is that in many models it determines both the threshold for the emergence of an epidemic and the expected final outcome of an outbreak [5, 6] . It has thus been widely used to gauge the degree of threat that a specific infectious agent will pose as an outbreak progresses [7, 8] .

However, over and over again differences have been observed between the predicted and observed sizes of epidemics, whereas in most cases the observed epidemic is much smaller than predicted [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] . Indeed, recent studies raise doubt about the validity of forecasting the outcome of epidemics using the Ro-based estimate [22, 23] .

Specifically, given an estimate of Ro, one can project in the ordinary differential equation (ODE) -based susceptibleinfected-recovered (SIR) and susciptible-infected-susceptible (SIS) models the future course of the epidemic. However, recent results have shown that these estimates are much larger than the true extent of the disease outcome [24] . This overestimate of late dynamics when using early dynamics has been raised by many authors in general cases, as well as specifically for the Ebola virus [25, 26] . We build upon these known observations to show that even in fully mixed models, * louzouy@math.biu.ac.il the early dynamics cannot be used to estimate the outcome of the late dynamics, and we propose an alternative approach.

This overestimate can be explained by a slow decrease in Ro. Ro can be reduced by human intervention, such as the removal of sick individuals from the society or removing carriers of the disease [27, 28] , or a limitation of movement either for the entire population or for people expressing clinical signs [29, 30] , which reduces the effective reproductive number [31] . Ro can also be reduced by passive vaccination of the population [32] or aggressive vaccination of the population [33, 34] .

All models mentioned above deal with external elements reducing Ro, and thus reducing the number of patients in steady state or the total number of removed individuals in SIR models. We argue here that the measure of Ro may not be indicative of future forecasts of the number of patients, even if no external factors or vaccinations are involved.

Consequently, additional information is required to estimate the number of people that will be affected by the disease in steady state, or when the epidemic is over. We propose here that the estimate of the change in the relative infection rate defined as I /I can be a good way to estimate the steady-state number of infected individuals in SIS models, and the total number of removed individuals in SIR models.

We introduce an epidemic model in which each of the susceptible individuals has a different probability to get infected and the same recovery probability. This model differs from uniform models, but also from network models. In the standard SIR and SIS models, the probability of infection is constant for everyone [1] [2] [3] [4] . However, in reality, different people have different tendencies to be clinically sick and infectious (e.g., elderly people, children, or immune and deficient individuals [35] [36] [37] [38] [39] [40] ).

In network SIR models, the infectivity of each node is a function of its degree (and thus it varies among nodes). However, network models assume a constant interaction pattern that may be realistic in sexually transmitted diseases, but it is not realistic for most noncontact (e.g., airborne or vehicleled) transmission. Moreover, in undirected networks, the probability to infect and to become infected are symmetrical [2] [3] [4] [41] [42] [43] , which is again mainly appropriate for sexually transmitted diseases, but not for airborne transmission.

We present evidence here that in a wide class of models, the variability in the probability to get infected can break the link between the early phase of the epidemics and the predicted outcome. We investigate the dynamical processes driving this result, its validity, and its consequences. We then compare the conclusions of this model to observed Ebola outbreaks. Such results have been presented in network models. However, in such models the connectivity is fixed over time [44] [45] [46] [47] . In the current analysis, we show that inhomogeneity has a crucial effect even in fully mixed systems.

The models used in this study are based on the SIS and SIR models. Each susceptible individual has a different probability to get infected, but all infected individuals have the same probability to infect other susceptible individuals. Specifically, individuals exist in two discrete states-"healthy" or "infected"-in the SIS model, and three discrete states-"healthy," "infected," or "recovered"-in the SIR model. At each time step, each susceptible (healthy) individual i is infected with rate β i . At the same time, infected individuals are cured and become susceptible again with rate γ in the SIS model, or recovered in the SIR model. The value of β i is only a function of the person getting infected and not the person infecting, and it does not vary over time for a given individual. Given a variable probability to get infected in the population, we define the number of susceptible individuals with a β i value between β and β + dβ as S(β)dβ and the number of infected individuals as I (β)dβ. We further define N(β) as N (β) = S(β) + I (β). Note that N (β) does not change over time and is thus equal to its initial condition. We have studied multiple distributions for this initial condition, as shall be further explained. S is the sum over all the susceptible individuals with different probability to get infected, I is the sum over all the infected individuals with different probability to get infected, and N is the sum over all the susceptible individuals and infected individuals with different probability to get infected. Formally,

β is the disease transmission rate (i.e., the probability that a person would be infected), and γ is the disease recovery rate. Note that here we use an integral approximation of the discrete sum. We will further show that this approximation is consistent with numerical simulations. The equations for the SIS model used here are

and N = I + S.

In the SIR model, we added the recovered individuals class. R(β) is defined as a recovered individual with a probability β to get infected. Formally,

None of the recovered individuals can be infected again in the SIR model. In this model, N (β) as N (β) = S(β) + I (β) + R(β). The equations of the SIR model are

where I is defined in Eq. (1), S is defined in Eq. (2), and R is defined above in Eq. (5) , when N = I + S + R.

In the different models studied here, we implicitly consider four different distributions for the probability of susceptible individuals to get infected:

(i) Constant infection rate for all susceptible individuals. This is equivalent to the mean-field SIS or SIR models. We used β = 2 × 10 −5 , γ = 1, and a population size of N = 100 000. Thus Ro = βN/γ = 2. All other models only differ in the distribution of β.

(ii) Uniform probability distribution of infection rates within a range 2 × 10 −11 < β < 4 × 10 −6 .

(iii) Normal probability distribution with mean μ and variance σ 2 , The probability of having infection rate β is

. We used μ = 2 × 10 −5 , and σ = 2μ.

(iv) Scale-free distribution with the probability of having infection rate β is P (β) = β −α , where α = 1.5, and β is limited to the range (2 × 10 −8 ,2 × 10 −4 ).

Monte Carlo simulations of the systems studied have been performed with a population size of N = 100 000. We assign a different infection rate to each individual using the four distributions above. The population is initiated with a small number of infected individuals, I 0 = 10. These 10 initial infected individuals are random individuals (i.e., they have β values randomly chosen from the population). Note that β is very low. However, β * N is high enough, allowing for the epidemics to spread, since N = 100 000. We assume full mixing and a mass action formalism, allowing for a very small number of infected individuals to spread the pathogen. Formally, we assign each susceptible individual at each iteration a probability β i t to be infected. Similarly, each infected individual is assigned in each iteration a probability of γ t = t to be removed and become susceptible again. The simulation updating is synchronous. The dynamics are simulated for different parameter values.

The ODEs were solved numerically using the MATLAB fourth-order Runge-Kutta method [48] , as applied in the MATLAB ode45 function assuming nonstiff equations [49] .

We study the behavior of an epidemic outbreak assuming that each susceptible individual has a different probability of getting infected. However, once it is infected, its contribution to the total force of infection is constant. In other words, P [S(β) + I → I (β) + I ] = β. This represents, for example, people with different levels of susceptibility to a given disease. Within this model, we study two possible models of recovery. Either the recovered hosts are immunized and then an SIR model is used, or the recovered hosts become susceptible again and then their β value does not change over time, and an SIS model is used. To study the effect of the distribution of β (the probability to be infected) on both initial and late dynamics, we assume one of four possible β value distributions: constant, uniform, Gaussian, or scale-free. In all cases, we maintain the expected value of β equal among models.

As will be further shown, the initial dynamics are only affected by the first moment of the distribution (the expected values of β), while the total number of infected individuals during the outbreak in the SIR model or the steady-state infected fraction in the SIS model can be strongly affected by the following moments. Thus, in some distributions, it is impossible to predict the "outcome" of the epidemics from the observed initial dynamics and the resulting estimate of Ro.

To examine the behavior of the infected class as a function of time, we developed a moment closure scheme, and we use the following notations:

For the SIR model, in addition to the above definitions, we define

We then use the following notation:

as a product of Eq. (8) with Eq. (1),

as a product of Eq. (9) with Eq. (2),

as a product of Eq. (10) with Eq. (5) , and

as a product of Eq. (7) with Eq. (3). For the SIS model in Eq. (4), the number of infected individuals can be estimated via (see Appendix A)

where I 1 is the first order of Eq. (11) and E N (β) is the first order of Eq. (7). An iterative equation can be developed to estimate the higher orders:

This would obviously lead to an infinite number of coupled equations. However, the number of equations can be limited using a simple moment closure method. The highest order of the set of ODEs is set to be zero, and the order below it is set to be a constant. Note that more advanced schemes could be used [50] [51] [52] . However, this simple scheme agrees well with simulation and is enough for the current analysis. For the SIR model in Eq. (6), the number of infected individuals can be estimated via

where I 1 is the first order of Eq. (11) , E N (β) is the first order of Eq. (7), and R 1 is the first order of Eq. (13). We developed a similar scheme for the SIR model, with Eq. (6). The iterative equation obtained for the number of susceptibles is

where S n is defined in Eq. (12) . The closure scheme above is also applied here. For example, if the stopping order is S 5 , we set S 5 = 0, and the previous order is set to

The moment closure scheme is consistent with the dynamics of the infected population simulated using a Monte Carlo simulation with the appropriate parallel distribution (Figs. 1-3) .

Increasing the number of moments for the solution of Eqs. (16) and (18) reduces the difference between Monte Carlo simulations and the moment closure results (data not shown). For example, in the SIS model, for the three cases of distribution of infection rates, an order of 20 coincides with the results obtained from the simulation. In the SIR model, only an eighth order was needed to reach a very high prediction. 20) for the SIR model. The results of a Monte Carlo simulation, the full moment closure scheme, and the approximation of the initial dynamics were compared, and a good fit was obtained. For all four distribution studied, N = 100 000, γ = 1, and E N (β) = 2 × 10 −5 . The rising line is the initial time approximation, and the two other, highly similar, lines are the full moment closure solution and the simulation results.

For the early dynamics of the SIS model in Eq. (15), we can approximate the dynamics by neglecting elements of the order of −I I 1 compared with other elements, since β 1 and (20) and the recovered individuals number as

Equation (19) for the SIS model and Eq. (20) for the SIR model are consistent with the simulations (Fig. 1) . Furthermore, as can be clearly seen, in both SIS and SIR, all distributions tested (uniform, Gaussian, and scale-free) have the same dynamics as the constant infectivity.

During the early dynamics, only the first moment of the sensitivity distribution affects the dynamics. Thus, the higher moments of the distribution cannot be estimated from the observed value of R 0 . However, these moments may affect the dynamics later in the epidemic. We thus examined if the term neglected in the previous subsection, (−I I 1 ), has a differential effect as a function of the higher moments of the infection rate distribution. The general solution of the infected class for Eq. (15) can be estimated to be (see Appendix A)

where

The ratio between the solution of Eq. (22) and the solution without the neglected term in Eq. (19) is

For t = 0, r = 1. For t ≈ 0, one can estimate

The same analysis can be performed for the SIR model. We examined again the effect of the neglected terms on the dynamics, with the general solution of the infected class in Eq. (17) approximated by

where

. The ratio between the solution without the neglected terms [Eq. (20) ] and the solution of Eq. (25), marked by r(t), is

with a similar approximation. One can clearly see that the difference between the models is affected by the value of E I (β)(t). Thus if E I (β)(t) differs among models, so would the resulting dynamics.

As can be seen in Fig. 2 uniform distributions, the values are close to those obtained with constant infectivity. However, in the SF distribution, r(t) deviates from the values obtained in the constant infectivity case. The source of the difference is in the drastic change in E I (β) in the SF model over time. This affects the denominator of Eq. (24) for the SIS and SIR models.

This difference results from the effect of rare events in the SF distribution. While in the uniform and Gaussian distribution the value of E I (β) is close to E S (β), this is not the case in the SF model, and early in the dynamics E I (β) E S (β). Thus in Eq. (24), r(t) is expected to decrease much faster than in the constant β scenario, as is indeed the case. The source of the difference between the SF and all other distributions is the presence of the very high β values in the population (even if they are rare). People with such high values of β are almost automatically infected, and they increase E I (β) sharply. Note that they also slightly decrease E S (β). In all other distributions, the variance of β is limited per definition by the requirement that all individuals have positive β values. Thus, while the normal distribution cannot have a large variance, in the case of the SF distribution, β could get in principle values between (0,∞). This is not the case in reality, since the sample is finite. Still, the upper bound is a few orders above the average (Fig. 2) .

The same happens for the SIR model, with the important distinction that the individuals with high β values are rapidly removed from the population, leading to a decrease in E I (β) following the initial sharp rise. Thus, while initially all distribution show dynamics purely determined by R 0 , the dynamics evolve differently as a function of the distribution of β.

One is thus led to ask whether the difference in the distribution can be estimated from the early dynamics. The conclusion from the results above are that in order to estimate the future dynamics, it is not enough to know R 0 , but the change in E I (β) should be estimated. While this cannot be estimated directly, we can directly quantify the effect this has on S * E S (β) through the disease dynamics. Specifically, a decrease in E S (β) is expected to lead to a parallel decrease in I /I . In the bottom graph of Fig. 3 , we calculated the difference in the number of infected individuals in every time step, divided by the total number of infected individuals in the same time. The total decrease in the expected value of β in the S compartment and the increase in the expected value of β in the I compartments are clear in the SF distribution, while it is more limited in all other distributions. This can be clearly seen in the measures of E S (β), in E I (β), and in I /I . The mechanism driving the difference is straightforward. In normal or uniform distribution the variance in β is limited, thus even infecting the most susceptible individuals does not significantly affect the values of E I (β) and E S (β). In the SF model, the expected value of β is strongly affected by the tail of the distribution.

We have shown in the previous subsections that the distribution of the infection rate has a drastic effect on the dynamics for SF distributions. For other distributions, the infected class dynamics are similar to the constant infectivity model, or slightly lower. At the steady state of the system, the difference is further enlarged. In the SIS model [Eq. (4)], we can compute the equilibrium frequency of infected individuals (see Appendix B) to be the solution of the following implicit equation: 

where C n are the Taylor series coefficients,

Expanding the terms [β − E N (β)] n leads to

As expected from the intermediate period, in the constant, uniform, and Gaussian distributions, the same number of infected individuals is obtained in the SIS model (Fig. 4) . However, in agreement with the intermediate period, the total number of infected individuals in equilibrium is much lower for the SF distribution. A similar analysis can be performed for the SIR model with similar results (Fig. 4) . The results obtained in both models in simulations and in Eq. (30) are equivalent (Fig. 4) . moments of β are important, the first moment approximation (that is only affected by the first moment of β) obviously fails to properly reproduce the number of infected individuals.

The intuitive explanation for the effect of the distribution of β on the number of infected individuals in equilibrium, and the resulting reduction in the number of people affected by the epidemics compared with the case of constant β, is the removal of individuals with high infection probability from the susceptible pool. This effect is directly quantifiable through the higher moments of β in the population and the resulting change in the first moment of β in the susceptible and infected pool. Only in the distribution where the higher moments can be important will there be a difference between the models with constant infectivity and the models with a larger variability.

To confirm the effect of heterogeneity in β on the outcome of observed epidemics, and that the steady-state number of infected individuals is much smaller than expected from the constant transmission rate model, we compared between observed epidemics and the analytical results described above. We studied the spread of the Ebola virus in three African countries.

The Ebola virus is one of five viruses of the Ebolavirus genus [53] . Four of the five known Ebola viruses, including EBOV, cause a severe and often fatal hemorrhagic fever in humans and other mammals, known as Ebola virus disease (EVD). The Ebola virus has caused the majority of human deaths from EVD, and it is the cause of the 2013-2015 Ebola virus epidemic in West Africa, which resulted in at least 28 424 suspected cases and 11 311 confirmed deaths [54] . The natural reservoir of the Ebola virus is believed to be bats, and it is primarily transmitted between humans and from animals to humans through bodily fluids [55] .

We analyzed the Ebola virus daily infection rates collected by health authorities in three African countries: Guinea, Liberia, and Sierra Leone. We computed the theoretical parameters providing the best fit to the epidemics in each of the three distributions described above as well as the standard SIR models. There are now known cases of reexposure to Ebola. We thus fit the SIR and not the SIS model. For example, the free parameters are as follows: for the SF distribution N is number of the population, α,β min ,β max ,I 0 ,γ ; for the Gaussian distribution, N,μ,σ,γ,I 0 ; for the uniform distribution, N,β min ,β max ,γ ,I 0 ; and for the constant model, N,β,γ,I 0 . Using the observed I 0 leads to suboptimal results for all models. This is probably the case since in the early dynamics, stochastic fluctuations can affect the results. We thus estimate I 0 from the dynamics later when enough cases are available.

In Fig. 5 , the top graph represents the observed infected class in Sierra Leone and the bottom graph represents the observed infected class in Guinea. We added the real times as labels. Both graphs are fitted to the models with the distribution above. While in Sierra Leone there is practically no difference between the different models (as can be observed in the F scores in Fig. 6 ), in the Guinea case there is a large difference. The scale-free and normal distributions provide a better fit than the constant case and uniform cases. . This rapid decrease even in the early stage of the disease suggests that highly infectious individuals are rapidly removed from the population. In the bottom graph, an F test was calculated between the best fit in the constant case and all the other distributions to incorporate the different number of free parameters. Three asterisks between the classic model and any of the distributions represent p < 0.001. In Guinea, all three distributions are much better than the constant case, but there is no significant difference between the quality of the nonconstant distributions. In Liberia, the SF and Gaussian distributions have a significantly better fit than the uniform distribution, and the same result is obtained for Sierra Leone.

The relative infection rate I /I can be used to detect significant deviations from the straightforward dynamics expected from these models quite early in the dynamics. The top of Fig. 6 represents the observed infected class in Sierra Leone for the relative infection rate I /I . We computed from the observed epidemics the number of new infection cases divided by the current number of infected individuals with a moving window of 3 days.

To determine the best-fitting model for multiple states in Africa, we calculated the sum of squared errors (SSE) of the optimal fit obtained for each distribution, where the error is the difference between the observed number of infected individuals and the solution of the theoretical model for the appropriate distribution. In the bottom graph, the SSE is plotted for three countries: Guinea, Liberia, and Sierra Leone. For all countries, the SSE of all distribution of transmission coefficient is smaller than the classic case.

An F test [56] [57] [58] was conducted to determine whether the reduced SSE is statistically significant. The F statistic is computed using one of two equations depending on the number of parameters in the models. If both models have the same number of parameters, the formula for the F statistic is F = σ SSE1 /σ SSE2 , where SSE1 is for the first model and SSE2 is for the second model. The p value of the results is computed using W -V degrees of freedom, where W is the number of data points and V is the number of parameters being estimated (one degree of freedom is lost per parameter estimated). The resulting F statistic can then be compared to an F distribution. If the models have different numbers of parameters, the formula becomes

where df is the number degree of freedom of every model. df1 is number of degree of model 1 and df2 is number of degree of model 2. In the bottom of Fig. 6 , the F test was calculated between the constant case and all the other distributions. Three asterisks between the classic model and any of the distributions represent a better fit of any of the distributions than the classic case. In Guinea, all three distributions are better than the constant case, but we cannot determine which of the three distributions is preferable. In Liberia, the SF and Gaussian distributions were shown to be preferable to the uniform distribution, and the same result was obtained for Sierra Leone. We also compared between the results of the const susceptible-exposed-infected-recovered (SEIR) model and other distributions of the SIR model by the F test, and the residuals of the SEIR model are practically identical to those of the SIR model, and in this case adding an extra parameter does not improve the fit to the real data.

We investigated the SIR model and the SIS model for the three cases of distribution of the transmission coefficient in the population and the classic case with a constant transmission coefficient. We found that in realistic cases, nonuniform models provide a better description of the observed epidemics than the model where there is a constant transmission coefficient for the entire population. In the early dynamics, only the first moment of β determines the dynamics, and there is no difference between the models as long as they have the same transmission rate. Later, the second moment of β starts affecting the dynamics for both the SIR and SIS models. During this period, the scale-free distribution behaves differently from other distributions and results in a lower number of infected individuals than expected. This reduction is the result of a difference between the first moment of β in the different compartments (S, I ). There is a sharp increase in the average value of β among infected compartments, accompanied by a smaller decrease in the susceptible compartment, compared with the initial value. Such an effect can only be observed in distributions where the second moment of β in the total population is large enough. The difference is then further enlarged in steady state where all distributions of the transmission coefficient in the population lead to a smaller number of infected individuals compared with the model with constant β, but the biggest difference is in the SF model, again explained by the high second moment of β in this distribution. We have then tested this conclusion by studying the spread of the Ebola virus in multiple African countries. An F test between the different distributions shows that they all produce a better fit than the constant β model. We attempted to detect early in the infection whether the total number of people affected by the epidemics will deviate significantly from the results expected from the early dynamics.

Classical models predict a large number of infected individuals in most epidemics with R 0 higher than 1. However, in reality, many epidemics end with a limited impact. The most classical example is perhaps the huge difference between the predictions and the observed amplitude of the JCD [9] [10] [11] [12] [13] [14] [15] [16] . Many models can explain this discrepancy, including, among others, nonlinear dynamics [59, 60] , delays [61, 62] , human intervention, passive vaccination [63, 64] , and small effective population [65] [66] [67] . We show here that even in the most standard SIR and SIS models, the initial dynamics cannot determine the total number of people affected by the epidemics. However, quite early in the dynamics, the relative infection rate I /I can be used to detect significant deviations from the straightforward dynamics expected from these models. The large difference is only expected if the distribution of β values is large. Such a large difference can be the result of intrinsic differences, but also the result of environmental differences, partial mixing, or subpopulation structure. Understanding this distribution in advance would improve our capacity to relate early and late dynamics. We now plan to develop methods to estimate this distribution from finer measures early in the disease dynamics.

We thank Miriam Beller for the language editing.

We solve Eq. (4) using iterative equations for the SIS model. The equation for the infected compartment from Eq. (4) can be integrated over all values of β as follows:

The term ∞ 0 βN(β)dβ is equal to E(β)N , and the term ∞ 0 βI (β)dβ can be defined as I 1 . Equation (A3) thus leads to dI dt

where it has a form similar to Eq. (15) . We can differentiate higher orders of I with respect to time to obtain equations similar to Eq. (A4). This can be written at the nth order of I (I n ) as the following general solution:

For the early dynamics for the SIS model, we dropped the term −I I 1 from Eq. (15) , as explained in the main text, and we obtained

In the second order for initial time, the term −I I 1 in Eq. (A4) is not neglected. We can write this term as −I 2 E I (β) using Eq. (11) . The new form of Eq. (A4) leads to

We denote −E N (β)N + γ = C 1 and −E I (β) = C 2 to obtain dI dt

The solution for that equation is

We get

One can change the constants to be E N (β)N − γ = ξ and E I (β) = ψ, and we obtain

The boundary condition is I (t = 0) = I o , leading to

The solution for the infected class in Eq. (A13) is

The ratio between the solution without the neglected terms Eq. (A8) and the current solution of Eq. (A16), marked by r(t), is 

For t ≈ 0, we get

and then we get 

where it has the same form as Eq. (17).

We drop the terms −I I 1 and −I R 1 for the early dynamics to obtain

and R is dR(β) dt = γ I (β).

We integrate again and obtain d dt 

which is

In the second order for initial time, the terms −I I 1 and −I R 1 in Eq. (C4) will not be neglected. The same development for the term I is identical to that in Appendix A.

We solve Eq. (6) with the closure scheme by integrating by β: If we differentiate S 1 , we get

and the equation will be

where we define S 2 = ∞ 0 β 2 S(β)dβ. In the same way, we can define in general

and the difference equation is

We use the moment closure method to derive the solution for these differential equations. The order where the differential equation stops will be zero, and the order before it will be const.

For example, if we want to stop at order S 5 , this order will be S 5 = 0, and order S 4 = const, where the const = S 4 (t = 0) = ∞ 0 β 4 S(β)dβ.

Infectious Diseases in Humans

Mathematical Tools for Understanding Infectious Disease Dynamics

Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation

Modeling Infectious Diseases in Humans and Animals

Generality of the final size formula for an epidemic of a newly invading infectious disease

Age-of-infection and the final size relation

The mathematics of infectious diseases

Pandemic potential of a strain of influenza A (H1N1): Early findings

Representations of mad cow disease

New phenotypes for new breeding goals in dairy cattle

Mad cow policy and management of grizzly bear incidents

No brainer-The USDA's regulatory response to the discovery of mad cow disease in the United States

Trust in food in the age of mad cow disease: A comparative study of consumers' evaluation of food safety in Belgium

Risk perception of the "mad cow disease" in France: Determinants and consequences

Villas-Boas, Consumer and market responses to mad cow disease

Tracking the human fallout from mad cow disease

Receptor recognition and cross-species infections of SARS coronavirus

Cross-reactive antibodies in convalescent SARS patients' sera against the emerging novel human coronavirus EMC (2012) by both immunofluorescent and neutralizing antibody tests

The genome sequence of the SARS-associated coronavirus

Clinical progression and viral load in a community outbreak of coronavirus-associated SARS pneumonia: A prospective study

Prion diseases: Update on mad cow disease, variant Creutzfeldt-Jakob disease, and the transmissible spongiform encephalopathies

Early estimation of the reproduction number in the presence of imported cases: Pandemic influenza H1N1-2009 in New Zealand

Epidemic models with uncertainty in the reproduction number

Rapid drop in the reproduction number during the Ebola outbreak in the Democratic Republic of Congo

Predicting the extinction of Ebola spreading in Liberia due to mitigation strategies

Spatiotemporal spread of the 2014 outbreak of Ebola virus disease in Liberia and the effectiveness of non-pharmaceutical interventions: A computational modeling analysis

Identifying transmission cycles at the humananimal interface: The role of animal reservoirs in maintaining gambiense human african trypanosomiasis

Transmission and control of African horse sickness in The Netherlands: A model analysis

Simulating school closure strategies to mitigate an influenza epidemic

World Health Organization Writing Group, Nonpharmaceutical interventions for pandemic influenza, national and community measures

Estimating the effective reproduction number for pandemic influenza from notification data made publicly available in real time: A multi-country analysis for influenza A/H1N1v

Vaccination and passive immunisation against Staphylococcus aureus

A vaccination model for a multi-city system

Modeling vaccination in a heterogeneous metapopulation system

Strategies for mitigating an influenza pandemic

Clinical recognition and diagnosis of Clostridium difficile infection

Ageing and infection

Case definitions of clinical malaria under different transmission conditions in Kilifi District

The management of community-acquired pneumonia in infants and children older than 3 months of age: Clinical practice guidelines by the Pediatric Infectious Diseases Society and the Infectious Diseases Society of America

Clinical practice guidelines for Clostridium difficile infection in adults: 2010 update by the Society for Healthcare Epidemiology of America (SHEA) and the Infectious Diseases Society of America (IDSA)

Infection dynamics on scale-free networks

Automata network SIR models for the spread of infectious diseases in populations of moving individuals

Epidemic outbreaks in complex heterogeneous networks

Slow epidemic extinction in populations with heterogeneous infection rates

Dynamics of person-to-person interactions from distributed RFID sensor networks

Griffiths Phases on Complex Networks

Epidemic spreading in metapopulation networks with heterogeneous infection rates

Numerical solutions of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes

Mastering MATLAB 5: A Comprehensive Tutorial and Reference

Moment closure and the stochastic logistic model

Novel moment closure approximations in stochastic epidemics

Moment-closure approximations for mass-action models

Proposal for a revised taxonomy of the family Filoviridae: Classification, names of taxa and viruses, and virus abbreviations

Ebola situation report

Killers in a Cell but on the Loose-Ebola and the Vast Viral Universe

The Analysis of Variance: Fixed, Random and Mixed Models

Mathematics for the Clinical Laboratory

Kenward-Roger approximate F test for fixed effects in mixed linear models

Mathematical Structures of Epidemic Systems

A stabilizability problem for a reaction-diffusion system modeling a class of spatially structured epidemic systems

Production and germination of conidia of Trichoderma stromaticum, a mycoparasite of Crinipellis perniciosa on cacao

Differential Equations and Applications in Ecology, Epidemics, and Population Problems

Effect of vaccination in environmentally induced diseases

Listeriosis: A model for the fine balance between immunity and morbidity

Temporal genetic samples indicate small effective population size of the endangered yellow-eyed penguin

Small effective population size and genetic homogeneity in the Val Borbera isolate

Effective population size is positively correlated with levels of adaptive divergence among annual sunflowers