Computational Modeling of Reaction-Diffusion COVID-19 Model Having Isolated Compartment

1  Introduction

In December 2019, the world encountered the most destructive disease named COVID-19. Humankind did not face such crises after World War II, which wreaked great havoc on the nations with the dwindling economies and inadequate health care facilities. It raised the causality rate by 5,534,735 and the fall in trade, movements, and traveling. According to a survey done by international traveling cooperation, about 4 billion people cannot travel due to traveling curbs. This pandemic ceased the government and private medical units as most countries were unprepared and were unaware of the multiplication of causative agents. Today, the coronavirus has produced a new chapter for mathematical researchers who found modeling the most appropriate tool for investigating the spread of a disease in a community.

According to the union health ministry, 653 omicron cases were detected across India in December 2021. However, a dramatic rise has been seen from 29 to 30 December 2021. On 30 December 2021, 13,187 new cases of omicron variant were reported, which were 77% of the total population of India, and the states like Mumbai, Delhi, Pune, Bengaluru, Chennai, Thane, Kolkata, and Ahmedabad were considered as the most susceptible ones. Mumbai showed an 80% rise in the cases daily with an approximation of 2510 cases, leading to a 400% rise weekly. Additionally, 923 cases were detected in New Delhi, which revealed a 600% rise weekly, along with Bengaluru 400 cases (90%), Chennai 294 (100%), and Mumbai (15%). These cases are of the omicron variant, the type of SARS-CoV-2 virus, and have been in destructive action in India since December 2021. The target of vaccination endorse 63% of adults.

The development of vaccines and medicines rectified the pandemic but did not reduce the chance of disease spreading, which has become a great challenge for scientists. This problem can be solved using equilibrium points of COVID-19. People have followed intense SOPs (standard operational procedures) to combat this variant, including quarantine, hand wash, suitable distance, and isolation.

NIH (National Institute of Health) reported that Pakistan had 75 cases since 27 December. Thirty-three were from Karachi, 17 from Islamabad, and 13 from Lahore. Individuals in these cases were found to be international travel. NIH revealed that these patients were immediately isolated and contacted their relatives to control the spread in a statement. On 13 December, The News, the largest-selling paper in Pakistan, coated the figures of NIH and revealed that the first case of omicron was from Karachi.

NIH further stressed the need to follow SOPs and advise people to get vaccinated. The government of Pakistan allowed all the vaccines to be administered as early as possible to get rid of the dire consequences of the omicron variant. NCOC (National Command and Operation Centre) urged people to administer vaccines and booster doses with criteria and precautions. NIH revealed that omicron is a lethal variant whose aftermaths have been multiplying in several countries.

Today, the NCOC issued the latest COVID-19 statistics for the last 24 h in Pakistan, which was 0.69%. Yesterday, 291 new COVID-19 cases were reported while 41,869 diagnostic tests were conducted, whereas three died due to the omicron variant.

Coronavirus was first reported in China city of Wuhan and named COVID-19 by WHO as it shares a common phylogenetic lineage with SARS COVID-II. It is mainly transferred from one person to another due to coughing, sneezing, or contact with infected individuals. Salivary droplets released by infected individuals are denser than air and immediately fall on the ground or nearby. WHO allocated some safety measures, including the distance of six feet, hand washing, wearing masks and gloves, and the quarantining [1,2]. Reportedly, 1,458,000 people were getting infected in more than 180 countries, raising the number of infected individuals by 4 million [3].

A model named as retrofitted state SIR model was proposed by [4] to predict and reckon the numeric of infected and susceptible individuals. Nesteruk focuses on the epidemic and calculates the number of infected individuals. Based on his assumptions, the mortality rate was higher than estimated. Nesteruk et al. [5–8] proved that infection multiplies quickly in densely packed areas, which moves the hypothesis towards applying social distancing and quarantining measures. In [9], authors have investigated the SIR model to estimate the primary cause of coronavirus spread.

Okhuese [10] demonstrated a framework of a stochastic model. The condition dictates extinction and persistence. Also, they debated the threshold of the stochastic model proposed when small or large noise occurs. A method of the potential besides misinformation transmission within the population is the basic reproduction number, mathematically disease-free threshold and stability are related to an epidemic peak and final size [11–15]. The purpose of the present article is to investigate the effects of diffusion on the spread of disease; the reproductive number of the model is given, which is the key element of such models, and the numerical results are obtained through some schemes such as NSFD, FTCS, and Crank Nicolson scheme.

A disease-free equilibrium’s local and global stability is linked to the calculation and epidemiological interpretation of this threshold parameter [16–18]. Following the spread of the disease, researchers have activated to speed innovative diagnostics and are in the throes of several vaccines to guard against COVID-19. Zeb et al. [19] considered the model comprised of five compartments. After the emergency of the acute syndrome coronavirus in 2002, which spread to 37 countries, COVID-19 is the third emerging human application purposes disease in the current century and the middle east respiratory syndrome coronavirus in 2012, which spread to 27 countries bilateral lung infiltration including dry cough, fever, trouble, breathing, fatigue and related symptoms caused by [20–26]. Alqarni et al. [27] constructed a new mathematical model for transmission dynamics of COVID-19. The model is based on the data from Saudi Arabia. The authors have discussed the concentration of the disease by introducing first-order ODE into the model. The whole population is divided into five compartments (SEIAR). Analysis of the model is performed by considering the basic reproduction number [28]. Deals fractional form of projectile motion (wind-influenced) is discussed in detail. Inverse singular spectral problems are discussed in [29] by Ozarslan et al. [30] discussed the Lewis model for the soybean drying process using fractional differential operators. The results are compared with Caputo fractional derivative.

One of the most effective ways to uncover the truth about diseases is through mathematical modeling. For the most part, determining differential equations is difficult and does not yield closed-form solutions. To accomplish this, we turned to a variety of numerical methods.

The detail of the rest of the sections is as follows. The Section 2 discusses the COVID-19 epidemic with an isolated compartment and presented the proposed model. Moreover, we prove the positivity of the model in this section also. The basic reproduction number and equilibrium points are discussed in Section 3. Furthermore, the stability of the equilibrium point is discussed in the same section. In the Section 4, numerical schemes are applied, and the stability and consistency of FTCS, Crank Nicolson, and NSFD schemes are investigated. Sections 5 and 6 are respectively presenting the results and conclusion.

2  Proposed Model for COVID-19

We consider the model for COVID-19, which comprises five compartments susceptible, exposed, infected, isolated, and recovered. The compartments are denoted by S,E,I,Q, and R, see [19].

Most infectious disease models comprise ordinary differential equations of the first order. Such models cannot give an accurate picture of disease because of the mobility of the population within the area; therefore, to study pandemic diseases, the spatial content cannot be ignored. Due to the mobility of individuals, the spread of disease may differ from one area to the other.

2.1 Why Diffusion?

The movement of people happens in special regions like countries and cities and generally in the local domain. Let’s look at the movement of people in the United States and China (densely populated regions). The mobility of people may be covering millions of square kilometers, whereas, in small countries, it may reduce to just several square kilometers. To what extent does mobility have larger effects on the variation in ecology. This may variate cultural values and health concerns.

Contrary to typical mathematical models, which deal with the spread and control of epidemics, the present model deals with the greater mobility of humans getting mixed rapidly. The specialty of this model is that it keeps the spatial content into consideration. In the case of pandemics, we cannot ignore the factor of mobility because the disease may spread faster in one area than others due to mobility [31–40]. The following flow map in Fig. 1 shows the different compartment SEIQR of the population and the factors affecting the said compartment.

Figure 1: Flow map for COVID-19 model

∂Sdt=δS∂2S∂x2+A−μS−βNS(E+I)∂Edt=δE∂2E∂x2+βNS(E+I)−πE−(μ+γ)E∂Idt=δI∂2I∂x2+πE−ςI−μI∂Qdt=δQ∂2Q∂x2+γE+ςI−θQ−μQ∂Rdt=δR∂2R∂x2+θQ−μR}.(1)

With initial condition as

S(x,0)≥0,E(x,0)≥0I(x,0)≥0,Q(x,0)≥0R(x,0)≥0}.(2)

As the first four equations are independent of R(t), so we modify the system as follows:

∂Sdt=δS∂2S∂x2+A−μS−βNS(E+I)∂Edt=δE∂2E∂x2+βNS(E+I)−πE−(μ+γ)E∂Idt=δI∂2I∂x2+πE−ςI−μI∂Qdt=δQ∂2Q∂x2+γE+ςI−θQ−μQ}.(3)

Consider the following assumptions

N=Aμ,s=SN,e=EN,i=IN,q=QN.(4)

With initial conditions

s(x,0)=s0≥0,e(x,0)=e0≥0i(x,0)=io≥0,q(x,0)=q0≥0}..(5)

We rewrite the system (3) as follows:

st=δs∂2s∂x2+μ−μs−βNs(e+i);et=δe∂2e∂x2+βNs(e+i)−πe−(μ+γ)e;it=δi∂2i∂x2+πe−ςi−μi;qt=δq∂2q∂x2+γe+ςi−θq−μq;}.(6)

Lemma.1

Under the initial conditions (5), all the solutions of system (6) are non-negative for t≥0.

Proof: By the I. C. s(0)=s0≥0,e(0)=e0≥0,i(0)=i0≥0,q(0)=q0≥0.

st|s=0=δs∂2s∂x2+μ≥0et|e=0=δe∂2e∂x2+βNs(e+i)≥0it|i=0=δi∂2i∂x2+πe≥qt|q=0=δq∂2q∂x2+γe≥0}.(7)

Eq. (7) proves the positivity of the model.

2.2 Basic Reproductive Number

It is the key element in the disease model whose value is a quick check to tell whether the disease has spread or not. The basic reproductive number for the model is given as

R0=βN(ς+μ+π(μ+π+γ)(ς+μ)).(8)

where δs=δe=δi=δq=0.

2.3 Existence and Stability of Equilibrium Points

Equilibrium points are those critical points that yield a constant solution of differential equations. These points may be stable or unstable depending upon the basic reproductive number R0. The system (6) has two equilibrium points as follows:

C0=(so,eo,io,qo),C1=(s1,e1,i1,q1).

where so=1,eo=io=qo=0

s1=ς(π+μ+γ)+μ(π+γ)+μ2βN(π+μ+γ)

e1=βN(π+μ+γ)−ς(π+μ+γ)−μ(ς+μ)βN(ς(π+μ+γ)+π(π+γ)+2πμ+γμ+μ2)

i1=μπ(βN(ς+μ+γ)−ς(π+μ+γ)−2μ(γ+μ)βN(ς2(π+μ+γ)+ςπ2+ςπγ+3ςπμ+2ςγμ+2ςμ2+π2μ+πγμ+2πμ2+γμ2+μ3

q1=μ[(βN(ς+μ+γ)−ς(π+μ+γ)−πμ−γμ−μ2](γπ+γμ+ςπ)βN(ς2(π+μ+γ)+ςπ(π+γ)+3ςπμ+2ςγμ+2ςμ2+π2μ+πγμ+2πμ2+γμ2+μ3(θ+μt)

2.4 System Stability at Equilibrium Point

For the stability of the system, we perturb (6) at the equilibrium point C1=(s1,e1,i1,q1) as under [28].

∂s1∂t=b11s1+b12e1+b13i1+b14q1+d1∂2s1∂x2∂e1∂t=b21s1+b22e1+b23i1+b24q1+d2∂2e1∂x2∂i1∂t=b31s1+b32e1+b33i1+b34q1+d3∂2i1∂x2∂q1∂t=b41s1+b42e1+b43i1+b44q1+d4∂2q1∂x2}.

where

δs=d1, δe=d2, δi=d3, δq=d4.

and

b11=−μ, b12=−βN, b13=−βN, b14=0,

b21=0, b22=−π−μ+βN, b23=βN, b24=0,

b31=0, b32=π, b33=−ς−μ, b34=0,

b41=0, b42=γ, b43=ς, b44=−θ−μ.

Suppose the above system possesses a Fourier solution

s1(x,t)=∑KskeλtCos(kx)e1(x,t)=∑KekeλtCos(kx)i1(x,t)=∑KikeλtCos(kx)q1(x,t)=∑KqkeλtCos(kx)}.

For k=nπ2 where n is a natural number, called the wave number at node n, by substituting the values of s1, e1, i1, q1 in the above equations. We get

λ∑ksk=b11∑Ksk+b12∑Kek+b13∑Kik−d1∑kk2skλ∑kek=b21∑Ksk+b22∑Kek+b23∑Kik−d2∑kk2ekλ∑kik=b31∑Ksk+b32∑Kek+b33∑Kik−d3∑kk2ikλ∑kqk=b41∑Ksk+b42∑Kek+b43∑Kik−d4∑kk2qk.

This leads to the following:

∑K(b11−λ−d1k2)sk+∑Kb12ek+∑Kb13ik+∑kb14qk=0∑Kb21sk+∑K(b22−λ−d2k2)ek+∑Kb23ik+∑kb24qk=0∑Kb31sk+∑Kb32ek+∑K(b33−λ−d3k2)ik+∑kb34qk=0∑Kb41sk+∑Kb42ek+∑Kb43ik+∑k(b44−λ−d4k2)qk=0.

The variational matrix V can be written as

V=(b11−d1k2b12b13b14b21b22−d2k2b23b24b31b32b33−d3k2b34b41b42b43b44−d4k2).

The characteristic polynomial for the above matrix can be written as

P(λ)=F0λ4+F1λ3+F2λ2+F3λ+F4=0.

where

F0=1

F1=−b11−b22−b33−b44+(d1+d2+d3+d4)k2

F2=−b12b21−b13b31−b23b32−b14b41−b24b42−b34b43+b33b44−(b44(d1+d2+d3)+b33(d1+d2+d4))k2+(d3d4+d2(d3+d4)+d1(d2+d3+d4))(k2)2+b22(b33+b44−(d1+d3+d4)k2)+b11(b22+b33+b44−(d2+d3+d4)k2)

F3=b11b23b32−b11b22b33+b14b22b41+b14b33b41−b14b21b42+b11b24b42+b24b33b42−b23b34b42−b14b31b43−b24b32b43+b11b34b43+b22b34b43−b11b22b44+b23b32b44−b11b33b44−b22b33b44+(−b24b42d1−b34b43d1+b33b44d1+b11b33d2−b14b41d2−b34b43d2+b11b44d2+b33b44d2−b14b41d3−b24b42d3+b11b44d3+b11b33d4−b23b32(d1+d4)+b22(b44(d1+d3)+b33(d1+d4)+b11(d3+d4)))k2−((b22d1+b11d2)d3+b44(d1d2+(d1+d2)d3)+(b22(d1+d3)+b11(d2+d3))d4+b33(d2d4+d1(d2+d4)))(k2)2+(d2d3d4+d1(d3d4+d2(d3+d4)))(k2)3+b13(−b21b32−b34b41+b31(b22+b44−(d2+d4)k2))−b12(b23b31+b24b41+b21(−b33−b44+(d3+d4)k2))

F4=b12b24b33b41−b12b23b34b41−b11b24b33b42+b11b23b34b42−b12b24b31b43+b11b24b32b43+b12b21b34b43−b11b22b34b43+b12b23b31b44−b11b23b32b44−b12b21b33b44+b11b22b33b44+(b22b34b43d1−b22b33b44d1+b11b34b43d2−b11b33b44d2+b12b21b44d3−b11b22b44d3+b24(b33b42d1−b32b43d1+(−b12b41+b11b42)d3)+(b12b21−b11b22)b33d4+b23(−b34b42d1−b12b31d4+b32(b44d1+b11d4)))k2+(−b34b43d1d2+(−b24b42d1+b44(b22d1+b11d2))d3−(b23b32d1+(b12b21−b11b22)d3)d4+b33(b44d1d2+(b22d1+b11d2)d4))(k2)2−(b44d1d2d3+(b33d1d2+(b22d1+b11d2)d3)d4)(k2)3+d1d2d3d4(k2)4+b14(b23(b32b41−b31b42)+b21(b33b42−b32b43)+((b33b41−b31b43)d2−b21b42d3)k2−b41d2d3(k2)2+b22(b31b43+b41(−b33+d3k2)))−b13(b24(b32b41−b31b42)+b21(b34b42−b32b44)+((b34b41−b31b44)d2+b21b32d4)k2+b31d2d4(k2)2−b22(b34b41+b31(−b44+d4k2))).

According to Routh Hurwitz’s criterion, the system is stable under the following condition:

F1>0,F1F2−F0F3>0,(F1F2−F0F3)F3−F12F4>0,F4>0.

Theorem: 1 If R0<1 then the system (6) is globally stable.

Proof: We construct the Lyapunov function as

L=(e−e0)+βNς+μ(i−i0)

L′=et+βNς+μit

L′≤βNe−(π+μ+γ)e+βN(ς+μ)(πe)+δe∂2e∂x2+βNς+μ+δi∂2i∂x2a

L′≤βNe−(π+μ+γ)e(ς+μ)(πe)≤(R∘−1)e

This implies that L′<1 when Ro<1, therefore, the system is globally stable.

3  Numerical Methods

In this section, we construct numerical methods for system (6).

st=δs∂2s∂x2+μ−μs−βNs(e+i);et=δe∂2e∂x2+βNs(e+i)−πe−(μ+γ)e;it=δi∂2i∂x2+πe−ςi−μi;qt=δq∂2q∂x2+γe+ςi−θq−μq;}.(9)

For t≥0, x∈[0,L] and the initial conditions are

s(x,0)=s0≥0,e(x,0)=e0≥0i(x,0)=io≥0,q(x,0)=q0≥0}.(10)

Dividing [0,L]×[0,T] into G×P grid points. Define the step size as h=LP & t=TG

xj=jh,j=0,1,2,…,Ptg=gt,g=0,1,2,…,G.

We denote the values sjn, ejn,ijn and qjn as finite-difference approximate values for s(jh,nt),e(jh,nt),i(jh,nt),q(jh,nt), respectively.

3.1 FTCS Scheme

The FTCS scheme for the above system can be written as

sjn+1=sjn+δsΔtΔx2(sj+1n−2sjn+sj−1n)+μΔt−μΔtsjn−βΔtNsjn(ejn+ijn)sjn+1=sjn+d1(sj+1n−2sjn+sj−1n)+μΔt−μΔtsjn−βΔtNsjn(ejn+ijn)sjn+1=sjn(1−2d1−μΔt−βNΔt(ejnijn))+d1(sj+1n+sj−1n)−μΔtejn+1=ejn(1−2d2+βNΔtsjn−πΔt−Δt(μ+γ))+d2(ej+1n+ej−1n)+ΔtβNsjnijnijn+1=ijn(1−2d3−ςΔt−−μΔt)+d3(ij+1n+ij−1n)+πΔtejnqjn+1=qjn(1−2d4+Δtθ−−μΔt)+d4(qj+1n+qj−1n)+γΔtejn+ςΔtijn}

where

d1=δsΔtΔx2,d2=δeΔtΔx2,d3=δiΔtΔx2,d4=δqΔtΔx2(11)

3.2 Stability of FTCS

We carry out stability using Von Neumann stability analysis.

Consider the following:

sjn+1=sjn+δsΔtΔx2(sj+1n−2sjn+sj−1n)+Δt(μ−sjn−βnsjn(ejn+ijn)).

Using the following values:

snj=eαnkeι~ϕjh,sn+1j=eα(n+1)keι~ϕjh,

snj+1=eαnkeι~ϕ(j+1)h,snj−1=eαnkeι~ϕ(j−1)h.

where i~=−1

eαk=1+d1(eι~φh−2+e−i~φh)−(μ+βNlm)Δt

|eαk|=|1+2d1(cos⁡φh−1)−(μ+βNlm)Δt|.

For

cos⁡φh=−1|eαk|=|1−4d1−Δt(μ+βNlm)|−1≤1−4d1−Δt(μ+βNlm)≤1.

After some simplification, we get

−μ+βNlm4≤d1≤2−μ−βNlm4.

Similarly,

ejn+1=ejn+d2(ej+1n−2ejn+ej−1n)+Δt(−(μ+γ)ejn−πejn+βNsjn(ejn+ijn)).Putting values of ejn,ejn+1,ej+1n,ej−1n,We get the following relation after some simplification:eαk=1+d2(2cos⁡φh−2)+(βNa−π−(μ+γ))Δt|eαk|=|1+2d2(cos⁡φh−1)+(βNa−π−(μ+γ))Δt||eαk|=|1−4d2−Δt(μ+γ)+βNa−π)|.

For zero linear constants

eαnk=|1−4d2−(π+μ+γ)|≤1.

Now similarly,

ijn+1=ijn+d3(ij+1n−2ijn+ij−1n)+(πejn−(ς+μ)ijn)Δt.Usingvaluesofijn+1,ij+1n,ijn,ij−1n,eαk=1+2d3(cos⁡φh−1)−(μ+ς)Δt|eαk|=|1+2d3(cos⁡φh−1)−(μ+ς)Δt|.

For

cos⁡φh=−1|eαk|=|1−4d3−Δt(ς+μ)|−1≤1−4d3−Δt(μ+ς)≤1.

After some simplification, we have

−(μ+ς)Δt4≤d3≤2−(μ+ς)Δt4.

qjn+1=qjn+d4(qj+1n−2qjn+qj−1n)+(γejn−ςijn−(θ+μ)qjn)Δt.Usingvaluesofqjn+1,qj+1n,qjn,qj−1n,etc.eαk=1+2d4(cos⁡φh−1)−(μ+θ)Δt|eαk|=|1+2d4(cos⁡φh−1)−(μ+θ)Δt|.

For

cos⁡φh=−1|eαk|=|1−4d4−Δt(θ+μ)|−1≤1−4d4−Δt(μ+θ)≤1.