A rigorous theoretical and numerical analysis of a nonlinear reaction-diffusion epidemic model pertaining dynamics of … – Nature.com

This section presents qualitative analysis of the reaction-diffusion COVID-19 compartmental epidemic model (1). We proceed to prove the basic mathematical properties of the model solution as follows.

One of the most important properties of an epidemic model is the solution boundedness. We take into consideration the approach described in28 in order to analyze the solution boundedness of the problem (1). The result is given in the following Theorem.

The solution of the model (1) i.e., ((S(.,tilde{t}), E(.,tilde{t}), I_1(.,tilde{t}), I_2(.,tilde{t}), I_3(.,tilde{t}), R(.,tilde{t}))) is bounded (forall) (tilde{t}ge 0).

In order to prove the desired result, add all equations in the model (1)

$$begin{aligned}&dfrac{partial }{partial tilde{t}}S(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}E(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}I_1(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}I_2(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}I_3(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}R(tilde{t},tilde{y}), \&= D_1dfrac{partial ^2}{partial tilde{y}^2}S(tilde{t},tilde{y})+D_2dfrac{partial ^2}{partial tilde{y}^2}E(tilde{y}, t)+D_3dfrac{partial ^2}{partial tilde{y}^2}I_1(tilde{t},tilde{y})+D_4dfrac{partial ^2}{partial tilde{y}^2}I_2(tilde{t},tilde{y})+D_5dfrac{partial ^2}{partial tilde{y}^2}I_3(tilde{t},tilde{y})\&quad +D_6dfrac{partial ^2}{partial tilde{y}^2}R(tilde{t},tilde{y})+Pi -zeta (S(tilde{t},tilde{y})+E(tilde{t},tilde{y})+I_1(tilde{t},tilde{y})+I_2(tilde{t},tilde{y})+I_3(tilde{t},tilde{y})+R(tilde{t},tilde{y}))\&quad -xi _{1}I_1(tilde{t},tilde{y})-xi _{2}I_2(tilde{t},tilde{y}). end{aligned}$$

Integrating over (Lambda), using the well known Divergence theorem29 and make use of no flux boundary conditions, which yield to

$$begin{aligned}&int _{Lambda }bigg {dfrac{partial }{partial tilde{t}}S(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}E(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}I_1(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}I_2(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}I_3(tilde{t},tilde{y})+dfrac{partial }{partial tilde{t}}R(tilde{t},tilde{y}) bigg }dtilde{y},\&=Pi |Lambda |-dint _{Lambda }bigg {(S(tilde{t},tilde{y})+E(tilde{t},tilde{y})+I_1(tilde{t},tilde{y})+I_2(tilde{t},tilde{y})+I_3(tilde{t},tilde{y})+R(tilde{t},tilde{y}))bigg }dtilde{y}\&quad -int _{Lambda }bigg {xi _{1}I_1(tilde{t},tilde{y})+xi _{2}I_2(tilde{t},tilde{y})bigg }dtilde{y},\&le Pi |Lambda |-zeta N(tilde{t}).\ dfrac{d}{dtilde{t}}N(tilde{t})&=Pi |Lambda |-zeta N(tilde{t}). end{aligned}$$

It gives

$$begin{aligned}0le N(tilde{t})le dfrac{Pi |Lambda |}{zeta }-exp (-zeta tilde{t}) N(0), quad ;forall ; tilde{t}ge 0.end{aligned}$$

Hence

$$begin{aligned} lim limits _{trightarrow +infty } N(t)le dfrac{Pi |Lambda |}{zeta }. end{aligned}$$

(square)

A positively invariant set for the system (1) is defined as follows:

$$begin{aligned} Phi = bigg {(S(tilde{t},tilde{y}), E(tilde{t},tilde{y}), I_1(tilde{t},tilde{y}), I_2(tilde{t},tilde{y}), I_3(tilde{t},tilde{y}), R(tilde{t},tilde{y})^{T}in mathbb {R}_+^6: N(tilde{t})le dfrac{Pi |Lambda |}{d}bigg }subset mathbb {R}_+^6. end{aligned}$$

For derivation of the threshold parameter known as the basic reproductive number, we used the well-known approached considered in30. Model (1) possess two equilibria i.e. the disease-free equilibrium (DFE) and the endemic equilibrium (EE) represented by (Xi _0) and (Gamma _{EE}) respectively such that:

$$begin{aligned} Gamma _0 = (S_0, E_0, I_{1_0}, I_{2_0}, I_{3_0}, R_0)= left( Pi /zeta , 0, 0, 0, 0, 0right) . end{aligned}$$

The EE is calculated as follows:

$$begin{aligned} Gamma _{EE} = (S^*, E^*, {I^*_{1}}, {I^*_{2}}, {I^*_{3}}, R^*), end{aligned}$$

with the following analytical values

$$begin{aligned} {left{ begin{array}{ll} \ S^{*} = dfrac{Pi text{c{N}}}{text{c{N}}zeta +mathfrak {R}^0-1}, \ \ E^{*} = dfrac{Pi (mathfrak {R}^0-1)}{alpha _1left( text{c{N}}zeta +mathfrak {R}^0-1 right) },\ \ I_1^{*} = dfrac{Pi mathfrak {R}_{01}(mathfrak {R}^0-1)}{beta left( text{c{N}}zeta +mathfrak {R}^0-1 right) },\ \ I_2^{*} = dfrac{Pi mathfrak {R}_{01}(mathfrak {R}^0-1)}{beta _Pleft( text{c{N}}zeta +mathfrak {R}^0-1 right) },\ \ I_3^{*} = dfrac{Pi mathfrak {R}_{03}(mathfrak {R}^0-1)}{beta psi left( text{c{N}}zeta +mathfrak {R}^0-1right) }, \ \ R^* = dfrac{Pi }{zeta }left[ eta _1dfrac{mathfrak {R}_{01}}{beta }+eta _2dfrac{mathfrak {R}_{02}}{beta _P}+eta _3dfrac{mathfrak {R}_{03}}{beta psi } right] dfrac{(mathfrak {R}^0-1)}{left( text{c{N}}zeta +mathfrak {R}^0-1right) }, end{array}right. } end{aligned}$$

and

$$begin{aligned} text{c{N}} = frac{1}{r+zeta }+left( 1+frac{eta _1}{zeta }right) frac{mathfrak {R}_{01}}{beta }+left( 1+frac{eta _2}{zeta }right) frac{mathfrak {R}_{01}}{beta _P}+left( 1+frac{eta _3}{zeta }right) frac{mathfrak {R}_{03}}{beta }. end{aligned}$$

Moreover, the basic reproductive number (mathfrak {R}^0) is computed as follows:

The infectious classes in the proposed model (1) are (E, I_1, I_2) and (I_3), Henceforth, the vectors below present the transmission of newborn infections and the transitions between various classes.

$$begin{aligned} {textbf {F}}= left( begin{array}{c} beta frac{(I_1+psi I_3)S}{N}+beta _Pfrac{ I_2 S}{N}\ 0\ 0\ 0\ end{array}right) , quad quad {textbf {V}} = left( begin{array}{c} (r+zeta )E\ -r k_1 E+(eta _1+zeta +zeta _{1})I_1\ -r k_2 E+(eta _2+zeta +zeta _{2})I_2\ -r(1-k_1-k_2) E+(eta _3+zeta )I_3 end{array} right) , end{aligned}$$

The Jacobian of above matrices are evaluated as:

$$begin{aligned} mathcal {F}= & {} left( begin{array}{cccc} 0 &{} beta &{} beta _mathcal {P} &{} psi beta \ 0 &{} 0 &{} 0 &{} 0\ 0 &{} 0 &{} 0 &{} 0\ 0 &{} 0 &{} 0 &{} 0\ end{array}right) ,\ mathcal {V}= & {} left( begin{array}{cccc} r+zeta &{} 0 &{} 0 &{} 0\ -r k_1 &{} eta _1+zeta _1+zeta &{} 0 &{} 0\ -r k_1 &{} 0 &{} eta _2+zeta _2+zeta &{} 0\ -r(1-k_1-k_2) &{} 0 &{} 0 &{} eta _3+zeta \ end{array}right) , end{aligned}$$

and (N = dfrac{Pi }{zeta }) in case of disease-free equilibrium. Thus, the associated next generation matrix is,

$$begin{aligned} mathcal {F}mathcal {V}^{-1} = left( begin{array}{cccc} frac{rbeta psi (1-k_1-k_2)}{(r+zeta )(eta _3+zeta )}+frac{rbeta k_1 }{(r+zeta )(eta _1+zeta _{1}+zeta )}+frac{rbeta _P k_2 }{(r+zeta )(eta _2+zeta _{2}+zeta )} &{} frac{beta }{(eta _1+zeta _{1}+zeta )} &{} frac{beta _P }{(eta _2+zeta _{2}+zeta )} &{} frac{beta psi }{(eta _3+zeta )}\ 0 &{} 0 &{} 0 &{} 0\ 0 &{} 0 &{} 0 &{} 0\ 0 &{} 0 &{} 0 &{} 0 end{array} right) . end{aligned}$$

The basic reproductive number which is the spectral radius of (FV^{-1}) and is given by:

$$begin{aligned} mathfrak {R}^0 = rho left( FV^{-1}right) = mathfrak {R}_{01}+mathfrak {R}_{02}+mathfrak {R}_{03}, end{aligned}$$

(5)

where

$$begin{aligned} mathfrak {R}_{01} =frac{rbeta k_1 }{(r+zeta )(eta _1+zeta +zeta _{1})}, mathfrak {R}_{02}= frac{rbeta _P k_2 }{(r+zeta )(eta _2+zeta +zeta _{2})} text {and} mathfrak {R}_{03}= frac{rbeta psi (1-k_1-k_2)}{(r+zeta )(eta _3+zeta )}. end{aligned}$$

The DFE (Gamma _0) is stable locally asymptomatically in (Phi) if (mathfrak {R}^0<1), otherwise it is unstable.

The Jacobian of (1) at the DFE point (Gamma _0) is

$$begin{aligned} J(Gamma _{0}) = left( begin{array}{cccccc} -zeta &{} 0 &{} -beta &{} -beta _P &{} -beta psi &{} 0\ 0 &{} -r-zeta &{} beta &{} beta _P &{} beta psi &{} 0\ 0 &{} r k_1 &{} -eta _{1}-zeta -zeta _1 &{} 0 &{} 0 &{} 0 \ 0 &{} r k_2 &{} 0 &{} -eta _{2}-zeta -zeta _2 &{} 0 &{} 0 \ 0 &{} r(1-k_{1}-k_{2}) &{} 0 &{} 0 &{} -zeta -eta _3 &{} 0 \ 0 &{} 0 &{} eta _1 &{} eta _2 &{} eta _3 &{} -zeta end{array}right) , end{aligned}$$

and the characteristic polynomial associated to above matrix is given as follows:

$$begin{aligned} P(varphi ) = (zeta +varphi )^2(varphi ^{3}+c_{2}varphi ^{2}+c_{1}varphi +c_{0}), end{aligned}$$

(6)

where,

$$begin{aligned} c_{2}&= b_1l_2(1-mathfrak {R}_{01})+l_1l_3(1-mathfrak {R}_{02})+l_1l_4(1-mathfrak {R}_{03})+l_2l_3+l_2l_4+l_3l_4,\ c_{1}&= l_1l_2(l_3+l_4)(1-mathfrak {R}_{01})+l_1l_3l_4(1-mathfrak {R}_{03}-mathfrak {R}_{02})+l_1l_4(l_3-l_2)mathfrak {R}_{03}\&quad -l_1l_2l_3mathfrak {R}_{02}, \ c_{0}&=(r+zeta )(eta _3+zeta )(eta _1+zeta _1+zeta )(eta _2+zeta _2+zeta ) left( 1-mathfrak {R}^0 right) >0 text {for} mathfrak {R}^0<1, end{aligned}$$

The values of (c_0, c_1, c_2) in (6) are claimed to be positive under the condition (mathfrak {R}^0<1). Further, (c_1c_2-c_3>0) confirming the Routh-Hurwitz conditions. Hence, the (texttt {DFE}) stable locally when (mathfrak {R}^0<1). (square)

To discuss the local stability of endemic equilibria (Gamma _{EE}), the model (1) is linearized at (Gamma _{EE} = (S^*, E^*, I^*_{1}, I^*_{2}, I^*_{3}, R^* )). For this purpose assume that

$$begin{aligned} left. begin{array}{ll} \ S(tilde{t},tilde{y}) &{}= bar{S}(tilde{t},tilde{y})+S^*,\ \ E(tilde{t},tilde{y}) &{}= bar{E}(tilde{t},tilde{y})+E^*,\ \ I_1(tilde{t},tilde{y}) &{}= bar{I}_{1}(tilde{t},tilde{y})+I^*_{1},\ \ I_2(tilde{t},tilde{y}) &{}= bar{I}_{2}(tilde{t},tilde{y})+I^*_{2},\ \ I_3(tilde{t},tilde{y}) &{}= bar{I}_{3}(tilde{t},tilde{y})+I^*_{3},\ \ R(tilde{t},tilde{y}) &{}= bar{R}(tilde{t},tilde{y})+R^*. end{array} right} end{aligned}$$

(7)

(bar{S}(tilde{t},tilde{y}), bar{E}(tilde{t},tilde{y}), bar{I}_{1}(tilde{t},tilde{y}), bar{I}_{2}(tilde{t},tilde{y}), bar{I}_{3}(tilde{t},tilde{y})) and (bar{R}(tilde{t},tilde{y})) are minimal perturbation. The linearized form of the problem (1) is given by (8),

$$begin{aligned} left. begin{array}{ll} \ dfrac{partial bar{S}}{partial tilde{t}} &{}= D_1dfrac{partial ^2 bar{S}}{partial tilde{y}^2}+g_{11}bar{S}+g_{12}bar{E}+g_{13}bar{I}_{1}+g_{14}bar{I}_2+g_{15}bar{I}_3+g_{16}bar{R},\ \ dfrac{partial bar{E}}{partial tilde{t}} &{}= D_2dfrac{partial ^2 bar{E}}{partial tilde{y}^2}+g_{21}bar{S}+g_{22}bar{E}+g_{23}bar{I}_{1}+g_{24}bar{I}_2+g_{25}bar{I}_3+g_{26}bar{R},\ \ dfrac{partial bar{I}_1}{partial tilde{t}} &{}= D_3dfrac{partial ^2 bar{I}_1}{partial tilde{y}^2}+g_{31}bar{S}+g_{32}bar{E}+g_{33}bar{I}_{1}+g_{34}bar{I}_2+g_{35}bar{I}_3+g_{36}bar{R},\ \ dfrac{partial bar{I}_2}{partial tilde{t}} &{}= D_4dfrac{partial ^2 bar{I}_2}{partial tilde{y}^2}+g_{41}bar{S}+g_{42}bar{E}+g_{43}bar{I}_{1}+g_{44}bar{I}_2+g_{45}bar{I}_3+g_{46}bar{R},\ \ dfrac{partial bar{I}_3}{partial tilde{t}} &{}= D_5dfrac{partial ^2 bar{I}_3}{partial tilde{y}^2}+g_{51}bar{S}+g_{52}bar{E}+g_{53}bar{I}_{1}+g_{54}bar{I}_2+g_{55}bar{I}_3+g_{56}bar{R},\ \ dfrac{partial bar{R}}{partial tilde{t}} &{}= D_6dfrac{partial ^2 bar{R}}{partial tilde{y}^2}+g_{61}bar{S}+g_{62}bar{E}+g_{63}bar{I}_{1}+g_{64}bar{I}_2+g_{65}bar{I}_3+g_{66}bar{R}, end{array} right} end{aligned}$$

(8)

such that

$$begin{aligned} begin{array}{ll} &{}g_{11} = -dfrac{beta }{N}left( I^*_{1}+psi I^*_{3}right) -dfrac{beta _P}{N}I^*_{2}-zeta , g_{12} = 0, g_{13} = -dfrac{beta }{N}S^*, g_{14} = -dfrac{beta _P}{N}S^*, \ &{}g_{15} = -dfrac{beta psi _{1}}{N}S^*,hspace{0.3cm} g_{16} = 0,hspace{0.3cm} g_{21} = dfrac{beta }{N}left( I^*_{1}+psi I^*_{3}right) +dfrac{beta _P}{N}I^*_{2}, hspace{0.3cm} g_{22} = -(r+zeta ), \ &{} g_{23} = dfrac{beta }{N}S^*, hspace{0.3cm} g_{24} = dfrac{beta _P}{N}S^*, hspace{0.3cm} g_{25} = dfrac{beta psi }{N}S^*, hspace{0.3cm} g_{26} = 0, hspace{0.3cm} g_{31} = 0, hspace{0.3cm} g_{32} = r k_1, \ &{}g_{33} = -(eta _{1}+zeta +zeta _1), hspace{0.3cm} g_{34} = 0, hspace{0.3cm} g_{35} = 0, hspace{0.3cm} g_{36} = 0, hspace{0.3cm} g_{41} = 0, hspace{0.3cm} g_{42} = r k_2,\ &{}g_{43} = 0, hspace{0.3cm} g_{44} = -(eta _{2}+zeta +zeta _2), hspace{0.3cm} g_{45}=0, hspace{0.3cm} g_{46}=0, hspace{0.3cm} g_{51} = 0, hspace{0.3cm} g_{56} = 0,\ &{}g_{52} = r(1-k_1-k_2), hspace{0.3cm} g_{53} = 0, g_{54} = 0, g_{55} = -(eta _{3}+zeta ), hspace{0.3cm}g_{61} = 0,\ &{}g_{62} = 0, hspace{0.3cm} g_{55} = -(eta _{3}+zeta ), hspace{0.3cm} g_{63} = eta _1, hspace{0.3cm} g_{64}=eta _2, hspace{0.3cm} g_{65}=eta _3, hspace{0.3cm} g_{66} = -zeta . end{array} end{aligned}$$

Given that the linearized system (8) has a solution in Fourier series form, then

$$begin{aligned} left. begin{array}{ll} \ bar{S}(tilde{t},tilde{y}) &{}= sum limits _{k}e^{lambda t}b_{1k}cos (ktilde{y}), \ \ bar{E}(tilde{t},tilde{y}) &{}= sum limits _{k}e^{lambda t}b_{2k}cos (ktilde{y}), \ \ bar{I}_{1}(tilde{t},tilde{y}) &{}= sum limits _{k}e^{lambda t}b_{3k}cos (ktilde{y}), \ \ bar{I}_{2}(tilde{t},tilde{y}) &{}= sum limits _{k}e^{lambda t}b_{4k}cos (ktilde{y}), \ \ bar{I}_{3}(tilde{t},tilde{y}) &{}= sum limits _{k}e^{lambda t}b_{5k}cos (ktilde{y}), \ \ bar{R}(tilde{t},tilde{y}) &{}= sum limits _{k}e^{lambda t}b_{6k}cos (ktilde{y}). end{array} right} end{aligned}$$

(9)

In above, (k =frac{npi }{2}), with (nin Z^{+}), indicates wave-number for the node n. Using (9) in (8) yield to,

$$begin{aligned} left. begin{array}{ll} \ sum limits _{k}(g_{11}-k^2D_1-lambda )b_{1k}+sum limits _{k}g_{12}b_{2k} +sum limits _{k}g_{13}b_{3k}+sum limits _{k}g_{14}b_{4k}+sum limits _{k}g_{15}b_{5k}+sum limits _{k}g_{16}b_{6k}&{}=0, \ \ sum limits _{k}g_{21}b_{1k}+sum limits _{k}(g_{22}-k^2D_{2}-lambda )b_{2k}+sum limits _{k}g_{23}b_{3k}+sum limits _{k}g_{24}b_{4k}+sum limits _{k}g_{25}b_{5k}+sum limits _{k}g_{26}b_{6k}&{}=0, \ \ sum limits _{k}g_{31}b_{1k}+sum limits _{k}g_{32}b_{2k}+sum limits _{k}(g_{33}-k^2D_{I_1}-lambda )b_{3k}+sum limits _{k}g_{34}b_{4k}+sum limits _{k}g_{35}b_{5k}+sum limits _{k}g_{36}b_{6k} &{}=0, \ \ sum limits _{k}g_{41}b_{1k}+sum limits _{k}g_{42}b_{2k}+sum limits _{k}g_{43}b_{3k}+sum limits _{k}(g_{44}-k^2D_4-lambda )b_{4k}+sum limits _{k}g_{45}b_{5k}+sum limits _{k}g_{46}b_{6k} &{}=0, \ \ sum limits _{k}g_{51}b_{1k}+sum limits _{k}g_{52}b_{2k}+sum limits _{k}g_{53}b_{3k}+sum limits _{k}g_{54}b_{4k}+sum limits _{k}(g_{44}-k^2D_5-lambda )b_{5k}+sum limits _{k}g_{56}b_{6k} &{}=0, \ \ sum limits _{k}g_{61}b_{1k}+sum limits _{k}g_{62}b_{2k}+sum limits _{k}g_{63}b_{3k}+sum limits _{k}g_{64}b_{4k}+sum limits _{k}g_{65}b_{5k}+sum limits _{k}(g_{66}-k^2D_{R}-lambda )b_{6k} &{}=0. end{array} right} end{aligned}$$

(10)

The matrix V representing the variational matrix of (8) is

$$begin{aligned} V = left( begin{array}{cccccc} -c_{11} &{} 0 &{} -g_{13}&{} -g_{14}&{} -g_{15}&{} 0\ g_{21} &{} -c_{22}&{} g_{23} &{} g_{24}&{} g_{25}&{} 0\ 0 &{}g_{23}&{}-c_{33}&{} 0&{} 0&{} 0\ 0 &{}g_{42}&{}0&{}-c_{44}&{} 0&{} 0\ 0 &{}g_{52}&{}0&{} 0&{}-c_{55}&{} 0\ 0 &{}0&{}g_{63}&{} g_{64}&{} g_{65}&{}-c_{66}\ end{array}right) , end{aligned}$$

(11)

where,

$$begin{aligned} c_{11}&= k^2D_1+g_{11}, \ c_{22}&= k^2D_{2}+g_{22}, \ c_{33}&= k^2D_3+g_{33}, \ c_{44}&= k^2D_4+g_{44},\ c_{55}&= k^2D_5+g_{55}, \ c_{66}&= k^2D_{R}+g_{66}. end{aligned}$$

The subsequent polynomial is,

$$begin{aligned} P(lambda ) = (lambda +c_{66})(lambda ^{5}+mathcal {A}_{4}lambda ^{4}+mathcal {A}_{3}lambda ^{3}+mathcal {A}_{2}lambda ^{2}+mathcal {A}_{1}lambda +mathcal {A}_{0}). end{aligned}$$

(12)

The relative coefficient values are;

$$begin{aligned} mathcal {A}_{4}&= c_{11}+k^{2}(D_2+D_3+D_4+D_5+D_6)+g_{11}+g_{22}+g_{33}+g_{44},\ mathcal {A}_{3}&= b_1+g_{22}left( g_{33}left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}right) +g_{44}left( 1-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +g_{55}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}right) right) ,\ mathcal {A}_{2}&=b_2+b_3+c_{11}g_{22}left( g_{33}left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}right) +g_{44}left( 1-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +g_{55}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}right) right) +g_{22}\&quad times g_{33}(D_5+D_4)k^{2}left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}right) +g_{22}g_{44}(D_5+D_3)k^{2}left( 1-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +g_{22}g_{55}k^{2}\&quad times (D_3+D_4)left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}right) +g_{22}g_{33}g_{44}left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0} right) +g_{22}g_{33}g_{55}\&quad times left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0} right) +g_{22}g_{44}g_{55}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +g_{21}g_{22}g_{33}dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}+ g_{21}g_{22}\&quad times g_{55}dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}+g_{21}g_{22}g_{44}dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}, \ mathcal {A}_{1}&= b_4+g_{22}g_{33}D_5D_4k^{4}left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}right) +c_{11}left( g_{22}g_{33}(D_5+D_4)k^{2}left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}right) right) \&quad +g_{22}g_{44}D_5D_3k^{4}left( 1-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +c_{11}left( g_{22}g_{44}(D_5+D_3)k^{2}left( 1-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) right) +c_{11}\&quad times left( g_{22}g_{55}(D_4+D_3)k^{2}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}right) right) +c_{11}g_{22}g_{33}g_{44}left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +c_{11} \&quad times g_{22}g_{55}left( g_{33}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}right) +g_{44}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) right) +g_{22}g_{33}g_{44}D_5k^{2}\&quad times left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +g_{22}g_{55}D_4D_3k^{4}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}right) +g_{22}g_{33}g_{55}D_4k^{2}\&quad times left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}right) +g_{22}g_{44}g_{55}D_3k^{2}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +g_{21}(c_{44}+c_{55})g_{22}\&quad times g_{33}dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}+g_{21}(c_{33}+c_{44})g_{22}g_{55}dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}+g_{21}(c_{33}+c_{55})g_{22}g_{44}dfrac{mathfrak {R}_{02}}{mathfrak {R}^0},\ mathcal {A}_{0}&= b_5 +c_{11}g_{22}g_{44}D_5D_3k^{4}left( 1-dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}right) +c_{11}g_{22} g_{55}D_4D_3k^{4}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}right) +c_{11}g_{22}\&quad times g_{33}g_{44}D_5k^{2}left( 1-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) +c_{11}g_{22}g_{33}g_{55}D_4k^{2}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}right) +c_{11}g_{22}g_{44}g_{55}D_3k^{2}\&quad times left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}right) +c_{11}g_{22}g_{33}g_{44}g_{55}left( 1-dfrac{mathfrak {R}_{03}}{mathfrak {R}^0}-dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) -c_{11}g_{22}g_{33}D_5D_4k^{4}dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}\&quad -c_{11}g_{22}g_{33}g_{44}k^{2}left( D_4dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}+D_3dfrac{mathfrak {R}_{02}}{mathfrak {R}^0}right) -c_{11}c_{55}g_{22}g_{44}^{2}dfrac{mathfrak {R}_{01}}{mathfrak {R}^0}. end{aligned}$$

The values of (b_i), where (i=1,ldots ,5) are given below:

$$begin{aligned} b_1&= c_{11} (c_{22}+c_{33}+c_{44}+c_{55})+g_{33}left( D_2k^{2}+c_{44}+c_{55}right) +g_{44}k^{2}left( D_2+D_3+D_5right) \&quad +g_{55}left( (D_2+D_3)k^{2}+c_{44}right) +k^{4}left( D_5(D_4+D_{I_E})+(D_5+D_4)(D_3+D_{I_E})right) \&quad +g_{22}k^{2}left( D_3+D_4+D_5right) ,\ b_2&=c_{11}left( k^{4}left( D_5(D_4+D_{I_E})+(D_5+D_4)(D_3+D_{I_E})right) +g_{44}k^{2}left( D_2+D_5+D_3right) right) \&quad +c_{11}g_{22}k^{2}left( D_3+D_4+D_5right) ,\ b_3&=c_{11}left( g_{44}((D_2+D_3)k^{2}+c_{44})+g_{33}(D_2k^{2}+c_{44}+c_{55})right) +(D_4D_5(D_2+D_3))k^{6}\&quad +(D_2D_3(D_5+D_4))k^{6}+g_{22}(D_4D_4+D_5(D_4+D_3))k^{2}+g_{33}g_{44}k^{2}(D_2+D_5)\&quad +g_{33}(D_5D_4+D_{I_E}(D_5+D_4))k^{2}+g_{44}k^{4}(D_5D_3+D_{I_E}(D_5+D_3))\&quad +g_{55}k^{4}(D_2D_4+D_3(D_2+D_4))+g_{44}g_{55}k^{2}(D_2+D_3)+g_{33}g_{55}(D_2k^{2}+c_{44}),\ b_4&=c_{11}left( D_5D_4(D_2+D_3)+D_2D_3(D_5+D_4)k^{6}+left( (g_{22}+g_{55})D_4D_3right) k^{4} right) \&quad +c_{11}left( (D_4+D_3)(g_{22}D_5+g_{55}D_2)k^{4}+g_{44}left( D_5D_3+D_{I_E}(D_5+D_3)right) k^{4} right) \&quad +c_{11}left( g_{44}left( D_5D_4+D_{I_E}(D_5+D_4)right) k^{4}+g_{55}(g_{33}(D_2+D_4)+g_{44}(D_2+D_5))k^{2} right) \&quad +c_{11}left( g_{44}g_{55}+2(D_2+D_3)k^{2}right) +g_{22}D_3D_5D_4k^{6}+c_{33}c_{44}c_{55}D_2k^{2},\ b_5&=g_{21}c_{22}c_{33}c_{44}c_{55}+c_{11}left( c_{33}k^{2}(g_{22}D_5D_5k^{2})+D_2(c_{44}+c_{55}) right) . end{aligned}$$

The coefficients (mathcal {A}_{i}), (i = 0,ldots , 4), of P given by (13) are positive, if (mathfrak {R}^0>1). Further, it also satisfies the Routh-Hurwitz stability conditions for a polynomial having degree five, i.e.,

$$begin{aligned}&mathcal {A}_{4}mathcal {A}_{3}mathcal {A}_{2}-mathcal {A}_{2}^{2}-mathcal {A}_{4}^{2}mathcal {A}_{1}>0,\&left( mathcal {A}_{4}mathcal {A}_{1}-mathcal {A}_{0}right) left( mathcal {A}_{4}mathcal {A}_{3}mathcal {A}_{2}-mathcal {A}_{2}^2-mathcal {A}_{4}^2mathcal {A}_{1}right) -mathcal {A}_{0}(mathcal {A}_{4}mathcal {A}_{3}mathcal {A}_{2})^2-mathcal {A}_{4}mathcal {A}_{0}^2>0. end{aligned}$$

Therefore, it is asserted that (Gamma _{EE}) is stable locally for (mathfrak {R}^{0}>1).

In the subsequent section, we will investigate the stability of the problem (1) in global case at the steady state (Gamma _{0}) using a nonlinear Lyapunov stability approach.

If (mathfrak {R}^0<1) then the DFE (Gamma _{0}) of the system (1) is globally stable in (Phi).

The Lyapunov-type function is define as

$$begin{aligned} V(tilde{t}) = int _{Lambda }biggl {E(tilde{t},tilde{y})+j_1I_1(tilde{t},tilde{y})+j_2I_2(tilde{t},tilde{y})+j_3I_3(tilde{t},tilde{y})biggr }dtilde{y}, end{aligned}$$

with

$$begin{aligned} j_1 = dfrac{beta }{eta _{1}+zeta +zeta _1}, j_2 = dfrac{beta }{eta _{2}+zeta +zeta _2}, hbox {and} j_3 = dfrac{psi beta }{eta _{1}+zeta }. end{aligned}$$

Differentiating (V(tilde{t},tilde{y})) with the solution of (1) as follows:

$$begin{aligned} dfrac{dV}{dtilde{t}}&= int _{Lambda }biggl {D_{2}dfrac{partial ^2 E(tilde{t},tilde{y})}{partial tilde{y}^2}+j_1D_{3}dfrac{partial ^2 I_1(tilde{t},tilde{y})}{partial tilde{y}^2}+j_2D_{4}dfrac{partial ^2 I_2(tilde{t},tilde{y})}{partial tilde{y}^2}+j_3D_{5}dfrac{partial ^2 I_2(tilde{t},tilde{y})}{partial tilde{y}^2}biggr }dx\&quad +int _{Lambda }biggl {j_1r k_1+j_2r k_2+j_3r(1-k_1-k_2)-(r+zeta ) biggr }E(tilde{t},tilde{y})dtilde{y}\&quad +int _{Lambda }biggl {beta frac{I_1(tilde{t},tilde{y})S(tilde{t},tilde{y})}{N(tilde{t},tilde{y})}+beta psi frac{I_3(tilde{t},tilde{y})S(tilde{t},tilde{y})}{N(tilde{t},tilde{y})}+beta _Pfrac{I_2(tilde{t},tilde{y})S(tilde{t},tilde{y})}{N(tilde{t},tilde{y})} biggr }dtilde{y}\&quad -int _{Lambda }biggl {j_1(eta _{1}+zeta +zeta _{1})I_1+j_2(eta _{2}+zeta +zeta _{2})I_2+j_3(eta _{3}+zeta )I_3biggr }dtilde{y}. end{aligned}$$

As, (S(tilde{t},tilde{y})le N(tilde{t},tilde{y})) for (tilde{y}in Lambda) and (tilde{t}ge 0), therefore,

$$begin{aligned} dfrac{dV}{dtilde{t}}&le int _{Lambda }biggl {D_{2}dfrac{partial ^2 E(tilde{t},tilde{y})}{partial tilde{y}^2}+D_{3}dfrac{partial ^2 I_1(tilde{t},tilde{y})}{partial tilde{y}^2}+D_{4}dfrac{partial ^2 I_2(tilde{t},tilde{y})}{partial tilde{y}^2}+D_{5}dfrac{partial ^2 I_2(tilde{t},tilde{y})}{partial tilde{y}^2}biggr }dtilde{y}\&quad +int _{Lambda }biggl {j_1r k_1+j_2r k_2+j_3r(1-k_1-k_2)-(r+zeta ) biggr }E(tilde{t},tilde{y})dtilde{y}\&quad +int _{Lambda }biggl {beta I_1(tilde{t},tilde{y})+beta psi I_3(tilde{t},tilde{y})+beta _PI_2(tilde{t},tilde{y}) biggr }dtilde{y}\&quad -int _{Lambda }biggl {j_1(eta _{1}+zeta +zeta _{1})I_1+j_2(eta _{2}+zeta +zeta _{2})I_2+j_3(eta _{3}+zeta )I_3biggr }dtilde{y}, \ dfrac{dV}{dtilde{t}}&le int _{Lambda }biggl {D_{2}dfrac{partial ^2 E(tilde{t},tilde{y})}{partial tilde{y}^2}+D_{3}dfrac{partial ^2 I_1(tilde{t},tilde{y})}{partial tilde{y}^2}+D_{4}dfrac{partial ^2 I_2(tilde{t},tilde{y})}{partial tilde{y}^2}+D_{5}dfrac{partial ^2 I_2(tilde{t},tilde{y})}{partial tilde{y}^2}biggr }dtilde{y}\&quad +int _{Lambda }biggl {j_1r k_1+j_2r k_2+j_3r(1-k_1-k_2)-(r+zeta ) biggr }E(tilde{t},tilde{y})dtilde{y}\&quad +int _{Lambda }biggl {(beta -j_1(eta _{1}+zeta +zeta _{1}))I_1(tilde{t},tilde{y})+(beta psi -j_2(eta _{2}+zeta +zeta _{2}))I_3(tilde{t},tilde{y})biggr }dtilde{y}\&quad +int _{Lambda }biggl {(beta _P-j_3(eta _{3}+zeta ))I_2(tilde{t},tilde{y})biggr }dtilde{y}. end{aligned}$$

Using, the criterions described in (2), we have

$$begin{aligned} dfrac{dV}{dtilde{t} }&le (r+zeta )int _{Lambda }biggl {j_1dfrac{r k_1}{(r+zeta )}+j_2dfrac{r k_2}{(r+zeta )}+j_3dfrac{r(1-k_1-k_2)}{(r+zeta )}-1 biggr }E(tilde{t},tilde{y})dtilde{y}\&le (r+zeta )(mathfrak {R}^0-1)int _{Lambda }E(tilde{t},tilde{y})dtilde{y}. end{aligned}$$

It is obvious (dfrac{dV}{dtilde{t} }<0) (forall) (tilde{t} ge 0) and (tilde{y}in Lambda) (Leftrightarrow) (mathfrak {R}^0<1). Further, (dfrac{dV}{dtilde{t} }=0) if and only if (E(tilde{t},tilde{y})rightarrow 0), then the model (1) implies that (I_1(tilde{t},tilde{y})rightarrow 0), (I_2(tilde{t},tilde{y})rightarrow 0), (I_2(tilde{t},tilde{y})rightarrow 0), (R( t,tilde{y})rightarrow 0) and (S(tilde{t},tilde{y})rightarrow dfrac{Pi }{zeta}) ((S, E, I_1, I_2, I_3, R) rightarrow left( dfrac{Pi }{zeta }, 0, 0, 0, 0, 0right)). In a result we concluded that the largest compact invariant set within the region (left{ (S, E, I_1, I_2, I_3, R): dfrac{dV}{dtilde{t}}=0 right}) is (Gamma _0). By employing the well-established LaSalles principle, (Gamma _0) is globally asymptotically stable under the condition (mathfrak {R}^0<1). (square)

View original post here:

A rigorous theoretical and numerical analysis of a nonlinear reaction-diffusion epidemic model pertaining dynamics of ... - Nature.com