Effects of distribution of infection rate on epidemic models.

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.

INTRODUCTION

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-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-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 susceptible-infected-recovered () and susciptible-infected-susceptible () 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, 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 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 can be a good way to estimate the steady-state number of infected individuals in models, and the total number of removed individuals in 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 and models, the probability of infection is constant for everyone [1-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-40]).

In network 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 vehicle-led) transmission. Moreover, in undirected networks, the probability to infect and to become infected are symmetrical [2-4,41-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-47]. In the current analysis, we show that inhomogeneity has a crucial effect even in fully mixed systems.

MODELS AND METHODS

model

The models used in this study are based on the and 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 model, and three discrete states—“healthy,” “infected,” or “recovered”—in the model. At each time step, each susceptible (healthy) individual is infected with rate . At the same time, infected individuals are cured and become susceptible again with rate in the model, or recovered in the model. The value of 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 value between and as and the number of infected individuals as . We further define as . Note that 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. is the sum over all the susceptible individuals with different probability to get infected, is the sum over all the infected individuals with different probability to get infected, and 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 model used here are and .

In the model, we added the recovered individuals class. is defined as a recovered individual with a probability to get infected. Formally, None of the recovered individuals can be infected again in the model. In this model, as . The equations of the model are where is defined in Eq. (1), is defined in Eq. (2), and is defined above in Eq. (5), when .

Distributions

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 or models. We used , , and a population size of . Thus . All other models only differ in the distribution of .

(ii) Uniform probability distribution of infection rates within a range .

(iii) Normal probability distribution with mean and variance , The probability of having infection rate is in the domain in (). We used , and .

(iv) Scale-free distribution with the probability of having infection rate is , where , and is limited to the range ().

Simulations

Monte Carlo simulations of the systems studied have been performed with a population size of . 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, . These 10 initial infected individuals are random individuals (i.e., they have values randomly chosen from the population). Note that is very low. However, is high enough, allowing for the epidemics to spread, since . 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 to be infected. Similarly, each infected individual is assigned in each iteration a probability of to be removed and become susceptible again. The simulation updating is synchronous. The dynamics are simulated for different parameter values.

Numerical solution of ODEs

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].

RESULTS

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, . 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 model is used, or the recovered hosts become susceptible again and then their value does not change over time, and an 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 model or the steady-state infected fraction in the 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.

Moment closure method

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 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 model in Eq. (4), the number of infected individuals can be estimated via (see Appendix A) where is the first order of Eq. (11) and 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-52]. However, this simple scheme agrees well with simulation and is enough for the current analysis.

For the model in Eq. (6), the number of infected individuals can be estimated via where is the first order of Eq. (11), is the first order of Eq. (7), and is the first order of Eq. (13). We developed a similar scheme for the model, with Eq. (6). The iterative equation obtained for the number of susceptibles is where is defined in Eq. (12).

The closure scheme above is also applied here. For example, if the stopping order is , we set , and the previous order is set to , where .

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).

FIG. 1.

Infected individuals as a function of time. Initial time dynamics based on Eq. (19) for the model and Eq. (20) for the 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, , , and . 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.

FIG. 3.

Plot of and as a function of time for the model and for the model for all distributions tested. While in the uniform and Gaussian distribution the value of is close to , in the SF model this is not the case, and early in the dynamics . The source of the difference between the SF and all other distributions is the presence of very high values in the population (even if they are rare). People with such high values of are almost automatically infected, and they increase sharply. The bottom plot represents the relative infection rate as a function of time. The total decrease in the expected value of in the compartment and the increase in the expected value of in the compartment induce a rapid decrease in this rate in the SF distribution.

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 model, for the three cases of distribution of infection rates, an order of 20 coincides with the results obtained from the simulation. In the model, only an eighth order was needed to reach a very high prediction.

Early dynamics

For the early dynamics of the model in Eq. (15), we can approximate the dynamics by neglecting elements of the order of compared with other elements, since and . The resulting solution of Eq. (15), when , is (see Appendix A) where is the initial number of infected individuals in the population. To check this approximation, we compared it to the results of a Monte Carlo simulation, and a good fit was obtained. For all four types of distributions, in the early time of the dynamics, Eq. (19) fits the full moment closure scheme as well as the values obtained from the Monte Carlo simulation.

For the model, we can use a similar approach in Eq. (17). We can approximate the dynamics by ignoring elements of the order of and compared with other elements, since and and . The resulting solution of Eq. (17), when , is (see Appendix C) and the recovered individuals number as Equation (19) for the model and Eq. (20) for the model are consistent with the simulations (Fig. 1). Furthermore, as can be clearly seen, in both and , all distributions tested (uniform, Gaussian, and scale-free) have the same dynamics as the constant infectivity.

Second-order analysis

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 . However, these moments may affect the dynamics later in the epidemic. We thus examined if the term neglected in the previous subsection, , 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 , . For , one can estimate

The same analysis can be performed for the 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 , is with a similar approximation. One can clearly see that the difference between the models is affected by the value of . Thus if differs among models, so would the resulting dynamics.

As can be seen in Fig. 2, the values obtained from the Monte Carlo simulations and the values estimated in Eq. (23) for the model and Eq. (26) for the model agree for all distributions. Note that in both the Gaussian and the uniform distributions, the values are close to those obtained with constant infectivity. However, in the SF distribution, deviates from the values obtained in the constant infectivity case. The source of the difference is in the drastic change in in the SF model over time. This affects the denominator of Eq. (24) for the and models.

FIG. 2.

[Eq. (23) for the model and Eq. (26) for the model]. The axis is time since the initial infections. These values agree for all distributions with the values obtained from Monte Carlo simulations. Note that in both the Gaussian and the uniform distributions, the values are close to those obtained with constant infectivity. However, in the SF distribution, deviates from the values obtained in the constant infectivity case. The source of the difference is in the drastic change in in the SF model over time. This affects the denominator of Eq. (24) for the and models. For all four cases, , , and .

This difference results from the effect of rare events in the SF distribution. While in the uniform and Gaussian distribution the value of is close to , this is not the case in the SF model, and early in the dynamics . Thus in Eq. (24), 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 sharply. Note that they also slightly decrease . 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 . 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 model, with the important distinction that the individuals with high values are rapidly removed from the population, leading to a decrease in following the initial sharp rise. Thus, while initially all distribution show dynamics purely determined by , 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 , but the change in should be estimated. While this cannot be estimated directly, we can directly quantify the effect this has on through the disease dynamics. Specifically, a decrease in is expected to lead to a parallel decrease in . 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 compartment and the increase in the expected value of in the 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 , in , and in . 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 and . In the SF model, the expected value of is strongly affected by the tail of the distribution.

Steady-state analysis

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 model [Eq. (4)], we can compute the equilibrium frequency of infected individuals (see Appendix B) to be the solution of the following implicit equation: One can develop the expansion of the integrand of Eq. (27) around the average infection rate, leading to an expression of the different distributions through the moments of : Eq. (27) then becomes where are the Taylor series coefficients, Expanding the terms 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 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 model with similar results (Fig. 4). The results obtained in both models in simulations and in Eq. (30) are equivalent (Fig. 4). For the uniform, constant, and normal distribution, where the higher moments of are of limited importance, the results of a first-order approximation of Eq. (30) are very similar to the full expansion and the simulation (about relative error). However, in the SF model where the higher 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.

FIG. 4.

Comparison between the analytic results and simulation results for the steady-state value in the four distributions studied (upper plot), and the total number of removed individuals in (lower plot). The SF model distribution has the lowest steady-state concentration of infected individuals in and the lowest number of removed individuals in . This results from 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.

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.

Application to observed epidemics

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 models. There are now known cases of reexposure to Ebola. We thus fit the and not the model. For example, the free parameters are as follows: for the SF distribution is number of the population, ; for the Gaussian distribution, ; for the uniform distribution, ; and for the constant model, . Using the observed 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 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 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.

FIG. 5.

Observed number of Ebola infections in Sierra Leone (upper plot) and Guinea (bottom plot) and a comparison of the dynamics with the best fit obtained using the different distributions.

FIG. 6.

Decrease of the relative infection rate as a function of time (top graph). 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 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 . 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 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 . 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 test [56-58] was conducted to determine whether the reduced SSE is statistically significant. The 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 statistic is , where is for the first model and is for the second model. The value of the results is computed using degrees of freedom, where is the number of data points and is the number of parameters being estimated (one degree of freedom is lost per parameter estimated). The resulting statistic can then be compared to an 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 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 model by the test, and the residuals of the SEIR model are practically identical to those of the model, and in this case adding an extra parameter does not improve the fit to the real data.

CONCLUSIONS

We investigated the model and the 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 and 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 (, ). 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 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 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-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-67]. We show here that even in the most standard and 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 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.