Explicit solution of a Lotka-Sharpe-McKendrick system involving neutral delay differential equations using the r-Lambert W function

1.   Introduction

An outbreak of pneumonia initially took place in Wuhan, China, which turned out to be coronavirus disease 2019 (COVID-19), caused by a novel coronavirus (SARS-CoV-2). The pandemic has brought about devastating public health, economic, and political crises. The outbreak has spread to 222 countries and territories, inflicted more than 116.17 million confirmed cases, and claimed more than 2.58 million lives, as reported on March 7th, 2021 ^[1]. People residing in the United States (US) have paid a high price; the pandemic has caused more than 29.05 million confirmed cases and claimed 527,726 lives, as reported on March 11th, 2021 ^[2,3]. The emergence of COVID-19 pandemic demonstrates the importance of research on accurately understanding and modeling the transmission and spread of the COVID-19 and SARS-CoV-2 ^[4,5,6]. It is essential to understand the mathematical mechanism and dynamics of COVID-19 to predict the pandemic better and evaluate the effects of intervention strategies and offer effective policies ^[7]. This paper is motivated by the need for in-depth investigations.

The COVID-19 is a communicable disease. An exposed person may interact with and transmit the virus to a susceptible individual. An interaction between the exposed person and the susceptible individual may result in the susceptible individual becoming exposed. An exposed individual may experience an incubation period before becoming ill, i.e., infected. Thus, an exposed person is infectious but not yet infected. An infected individual may be lucky enough to recover, at least to some extent, after a period of time. Alternatively, they will die. This process can be described by a dynamical system of susceptible-exposed-infected-recovered-deceased (SEIRD) differential equations, which can connect the numbers of susceptible, exposed, infected, recovered, and deceased individuals. Figure 1 shows the SEIRD model flow diagram of the transmission of COVID-19.

Figure 1.  The SEIRD model flow diagram of the transmission of COVID-19.

DownLoad: Full-Size Img PowerPoint

In the Figure 1, $\alpha (t)$ is the contact rate between susceptible and exposed people, $\beta (t)$ is the incidence rate of the exposed people becoming infected, $\mu (t)$ is the recovery rate of the infected individuals, and $\delta (t)$ is the mortality rate of the infected individuals. In the literature, a system of susceptible-exposed-infected-recovered (SEIR) differential equations was developed to describe diseases like chickenpox and vector-borne diseases such as dengue hemorrhagic fever ^{[8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]}. These diseases have a long incubation period, during which time the exposed individual cannot yet transmit the pathogen to others. The SEIR assumes that only an infected individual is infectious and can transmit the virus to a susceptible person and an exposed person is not infectious, as with susceptible-infected-recovered (SIR) models. For the very dangerous COVID-19, an exposed individual is infectious and can transmit the virus to susceptible people rapidly. An infected person is infectious but may be sent to hospitals quickly and separated from the susceptible people, making them unlikely to transmit the virus to susceptible people. Therefore, the existing SEIR and SIR systems cannot be used to analyze the COVID-19, and we need to build new models to fill the gap.

To get an intuitive idea of these rates for the COVID-19 transmission process, we obtain population data from US Census Bureau and COVID-19 daily reports from the Johns Hopkins University website through a GitHub repository ^[4,5,6]. The data are then combined for a unified analysis. We then calculate the empirical incidence, mortality, and recovered rates. Figure 2 provides trajectories of total numbers of people tested, confirmed cases, deaths, and recovered as well as related linear regression lines for COVID-19 in the US from April 16th to June 30th, 2020.

Figure 2.  The trajectories of total numbers of people tested, confirmed cases, deaths, and recovered as well as related linear regression lines for COVID-19 in the United States.

DownLoad: Full-Size Img PowerPoint

In Figures 3 and S1, we show trajectories of empirical testing rate, incidence rate, mortality rate, and recovery rate as well as related quadratic and linear regression lines for COVID-19 based on the results of Figure 2, respectively. It can be seen from the Figures 3 and S1 that both the testing rate and recovery rate increase, the incidence rate decreases, and the mortality rate increases at the beginning and then remains relatively stable before decreasing. The features of these rate may be due to people's being inexperienced at the beginning of the pandemic but gradually becoming more experienced at dealing with it and getting better prepared.

Figure 3.  The trajectories of testing rates, incidence rates, mortality rates, and recovery rates as well as related quadratic and linear regression lines for COVID-19 in the United States.

DownLoad: Full-Size Img PowerPoint

The empirical results shown in Figures 1, 3 and S1 provide us with a rationale for building the models. One noticeable feature is that these rates are not constants but time-dependent. Figures 1, 3 and S1 are based on COVID-19 data in the US from April 16th to June 30th, 2020, in which the sample size is large, and the deterministic SEIRD system is appropriate. The results of Figures 1, 3 and S1 provide a good estimate of the incidence, mortality, and recovery rates of COVID-19. Based on the empirical results, we can consider a few cases of contact, incidence, recovery, and mortality rates to model dynamics of COVID-19.

In this article, we build a system of SEIRD ordinary differential equations to describe relationship among the numbers of susceptible, exposed, infected, recovered, and deceased individuals in a town, state or country based on the empirical rates. The organization of the paper is as follows. In section 2, we build SEIRD models in details. We perform empirical analysis for the SEIRD models in section 3. In section 4, we provide some discussion about the property of the SEIRD models and some thoughts about future research direction.

2.   The SEIRD models

Consider an independent time variable $t$ which is measured in days. At time $t$, let $s(t)$ be the number of susceptible individuals in a town, state or country; $e(t)$, the number of exposed individuals that have transmitted the virus; $i(t)$, the number of exposed individuals that are infected; $r(t)$, the number of infected individuals that have recovered; and $d(t)$, the number of infected individuals that have died. To describe the transmission process of COVID-19, consider a SEIRD system of ordinary differential equations

with initial conditions $s(0) = S_0, e(0) = E_0, i(0) = I_0, r(0) = 0$, and $d(0) = 0$. In the system (2.1), $\alpha (t)$ is contact rate, $\beta (t)$ is incidence rate, $\mu (t)$ is recovery rate, and $\delta (t)$ is mortality rate displayed in the Figure 1. Note all these rates are non-negative. The first equation in the system (2.1) is called susceptible equation, the second equation in the system (2.1) is called exposed equation, the third is called infected equation, the fourth is called recovery equation, and the fifth is called deceased equation, respectively. The system (2.1) is different from those in the literature of SEIR system as the susceptible and exposed equations are different. In the SEIR models, it is assumed that only an infected individual can transmit the virus to a susceptible person like SIR models and so the first two SEIR equations are

For COVID-19, however, an exposed individual is infectious and may interact with a susceptible individual, thus transmitting the virus. This transmission is described by the susceptible equation in system (2.1). The exposed equation in system (2.1) implies that an exposed individual may contract coronavirus and experience an incubation period before becoming ill, i.e., infected. Once contracting the virus, an exposed person can transmit the virus to a susceptible individual. The infected can either recover with a recovery rate $\mu(t)$ or die with a mortality rate $\delta (t)$.

Let us denote $N = S_0 + E_0 + I_0$. Note

It follows that $s(t) + e(t) + i(t) + r(t) + d(t) = N$ is a constant. Since $\frac {ds(t)}{dt} = - \alpha (t) s(t)e(t) \le 0$, the susceptible function $s(t)$ is a decreasing function. Similarly, the recovery function $r(t)$ and deceased function $d(t)$ are increasing functions since $\frac {dr(t)}{dt} = \mu(t) i(t) \ge 0$ and $\frac {dd(t)}{dt} = \delta(t) i(t) \ge 0$, respectively. In the following, we impose conditions on incidence rate $\beta (t),$ recovery rate $\mu (t)$, and mortality rate $\delta (t)$ based on the empirical results of the Figures 3 and S1 as follows

As shown in Figure 3, a declining incidence rate implies that the proportion of infected people is high at the beginning but decreases as time goes by. The empirical mortality rate is assumed to be a quadratic curve. It implies that the mortality rate increases initially, after which it remains stable and eventually decreases. One possible reason for this is that the virus is strongly virulent at the beginning but gradually weakens. Also, healthcare for infected people improves, and this contributes to an eventual decline in mortality. On the other hand, an increasing recovery rate is most likely due to improving healthcare and the fact that people are more experienced at handling the disease as time progresses. We assume that the recovery rate is linearly increasing from the empirical result of COVID-19 in the US, as shown in Figure 3.

For the contact rate $\alpha (t)$, we consider an exponential declining rate: $\alpha (t) \; = \; \alpha_0 \exp(-\alpha_1 t), \; \alpha_0 > \; 0, \; \alpha_1 > 0.$ The rationale for a declining contact rate is that people may contact each other as normal at the beginning stage of the pandemic but reduce contact quickly as time moves on. The exponential decline means that the contact reduces to a minimum fast.

From the susceptible equation and the exposed equation in the system (2.1), we have

Integrating the above equation gives

If $\alpha (t) = \alpha_0$ and $\beta (t) = \beta_0$ are constants, the above equation provides an exact solution

In addition, $\frac {de(t)}{ds(t)} = -1 + \frac {\beta_0}{\alpha_0 s(t)} = 0$ implies that $e(t)$ reaches its maximum when $s(t) = 1/q$. Hence, the maximal number of exposed individual is given by

3.   Numerical results

Under an assumption of constant contacting rate $\alpha (t) = \alpha_0$ and constant incidence rate $\beta (t) = \beta_0$, we obtain the exposed function $e(t)$ and the maximal number of exposed individual $e_{\mbox{max}}$ in relations (2.3) and (2.4). In general, the contact rates, incidence rates, recovery rates, and mortality rates are functions of time variable $t$. An exact solutions of SEIRD system (2.1) is not possible. Instead, a numerical approximation procedure can provide us a solution via Euler's method ^[33,34,35]. Taking a small time interval $\Delta t$ and increment $\Delta s(t) = s(t + \Delta t) - s(t)$, we may make an approximation $\frac {d s(t)}{dt} = \frac {\Delta s(t)} {\Delta t}$. Similarly, we may approximate $\frac {d i(t)}{dt} = \frac {\Delta i(t)} {\Delta t}$, where $\Delta i(t) = i(t + \Delta t) - i(t)$, etc. Then, we may obtain a system of recurrence equations as follows to approximate the numbers of susceptible, exposed, infected, recovered, and deceased

Using the recurrence system (3.1) and initial conditions $s(0) = S_0, e(0) = E_0, i(0) = I_0, r(0) = 0$, and $d(0) = 0$, we perform numerical approximations for a few cases of the rates defined in relations (2.2).

By using the rate functions provided in the Figure 3, we perform numerical calculation for approximate numbers of susceptible, exposed, infected, and deceased. The incidence, recovery, and mortality rates are from Figure 3: $\beta(t) = 0.01(21.0 -0.30t +0.0018 t^2),$ $\mu (t) = 0.01(14.65 + 0.188 t),$ and $\delta (t) = 0.01(5.25 + 0.045t -0.00067 t^2)$. The initial conditions are given by $S_0 = 10^6, E_0 = 10$, and $I_0 = 2$. We consider four exponential decline contact functions $\alpha (t) = 1.2 \times 10^{-6} \exp (\alpha_1 t)$: (a) $\alpha_1 = -0.01$, (b) $\alpha_1 = -0.025$, (c) $\alpha_1 = -0.05$, and (d) $\alpha_1 = -0.065$, respectively. Thus, the exponential declining rates of the four contact functions decrease from $-0.01$ to $-0.065$. The results are plotted in Figure 4.

Figure 4.  Approximate numbers of susceptible, exposed, infected, and deceased using exponential decline contact function $\alpha (t) = 1.2 \times 10^{-6} \exp (\alpha_1 t)$, (a) $\alpha_1 = -0.01$, (b) $\alpha_1 = -0.025$, (c) $\alpha_1 = -0.05$, and (d) $\alpha_1 = -0.065$, when other rates are fixed: $\beta(t) = 0.01(21.0 -0.30t +0.0018 t^2), \; \mu (t) = 0.01(14.65 + 0.188 t),$ and $\delta (t) = 0.01(5.25 + 0.045t -0.00067 t^2)$. The initial conditions are given by $S_0 = 10^6, E_0 = 10$, and $I_0 = 2$.

DownLoad: Full-Size Img PowerPoint

It can be seen from Figure 4(a), (b) that the number of susceptible individuals rapidly decreases to zero when $\alpha_1 = -0.01$ and $\alpha_1 = -0.025$. If $\alpha_1 = -0.05$, the number of susceptible individuals decreases and stabilizes around 100,000 when the time t is about 30 days (Figure 4 (c)); if $\alpha_1 = -0.065$, the number of susceptible individuals decreases and stabilizes around 500,000 when the time t is about 30 days (Figure 4 (d)). For each plot, the numbers of exposed and infected start low and then increase. Once they reach their peak, they decrease to low numbers again. When $\alpha_1 = -0.01$ and $\alpha_1 = -0.025$, both the numbers of exposed and infected reach high levels when the time $t$ is about 15 days, while fewer people are exposed and infected when $\alpha_1 = -0.05$ and $\alpha_1 = -0.065$, respectively. The numbers of recovered and deceased individuals increase. The numbers of recovered and deceased decrease as the exponential declining rate decreases exponentially from –0.01 to –0.065, e.g., the curve of recovered in Figure 4(a) is higher than that in Figure 4(d).

In the supplementary materials, we consider another set of incidence, recovery, and mortality rates from Figure S1: $\beta(t) = 0.01(19.25 -0.17t), \; \mu (t) = \; 0.01 (14.65 + 0.20 t),$ and $\delta (t) = 0.01 (5.25 + 0.045t -0.001 t^2)$. In this set of rates, the incidence rate is fitted by linear regression, instead of a quadratic regression in Figure 3. Again, the initial conditions are given by $S_0 = 10^6, E_0 = 10$, and $I_0 = 2$. Using the rates, approximate numbers of susceptible, exposed, infected, and deceased are given in Figure S2, which provides similar results as Figure 4.

In Figure 5, we use another set of rates to calculate approximate numbers of susceptible, exposed, infected, and deceased: $\beta(t) = 0.01 (14.5 -0.25t), \mu (t) = 0.01 (8.5 + 0.30 t),$ and $\delta (t) = 0.01 (4.5 -0.05t -0.001 t^2)$. The features of the numbers of susceptible, recovered, and deceased are similar to that of Figure 4. However, the number of exposed is still high and does no reach a low level when time $t$ is around 50 days, and the number of infected does not stabilize when time $t$ is around 50 days. Hence, this set of rate functions leads to a good fit for the transmissions and progressions of US COVID-19 pandemic.

Figure 5.  Approximate numbers of susceptible, exposed, infected, and deceased using exponential decline contact function $\alpha (t) = 1.2 \times 10^{-6} \exp (\alpha_1 t)$, (a) $\alpha_1 = -0.01$, (b) $\alpha_1 = -0.05$, (c) $\alpha_1 = -0.065$, and (d) $\alpha_1 = -0.075$, when other rates are fixed: $\beta(t) = 0.01 (14.5 -0.25t), \mu (t) = 0.01 (8.5 + 0.30 t),$ and $\delta (t) = 0.01 (4.5 -0.05t -0.001 t^2)$. The initial conditions are given by $S_0 = 10^6, E_0 = 10$, and $I_0 = 2$.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

In this article, we build SEIRD differential equations to describe relationships among the numbers of susceptible, exposed, infected, recovered, and deceased individuals. Based on the empirical rates, empirical analyses for the SEIRD models are carried out to explore the applications and implications of the models. For COVID-19, a susceptible individual may interact with an exposed person and contract the virus. Hence, an exposed individual is infectious and can transmit the virus to susceptible people rapidly. The SEIRD system differs from SEIR and SIR differential equations in the literature. In these literature models, only an infected individual is infectious and can transmit the virus to a susceptible person. The proposed SEIRD system can directly describe the transmission of COVID-19 in a susceptible population.

Using the US empirical incidence, mortality, and recovery rate functions of COVID-19, we calculate approximate numbers of susceptible, exposed, infected, recovered, and deceased. The incidence, recovery, and mortality rate functions are functions of time in days and are not constant. We consider four exponential decline contact functions to investigate the impact of contact rate on the transmission of COVID-19. An exponential declining contact rate implies that people have contact with each other as usual at the beginning, with a reduction in contact quickly later on. Our empirical results show that a very rapid decline in contact significantly reduces the numbers of exposed, infected, and deceased. Therefore, it makes sense to wear masks and reduce social activities during the pandemic.

In the literature, the SEIR models were extended to be a susceptible-exposed-infected-recovered-susceptible (SEIRS) models, where recovered people may become susceptible again, i.e., recovery does not confer lifelong immunity ^[25]. For example, rotavirus and malaria are diseases with long incubation periods, and recovery only confers temporary immunity. In the SEIRD system (2.1), we assume immunity for COVID-19 and that recovery confers lifelong immunity. The models can be extended to remove this assumption. As time progresses, we may have empirical data supporting the fact that recovery does not confer lifelong immunity for COVID-19; SEIRD system (2.1) can then be extended to accommodate the finding. Another extension of SEIRD system (2.1) allows for immunity by vaccination. Because vaccine rollout began in December 2020, its effectiveness and influence can be evaluated and provide us with data to help us develop an extension of the SEIRD system that will accommodate immunity by vaccination.

In this article, we only consider a deterministic system of SEIRD differential equations which can not capture random fluctuations of the numbers of susceptible, exposed, infected, recovered, and deceased individuals. The deterministic models can be extended to include stochastic components to describe random fluctuations of the virus transmissions of COVID-19 and viral phylodynamics of SARS-CoV-2 ^{[36,37,38,39,40,41,42,43,44,45]}. More research is needed for further investigations.

Acknowledgments

This study was supported by U.S. National Science Foundation grant DMS-1915904 (Ruzong Fan). One anonymous reviewer and Dr. Jianhong Wu, provided very good and insightful comments for us to improve the manuscript. Mr. Marion L. Hartley helped us to edit the paper. Ruzong Fan conceives the project; Shuqi Wang, Wen Tang, Liyan Xiong, Menyu Fang, and Chi-Yang Chiu provide technical support for the project.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Supplementary

In the supplementary materials, we consider another set of incidence, recovery, and mortality rates from Figure S1: $\beta(t) = 0.01(19.25 -0.17t), \mu (t) = 0.01 (14.65 + 0.20 t),$ and $\delta (t) = 0.01 (5.25 + 0.045t -0.001 t^2)$. In this set of rates, the incidence rate is fitted by linear regression, instead of a quadratic regression in Figure 3. Again, the initial conditions are given by $S_0 = 10^6, E_0 = 10$, and $I_0 = 2$. Using the rates, approximate numbers of susceptible, exposed, infected, and deceased are given in Figure S2, which provides similar results as Figure 4.

Figure S1.  The testing rates, incidence rates, mortality rates, and recovery rates as well as related quadratic and linear regression lines for COVID-19 in the United States.

DownLoad: Full-Size Img PowerPoint

Figure S2.  Approximate numbers of susceptible, exposed, infected, and deceased using exponential decline contact function $\alpha(t) = 1.2 \times 10^{-6} \exp (\alpha_1 t)$, (a) $\alpha_1 = -0.01$, (b) $\alpha_1 = -0.025$, (c) $\alpha_1 = -0.05$, and (d) $\alpha_1 = -0.065$, when other rates are fixed: $\beta(t) = 0.01(19.25 -0.17t), \mu (t) = 0.01 (14.65 + 0.20 t),$ and $\delta (t) = 0.01 (5.25 + 0.045t -0.001 t^2)$. The initial conditions are given by $S_0 = 10^6, E_0 = 10$, and $I_0 = 2$.

DownLoad: Full-Size Img PowerPoint