An Epidemic Model with Time Delay Determined by the Disease Duration

Abstract

Immuno-epidemiological models with distributed recovery and death rates can describe the epidemic progression more precisely than conventional compartmental models. However, the required immunological data to estimate the distributed recovery and death rates are not easily available. An epidemic model with time delay is derived from the previously developed model with distributed recovery and death rates, which does not require precise immunological data. The resulting generic model describes epidemic progression using two parameters, disease transmission rate and disease duration. The disease duration is incorporated as a delay parameter. Various epidemic characteristics of the delay model, namely the basic reproduction number, the maximal number of infected, and the final size of the epidemic are derived. The estimation of disease duration is studied with the help of real data for COVID-19. The delay model gives a good approximation of the COVID-19 data and of the more detailed model with distributed parameters.

1. Introduction

Mathematical modeling of infectious diseases attracts much attention due to successive epidemics, such as HIV, emerging in the 1980s and still continuing [1,2], SARS epidemic in 2002–2003 [3,4], H5N1 influenza in 2005 [5,6] and H1N1 in 2009 [7,8], Ebola in 2014 [9,10]. The ongoing COVID-19 pandemic has stimulated unprecedented efforts of mathematical modeling in epidemiology. A wide variety of mathematical approaches are developed to study epidemiological problems. However, sufficiently simple and validated models still remain in the focus of mathematical modeling in epidemiology.

Modern studies in mathematical epidemiology begin with the SIR model developed in the works by W. O. Kermack and A. G. McKendrick [11,12,13], stimulated by the Spanish flu epidemic in 1918–1919. Among many developments of such models, we can cite multi-compartment models [14,15,16], models with a time-varying or nonlinear disease transmission rate [17,18], multi-patch models [19,20,21], multi-group models incorporating the effect of the heterogeneity of the population [22], and epidemic models with vaccination and other control measures [23,24]. Random movement of individuals in the population is considered in spatiotemporal models in order to describe spatial distributions of susceptible and infected individuals [25,26]. A more detailed literature review can be found in the monographs [27,28] and review articles [29,30].

The conventional SIR model, which includes susceptible (S), infected (I), and recovered (R) compartments, and similar models assume that recovery and death rates at time t are proportional to the number of actively infected individuals $I\left(t\right)$ at the same moment of time. This assumption does not take into account disease duration, and it can lead to a large error. In our previous work [31], we showed that this assumption leads to an overestimation of actual recoveries and deaths. Instead, if we use distributed recovery and death rates, properly chosen from real data, the description of the epidemic progression becomes more precise. However, since distributed recovery and death rates are not easily available, we develop a simpler delay model in this work. It gives close results, but it does not require precise immunological data. The model considered in [31] involves distributed recovery and death rates. The model considered in [32] is an extension of the model by incorporating the vaccinated compartment, and the resulting model is an immuno-epidemic model. The delay model, with disease duration delay, considered here is derived from the model proposed in [31] with an appropriate assumption on the recovery and death rate functions. The present model is quite different from the delay model considered in [33] as the delay parameter involved with the earlier model was the measure of the incubation period and the departure of infected individuals from the infected compartment was due to the imposition of quarantine measure. The model proposed and analyzed here is solely dependent on the disease duration period and without imposed quarantine.

Most of the existing delay epidemic models consider time delay either in the disease incidence function or in the susceptible recruitment function (Appendix A). The delay in the recovery and death rates has not been studied yet thoroughly. In this work, we introduce time delay in recovery and death rates with the average disease duration considered as the delay parameter.

In Section 2, we discuss the distributed model and derive the delay model where the discrete time delay estimates average disease duration. We obtain epidemic characterization of the delay model in Section 3. Then, in Section 4, we perform a numerical comparison among the distributed model, delay model, and conventional SIR model with the equivalent parameter values. Next, we discuss a method to estimate the value of disease duration using the real data of disease incidence in Section 5. In Section 6, we validate our delay model with epidemiological data collected during the COVID-19 epidemic. The main outcomes of the proposed model and its epidemiological implications are discussed in the concluding section.

2. Model Formulation

In contrast with the existing compartmental epidemiological models, we start the model derivation by the introduction of the class of newly infected individuals instead of the total number of infected individuals. This approach is appropriate to evaluate daily recovery and death rates. We recall in this section the model with a distributed recovery and death rate [31,33]. We then use this model to derive the delay model and study its properties in the next sections.

2.1. Model with Distributed Parameters

The number of newly infected individuals $J\left(t\right)$ is determined by the rate of decrease of the number of susceptible individuals, $J\left(t\right)=-S^{\prime }\left(t\right)$. Assuming that

$$N=S\left(t\right)+I\left(t\right)+R\left(t\right)+D\left(t\right)$$

(1)

is constant, where $I\left(t\right)$ is the total number of infected at time t and $R\left(t\right)$ and $D\left(t\right)$ denote, respectively, recovered and dead, we can write

$$I\left(t\right)={\int }_0^{t}J\left(\eta \right)d\eta -R\left(t\right)-D\left(t\right).$$

(2)

Following conventional epidemiological models, we set

$$\frac{dS\left(t\right)}{dt}=-\beta \frac{S\left(t\right)}{N}I\left(t\right),$$

where $\beta$ is the disease transmission rate. Let $r\left(\eta \right)$ and $d\left(\eta \right)$ be the recovery and death rates depending on the time-since-infection $\eta$. Then, the number of infected individuals who will recover at time t is given by the expression:

$$\frac{dR\left(t\right)}{dt}={\int }_0^{t}r(t-\eta )J\left(\eta \right)d\eta ,$$

and the number of infected individuals who will die at time t:

$$\frac{dD\left(t\right)}{dt}={\int }_0^{t}d(t-\eta )J\left(\eta \right)d\eta .$$

Differentiating Equality (1), we obtain

$$\frac{dI\left(t\right)}{dt}=\beta \frac{S\left(t\right)}{N}I\left(t\right)-{\int }_0^{t}r(t-\eta )J\left(\eta \right)d\eta -{\int }_0^{t}d(t-\eta )J\left(\eta \right)d\eta .$$

Thus, we obtain the following integro-differential equation model:

$$\begin{array}{ccc}\hfill \frac{dS\left(t\right)}{dt}& =& -\beta \frac{S\left(t\right)}{N}I\left(t\right)(=-J\left(t\right)),\hfill \end{array}$$

(3a)

$$\begin{array}{ccc}\hfill \frac{dI\left(t\right)}{dt}& =& \beta \frac{S\left(t\right)}{N}I\left(t\right)-{\int }_0^{t}r(t-\eta )J\left(\eta \right)d\eta -{\int }_0^{t}d(t-\eta )J\left(\eta \right)d\eta ,\hfill \end{array}$$

(3b)

$$\begin{array}{ccc}\hfill \frac{dR\left(t\right)}{dt}& =& {\int }_0^{t}r(t-\eta )J\left(\eta \right)d\eta ,\hfill \end{array}$$

(3c)

$$\begin{array}{ccc}\hfill \frac{dD\left(t\right)}{dt}& =& {\int }_0^{t}d(t-\eta )J\left(\eta \right)d\eta ,\hfill \end{array}$$

(3d)

with the initial condition $S\left(0\right)=N$, $I\left(0\right)=I_0>0$ ($I_0$ is sufficiently small as compared to N), $R\left(0\right)=0$, $D\left(0\right)=0$.

2.2. Reduction to SIR Model

In a particular case, if we assume the uniform distribution of recovery and death rates:

$$r(t-\eta )=\left\{\begin{array}{ccc}{r}_0& ,& t-\tau <\eta \le t\\ 0& ,& \eta <t-\tau \end{array}\right.,d(t-\eta )=\left\{\begin{array}{ccc}{d}_0& ,& t-\tau <\eta \le t\\ 0& ,& \eta <t-\tau \end{array}\right.,$$

where $\tau >0$ is disease duration, $r_0$ and $d_0$ are some constants, and if $r_0$ and $d_0$ are small enough, then the model (3) can be reduced to the conventional SIR model (see [31]):

$$\begin{array}{ccc}\hfill \frac{dS\left(t\right)}{dt}& =& -\beta \frac{S\left(t\right)}{N}I\left(t\right),\hfill \end{array}$$

(4a)

$$\begin{array}{ccc}\hfill \frac{dI\left(t\right)}{dt}& =& \beta \frac{S\left(t\right)}{N}I\left(t\right)-(r_0+d_0)I\left(t\right),\hfill \end{array}$$

(4b)

$$\begin{array}{ccc}\hfill \frac{dR\left(t\right)}{dt}& =& r_{0}I\left(t\right),\frac{dD\left(t\right)}{dt}=d_{0}I\left(t\right).\hfill \end{array}$$

(4c)

However, the approximation of a uniform distribution of recovery and death rates may not be precise since infected individuals are unlikely to recover or die at the beginning of the disease. In the following subsection, we consider another choice of recovery and death rates, which will approximate the real scenario of recovery and death more accurately.

2.3. Delay Model

Let us assume that disease duration is $\tau$, and the individuals $J(t-\tau )$ infected at time $t-\tau$ recover or die at time t with certain probabilities. This assumption corresponds to the following choice of the functions $r(t-\eta )$ and $d(t-\eta )$:

$$r(t-\eta )=r_0\delta (t-\eta -\tau ),d(t-\eta )=d_0\delta (t-\eta -\tau ),$$

where $r_0$, $d_0$ are constants, $r_0+d_0=1$, and $\delta$ is the Dirac delta-function. Then,

$$\frac{dR\left(t\right)}{dt}={\int }_0^{t}r(t-\eta )J\left(\eta \right)d\eta =r_{0}J(t-\tau ),$$

$$\frac{dD\left(t\right)}{dt}={\int }_0^{t}d(t-\eta )J\left(\eta \right)d\eta =d_{0}J(t-\tau ).$$

Note that $J\left(t\right)$ is the number of newly infected individuals appearing at time t. If we assume that the first infected case was reported at time $t=0$, then we can set $J\left(t\right)=0$ for all $t<0$. Now, integrating the above two equations from 0 to t and assuming that $R\left(0\right)=D\left(0\right)=0$, we obtain

$$R\left(t\right)=r_0{\int }_0^{t}J(s-\tau )ds=r_0{\int }_{-\tau }^{t-\tau }J\left(y\right)dy=r_0{\int }_0^{t-\tau }J\left(y\right)dy,$$

$$D\left(t\right)=d_0{\int }_0^{t}J(s-\tau )ds=d_0{\int }_{-\tau }^{t-\tau }J\left(y\right)dy=d_0{\int }_0^{t-\tau }J\left(y\right)dy.$$

Then, instead of (2), we have

$$I\left(t\right)={\int }_{t-\tau }^{t}J\left(s\right)ds,$$

(5)

and from (3b),

$$\frac{dI\left(t\right)}{dt}=\beta \frac{S\left(t\right)}{N}I\left(t\right)-J(t-\tau ).$$

(6)

From (5) we obtain

$$J\left(t\right)=\beta \frac{S\left(t\right)}{N}{\int }_{t-\tau }^{t}J\left(s\right)ds,\frac{dS\left(t\right)}{dt}=-J\left(t\right).$$

(7)

Hence,

$$\frac{dS\left(t\right)}{dt}=-\beta \frac{S\left(t\right)}{N}{\int }_{t-\tau }^{t}J\left(s\right)ds=\beta \frac{S\left(t\right)}{N}{\int }_{t-\tau }^t\frac{dS\left(s\right)}{ds}ds=-\beta \frac{S\left(t\right)}{N}(S(t-\tau )-S\left(t\right)).$$

(8)

Once we obtain the solution $S\left(t\right)$ from (8), then we can find $I\left(t\right)$ using the following relation:

$$I\left(t\right)={\int }_{t-\tau }^{t}J\left(s\right)ds=-{\int }_{t-\tau }^t\frac{dS\left(s\right)}{ds}ds=S(t-\tau )-S\left(t\right).$$

Hence, System (3) is reduced to the following delay model:

$$\begin{array}{ccc}\hfill \frac{dS\left(t\right)}{dt}& =& -J\left(t\right),\hfill \end{array}$$

(9a)

$$\begin{array}{ccc}\hfill \frac{dI\left(t\right)}{dt}& =& J\left(t\right)-J(t-\tau ),\hfill \end{array}$$

(9b)

$$\begin{array}{ccc}\hfill \frac{dR\left(t\right)}{dt}& =& r_{0}J(t-\tau ),\hfill \end{array}$$

(9c)

$$\begin{array}{ccc}\hfill \frac{dD\left(t\right)}{dt}& =& d_{0}J(t-\tau ),\hfill \end{array}$$

(9d)

$$\begin{array}{ccc}\hfill J\left(t\right)& =& \beta \frac{S\left(t\right)}{N}I\left(t\right),\hfill \end{array}$$

(9e)

with $J\left(t\right)=0$ for all $t<0$. A similar model was proposed in [33] without derivation from the distributed model.

3. Epidemic Characteristics

In this section, we determine the basic reproduction number, the final size of the epidemic, and maximum number of infected individuals for the delay model (9).

3.1. Basic Reproduction Number

To determine the initial exponential growth rate and the basic reproduction number ${\mathcal{R}}_0$, we use the equations for the infected compartments $I\left(t\right)$ considered in the form:

$$\frac{dI\left(t\right)}{dt}=\beta \frac{S\left(t\right)}{N}I\left(t\right)-\beta \frac{S(t-\tau )}{N}I(t-\tau ).$$

(10)

Suppose that, at the beginning of the epidemic, $S\left(t\right)\approx S_0$ and $S(t-\tau )=S_0.$ Then, from (10), we have

$$\frac{dI\left(t\right)}{dt}=\beta \frac{S_0}{N}(I\left(t\right)-I(t-\tau )).$$

Substituting $I\left(t\right)=I_0{e}^{\lambda t}$, we have

$$\lambda =\beta \frac{S_0}{N}\left(1-e^{-\lambda \tau }\right).$$

(11)

Let $\mathcal{G}\left(\lambda \right)=\beta \frac{S_0}{N}\left(1-e^{-\lambda \tau }\right).$ Clearly, (11) has a solution $\lambda =0$ and a non-zero solution with the sign determined by ${\mathcal{G}}^{\prime }\left(0\right).$ Denote

$${\mathcal{R}}_0={\mathcal{G}}^{\prime }\left(0\right)=\beta \tau \frac{S_0}{N}.$$

Let ${\mathcal{R}}_0>1.$ This implies, ${\mathcal{G}}^{\prime }\left(0\right)>1.$ We observe that $\mathcal{G}\left(\lambda \right)$ is an increasing function of $\lambda$ and $\mathcal{G}\left(\lambda \right)\to \beta \frac{S_0}{N}$ as $\lambda \to \infty .$ This implies that Equation (11) has a positive solution, i.e., there exists ${\lambda }_{*}>0$ such that $\mathcal{G}\left({\lambda }_{*}\right)={\lambda }_{*}.$ If ${\mathcal{R}}_0<1,$ i.e., ${\mathcal{G}}^{\prime }\left(0\right)<1,$ then

$${\mathcal{G}}^{\prime }\left(\lambda \right)=\beta \tau \frac{S_0}{N}{e}^{-\lambda \tau }<1$$

for all $\lambda \ge 0.$ Therefore, the equation $\mathcal{G}\left(\lambda \right)=\lambda$ has no positive solution.

Hence, the basic reproduction number is given by the following expression:

$${\mathcal{R}}_0=\beta \tau \frac{S_0}{N}.$$

(12)

Letting $\widehat{\lambda }=\lambda /\left(\beta \frac{S_0}{N}\right)$, we can write Equation (11) in the following form:

$$\widehat{\lambda }=1-e^{-{\mathcal{R}}_0\widehat{\lambda }}.$$

The solution of this equation determines the initial exponential growth rate, which depends on the single parameter ${\mathcal{R}}_0$.

3.2. Final Size of the Epidemic

Next, we determine the final size of the susceptible compartment $S_f={lim}_{t\to \infty }S\left(t\right)$. From (3a), we obtain:

$$\frac{dS\left(t\right)}{dt}=-\beta \frac{S\left(t\right)}{N}{\int }_{t-\tau }^{t}J\left(\eta \right)d\eta .$$

Then, integrating from 0 to ∞, we obtain

$${\int }_0^{\infty }\frac{dS}{S}=-\frac{\beta }{N}{\int }_0^{\infty }\left({\int }_{t-\tau }^{t}J\left(\eta \right)d\eta \right)dt.$$

Changing the order of integration, we get

$$\begin{array}{cc}\hfill ln\frac{S_0}{S_f}& =\frac{\beta }{N}\left[{\int }_{-\tau }^0\left({\int }_0^{\eta +\tau }J\left(\eta \right)dt\right)d\eta +{\int }_0^{\infty }\left({\int }_{\eta }^{\eta +\tau }J\left(\eta \right)dt\right)d\eta \right]\hfill \\ \hfill & =\frac{\beta }{N}\left[{\int }_{-\tau }^0(\eta +\tau )J\left(\eta \right)d\eta +{\int }_0^{\infty }\tau J\left(\eta \right)d\eta \right].\hfill \end{array}$$

Now, integrating (9a) from 0 to ∞, we have:

$${\int }_0^{\infty }J\left(\eta \right)d\eta =S_0-S_f.$$

(13)

Since $J\left(t\right)=0$ for all $t\in [-\tau ,0]$, then

$$\begin{array}{ccc}\hfill ln\frac{S_0}{S_f}& =& \frac{\beta }{N}\left[{\int }_0^{\infty }\tau J\left(\eta \right)d\eta \right]=\frac{\beta }{N}\tau \left(S_0-S_f\right).\hfill \end{array}$$

(14)

Thus, the final size can be obtained from the equation:

$$ln\frac{S_0}{S_f}=\frac{\beta }{N}\tau \left(S_0-S_f\right).$$

(15)

Integrating (9c) and (9d) from 0 to ∞ and using (13), we obtain:

$$R_f\equiv \underset{t\to \infty }{lim}R\left(t\right)=r_0{\int }_0^{\infty }J(s-\tau )ds=r_0{\int }_0^{\infty }J\left(\eta \right)d\eta =r_0(S_0-S_f),$$

$$D_f\equiv \underset{t\to \infty }{lim}D\left(t\right)=d_0{\int }_0^{\infty }J(s-\tau )ds=d_0{\int }_0^{\infty }J\left(\eta \right)d\eta =d_0(S_0-S_f).$$

Final size $S_f$ obtained from the numerical simulation and using the Formula (15) is verified for three different values of $\tau$ and a range of values of $\beta$ (Figure 1). Observe that Equation (15) can be written as:

$$ln\frac{S_0}{S_f}={\mathcal{R}}_0\left(1-\frac{S_f}{S_0}\right).$$

(16)

Note that the SIR model (4) gives the same equation for the final size, and also in the SIR model (4), if we assume that $r_0+d_0=1/\tau$, then the expression of ${\mathcal{R}}_0$ is equivalent for both the SIR model (4) and the delay model (9).

3.3. Maximum Number of Infected Individuals

We now derive an approximate formula for the maximal number of infected individuals for the delay model (9). We have $I\left(t\right)=S(t-\tau )-S\left(t\right)$. Suppose that $I\left(t\right)$ attains its maximum at $t=t_m$. Set $I_m=I\left(t_m\right)$, $S_m=S\left(t_m\right)$. From the equality $I^{\prime }\left(t_m\right)=0$, we obtain $S^{\prime }(t_m-\tau )=S^{\prime }\left(t_m\right)$. This implies

$$\frac{\beta S(t_m-\tau )I(t_m-\tau )}{N}=\frac{\beta S\left(t_m\right)I\left(t_m\right)}{N}.$$

(17)

Substituting the relation $I\left(t_m\right)=S(t_m-\tau )-S\left(t_m\right)$ in (17), we obtain

$$I(t_m-\tau )=\frac{S_m{I}_m}{S_m+I_m}.$$

(18)

From (8), we have

$$\frac{dS}{dt}=-\frac{\beta S\left(t\right)}{N}\left(S(t-\tau )-S\left(t\right)\right).$$

(19)

Integrating (19) from 0 to $t_m$ and changing the variable inside the first integral of the right-hand side, we obtain:

$$\begin{array}{ccc}\hfill {\int }_0^{t_m}\frac{dS}{S}& =& -\frac{\beta }{N}\left({\int }_{-\tau }^{0}S\left(t\right)dt-{\int }_{t_m-\tau }^{t_m}S\left(t\right)dt\right).\hfill \end{array}$$

(20)

We assume that $S\left(t\right)=S_0$ for all $t\in [-\tau ,0]$ and use the approximation (Figure 2):

$${\int }_{t_m-\tau }^{t_m}S\left(t\right)dt\approx \frac{\tau }{2}{I}_m+\tau S_m.$$

Then from (20) we have

$$\begin{array}{cc}\hfill ln\frac{S_m}{S_0}& =-\frac{\beta }{N}\left(\tau S_0-\frac{\tau }{2}{I}_m-\tau S_m\right).\hfill \end{array}$$

Let $x=\frac{I_m}{S_0}$, $y=\frac{S_m}{S_0}$ and $\frac{S_0}{N}\approx 1$. Then we have

$$lny=-\beta \tau \left(1-\frac{1}{2}x-y\right).$$

(21)

Again, integrating (19) from $t_m-\tau$ to $t_m$, we obtain:

$${\int }_{t_m-\tau }^{t_m}\frac{dS}{S}=-\frac{\beta }{N}\left({\int }_{t_m-\tau }^{t_m}S(t-\tau )dt-{\int }_{t_m-\tau }^{t_m}S\left(t\right)dt\right).$$

Changing the variable inside the first integral of the right-hand side, we get:

$$ln\frac{S_m}{S(t_m-\tau )}=-\frac{\beta }{N}\left({\int }_{t_m-2\tau }^{t_m-\tau }S\left(t\right)dt-{\int }_{t_m-\tau }^{t_m}S\left(t\right)dt\right).$$

Now, using the approximation described in Figure 2, we conclude that

$$ln\frac{S_m}{S(t_m-\tau )}=-\frac{\beta }{N}\left(\tau S(t_m-\tau )+\frac{\tau }{2}I(t_m-\tau )-\tau S_m-\frac{\tau }{2}{I}_m\right).$$

(22)

Using the relation $I\left(t_m\right)=S(t_m-\tau )-S\left(t_m\right)$, after some transformations, (22) can be written as:

$$ln\frac{S_m}{S_m+I_m}=-\frac{\beta }{N}\frac{\tau }{2}\left(I_m+I(t_m-\tau )\right).$$

Using (18), we obtain:

$$ln\frac{S_m}{S_m+I_m}=-\frac{\beta }{N}\frac{\tau }{2}\left(I_m+\frac{S_m{I}_m}{S_m+I_m}\right).$$

Substituting $x=\frac{I_m}{S_0}$, $y=\frac{S_m}{S_0}$ and $\frac{S_0}{N}\approx 1$, we have

$$ln\frac{y}{x+y}=-\frac{\beta \tau }{2}\left(x+\frac{xy}{x+y}\right).$$

(23)

Solving (21) and (23), we can find x, y and, consequently, $I_m$, $S_m$.

In Figure 3, we show a comparison between the maximum number of infected obtained by Equations (21) and (23) and the maximum number of infected obtained by direct numerical simulation of the delay model (9). From Figure 3, we can observe that the approximation gives a very close upper bound to the maximum number of infected.

4. Comparison of Models (3) and (9) and SIR (4)

We compare the results obtained from the distributed model (3), the delay model (9), and the conventional SIR model (4). In [31], we showed that Model (3) can describe the epidemic spread more precisely, and it can exactly capture the recovery and death dynamics by using suitable distributed recovery and death rates. However, the main constraint in using the distributed model (3) is the availability of distributed recovery and death rates. Instead, the average duration of the disease is easier to determine. We show that, in the delay model (9), if we take $\tau$ as the mean of the recovery and death distribution functions involved in the distributed model (3), then the delay model (9) gives a close solution to the distributed model (3).

To estimate the recovery and death distributions in the model (3), we use the function fitdist(:,gamma) in MATLAB. This function is used to fit a vector of data $X=(x_1,x_2,\cdots ,x_n)$ by a gamma distribution of the form $\frac{1}{b^a\Gamma \left(a\right)}{x}^{a-1}{e}^{-x/b}$, where a and b are the shape and scale parameters. This function gives the maximum likelihood estimators of a and b for the gamma distribution, which are the solutions of the simultaneous equations

$$log\widehat{a}-\Psi \left(\widehat{a}\right)=log\left(\overline{X}/{(\prod _{i=1}^n{x}_i)}^{1/n}\right),$$

$$\widehat{b}=\overline{X}/\widehat{a},$$

where $\overline{X}$ is the sample mean of the data X and $\Psi$ is the digamma function given by

$$\Psi \left(x\right)={\Gamma }^{\prime }\left(x\right)/\Gamma \left(x\right).$$

The function fitdist(:,gamma) estimates the shape and scale parameters with the $95\%$ confidence interval. Using this statistical technique, we estimated the recovery and death distributions for the model (3) for the data used in [31] and given by the formulas

$$r\left(t\right)=\frac{p_0}{b_1^{a_1}\Gamma \left(a_1\right)}{t}^{a_1-1}{e}^{-\frac{t}{b_1}},d\left(t\right)=\frac{(1-p_0)}{b_2^{a_2}\Gamma \left(a_2\right)}{t}^{a_2-1}{e}^{-\frac{t}{b_2}},$$

(Figure A1 in Appendix C), where the estimated parameter values are given in

Table 1. Estimated parameter values (gamma distributions).

Parameters	Estimated Value	95% Confidence Interval
$a_1$	8.06275	[6.15314, 10.565]
$b_1$	2.2140	[1.67523, 2.92623]
$a_2$	6.00014	[3.69566, 9.74161]
$b_2$	2.19887	[1.32639, 3.64526]

. Here, $p_0$ is the survival probability, and its estimated value, from the data, is $p_0=0.97$ [34].

The estimated average time to recovery is $17.85$ days and to death is $13.19$ days. The value of survival probability $p_0$ is $0.97$ [34], 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 [34,35]. Thus, we can take average disease duration $\tau$ as $17.7$ days. The corresponding value in the SIR model is $r_0+d_0\approx 1/17.7$ days${}^{-1}$.

It is essential to mention here that, instead of a gamma distribution, we can use the Erlang distribution $E_k=\frac{{\lambda }^k{x}^{k-1}{e}^{-\lambda x}}{(k-1)!},x,\lambda \ge 0$ and $k\in \mathbb{N}.$ Since, for lower values of k, e.g., $E_1$, $E_2$ distributions, with a reasonable choice of $\lambda$, give a significant probability of recovery or death at the beginning of infection, we should look for an Erlang distribution with higher values of k.

We found that the Erlang-8 distribution

$$E_8=\frac{{\lambda }^8{x}^7{e}^{-\lambda x}}{\left(7\right)!},$$

with $\lambda =0.4545$, and the Erlang-6 distribution

$$E_6=\frac{{\lambda }^6{x}^5{e}^{-\lambda x}}{\left(5\right)!},$$

with $\lambda =0.4548$, closely match with the gamma distributions for recovery and death rates respectively, as estimated above. The estimated values of the shape parameters $a_1$ and $a_2$ associated with the gamma distributions can be approximated as $a_1\approx 8$ and $a_2\approx 6$ for the recovery and death rate distributions, respectively. Thus, one can use the $E_8$ and $E_6$ distributions instead of the gamma distributions for recovery and death rates. However, the result will be similar in both cases, and we continue our numerical simulation with the gamma distributions.

Though the parameters of the three models correspond to each other, the distributed model (3) and the delay model (9) both give far different dynamics as compared to the SIR model (4), whereas the distributed model (3) and the delay model (9) give reasonably close dynamics to each other (Figure 4).

The comparison of the final size of epidemic $S_f$, maximal number of infected $I_m$, and the time to the maximal number of infected $t_m$ of the distributed model (3) with the gamma distribution, the delay model (9), and the SIR model (4) is shown in Figure 5 for different values of the transmission rate $\beta$. As before, the maximal numbers of infected individuals $I_m$ in the models (3) and (9) are sufficiently close, but they are much higher than for the SIR model (4). The times to maximum infected $t_m$ in the models (3) and (9) are reasonably close, but less than for the SIR model (4). The final size of the epidemic $S_f$ is more or less the same for all three models. Similar properties are observed if the gamma distribution is replaced by the Erlang distribution with the corresponding parameters.

This result indicates the relevance of the delay model. If we do not have sufficient individual-level data to estimate the recovery and death distributions, but we have an approximate value of disease duration, then we can describe the epidemic progression in a sufficiently precise way.

5. Determination of Disease Duration from Data

In this section, we determine the disease duration $\tau$ from the epidemiological data for the daily number of infected $J\left(t\right)$ and the total number of infected individuals $I\left(t\right)$, using the equation

$$\frac{dI\left(t\right)}{dt}=J\left(t\right)-J(t-\tau ).$$

Let $I\left(t\right)$ have the maximum at $t=t_m$. Set $I\left(t_m\right)=I_m$. Then, $J\left(t_m\right)=J(t_m-\tau )$, i.e., the daily number of infected is the same at two different time points $t=t_m$ and $t=t_m-\tau$. From the real data of the infected individuals $I\left(t\right)$, we can find the day on which the daily number of active cases is maximal, and it determines $t_m$. From the data of daily reported cases $J\left(t\right)$, we can observe that $J\left(t\right)$ crosses its maximum at some time before $t_m$. Now, we have to find the value of $J\left(t\right)$ such that $J\left(t_m\right)$ will be equal to $J(t_m-\tau )$, which in turn determines the disease duration $\tau$. Hence, considering the delay model, using the real data of daily new cases $J\left(t\right)$ and active cases $I\left(t\right)$ around a peak, we can find the disease duration $\tau$.

We illustrate this method using the data of $J\left(t\right)$ and $I\left(t\right)$ taken from [35] for COVID-19 in Italy. We collected the daily new reported data $J\left(t\right)$ and active case data $I\left(t\right)$ for Italy from 21 February 2020 to 31 May 2021 (which capture the first three peaks in Italy) and from 10 November 2021 to 28 February 2022 (which capture the peak due to Omicron in Italy). To have smoother data, we used the 7-day moving average, the data on the j-th day replaced by the average data from the $(j-3)$th day to the $(j+3)$th day. As the concerned method is focused on the peaks, the error at the beginning and end of the time interval is not essential. In Italy, during the first peak (in April 2020), the peak of $I\left(t\right)$ is attained at $t_m=51$ (Figure 6a) and the peak of $J\left(t\right)$ is attained before on $t=51$, which is less than $t_m$ (Figure 6b). First, we find that $J(t_m=51)=4.17\times {10}^3$ and then find $J\left(32\right)=4.15\times {10}^3\approx J(t_m=51)$. This implies $J(t_m-\tau )=J\left(32\right)$, and consequently, we can calculate $\tau =19$ as the disease duration during the first peak. Similarly, during the second peak (in November 2020) and third peak (in March 2021), we estimated the disease duration as $\tau =20$ days and $\tau =14$ days, respectively. The peak of epidemic progression due to Omicron in January 2022 is shown in Figure 6c,d, and we estimated the value of $\tau$ as 11 days. Similarly, the value of $\tau$ is estimated for some other countries (

Table 2. Estimated value of $\tau$ for different countries during different outbreaks of the COVID-19 epidemic. The months indicated in the table correspond to the time when the corresponding peak appeared.

Country	Estimated Value	Estimated Value	Estimated Value	Estimated Value
	of $\mathit{\tau }$ (in Days)	of $\mathit{\tau }$ (in Days)	of $\mathit{\tau }$ (in Days)	of $\mathit{\tau }$ (in Days)
	during Peak 1	during Peak 2	during Peak 3	during Peak 4
Italy	19 (April 2020)	20 (November 2020)	14 (March 2021)	11 (January 2022)
Russia	25 (May 2020)	24 (January 2021)	26 (November 2021)	9 (February 2022)
China	16 (February 2020)	-	-	-
Romania	16 (November 2020)	14 (March 2021)	18 (October 2021)	12 (February 2022)
Sweden	20 (July 2020)	20 (December 2020)	19 (April 2021)	14 (February 2022)
Iran	14 (December 2020)	24 (May 2021)	28 (August 2021)	10 (February 2022)

).

Remark 1.

It is important to note that in the case of the combination of strains, the estimated value of τ corresponds to the weighted average of the strain-specific disease duration. Furthermore, note that this method depends on how the daily cases of infection were reported. However, in many cases, there is a dominant variant during epidemic outbreaks. Thus, the estimated value of τ can be considered as the disease duration corresponding to the dominant variant during a specific epidemic wave.

Remark 2.

Note that the delay model (9) does not take into account the difference in duration to recovery and the duration to death. If we assume that the death cases are relatively rare (e.g., approximately $\le 2\%$ in the case of COVID-19), then this difference may not be very essential.

Remark 3.

Note that this method of the estimation of τ remains the same for time-varying $\beta =\beta \left(t\right)$. Thus, the method discussed above does not depend on whether β is a constant or β varies with respect to time. Since the time-varying $\beta \left(t\right)$ is capable of incorporating the effect of all possible infected compartments (i.e., exposed $E\left(t\right)$, asymptomatic $I_A\left(t\right)$, symptomatic $I_S\left(t\right)$, hospitalized $H\left(t\right)$, etc.) on disease transmission, we can use this method to obtain an estimation of τ for any infectious disease with available data.

6. Model Validation with Epidemiological Data

In order to validate the delay model (9) and the method of estimation of $\tau$, we compared the results with the epidemiological data. We used the estimated value of $\tau$ obtained by the method described above during different peaks of COVID-19.

In Italy, the estimated value of $\tau$ is 19 days, 20 days, 14 days, and 11 days during the first peak (April 2020), the second peak (November 2020), the third peak (March 2021), and the fourth peak (January 2022), respectively. We assumed that the disease duration $\tau$ remains $\tau =19$ days from 21 February 2020 to 15 August 2020; $\tau =20$ days from 16 August 2020 to 19 February 2021; $\tau =14$ days from 20 February 2021 to 21 May 2021; and 11 days from 10 November 2021 to 28 February 2022 (during Omicron in Italy).

Once these parameter values were determined, we took the number $J\left(t\right)$ of daily infected individuals from the epidemiological data [35] and found the sum of daily recoveries and deaths from the expression

$$\Sigma \left(t\right)=J(t-\tau ).$$

(24)

These results were compared with the sum of recoveries and deaths in the data. Figure 7 shows the result of such a comparison for Italy from 21 February 2020 to 21 May 2021 and from 10 November 2021 to 28 February 2022, with the data from [35] (7-day moving average). Recoveries and deaths can also be determined as a proportion of active cases $\sigma \left(t\right)=(r_0+d_0)I\left(t\right)$, as is done in the SIR model. Here, $I\left(t\right)$ is taken from the data and $r_0+d_0=1/\tau$, and we observed that the SIR model overestimates the sum of recovered and dead. Thus, the delay model (9) gives a good description of the recovery and death compared with the epidemiological data, while the SIR model overestimates them.

7. Discussion

We proposed a delay model under the assumption that the infected individuals recover or die exactly after an average disease duration $\tau$. Generally, in the case of COVID-19, we observed that the recovery and death distributions follow unimodal or bimodal gamma distributions [34,36]. Estimating these gamma distributions requires individual-level data with the date of onset of disease and the date of recovery or death, which may be very difficult to gather for every country or province. Instead, the only information about the disease duration can help us to obtain sufficiently good results using the delay model. Furthermore, we developed a method to estimate the disease duration from the epidemiological data. It is important to mention here that one can use the Erlang distributions, instead of the gamma distributions, for a compartmental epidemic model with multi-phase disease transition [37].

Let us note that we consider only symptomatic individuals in the model. The influence of asymptomatic individuals is widely discussed in the COVID-19 literature. According to some estimates, they can constitute between $25\%$ and $50\%$ of the total number of cases [38,39]. On the other hand, the infectivity of asymptomatic individuals is much lower than the infectivity of symptomatic individuals because infectivity is proportional to the viral load in the upper respiratory tract [40] and symptoms correlate with viral load. Hence, in the first approximation, we can consider only symptomatic individuals. Further studies are needed to take into account asymptomatic individuals more precisely.

We noticed that during different peaks of COVID-19, the estimated value of $\tau$ was different. This difference might be due to different strains or the change of proportion of different infected compartments (such as asymptomatic, hospitalized) that we counted in the same compartment.

The presence of exposed individuals can be taken into account by means of time-dependent infectivity rate $\beta (t-\eta )$. For the individuals infected at time $\eta$, their infectivity at time t depends on the difference $t-\eta$. The function $\beta (t-\eta )$ is small if the difference $t-\eta$ is small, which is the case of exposed individuals. This case was studied in [32]. We did not consider the time-dependent infectivity rate in this work since its main objective is to compare the model with distributed recovery and death rates with the delay model.

We compared the final size of the epidemic and maximal number of infected obtained in the Formulas (15), (21), and (23), respectively. Note that these formulae depend on $\tau$ and $\beta$. Since the outbreak due to Omicron is the only one when the social distancing measures such as lockdowns or isolation were not strictly imposed, we can assume that the transmission rate $\beta$ is approximately constant. We took the data of the Omicron outbreak in Italy from [35] from 10 November 2021 to 10 December 2021, and we fit the delay model to these 40 days of data and estimated the disease transmission rate $\beta$ as $0.118$. We took the value of disease duration $\tau =11$ days, which was obtained by using the method discussed in Section 5. Then, using the formula (15), we calculated the final size of the epidemic as $3.387\times {10}^7$. Using the formulae (21) and (23), we calculated the maximal number of infected as $I_m=3.4139\times {10}^6$, whereas the maximal number of infected was $I_m=2.7317\times {10}^6$ in the data. Thus, the formulae (21) and (23) give a reasonable estimate of the maximal number of infected. Similarly, we obtained an accurate estimate of the maximal number of infected for some other countries (not shown).

Let us finally note that the delay model presented in this work is simple and generic. It describes epidemic progression with two parameters $\beta$ and $\tau$, which can be easily estimated from the data. Our next goal will be to apply the proposed modeling approach to multi-compartment models consisting of different groups of susceptible and/or infected and to immuno-epidemic models with time-varying recovery and death rates [32]. It is also interesting to check the applicability of the proposed model to other transmissible diseases.