A generalized distributed delay model of COVID-19: An endemic model with immunity waning

1.   Introduction

The first outbreak of the novel coronavirus disease 2019 (COVID-19) occurred in December 2019 in the city of Wuhan, Hubei province, China ^[1]. As of July 2022, the disease is still spreading around the world with over 500 million reported cases worldwide and over 6.1 million reported deaths ^[2]. The COVID-19 disease is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which is believed to have originated from bats ^[3]. It is still unclear how the virus made its way into human population. In human populations, the virus is transmitted from one person to another through the inhalation of infectious aerosols that are dispersed when an infected individual talks, sneezes or coughs. The virus can also be transmitted when a susceptible individual has in-person contact with an infected individual (e.g., hugs, handshakes, kisses, etc.) or touches a contaminated surface and uses the hand to touch their eyes, nose or mouth ^[4,5,6,7]. SARS-CoV-2 is similar to other coronaviruses, such as the severe acute respiratory syndrome (SARS) and the Middle East respiratory syndrome coronavirus (MERS). However, it is believed to spread more easily than these other viruses ^[8,9].

The world health organization (WHO) declared COVID-19 as a public health emergency on January 20, 2020 ^[10] and a pandemic on March 11, 2020 ^[11]. Before this declaration, many governments worldwide have already implemented non-pharmaceutical interventions (NPIs), such as physical distancing, wearing of facemasks, closure of schools and businesses, and travel restrictions, among others, to reduce the spread of the disease. Although implementing these NPIs helped slow the spread of COVID-19, it does not prevent its continuous spread from one region to another. To further control the spread of SARS-CoV-2, there was an urgent need for safe and effective vaccines against the virus ^[12,13]. Vaccination was believed to be the only way to bring the pandemic under control without causing much damage to our physical and mental health and also without impacting the world economy much further than it has already been affected ^[14,15,16]. The first set of COVID-19 vaccines became available towards the end of 2020 ^[17]. These vaccines provide significant protection against the original strain of SARS-CoV-2 virus ^[18,19,20]. However, the emergence of highly transmissible variants of the virus, such as the Delta and Omicron variants ^[21,22] has led to the continuous spread of the disease. In addition, the waning of both vaccine-induced and natural immunities has led to conversations about the endemicity of COVID-19 ^[23,24].

Following the outbreak of COVID-19, many mathematical models have been developed to study the transmission dynamics of the disease ^{[25,26,27,28,29]}. In addition, mathematical models have been used to study the effectiveness of intervention strategies (pharmaceutical and non-pharmaceutical) implemented during the early stages of the pandemic ^{[30,31,32,33,34,35,36]}. One of the earliest mathematical models proposed for studying the dynamics of COVID-19 includes the model of ^[29], which focuses on the transmission of the disease by super-spreaders. In ^[28], an age-structured model was used to study the transmission dynamics of COVID-19 in British Columbia, Canada. The authors of this article proposed a framework that incorporates empirical contact surveys into disease transmission models. They showed that time-dependent contact rates are required to model the dynamics of COVID-19 adequately. The work of ^[32] introduced a Bayesian epidemiological model of COVID-19, which explicitly accounts for individuals willing and able to participate in physical distancing. By fitting their model to the reported cases of COVID-19 in British Columbia (BC), Canada, they estimated the effect of physical distancing on contact rate and studied the projected effect of relaxing this measure on the transmission of COVID-19 in BC. Their results show that physical distancing has a strong impact on COVID-19 transmission in BC, and it is consistent with a decline in reported cases and hospitalization.

Many of the mathematical models used to study the dynamics of COVID-19 considered an epidemic modeling framework. Towards the middle of the year 2020, there have been several reports and discussions on COVID-19 reinfection ^[37,38,39]. By mid-2021, there were reports of breakthrough infections (vaccinated individuals getting infected) ^[40,41,42] and waning of both vaccine-induce and natural immunities ^[43,44,45]. These developments triggered discussions about the endemicity of COVID-19 and the use of endemic models to study the long term behavior of the disease ^{[46,47,48,49,50,51]}. In this study, we develop an endemic model of COVID-19 that incorporates the waning of both vaccine-induced and natural immunities. Many of the endemic models proposed to study COVID-19 assume that the waning rate of both vaccine-induced and natural immunities follow exponential distributions. However, it is well-known that immunity waning happens differently between individuals and over a period of time in the population ^[23,24,52]. Our modeling framework uses distributed delay equations to capture the dynamics of immunity waning in a population for both vaccine-induced and infection-induced (natural) immunities. These distributed delay equations assume that the time to lose immunity in a population follows a Gamma distribution and immunity loss happens gradually and not instantaneously as in the case of an exponential distribution. A similar approach was used to model immunity waning in the susceptible-infected-hospitalized-recovered-deceased (SIHRD) model of ^[53] and the immuno-epidemiological model of ^[54].

We present the derivation of our distributed delay model in Section 2 and use the linear chain trick method to derive a nonlinear ODE system from the integro-differential equation model in Section 2.1. We show that the reduced ODE system is biologically well-posed by proving the non-negativity and boundedness of its solutions in Section 3. In Section 4, we study the stability properties of the ODE system and performed bifurcation analysis on the endemic equilibrium. Furthermore, we study the global stability of the disease-free equilibrium in Section 4.2 and that of the endemic equilibrium in Section 4.3. In Section 5, we present numerical simulations of the ODE system for hypothetical scenarios. Lastly, we discuss our modeling results in Section 6.

2.   Mathematical model

We develop an endemic model of COVID-19 that uses distributed delay differential equation to model the waning of both vaccine- and infection-induced immunities. The model consists of six compartments: susceptible ($S$), vaccinated $(V)$, exposed ($E$), asymptomatic infectious $(I_A)$, symptomatic infectious $(I_S)$, and recovered $(R)$. There is a constant inflow of individuals into the susceptible population (due to birth and immigration) at a rate $\Lambda$. A susceptible individual becomes infected at a rate $\beta$ upon coming in contact with an infectious individual in either $I_A$ or $I_S$, and transitions to the exposed population ($E$). In addition, susceptible individuals transition to the vaccinated compartment $(V)$ at a vaccination rate $\pi$. Vaccinated individuals can become infected when they come in contact with infectious individuals (breakthrough infection) or lose their vaccine-induced immunity and return to the susceptible compartment. Following a similar approach used in ^[55], we assume that the time it takes for vaccine-induced immunity to wane in individuals differ, and that waning happens gradually at the population level ^[23,24,52].

Let $p(k)$ be the probability density function of the random variable $k$, which models the time it takes for an individual to lose their COVID-19 vaccine-induced immunity, starting from when they got vaccinated. Then the probability that an individual loses his/her vaccine-induced immunity between $0$ and $\tau_1$ units time after obtaining the vaccine is $\int_0^{\tau_1} p(k) \; \text{d}k$. Thus, the probability that a vaccinated individual does not loose their vaccine-induced immunity between $0$ and $\tau_1$ units time after vaccination is given by

which is the same as the probability that they lose their immunity after $\tau_1$. We assume that at time $t - \tau_1$, $\pi\, S(t - \tau_1)$ susceptible individuals were vaccinated and transitioned to the vaccinated compartment. Vaccinated individuals that lose their vaccine-induced immunity transition back to the susceptible population, from which they can be infected. The probability that a vaccinated individual is still alive or has not left the community under consideration at $\tau_1$ time units after vaccination is given by $e^{-\mu\, \tau_1}$, where $\mu$ is the per capital death/emigration rate for individuals in the population. Due to imperfection in the vaccine, we assume that vaccinated individuals get infected with COVID-19 at a rate $\xi \equiv \xi(I_A, I_S)$, where $I_A$ and $I_S$ are the asymptomatic and symptomatic infectious populations, respectively, ($\xi = 0$ for a perfect vaccine), and define the probability that a vaccinated individual is not yet infected at time $\tau_1$ as $e^{-\xi\, \tau_1}$. Therefore, the population of vaccinated individuals at time $t$ is given by

Using the change of variable $\hat{\tau_1} = t - \tau_1$, (2.2) becomes

Upon differentiating both sides of (2.3) with respect to $t$ and using Leibniz integral rule, we derive the differential equation for the vaccinated population as

The first term on the right hand side of (2.4) represents the newly vaccinated individuals transitioning to the vaccinated compartment from the susceptible compartment, and the second term accounts for the vaccinated individual that got infected (breakthrough infection). The integral in (2.4) represents the vaccinated individuals losing their vaccine-induced immunity over a period of time. These individuals eventually return to the susceptible population. The last term accounts for emigration and death of vaccinated individuals.

A susceptible or vaccinated individual become infected at some rate after coming in contact with an infectious individual, and moves to the exposed compartment before transitioning to become asymptomatic or symptomatic infectious at a constant rate $\alpha$. It is important to note that infected individuals in the exposed compartment are not infectious and do not transmit the disease. A fraction $\eta$ of these individuals transitioning from the exposed to the infectious stage of the disease become asymptomatic, while the remaining $1-\eta$ are symptomatic. Asymptomatic and symptomatic infectious individuals recover from the disease at the constant rates $\gamma_A$ and $\gamma_S$, respectively, and transition to the recovered compartment $(R)$ upon recovery.

We also assume that recovered individuals can lose their infection-induced (natural) immunity, and model this process using a distributed delay equation. In this case, we let $l$ be a random variable that models the time it takes for an individual to lose his/her infection-induced immunity starting from the time after recovering from COVID-19, and let $q(l)$ be the probability density function for $l$. Therefore, the probability that an individual loses his/her immunity between $0$ and $\tau_2$ units time after recovery is $\int_0^{\tau_2} q(l) \; \text{d}l$. On the other hand, the probability that a recovered individual does not lose their infection-induced immunity between $0$ and $\tau_2$ units time after recovery is

Assuming that at time $t - \tau_2$, $\gamma_A\, I_A(t - \tau_2)$ and $\gamma_S\, I_S(t - \tau_2)$ individuals recovered from the disease, thereby, transitioning from the asymptomatic and symptomatic infectious compartments, respectively, to the recovered compartment $(R)$. Since the probability that a recovered individual is still alive at $\tau_2$ time units after recovery is given by $e^{-\mu\, \tau_2}$, therefore, the total number of individuals in the recovered compartment at time $t$ is given by

Similar to (2.2), we apply a change of variable to (2.6) and differentiate both sides of the resulting equation. Using Leibniz integral rule, we obtain the differential equation for $R(t)$ given by

Here, the first term on the right hand side of (2.7) accounts for the individuals that just recovered from COVID-19, while the integral accounts for those losing their infection-induced immunity. These individuals return to the susceptible population gradually over a period of time. The last term in this equation accounts for emigration and death.

In both (2.4) and (2.7), we assume that the random variables $k$ and $l$ follow the Gamma distribution, so that $p(k)$ and $q(l)$ are probability density functions for the Gamma distribution, given by

where $n, m = 1, 2\dots$ and $\theta, \phi > 0$ are the shape and rate parameters, respectively. Based on the derivation described above, the differential equations for our distributed delay model are

where $\pi$ is the vaccination rate, $\alpha$ is the rate of transitioning from exposed to infectious, $\eta$ is the fraction of exposed individuals that become asymptomatic, and $\gamma$ is the recovery rate. Here, $\Lambda$ is the birth/immigration rate of susceptible individuals into the population, $\mu$ is the natural death/ immigration rate, and $\delta_A$ and $\delta_S$ are the COVID-19 induced death rates for the asymptomatic and symptomatic infectious populations, respectively. In (2.9), $\lambda \equiv \lambda(I_A, I_S)$ is the force of infection, which determines the disease spread in the population. It is defined as

where $\beta_A$ and $\beta_S$ are the disease transmission rates for the asymptomatic and symptomatic infectious populations, respectively, and $N(t) = S(t) + V(t) + E(t) + I_A(t) + I_S(t) + R(t)$ is the total population. It is important to remark that we have used the explicit definition for the rate at which vaccinated individual get infected, $\xi = (1-\varepsilon) \lambda$ in (2.9), where $0 \leq \varepsilon < 1$ is the vaccine efficacy and $\lambda$ is the force of infection (2.10). Here, $\varepsilon = 0$ implies that the vaccine is not effective in blocking COVID-19 infection at all, while $\varepsilon = 1$ implies that the vaccine perfectly blocks the acquiring of COVID-19 infection.

2.1.   Ordinary differential equation (ODE) derivation

Next, we use the linear chain trick method ^{[55,56,57,58,59]} to derive a nonlinear ordinary differential equation (ODE) model from the distributed delay differential equation (2.9). Our goal is to study the stability properties of the ODE model, and the long-term dynamics of COVID-19.

We begin by considering the integrals in (2.9). We let

where $p_{\theta, j}(\tau_1)$ for $j = 1, 2, \dots, n$ ($n = 1, 2, \dots$) and $\theta > 0$, is the Gamma distribution function given in (2.8). Here, $V_w$ represents vaccinated individuals losing their vaccine-induced immunity. We want to derive an ODE representation of this integral for an arbitrary $n$ following the approach of ^[59]. Using the change of variable $\tau = t - \tau_1$, we have

Differentiating both side of (2.12) with respect to $t$ and using the Leibniz integral rule, we obtain

From (2.8), we derive that the Gamma distribution satisfies the following initial value problems

Upon using the initial value problems (2.14) in (2.13), we derive an ODE system for $V_{w \, j}(t)$, given by

Similarly, let $R_{w \, j}$ represent the recovered individuals losing their infection-induced immunity and set

where $q_{\phi, j}(\tau_1)$ for $j = 1, 2, \dots, m$ ($m = 1, 2, \dots$) and $\phi > 0$, is the Gamma distribution function given in (2.8). We use the change of variable $\hat{\tau} = t - \tau_2$ in (2.16) to obtain

Upon differentiating both side of (2.17) with respect to $t$ and following a similar procedure used in deriving the ODEs for $V_{w \, j}$, we derive the following ODEs for $R_{w \, j}$

Next, we couple the system of ODEs for $V_{w\, j}(t)$ and $R_{w \, j}(t)$ as defined in (2.15) and (2.18), respectively, to (2.9) to obtain our reduced ODE model given by

with the initial conditions $S(0) > 0, V(0) \geq 0, E(0) \geq 0, I_A(0) \geq 0, I_S(0) \geq 0, R(0) \geq 0$,

Here, the total population is given by $\hat{N}(t) = S(t) + V(t) + E(t) + I_A(t) + I_S(t) + R(t) + \sum_{j = 1}^n V_{w \, j}(t) + \sum_{j = 1}^m R_{w \, j}(t)$, where the ODEs for $V_{w\, j}(t)$ for $j = 1, \dots, n$ and $R_{w \, j}(t)$ for $j = 1, \dots, m$ are as given in (2.15) and (2.18), respectively. The model variables and their descriptions are given in Table 1.

In Figure 1, we present the schematic diagram for the nonlinear ODE system (2.19), where the black arrows show the progression of individuals through the different compartment of the model and disease stages. The blue arrow shows birth and immigration into the susceptible population, the red arrows show both natural and COVID-19 induced deaths, while the green arrows show immunity waning. Immunity is lost in stages and over a period of time. In the schematic diagram, $V_{w \, j}$ represents the vaccinated individuals in the $j^{th}$ stage of losing their vaccine-induced immunity, and $R_{w \, j}$ represents recovered individuals in the $j^{th}$ stage of losing their infection-induced immunity.

3.   Non-negativity and boundedness of solution

In this section, we consider the ODE model (2.19) for $n = m = 2$, and show the non-negativity and boundedness of its solution. We consider the ODEs (2.15) and (2.18) for $j = 1, 2$, which implies that an individual in the process of losing their immunity will pass through two intermediate immunity waning stages before eventually losing the immunity, namely, $V_{w \, 1}$ and $V_{w \, 2}$ for the individuals losing their vaccine-induced immunity, and $R_{w \, 1}$ and $R_{w \, 2}$ for those losing their infection-induced immunity.

We consider (2.11) for $j = 1$ and substitute $p_{\theta, 1}(\tau_1) = \theta \, e^{-\theta\, \tau_1}$ to obtain

Similarly, we substitute $p_{\theta, 1}(\tau_1) = \theta \, e^{-\theta\, \tau_1}$ into $V(t)$ as defined in (2.2) and evaluate the inner integral to obtain

Upon comparing (3.1) and (3.2), we deduce that $V_{w \, 1}(t) = \theta \, V(t)$. Following a similar procedure for $R(t)$ and $R_{w, 1}(t)$ as defined in (2.6) and (2.16), respectively, we obtain $R_{w \, 1}(t) = \phi \, R(t)$. Using this relations in (2.19), we derived the corresponding ODE system for the case of $n = m = 2$, given by

}where $V_{w \, 2}(t) = \theta \, \hat{V}(t)$ and $R_{w\, 2}(t) = \phi\, \hat{R}(t)$, with $V_{w \, 2}(t)$ and $R_{w\, 2}(t)$ satisfying the ODEs in (2.15) and (2.18), respectively, for $j = 2$. Here, the total population is given by

and the initial conditions are $S(0) \geq 0, V(0) \geq 0, \hat{V}(0) \geq 0, E(0) \geq 0, I_A(0) \geq 0, I_S(0) \geq 0, R(0) \geq 0$, and $\hat{R}(0) \geq 0$. Observe that we have eliminated the equations for the first stage of immunity waning $V_{w\, 1}(t)$ and $R_{w\, 1}(t)$ in (3.3), using the properties $V_{w\, 1}(t) = \theta\, V(t)$ and $R_{w\, 1}(t) = \phi \, R(t)$, respectively.

Now, we show the non-negativity and boundedness of the solutions of the ODE system (3.3), which guarantees that the model is biologically well-posed. To show that the solution of this system remains positive starting with a positive initial condition, we start by considering the equation for the susceptible population in (3.3), which gives

Observe that (3.5) is a first-order linear ODE with integrating factor

where $\lambda(s) \equiv \lambda(S, V, E, I_A, I_S, R, \hat{V}, \hat{R})$. Upon solving (3.5) using this integrating factor and with equality, we obtain

Evaluating the integral on the right hand side of this equation and using the fact that $S(0) \geq 0$, we get

which can easily be simplified to

Since $S(0) > 0$ and $\Lambda > 0$, from (3.6), it is guaranteed that the solution $S(t)$ is non-negative for all $t > 0$. From the equation for $V(t)$ in (3.3), we have $V'(V = 0) = \pi\, S(t)$. Since $S(0) > 0$ and $S(t) \geq 0$, we conclude that $V(t) \geq 0$ for all $t > 0$. Following a similar procedure, it can easily be shown that $E(t), I_A(t), I_S(t), R(t), \hat{V}(t)$ and $\hat{R}(t)$ are non-negative for all $t > 0$.

Next, to show that the solutions of (3.3) are bounded for all $t > 0$, we differentiate $\hat{N}(t)$ as given in (3.4) with respect to $t$ to obtain

Note that we have used the fact that $I_A \geq 0$ and $I_S \geq 0$ in deriving the inequality in (3.7). The solution of this inequality is given by

Taking the limit of both sides as $t \to \infty$ gives

This inequality shows that the total population $\hat{N}(t)$ is bounded by $\Lambda/\mu$ for all $t > 0$. In addition to the boundedness of the solutions, since all the state variables are non-negative for all $t > 0$, we are guaranteed that the ODE system (3.3) is well-posed biologically. Although the proofs presented here are for $n = m = 2$ in (2.19), they can easily be extend to the case of finite $n \in \mathbb{Z}^+$ and $m \in \mathbb{Z}^+$ to show that the ODE system (2.19) is biologically well-posed.

4.   Stability analysis

Next, we study the stability properties of the ODE model (3.3). We first derive the disease-free steady-state of the system, after which we derive the control reproduction number ($R_c$). In addition, we perform bifurcation analysis on the endemic equilibrium of a limiting model and show that the model exhibits a transcritical bifurcation at $R_c = 1$, which can be either backward or forward depending on the immunity waning rates. Lastly, we proof the global stability of both the disease-free and endemic equilibria.

We derive the disease-free equilibrium (DFE) of the ODE system (3.3) as

where $\Phi = \mu \big[(\mu + \theta)^2 + v\, (\mu + 2\theta) \big]$. The total population at the DFE is given by $\hat{N}_e = S_e + V_e + \hat{V}_e = \Lambda/\mu$. Using the next generation matrix approach ^[60,61,62], we derive the control reproduction number of the ODE model (3.3) as

where $\psi_1 = \gamma_A + \mu + \delta_A$ and $\psi_2 = \gamma_S + \mu + \delta_S$. The control reproduction number, $R_c$ denotes the average number of new COVID-19 infection generated by a single infection introduced into a susceptible population where a fraction of them are vaccinated. We establish the following theorem.

Theorem 4.1. The DFE of the model equation (4.1) is locally asymptotically stable (LAS) if $R_c <$ 1, and unstable if $R_c >$ 1.

The results in Theorem 4.1 implies that COVID-19 can be eliminated from the population when the control reproduction number $R_c < 1$ if the initial size of the population is under the basin of attraction of the DFE, $\Theta$. The proof of 4.1 is elementary and can be established using Theorem 2 of ^[61].

4.1.   Bifurcation analysis

We consider a limiting ODE model where vaccinated individual can only get infected with COVID-19 after completely losing their vaccine-induced immunity. We perform bifurcation analysis on the endemic equilibrium of the limiting model. In this limit, the ODE system (3.3) reduces to

The disease-free equilibrium of the model (4.3) is the same as that of the ODE system (3.3), given in (4.1). Based on the assumptions used to derived the model (4.3), the control reproduction number (4.2) reduces to

where $\psi_1 = \gamma_A + \mu + \delta_A$ and $\psi_2 = \gamma_S + \mu + \delta_S$.

We begin our analysis by defining the endemic equilibrium of (4.3) as

and the force of infection at the endemic equilibrium point as

where $N_d(t) = S_d(t) + V_d(t) + \hat{V}_d(t) + E_d(t) + I_{A d}(t) + I_{S d}(t) + R_d(t) + \hat{R}_d(t)$ is the total population at the endemic equilibrium. At the endemic equilibrium, the state variables of (4.3) satisfy

where

Note that $\psi_1$ and $\psi_2$ are as defined in (4.2). Substituting the parameters in (4.8) into (4.7) and writing all the state variables in terms of $\lambda_d$, it can be shown that the non-zero equilibria of (4.3) satisfy the quadratic equation

where the coefficients $T_1, T_2$ and $T_3$ satisfy

Here, the parameters $A_0, \dots, A_6$ are defined as follows

It can easily be shown that $T_1 > 0$ since $A_1(\phi+\mu)-\phi^2A_0 > 0$, and $T_3\geq 0$ if and only if $\hat{R}_c\leq 1$. The endemic equilibrium (equilibria) of the model can be obtained by solving for the positive values of $\lambda_d$ in (4.9) and substituting it into (4.7). These equilibria exist for $\hat{R}_c > 1$ and are summarized in the following theorem.

Theorem 4.2. The ODE model (4.3) has

1) a unique endemic equilibrium if $T_3 < 0 \iff \hat{R}_c > 1$;

2) a unique endemic equilibrium if $T_2 < 0$ and $T_3 = 0$ or $T_2^2-4T_1T_3 = 0$;

3) two endemic equilibria if $T_3 > 0$, $T_2 < 0$ and $T_2^2-4T_1T_3 > 0$;

4) no endemic equilibrium otherwise.

In Theorem 4.2, case (1) shows the existence of the endemic equilibrium point $\Upsilon$ whenever $\hat{R}_c > 1$, while case (3) shows the possibility of a backward bifurcation. The phenomenon of backward bifurcation occurs when an asymptotically stable disease-free equilibrium co-exists with an asymptotically stable endemic equilibrium for some $\hat{R}_c < 1$ ^[33,63,64]. To prove this fact, a critical value of $\hat{R}_c$ denoted by $R_{critical}$ can be obtained from the discriminant equation $T_2^2-4T_1T_3 = 0$. Upon substituting (4.10) into this equation and simplifying, we obtain

where $A_1$ is as defined in (4.11), $T_1$ and $T_2$ are as defined in (4.10), and $\Psi_0$ and $\Psi_2$ are defined in (4.8). Thus, a backward bifurcation occurs at $\hat{R}_c = 1$, giving rise to a branch of unstable endemic equilibrium solutions for $R_{\text{critical}} < \hat{R}_c < 1$. Another bifurcation occurs at $\hat{R}_c = R_{\text{critical}}$, where the unstable branch of endemic solutions switch stability and become a stable branch of endemic equilibrium solutions for $R_{\text{critical}} < \hat{R}_c$. The biological interpretation of having a backward bifurcation is that the control reproduction number $\hat{R}_c$ being less than unity is no longer sufficient to guarantee the disease eradication. More control measures would be required to eradicate the disease from the population.

In Figure 2, we show that (4.3) can exhibit either a forward or backward bifurcation, depending on the immunity waning rates. When we fixed the waning rate of vaccine-induced immunity as $\theta = 0.01$ and used an infection-induced immunity waning rate of $\phi = 0.01$, the model exhibits a forward bifurcation, which occurs at $\hat{R}_c = 1$ (left panel of Figure 2). The bifurcation changes to a backward bifurcation when $\phi$ is increased to $\phi = 0.02$ (right panel of Figure 2). The parameters used are $\beta_A = 0.2189$, $\beta_S = 0.3521$, $\pi = 0.975$, $\eta = 0.25$, and as given in Table 2. Having a forward bifurcation implies that COVID-19 can be eradicated from the population when $\hat{R}_c < 1$. On the other hand, a backward bifurcation implies that $\hat{R}_c < 1$ is not sufficient to eradicate COVID-19 from the population. Observe that a forward bifurcation is obtained with a smaller immunity waning rate, which implies that disease eradication is easier when immunity waning happens over a longer time compared to when it happens faster.

In Figure 3, we study the effect of immunity waning on the bifurcation types. In other words, we study the switching of the bifurcation at $\hat{R}_c$ from forward to backward. For a fixed value of $\theta$, there is a threshold value of $\phi$ where the bifurcation at $\hat{R}_c = 1$ switches from a forward transcritical bifurcation to a backward bifurcation. This threshold value depends on the value of $\theta$ and increases as $\theta$ increases as illustrated in Figure 3. When $\theta = 0.01$, the bifurcation switches at $\phi \sim 0.0175$ (left panel of Figure 3). Increasing the vaccine-induced immunity waning rate to $\theta = 0.015$, the threshold value for the infection-induced immunity where the bifurcation switches forward to backward also increases to $\phi \sim 0.024$ (right panel of Figure 3). It is important to mention that for a fixed value of $\phi$, changes is the value of $\theta$ does not affect the bifurcation type.

4.2.   Global stability analysis of the disease-free equilibrium

Next, we establish the global stability of the disease-free equilibrium (DFE) of the ODE model (4.3) under the condition that the disease-induced mortality rates are negligible, i.e., $\delta_A = \delta_S = 0$. We used the Castillo-Chavez approach ^[63], which uses the Jacobian matrices of both the disease-free and endemic equilibria to establish the stability of the DFE for $\hat{R}_c < 1$. The approach has been used in several articles including ^[64,65,66].

By re-writing the model (4.3) in vector form, we have

where the vector $\mathcal{U} = (S, V, \hat{V}, R, \hat{R})$ denotes the uninfected compartments of the system and $\mathcal{V} = (E, I_A, I_S)$ denotes the infected compartments. Here, $H(\mathcal{U}, \pmb{0}) = \pmb{0}$. For global stability to be ensured, the following conditions must be satisfied:

$1)$ $\dot{\mathcal{\mathcal{U}}}(t) = G(\mathcal{U}^{0}, \pmb{0}) = \pmb{0}$, where $\mathcal{U}^{0}$ is the DFE of (4.3).

$2)$ $\hat{H}(\mathcal{U}, \mathcal{V}) = J\mathcal{V}-{H}(\mathcal{U}, \mathcal{V})$, with $\hat{H}(\mathcal{U}, \mathcal{V})\geq0$, where $J$ is the Jacobian matrix defined by $J = \frac{\partial H}{\partial \mathcal{V}}(\mathcal{U}^{0}, \pmb{0})$.

From the ODE model (4.3), we construct the Jacobian matrix of the infected compartments at the DFE as

Where $S_e, V_e$, and $\hat{V_e}$ are as given in (4.1) and $\hat{N_e}$ is the total population at the DFE. We construct the vector $\mathcal{V}$ of the infected compartments and compute the matrix-vector multiplication $J\mathcal{V}$ to obtain

Furthermore, from (4.3), we construct the vector $H(\mathcal{U}, \mathcal{V})$ as

Using (4.14) and (4.15), we compute $\hat{H}(\mathcal{U}, \mathcal{V}) = J\mathcal{V}-H(\mathcal{U}, \mathcal{V})$ as

Since $\frac{S_e}{\hat{N}_e}\geq \frac{S}{\hat{N}}$,

we conclude that $\hat{H}(\mathcal{U}, \mathcal{V}) \geq 0$. It can easily be shown that $\lim_{t\to\; \infty} \mathcal{U}(t) = \mathcal{U}^{0}$ and that $J$ is an $M-$Matrix. Therefore, $\mathcal{U}^0$ is globally asymptotically stable whenever $\hat{R}_c < 1$. This result is summarized in the following theorem:

Theorem 4.3. For the special case where $\delta_A = \delta_S = 0$, the unique positive disease-free equilibrium (DFE) of the model (4.3) is globally-asymptotically stable if the control reproduction number $\hat{R}_c < 1$.

4.3.   Global stability analysis of the endemic equilibrium

Lastly, we establish the global stability of the unique endemic equilibrium of the model (4.3), for the special case $\delta_A = \delta_S = 0$, that is, no disease-induced death. For this analysis to be possible, the control reproduction number must be greater than unity, since this endemic equilibrium only exists for $\hat{R}_c > 1$. We establish the following theorem.

Theorem 4.4. For the special case where $\delta_A = \delta_S = 0$, the unique positive endemic equilibrium of the model (4.3) is globally-asymptotically stable if the control reproduction number $\hat{R}_c > 1$.

Proof. Consider the following nonlinear Lyapunov function

where $H_k$ is a constant and $F_k$ is given by

Here, $f_{d\, k}$ for $k = 1, \dots, 8$ are the elements of the vector $f_{d} = (S_d, V_d, \hat{V}_d, E_d, I_{Ad}, I_{Sd}, R_{d}, \hat{R}_d)$. From (4.3), we have

whose time derivative is given by

Since the disease-induced death rates $\delta_A = \delta_S = 0$, from (3.7), we have

whose solution is given by

Substituting for $\hat{N}$ in (4.6) gives $\lambda_d = q\, (\beta_A\, I_{Ad}+\beta_S\, I_{Sd})$, where $q = \mu/\Lambda$. Substituting the state variables from (4.7) into the expressions in (4.16), we obtain

where $\chi_1, \chi_2, \chi_3, \bar{\chi_4}, \chi_5, \bar{\chi_6}$, and $\chi_7$ are defined as follows:

The constants $H_1$ to $H_6$ are also given by

Upon substituting the expressions in (4.18) into (4.17), the first term gives

The addition of the second and third terms gives

Furthermore, the addition of the forth and fifth is given by

The sum of the sixth and seventh terms in (4.17) is given by

and the last term in (4.17) is simplified to obtain

Combining, simplifying and factorizing the above expressions, we obtain

From (4.19), we have the following

Since $S\geq 0, V\geq 0, \hat{V}\geq 0, E\geq 0, I_A\geq 0, I_S\geq 0, R\geq 0$ and $\hat{R}\geq 0$, and the provision of Theorem 4.4 is satisfied, then $\dot{F}\leq 0$ since all the model parameters are positive. Furthermore, $\dot{F} = 0$, if $S = S_d, V = V_d, \hat{V} = \hat{V}_d, E = E_d, I_A = I_{Ad}, I_S = I_{Sd}, R = R_d$, and $\hat{R} = \hat{R}_d$ Thus, $F$ is a Lyapunov function for the ODE model (4.3). By LaSalle's Invariance Principle ^[67], the model has an endemic equilibrium point that is globally asymptotically stable. This result shows that when the COVID-19-induced death rate is negligible, the virus will persist in the environment if the control reproduction number is greater than unity.

5.   Numerical simulation

We study the ODE model (3.3) numerically for different vaccination scenarios and immunity waning rates. This study helps us understand the effect of vaccination and other important parameters in the model on the disease dynamics. The model parameters and descriptions are given in Table 2.

We begin by studying the effect of the model parameters on the control reproduction number, $R_c$ (4.2). Of particular interest are the vaccination and vaccine-induced immunity waning rates. Our results for this study are presented in the form of contour plots. In the top left panel of Figure 4, we have the contour plot for $R_c$ in terms of the fraction of exposed individuals that become asymptomatic infectious $(\eta)$ and the waning rate of vaccine-induced immunity ($\theta$). This result shows that for any given values of $\eta$, $R_c$ increases as the immunity waning rate increases. Conversely, for any given value of $\theta$, the control reproduction number decreases as $\eta$ increases. For this result, the maximum values of $R_c$ occurs when the asymptomatic fraction is smallest and the waning rate is highest. It is important to remark that this result is influence by the fact that our asymptomatic infectious transmission rate ($\beta_A$) is smaller then our symptomatic infectious transmission rate ($\beta_S$). A similar result is shown in the top right panel of Figure 4 for the transmission rate of asymptomatic infectious population ($\beta_A$) and the waning rate of vaccine-induced immunity ($\theta$). As one would expect, $R_c$ increase with respect to both $\beta_A$ and $\theta$. Although, for a fixed transmission rate, $R_c$ increases more rapidly as the vaccine-induced immunity waning rate increases for high values of $\beta_A$ compared to lower values. This shows that when $\beta_A$ is small, changes in the waning rate of vaccine-induced immunity has little effect on the spread of the disease. However, as the transmission rate increases, the effect of immunity waning on the disease spread also increases.

In the bottom left panel of the same figure, we present the result for the vaccination rate ($\pi$) and the vaccine-induced immunity waning rate ($\theta$). Similar to the other results, we show that $R_c$ increases as $\theta$ increase for any fixed values of the vaccination rate. On the other hand, $R_c$ decreases as the vaccination rate increases for any fix values of the vaccine-induced immunity waning rate ($\theta$). Lastly, in the bottom right panel of Figure 4, we present the result for vaccine efficacy ($\varepsilon$) and the vaccine-induced immunity rate $\theta$. We observe from this result that $\theta$ has little effect on the control reproduction number when the vaccine efficacy is small and its effect on $R_c$ increases as the vaccine efficacy increases. This result shows that when a vaccine is not effective, its waning rate has little effect on the spread of COVID-19. On the other hand, the waning rate of more effective vaccines have significant effect on the spread of the disease.

Next, we study the effect of the model parameters on the disease dynamics. Throughout this study, we assume a total population of 1.5 million, based on the typical population of a large city in North America ^[81], and assume a daily birth and immigration (population increase) rate of 1, 150 per day ^[68]. We assumed that $5\%$ of the population had already recovered from COVID-19 at the beginning of the modeling period. This is used in computing the initial susceptible population using the formula

In addition, we assume that there are no vaccinated individuals, and individuals losing both vaccine- and infection-induced immunities at the beginning of the study. We assume that there are 1000 exposed and 2500 infectious individuals in the population at the beginning of our simulations, and that $20\%$ of the infectious individuals are asymptomatic while the remaining $80\%$ are symptomatic. A summary of the initial conditions is given in Table 3.

We study the effect of both vaccine- and infection-induced immunity waning on the disease dynamics. The results of this study are presented in Figure 5. In the left panel of this figure, we show the dynamics of the infectious population (asymptomatic and symptomatic) over time. We used an infection-induced immunity waning rate of $\phi = 0.01$ and vary the waning rate of vaccine-induced immunity, using $\theta = 0.05, 0.3, 0.5$ and $0.75$. These values of $\theta$ correspond to the control reproduction numbers, $R_c = 0.8065, 1.0773, 1.2402$ and $1.3993$, respectively, computed with (4.2). We observe that the control reproduction number, $R_c < 1$ for only $\theta = 0.05$. This is the only immunity waning rate for which the disease is expected to be eradicated from the population. As $\theta$ increases, so does the control reproduction number, giving rise to $R_c > 1$. For these scenarios, the disease becomes endemic in the population. Similar results are given in the right panel of Figure 5 for $\phi = 0.05$. Since $R_c$ is independent of the waning rate of infection-induced immunity $(\phi)$, the computed $R_c$ for this scenario are the same as those computed for the results in the left panel. However, an increase in $\phi$ leads to more infections in the population.

In Figure 6, we investigate the effect of vaccination and vaccine efficacy on the disease dynamics. We consider different waning rates for vaccine-induced immunity. In the top left panel of this figure, we used a vaccination rate of $\pi = 0.65$ and a vaccine efficacy rate of $\varepsilon = 0.65$. For this scenario, the computed control reproduction numbers are $R_c = 0.7879$ ($\theta = 0.05$), $R_c = 0.9882$ ($\theta = 0.30$), $R_c = 1.1087$ ($\theta = 0.50$), and $R_c = 1.2264$ ($\theta = 0.75$). As expected, the disease is eradicated from the population for $\theta = 0.05$ and $\theta = 0.30$ since the control numbers for these cases are less than 1. On the other hand, the disease continues to spread in the population for the waning rate $\theta = 0.50$ and $\theta = 0.75$ since their $R_c$ are greater than 1. In the top right panel of Figure 6, the vaccination rate is increase to $\pi = 0.80$. We observe from the results in this panel that the spread of COVID-19 in the population is relatively lower compared to when the vaccination rate is lower (top left panel). This shows the impact of an increase in the vaccination rate. The computed control reproduction numbers for this scenario are $R_c = 0.7770, 0.9473, 1.0545$ and $1.1630$ for $\theta = 0.05, 0.30, 0.50$ and $0.75$, respectively.

Next, we consider a scenario where the vaccination rate, $\pi = 0.65$ and vaccine efficacy, $\varepsilon = 0.75$ (bottom left panel of Figure 6). Comparing these results to those in the top left panel of the same figure (for $\pi = 0.65$ and $\varepsilon = 0.75$), we notice that even though the vaccination rates in the two scenarios are the same, there is a decrease in the spread of COVID-19 in the bottom left panel, which can be attributed to an increase in the vaccine efficacy. The computed control reproduction for the scenario in the bottom left panel are $R_c = 0.5893, 0.8204, 0.9594$ and $1.0952$ for $\theta = 0.05, 0.30, 0.50$ and $0.75$, respectively. We observe that the computed $R_c$ for $\theta = 0.50$ in this scenario is less than 1, unlike the case of the same immunity waning rate in the top panel. This suggests that having a safe vaccine with a higher efficacy is important in eradicating the spread of COVID-19. Lastly, by increasing both the vaccination rate and vaccine efficacy to $\pi = 0.80$ and $\varepsilon = 0.75$, respectively, we are able to bring the spread of COVID-19 under control (bottom right panel of Figure 6). For this scenario, the compute control reproduction numbers are $R_c = 0.5768 \; (\theta = 0.05$), $0.7732 \; (\theta = 0.30$), $0.8969 \; (\theta = 0.50)$ and $1.0221 \; (\theta = 0.75$). We observe from the results in this plot that there is a decline in the infectious population over time for all the four cases, and that we are able to eradicate COVID-19 in three of the four cases (where $R_c < 1$). For the last case, where $\theta = 0.75$, although, $R_c > 1$, the infectious population continue to decrease and become relatively small.

6.   Discussion

We have developed an endemic SEIR model of COVID-19 that incorporate the waning of both vaccine- and infection-induced immunities. Our model assumes that the waning of both immunities happens gradually at the population level and describes them using distributed delay equations. The model consists of six major compartments: susceptible $(S)$, vaccinated $(V)$, exposed $(E)$, asymptomatic infectious $(I_A)$, symptomatic infectious $(I_S)$ and recovered $(R)$. We represent the distribution of individuals at different stages of losing their COVID-19 immunity with integrals, and assume that the time it takes to lose both vaccine-induced and infection-induced immunities follow Gamma distributions. The linear chain trick method was used to derive a nonlinear system of ODEs from the distributed delay model, which gave rise to additional compartments in the model. These additional compartments are used to account for the individuals at different stages of losing their immunity.

We proved the non-negativity and boundedness of the reduced ODE model solutions for the case where there are only two stages in losing immunity. These proofs can be extended for the general case with finite number of immunity waning stages, which will guarantee the well-posedness of the ODE model. We computed the control reproduction number of the ODE model and performed bifurcation analysis on the endemic equilibrium of a limiting model. In the limiting ODE model, we assume that a vaccinated individual can only be infected with COVID-19 after completely losing their vaccine-induced immunity. The bifurcation analysis was used to study the stability properties of the endemic equilibrium solution branch of the model. We showed that the limiting ODE model can exhibit both forward and backward transcritical bifurcation at $R_c = 1$, depending on the waning rates of both vaccine-induced and infection-induced immunities. Although the change in bifurcation type depends largely on the waning rate of infection-induced immunity, the threshold value of the waning rate of infection-induced immunity for which the change in bifurcation type occurs changes with respect to the vaccine-induced immunity waning rate. This result suggests that infection-induced immunity waning influences the bifurcation type more than vaccine-induced immunity waning. Furthermore, the existence of backward bifurcation implies that $R_c < 1$ does not guarantee the eradication of COVID-19 in the population, and that additional measures are required to eradicate the disease. This emphasizes the difficulty in eradicating COVID-19, which is as a result of immunity waning. We also performed global stability analysis on both the disease-free and endemic equilibria of the limiting ODE model.

Furthermore, we investigated the effect of model parameters on the control reproduction number. Of particular interest is the effect of the vaccine-induced immunity waning rate. Our results are presented in the form of contour plots for the control reproduction number as a function of the vaccine-induced immunity waning rate and another model parameter such as the fraction of exposed individuals that become asymptomatic, the transmission rate of asymptomatic infectious individuals, vaccination rate, and vaccine efficacy. For all the scenarios considered, our results show that the control reproduction number increases as the immunity waning rate increases, as expected. In addition, we showed that the increase in the control reproduction number as the immunity waning rate increases is more rapid when the asymptomatic transmission rate is higher compared to when it is lower. When COVID-19 spreads at a low rate, the waning rate of the vaccine-induced immunity will have little affect on the disease spread. But when the transmission rate is higher, an increase in the waning rate will lead to significant increase in the spread of the disease. Similar, for the vaccination rate and vaccine efficacy, we noticed that the increase in the control reproduction number as the immunity waning rate increase is more rapid when both parameters are higher compared to when they are low. When the vaccination rate is low, the waning of vaccine-induced immunity has no significant impact on the spread of the disease since only a small fraction of the population has protection from the disease due to vaccination. However, as the vaccination rate increases, vaccine protection against COVID-19 in the population also increases, and as a result of this, an increase in the waning rate of vaccine-induced immunity will have a significant effect on the spread of the disease in the population. A similar explanation applies to the case of vaccine efficacy.

We also studied the effect of the model parameters on the disease dynamics. We started by investigating how the immunity waning rates affect the disease dynamics. As expected, an increase in the waning rate of vaccine-induced immunity leads to an increase in the spread of COVID-19. Recall that the control reproduction number is independent of the waning rate of infection-induced immunity. However, the values of this parameter affects the dynamics of the disease in the population. We observed an increase in infections in the population as this waning rate increases. The control reproduction number is computed at the beginning of our study period and considers new infections and the transfer of infections. It does not depend on the waning rate of infection-induced immunity because this immunity waning happens in the recovered compartment and does not happen at the beginning of the study period. However, when recovered individual lose the infection-induced immunity, they return to the susceptible population and can be reinfected. Lastly, we investigated the effect of vaccination and vaccine efficacy on the spread of COVID-19 in the population by solving the ODE model with various vaccination rates and vaccine efficacies. We started with baseline values for these parameters and then increase them gradually while studying the changes in the infectious population. Our results show that increase in both the vaccination rate and the vaccine efficacy lead to a decrease in both the control reproduction number and the infectious population. In addition, we showed that having a high vaccination rate with a moderately effective vaccine is sufficient to eradicate COVID-19 in the population.

We have proposed a distributed delay framework for modeling COVID-19 and studied the stability properties of the model. An interesting and important next step to this study would involve fitting our reduced ODE model to the reported case of COVID-19 in a particular city, where the estimated parameters would then be used to study the endemicity of COVID-19 in the city. We are currently working on this extension for some selected cities around the world. In this study, we focused on analyzing the reduced ODE system. It would be worthwhile to analyze the complete distributed delay model. Another interesting extension of our modeling framework includes stratifying the population into age groups so that the dynamics of each a group can be studies more extensively. Although, we have presented our modeling framework for a single homogeneous population, the idea can be extended to a scenario with multiple populations, where the different populations interact with each other through the movement of individuals between the populations. Such a model can be used to study the spread of COVID-19 between regions.

Acknowledgments

JDK acknowledges supports from Canada's International Development Research Centre (IDRC) (Grant No. 109981) and the Swedish International Development Cooperation Agency (SIDA) (Grant No. 109559-001). JDK equally acknowledges support from NSERC Discovery Grant (Grant No. RGPIN-2022-04559), NSERC Discovery Launch Supplement (Grant No. DGECR-2022-00454) and New Frontier in Research Fund-Exploratory (Grant No. NFRFE-2021-00879).

Conflict of interest

The authors declare that there is no competing interest.