Mathematical Modelling of the Spatial Distribution of a COVID-19 Outbreak with Vaccination Using Diffusion Equation

Abstract

The formulation of mathematical models using differential equations has become crucial in predicting the evolution of viral diseases in a population in order to take preventive and curative measures. In December 2019, a novel variety of Coronavirus (SARS-CoV-2) was identified in Wuhan, Hubei Province, China, which causes a severe and potentially fatal respiratory syndrome. Since then, it has been declared a pandemic by the World Health Organization and has spread around the globe. A reaction–diffusion system is a mathematical model that describes the evolution of a phenomenon subjected to two processes: a reaction process, in which different substances are transformed, and a diffusion process, which causes their distribution in space. This article provides a mathematical study of the Susceptible, Exposed, Infected, Recovered, and Vaccinated population model of the COVID-19 pandemic using the bias of reaction–diffusion equations. Both local and global asymptotic stability conditions for the equilibria were determined using a Lyapunov function, and the nature of the stability was determined using the Routh–Hurwitz criterion. Furthermore, we consider the conditions for the existence and uniqueness of the model solution and show the spatial distribution of the model compartments when the basic reproduction rate ${\mathcal{R}}_0<1$ and ${\mathcal{R}}_0>1$. Thereafter, we conducted a sensitivity analysis to determine the most sensitive parameters in the proposed model. We demonstrate the model’s effectiveness by performing numerical simulations and investigating the impact of vaccination, together with the significance of spatial distribution parameters in the spread of COVID-19. The findings indicate that reducing contact with an infected person and increasing the proportion of susceptible people who receive high-efficacy vaccination will lessen the burden of COVID-19 in the population. Therefore, we offer to the public health policymakers a better understanding of COVID-19 management.

1. Introduction

The use of differential equations, especially partial differential equations, has become very interesting to predict many evolutionary problems, in particular, the evolution of some biological phenomena. They are essential in fields such as aeronautical simulation, finance, weather prediction, and disease forecasting [1]. One of the most important classes of Partial Differential Equations is the class of reaction–diffusion equations. A reaction–diffusion system is a mathematical model that describes the evolution of a phenomenon subjected to two processes: a reaction process, in which different substances are transformed, and a diffusion process, which causes their distribution in space. Many research papers in epidemiology have proposed a modelling approach using real datasets from affected countries and have identified different characteristics controlled according to various parameters of the epidemic and to the effects of intervention strategies in the different countries concerned, depending on their situation. Ian Cooper and co-authors studied an SIR model of COVID-19 in diverse communities in Asia and North America [2]. In their study, they did not consider migration and death of individuals, so the population size remained constant during the study. In order to compensate for this deficiency, Ref. [3] conducted a study where the population size is not constant and the death rates of all compartments are the same. This led to more general results. Ref. [3] also studied the factors determining the spread of COVID-19 and used statistical modelling to propose strategies to prevent future accelerated viral infection similar to that observed in COVID-19. Furthermore, in [4], the authors presented a novel Susceptible–Infectious–Goneanewsusceptible–Recovered (SIGR) model to study the influence of vaccination at the sub-population level in the spread of COVID-19 pandemic, while in [5], authors proposed Susceptible–Exposed-Infectious–Recovered (SEIR) epidemic model with a convex incidence rate incorporated with a time delay in order to study the influence of delay in the dynamical system.

Other studies (such as [6,7,8,9,10,11,12]) have proposed spatial epidemic models. In these studies, reaction–diffusion equations were used to explain both the temporal and spatial evolution of the spread of diseases. Most of them use a continuous diffusion approach, but some others, such as Mimura’s team’s article, adopt the spatial SEIR model, in which individuals move randomly on a two-dimensional lattice with periodic boundary conditions [7]. Many authors used this class of equations to understand the behaviour of hepatitis C virus. The reaction term is the process of change of individuals involved in the interactions between species in the absence of diffusion, and the diffusion term describes the spatial movement of individuals [13,14,15,16]. In this paper, Ref. [17] studied a model of an SI diffusion reaction in which a pathogen is active in a population with two subgroups: healthy individuals who are susceptible to infection and already infected individuals who can transmit the disease to healthy individuals. To generalize this result to a heterogeneous population, Refs. [18,19] studied a reaction–diffusion Susceptible–Vaccinated–Infected–Recovered (SEIRV) model in a spatially heterogeneous environment with Dirichlet boundary conditions.

Furthermore, in [20], the authors developed a reaction–diffusion epidemic model on human mobility networks to characterize the spatio-temporal propagation of the COVID-19 pandemic, and a novel time-dependent function was incorporated into the model to describe the effects of human interventions. In [21], the authors studied an optimal control problem for a generalized multi-group reaction–diffusion SIR epidemic model, with heterogeneous nonlinear incidence rates, which is an extension of the study of an optimal control problem to a large class of reaction–diffusion multigroup epidemic models. In [22], the research carried out by the authors explored a novel SEIR-A reaction–diffusion COVID-19 epidemic system with direct and aerosol transmission. The effects of three strategies, including vaccination, receiving treatment, and wearing a mask, were evaluated numerically by authors. The findings of their research suggest that the three strategies can effectively control the peak and final scale of infection and shorten the duration of the COVID-19 epidemic.

The present work aims to carry out a mathematical study of a model of infection of COVID-19 by the bias of reaction–diffusion equations. In this current work, we want to make some additions to the work carried out earlier by researchers, by considering the exposed population (not still infectious) and the diffusion phenomena and by showing the existence, the uniqueness, the positiveness, and the boundedness of the solutions of the SEIRV diffusion reaction model considered. A study on the existence and nature of the equilibrium states is carried out according to the basic reproduction number ${\mathcal{R}}_0$ whether it is less than or greater than 1. We were also able to find a Lyapunov function to investigate the global stability of the disease-free and endemic equilibria. For the confirmation of these properties, a numerical simulation is presented at the end of the study with Python software.

The remainder of the article is as follows: in Section 2, we present the mathematical model we formulated alongside parameters and their biological signification. In Section 3, an extensive mathematical analysis is presented to understand the dynamical behaviour of the system. In Section 4, we present a numerical simulation and sensitivity analysis of ${\mathcal{R}}_0$ in order to identify parameters sensitive to the disease spread, and lastly in Section 5 and Section 6, we present some perspective, discussion and key conclusions derived from our research approach to the spatial distribution of COVID-19 outbreak with vaccination using diffusion equation.

2. Mathematical Model

Let us consider the following model:

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\partial S(x,t)}{\partial t}}\hfill & =& d_1\Delta S+\mathsf{\Lambda }-d_{S}S-{\beta }_{1}SI-\lambda S,\hfill \\ {\displaystyle \frac{\partial E(x,t)}{\partial t}}\hfill & =& d_2\Delta E+{\beta }_{1}SI+\beta VI-\sigma E-d_{E}E\hfill \\ {\displaystyle \frac{\partial I(x,t)}{\partial t}}\hfill & =& d_3\Delta I+\sigma E-d_{I}I-\theta I,\hfill \\ {\displaystyle \frac{\partial R(x,t)}{\partial t}}\hfill & =& d_4\Delta R+\gamma V-d_{R}R+\theta I,\hfill \\ {\displaystyle \frac{\partial V(x,t)}{\partial t}}\hfill & =& d_5\Delta V+\lambda S-\beta VI-d_{V}V-\gamma V.\hfill \end{array}\right.$$

(1)

where S, E, I, R, and V are the population of susceptible, exposed, infected, recovered (and definitely immunized), and vaccinated at the position x at time t, respectively. The susceptible population reproduces at a constant rate $\mathsf{\Lambda }$ and is infected at a rate ${\beta }_{1}SI$, where ${\beta }_1$ is the contact rate per day between the susceptible and infected populations and $d_S$, $d_E$, $d_I$, $d_R$, and $d_V$ are, respectively, the mortality rates of susceptible, exposed, infected, recovered, and vaccinated people in the studied region. The parameter $\lambda$ is the vaccination rate against COVID-19 in this population. The recovered population is produced from the infected at a rate $\theta I$. We assume that the population moves in the region $\mathsf{\Omega }$ according to Fick’s law [1], and $d_i$’s being the diffusion coefficients and $\Delta$ the Laplacian operator. In this work, we consider the system (1) with initial conditions as follows

$$S(x,0)={\phi }_1\left(x\right),E(x,0)={\phi }_2\left(x\right),I(x,0)={\phi }_3\left(x\right),V(x,0)={\phi }_4\left(x\right),R(x,0)={\phi }_5\left(x\right)x\in \mathsf{\Omega },$$

(2)

where ${\phi }_i\in C^2\left(\mathsf{\Omega }\right)\cap C\left(\overline{\mathsf{\Omega }}\right)$ and homogeneous Neumann boundary conditions

$$\frac{\partial S}{\partial \nu }=\frac{\partial E}{\partial \nu }=\frac{\partial I}{\partial \nu }=\frac{\partial R}{\partial \nu }=\frac{\partial V}{\partial \nu }=0x\in \partial \mathsf{\Omega },t>0,$$

(3)

where $\mathsf{\Omega }$ is an open bounded subset of ${\mathbb{R}}^n$ with smooth boundary $\partial \mathsf{\Omega }$, $\nu$ being the unit outer normal to $\partial \mathsf{\Omega }$. The interaction graph of the system (1) is presented in Figure 1 while the signification of parameters is presented in

Table 1. Parameters and their signification.

Parameters	Biological Signification of Parameters
$\mathsf{\Lambda }$	Number of incoming susceptible per day
$d_S$	Mortality rate of susceptible per day
${\beta }_1$	Contact between S and I
$\beta$	Contact between V and I
$\lambda$	Vaccination rate
$d_E$	Mortality rate of exposed population per day
$\sigma$	Infection rate
$d_I$	Mortality rate of infected population per day
$d_R$	Mortality rate of recovered population per day
$\theta$	Recovery rate
$d_V$	Mortality rate of vaccinated population per day

.

3. Dynamical Behaviour of the System and Qualitative Analysis

In this section, we study the dynamical behaviour of system (1), such as the existence and uniqueness of positive solutions and existence of equilibria, and their basic reproduction number, local stability, and global stability.

3.1. Existence, Uniqueness, and Positivity

Definition 1. 

Let $(\widehat{S},\widehat{E},\widehat{I},\widehat{V},\widehat{R})$ and $(\tilde{S},\tilde{E},\tilde{I},\tilde{V},\tilde{R})$ in $C(\overline{\mathsf{\Omega }}\times [0,\infty ))\cap C^{1,2}(\mathsf{\Omega }\times [0,\infty ))$ are a pair of upper and lower solution to the problem (1), if $\tilde{S}\le \widehat{S}$, $\tilde{E}\le \widehat{E}$, $\tilde{I}\le \widehat{I}$, $\tilde{V}\le \widehat{V}$, $\tilde{R}\le \widehat{R}$ in $\overline{\mathsf{\Omega }}\times [0,\infty )$, and the following differential inequalities hold:

$$\begin{array}{lll}\frac{\partial \tilde{S}}{\partial t}\le d_1\Delta \tilde{S}+\mathsf{\Lambda }-{\beta }_1\tilde{S}\tilde{I}-(d_S+\lambda )\tilde{S}\hfill & ;& \frac{\partial \widehat{S}}{\partial t}\ge d_1\Delta \widehat{S}+\mathsf{\Lambda }-{\beta }_1\widehat{S}\widehat{I}-(d_S+\lambda )\widehat{S}\hfill \\ \frac{\partial \tilde{E}}{\partial t}\le d_2\Delta \tilde{E}+{\beta }_1\tilde{S}\tilde{I}+\beta \tilde{V}\tilde{I}-(\sigma +d_E)\tilde{E}& ;& \frac{\partial \widehat{E}}{\partial t}\ge d_2\Delta \widehat{E}+{\beta }_1\widehat{S}\widehat{I}+\beta \widehat{V}\widehat{I}-(\sigma +d_E)\widehat{E},\\ \frac{\partial \tilde{I}}{\partial t}\le d_3\Delta \tilde{I}+\sigma \tilde{E}-(d_1+\theta )\tilde{I}& ;& \frac{\partial \widehat{I}}{\partial t}\ge d_3\Delta \widehat{I}+\sigma \widehat{E}-(d_1+\theta )\widehat{I},\\ \frac{\partial \tilde{R}}{\partial t}\le d_4\sigma \tilde{R}+\gamma \tilde{I}-d_R\tilde{R}+\theta \tilde{I}& ;& \frac{\partial \widehat{R}}{\partial t}\ge d_4\Delta \widehat{R}+\gamma \widehat{I}-d_R\widehat{R}+\theta \widehat{I},\\ \frac{\partial \tilde{V}}{\partial t}\le d_5\Delta \tilde{R}+\lambda \tilde{S}-\beta \tilde{V}\tilde{I}-(d_V+\gamma )\tilde{V}& ;& \frac{\partial \widehat{V}}{\partial t}\ge d_5\Delta \widehat{R}+\lambda \widehat{S}-{\beta }_1\widehat{V}\widehat{I}-(d_V+\gamma )\widehat{V}.\end{array}$$

for $(x,t)\in \mathsf{\Omega }\times (0,\infty )$ and

$$\begin{array}{llll}\frac{\partial \tilde{S}}{\partial \nu }\le 0\le \frac{\partial \widehat{S}}{\partial \nu },\frac{\partial \tilde{E}}{\partial t}\le 0\le \frac{\partial \widehat{E}}{\partial \nu },& \frac{\partial \tilde{I}}{\partial \nu }\le 0\le \frac{\partial \widehat{I}}{\partial \nu },\frac{\partial \tilde{V}}{\partial \nu }\le 0\le \frac{\partial \widehat{V}}{\partial t},\hfill & \frac{\partial \tilde{R}}{\partial \nu }\le 0\le \frac{\partial \widehat{R}}{\partial \nu },\mathit{for}(x,t)\in \partial \mathsf{\Omega }\times (0,\infty )\hfill & \\ \tilde{S}(x,t)\le {\phi }_1(x,t)\le \widehat{S}(x,t),& \tilde{E}(x,t)\le {\phi }_2(x,t)\le \widehat{E}(x,t),& \tilde{I}(x,t)\le {\phi }_3(x,t)\le \widehat{I}(x,t),\\ \tilde{V}(x,t)\le {\phi }_4(x,t)\le \widehat{V}(x,t),& \tilde{R}(x,t)\le {\phi }_5(x,t)\le \widehat{R}(x,t),& \mathit{for}(x,t)\in \partial \mathsf{\Omega }\times (0,\infty )\end{array}$$

It is easy to see that $\underline{\mathbf{0}}$ $=(0,0,0,0,0)$ and $\mathbf{K}$ $=(K_1,K_2,K_3,K_4,K_5)$ are a pair of coupled lower-upper solutions to problem (1), where

$$\begin{array}{ccc}{K}_1\hfill & =\hfill & max\{\frac{\mathsf{\Lambda }}{d},\parallel {\phi }_1{\parallel }_{C(\overline{\mathsf{\Omega }},\mathbb{R})}\}\hfill \\ K_2\hfill & =\hfill & max\{\frac{\mathsf{\Lambda }}{d},\parallel {\phi }_2{\parallel }_{C(\overline{\mathsf{\Omega }},\mathbb{R})}\}\hfill \\ K_3\hfill & =\hfill & max\{\frac{\sigma \mathsf{\Lambda }}{d^2},\parallel {\phi }_3{\parallel }_{C(\overline{\mathsf{\Omega }},\mathbb{R})}\}\hfill \\ K_4\hfill & =\hfill & max\{\frac{d\gamma \mathsf{\Lambda }+\gamma \mathsf{\Lambda }\theta }{d^3},\parallel {\phi }_4{\parallel }_{C(\overline{\mathsf{\Omega }},\mathbb{R})}\}\hfill \\ K_5\hfill & =\hfill & max\{\frac{\lambda \mathsf{\Lambda }}{d^2},\parallel {\phi }_5{\parallel }_{C(\overline{\mathsf{\Omega }},\mathbb{R})}\},\hfill \end{array}$$

and $d=min\{d_S+\lambda ;\sigma +d_E;d_I+\theta ;d_R\}$. Using the following lemma provided by Redinger [13], we obtain the existence and uniqueness of the solution.

Lemma 1. 

Let $\widehat{U}$ and $\tilde{U}$ be a pair of upper and lower solutions for problem (1) and suppose that the initial functions ${\phi }_i$$(i=1,2\dots ,5)$ are continuous in $\overline{\mathsf{\Omega }}$. Then problem (1) has exactly one regular solution $U(x,t)=\left(S\right(x,t),E(x,t),I(x,t),V(x,t),R(x,t\left)\right)$ satisfying $\tilde{U}\le U\le \widehat{U}$ in $\overline{\mathsf{\Omega }}\times [0,\infty )$.

Hence, $0\le S(x,t)\le K_1$, $0\le E(x,t)\le K_2$, $0\le I(x,t)\le K_3$, $0\le V(x,t)\le K_4$, $0\le R(x,t)\le K_5$. Furthermore, also, by the maximum principle, if ${\phi }_i\left(x\right)\ne 0$, we have $S(x,t)>0,E(x,t)>0,I(x,t)>0,V(x,t)>0,R(x,t)>0$ for all $t>0,x\in \overline{\mathsf{\Omega }}$.

3.2. Equilibria and Basic Reproduction Number

3.2.1. Disease-Free Equilibrium and Basic Reproduction Number

It is easy to verify that system (1) always has a disease-free equilibrium given by $(S_0,0,0,0,V_0)$ where

$$\begin{array}{ccc}{S}_0\hfill & =\hfill & \frac{\mathsf{\Lambda }}{d_S+\lambda }\hfill \\ V_0\hfill & =\hfill & \frac{\lambda }{d_V+\gamma }\frac{\mathsf{\Lambda }}{d_S+\lambda }.\hfill \end{array}$$

In order to find the basic reproduction number ${\mathcal{R}}_0$ of the system (1), we obtain the following linear system at disease-free equilibrium ${\mathcal{E}}_0$:

$$\left\{\begin{array}{cc}{\displaystyle \frac{\partial E(x,t)}{\partial t}}=d_2\Delta E+{\beta }_{1}SI+\beta VI-\sigma E-d_{E}E,\hfill & x\in \mathsf{\Omega },t>0\hfill \\ {\displaystyle \frac{\partial I(x,t)}{\partial t}}=d_3\Delta I+\sigma E-d_{I}I-\theta I,\hfill & x\in \mathsf{\Omega },t>0\hfill \\ {\displaystyle \frac{\partial R(x,t)}{\partial t}}=d_4\Delta R+\gamma V-d_{R}R+\theta I,\hfill & x\in \mathsf{\Omega },t>0\hfill \\ {\displaystyle \frac{\partial E}{\partial \nu }}={\displaystyle \frac{\partial I}{\partial \nu }}={\displaystyle \frac{\partial R}{\partial \nu }}=0\hfill & x\in \partial \mathsf{\Omega },t>0.\hfill \end{array}\right.$$

(4)

By following the idea of Wendi Wang [23], consider the vectors $\mathcal{F}$ and $\mathcal{V}$ given by

$$\mathcal{F}=\left(\begin{array}{c}{\beta }_{1}SI+\beta VI\\ \sigma E\\ \gamma V+\theta I\end{array}\right),\mathcal{V}=\left(\begin{array}{c}\sigma E+d_{E}E\\ d_{I}I+\theta I\\ d_{R}R\end{array}\right).$$

where $\mathcal{F}$ is the input rate vector and $\mathcal{V}$ the transfer rate vector for each compartment. In order to find the next generation matrix $\mathbb{K}$, we have

$$J\mathcal{F}\left({\mathcal{E}}_0\right)=\left(\begin{array}{ccc}0& \frac{\mathsf{\Lambda }}{d_S+\lambda }\left({\beta }_1+\frac{\beta \lambda }{d_V+\gamma }\right)& 0\\ \sigma & 0& 0\\ 0& \theta & 0\end{array}\right),J\mathcal{V}\left({\mathcal{E}}_0\right)=\left(\begin{array}{cc}\sigma +d_E& 0\\ 0& d_I+\theta & 0\\ 0& 0& d_R\end{array}\right).$$

Therefore, $\mathbb{K}$ is given by

$$\mathbb{K}=J\mathcal{F}\left({\mathcal{E}}_0\right)J{\mathcal{V}}^{-1}\left({\mathcal{E}}_0\right)=\left(\begin{array}{ccc}0& \frac{{\beta }_1{S}_0+\beta V_0}{d_I+\theta }& 0\\ \frac{\sigma }{\sigma +d_E}& 0& 0\\ 0& 0& 0\end{array}\right).$$

Finally, the basic reproduction number is given by

$${\mathcal{R}}_0=\rho \left(\mathbb{K}\right)=\frac{\sigma \mathsf{\Lambda }({\beta }_1(d_V+\gamma )+\beta \lambda )}{(d_S+\lambda )(d_V+\gamma )(\sigma +d_E)(d_I+\theta )}.$$

(5)

3.2.2. Existence and Uniqueness of the Endemic Equilibrium

For the second equilibrium, we are looking for ${\mathcal{E}}^{*}=(S^{*},E^{*},I^{*},R^{*},V^{*})$ all different from 0 such that

$$\begin{array}{ccc}\hfill \mathsf{\Lambda }-d_{S}S-{\beta }_{1}SI-\lambda S& =\hfill & \hfill 0\\ \hfill {\beta }_{1}SI+\beta VI-\sigma E-d_{E}E& =\hfill & \hfill 0\\ \hfill \sigma E-(d_I+\theta )I& =\hfill & \hfill 0\\ \hfill \gamma V-d_{R}R+\theta I& =\hfill & \hfill 0\\ \hfill \lambda S-\beta VI-d_{V}V-\gamma V& =\hfill & \hfill 0.\end{array}$$

(6)

From these equations, we come out with

$$\begin{array}{ccc}\hfill S=f_S\left(I\right)& =& \frac{\mathsf{\Lambda }}{d_S+\lambda +{\beta }_{1}I}\hfill \\ \hfill E=f_E\left(I\right)& =& \frac{I(d_I+\theta )}{\sigma }\hfill \\ \hfill V=f_V\left(I\right)& =& \frac{\lambda f_S\left(I\right)}{\beta I+d_V+\gamma }\hfill \\ \hfill R=f_R\left(I\right)& =& \frac{\gamma f_V\left(I\right)+\theta I}{d_R}.\hfill \end{array}$$

(7)

By using (6) and (7), we will prove the existence of solution of equation

$$G\left(I\right)={\beta }_{1}I{f}_S\left(I\right)+\beta I{f}_V\left(I\right)-(\sigma +d_E)f_E\left(I\right)=0.$$

(8)

It is easy to see that $G\left(0\right)=0$ and $\underset{I\to +\infty }{lim}G\left(I\right)=-\infty$. Therefore, the previous equation has a solution if $G^{\prime }\left(0\right)>0$. However, we have

$$\begin{array}{ccc}\hfill G^{\prime }\left(0\right)& =& {\beta }_1{f}_S\left(0\right)+\beta h\left(0\right)-(\sigma +d_E)g^{\prime }\left(0\right)\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& =& \frac{{\beta }_1\mathsf{\Lambda }}{d_S+\lambda }+\frac{\beta \lambda \mathsf{\Lambda }}{(d_V+\gamma )(d_S+\lambda )}-\frac{(\sigma +d_E)(d_I+\theta )}{\sigma }.\hfill \end{array}$$

(10)

It follows that $G^{\prime }\left(0\right)>0$ if and only if

$$0<\frac{{\beta }_1\mathsf{\Lambda }}{d_S+\lambda }+\frac{\beta \lambda \mathsf{\Lambda }}{(d_V+\gamma )(d_S+\lambda )}-\frac{(\sigma +d_E)(d_I+\theta )}{\sigma }.$$

This means

$$\begin{array}{ccc}\hfill \frac{(\sigma +d_E)(d_I+\theta )}{\sigma }& <\hfill & \frac{{\beta }_1\mathsf{\Lambda }}{d_S+\lambda }+\frac{\beta \lambda \mathsf{\Lambda }}{(d_V+\gamma )(d_S+\lambda )}\hfill \\ \hfill 1& <\hfill & \frac{\mathsf{\Lambda }\sigma \left({\beta }_1(d_V+\gamma )+\beta \lambda \right)}{(\sigma +d_E)(d_I+\theta )(d_S+\lambda )(d_V+\gamma )}\equiv {\mathcal{R}}_0.\hfill \end{array}$$

To prove uniqueness, we are looking for a solution $I\ne 0$ which verifies (8). Since $I\ne 0$, it can divided through, and then using (7), we obtain

$$h\left(I\right)=\frac{{\beta }_1\mathsf{\Lambda }}{d_S+\lambda +{\beta }_{1}I}+\frac{\beta \lambda \mathsf{\Lambda }}{(d_S+\lambda +{\beta }_{1}I)(\beta I+d_V+\gamma )}-\frac{(\sigma +d_E)(d_I+\theta )}{\sigma },$$

and it follows that

$$h^{\prime }\left(I\right)=-\frac{{\beta }_1^2\mathsf{\Lambda }}{{(d_S+\lambda +{\beta }_{1}I)}^2}-\frac{\beta \mathsf{\Lambda }\lambda \left({\beta }_1(\beta I+d_V+\gamma )+\beta (d_S+\lambda +{\beta }_{1}I)\right)}{{(d_S+\lambda +{\beta }_{1}I)}^2{(\beta I+d_V+\gamma )}^2}<0.$$

(11)

Finally, the function G is injective. Therefore, if $I_1$ and $I_2$ are two solutions of (8), then $G\left(I_1\right)=G\left(I_2\right)=0$ and by injectivity of G, we have $I_1=I_2$. This allows us to conclude the uniqueness of the endemic equilibrium. We come out with this following theorem.

Theorem 1. 

If ${\mathcal{R}}_0>1$, then the infected equilibrium ${\mathcal{E}}^{*}=(S^{*},E^{*},I^{*},R^{*},V^{*})$ exists and is unique.

3.3. Stability of Equilibrium

Let $0={\mu }_0<{\mu }_i<{\mu }_{i+1}$ for $i=1,2,\dots$ be the eigenvalues of $-\Delta$ on $\mathsf{\Omega }$ with homogeneous Neumann boundary conditions, $\mathcal{E}\left({\mu }_i\right)$ is the space of eigenfunctions corresponding to ${\mu }_i$’s and $\{{\varphi }_{ij}:j=1,2,\dots ,dim\mathcal{E}\left({\mu }_i\right)\}$ an orthogonal basis of $\mathcal{E}\left({\mu }_i\right)$. Then, $\mathbb{X}={\left[C^1\left(\overline{\mathsf{\Omega }}\right)\right]}^3$ can be decomposed as

$$\mathbb{X}=\underset{i=1}{\overset{\infty }{⨁}}{\mathbb{X}}_i,{\mathbb{X}}_i=\underset{i=1}{\overset{dim\mathcal{E}\left({\mu }_i\right)}{⨁}}{\mathbb{X}}_{ij}$$

where ${\mathbb{X}}_{ij}=\{c{\varphi }_{ij}:c\in {\mathbb{R}}^3\}$. Then we can prove the local stability of equilibrium as in [2,24].

Theorem 2. 

If ${\mathcal{R}}_0<1$, then the disease-free equilibrium ${\mathcal{E}}_0$ of system (1) is locally asymptotically stable.

Proof. 

The linearisation of system (1) at ${\mathcal{E}}_0$ can be expressed by

$$\frac{\partial Z(x,t)}{\partial t}=\mathcal{D}\Delta Z(x,t)+\mathcal{A}Z(x,t),$$

where $\Delta$ is the Jacobian matrix of system (1) and $\mathcal{A}$ is the residual non-linear operator $Z=(S,E,I,R,V)$, $\mathcal{D}=diag(d_1,d_2,d_3,d_4,d_5)$, and if J is the Jacobian matrix of system (1), we have:

$$-J+\mathcal{D}\Delta =\left(\begin{array}{ccccc}-(d_S+\lambda )-\beta I& 0& -{\beta }_{1}S& 0& 0\\ {\beta }_I& -(\sigma +d_E)& {\beta }_{1}S+\beta V& 0& \beta I\\ 0& \sigma & -(d_I+\theta )& 0& 0\\ 0& 0& \theta & -d_R& \gamma \\ \lambda & 0& -\beta V& 0& -\beta I-(d_I+\gamma ).\end{array}\right)$$

(12)

Therefore, the characteristic equation of −J at ${\mathcal{E}}_0$ is

$$(-d_R-{\mu }_i{d}_4-x)(-(d_S+\lambda )-{\mu }_i{d}_1-x)(-(d_I+\gamma )-{\mu }_i{d}_3-x)P\left(x\right),$$

(13)

where

$$P\left(x\right)=x^2+(\sigma +d_E+d_I+\theta +{\mu }_i(d_2+d_5))x+(\sigma +d_E)(d_I+\theta )-\sigma ({\beta }_1{S}_0+\beta V_0).$$

It is obvious that (13) has as eigenvalues

$$x_1=-d_R-{\mu }_i{d}_4,x_2=-(d_S+\lambda )-{\mu }_i{d}_1,x_3=-(d_I+\gamma )-{\mu }_i{d}_3.$$

The two other ones are $x_4$ and $x_5$ such that

$$\begin{array}{ccc}\hfill x_4+x_5& =\hfill & -(\sigma +d_E+d_I+\theta +{\mu }_i(d_2+d_5))\hfill \\ \hfill x_4{x}_5& =\hfill & (\sigma +d_E)(d_I+\theta )-\sigma ({\beta }_1{S}_0+\beta V_0).\hfill \end{array}$$

They are negative if and only if $x_4{x}_5>0$; that is,

$$0<(\sigma +d_E)(d_I+\theta )-\sigma ({\beta }_1{S}_0+\beta V_0).$$

(14)

By using the values of $S_0$ and $V_0$ given by (6), Equation (14) is equivalent to

$$\sigma ({\beta }_1{S}_0+\beta V_0)<(\sigma +d_E)(d_I+\theta ).$$

This implies

$$\frac{\sigma \mathsf{\Lambda }}{d_S+\lambda }\left(\frac{\beta \lambda }{d_V+\gamma }+{\beta }_1\right)<(\sigma +d_E)(d_I+\theta ),$$

and it follows that $\frac{\sigma \mathsf{\Lambda }}{d_S+\lambda }\frac{{\beta }_1(d_V+\gamma )+\beta \lambda }{d_V+\gamma }<(\sigma +d_E)(d_I+\theta )$. Finally, we come out with

$${\mathcal{R}}_0\equiv \frac{\sigma \mathsf{\Lambda }({\beta }_1(d_V+\gamma )+\beta \lambda )}{(d_S+\lambda )(d_V+\gamma )(\sigma +d_E)(d_I+\theta )}<1,$$

and the result follows. □

Theorem 3. 

If ${\mathcal{R}}_0>1$, then the endemic equilibrium ${\mathcal{E}}^{*}$ of system (1) is locally asymptotically stable.

Proof. 

The linearization of system (1) at ${\mathcal{E}}^{*}$ can be expressed by

$$\frac{\partial Z(x,t)}{\partial t}=\mathcal{D}\Delta Z(x,t)+\mathcal{B}Z(x,t),$$

where

$$-J+\mathcal{D}\Delta =\left(\begin{array}{ccccc}-\left(d_S+\lambda +{\beta }_1{I}^{*}\right)& 0& -{\beta }_1{S}^{*}& 0& 0\\ {\beta }_1{I}^{*}& -(\sigma +d_E)& {\beta }_1{S}^{*}+\beta V^{*}& 0& \beta I^{*}\\ 0& \sigma & -(d_I+\theta )& 0& 0\\ 0& 0& \theta & -d_R& \gamma \\ \lambda & 0& -\beta V^{*}& 0& -\left(\beta I^{*}+d_I^{*}+\gamma \right)\end{array}\right),$$

and $\mathcal{D}=(d_1,d_2,d_3,d_4,d_5)$.

Thus, the characteristic equation at ${\mathcal{E}}^{*}$ is

$$\left(-d_R-{\mu }_i{d}_4-x\right)\left(x^4+{\lambda }_3{x}^3+{\lambda }_2{x}^2+{\lambda }_{1}x+{\lambda }_0\right)=0,$$

(15)

where

$$\begin{array}{ccc}\hfill {\lambda }_3& =\hfill & (\left(d_S+\lambda +{\beta }_1{I}^{*}\right)+{\mu }_i{d}_1+\left(\beta I^{*}+d_V+\gamma \right)+{\mu }_i{d}_5+(\sigma +d_E)+{\mu }_i{d}_2+(d_I+\theta )+{\mu }_i{d}_3)>0\hfill \\ \hfill {\lambda }_2& =\hfill & -\left(d_S+\lambda +{\beta }_1{I}^{*}\right)+{\mu }_i{d}_1\left(\beta I^{*}+d_V+\gamma +{\mu }_i{d}_5\right)\sigma {\beta }_1{S}^{*}+\beta V^{*}(\sigma +d_E+{\mu }_i{d}_2)(d_I+\theta +{\mu }_i{d}_3)>0\hfill \\ \hfill {\lambda }_1& =\hfill & \left(\sigma {\beta }_1{S}^{*}+\beta V^{*}\right)\left(1+\left(d_S+\lambda +{\beta }_1{I}^{*}\right)+{\mu }_i{d}_1+\left(\beta I^{*}+d_V+\gamma \right)+{\mu }_i{d}_5\right)+\beta {\beta }_1\sigma I^{*}(S^{*}-V^{*})\hfill \\ \hfill {\lambda }_0& =\hfill & \left(d_S+\lambda +{\beta }_1{I}^{*}+{\mu }_i{d}_1\right)\left(\beta I^{*}+d_V+\gamma +{\mu }_i{d}_5\right)(\sigma +d_E+{\mu }_i{d}_2)(d_I+\theta +{\mu }_i{d}_3)\hfill \\ \hfill & \hfill & -\beta {\beta }_1\sigma I^{*}\left(\lambda S^{*}-\left(\beta I^{*}+d_V+\gamma +{\mu }_i{d}_5\right)S^{*}-V^{*}\left(d_S+\lambda +{\beta }_1{I}^{*}+{\mu }_i{d}_1\right)+\right)\hfill \\ \hfill & \hfill & -\sigma \left(\beta I^{*}+d_V+\gamma +{\mu }_i{d}_5\right){\beta }_1{S}^{*}+\beta V^{*}(A+1).\hfill \end{array}$$

Hence, ${\lambda }_0,{\lambda }_2,{\lambda }_3>0$ and ${\lambda }_1\left({\lambda }_3{\lambda }_2-{\lambda }_1\right)>{\lambda }_3^2{\lambda }_0$ whenever ${\mathcal{R}}_0>1$. Then, using the Routh–Hurwitz criterion, we claim that all eigenvalues of (15) have negative real parts. Thus, the endemic equilibrium ${\mathcal{E}}^{*}$ of system (1) is locally asymptotically stable when ${\mathcal{R}}_0>1$. This completes the proof. □

3.4. Global Stability of the Disease-Free and Endemic Equilibria by Means of a Lyapunov Function

In this section, we investigate the global stability of the disease-free equilibrium ${\mathcal{E}}_0$ for system (1). We consider a Lyapunov functional based on the Volterra function $\mathsf{\Phi }\left(x\right)=x-1-lnx$. Clearly, $\mathsf{\Phi }\left(x\right)\ge 0$ for all $x>0$ and the equality holds if and only if $x=1$. In the presence of diffusion, the aim is to show that every solution of system (1) with a positive initial value that is different from the equilibrium point will converge to the equilibrium.

Theorem 4. 

If ${\mathcal{R}}_0<1$, then the disease-free equilibrium ${\mathcal{E}}_0$ of system (1) is globally asymptotically stable in the feasible region Γ. If ${\mathcal{R}}_0>1$, then ${\mathcal{E}}_0$ is unstable.

Proof. 

Define a Lyapunov function

$$\mathcal{L}\left(t\right)={\int }_{\mathsf{\Omega }}{\mathcal{L}}_1(x,t)dx,$$

where

$${\mathcal{L}}_1(x,t)=aE(x,t)+bI(x,t)$$

and $a,b$ are positive constants to be determined later. Then, along with the solutions of system (1), we have

$$\begin{array}{ccc}\hfill \frac{\partial {\mathcal{L}}_1}{\partial t}& =\hfill & a\frac{\partial E}{\partial t}+b\frac{\partial I}{\partial t}\hfill \\ \hfill & =\hfill & a({\beta }_{1}SI+\beta VI-\sigma E-d_{E}E)+b(\sigma E-(d_I+\theta )I)+a{d}_2\Delta E+b{d}_3\Delta I\hfill \\ \hfill & \le \hfill & a({\beta }_1{S}_{0}I+\beta V_{0}I-\sigma E-d_{E}E)+b(\sigma E-(d_I+\theta )I)+a{d}_2\Delta E+b{d}_3\Delta I\hfill \\ \hfill & =\hfill & I\left(a{\beta }_1{S}_0+a\beta V_0-b(d_I+\theta )\right)+\left(b\sigma -a(\sigma +d_E)\right)E+a{d}_2\Delta E+b{d}_3\Delta I\hfill \end{array}$$

By choosing $a=\sigma$ and $b=\sigma +d_E$, we obtain

$$\frac{\partial {\mathcal{L}}_1}{\partial t}\le (d_I+\theta )(\sigma +d_E)({\mathcal{R}}_0-1)I+a{d}_2\Delta E+b{d}_3\Delta I.$$

Using Green’s formula and the Neumann boundary conditions in (3), we obtain

$${\int }_{\mathsf{\Omega }}\Delta Idx={\int }_{\mathsf{\Omega }}\Delta Edx=0.$$

Using the above conditions, we have

$$\frac{\partial \mathcal{L}}{\partial t}\le {\int }_{\mathsf{\Omega }}\left((d_I+\theta )(\sigma +d_E)({\mathcal{R}}_0-1)I\right)dx.$$

(16)

Therefore, $\frac{\partial \mathcal{L}}{\partial t}\le 0$ whenever ${\mathcal{R}}_0<1$. Furthermore, $\frac{\partial \mathcal{L}}{\partial t}=0$ if and only if $I=0$. It follows that the largest invariant subset $\{(S,E,I,V,R)\mathrm{such}\mathrm{that}\dot{\mathcal{L}}=0\}$ when ${\mathcal{R}}_0<1$ is reduced to the singleton ${\mathcal{E}}_0$. By LaSalle’s Invariance Principle [25], the infection-free equilibrium of system (1) is globally asymptotically stable when ${\mathcal{R}}_0<1$ and if ${\mathcal{R}}_0>1$, then ${\mathcal{E}}_0$ is unstable. □

Theorem 5. 

Consider a Lyapunov function

$$\mathcal{H}\left(t\right)={\int }_{\mathsf{\Omega }}{\mathcal{H}}_1(x,t)dx,$$

with

$${\mathcal{H}}_1(x,t)=S^{*}\mathsf{\Phi }\left(\frac{S}{S^{*}}\right)+I^{*}\mathsf{\Phi }\left(\frac{I}{I^{*}}\right)+R^{*}\mathsf{\Phi }\left(\frac{R}{R^{*}}\right).$$

where

$$\mathsf{\Phi }\left(x\right)=x-1-lnx$$

Then, $\mathcal{H}$ is non-negative and is strictly minimized at the unique equilibrium $(S^{*},E^{*},I^{*},R^{*},V^{*})$; i.e., it is a valid Lyapunov function. Hence, ${\mathcal{E}}^{*}=(S^{*},E^{*},I^{*},R^{*},V^{*})$ is globally asymptotically stable.

Proof. 

According to (1), we have

$$\begin{array}{ccc}\hfill \frac{\partial \mathcal{H}(x,t)}{\partial t}& =\hfill & \left(1-\frac{S^{*}}{S}\right)\frac{\partial S}{\partial t}+\left(1-\frac{I^{*}}{I}\right)\frac{\partial I}{\partial t}+\left(1-\frac{R^{*}}{R}\right)\frac{\partial R}{\partial t}\hfill \\ \hfill & =\hfill & \left(1-\frac{S^{*}}{S}\right)\left(d_1\Delta S+\mathsf{\Lambda }-d_{S}S-{\beta }_{1}SI-\lambda S\right)+\left(1-\frac{I^{*}}{I}\right)\left(d_3\Delta I+\sigma E-d_{I}I-\theta I\right)\hfill \\ \hfill & +\hfill & \left(1-\frac{R^{*}}{R}\right)\left(d_4\Delta R+\gamma V-d_{R}R+\theta I\right)\hfill \\ \hfill & =\hfill & -\frac{{(S^{*}-S)}^2(d_S+\lambda )}{S}+\left(1-\frac{S^{*}}{S}\right)d_1\Delta S+{\beta }_1{S}^{*}{I}^{*}\left(1-\frac{SI}{S^{*}{I}^{*}}-\frac{S^{*}}{S}+\frac{I}{I^{*}}\right)\hfill \\ \hfill & +\hfill & \sigma \frac{(I-I^{*})(E-E^{*})}{I}-\frac{(d_I+\theta ){(I-I^{*})}^2}{I}+\left(1-\frac{I^{*}}{I}\right)d_3\Delta I+\left(1-\frac{R^{*}}{R}\right)\frac{\partial R}{\partial t}\hfill \\ \hfill & =\hfill & -\frac{{(S^{*}-S)}^2(d_S+\lambda )}{S}+\left(1-\frac{S^{*}}{S}\right)d_1\Delta S-{\beta }_1{S}^{*}{I}^{*}\left(\mathsf{\Phi }\left(\frac{SI}{S^{*}{I}^{*}}\right)+\mathsf{\Phi }\left(\frac{S^{*}}{S}\right)-\mathsf{\Phi }\left(\frac{I}{I^{*}}\right)\right)\hfill \\ \hfill & +\hfill & \left(1-\frac{I^{*}}{I}\right)d_3\Delta I-\sigma E^{*}{I}^{*}\left(\mathsf{\Phi }\left(\frac{I}{I^{*}}\right)+\mathsf{\Phi }\left(\frac{E}{E^{*}}\right)-\mathsf{\Phi }\left(\frac{EI}{E^{*}{I}^{*}}\right)\right)\hfill \\ \hfill & -\hfill & \frac{d_R}{R}{(R^{*}-R)}^2-R^{*}{I}^{*}\left(\mathsf{\Phi }\left(\frac{I}{I^{*}}\right)+\mathsf{\Phi }\left(\frac{R}{R^{*}}\right)-\mathsf{\Phi }\left(\frac{RI}{R^{*}{I}^{*}}\right)\right)+\left(1-\frac{R^{*}}{R}\right)d_4\Delta R.\hfill \end{array}$$

Using Green’s formula and the Neumann boundary conditions in (3), we obtain

$${\int }_{\mathsf{\Omega }}\left(1-\frac{S^{*}}{S}\right)d_1\Delta Sdx=-d_1{\int }_{\mathsf{\Omega }}\nabla \left(1-\frac{S^{*}}{S}\right)\nabla Sdx=-d_1{\int }_{\mathsf{\Omega }}\frac{S^{*}}{S^2}{\left|\nabla S\right|}^{2}dx\le 0,$$

$${\int }_{\mathsf{\Omega }}\left(1-\frac{V^{*}}{V}\right)d_1\Delta Vdx=-d_1{\int }_{\mathsf{\Omega }}\nabla \left(1-\frac{V^{*}}{V}\right)\nabla Vdx=-d_1{\int }_{\mathsf{\Omega }}\frac{V^{*}}{V^2}{\left|\nabla V\right|}^{2}dx\le 0,$$

$${\int }_{\mathsf{\Omega }}\left(1-\frac{R^{*}}{R}\right)d_1\Delta Rdx=-d_1{\int }_{\mathsf{\Omega }}\nabla \left(1-\frac{R^{*}}{R}\right)\nabla Rdx=-d_1{\int }_{\mathsf{\Omega }}\frac{R^{*}}{R^2}{\left|\nabla R\right|}^{2}dx\le 0.$$

Using the above conditions, we conclude that

$$\frac{d\mathcal{H}\left(t\right)}{dt}\le {\int }_{\mathsf{\Omega }}\left[-\frac{{(S^{*}-S)}^2(d_S+\lambda )}{S}-\frac{d_R}{R}{(R^{*}-R)}^2-\frac{(d_I+\theta ){(I-I^{*})}^2}{I}\right].$$

(17)

Furthermore, we have $\frac{d\mathcal{H}\left(t\right)}{dt}=0$ only at steady-state ${\mathcal{E}}^{*}=(S^{*},E^{*},I^{*},R^{*},V^{*})$. Therefore, by Lyapunov’s direct method, the steady state solution ${\mathcal{E}}^{*}$ is globally asymptotically stable. □

4. Numerical Simulation and Sensitivity Analysis of ${\mathcal{R}}_{\mathbf{0}}$

4.1. Numerical Simulation

4.1.1. Experiment 1: Numerical Simulation When ${\mathcal{R}}_0<1$

Here, $n=1$ (spatial dimension). We choose

$$(d_1,d_2,d_3,d_4,d_5)=(1000,900,20,900,1200){\mathrm{m}}^2{\mathrm{day}}^{-1}.$$

(18)

With the values in

Table 2. Values of the parameters when ${\mathcal{R}}_0<1$.

Parameters	Values	References
$\mathsf{\Lambda }$	350	Estimated
$d_S$	$0.03324588$	[26]
${\beta }_1$	$0.0000051$	[26]
$\lambda$	$0.000091$	[26]
$\beta$	$0.00001611$	[26]
$\sigma$	$17/100$	[26]
$d_E$	$0.003324588$	[26]
$d_I$	$0.06184588$	[26]
$\theta$	$0.1109289$	[26]
$d_R$	$0.03324588$	[26]
$d_V$	$0.003324588$	[26]
$\gamma$	$0.15$	Estimated

and simulation, we obtain ${\mathcal{R}}_0=2.{94}^{-6}<1$ and at the disease-free equilibrium, we show $S_0$ to be 1800 and $V_0$ to be 8000. Figure 2 displays the solution curves of the model where they tend to the stability of the disease-free equilibrium point with different initial histories. Hence, numerical simulations of Experiment 1 confirm the qualitative results.

4.1.2. Experiment 2: Numerical Simulation When ${\mathcal{R}}_0>1$

Here, $n=1$ (spatial dimension), and we use the same diffusion coefficient as presented in experiment 1.

With the values in

Table 3. Values of the parameters when ${\mathcal{R}}_0>1$.

Parameters	Values	References
$\mathsf{\Lambda }$	15,000	Estimated
$d_S$	$0.003324588$	[26]
${\beta }_1$	$0.0000051$	[26]
$\lambda$	$17/100$	Estimated
$\beta$	$0.00000091$	[27]
$\sigma$	$0.00001611$	Estimated
$d_E$	$0.003324588$	[26]
$d_I$	$0.06184588$	[26]
$\theta$	$0.1109289$	[26]
$d_R$	$0.03324588$	[26]
$d_V$	$0.003324588$	[26]
$\gamma$	$0.95$	[26]

, the endemic equilibrium value is ${\mathcal{E}}^{*}=(4000,1.25,6000,$$3500,{10}^5)$, and ${\mathcal{R}}_0=2.05>1$. The result is presented in Figure 3 and is consistent with the stability of endemic equilibrium presented in the qualitative analysis.

4.1.3. Experiment 3: Numerical Simulation When ${\mathcal{R}}_0<1$ and $\lambda =0$

In this section, all parameters used are the same as the ones used for experiment 1, except that in this case $\lambda =0$. The result of the simulation is presented in Figure 4.

4.1.4. Experiment 4: Numerical Simulation When ${\mathcal{R}}_0>1$ and $\lambda =0$

In this section, all parameters used are the same as those used in experiment 2, except that in this case $\lambda =0$. The result of the simulation is presented in Figure 5.

4.2. Local Sensitivity of ${\mathcal{R}}_{\mathbf{0}}$

To know the effect of the parameters in system (1) on the control reproduction number ${\mathcal{R}}_0$, we perform sensitivity analysis on ${\mathcal{R}}_0$. This analysis is investigated analytically by computing $\frac{\partial {\mathcal{R}}_0}{\partial \mathcal{Z}}$, where $\mathcal{Z}=(\sigma ,\mathsf{\Lambda },{\beta }_1,\beta ,d_S,\lambda ,d_V,\gamma ,d_I,\theta ,d_E).$ The sensitivity of ${\mathcal{R}}_0$ to each parameter is as follows:

$$\begin{array}{ccc}\frac{\partial {\mathcal{R}}_0}{\partial \sigma }\hfill & =\hfill & \frac{\mathsf{\Lambda }\left({\beta }_1(d_V+\gamma )+\beta \lambda \right)}{(d_S+\lambda )(d_V+\gamma )(d_I+\theta )}\frac{d_E}{{(\sigma +d_E)}^2}>0\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial \mathsf{\Lambda }}\hfill & =\hfill & \frac{\sigma \left({\beta }_1(d_V+\gamma )+\beta \lambda \right)}{(d_S+\lambda )(d_V+\gamma )(\sigma +d_E)(d_I+\theta )}>0\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial {\beta }_1}\hfill & =\hfill & \frac{\sigma \mathsf{\Lambda }(d_V+\gamma )}{(d_S+\lambda )(d_V+\gamma )(\sigma +d_E)(d_I+\theta )}>0\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial \beta }\hfill & =\hfill & \frac{\sigma \mathsf{\Lambda }\lambda }{(d_S+\lambda )(d_V+\gamma )(\sigma +d_E)(d_I+\theta )}>0\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial d_S}\hfill & =\hfill & -\frac{\mathsf{\Lambda }\left({\beta }_1(d_V+\gamma )+\beta \lambda \right)}{(d_V+\gamma )(d_I+\theta )}\frac{1}{{(d_S+\lambda )}^2}<0\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial \lambda }\hfill & =\hfill & \frac{\sigma \mathsf{\Lambda }\left(\beta d_S-{\beta }_1(d_V+\gamma )\right)}{(d_V+\gamma )(\sigma +d_E)(d_I+\theta )}\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial d_V}\hfill & =\hfill & -\frac{\sigma \mathsf{\Lambda }}{(d_S+\lambda )(\sigma +d_E)(d_I+\theta )}\frac{\beta \lambda }{{(d_V+\gamma )}^2}<0\hfill \end{array}$$

$$\begin{array}{ccc}\frac{\partial {\mathcal{R}}_0}{\partial \gamma }\hfill & =\hfill & -\frac{\sigma \mathsf{\Lambda }\left({\beta }_1(d_V+\gamma )+\beta \lambda \right)}{(d_S+\lambda )(d_I+\theta )(\sigma +d_E)}\frac{1}{{(d_V+\gamma )}^2}<0\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial d_I}\hfill & =\hfill & -\frac{\sigma \mathsf{\Lambda }\left({\beta }_1(d_V+\gamma )+\beta \lambda \right)}{(d_S+\lambda )(d_V+\gamma )(\sigma +d_E)}\frac{1}{{(d_I+\theta )}^2}<0\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial \theta }\hfill & =\hfill & -\frac{\sigma \mathsf{\Lambda }\left({\beta }_1(d_V+\gamma )+\beta \lambda \right)}{(d_S+\lambda )(d_V+\gamma )(\sigma +d_E)}\frac{1}{{(d_I+\theta )}^2}<0\hfill \\ \frac{\partial {\mathcal{R}}_0}{\partial d_E}\hfill & =\hfill & -\frac{\sigma \mathsf{\Lambda }\left({\beta }_1(d_V+\gamma )+\beta \lambda \right)}{(d_S+\lambda )(d_V+\gamma )(d_I+\theta )}\frac{1}{{(d_E+\sigma )}^2}<0.\hfill \end{array}$$

We have the plot of ${\mathcal{R}}_0$ in terms of their different parameters in Figure 6. The values of other parameters used for plotting while fixing the parameter that is plotted against ${\mathcal{R}}_0$ for each plot are taken from Table 2.

The sensitivity index technique will help to measure the most sensitive parameters for the fundamental reproductive number ${\mathcal{R}}_0$. The fundamental reproduction number’s normalised sensitivity index is provided by $S_{\mathcal{Z}}^{{\mathcal{R}}_0}=\frac{\partial {\mathcal{R}}_0}{\partial \mathcal{Z}}.\frac{\mathcal{Z}}{{\mathcal{R}}_0}$, where $\mathcal{Z}$ is a parameter, as defined earlier. We obtain

$$\begin{array}{cc}\hfill S_{\sigma }^{{\mathcal{R}}_0}& =d_E/(d_E+\sigma ),\hfill \\ \hfill S_{\mathsf{\Lambda }}^{{\mathcal{R}}_0}& =1,\hfill \\ \hfill S_{{\beta }_1}^{{\mathcal{R}}_0}& =\frac{{\beta }_1(d_V+\gamma )}{{\beta }_1(d_V+\gamma )+\beta \lambda },\hfill \\ \hfill S_{\beta }^{{\mathcal{R}}_0}& =\frac{\beta \lambda }{{\beta }_1(d_V+\gamma )+\beta \lambda },\hfill \\ \hfill S_{d_S}^{{\mathcal{R}}_0}& =-d_S/(d_S+\lambda ),\hfill \\ \hfill S_{\lambda }^{{\mathcal{R}}_0}& =\frac{\lambda (-{\beta }_1(d_V+\gamma )-\beta \lambda +\beta (d_S+\lambda ))}{(d_S+\lambda )({\beta }_1(d_V+\gamma )+\beta \lambda )},\hfill \\ \hfill S_{d_V}^{{\mathcal{R}}_0}& =-\frac{d_V\beta \lambda }{(d_V+\gamma )({\beta }_1(d_V+\gamma )+\beta \lambda )},\hfill \\ \hfill S_{\gamma }^{{\mathcal{R}}_0}& =-\frac{\gamma \beta \lambda }{(d_V+\gamma )({\beta }_1(d_V+\gamma )+\beta \lambda )},\hfill \\ \hfill S_{d_I}^{{\mathcal{R}}_0}& =-d_I/(d_I+\theta ),\hfill \\ \hfill S_{\theta }^{{\mathcal{R}}_0}& =-\theta /(d_I+\theta ),\hfill \\ \hfill S_{d_E}^{{\mathcal{R}}_0}& =-d_E/(d_E+\sigma ).\hfill \end{array}$$

The sensitivity indices obtained by using the parameters values in Table 2 are presented in

Table 4. The sensitivity index of ${\mathcal{R}}_0$ with respect to the parameters $\mathcal{Z}$ of the system (1).

Parameter	Sensitivity Index
$\sigma$	0.019
$\mathsf{\Lambda }$	1
${\beta }_1$	0.99
$\beta$	0.00187597
$d_S$	−0.99727
$\lambda$	−0.00085543
$d_V$	−0.00004068
$\gamma$	−0.00183557
$d_I$	−0.3579566
$\theta$	−0.642
$d_E$	−0.019

. Four of the sensitivity indices are positive, while the others are negative, as can be seen in Table 4. We notice a difference in sensitivity between the parameters’ contact between V and I ($\beta$) and contact between S and I (${\beta }_1$). This difference is minimal for $\beta$ because the vaccinated population has developed some kind of immune system when in contact with the infected population, which is due to the fact that significant proportion of the population have been vaccinated, and maximal for ${\beta }_1$ because the parameter determines the spread of the disease during contact between the susceptible population and the infected population, hence the high sensitivity of this parameter. We conclude that increasing the recovery rate and vaccination rate will aid in decreasing the ${\mathcal{R}}_0$, which affirms the effect of vaccination, and by extension, it is important to encourage a significant chunk of the population to get vaccinated, which will help to combat the spread of the virus.

5. Discussion

To demonstrate the influence of the spatial diffusion in the model, we used the method presented in [28] by using implicit–explicit finite differences, where the Laplacian term is discretised implicitly (for numerical stability, allowing for a relatively large time stepping), whereas nonlinear terms are handled explicitly. We assume $\mathsf{\Lambda }=0$ and the diffusion term $d_1\Delta =d_2\Delta =0.00005$. Furthermore, we assume $S_0\left(x\right)=1.3+cos\left(3\pi x\right)$ and $I_0=0.01exp(-1000x)$ with $E_0=0,V_0=1,I_0=0$ and $I_0\left(1\right)=0.001$ in order to have a wave-like or exponential propagation of the disease. When varying $\beta$, we chose ${\beta }_1=1$ and $\lambda =0.7$ and when varying $\lambda$, we chose $\beta ={\beta }_1=1$. Other parameters were chosen from Table 2.

For the sensitivity analysis, we discovered that some of the parameters, i.e., $\mathsf{\Lambda },{\beta }_1,\beta ,$ and $\sigma$, cause ${\mathcal{R}}_0$ to increase, while other parameters, corresponding to $d_S,d_E,\theta ,d_I,d_V$, $\lambda$ and $\gamma$, cause it to decrease.

When ${\mathcal{R}}_0<1$, all compartments tend to zero except for the susceptible population, which increases with rate depending on $\mathsf{\Lambda }$. Moreover, the extinction of other compartments is very fast (less than 50 days). For ${\mathcal{R}}_0>1$, the susceptible and vaccinated populations tend towards zero, and the majority of the population becomes infected. Whenever ${\mathcal{R}}_0<1$, the behaviour of the disease in the population does not change much with $\lambda$. However, when ${\mathcal{R}}_0>1$, if we increase $\lambda$, the infected population will disappear, and the majority of the population becomes vaccinated. Hence, we cannot take $\lambda <0$ since the parameters are positive.

Figure 3 shows both (i) the influence of a nonzero $\mathsf{\Lambda }$ flow of entrants into the susceptible population: the asymptotic value of $S(x,\infty )$ is not zero, contrary to what we observe in Figure 7, where $\mathsf{\Lambda }=0$ and $S(x,\infty )=0$, and (ii) the influence of the vaccination parameter $\lambda$, which, when it increases, increases the ratio $V(x,\infty )/S(x,\infty )$ (it goes from the value 4.4 to the value 25 between Figure 2, where $\lambda =9.1\times {10}^{-5}$ and $\mathsf{\Lambda }=350$, and Figure 3, where $\lambda =0.17$ and $\mathsf{\Lambda }=\mathrm{15,000}$). When $\lambda =0$ (Figure 4 and Figure 5), we can see the influence of ${\mathcal{R}}_0:S(x,\infty )=0$, if ${\mathcal{R}}_0>1$, even if $\mathsf{\Lambda }$ is high, equal to 15,000 (Figure 5), and $S(x,t)$ remains increasing with t if ${\mathcal{R}}_0<1$, even if $\mathsf{\Lambda }$ is small, equal to 350 (Figure 4).

Figure 7a,b shows the influence of the parameter $\beta$ and that of the diffusion: when $\beta$ goes from 0.1 to 0.2, the asymptotic values $S(x,\infty )=0$ and $V(x,\infty )=0$ are reached later, and the initial shape of $S(x,0)$ is kept longer. Figure 7c,d shows the influence of the parameter $\lambda$ and that of the diffusion: when $\lambda$ goes from 0.1 to 0.2, the asymptotic values $S(x,\infty )=0$ and $V(x,\infty )=0$ are also reached later, and the initial shape of $S(x,0)$ is kept longer. If $\lambda$ is small, that is, less than 0.5, the trend of $S(x,300)$ will be close to linear. This shows, as expected, that we must wait a long time to obtain the asymptote for high values of $\beta$ and $\lambda$, probably 500 or more time steps. Furthermore, for increasing $\lambda$, the disappearance of S would correspond to the appearance of a high asymptotic value of the evolution of V.

If the epidemic process starts with a non-homogeneous spatial initial condition for ${\beta }_1{S}_0(x,t)$, i.e., for the initial flow of exposed individuals resulting from a not spatially uniform contact rate per day between the susceptible and infected population, then Figure 7 shows the evolution of the susceptible and exposed sub-populations until their asymptotic state is reached. The susceptible individuals disappear progressively depending on both the exposition rate ${\beta }_1$ and the vaccination rate $\lambda$, while their number seems to remain constant for high values of $\lambda$.

6. Conclusions

In this paper, we have been able to develop a mathematical modelling of COVID-19 using diffusion equations built from some previous works to construct a more realistic model of COVID-19. Using this model, we were able to prove the existence, uniqueness, and local stability of constant stationary solutions. Another mathematical perspective we contributed was to build an appropriate Lyapunov function for these constant solutions in order to prove the global stability of the model. Then, we performed sensitivity analysis of ${\mathcal{R}}_0$ in order to understand the dynamics of each parameter, which gives a better insight into how to control the evolution of the disease in the population, which is of interest to policy makers and public health experts. The numerical analysis showed agreement with the results of the qualitative analysis. Moreover, the proof of existence and stability of non-spatially homogeneous solutions and application of the model to time series data from different countries will be the subject of a future investigation.