An epidemic model with time delays determined by the infectivity and disease durations

1.   Introduction

Mathematical modeling of infectious diseases is a valuable tool that has been used to understand the dynamics by which infections spread, to predict the future progress of an outbreak and to evaluate strategies to control an epidemic ^[1]. This tool attracts much attention due to successive epidemics of viral infections, such as HIV, emerging in the 1980s and still continuing ^[2,3], SARS epidemic in 2002–2003 ^[4,5], H5N1 influenza in 2005 ^[6,7] and H1N1 in 2009 ^[8,9], Ebola in 2014 ^[10,11].

Mathematical epidemiology was put up to a new level by the SIR model which appeared in the works published by Kermack and McKendrick ^[12,13], stimulated by the Spanish influenza epidemic in 1918–1919. Many epidemic models have been evolved from that innovative model, such as multi-compartment models ^[14,15], models with a nonlinear disease transmission rate ^[16,17], multi-patch models ^[18,19], multi-group models incorporating the effect of the heterogeneity of the population ^[20], and epidemic models with vaccination and other control measures ^[21,22]. Random movement of individuals in the population is considered in spatiotemporal models in order to describe spatial distributions of susceptible and infected individuals ^[23,24]. A more detailed literature review can be found in the monographs ^[25,26] and review articles ^[27,28].

The compartmental epidemiological models, like the conventional SEIRD model (Susceptible, Exposed, Infectious, Recovered, Dead), are based on the assumptions that newly exposed individuals at time $t$ depend on the rate proportional to the product of the numbers of susceptible individuals $S(t)$ and infectious 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 because it does not take into account infectivity period and disease duration, which can lead to a large error. Assuming average disease duration $\tau$, the recovery and death rates at time $t$ are determined by the number of newly infected individuals at time $t-\tau$, which can be very different from the number of infected individuals at time $t$. To take into account this factor, we consider the recovery and death functions as functions of the time since infection. However, since distributed recovery and death rates are not easily available, we develop a simpler delay model which does not require precise immunological data.

Let us recall that exposed individuals are those for whom viral infection develops in the organism, but they do not yet transmit infection to other individuals. Infectivity of respiratory viral infections is determined by the level of viral load in the upper respiratory tract. Exposed individuals become infectious, that is, they transmit the disease to uninfected individuals, when the viral load becomes sufficiently high. Transition from exposed to infectious compartments for SARS-CoV-2 is from 2 to 5 days, depending on the virus variant. We will call this time period infectivity delay. On the other hand, appearance of symptoms also occurs several days post-infection, and it corresponds to the incubation period. Emergence and severity of symptoms is determined by several factors, including virus virulence and individual immune response. Strictly speaking, incubation period is different from the infectivity delay, though their duration can be quite close.

In this work, we study the influence of the infectivity and disease duration on the epidemic progression. We propose a model with two time delays. At every moment $t$, a number of newly exposed people will leave the susceptible category. After the first time delay, this group of people will become infectious. Then, after the disease duration, the members of this group will either recover or die.

Let us note that the model considered in this work is different from the model considered in ^[30] since time delay there was determined by the incubation period and quarantine measures. The model proposed in this paper depends on the infectivity period rather than incubation period, and disease duration, without imposed quarantine. Furthermore, in addition to ^[31], this model takes into account exposed individuals and two time delays (infectivity and disease durations).

It is noteworthy to mention the impact of time delays on Turing instability and indicate the role of diffusion in epidemic progression ^[32]. Stability and bifurcation analysis of equilibrium points can be useful to reveal the epidemic state. Moreover, the threshold value of time delay provides a way to control the periodic behavior important for the study of disease progression.

The contents of the paper are as follows. First, we introduce the model with distributed recovery and death rates and study the existence, uniqueness and positiveness of its solutions. Next, we show how we derive the delay model from the distributed model, and we study the existence, uniqueness and positiveness of its solutions. After that, we calculate some characteristics of epidemic progression in the delay model. Next, we apply this model to the investigation of Omicron variant of the SARS-CoV-2 infection, and compare the delay model with the conventional ODE model. After that, we present an analytical comparison between the growth rates of ODE and delay models. Finally, we present conclusions and further perspectives.

2.   Model with distributed parameters

Let us start the model derivation by the introduction of the class of newly exposed individuals in the model with distributed recovery and death rates previously developed in ^[31]. This approach is appropriate to evaluate daily recovery and death rates. We will study the properties of this model, and we will reduce it to the delay model in the next section. This model can also be reduced to the SEIRD model under certain assumptions on the recovery and death rates.

2.1.   Model derivation

The size of newly exposed individuals $J(t)$ is determined by the rate of decrease of the size of susceptible individuals

Assuming that

is constant, where $E(t)$ is the total size of exposed compartment, $I(t)$ is the total size of infectious compartment at time $t$ and $R(t)$ and $D(t)$ denote, respectively, recovered and dead. Let $\tau_0$ be the infectivity period. Then $J(t-\tau_0)$ is the size of newly infectious at time $t$. Thus,

Now, differentiating the Eq (2.2), we get:

Therefore,

Let $\rho(\eta)$ and $\mu(\eta)$ be the recovery and death rates depending on the time since being infectious $\eta$. Then $\rho(t-\eta)$ and $\mu(t-\eta)$ are the recovery and death rates at time $t$ of the individuals became infectious at time $\eta$. Thus, the rate of change of the size of infected individuals who will recover at time t is given by the expression:

and the rate of change of the size of infected individuals who will die at time t:

Here $J(\eta-\tau_0)$ is the size of individuals who became exposed at time $\eta-\tau_0$ and infectious at time $\eta$. Thus, we obtain the following integro-differential system of equations:

2.2.   Existence and uniqueness of solution

We will prove the existence and uniqueness of solution of system (2.8) for $t \in [0, \mathcal{T}_f]$ where $\mathcal{T}_f \in (0, \infty)$, with the initial conditions

In what follows we assume that the recovery and death rates satisfy the following inequality

for any $\eta$ and $t_0 > \eta$, and $\rho(\eta) = \mu(\eta) = 0$ for $\eta < 0$. This condition means that the total size of recovered and dead individuals among those infected at time $\eta$ remains less than $J(\eta)$. We assume that the coefficients are continuous and positive. Note that if $J(t)$ is uniquely determined, then the Eqs (2.8b), (2.8d) and (2.8e) have unique solutions. Hence, it is sufficient to prove the existence and uniqueness of solution for the equations (2.8a) and (2.8c).

Before proving the existence and uniqueness of solution, we will verify that the solutions of system (2.8) with initial conditions (2.9) are positive and bounded.

Lemma 1. If condition (2.10) is satisfied, then any solution of system (2.8) with initial condition (2.9) satisfies the inequality $0 \leq A \leq S_0 + I_0$, where $A \in \{ S(t), E(t), I(t), R(t), D(t)\}$.

Proof. From (2.8a) we observe that if for some $t_* > 0$, $S(t_*) = 0$, then $\frac{dS(t)}{dt}|_{t = t_*} = 0$. This implies that $S(t) \geq 0$ for $t > 0$. From (2.8d) and (2.8e), we conclude that $R(t)$ and $D(t)$ also remain positive for all t.

Using the initial conditions (2.9), we have $J(t) = 0$ for $-\tau \leq t < 0$ and $J(0) = \frac{\beta S_0 I_0}{N} > 0$. Suppose $t_0 > 0$ be such that $J(t) \geq 0$ for $0 \leq t < t_0$. Next, we have

Integrating (2.8d) and (2.8e) with respect to $t$ from $0$ to $t_0$ and taking into account the initial conditions (2.9), we get

By changing the order of integration and taking into account inequality (2.10), we find:

Together with (2.11), this gives $I(t_0) \geq 0$. This implies $J(t_0) = \frac{\beta S(t_0) I(t_0)}{N} \geq 0$. This proves that $I(t)$ and $J(t)$ remains non-negative for all $t \geq 0$. On the other hand, integrating (2.8b) from $0$ to $t$, we get

Furthermore, $S(t) + E(t) + I(t) + R(t) + D(t) = S_0 + I_0$. Thus, any solution of system (2.8) lies between 0 and $S_0 + I_0$. □

We now proceed to the proof of the existence and uniqueness theorem.

Theorem 2. There exists a unique solution $(S(t), I(t))$ of system (2.8a) and (2.8c) in the domain $\Omega^2$, where $\Omega$ is defined by

To prove this theorem, we need a mathematical setup of complete metric space, which is defined in the following lemma.

Lemma 3. $(\Omega, d)$ is a complete metric space with respect to the metric $d(T_1, T_2)$ defined by the equality

where $\gamma \geq 0$ is some constant.

Proof. First, we prove that $\Omega$ is a complete metric space with respect to the supremum metric given by the equality

Consider a Cauchy sequence $\big\{T_n(t)\big\}$ in $\Omega$. Then for any $\epsilon > 0$, there exists $N_0 \in \mathbb{N}$ such that

Therefore, for all $t \in [-\tau_0, \mathcal{T}_f]$, $\big\{T_n(t)\big\}$ is a Cauchy sequence in $\mathbb{R}$, and hence converges to a real number denoted by $T(t)$. Choose any $t\in [-\tau_0, \mathcal{T}_f]$. Hence, there exists $P_t \in \mathbb{N}$ such that if $p > P_t$ then $|T_p(t)-T(t)| < \epsilon /2.$

Furthermore, since $\big\{T_n\big\}$ is a Cauchy sequence in $(\Omega, d_{sup})$, there exists $N_1$ such that

Next, choose $p > max\{ N_1, P_t\}$. Then for all $n \geq N_1$

Taking supremum over $[-\tau_0, \mathcal{T}_f]$ in both sides of the above inequality, we get

It remains to show that $T\in \Omega$. It is clear that for all $n \in \mathbb{N}$, $0 \leq T_n(t) \leq S_0+I_0$, for all $t \in [0, \mathcal{T}_f]$. Taking limit as $n \to \infty$, we get $0 \leq T(t) \leq S_0+I_0$, for all $t \in [0, \mathcal{T}_f]$. Similarly, for all $n \in \mathbb{N}$, $T_n(t) = 0$, $\forall t \in [-\tau_0, 0)$ implies $T(t) = 0$, $\forall t \in [-\tau_0, 0)$.

Take any $t_0 \in [0, \mathcal{T}_f]$. Then

which proves that $T$ is continuous at $t_0$. Thus, $T \in \Omega$, and hence $(\Omega, d_{sup})$ is a complete metric space.

Next, we have the following relation between the two metrics $d$ and $d_{sup}$ on $\Omega$:

which implies that $d$ and $d_{sup}$ are equivalent metrics. This proves that $(\Omega, d)$ is a complete metric space.

□

We now proceed to prove the existence and uniqueness of solution of system (2.8a) and (2.8c) in the metric space $(\Omega, d)$. For any given function $T(t) \in \Omega$, equation

with initial condition $S(0) = S_0 > 0$ has a unique solution

Note that subscript $T$ is used to denote the unique solution of Eq (2.12) for a given function $T(t) \in \Omega$. Set $J_T(t) = \frac{\beta}{N} S_T(t)T(t)$. Under the assumption that $\rho(t) = 0 = \mu(t)$ for $t < 0$, equation

with $I(0) = I_0 > 0$ also has a unique solution which can be written in the form

where

Let us now consider the map $L:(\Omega, d)\rightarrow (\Omega, d)$ defined by the equality

where $G(\eta, T)$ is given in (2.16). Before proceeding further, we verify that $L$ maps $(\Omega, d)$ into itself.

Lemma 4. The map $L:(\Omega, d)\rightarrow (\Omega, d)$ given by equality (2.17) is well defined.

Proof. From (2.13) we obtain

Substituting this relation into (2.16), we can write

Next,

Changing the order of integration in the right-hand side, we get

Note that $\frac{dS_T(\xi)}{d\xi} \leq 0$ and following condition (2.10), we conclude that

This implies $L(T(t)) = I_0+\int_{0}^{t}G(\eta, T)d\eta$ lies between $0$ and $S_0+I_0$.

Let us also note that if $T_1(t), T_2(t) \in \Omega$ and $T_1(t) = T_2(t)$, then $S_{T_1}(t) = S_{T_2}(t)$, and consequently $G(\eta, T_1) = G(\eta, T_2)$. Hence the map $L$ is well defined. □

Next, we prove that the map $L:(\Omega, d)\rightarrow (\Omega, d)$ defined in (2.17) is a contraction.

Lemma 5. The map $L:(\Omega, d)\rightarrow (\Omega, d)$ defined in (2.17) is a contraction map.

Proof. For any two functions $T_1(t), T_2(t)\in \Omega$, we obtain the inequality

Then we have the following estimate:

Therefore,

Using the inequalities $|e^{-x}-e^{-y}|\leq|x-y|$ and $|e^{-x}|\leq1$, for any $x, y\geq0$, we get

Since $T_j(t)\leq S_0+I_0$, $\forall t \in [0, \mathcal{T}_f]$, then $T_j(t) \leq M, \; j = 1, 2$, where $M = S_0+I_0$. Next, $|T_1(\eta)-T_2(\eta)| \leq e^{\gamma \eta} d(T_1, T_2)$. Using this inequality and condition (2.10), we can write

Thus,

This implies the estimate

Since $\int_{0}^{t}e^{\gamma(\eta-\tau_0)}d\eta \leq \int_{0}^{t}e^{\gamma \eta}d\eta$, we get

Therefore,

Taking supremum of both sides, we get

We choose the value of $\gamma > 0$ large enough such that $\frac{\beta S_0}{N}\bigg(\frac{1}{\gamma}+\frac{N+\beta M}{N \gamma^2}+\frac{\beta M}{N \gamma^3} \bigg) < 1$. Consequently, $L:(\Omega, d)\rightarrow (\Omega, d)$ is a contraction map on the complete metric space $(\Omega, d)$. □

To finish the proof of the existence of solution, we use the following theorem ^[33].

Theorem 6. Let $(X, d)$ be a complete metric space and let $T : X\rightarrow X$ be a contraction mapping on X. Then T has a unique fixed point $x \in X$ (such that $T(x) = x$).

It follows from this theorem that the map $L$ has a unique fixed point. Thus, there exists a unique function $T_u \in \Omega$ satisfying the equality $T_u(t) = I_0+\int_{0}^{t}G(\eta, T_u)d\eta$, where $G(\eta, T_u)$ is given in (2.16). Besides, we note that $G(\eta, T)$ is a continuous function on $[0, \mathcal{T}_f]$. Hence, the derivative $\frac{dT_u(t)}{dt}$ exists. This completes the proof of the existence and uniqueness of solution of system (2.8).

2.3.   Reduction to the SEIRD Model

Let us assume the uniform distribution of the recovery and death rates:

where $\tau > 0$ is the disease duration, $\rho_0$ and $\mu_0$ are small enough constants, and consider the initial conditions: $S(0) = N$, $E(0) = 0$, $I(0) = I_0 > 0$, where $I_0$ is sufficiently small compared to $N$, $R(0) = 0$, $D(0) = 0$.

From (2.8b), we get

Therefore, (2.8b) can be written as follows:

where $\lambda = 1/\tau_0$.

Substituting functions (2.18) and (2.19) into (2.8d) and (2.8e) respectively, we get

Integrating these equalities from $t-\tau$ to $t$, we obtain

Therefore, from (2.5)

Hence,

or

Since $(\rho_0 + \mu_0)$ is assumed to be small enough, we neglect the term involving $(\rho_0+\mu_0)^2$. Hence, we obtain

In this case, system (2.8) is reduced to the conventional SEIRD model

3.   Delay model

3.1.   Reduction to the delay model

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

where $\rho_1 + \mu_1 = 1$, and $\delta$ is the Dirac delta-function. Then

and

Note that $J(t)$ is the number of newly exposed individuals appearing at time $t$. If we assume that the first infected case was reported at time t = 0, then we can set $J(t) = 0$ for all $t < 0$.

Integrating the Eqs (3.1) and (3.2) from 0 to $t$ and assuming that $R(0) = D(0) = 0$, we get:

and

Therefore,

Since $\rho_1 + \mu_1 = 1,$ then,

Finally, we obtain

Let $y = x-\tau_0$. Then

We obtain

Then the delay model becomes as follows:

Here $J(t)$ is the size of newly exposed at time $t$, $J(t-\tau_0)$ is the size of newly infectious at time $t$, $J(t-\tau_0 -\tau)$ is the size of newly recovered and dead at time $t$. The term in the right-hand side of Eq (3.6a) describes the decrease of the size of susceptible due to infection, with the same term with sign plus in the next equation. The second term in the right-hand side of Eq (3.6b) corresponds to the decrease of the size of exposed due to their infectiousness, the second term in the right-hand side of Eq (3.6c) corresponds to the decrease of the size of infectious due to their recovery or death. The term in the right-hand side of (3.6d) represents the size of newly recovered individuals, and the term in the right-hand side of (3.6e) represents the size of newly dead individuals. System (3.6) is completed by the initial conditions

3.2.   Epidemic characteristics for the delay model

3.2.1.   Basic reproduction number

From (3.4), we get:

Then

At the beginning of epidemic, we have $S(t-\tau_0) = S(t-\tau_0-\tau) = S_0$. Then

Substituting $I(t) = I_0 e^{xt}$, we get

Hence,

This equation has a positive solution if and only if $\Re_0 > 1$, where $\Re_0 = \tau \beta \frac{S_0}{N}$ is the basic reproduction number.

3.2.2.   Final size of epidemic

Let us determine the final size of the susceptible class, $S_f = \lim_{x\to\infty} S(x)$. Integrating (3.6a) from 0 to $\infty$, we get

By changing the order of integration,

where $J(t) = 0$ if $t < 0$. Therefore,

Since $J(x) = -S'(x)$, then

where $w = \frac{S_f}{S_0}.$ This equation is the same as for model (2.20), but the basic reproduction number is different.

Integrating (3.6d) and (3.6e) and taking the limits as $t\rightarrow \infty$, we obtain the final size of recovered and dead populations:

3.3.   Growth rates for SEIRD and delay model

Linearizing the delay model about the disease free equilibrium, we obtain the equation for the principal eigenvalue $\nu$ which determines infection growth rate:

A similar equation for the SEIRD model is as follows:

where incubation rate $\lambda = \frac{1}{\tau_0}$ and recovery and death rates of model (2.20) satisfy $\rho_0+\mu_0 = \frac{1}{\tau}$.

The delay model (3.6) has higher growth rate (larger $\nu$) for the same basic reproduction number as compared to the ODE basic model (2.20). A higher growth rate can lead to higher peak of epidemic achieved much earlier. Thus, to describe the epidemic progression more precisely, basic reproduction number is not sufficient, rather the estimate of the growth rate can give better description of the epidemic progression. Note that the curves of growth rate are above the $\beta-$axis for $\beta > \frac{1}{\tau}$ which corresponds to $\Re_0 > 1$.

4.   Numerical simulation and model comparison

In this section we present some numerical simulations to validate the proposed model and to compare its results with real data and with the SEIRD model. We will compare modelling results with reported active cases for the Omicron variant of the SARS-CoV-2 infection. The Omicron variant was reported to the World Health Organization (WHO) by the Network for Genomics Surveillance in South Africa on 24 November 2021 ^[34,35]. It was first detected in Botswana and has spread to become the predominant variant in circulation around the world.

The Omicron variant has a shorter incubation period, compared to other variants, from 1 to 4 days ^[36]. While BA.5, like previous Omicron subvariants, seems to spread more easily than other COVID-19 variants, the symptoms have generally been milder and have a shorter duration of six to seven days ^[37]. These data allow us to estimate the values of time delay, namely, infectivity period and disease duration. Furthermore, recovery and death rates are also evaluated from the epidemiological data ^[38]. The value of $\beta$ is chosen to fit the data. The values of parameters can differ for different countries.

The parameters of the SEIRD model (2.20) can be obtained from the corresponding parameters of the delay model (3.6):

The initial conditions for the delay model are as follows:

where $\theta = \tau+\tau_0$. The values $I^-$ and $I_0$ are determined from the epidemiological data, the former as an average value for this time period, the latter is the number of daily new cases in the beginning of simulations. Though the data are not usually precise, and the relative error can be quite large in the beginning of outbreak, it appears that the results of the simulations and the best-fit value of $\beta$ are not very sensitive to the initial conditions.

The results of numerical simulations for different countries are shown in Figure 2. The delay model gives a good approximation of the data, while the SEIRD model is less accurate. We can notice that the maximum number of infected individuals in the delay model (3.6) is much higher than for the conventional ODE model (2.20). Also, the time to maximum infected in the delay model is less than for the conventional ODE model. As it is follows from the derivation of the SEIRD model from the distributed delay model, it overestimates the recovery and death rates witch leads to an underestimation in the active cases. The variation of $\beta$ (Figure 3) does not essentially improve the description of the data by the SEIRD model.

5.   Discussion

In this work we develop an epidemiological model with distributed recovery and death rates and with exposed (infected but not infectious) individuals which were not considered in the previous study ^[29,39]. This distributed delay model represents a system of integro-differential equations. It is an appropriate tool to study epidemic progression which allows an accurate desciption of its dynamics. A disadvantage of this model is that it is relatively complex and that it requires the knowledge of distributed recovery and death rates which may not be available in the literature.

In order to compensate this disadvantage of the distributed model, we derive from it two simplified models corresponding to different limiting cases. In the first one, where it is assumed that recovery and death rates are uniformly distributed, we obtain conventional compartmental SEIRD model. In the second case, where these distributions are delta-functions, we obtain a new delay model, which was not previously considered in the epidemiological literature. Since distributed recovery and death rates are described by gamma-distributions ^[29,39], the approximation by the delta-function can be more appropriate than by a uniform distribution which overestimates recoveries and deaths during first several days post-infection.

The point-wise delay model is quite simple, it has a clear biological meaning, and it is determined by two main parameters: time delay before infected individual becomes infectious and disease duration. Both of them can be easily determined from the clinical data for each particular viral infection (or virus variant). Two other parameters, recovery and death rates are also known. The only unknown parameter, disease transmission rate $\beta$, is determined by the comparison with the data on new daily infections. It appears that this model gives a good description of the COVID-19 epidemic progression in different countries.

Since the SEIRD model is obtained under the assumption of uniform distribution of recovery and death rates, then it overestimates recoveries and deaths and underestimates the number of infectious individuals. As a result, epidemic progression is slower than in the delay model and in the data. We can observe that the maximum number of infected individuals in the delay model is much higher than for the conventional ODE model. Also, the time to maximum infected in the delay model is less than for the conventional ODE model. Furthermore, the delay model (3.6) has a higher growth rate in comparison with the SEIRD model (2.20). A higher growth rate can lead to a higher peak of epidemic achieved much earlier. Therefore, to describe the epidemic progression more precisely, it is not sufficient to rely only on the basic reproduction number, but also on the estimation of the growth rate which can give a better description of the epidemic progression.

As it is indicated above, the SEIRD model overestimates recoveries and deaths during first time interval post-infection. Therefore, in order to describe the data with this model, we need to increase the disease transmission rate $\beta$ but then also to increase even more recovery and death rate. In the example shown in Figure 3(b), the average disease duration becomes 4 days for the SEIRD model, while it is 6 days for the delay model, in agreement with the clinical data. Altogether, this gives the value of the basic reproduction number $1.6$ for the SEIRD model instead of $1.32$ for the delay model.

Thus, SEIRD and delay models represent two different limiting cases of the distributed delay model. Both of them can describe the epidemic progression, but the delay model seems to be more precise from the point of view of parameter estimation.

Finally, let us note that the delay model presented in this work is generic but a little bit more complex than the delay model presented in ^[31] because it has one additional compartment. It describes epidemic progression with three parameters $\beta$, $\tau_0$ and $\tau$, which can be easily estimated from the data. This approach opens further applications to more complex multi-compartment models consisting of different groups of susceptible and/or infected and to immuno-epidemic models with time-varying recovery and death rates ^[39]. It is also interesting to check the applicability of the proposed model to other transmissible diseases.

Acknowledgments

This work is supported by the Ministry of Science and Higher Education of the Russian Federation (project number FSSF-2023-0016).

Conflict of interest

Vitaly Volpert is an Guest Editor for Mathematical Biosciences and Engineering and was not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.

Appendix: Existence and uniqueness of solution of the delay model

We will prove the existence and uniqueness of solution of system (3.6) for $t \in [0, \mathcal{T}_f]$ where $\mathcal{T}_f \in (0, \infty)$, with the initial conditions (3.7).

Note that if $J(t)$ is uniquely determined, then Eqs (3.6b), (3.6d) and (3.6e) have unique solutions. Hence, it is sufficient to prove the existence and uniqueness of solution for the equations (3.6a) and (3.6c). Before proving the existence and uniqueness of solution, we will verify that the solutions of system (3.6) with initial conditions (3.7) are positive and bounded.

Lemma 7. Any solution of system (3.6) with initial condition (3.7) satisfies the inequality $0 \leq A \leq S_0 + I_0$, where $A \in \{ S(t), E(t), I(t), R(t), D(t)\}$.

Proof. From (3.6a) we observe that if for some $t_* > 0$, $S(t_*) = 0$, then $\frac{dS(t)}{dt}|_{t = t_*} = 0$. This implies that $S(t) \geq 0$ for $t > 0$. From (3.6d)and (3.6e), we conclude that $R(t)$ and $D(t)$ also remain positive for all t. Integrating (3.6b) from $0$ to $t$ we get

Next, integrating (3.27c) from $0$ to $t$, we have

Since $J(t) = 0$ for $t < 0$, then

This gives $I(t) \geq 0$. Furthermore, $S(t) + E(t) + I(t) + R(t) + D(t) = S_0 + I_0$. Thus, any solution of system (3.6), lies between 0 and $S_0 + I_0$. □

Let us proceed to the proof of the existence and uniqueness theorem.

Theorem 8. There exists a unique solution $(S(t), I(t))$ of system (3.6a) and (3.6c) in the domain $\Omega^2$, where $\Omega$ is defined by

To prove this theorem, we use a mathematical setup of complete metric space, which is defined properly in the following lemma.

Lemma 9. $(\Omega, d)$ is a complete metric space with respect to the metric $d(T_1, T_2)$ defined by

and $\gamma \geq 0$ is a constant.

Proof. The proof is similar to that of Lemma 3. □

We now proceed to prove the existence and uniqueness of solution of system (3.6a) and (3.6c) in the metric space $(\Omega, d)$. For any given function $T(t) \in \Omega$, equation

with initial condition $S(0) = S_0 > 0$ has a unique solution

Note that subscript $T$ is used to denote the unique solution of Eq (A.1) for a given function $T(t) \in \Omega$. Let us denote $J_T(t) = \frac{\beta}{N} S_T(t)T(t)$, then equation

with $I(0) = I_0 > 0$ also has a unique solution which can be written in the form

where

Let us now consider the map $L:(\Omega, d)\rightarrow (\Omega, d)$ defined by the equality

where $G(\eta, T)$ is given in (A.5). Before proceeding further, we verify that $L$ maps $(\Omega, d)$ into itself.

Lemma 10. The map $L:(\Omega, d)\rightarrow (\Omega, d)$ defined in (A.6) is well-defined.

Proof. From (A.2) we obtain

Substituting this relation into (A.5), we can write

Next,

Thus,

This implies $L(T(t)) = I_0+\int_{0}^{t}G(\eta, T)d\eta$ lies between $0$ and $S_0+I_0$.

Let us also note that if $T_1(t), T_2(t) \in \Omega$ and $T_1(t) = T_2(t)$, then $S_{T_1}(t) = S_{T_2}(t)$, and consequently $G(\eta, T_1) = G(\eta, T_2)$. Hence, the map $L$ is well-defined. □

Next, we prove that the map $L:(\Omega, d)\rightarrow (\Omega, d)$ defined in (A.6) is a contraction.

Lemma 11. The map $L:(\Omega, d)\rightarrow (\Omega, d)$ defined in (A.6) is a contraction map.

Proof. For any two functions $T_1(t), T_2(t)\in \Omega$,

Then we have the following estimate:

Therefore,

Using the inequalities $|e^{-x}-e^{-y}|\leq|x-y|$ and $|e^{-x}|\leq1$, for any $x, y\geq0$, we get

Since $T_j(t)\leq S_0+I_0$, $\forall t \in [0, \mathcal{T}_f]$, we get $T_j(t) \leq M, \; j = 1, 2$, where $M = S_0+I_0$.

Next, $|T_1(\eta)-T_2(\eta)| \leq e^{\gamma \eta} d(T_1, T_2)$.

Using the inequality, we can write

Thus,

This implies the estimate

Since $\int_{0}^{t}e^{\gamma(\eta-\tau_0)}d\eta \leq \int_{0}^{t}e^{\gamma \eta}d\eta$, we get

This implies

Taking the supremum of both sides, we get

We choose the value of $\gamma > 0$ large enough such that $2\frac{\beta S_0}{N}\bigg(\frac{1}{\gamma}+\frac{\beta M}{N \gamma^2} \bigg) < 1$. Consequently, $L:(\Omega, d)\rightarrow (\Omega, d)$ is a contraction map on the complete metric space $(\Omega, d)$. □

Finally, we use Theorem 6, which completes the proof of the existence and uniqueness of solution of system (3.6).