Spatiotemporal dynamics of a diffusive predator-prey system incorporating social behavior

1.   Introduction

Coronavirus disease (COVID-19) pandemic is considered the biggest global threat worldwide because of thousands of confirmed infections, accompanied by thousands of deaths ^[1]. COVID-19 was identified and named by the World Health Organisation (WHO) on January 10, 2020 following an earlier viral infection episode in Wuhan, China in December 2019, and was declared by the WHO to be a public health emergency of international concern in 2020 ^[2,3]. COVID-19 by its nature is a very contagious disease that spreads easily from person-to-person through direct contact with objects or surfaces that are contaminated with the virus. Moreover, inhalation of respiratory droplets from both asymptomatic and symptomatic infectious individuals cause transmission ^[6]. The symptoms of COVID-19 appear 2-14 days after exposure and may include fever, dry cough, muscle pain, fatigue, and shortness of breath ^[7]. The symptoms are mild in 85% of cases, and vary from severe in 10% to critical in 5% of those infected, but a larger proportion of infected individuals exhibit mild or no symptoms ^[3]. The severity and progression of COVID-19 are known to be exacerbated by the presence of co-morbidities such as diabetes, hypertension and cardio/cerebrovascular diseases ^[8]. It has also been observed that COVID-19 mortality risk is highly concentrated within the elderly population ^[9].

Available scientific evidence classify COVID-19 infected individuals into three broad categories; individuals who manifest severe symptoms, individuals who manifest mild symptoms and individuals who do not manifest any COVID-19 symptoms (asymptomatic) and yet remain infectious undetected. The non-manifestation of COVID-19 symptoms in some infected people complicates the epidemiology of the COVID-19 pandemic. Firstly, asymptomatic individuals are unlikely to seek medical care or self-quarantine given that they cannot tell whether they have the disease unless detected through testing or contact tracing. Secondly, they will continue interacting with healthy people thereby, spreading the virus. Although the asymptomatic categories form a large proportion of COVID-19 infections, it is not yet known to what extent they spread the virus relative to categories with severe symptoms which constitute a small proportion of COVID-19 infections. The size of this delay may play an important role in minimizing the spread of the disease in the community. It is therefore essential to gain a better and more comprehensive understanding of the effects of time delay on COVID-19 transmission and control.

Mathematical models have proved to be essential guiding tools for epidemiologists, biologists as well as policymakers. Models can provide solutions to phenomena which are difficult to measure practically. Recently, a number of mathematical models have been proposed to study the spread and control of COVID-19 (see, for example ^{[1,2,6,7,8,9,19,20,21,22,23,24]}, and references therein) have certainly produced many useful results and improved the existing knowledge on COVID-19 dynamics. In ^[4] a discrete fractional Susceptible-Infected-Treatment-Recovered-Susceptible (SITRS) model for simulating the coronavirus (COVID-19) pandemic was proposed by taking into account the possibility that people who have been infected before can lose their temporary immunity and get reinfected. In ^[5] a novel reaction-diffusion coronavirus (COVID- 19) model was employed to investigate the effect of random movements of individuals from different compartments in their environments. A limitation of these studies, however, is the non-inclusion of the time taken before an infectious human is detected and quarantined, despite the fact that in many countries where the disease is endemic, lack of financial and human resources often results on delay in detection and quarantining of infectious individuals. In addition epidemic models with time delay often exhibit periodic solutions and as a consequence understanding the nature of these periodic outbreaks plays a crucial role in designing policies that can successfully control the disease (see, for example ^[10]). In ^[11] a mathematical model with time delay was proposed to describe the outbreak of 2019-nCoV in China to show that the novel dynamic system can well predict the outbreak trend of the disease. In ^[18] Pei and Zhang constructed a SIRD epidemic model (S-Susceptible, I-Infected, R-Recovered, D-Dead) which is a non-autonomous dynamic system with an incubation time delay to study the evolution of the COVID-19 in Wuhan City, Hubei Province and China Mainland. In ^[25] a system of ordinary differential equations with delays was utilized to describe the evolution of the COVID-19 pandemic.

Although this study is not the first to incorporate discrete delay in studying COVID-19 transmission, the main goal of this article is to explore the dynamics and stability analysis of a COVID-19 model with discrete delay. Hence, we formulated a mathematical model that incorporates a discrete delay that represents the incubation period. In addition, we investigate the impact of the time taken to detect and quarantine infectious individuals on the disease dynamics. The rest of the paper is organised as follows. In Section 2, the propose model present the analytical results. In Section 3, numerical simulations are done to verify the theoretical results presented in the study. Finally, a concluding remark rounds up the paper.

2.   Methods and results

The COVID-19 pandemic remains a major global threat worldwide. This is mainly attributed to several challenges associated with effective control which range from the inadequate use of control measures such as wearing of disposable surgical face masks, regular hand-washing with plenty of soap under running water, the use of alcohol-based hand sanitizers in the absence of soap and water, vaccines among others as recommended by the WHO ^[16,17]. Although these interventions have succeeded greatly in many countries the major problem in the spread of COVID-19 is human-to-human transmission in a heterogeneous community. The implementation of interventional strategies such as quarantine/isolation during infection remains a big challenge in the fight against the disease because of hunger, poverty, and poor health facilities, especially in developing countries in sub-Saharan Africa where governments lack social securities. Furthermore, these challenges often lead to delay in detection and quarantine/isolation of infectious individuals.

2.1.   Mathematical model

In this article, we propose a model to analyze the impact of delay in treatment and the time needed to detect/diagnose and quarantine individuals infected with COVID-19. Assume that the infected individuals are in two different categories, those with mild symptoms and those with severe illness. Anyone can have mild to severe symptoms. We subdivide the total population, $N(t)$, at time $t$, into six compartments namely; susceptible individuals $S(t)$, asymptomatic (undetected infectious) individuals $I(t)$, infectious detected and quarantined $Q(t)$ and recovered individuals $R(t).$ The recovered population $R(t)$ is made up of individuals who have successfully recovered from the infection either naturally or through various health support mechanisms (since the disease has no treatment). The two additional compartments $H_m(t)$ and $H_s(t)$ represent the symptomatic (hospitalized) individuals who develop symptoms. The distinction between the two categories of hospitalized individuals represent that $H_m$ are hospitalized individuals who develop mild symptoms while $H_s(t)$ represent hospitalized individuals who develop severe symptoms. The human population at any given time $t$, is given by $N(t) = S(t)+I(t)+Q(t)+H_m(t)+H_s(t)+R(t).$ The proposed COVID-19 model with a time delay factor is given by:

where $\Lambda$ is the recruitment rate, $\delta$ denotes natural mortality rate, $\beta_1$ denotes the contact rate of asymptomatic (undetected infectious) and susceptible humans, $\beta_2$ denotes the contact rate of symptomatic (hospitalized) and susceptible humans, $\alpha$ is the detection rate of asymptomatic (undetected infectious) patients, $\sigma_1$ denotes the recovery rate of asymptomatic (undetected infectious) and quarantined individuals, $\sigma_2$ is the recovery rate of hospitalized individuals, and $\mu$ represents the disease-induced death rate of symptomatic humans. Proportions $0 < p < 1$ account for hospitalized individuals with severe symptoms, while the remainder $(1-p)$ accounts for hospitalized individuals with mild symptoms. The quarantined individuals are hospitalized after $1/\gamma$ days. $\tau$ is a discrete time delay representing the latent period. The model diagram is depicted in Figure 1.

2.2.   Initial conditions

The initial conditions for the model (2.1) are given as follows:

Here, $C$ is the Banach space $C([-\tau, 0], \mathbb{R}^6)$ of continuous functions mapping the interval $[-\tau, 0]$ into $\mathbb{R}^6$ with the sup-norm $\| \varphi\| = \sup_{\theta\in[-\tau, 0]}\|\varphi_i\|$, for $i = 1, 2, 3, 4, 5, 6$ for $\varphi\in C.$ The non-negative cone of $C$ is defined as $C_+ = C([-\tau, 0], \mathbb{R}^6).$

2.3.   Positivity and boundedness

Equation (2.3) in theorem (2.1) shows that the model formulated in this study is biologically meaningful. Precisely, the theorem demonstrates that for non-negative initial conditions, the solutions of the proposed model are non-negative and bounded for all $t > 0$.

Theorem 2.1. There exists a unique solution for the COVID-19 model (2.1). Furthermore, the solution is non-negative for all $t > 0$ and lies in the set:

Proof. To prove the positivity of the model system (2.1), we investigate the direction of the vector field given by the right-hand side of the model system (2.1) on each space and note whether the vector field points to the interior of $\mathbb{R}^6_+$ or is tangent to the coordinate space. We observe that:

It follows that the vector field given by the right-hand side of the model system (2.1) on each coordinate plane is either tangent to the coordinate plane or points to the interior of $\mathbb{R}^6_+$. Hence, the positivity of the solutions starting in the interior of $\mathbb{R}^6_+$ is assured. $\mathbb{R}^6_+$ is a positively invariant set of the $SIQH_m H_s R$ model system (2.1). Moreover, if the initial conditions $\varphi_i \geq$ 0, (i = 1, 2, 3, 4, 5, 6.) are, therefore, the corresponding solutions of the model system (2.1).

Theorem 2.2. The solutions of the $SIQH_m H_s R$ system (2.1) with the initial conditions of (2.2) are uniformly bounded in the region $\Gamma$.

Proof. To prove the boundedness of the model system (2.1), we add all model equations, which gives:

Since $N'(t) \leq \Lambda-\delta N(t)$ for $0\leq \tau < t,$ it follows by applying the standard comparison Theorem in ^[27] that $N(t) \leq N(0)e^{-\delta t}+(\frac {\Lambda}{\delta})(1-e^{-\delta t}).$ In particular, we have $N(t) \leq (\frac {\Lambda}{\delta})$ if $N(0) \leq (\frac {\Lambda}{\delta}).$ Therefore, we conclude that the population is bounded. Hence, all solutions in $\mathbb{R}_{+}^{6}$ eventually enter $\Gamma$.

2.4.   The basic reproduction number and stability analysis

Since the last equation in the model system (2.1) is independent of the other equations, system (2.1) may be reduced to the following system:

where $m_1 = (\alpha+\delta+\sigma_1)$, $m_2 = (\gamma+\delta+\sigma_1)$ and $m_3 = (\sigma_2+\mu+\delta)$.

It can easily be verified that in the absence of the disease in the community, system (2.6)-(2.10) admit a disease-free equilibrium given by ${\cal E}^0 = (S^0 = \frac{\Lambda}{\delta}, 0, 0, 0, 0)$. In addition, from equations (2.8), (2.9) and (2.10) we have:

Substituting these results into equation (2.7) gives:

From (2.12) we have:

It follows that:

Substituting the value of $S^*$ into equation (2.6) gives:

Substituting (2.15) into (2.11) yields:

From the computations in equation (2.16), we observe that $I^*$, $Q^*$, $H^*_M$ and $H^*_S$ makes biological sense whenever:

Therefore, if we let the basic reproduction number of model (2.6)-(2.10) be:

It follows that models (2.6)-(2.10) has a second equilibrium ${\cal E}^* (S^*, I^*, Q^*, H^*_m, H^*_s)$ point known as the endemic equilibrium which exists whenever ${\cal R}_0 > 1.$

Biologically, the basic reproduction number ${\cal R}_0$ represents the average number of new or secondary COVID-19 infections caused by the introduction of an infectious individual into a totally susceptible population. In fact,

● the term $\frac{\beta_1\Lambda}{m_1\delta}$ is the average number of secondary infections generated as result of contact between susceptible individuals and one asymptomatic (Undetected infectious) COVID-19 patient,

● the term $\frac{\beta_2(1-p)\alpha\gamma \Lambda}{m_1m_2m_3\delta}$ represents the average number of new COVID-19 cases generated when susceptible individuals come into contact with a hospitalized patient of class $H_m$,

● the term $\frac{\beta_2 p\alpha\gamma \Lambda}{m_1m_2m_3\delta}$ gives the average number of secondary COVID-19 infections which occur in the community when susceptible individuals come into contact with a hospitalized patient of class $H_s$.

Then we have the following results:

Theorem 2.3. If ${\cal R}_0\leq 1$, then the disease-free equilibrium ${\cal E}^0$ is globally asymptotically stable.

The detailed proof process can be obtained in Appendix A.

Next, we investigate the global stability of the endemic equilibrium point ${\cal E}^*$ of models (2.6)-(2.10) when ${\cal R}_0 > 1.$

Theorem 2.4. If ${\cal R}_0 > 1$, then model (2.6)-(2.10) has a globally asymptotically stable endemic equilibrium point.

The proofs of Theorem 2.4 is given in Appendix B.

3.   Numerical results and discussions

In this section, we perform numerical analysis to explore the behavior of the model system (2.1) and illustrate the stability of the equilibria solutions. We numerically solve the model system (2.1) using dde23 ^[14] based on Runge-Kutta methods through MATLAB software and parameters values adopted from Table 1, and the initial population levels were assumed as follows: $S(0) = 10$, and $I(0) = Q(0) = H_m(0) = H_s(0) = 2$.

The numerical results in Figures 2–4 illustrate the dynamical solutions of the model system (2.1) for different values of $\tau$ at the endemic equilibrium point of ${\cal R}_0 = 3.7095$. The results were obtained using the parameter values in Table 1, with $\delta = 0.00000357$, $\alpha = 0.001$, $\sigma_1 = 0.00124$, $\sigma_2 = 0.0182$, $\gamma = 0.001$, and $p = 5.4\times 10^{-10}$ coupled with delay values $\tau = 5$, $\tau = 15$, and $\tau = 25$ for Figures 2, 3, and 4 respectively. To improve the clarity of the results, the solution for all populations were zoomed in. In all cases for certain parameter values and initial population levels, the model system (2.1) exhibits some periodic oscillation. Precisely, we note that the infected population of class $I(t)$ and $Q(t)$ oscillates with a reduced amplitude from the start for a considerable time frame, thereafter the oscillations dies off and converges to the endemic equilibrium point. Similar patterns are observed for compartment $S(t)$ in Figure 2(a), 3(a), and 3(a). We can also note that the intensity and amplitude of oscillation at $\tau = 15$ are high compared to that at $\tau = 5$ and $25.$ In addition, the implication of these results is that the inclusion of the time delay factor destabilizes the endemic equilibrium point for a certain period of time, leading to periodic oscillations which arise due to the existence of Hopf bifurcations. These results agree with the analytical analysis of the global stability for the endemic equilibrium point in Theorem 2.4.

The results in Figure 5 demonstrate the existence of Hopf bifurcations that arise due to the inclusion of the time delay factor $\tau = 15$ in the model system (2.1). The results were obtained by using parameter values in Table (1), with $\delta = 0.00000357$, $\alpha = 0.001$, $\sigma_1 = 0.00124$, $\sigma_2 = 0.0182$, $\gamma = 0.001$, and $p = 5.4\times 10^{-10}$, leading to ${\cal R}_0 = 3.7095$. Overall, these results are in agreement to those depicted in Figures 2–4. Thus, the inclusion of the time delay factor leads to the existence of Hopf bifurcations.

The results in Figure 6 demonstrate the effects of $\sigma_1$ (recovery rate of asymptomatic and quarantined individuals) and $\sigma_2$ (recovery rate of hospitalized individuals) on the dynamics of the disease in the population. The parameters and initial values are fixed and provided in Table 1. Overall, we note that increasing the recovery rate of asymptomatic and quarantined individuals (modeled by parameter $\sigma_1$) reduces the spread of the disease in the population. In particular, when $\sigma_1$ is greater than $50\%$, the disease dies in the population and persists when less than $50\%.$

The results in Figures 7 and 8 demonstrate the solution profile for mild ($H_m$) and severe ($H_s$) hospitalized individuals for different values of $\tau$ at the endemic equilibrium point of ${\cal R}_0 = 3.7095$. The results were obtained using parameters values in Table 1, with $\delta = 0.00000357$, $\alpha = 0.001$, $\sigma_1 = 0.00124$, $\sigma_2 = 0.0182$, $\gamma = 0.001$, and $p = 5.4\times 10^{-10}$ coupled with delay values $\tau = 5$, and $\tau = 15$ for Figures 7 and 8 respectively. Overall, for certain parameter values and initial population levels, the model system (2.1) does not exhibit periodic oscillation. Precisely, we note that the infected population of class $H_m$ and $H_s$ decrease gradually from the start for a considerable time frame, thereafter the solutions converge to the endemic equilibrium point. This has the implication that the inclusion of the time delay factor has less effects on hospitalized individuals. These results agree with the analytical analysis of global stability for the endemic equilibrium point in Theorem 2.4.

Figure 9 demonstrates the convergence of the solution profile of model system (2.1) to the disease-free-equilibrium for ${\cal R}_0 < 1.$ The parameter values used in simulations are in Table (1), with $\beta_1 = 5.33\times 10^{-7}$, $\beta_2 = 0.089$, $\alpha = 0.0001$, $\sigma_1 = 0.00124$, $\sigma_2 = 0.0182$, $\gamma = 0.01$, and $p = 5.4\times 10^{-10}$, leading to ${\cal R}_0 = 0.171$. We observe that the variable for epidemiological classes $S(t)$ and $I(t)$ for $t\leq 50$ all solutions in Figures 9(a) and 9(b) respectively decrease at the beginning and finally attain stability to the disease-free-equilibrium point. In addition, the variable for quarantine individual in Figure 9(a) increase rapidly during the first $50$ days, followed by a gradual decline and stability of solutions at the disease-free equilibrium point. In particular, the disease dies out in the population after 50 days which is in agreement with the analytical results summarized the Theorem 2.3. The results in Figure 10 demonstrate the existence of the bifurcations that arise due to the inclusion of the time delay factor $\tau = 15$ in the model system (2.1). The results were obtained by using parameter values in Table (1), with $\beta_1 = 5.33\times 10^{-7}$, $\beta_2 = 0.089$, $\alpha = 0.0001$, $\sigma_1 = 0.00124$, $\sigma_2 = 0.0182$, $\gamma = 0.01$, and $p = 5.4\times 10^{-10}$, leading to ${\cal R}_0 = 0.171$. We can note that the inclusion of the time delay factor leads to the existence of bifurcations in the model system.

4.   Concluding remarks

In this article, we have developed and analyzed a mathematical model for COVID-19 that incorporates a discrete delay that accounts for the latent period. We compute the basic reproduction number and demonstrate that it is an important threshold quantity for the stability of equilibria. By constructing suitable Lyapunov functionals, it is shown that the model has a globally asymptotically stable infection-free equilibrium whenever the reproduction number is less than unity. Furthermore, whenever the reproduction number is greater than the unity then the model has a unique endemic equilibrium point which is globally asymptotically stable. Numerical simulations are carried out to illustrate the main results. Although quarantine/isolation of an asymptomatic individual is a relatively easy strategy to implement, some studies suggest that quarantine/isolation of asymptomatic, symptomatic and susceptible individuals maybe more effective (see for example ^[26]). The rationale being that by decreasing host density, the number of contacts per unit time between humans is low, thereby reducing disease transmission. In ^[26] it was demonstrated that quarantine/isolation of both asymptomatic and symptomatic individuals only can be effective whenever the number of infected hosts is above a certain critical level ^[26]. We expect to improve this study in the future by developing (COVID-19) model(s) with a time delay that will enable the comparison of the aforementioned aspects. In addition the bifurcation analysis of epidemic models with more compartments and parameters will be more complex and this is a major challenge for the future.

Acknowledgment

All authors are grateful to their respective institutions for their support during the preparation of the manuscript. Paride O. Lolika acknowledges the support from the University of Juba, South Sudan.

Conflict of interest

The authors declare that there are no conflicts of interest.

Supplementary

Appendix A

Proof of Theorem 2.3. We denote by $x_t$ the translation of the solution of models (2.6)-(2.10), that is:

and consider the Lyapunov function:

Then, the time derivative of $V(t)$ along solutions of models (2.6)-(2.10):

Since $S(t)\leq S^0\; (S^0 = \frac{\Lambda}{\delta})$ for $t\geq 0$, we have:

Therefore, $\dot{V}(t) < 0$ holds if ${\cal R}_0 < 1$. Furthermore, $\dot{V}(t) = 0$ if ${\cal R}_0 = 1$. Thus, the largest invariant set of $\dot{V}(t)$ is a singleton such that $S(t) = S^0$, $I(t) = Q(t) = H_m(t) = H_s(t) = 0.$ From the LaSalle invariance principle ^[12], the disease-free equilibrium of models (2.6)-(2.10) denoted by ${\cal E}^0$ is globally asymptotically stable whenever ${\cal R}_0\leq 1$. This completes the proof.

Appendix B

Proof of Theorem 2.4. Let us consider the Lyapunov function:

Here,

The derivatives of ${\cal W}_1(t)$ are given by:

Substituting the appropriate differentials from (2.6)-(2.10), we have:

At endemic equilibrium, we have:

Using the above constants, we have:

The derivatives of ${\cal W}_2^+$ are given by:

Combining the derivatives $\dot{\cal W}_1(t)$ and $\dot{\cal W}_2(t)$, we have:

Since the arithmetic mean is greater than or equal to the geometric mean, we have

and it follows that

for all $S(t) > 0$, because the arithmetic mean is greater than or equal to the geometric mean.

Further, note that a continuous and differentiable function $G(t) = 1-g(t)+\ln g(t)$ is always non positive for any function $g(t) > 0$, and $g(t) = 0$ if and only if $g(t) = 1$. Thus we note that

Hence, it follows that ${\cal W}(t)\leq0$ and consequently, $\dot{\cal W}(t)\leq0$. Moreover, the largest invariant set of $\dot{\cal W}(t) = 0$ is a singleton where $S(t)\equiv S^*$, $I(t)\equiv I^*$, $Q(t)\equiv Q^*$, $H_m(t)\equiv H_m^*$, and $H_s(t)\equiv H_s^*$. Using LaSalle's invariance principle ^[12], we conclude that the endemic equilibrium point ${\cal E}^*$ is globally asymptotically stable if ${\cal R}_0 > 1.$