Time-varying and state-dependent recovery rates in epidemiological models.

Differential equation models of infectious disease have undergone many theoretical extensions that are invaluable for the evaluation of disease spread. For instance, while one traditionally uses a bilinear term to describe the incidence rate of infection, physically more realistic generalizations exist to account for effects such as the saturation of infection. However, such theoretical extensions of recovery rates in differential equation models have only started to be developed. This is despite the fact that a constant rate often does not provide a good description of the dynamics of recovery and that the recovery rate is arguably as important as the incidence rate in governing the dynamics of a system. We provide a first-principles derivation of state-dependent and time-varying recovery rates in differential equation models of infectious disease. Through this derivation, we demonstrate how to obtain time-varying and state-dependent recovery rates based on the family of Pearson distributions and a power-law distribution, respectively. For recovery rates based on the family of Pearson distributions, we show that uncertainty in skewness, in comparison to other statistical moments, is at least two times more impactful on the sensitivity of predicting an epidemic's peak. In addition, using recovery rates based on a power-law distribution, we provide a procedure to obtain state-dependent recovery rates. For such state-dependent rates, we derive a natural connection between recovery rate parameters with the mean and standard deviation of a power-law distribution, illustrating the impact that standard deviation has on the shape of an epidemic wave.

Introduction

Compartmental models of infectious-disease transmission have proven to be an invaluable tool for the prediction of disease progression and the evaluation of public health policies and interventions. In general, compartmental models describe how an infectious disease propagates throughout a population by characterizing the rates of transition of a populace from states of susceptible to, infected with, and finally recovered from disease. The rate of transition from the susceptible state to the infected state is often called the force of infection or incidence rate (Liu, Hethcote, & Levin, 1987), and the rate of transition from the state of infected with the disease to recovered from disease is called the recovery rate. While classically the incidence rate is taken as a bilinear term, there has been extensive theoretical work regarding alternative nonlinear incidence rates. These alternative nonlinear incidence rates are normally developed to account for behavioral characteristics, such as the crowding of infected individuals and the avoidance of exposure to infection (Alexander and Moghadas, 2005, Liu et al., 1987, Liu et al., 1986, Ruan and Wang, 2003, van den Driessche and Watmough, 2000) or to account for multi-stage infections (Krylova & Earn, 2013). In fact, formulations of nonlinear incidence rates can be justified through a first-principles derivation (Ponciano & Capistrán, 2011).

Like incidence rates, recovery rates have received significant attention in regards to their formulations in integral equations (Feng et al., 2007, Fowler and Hollingsworth, 2015, Hethcote and Tudor, 1980) and stochastic epidemic models (Ball, 1983, Ball, 1986, Ball, 1991, Britton, 2010, Clancy, 1999, Clancy, 2014). Recovery rates have also been generalized in ordinary differential equation (ODE) models, mainly through the method of stages, whereby the infected state is broken into multiple stages, where each stage has a constant recovery rate (Krylova and Earn, 2013, Lloyd, 2001a, Lloyd, 2001b). However, theoretical extensions of recovery rates beyond multiple constants in ODEs have only started to be developed, despite the fact that constant recovery rates based on the mean value of an exponentially distributed infectious period are epidemiologically unrealistic (Bailey et al., 1975, Gough, 1977, Keeling and Grenfell, 1999, Keeling and Rohani, 2007, Lloyd, 2001a, Lloyd, 2001b).

Here we derive time-varying and state-dependent recovery rates for ODE models of infectious diseases. To obtain such novel recovery rates, we use the connection between integral equation models, ODEs and survival functions. First we show how an integral equation model for an infectious disease (Hethcote & Tudor, 1980) is related to more typical ODEs. A key component of this relationship is the probability distribution of the time spent in the infectious class. Most traditional models assume this time is exponentially distributed, and so we review this distribution and its connection to stochastic processes. Next, we show how an alternative probability distribution of the time spent in the infectious class leads to different forms of recovery rates in ODE models.

We validate our work by modeling historical measles outbreaks in Iceland (Cliff, Haggett, Ord, & Versey, 1981). First, we develop a compartmental model with a time-varying recovery rate based on infectious periods that follow the family of Pearson distributions (Pearson, 1893, Pearson, 1895). The family of Pearson distributions span all possible values of skewness and kurtosis, so we evaluate the sensitivity of model predictions of the timing and magnitude of an epidemic peak with such time-varying recovery rates, relative to the uncertainty in the skewness and kurtosis derived from infectious period data.

Finally, we turn to state-dependent recovery rates, where we compare recovery rates based on power-law distributions to constant recovery rate models. Through this comparison, we demonstrate advantages of recovery rates based on a power-law distribution, including the potential for scaling invariance, the finite contribution of infected individuals to disease spread and the potential for deterministic disease burnout.

Methods

We now outline how to obtain state-dependent and time-varying recovery rates in ODE models. To accomplish this, we illustrate the connection between integral equation models to survival functions and finally ODE models.

Integral equation models of infectious diseases

We considered a population divided into the proportion of susceptible , infected and recovered individuals, with , which follow an integral equation formulation of an SIR compartmental model:Here is the force of infection, is the population birth rate minus the population death rate, and is the waiting-time distribution for the time spent in the infectious class before recovery, where is the time of infection.

As characterizes the time waiting in the infectious class before recovery, it also may be interpreted as a survival function. Thereby, possesses all the properties of a survival function, including corresponding to cumulative distribution function, being a hazard function, and as the average waiting-time (Arino and van den Driessche, 2006, Rodriguez, 2007).

Note, as , we omit entirely the integral equation, integro-differential equation (IDE) and ODE formulations of for the ease of presentation.

The probability of remaining infected

To generalize constant recovery rates, we now provide insight into the formulation of . We begin by expressing in terms of the recovery rate. Taking to be the probability of recovery during the time interval , it follows thatwhere is a time-dependent and/or state-dependent recovery rate. Thus, if we want to know , where is some small increment of time, then

If is constant and time invariant, then , andwhere is the constant recovery rate.

By rearranging the terms in (3) and letting , the ODE formulation of is obtained:and thus,

Alternatively, if is dependent on factors, such as the state of the epidemic, or time, then where is the recovery rate at time . It once again follows from (2) that

Rearranging terms and letting yields the ODE formulation of :and thus,

Equation (6) is not the only alternative to the waiting-time distribution (4). For instance, takingwhere is the constant waiting-time in the infected state, is a common assumption used to motivate delay-differential equations (Arino & van den Driessche, 2006).

Differential equation models of infectious diseases

To obtain ODE models of infectious diseases, we start by differentiating (1) to obtain the IDEs:

From system (7), we obtain ODEs by imposing assumptions on , and then using system (1) to eliminate the integral terms.

Constant recovery rates

To derive ODE models with a constant recovery rate from integral equation models, we make use of the standard assumptions that Thus it follows that Thereby system (7) simplifies to

From system (8), the standard ODE formulation of an SIR compartmental model is obtained by further substituting system (1), which yields:

State-dependent and time-varying recovery rates

Alternative assumptions on leads to formulation of other recovery rates. For instance, as in (5). To obtain state-dependent and time-varying recovery rates, first note that (6) reduces the IDEs from system (7) to

The IDEs from system (10), through the substitution of system (1), reduce to the ODEs,with recovery rate .

The reproductive numbers

The next-generation method is a standard technique to determine the basic reproduction number, (Diekmann et al., 1990, Heffernan et al., 2005, Liu et al., 1987, van den Driessche and Watmough, 2000, van den Driessche and Watmough, 2002). Traditionally, the next-generation method is applied to system (10), with and linearizing the flow of newly infected individuals, , and the flow from the infected state , at the disease-free equilibrium From the linearized flows,

The largest eigenvalue of the next-generation matrix, is taken as , which in this case corresponds to .

In general, system (9) with arbitrary and need not be linearizable (Greenhalgh, Galvani, & Medlock, 2015), which impedes the use of the next-generation method. Thus, we consider computing the effective reproduction number directly asand the basic reproduction number as,

If the assumption of a constant is still imposed, it follows that the force of infection is of the form:

For the force of infection (12) to be biologically meaningful (Liu et al., 1987), it follows that

Such a condition ensures the infectious period is well-posed near the disease-free equilibrium and provides a natural connection to the null probability of transmission in the absence of infection (Korobeinikov and Maini, 2005, Ponciano and Capistrán, 2011). It also follows thatandHere, such conditions placed on the force of infection amount to bounds on how quickly changes relative to the infected proportion.

Novel recovery rates

To demonstrate model predictions with time-varying and state-dependent recovery rates, we model measles outbreaks in Reykjavik Iceland. To do this, we use a classical data set on the infectious period of measles (Simpson, 1952), time series data of measles incidence in Reykjavik Iceland (Cliff et al., 1981) from 1924 to 1928, and parameter estimates from the literature (Table 1).

Table 1

Parameter values.

Definition	Parameter	Point estimate	Reference
Transmission rate (time-varying)	β_TV	0.2374/day	Fit
Rate of demographic turnover	b˜	0.027/year	Statistics Iceland
Average infectious period	μ	8.00 days	(Cliff et al., 1981, Simpson, 1952)
Standard deviation of infectious period	σ	2.26 days	(Cliff et al., 1981, Simpson, 1952)
Skewness (μ_3) and kurtosis (μ_4) of infectious period distribution:			(Cliff et al., 1981, Simpson, 1952)
Estimates from data
	μ_3μ_4	0.39days2.92days
Best fits of Pearson distributions
Type I	μ_3μ_4	0.23days2.61days
Type III	μ_3μ_4	02.1days09.6days
Type IV	μ_3μ_4	00.1days3.03days
Type V	μ_3μ_4	0.66days3.82days
Type VI	μ_3μ_4	4.13days46.5days

First, we apply our methodology to obtain time-varying recovery rates that capture the first four statistical moments (mean, standard deviation, skewness and kurtosis) of infectious period data (Schuster, 2016, pp. 83–197). To do this, we assume the infectious period distribution belongs to the family of Pearson distributions. The family of Pearson distributions consists of types I, III, IV, V, and VI, which encompass nearly all classical continuous probability distributions, including normal, Student's t, exponential, gamma, beta, inverse gamma, beta prime, and F distributions. The family of Pearson distributions also has the distinction of being the first to fit the higher statistical moments of skewness and kurtosis arbitrarily well. Importantly, this ability implies that one can specify any mean, standard deviation, skewness and kurtosis and find a Pearson distribution that describes that probability density (Pearson, 1893, Pearson, 1895, Rhind, 1909) (Fig. 1).

Fig. 1

Square of skewness and kurtosis plane. The region for each member from the family of Pearson distributions is labelled, in addition to the impossible zone () and sampling region for the variance-based sensitivity analysis.

Using the family of Pearson distributions, we develop recovery rates that account for any mean, standard deviation, skewness and kurtosis of an infectious period distribution. We evaluate how uncertainty in these moments impacts predictions on the timing and magnitude of the epidemic peak through a variance-based sensitivity analysis (Sobol, 2001). To do this, we make use of a data set on the infectious period of measles (Simpson, 1952), and determine a sampling region for skewness and kurtosis based on the best fits of each member of the family of Pearson distributions (Fig. 2) to the observed values of skewness and kurtosis from the data.

Fig. 2

Best fit to Pearson distribution. Least-squares fit of the family of Pearson distributions to infectious period data of measles.

Next, we demonstrate how to obtain state-dependent recovery rates based on a power-law distribution, outlining how such a recovery rate leads to a finite contribution to the spread of disease from infected individuals. We then demonstrate how different power-law exponents affect the trajectory of disease, through a comparison to the traditional model with a constant recovery rate.

Time-varying recovery rates associated to the Pearson distribution

To obtain time-varying recovery rates for ODEs, we take and impose that the infectious period distribution is distributed according to the family of Pearson distributions, which are given by the relation:whereand is a scale parameter. Note, the parameters and correspond to the standard deviation, square of skewness and kurtosis of the distribution respectively.

Each member of the family of Pearson distributions corresponds to particular assumptions placed on parameters and . For each case of these assumptions we fit the resulting cumulative distribution function (), through a least-squares procedure, to infectious period data of measles (Fig. 2). After obtaining the best fit for each member of the family of Pearson distributions, we estimate the conditional waiting-time distribution given time bywhere is a random variable for the time spent in the infectious class, is the time spent in the infectious class and corresponds to the bi-weekly measles incidence data from Reykjvik Iceland.

As equation (14) is a survival function, integrating over all durations of infectiousness yields a series of constant residual waiting-times (Finkelstein, 2002):

Approximating for all through linear interpolation between each , we finally determine from (Finkelstein, 2002)

The sensitivity of model predictions to an infectious period distribution

To evaluate the sensitivity of model predictions to the uncertainty in the mean, standard deviation, skewness and kurtosis from an infectious period distribution, we calculate the variance-based sensitivity indices for predicting the time to, and magnitude of the measles epidemic peak.

We now outline the procedure to conduct the analysis. First, we estimate a skewness of and a kurtosis of directly from the infectious period data. We sample the ball of radius centered at in the − plane (Table 1). The value of is the average distance from to the point corresponding to the best fit parameters of each Pearson distribution type (Table 1, Fig. 1). Additional details of the sample distributions and parameters are available in Table 1. Next, for evaluating the sensitivity of model predictions relative to the uncertainty in model inputs, we measure 1) the predicted time to the measles epidemic peak, 2) the predicted magnitude of the measles epidemic peak, and 3) the combined time to, and magnitude of, the measles epidemic peak, as given by:where and is the total population size.

The calculation of the variance-based sensitivity indices yields that uncertainty in the skewness of the infectious period distribution is the largest contributor to sensitivity in both predictions of the magnitude of the measles epidemic and the combined scenario (Table 2).

Table 2

First-order sensitivity indices.

Scenario	Mean	Standard deviation	Skewness	Kurtosis
Magnitude of epidemic	0.013	0.013	0.244	0.071
Timing of epidemic	0.008	0.006	0.287	0.091
Timing and magnitude of epidemic	0.066	0.019	0.241	0.138

State-dependent recovery rates associated to the power-law distribution

To obtain state-dependent recovery rates we follow a similar approach to one that connects constant rates to exponential distributions. In particular, one assumes that no inflow into the infectious class occurs given a certain infectious proportion at some initial time (Martcheva, 2015, pp. 9–31). Under such conditions with an exponentially distributed infectious period, the infectious proportion from system (11) behaves according to a linear ODE with coefficient . The solution to this ODE corresponds to the exponential decay, and as , with for any . The waiting-time distribution stems from normalizing by , which leads to . Note, due to the exponential decay for any .Physically speaking, having for all finite infectious periods is unrealistic, even in the most extreme of conditions.

Under the same scenario of no inflow into the infectious class, an alternative to the exponential distribution that avoids the pitfall of for all finite infectious periods is taking to correspond to the power-law distribution,Here and . Once again, making use of the assumption we obtain

In addition, the mean and standard deviation for such a , are:

Similarly to the exponential distribution case we have as . However, now for .

To determine the recovery rate that corresponds to this choice of , we isolate for yielding,

Making the substitution and rearranging the factor yields

Through differentiation and simplification, it follows thatwhere the state-dependent recovery rate is

The epidemic wave from an infectious period that is power-law distributed

We now demonstrate how a power-law distributed infectious period affects the shape of the predicted epidemic wave. As we are only concerned with effects of a single epidemic wave, we take with an average infectious period of (Table 1). It follows from the formulas for the mean and standard deviation (17) thatand

For to be biologically valid requires . Violation of this interval results in 1) being undefined at when , or 2) the associated survival function (16) taking positive values before the initial moment of infection , as forces . In addition, for a state-dependent recovery rate to be well-posed over the course of an entire epidemic the parameters of the power-law distribution need to be selected to account for the largest occurring infected proportion . instead of just any infectious proportion at some time .

This proportion is precisely the endemic equilibrium value of , or when no such equilibrium exists, the value of along a trajectory when and thus .

Finally, taking as in (12), yields

Because demographic factors are not included it follows that

As the peak of infection occurs for , as this enforces , the magnitude of the peak of infection is

Further imposing reduces system (11) to a linear ODE system. From the solution of the linear ODE system we determine the time of the peak of infection,

With (20), (21) and the measles incidence in Reykjavik Iceland (Cliff et al., 1981) from 1924 to 1928, we determine so that the predicted epidemic peak corresponds to the peak of infection for and standard deviations of an infectious period that is power-law distributed (Fig. 3).

Fig. 3

The shape of the epidemic wave. Measles time series (black dashed line) and measles incidence for , with a) , b) and c) .

Discussion

We derive infectious disease models with both time-varying and state-dependent recovery rates. To do this, we illustrate how integral equations reduce to ODE models through their connection with survival functions. We first demonstrated the novelty of time-varying recovery rates by quantifying the contribution of uncertainty in an infectious period distribution's moments to the sensitivity in predicting the epidemic peak of Iceland's 1924 measles outbreak. Our findings show that skewness had a sensitivity index at least double the other moments (Table 2). This suggests that parameter estimation techniques that rely on fitting infectious period distributions should pay close attention to the resulting value of skewness if the prediction of epidemic peaks is the ultimate goal. Secondly, we demonstrated state-dependant recovery rates derived from a power-law distribution. For such rates, we illustrated how the parameters associated with the infectious-period distribution influence the shape of an epidemic wave.

To explore the impact of uncertainty on model predictions, we used the family of Pearson distributions because it spans all possible skewness and kurtosis values (Pearson, 1893, Pearson, 1895). However, alterative distribution families, such as the exponentially generalized Beta distribution, which is often used in estimation of hazard functions, or distributions associated with nonstandard parametric survival functions (Lawless, 2002) also provide further avenues for exploring the overall effects of uncertainty on model predictions.

An important extension of our approach in developing time-varying and state-dependent rates is in the evaluation of deterministic models that compare treatment outcomes, such as drug treatments for malaria (Greenhalgh et al., 2015b, Tarning et al., 2012) or vaccine interventions. These models often reduce collected treatment and vaccination data to mean values in order to determine rate parameters. While models that compare treatment outcomes in such a fashion are useful, incorporating higher moments, like skewness and kurtosis, yields models that use more information and thus stand to provide more powerful predictions.

Many physically motivated functions merit use as a state-dependent recovery rate. We demonstrate how to obtain a state-dependent recovery rate from a power-law distribution, however other distributions with inverse cumulative density functions, such as the Weibull and Kumaraswamy distributions (Jones, 2009), would naturally lead to novel state-dependent recovery rates. Specifically, state-dependent recovery rates based on distributions that capture the recovery process, such as those with a finite infectious period or those that more accurately capture the concentration around the mean infectious period, provide a more powerful description of how infected individuals contribute to the spread of disease. Such state-dependent recovery rates would improve the modeling of quarantine, treatment or other policies that curtail the ability of infected individuals to spread disease.

Another benefit of state-dependent and time-varying recovery rates is the potential for new exactly solvable systems. Exactly solvable systems, such as the typical Susceptible-Infected models, are exceptionally rare and exceptionally useful in disease modeling. Thus, any new such system achieved from state-dependent and time-varying recovery rates would open the door to a plethora of new analytical insights.

The development of time-varying and state-dependent recovery rates stands to benefit the study of disease elimination and eradication. Disease elimination and eradication are currently important issues to world health authorities (Hopkins, 2013), particularly because of the partial success of the Global Malaria Eradication program (Nájera, González-Silva, & Alonso, 2011) and the successful eradication of both smallpox (Fenner, Henderson, Arita, Jezek, & Ladnyi, 1988) and rinderpest (Normile, 2008). For diseases slated for eradication, such as measles, yaws and polio, models that accurately capture the physical process of disease spread stand to provide quantitative information on policies and interventions designed to expedite the elimination and eradication process. As time-varying and state-dependent recovery rates better account for the contribution of infected individuals to disease spread and also include the potential for disease burnout (Greenhalgh, Galvani, et al., 2015), they thereby provide invaluable information absent from more traditional modeling formulations.

Conclusions

Approaches to accurately account for non-constant rates in deterministic epidemic models is recognized as an open challenge in infectious disease modeling (Roberts, Andreasen, Lloyd, & Pellis, 2015). Here, we provided an approach to account for time-varying and state-dependent rates developed from infectious period data. By basing time-varying and state-dependent recovery rates on infectious period distributions derived from such data, ODE models of infectious diseases are able to more accurately account for the dynamics of transmission. While we developed ODE models of infectious diseases with non-constant rates based on such infectious period data, our approach is easily generalizable to novel rates from alternate sources, including viral load, latent period and pharmacokinetic data. In summary, our work enables disease and compartmental modellers to shed the shackles imposed by exponentially distributed waiting times and constant rates by enabling the discovery and use of new non-constant rates that more accurately capture the dynamics of disease.