An Epidemic Model with Time-Distributed Recovery and Death Rates.

A compartmental epidemiological model with distributed recovery and death rates is proposed. In some particular cases, the model can be reduced to the conventional SIR model. However, in general, the dynamics of epidemic progression in this model is different. Distributed recovery and death rates are evaluated from COVID-19 data. The model is validated by the epidemiological data for different countries, and it shows better agreement with the data than the SIR model. The time-dependent disease transmission rate is estimated.

Introduction

The world population undergoes successive epidemics of viral infections with important health, social and economical consequences. During the last decades these were SARS epidemic in (Anderson et al. 2004; Lam et al. 2003), H5N1 influenza in 2005 (Chen et al. 2006; Kilpatrick et al. 2006), H1N1 influenza in 2009 (Girard et al. 2010; Jain et al. 2009), Ebola in 2014 (Frieden et al. 2014; WHO Ebola Response Team 2014), and currently, COVID-19 pandemic which continues already two years and further evolution of which remains unpredictable.

Nowadays, importance of mathematical modeling in epidemiology is generally accepted, but the outcome of this modeling is controversial. On one hand, classical compartmental models allow the evaluation of the main tendencies of epidemic progression. There are various developments of the epidemic models to multicompartment models (see, e.g., Brauer 2008; Giordano et al. 2020; Sharma et al. 2020), models with time varying or nonlinear disease transmission rate (d’Onofrio et al. 2020; Sun et al. 2008). Multipatch models (Bichara and Iggidr 2018; Lahodny and Allen 2013; McCormack and Allen 2007), multigroup models (Elbasha and Gumel 2021), spatiotemporal models (Ahmed et al. 2019; Filipe and Maule 2004) have been formulated and studied to understand various aspects of epidemic growth (see the monographs (Brauer et al. 2019; Martcheva 2015) and review articles (Hethcote 2000; Hurd and Kaneene 1993) for further details). However, the main questions about the prediction of epidemic outbreaks and their efficient managing remain unsolved. This can be partially explained by the unpredictable emergence of new viruses or virus variants, but the lack of understanding of epidemic progression and its economical consequences in a complex multiconnected modern society leads to an empirical try and error method clearly illustrated during COVID-19 pandemic (Supino et al. 2020).

Ongoing pandemic stimulated important modeling efforts directed to the application of the existing models and to their critical rethinking. Compartmental epidemiological models, like the classical SIR model, are based on the assumptions that newly infected individuals at time t appear with the rate proportional to the product of the numbers of susceptible individuals S(t) and infected individuals I(t) and that the recovery and death rates are proportional to the number of infected individuals. The first assumption is justified for homogeneous populations, but the second assumption has a limited applicability. Indeed, assuming that an average disease duration is , we conclude that the recovery and death rates at time t are determined by the number of infected individuals at time (disease onset), which can be very different from I(t), unless the epidemic progression is slow (basic reproduction number is close to 1). In a more detailed description, we do not consider a fixed disease duration but take into account that the recovery and death rates depend on the disease status of infected individuals, that is, on time passed after the disease onset.

The recovery and death rates can significantly vary depending on the individual disease progression (Github 2022). These factors are rarely taken into account in mathematical models (Feng et al. 2007; Hethcote and Tudor 1980), and further studies are needed to enlighten the significance of immunological factors to capture the incubation period (Culshaw et al. 2003; Leclerc et al. 2014; Vargas-De-León 2012), time-dependent immunity (Kyrychko and Blyuss 2005; Taylor and Carr 2009; Yuan and Bélair 2014) and other factors. Continuous dependence of the disease transmission rate on immunological parameters (e.g., instantaneous viral load) is also incorporated and studied using continuous time delay models (Gilchrist and Sasaki 2002).

In this work, we continue to study the influence of the disease time course on the epidemic progression. We propose a compartmental model based on integro-differential equations where we take into account that recovery and death rates at time t depend on the time interval from the disease onset for the individuals infected at time . We determine the main features of epidemic progression and show that they are different in comparison with conventional compartmental models. Further, we illustrate the application of this modeling approach to the COVID-19 data.

The contents of the paper are as follows. In Sect. 2, we introduce the model and study the positiveness of solutions. Next, we show how it can be reduced to the conventional SIR model in some particular cases. In order to apply this model to the investigation of COVID-19 epidemic, we determine time-dependent recovery and death rates from the available data (Sect. 3). We compare the characteristics of epidemic progression in the data and in different models in Sects. 4 and 5. Time-dependent disease transmission rate is estimated in Sect. 6.

Model with Distributed Recovery and Death Rates

Recovery and death rates of infected individuals depend on time after the disease onset. In this section, we will derive a model based on the number of newly infected individuals and their recovery and death rates depending on time after infection. We will study some properties of this model and will show that conventional SIR model can be obtained from it under some particular assumptions.

Model Formulation

We propose an integro-differential equation model where the recovery and death rates depend on the time-since-infection of the infected individuals. Let J(t) be the number of newly infected individuals at time t, while S(t), I(t), R(t) and D(t) denote the total numbers of susceptible, infected, recovered and dead individuals at time t. The total number of infected at time t is given by the following expression:We assume that the total population size remains constant, , that is, natural natality and mortality rates are assumed to be equal to each other. Using the equality and differentiating equality (1), we obtain: . On the other hand, the rate of change of the susceptible population is given by the equationwhere is the disease transmission rate.

Let and be the recovery and death rates depending on the time-since-infection . Then the number of infected individuals who will recover at time t is given by the expression:Similarly, we determine the number of infected individuals who will die at time t:Thus, the rate of change of the infected compartment I(t) is given by the following equation:The rates of change of the recovered R(t) and the death compartment D(t) are given, respectively, by the equations:Hence, we obtain the following model: where . This system of equations should be completed by the initial condition , , , and for . We will study below some properties of this model and will apply it to assess the epidemic progression.

The proposed model is capable of capturing the features of multicompartment models consisting of symptomatic and asymptomatic compartments implicitly. Their explicit consideration assumes that the individuals belonging to two compartments have different strength of infectivity and difference in time required to recovery. We explain below that r(t) and d(t) follow gamma distribution. Without any loss of generality, we can assume that asymptomatic individuals can recover much earlier than symptomatic individuals. The distributed recovery rate takes care of the time difference between the recovery of individuals belonging to two different compartments. Multicompartment epidemic models for COVID-19 also include exposed compartments and they are less infectious than the infected individuals. This aspect is taken into account by calculating the rate of infectivity from the time series of daily infected. Available data for COVID-19 infection do not differentiate between exposed, symptomatic and asymptomatic infected individuals; hence, we can consider them as a single compartment (Ghosh et al. 2022).

Positiveness of Solutions

Since Eq. (2b) contains negative integral terms, we should verify that the solution of system (2a)–(2d) remains positive. From (2a), we observe that, if at some time then . This shows that for all . From (2c), (2d) we get that R(t), D(t) are increasing functions. Hence, R(t) and D(t) also remain positive for all t. Next, we prove that . Take some . Then from (1) we haveIntegrating (2c), (2d) from 0 to with and taking their sum, we get the equalityChanging the order of integration, we obtainSince the integral gives the proportion of recovered and dead individuals from time to among those infected at time , it follows that it is less than 1. Consequently,Together with (3), this inequality gives that . Therefore, I(t) remains positive for all t. This conclusion completes the proof of positiveness of solution of system (2a)–(2d).

Reduction to the SIR Model

In this section, we show that model (2a)–(2d) can be reduced to conventional SIR model under some assumptions. Consider the recovery and death rates in the formwhere is disease duration and and are some constants. Substituting these functions in (2c) and (2d), we getIntegrating these equalities from to t, we obtainSince we assume that the disease duration is , then (1) can be written as follows:where and represent the number of recovered and dead during the time interval , respectively. Hence, from (7), we haveNow, from (2b) and (8),Assuming that is small enough, we neglect the term involving . Hence, we obtainIn this case, system (2a)–(2d) is reduced to conventional SIR model Thus, assuming uniform distribution of recovery and death rates (5) and that they are small enough, we can reduce model (2a)–(2d) to the classical SIR model. However, in general, these assumptions do not hold, and we need to take into account more realistic recovery and death rate distributions.

Estimation of Recovery and Death Rate Functions

Gamma Distribution

In this section, we estimate the recovery r(t) and death d(t) rate functions in the case of COVID-19 epidemic. The data for 120 recovered patients and 31 dead individuals from Ref. (Github 2022; Verity et al. 2019) were used to fit a gamma distribution. Note that there are no recovery or death during the first two days after infection. The maximums of these distributions are reached between 13 and 18 days for recovery and 10-15 days for death (Fig. 2). For some individuals, the recovery time is quite long. These distribution functions for recovery and death take into account asymptomatic, symptomatic and hospitalized compartments. Individuals recovered within 7 to 10 days from infection can be considered as asymptomatic either due to less viral load or due to strong immune response. On the other hand, death after significant time period from the day of infection can be assumed to be contribution from the hospitalized compartments. Further, in the literature on epidemic modeling, the choice of gamma function to model distributed recovery period is well known (Bailey 1954; Chowell et al. 2009; Lloyd 2001).

Fig. 2

Probability distribution of recovery a and death b as a function of time (in days) after the onset of infection. The red curves show the best fit (the values of parameters are given in the text) (Color figure online)

We estimate the mean time from the disease onset to recovery as 17.85 days and the mean time to death as 13.19 days. The best-fitted gamma distribution corresponding to recovery, which is shown by the red curve in Fig. 1a is given by the expressionwith the estimated parameter values and . Similarly, the best-fitted gamma distribution corresponding to death shown by the red curve in Fig. 1b is given by:where and .

Fig. 1

Probability distribution of recovery a and death b as a function of time (in days) after the onset of infection. The red curves show the best fit gamma distributions (the values of parameters are given in the text) (Color figure online)

These functions are normalized in such a way that the total probability of recovery and death equals 1. We set , where is the survival probability. Its value is estimated from the data as (Paul and Lorin 2021).

Bimodal Gamma Distribution

Instead of the gamma distribution, some other distribution functions can be used to describe the recovery and death rates. It is observed that in some cases there are two groups of recovered (dead) individuals, where one group has a shorter time interval to recovery (death) and another group a longer time period. In such cases, to obtain a better parametrization of the recovery and death rate functions, we can consider a bimodal gamma distribution, that is, a linear combination of two different gamma distributions. In (Paul and Lorin 2021), the distribution of recovery and death as functions of onset-to-recovery and onset-to-death are estimated using the COVID-19 data in Canada. The corresponding data are shown in Fig. 2 by the blue bars. We have fitted these data by bimodal gamma distributions and (red curves in Fig. 2) corresponding to recovery and death rate functions, respectively:where the best-fitted parameter values are as follows: , , , and , , ,

Sensitivity Analysis

Parameters in the recovery and death distributions presented above are estimated from the individual data which can vary depending on country, time period, and on the virus variant. We will estimate the sensitivity of the model outcomes (maximal number of infected and time to the maximal number of infected ) to the shape and scale parameters , , , . For this purpose, we use variance-based sensitivity analysis with the Monte Carlo numerical procedure described in (Saltelli et al. 2008) for computing the full set of first-order sensitivity indices for corresponding to the parameters , , and , respectively. We have estimated the parameters , , and from the individual level data given in (Github 2022). Then we use a set of sample points obtained by using Latin hyper-cube sampling in the neighborhood of these estimated parameter values and perform numerical simulation as described in (Saltelli et al. 2008).

The first-order sensitivity indices are shown in Fig. 3 and summarized in Table 1. This sensitivity analysis shows that the model outcomes (maximal number of infected) and (time to the maximal number of infected) are most sensitive to the scale parameter in the gamma distribution for the recovery rate.

Fig. 3

The left panel shows the first-order sensitivity indices corresponding to the model outcome and the right panel to in the case of gamma distribution (Color figure online)

Table 1

First-order sensitivity indices (gamma distribution)

Parameters	a_1	a_2	b_1	b_2
Sensitivity of I_m	0.1940	0.0068	0.7905	0.0059
Sensitivity of t_m	0.0922	0.1312	0.4971	0.2208

A similar method is used to perform the sensitivity analysis for the parameters , , , , , , and involved in the bimodal gamma distribution. The corresponding first-order sensitivity indices are shown in Fig. 4 and summarized in Table 2. We can observe that is the most sensitive parameter to and as compared to other parameters.

Fig. 4

The left panel shows the first-order sensitivity indices corresponding to the model outcome and the right panel to in the case of bimodal gamma distribution (Color figure online)

Table 2

First-order sensitivity indices (bimodal gamma distribution)

Parameters	a_1	a_2	c_1	c_2	b_1	b_2	d_1	d_2
Sensitivity of I_m	0.0807	0.0188	0.0195	0.0191	0.8051	0.0244	0.0193	0.0194
Sensitivity of t_m	0.0480	0.0646	0.0304	0.0413	0.6291	0.0816	0.0435	0.0625

Comparison with the SIR Model

We showed in Sect. 2 that classical SIR model can be obtained as a particular case of distributed model (2a)–(2d). We will compare dynamics of epidemic progression in the two models using the estimated recovery and death rates.

Since the estimated average time to recovery is 17.85 days and to death 13.19 days, we take average disease duration as 16 days. The corresponding value in SIR model is . We set , that is, out of 100 infected individuals, 97 infected will recover. This estimate matches with most of the COVID-19 epidemic data from various countries (Worldometer 2022; Paul and Lorin 2021).

Though the parameters of the two models correspond to each other, system (2a)–(2d) and SIR model (10a)–(10c) give different dynamics of epidemic progression (Fig.  5). We notice that the maximal number of infected individuals I(t) is much higher for system (2a)–(2d) as compared to the SIR model (10a)–(10c), while time to the maximal number of infected is less.

Fig. 5

Comparison between the solutions of the system (2a)–(2d) (blue curves) and SIR model (red curves): a the number of susceptible individuals S(t), b the number of infected individuals I(t). In both models . The values of other parameters for the SIR model: , , ; and for model (2a)–(2d): , , , , , , , (Color figure online)

Comparison of the final size of epidemic , maximal number of infected and the time to the maximal number of infected between system (2a)–(2d) with gamma distribution and the SIR model are shown in Fig. 6 for different values of parameters. As before, the maximal number of infected individuals in model (2a)–(2d) is much higher than for the SIR model (10a)–(10c), time and the final size are less for the distributed model. This difference can be explained by the fact that the recovery and death rates are uniformly distributed during the disease duration for the SIR model (Sect. 2), contrary to the gamma distribution in (2a)–(2d). Therefore, there is a shift to earlier recovery and death for the SIR model.

Fig. 6

Comparison of the maximal number of infected individuals a, time to reach the maximal number b and the final size of epidemic c between the system (2a)–(2d) (blue curves) and the SIR model (red curves) for different values of . The values of parameters: , , for the SIR model ; and for the system (2a)–(2d): , , , , (Color figure online)

Similarly, we compare system (2a)–(2d) with bimodal gamma distribution with the conventional SIR model. In this case, the estimated mean time from onset-to-recovery is 22.63 and the mean time from onset-to-death is 21.80. Hence, the average disease duration is taken approximately 22.2 days, and . All other parameters are kept the same as above. The properties of the final size of epidemic, maximal number of infected individuals, time to maximum are similar to the previous case and not shown here.

Model Validation with Epidemiological Data

In order to validate the model with distributed recovery and death rates, we compare the results of modeling with the epidemiological data. Distributed recovery and death rates are estimated in Sect. 3.1 from the data in China in February 2020 (Github 2022; Verity et al. 2019). Once these functions are determined, we take the number J(t) of daily infected individuals from the epidemiological data and find the sum of daily recoveries and deaths from the expressionThese results are compared with the sum of recoveries and deaths in the data. Figure 7 shows the result of such comparison for China from January 23, 2020, to April 15, 2020, with the data from (Worldometer 2022) (7-day moving average).

Fig. 7

In the left panel, the blue curve shows the number of recovered and dead in the distributed model, the magenta curve corresponds to in the SIR model, and the black dots correspond to the 7-days moving average of daily recoveries and death in China. The right panel shows the corresponding cumulative recovery and death (Color figure online)

Recoveries and deaths can also be determined as a proportion of active cases as it is done in the SIR model. Here I(t) is taken from the data and . In agreement with the results of the previous section, SIR model overestimates the sum of recovered and dead.

A similar comparison is done for other countries (Fig. 8). It is important to mention here that we used the same gamma distribution as determined before from the data for China (Sect. 3.1).

Fig. 8

, are plotted for different countries, using the gamma distributions for recovery and death rates as estimated in Sect. 3.1. In the left panel, the blue curves correspond to , the magenta curves correspond to and the black dots correspond to the 7-day moving average of daily recovery and death in different countries. The right panel shows the corresponding cumulative recovery and death. a, b: France; c, d: Italy; e, f: Sweden (Color figure online)

Next, we consider the bimodal gamma distribution determined above (Sect. 3.2) and calculate for a longer period of time from March 10, 2020, to June 16, 2020 (Fig. 9), than used for the determination of the distribution parameters. As before, we compare the results with the SIR model and observe that it overestimates the total recovery and death.

Fig. 9

In the left panel, the blue curve corresponds to , the magenta curve corresponds to and the black dots correspond to the 7-day moving average of daily recovery and death in Canada. The right panel shows the corresponding cumulative recovery and death (Color figure online)

Thus, the model with gamma distribution gives a good description of the recovery and death in different countries compared with the epidemiological data, while the SIR model overestimates it.

Estimation of the Time-Dependent Disease Transmission Rate

In this section, we will consider time-dependent transmission rate, and will estimate it from the COVID-19 epidemiological data. Dynamics of the transmission rate can help in the understanding of epidemic progression (Mummert 2013).

Theorem 1

For the model (2a)–(2d), the time-dependent transmission function is given by the following expression:

Proof

We haveNow, we also know thatUsing in the previous equation, we getSubstituting this expression into (13), we obtain (12).

As an illustration of this theorem, we consider the COVID-19 data taken from (Worldometer 2022) for new daily cases and total active cases. In order to decrease the irregularity of data, we take the 7-day moving average of J(t) and I(t). Note that represents the cumulative number of infected. Consequently, we can determine the function using formula (12).

We consider the COVID-19 infection data for a span of approximately 450 days. We use the data for India from March 7, 2020, for France from March 2, 2020, for Italy from February 21, 2020, and for Sweden from March 23, 2020, and up to May 20, 2021, for all the four countries (Worldometer 2022). Then we plot for four countries with the help of (12) and plot in Fig. 10. The initial date corresponds to the first reported case in a given country (marked with vertical dashed lines). Initial transient observed in case of India may be related to the inaccuracy of the reporting strategy. Growth and decline in at different time correspond to various social restrictions as well as onset of a new outbreak. It is interesting to note that the declining pattern for two neighboring European countries, France and Italy, are similar in the beginning. However, an increasing peak for may indicate that the relaxation of lockdown restriction in France was more rapid compared to other countries.

Fig. 10

Time varying for COVID-19 in France, India, Italy, Sweden calculated by formula (12) for the model with gamma distribution and the parameter values: , , , , (Color figure online)

We can note from the presented results that oscillates according to increasing or decreasing epidemic waves. Furthermore, average value of this function is different in different countries. As such, it is about 0.1 in India and about 0.05 in France. In order to interpret dynamics of time-dependent transmission rate , we simplify expression (12) assuming that . This assumption is justified since the total number of infected remains in most countries much less than the total population. Then . The same expression for can be obtained from the SIR model if .

In order to give further estimates of , suppose that disease duration is . Then , that is, the individuals infected at time recover or die at time t but not before. Hence, we obtain approximate formulaSet and then substituting in above equation, we find . If , then , that is, the disease transmission rate is inversely proportional to the disease duration. This estimate is in agreement with an average disease duration 16 days determined in Sect. 2. For , characterizes the rate of growth of newly infected individuals, and for negative , the rate of decay.

Discussion

Ongoing COVID-19 pandemic has stimulated scientific research in various disciplines ranging from economy to education (Volpert et al. 2020). A wide variety of modeling approaches are considered in the recent literature (see, e.g., Rahimi et al. 2021; Sharma et al. 2020 for more detail). However, validation of these models is complicated by the uncertainty of the data, especially for asymptomatic cases.

Another shortcoming of conventional epidemiological model is that they consider recovery and death rate as a given proportion of active infected cases at the same moment of time. Clearly, this is a strong assumption which can lead to a large error in the evaluation of epidemic progression. In order to overcome this issue, we propose in this work a new type of immunoepidemiological models based on the daily number of infected individuals and their time-distributed recovery and death rates. Distributed recovery and death rates are evaluated from the data from China and Canada. They give a reliable description of data for different countries and time periods. We note that the parameters of gamma distribution can depend on the virus variant. This question can be addressed in the future studies when the time-dependent recovery rate is available in particular for the Omicron variant.

We compare this approach with the SIR model with appropriate recovery and death rates. It is clearly seen that the SIR model overestimates the daily recoveries and deaths which, in turn, underestimates the daily number of infected, maximal and total numbers of infected individuals. The recovery and death rate functions are estimated with the limited real data from China and Canada on the number of days spent since infection before recovery and death (Github 2022; Paul and Lorin 2021). These data sets are used to estimate the parameters involved with two different parametrization of the recovery and death rate functions. Time-since-infection-based recovery and death rates implicitly take into account mild and severe infection which can be considered as symptomatic and asymptomatic compartments. Numerical validation of the proposed model with the COVID-19 epidemic data from five different countries indicates that the parametrization of recovery and death rates with gamma function effectively captures the daily and cumulative recoveries and deaths, although the real data show large irregularity.

It is important to mention here that the proposed modeling approach can be used to predict the disease progression accurately if we have specific data for the first days in order to estimate the parameters involved in the recovery and death rate functions. Availability of such kind of data is a challenging issue in the beginning of epidemic. However, we should highlight that the estimates of r(t) and d(t) with the data from China during the onset of COVID-19 epidemic works well to study the disease progression in France, Italy and Sweden. Having these information, the proposed modeling approach can be used to predict the maximal number of infected and the time to maximal number of infected. This predictions can be used to estimate the required number of hospital beds and readiness of medical facilities based upon the appropriate information about the rate of hospitalization and severity of the viral strain.

We have also described a method to calculate time-dependent infectivity rate based upon the daily incidence data. It clearly indicates that the rate of transmission of infection from one individual to another not only depends upon the fixed transmission rate, but rather it is solely related with the time period over which one individual remain infected. The next challenging issue will be to estimate the rate of infection transmission depending on time-dependent viral load for different virus variants.