A numerical study of COVID-19 epidemic model with vaccination and diffusion

1.   Introduction

The COVID-19 pandemic is a global outbreak, caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) ^[1]. This infection is highly contagious and is transmittable from person to person via respiratory droplets, produce as a consequence of coughing and sneezing of an infectious person. Although the severity and mortality rate of this novel infection is notably reduced, still it is a huge panic for human health as well as economy around the globe. This infection is still the main cause of restricting lives of humans in many regions. Despite huge research on this novel disease, scientists are focusing to explore its different aspects. According to recent statistics, more than 621.8 million infected cases are reported globally. Out of the total reported infected cases about 6.5 million death cases and about 601.9 million recovered cases are recorded ^[1]. Many pharmaceutical and non-pharmaceutical preventing strategies were suggested and implemented to overcome the infection. For instance, taking vaccines, staying at home, usage of facemasks specifically in public areas, avoiding crowded places, etc ^[1,2].

One of the main causes of the infection spread is the spatial movement of population from one region to another. The assessment of geographical spreading patterns as well as the impact of those factors causing the transmission and eradication of infection can be carried out using mathematical models ^[3,4,5,6,7]. In existing literature, most of COVID-19 epidemic models have been constructed via ordinary differential equations (ODEs). In contrast to ODEs models, the mathematical models based on PDEs, particularly the reaction-diffusion is more appropriate and have promising results describing the dynamics of an infectious disease. Numerous mathematical spatio-temporal epidemic models have been developed to explore the transmission patrons of various infectious diseases including COVID-19. These problems were solved numerically using different numerical techniques to analyze the dynamics of different infectious and viral diseases. Such as, in ^[8] the authors formulated a reaction-diffusion deterministic model to analyze the spatio-temporal dynamics of influenza. They further utilized the finite-difference scheme in order to visualize the disease dynamics. In ^[9], the authors developed a diffusive vaccine epidemic model with variable transmission coefficient addressing the spatio-temporal dynamical aspects of influenza disease. Moreover, in ^[9] the authors presented a detailed analysis of the model. A compartmental spatio-temporal epidemic model is considered by Jawaz et al. ^[10], for investigating the complex dynamics of HIV/AIDS and proposed the non-standard finite difference numerical technique for the solution due to its positivity preserving property. Ahmed et al. ^[11], discussed the numerical solution of the whooping cough epidemic model by using an operator-splitting finite difference scheme. The numerical techniques, operator-splitting based finite difference and meshless prodedures have been used for solving the spatio-temporal epidemic models in Refs. ^[12,13]. Asif et al. ^[14] discussed the approximate solution of two spatio-temporal biological models by using meshless and finite difference procedures based on operator splitting techniques. A structure-preserving non-standard finite difference operator-splitting procedure has been implemented by Ahmed et al. ^[15], to obtain the numerical solution of the reaction-diffusion epidemic model. Sokolovsky et al. ^[16] examined analytic and numerical solution of spatio-temporal epidemiological models and analyze different stages of the epidemic. A novel reaction-diffusion model describing the dynamics of COVID-19 is developed in ^[17]. A reaction-diffusive spatial model addressing the dynamics of pandemic in Greece and Andalusia as a case study is presented in ^[18].

In this paper, we present a spatio-temporal vaccine mathematical model to address the dynamics of COVID-19 in a spatially heterogeneous environment. The model is actually the spatial extension of ^[19]. Initially, we present some basic and necessary mathematical properties of the proposed spatio-temporal model. The model is then solved numerically using an efficient numerical scheme to present the simulation results. The paper outlines are categorized as: Section 2 describes the formulation of spatio-temporal model. A basic mathematical analysis of the problem is presented in Section 3. Numerical solution of the model and detailed simulation with and with diffusion is investigated in Section 4. Finally, the concluding remarks are summarized in Section 5.

2.   Spatio-temporal COVID-19 model description

In this section, a compartmental epidemic reaction-diffusion model is formulated in order to analyze temporal and spatial dynamics of COVID-19 outbreak. The total population at time $t\geq0$ and spatial point $x\in[-2, 2]$ is denoted by $\mathcal{N}(x, t)$ and assumed to be constant i.e., in the population birth and deaths rates are equal. The population $N(x, t)$ is further categorized into sub-classes denoted by $S(x, t), E(x, t), I_S(x, t), I_A(x, t), V(x, t)$ and $R(x, t)$ describing the susceptible, latent, symptomatically, asymptomatically, vaccinated and recovered population groups respectively. The spatial distribution of total available population at time $t\geq0$ and $x\in \Xi = [-2, 2]$ is

The following assumptions are taken into the account to construct the desired mathematical model:

● Each newborn is susceptible i.e., having the capability of getting infection.

● The population in $I_S$ and $I_A$ are capable to transmit the infection.

● Vaccination is not permanent. The vaccinated individuals can wane vaccine immunity and become susceptible again.

The desired nonlinear spatio-temporal COVID-19 epidemic model is presented as follows:

subjected to non-negative initial conditions (ICs):

where $S(0, x) = S_0 \geq 0, E(0, x) = E_0\geq 0, I_S(0, x) = {I_S}_0 \geq 0, I_A(0, x) = {I_A}_0 \geq 0, V(0, x) = R(0, x) = 0$ and $x\in [-2, 2]$ is the computational domain for the proposed model (2.1). The biological description of embedded parameters is shown in the Table 2. Further, $D_S, D_E, D_{I_S}, D_{I_A}, D_{V}$ and $D_R$ are the coefficients of diffusivity for the respective population groups. Furthermore, no flux boundary conditions are being imposed for the proposed model given by:

3.   Theoretical analysis of model

In the current section, we establish some of the important theoretical analyses of the proposed COVID-19 reaction-diffusion model (2.1). The well-known threshold parameter of the epidemiological model is derived. Moreover, the disease free and endemic equilibria of model (2.1) are computed. In addition, the stability results of the model are accomplished.

3.1.   Well-posedness of model (2.1)

The well-posedness of partial differential equations (PDEs) with some initial and boundary conditions is important to verify in order to claim the existence and uniqueness of its solution. For this purpose, the semi-group theory approach can be utilized to check when a problem is well-posed. To proceeds with this analysis, a Banach space $X = C(\Xi; \mathbb{R})$ is considered, which is a space of real-valued continuous function $\Upsilon$ define on $\Xi$ and $X_{+}$ is a positive cone. The space is combined with norm $\|\Upsilon\|_{X} = \sup_{x\in\Xi}|\Upsilon(x)|$. Assume that $\mathcal{A}_{\xi}$ is a linear operator defined by

where $\xi$ represent the state variables $S, E, I_S, I_V, V, R$ and $\Delta$ is the second order differential operator $\frac{\partial^2}{\partial x^2}$. Further,

After utilizing the well-known facts regarding the operator $\mathcal{A}_{\xi}$, i.e., $\mathcal{A}_{\xi}$ is an infinitesimal generator of a strongly continuous semi-group $\{e^{t\mathcal{A}_{\xi}}: t\geq0 \}$ of linear operators in $X$ ^[20]. On other hand the operator,

is an infinitesimal generator of a strongly continuous semi-group $\{e^{t\mathcal{A}}: t\geq0 \}$ of linear operators in $Y = X^{6}$ ^[20], where $Y$ is a Banach space with norm define by,

and $Y_{+} = X_{+}^{6}\subset Y$ be a positive cone in $Y$. Let $F$ be a nonlinear operator define on $Y$ as:

for $(\phi_{1}, \phi_{2}, \phi_{3}, \phi_{4}, \phi_{5}, \phi_{6}) \in Y_{+}$. Thus model (2.1) can be expressed in compact form as:

where

Clearly, $F$ is Lipschitz continuous on $Y_{+}$. Since the following inequality holds,

where $L = \mu+\mu_1+2(\psi_{V}+\kappa+\gamma_1+\gamma_2+\omega_V)$. Hence, by Theorem $3.3.3$ in ^[21], the following result is concluded.

Proposition 3.1. Assume that $\mathcal{A}$ be an operator define by (3.1) and for each $u_{0} = (S_0, E_0, {I_S}_0, {I_S}_0, V_0, R_0)^T \in Y_{+}$ there exist a unique continuously differentiable solution $u(t)$ of (3.2) define on some maximum interval of existence $[0, \tau]$ such that

3.2.   Positivity and boundedness

To show that the solution of model (2.1) is non-negative, we follow the lemma provided in ^[22,23,24].

Lemma 3.2. Assume that $u \in C(\bar{\Xi}\times [0, \tau))\bigcap C^{2, 1}(\Xi\times(0, \tau))$ satisfy

where $k(x, t) \in C(\bar{\Xi}\times [0, \tau))$. Then $u(x, t) \geq 0$ 0n $\bar{\Xi}\times [0, \tau)$ and $u(x, t) > 0 \ \hbox{or} \ u(x, t)\equiv 0$ in $\Xi\times [0, \tau)$

In view of the above Lemma 3.2, the following result is obtained.

Proposition 3.3. The solution $(S(., t), E(., t), I_S(., t), I_A(., t), V(., t), R(., t))$ of model (2.1) is non-negative on $\Xi\times [0, \tau)$ provided that the initial condition $(S(.), E(.), I_S(.), I_A(.), V(.), R(.))$ is non-negative.

Proof. Suppose $\mathcal{K}(x, t) = \beta(I_S(x, t)+\beta_1 I_A(x, t))\frac{1}{N}+(\mu+\omega_{V})$. Then from first equation of model (2.1)

Clearly, (3.3) fulfill the conditions of Lemma 3.2. So, by direct application of Lemma 3.2 $S(x, t)\geq 0$ on $\Xi\times [0, \tau)$. The non-negativity of the rest of state variables $E(x, t), I_S(x, t)$, $I_A(x, t)$, $V(x, t)$ and $R(x, t)$ can be prove in similar manner.

Next, the boundedness of solution of model (2.1) is demonstrated in order to prove the time existence of global solution. For this the following result is derived.

Theorem 3.4. The solution $(S(., t), E(., t), I_S(., t), I_A(., t), V(., t), R(., t))$ of model (2.1) is bounded for all $t\geq 0.$

Proof. By sum up all equations of model (2.1), we obtained

Integrating over $\Xi = [-2, 2]$ yields to,

By utilizing the no flux boundary conditions imposed to the problem we have

Therefore,

It follows that

Thus,

Invariant region

The transmission dynamics of the problem is studied in the biological feasible region given by,

where,

3.3.   Steady-state of the model

The proposed spatio-temporal vaccine model (2.1) has two steady-states. The disease free steady-state is

The endemic equilibrium is stated as follows: ${\mathcal{E}}_1(S^{***}, E^{***}, I_S^{***}, I^{***}_A, V^{***}, R^{***})$ with

where,

Substituting values from (3.4) into the following force of infection:

After some manipulations, the non-zero equilibria of the model satisfies the following equation

with the coefficients

The basic reproduction number $\mathcal{R}_{0}$ is the epidemiological threshold parameter called the basic reproductive number. By utilizing the next generation technique $\mathcal{R}_0$ is computed as fallow:

3.4.   Stability of the model equilibria

In this subsection, we provide the stability analysis of the spatio-temporal COVID-19 model (2.1). In order to investigate the local stability of the model at $\mathcal{E}_0$, we follow the Fourier series expansion technique discussed in ^[15]. The variational matrix obtained after linearizing the spatio-temporal COVID-19 model (2.1) is given as follows:

where,

The corresponding characteristic equation of the variational matrix ${\bf{V}}$ is obtained as:

The eigenvalue $-(\mu+D_{S}k^2)$ has negative real part while the following coefficients of (3.10) are obtained

After substituting values from (3.9), it can be shown that $a_i > 0$ for $i = 1, 2, ..., 5$ under the condition $\mathcal{R}_{0} < 1$. Furthermore, after some manipulation, the Routh-Hurwitz criteria for characteristic polynomial (3.10) arrived when $\mathcal{R}_{0} < 1$. Thus, the spatio-temporal COVID-19 (2.1) is locally asymptotically at $\mathcal{E}_0$ whenever, $\mathcal{R}_{0} < 1$.

Moreover, to inveterate the global asymptotical stability of the spatio-temporal model (2.1) at $\mathcal{E}_0$, the well-known Lyapunov function approach is implemented. The subsequent theorem yield the desired result.

Theorem 3.5. If $\mathcal{R}_0\leq1$, the equilibrium point $\mathcal{E}_{0}$ of model (2.1) is globally asymptotically stable.

Proof. To prove the result, consider Lyapunov functional as follows:

where $g_1 = \dfrac{\mu+\psi_{V}}{\mu+\psi_{V}+\omega_{V}}$, $g_2 = \dfrac{\beta g_1}{(\mu+\mu_1+\gamma_1)(\mu+\gamma_2)}$, $g_2 = \dfrac{\beta \beta_1g_1}{\mu+\gamma_2}$ and $\tau$ is the maximum time.

Then, the time derivative $\mathcal{F}(t)$ is given by

By utilizing the no flux boundary conditions implies,

Since,

Therefore,

It is clear that $\dfrac{d}{dt}\mathcal{F}(t) \leq 0$ for all $t\geq 0$ and $x\in\Xi$ if and only if $\mathcal{R}_0\leq1$. Moreover $\dfrac{d}{dt}\mathcal{F}(t) = 0$ whenever, $E(x, t) \rightarrow 0$, for all $t\geq 0$ and $x\in \Xi$. It follows from the well-known LaSalle's invariance principle presented in ^[25], that the equilibrium $\mathcal{E}_0$ of spatio-temporal model (2.1) is globally asymptotically stable if $\mathcal{R}_0\leq 1$.

4.   Numerical method and simulation

This section presents the numerical scheme of the proposed spatio-temporal COVID-19 model (2.1). The operator-splitting finite difference approximation technique is utilized for the numerical solution and is described in detail as follows:

4.1.   Numerical scheme

The operator splitting approach is a well-known and efficient numerical technique. It has been applied to solve nonlinear partial differential equations arising from nonlinear and complex problems. In this scheme, we split the differential operator into sub-operators and split the problems under consideration for solution into sub-problems, where each of these problems correspond to a particular physical phenomenon. Further, the operator-splitting technique is implemented to solve the proposed reaction-diffusion epidemic model (2.1). Each differential equation in the model (2.1) denotes the physical process consisting of two individual processes: the population with different disease status interacting with each other and diffusing spatially in one direction, required different time steps. It is therefore important to split the time operator and shift model (2.1) into two sub-systems: The resulting non-linear reaction system used for time step $t_0$ to $\frac{1}{2}dt$ is,

and the linear diffusion system used for time step $\frac{1}{2}dt$ to $t^{n}$ is,

using conventional finite difference approximations, the first order time derivative in (4.1) and (4.2) is approximated by first order-forward difference,

while the second order spatial derivative in (4.2) is approximated by second order central finite difference given by,

where $\xi$ represent each of the variables $S, E, I_S, I_A, V$ and $R$. The sub-systems (4.1) and (4.2) in iterative form can be written as fallow:

and

4.2.   Simulation and discussions

In this section, numerical simulations of the proposed reaction-diffusion epidemic model (2.1) is performed using the iterative scheme derived in (4.5) and (4.6). The purpose of the model simulation is to analyze the impact of diffusion phenomena under several control strategies. Model (2.1) is simulated in two ways: with and without diffusion in parallel with other control measures such as, the effective contact rate $\beta$, the transmissibility of infection relative to asymptomatically infected individuals $\beta_{1}$, the vaccine waning rate $\psi_{V}$ and the effective vaccination rate $\omega_{V}$. The simulation is performed based on initial conditions (2.2) and estimated parameters given in Tabe 2, while the diffusion coefficients are assumed as $D_S = 0.00005, D_E = 0.0005, D_{I_S} = 0.0001, D_{I_A} = 0.001, D_V = 0.0001, D_R = 0$ so that the quantity

where $\xi$ represents each of the state variables. In order to perform simulation, the spatial step is taken as $dx = 0.06$ and the time step $dt = 0.02$ days based on Von Neumann stability criteria ^[26]. Figure 1 demonstrate in detail the dynamics of exposed, symptomatically and asymptomatically infected individuals, with and without diffusion using the baseline values of the model (2.1) embedded parameters. From the corresponding plots, it is observed that the number of infected individuals in all categories reduces significantly, in the case of diffusion as compared to without diffusion. Physically, it means that restricting public gathering plays an important role in reducing the transmission of infection. In the subsequent sections, we present the impact of different control interventions as mentioned before on the dynamics of exposed, symptomatically and asymptomatically infected individuals in case of diffusion and without diffusion. In Table 3, we compare the results of the present technique with fully explicit Euler scheme in terms of CPU time.

4.2.1.   Dynamics without diffusion at $x=0$

In this case, the proposed model (2.1) is simulated without diffusion by assuming $D_S = D_E = D_{I_S} = D_{I_A} = D_V = D_R = 0$. In addition, the effects of the above control interventions are analyzed considering different situations. The dynamics of exposed, symptomatically and asymptomatically infected individuals, in this case are graphically illustrated for the time period $0$ to $700$ days. A detailed interpretation is presented in Figures 2–4. Figure 2 describes the influence of control measure $\beta$ on the dynamics of COVID-19. The impact is observed for the estimated baseline value of the effective-contact rate $\beta$, without diffusion. Simulations are performed by reducing the baseline value of $\beta$ to 10, 20 and $30\%$. It is found that with a $20\%$ decrease in social contacts ($\beta$), $53\%$ reduction is observed in infected individuals, while $30\%$ reduction in social contacts, i.e., strengthening lock-down yield to $75\%$ decrease in the number of exposed, symptomatically and asymptomatically infected individuals. Therefore, this analysis shows that enforcing a quarantine policy is beneficial and helps control the spread of infection. The effect of the parameter $\beta_{1}$ on the dynamics of exposed asymptomatic and symptomatic infected individuals is shown in Figure 3. The corresponding dynamics are observed without diffusion for 10, 20 and $30\%$ reduction in $\beta_{1}$ to its estimated baseline value given in Table 2. Simulations show that the number of cases in each category decreased by 23, 45 and $64\%$ respectively. The model (2.1) is simulated without diffusion to investigate the impact of various vaccination rates on the dynamics of exposed, symptomatically and asymptomatically infected individuals. Parameter $\psi_{\nu}$ is reduced by 10, 20 and $30\%$ to its estimated baseline value. The resulting plots are shown in Figure 4.

4.2.2.   Dynamics with diffusion at $x=0$

This section presents a simulation of the model (2.1) with diffusion. The values of corresponding diffusion coefficients are assumed as $D_S = 0.00005, D_E = 0.0005, D_{I_S} = 0.0001, D_{I_A} = 0.001, D_V = 0.0001, D_R = 0$. Figure 5 demonstrates the dynamics of the infected population with variation in $\beta$ in the presence of diffusion. The parameter $\beta$ is reduced with the same rate discussed in the previous section. By implementing an isolation policy with diffusion, one can analyze that $20\%$ reduction in social contacts yields to $75\%$ decrease in infected individuals. Moreover, with $30\%$ reduction in parameter $\beta$ the infected individuals are decreased to $93\%$. Thus, the isolation policy with diffusion is more effective and plays a significant role in controlling the prevalence of infection.

Figure 6, describes the effect of control intervention $\beta_{1}$ on the disease incidence with diffusion. The simulation results are performed by reducing $\beta_{1}$ to 10, 20 and $30\%$. One can observe that in this case with a $30\%$ reduction in $\beta_{1}$ the peaks of infected population curvesis reduced to $95.1\%$. It shows that using the suggested control measure $\beta_{1}$ with diffusion is more beneficial in controlling infection prevalence as compared to without diffusion.

Figure 7 describes the impact of variation in parameter $\psi_{\nu}$ on the dynamics of infected individuals. The dynamics are observed in the presence of diffusion. We reduce the estimated value of $\psi_{V}$ to 10, 20, and $30\%$ and analyze the behavior of exposed, symptomatically, and asymptomatically infected individuals. Simulation results show that reduction in $\psi_{\nu}$ to $30\%$ decreased the population in the respective infected classes to $72.1\%$. These details are provided in Table 4 and 6.

4.2.3.   Dynamics with diffusion at $x=1$

This section presents a simulation of the model (2.1) with diffusion at $x = 1$. The corresponding plots are shown in Figures 8–10. The time evolutionary trajectories of exposed, symptomatically and asymptomatically infected individuals at $x = 1$ with a reduction in $\beta$ at different rates are shown in Figure 8. Initially, an increase is observed in the respective curves for the baseline value of $\beta$. It is due to the population is diffuses from a higher concentration place to a low concentration region. It has been observed that with diffusion in case of strict lock-down (with a $30\%$ reduction in $\beta$) lowest projected peaks are noticed. Figure 9 depicts the influence of $\beta_{1}$ on the dynamics of exposed, symptomatically and asymptomatically infected individuals with diffusion at $x = 1$. The simulation is obtained for a reduction in $\beta_1$ with various rates. According to initial condition (2.2) the population is concentrated at the origin and decreases exponentially at $x = 1$. As the population diffuses the concentration increase for $\beta_1 = 0.7459$ at $x = 1$ and significantly reduces with a $30\%$ reduction in $\beta_1$. Finally, the impact of $\psi_{\nu}$ (vaccine waning rate) on the dynamics of exposed, symptomatically and asymptomatically infected individuals with diffusion at $x = 1$ is described in Figure 10. By reducing the vaccine waning rate to 10, 20 and $30\%$, a reasonable decrease is analyzed in the respective compartments of infected individuals.

Figure 11 describes mesh plots of model (2.1), which represent the Spatio-temporal evolution of exposed, symptomatically and asymptomatically infected individuals for all time and spatial points in $[5, −5]$. As long as the stability condition given by (4.7) fulfills the proposed scheme produces a stable solution. Moreover the coefficients of diffusivity $D_{\xi}$, where $\xi$ represents each of the state variables $S, E, I_S, I_A, V$ and $R$ are chosen so that (4.7) is satisfied. The corresponding mesh plots show the consistent behavior of the proposed operator-splitting finite-difference numerical scheme. Clearly, it preserves the positivity property of the solution, i.e., the solution stays positive for all $t>0$ and spatial points within the domain of definition $[5, −5]$.

5.   Conclusions

In this study, a new mathematical model based on PDEs is constructed to explore the dynamics of COVID-19 outbreak in a spatially heterogeneous environment. The primary focus is to analyze the impact of various significant measures with and without diffusion on the dynamics and control of the pandemic. The basic mathematical properties including existence, uniqueness, positivity, and boundedness of the diffusive vaccinated model are presented. Moreover, the spatio-temporal COVID-19 mathematical model is numerically solved using the finite difference operator-splitting scheme with uniform and non-uniform initial conditions. Based on the derived iterative scheme the model is simulated for various key parameters such as $\beta$, $\beta_1$, $\psi_V$ with and without diffusion. Further, the dynamics are observed spatially at $x = 0.0$ and $x = 1.0$. The respective results are also shown in tabular form. These results revealed that with diffusion (and at $x = 1$) the reduction in disease transmission coefficients $\beta$ and $\beta_1$ with a 30% rate, the projected peaks of infected individuals reduce upto 98 and 78% respectively as shown in Table 5. The corresponding evolutionary trajectories showed a clear picture of diffusion phenomena and the number of infective individuals significantly reduced with diffusion as compared to without diffusion.

Acknowledgments

This Project was funded by the Denship of Scientific Research (DSR)at King Abdulaziz University, Jeddah, Saudi Arabia under grant no. (G: 394-130-1443). The authors, therefore, acknowledge with thanks DSR for technical and financial support.

Conflict of interest

The authors declare there is no conflict of interest.

Data availability statement

The authors declare that the all data supporting the findings of this study are available within the article.