Stability and bifurcation analysis of a nested multi-scale model for COVID-19 viral infection

Let S ( 0 ) ≥ 0 S\left(0)\ge 0 , E ( 0 ) ≥ 0 E\left(0)\ge 0 , and I ( 0 ) ≥ 0 I\left(0)\ge 0 , and then, the solutions S ( t ) , E ( t ) S\left(t),E\left(t) , and I ( t ) I\left(t) of Systems (2.18)–(2.20) are non-negative for all t ≥ 0 t\ge 0 .

Let x ( t ) = ( S ( t ) , E ( t ) , and I ( t ) ) x\left(t)=(S\left(t),E\left(t),\text{and}I\left(t)) be the solution of the system with the initial condition x 0 = ( S ( 0 ) , E ( 0 ) , and I ( 0 ) ) {x}_{0}=(S\left(0),E\left(0),\text{and}I\left(0)) . We assume that all the parameters of the model are positive. Because the quantity N h {N}_{h} measures the total amount of COVID-19 virus produced by an infected human, we consider N h {N}_{h} to be a positive parameter.

Now, using the continuity of solution, for all of S ( t ) , I ( t ) S\left(t),I\left(t) , and R ( t ) R\left(t) that have positive initial values at t = 0 t=0 , we have the existence of the solution on an interval ( 0 , t 0 ) \left(0,{t}_{0}) such that S ( t ) ≥ 0 S\left(t)\ge 0 , I ( t ) ≥ 0 I\left(t)\ge 0 , and R ( t ) ≥ 0 R\left(t)\ge 0 for 0 < t < t 0 0\lt t\lt {t}_{0} .

Now, let us say that we have t = t 1 ≥ 0 t={t}_{1}\ge 0 such that S ( t 1 ) = 0 S\left({t}_{1})=0 and E ( t 1 ) ≥ 0 E\left({t}_{1})\ge 0 , I ( t 1 ) ≥ 0 I\left({t}_{1})\ge 0 . Then, from System (2.18), we find that at t = t 1 t={t}_{1} d S d t t = t 1 = Λ ≥ 0 . {\left.\frac{{\rm{d}}S}{{\rm{d}}t}\right|}_{t={t}_{1}}=\Lambda \ge 0. This positive derivative ensures that at any time, the solution reaches the axis, its derivative increases, and the function S ( t ) S\left(t) does not cross to negative part. This shows that S ( t ) S\left(t) remains non-negative for all t ≥ 0 t\ge 0 .

Similarly, using the same analysis, we can show that d E d t t = t 1 = β N h S ( t 1 ) I ( t 1 ) ≥ 0 , {\left.\frac{{\rm{d}}E}{{\rm{d}}t}\right|}_{t={t}_{1}}=\beta {N}_{h}S\left({t}_{1})I\left({t}_{1})\ge 0, d I d t t = t 1 = π E ( t 1 ) ≥ 0 . {\left.\frac{{\rm{d}}I}{{\rm{d}}t}\right|}_{t={t}_{1}}=\pi E\left({t}_{1})\ge 0.

Thus, we find that x ( t ) x\left(t) never crosses the axes S = 0 , I = 0 S=0,I=0 , and R = 0 R=0 when it touches them. Thus, for any positive initial conditions, the solutions S ( t ) , E ( t ) S\left(t),E\left(t) , and I ( t ) I\left(t) of Systems (2.18)–(2.20) remain non-negative for all t ≥ 0 t\ge 0 .

Boundedness ̲ : \underline{{\bf{Boundedness}}}{\boldsymbol{:}}

Let N ( t ) = S ( t ) + E ( t ) + I ( t ) N\left(t)=S\left(t)+E\left(t)+I\left(t)

Now, d N d t = d S d t + d E d t + d I d t = Λ − μ ( S ( t ) + E ( t ) + I ( t ) ) − γ 1 E − γ 2 I − d N h ≤ Λ − μ N ( t ) . \begin{array}{rcl}\frac{{\rm{d}}N}{{\rm{d}}t}& =& \frac{{\rm{d}}S}{{\rm{d}}t}+\frac{{\rm{d}}E}{{\rm{d}}t}+\frac{{\rm{d}}I}{{\rm{d}}t}\\ \hspace{1em}& =& \Lambda -\mu \left(S\left(t)+E\left(t)+I\left(t))-{\gamma }_{1}E-{\gamma }_{2}I-d{N}_{h}\\ & \le & \Lambda -\mu N\left(t).\end{array} Therefore, d N d t + μ N ( t ) ≤ Λ . \frac{{\rm{d}}N}{{\rm{d}}t}+\mu N\left(t)\le \Lambda . The integrating factor is given by e μ t . {e}^{\mu t}. Therefore, after integration, we obtain

N ( t ) ≤ Λ μ + c e − μ t N\left(t)\le \frac{\Lambda }{\mu }+c{e}^{-\mu t} , where c c is a constant. Now, as t → ∞ t\to \infty , we obtain lim sup N(t) ≤ Λ μ . \hspace{0.1em}\text{lim sup N(t)}\hspace{0.1em}\le \frac{\Lambda }{\mu }. Thus, we have shown that the solutions of Systems (2.18)–(2.20) are ultimately bounded. Therefore, the biologically feasible region of Systems (2.18)–(2.20) is given by the following set: □ Ω = ( S ( t ) , E ( t ) , I ( t ) ) ∈ R + 3 : S ( t ) + E ( t ) + I ( t ) ≤ Λ μ , t ≥ 0 . \Omega =\left\{(S\left(t),E\left(t),I\left(t))\in {{\mathbb{R}}}_{+}^{3}:S\left(t)+E\left(t)+I\left(t)\le \frac{\Lambda }{\mu },\hspace{1em}t\ge 0\right\}.

The set Ω = ( S ( t ) , E ( t ) , I ( t ) ) ∈ R + 3 : S ( t ) + E ( t ) + I ( t ) ≤ Λ μ , t ≥ 0 \Omega =\left\{(S\left(t),E\left(t),I\left(t))\in {{\mathbb{R}}}_{+}^{3}:S\left(t)+E\left(t)+I\left(t)\le \frac{\Lambda }{\mu },\hspace{1em}t\ge 0\right\} is a positive invariant and an attracting set for Systems (2.18)–(2.20).

(Existence of solution) Let D be the domain defined above such that (2.22) holds. Then, there exists a unique solution of model system of equations (2.18)–(2.20), which is bounded in the domain D.

Let (2.23) f 1 = Λ − β N h S ( t ) I ( t ) − μ S ( t ) , {f}_{1}=\Lambda -\beta {N}_{h}S\left(t)I\left(t)-\mu S\left(t), (2.24) f 2 = β N h S ( t ) I ( t ) − ( μ + π + γ 1 ) E ( t ) , {f}_{2}=\beta {N}_{h}S\left(t)I\left(t)-\left(\mu +\pi +{\gamma }_{1})E\left(t), (2.25) f 3 = π E ( t ) − ( μ + γ 2 ) I ( t ) − d N h I ( t ) . {f}_{3}=\pi E\left(t)-\left(\mu +{\gamma }_{2})I\left(t)-d{N}_{h}I\left(t). We will show that ∂ f i ∂ x j , i , j = 1 , 2 , … , n , \frac{\partial {f}_{i}}{\partial {x}_{j}},\hspace{1em}i,j=1,2,\ldots ,n, is continuous and bounded in the domain D.

From equation (2.23), we have ∂ f 1 ∂ S = − β N h I − μ , ∂ f 1 ∂ S = ∣ − β N h I − μ ∣ < ∞ , ∂ f 1 ∂ E = 0 , ∂ f 1 ∂ E < ∞ , ∂ f 1 ∂ I = − β N h S , ∂ f 1 ∂ I = ∣ − β N h S ∣ < ∞ . \begin{array}{rcl}\frac{\partial {f}_{1}}{\partial S}& =& -\beta {N}_{h}I-\mu ,\hspace{1em}\left|\frac{\partial {f}_{1}}{\partial S}\right|=| -\beta {N}_{h}I-\mu | \lt \infty ,\\ \frac{\partial {f}_{1}}{\partial E}& =& 0,\hspace{1em}\left|\frac{\partial {f}_{1}}{\partial E}\right|\lt \infty ,\\ \frac{\partial {f}_{1}}{\partial I}& =& -\beta {N}_{h}S,\hspace{1em}\left|\frac{\partial {f}_{1}}{\partial I}\right|=| -\beta {N}_{h}S| \lt \infty .\end{array}

Similarly, from equation (2.24), we have ∂ f 2 ∂ S = β N h I , ∂ f 2 ∂ S = ∣ β N h I ∣ < ∞ , ∂ f 2 ∂ E = − μ − π − γ 1 , ∂ f 2 ∂ E = ∣ − ( μ + π + γ 1 ) ∣ < ∞ , ∂ f 2 ∂ I = β N h S , ∂ f 2 ∂ I = ∣ β N h S ∣ < ∞ . \begin{array}{rcl}\frac{\partial {f}_{2}}{\partial S}& =& \beta {N}_{h}I,\hspace{1em}\left|\frac{\partial {f}_{2}}{\partial S}\right|=| \beta {N}_{h}I| \lt \infty ,\\ \frac{\partial {f}_{2}}{\partial E}& =& -\mu -\pi -{\gamma }_{1},\hspace{1em}\left|\frac{\partial {f}_{2}}{\partial E}\right|=| -\left(\mu +\pi +{\gamma }_{1})| \lt \infty ,\\ \frac{\partial {f}_{2}}{\partial I}& =& \beta {N}_{h}S,\hspace{1em}\left|\frac{\partial {f}_{2}}{\partial I}\right|=| \beta {N}_{h}S| \lt \infty .\end{array}

Finally, from (2.25), we have ∂ f 3 ∂ S = 0 , ∂ f 3 ∂ S < ∞ , ∂ f 3 ∂ E = π , ∂ f 3 ∂ E = ∣ π ∣ < ∞ , ∂ f 3 ∂ I = − ( μ + γ 2 + d N h ) , ∂ f 3 ∂ I = ∣ − ( μ + γ 2 + d N h ) ∣ < ∞ . \begin{array}{rcl}\frac{\partial {f}_{3}}{\partial S}& =& 0,\hspace{1em}\left|\frac{\partial {f}_{3}}{\partial S}\right|\lt \infty ,\\ \frac{\partial {f}_{3}}{\partial E}& =& \pi ,\hspace{1em}\left|\frac{\partial {f}_{3}}{\partial E}\right|=| \pi | \lt \infty ,\\ \frac{\partial {f}_{3}}{\partial I}& =& -\left(\mu +{\gamma }_{2}+d{N}_{h}),\hspace{1em}\left|\frac{\partial {f}_{3}}{\partial I}\right|=| -\left(\mu +{\gamma }_{2}+d{N}_{h})| \lt \infty .\end{array}

Hence, we have shown that all the partial derivatives are continuous and bounded. Therefore, Lipschitz condition (2.22) is satisfied. Hence, by Theorem 2.1, there exists a unique solution of Systems (2.18)–(2.20) in the region D D .□

The infection-free equilibrium point E 0 {E}_{0} of Systems (2.18)–(2.20) is locally asymptotically stable provided R 0 < 1 {R}_{0}\lt 1 . If R 0 {R}_{0} crosses unity, E 0 {E}_{0} loses its stability and becomes unstable.

The infection-free equilibrium point E 0 {E}_{0} of Systems (2.18)–(2.20) is globally asymptotically stable whenever R 0 < 1 . {R}_{0}\lt 1.

We will prove the global stability of E 0 = ( ω μ , 0 , 0 ) {E}_{0}=\left(\frac{\omega }{\mu },0,0) of Systems (2.18)–(2.20) by showing that Systems (2.18)–(2.20) can be written in the general form stated in Theorem A.2 (Appendix) and both the conditions A 1 {A}_{1} and A 2 {A}_{2} of Theorem A.2 are satisfied.

Comparing the general system to Systems (2.18)–(2.20), the functions F F and G G are given by

F ( X , Y ) = Λ − β N h S I − μ S , F\left(X,Y)=\Lambda -\beta {N}_{h}SI-\mu S, G ( X , Y ) = ( β N h S I − ( m u + π + γ 1 ) E , π E − ( μ + γ 2 + d N h ) I ) , G\left(X,Y)=(\beta {N}_{h}SI-\left(mu+\pi +{\gamma }_{1})E,\hspace{0.33em}\pi E-\left(\mu +{\gamma }_{2}+d{N}_{h})I), where X = S X=S and Y = ( E , I ) Y=\left(E,I) .

The infection-free equilibrium point is U 0 = ( X 0 , 0 ¯ ) , {U}_{0}=\left({X}_{0},\bar{0}), where X 0 = Λ μ and 0 ¯ = ( 0 , 0 ) . {X}_{0}=\frac{\Lambda }{\mu }\hspace{1.0em}\hspace{0.1em}\text{and}\hspace{0.1em}\hspace{1.0em}\bar{0}=\left(0,0). From the stability analysis of E 0 {E}_{0} , we know that U 0 {U}_{0} is locally asymptotically stable iff R 0 < 1 {R}_{0}\lt 1 . Clearly, we see that G ( X , 0 ¯ ) = ( 0 , 0 ¯ ) G\left(X,\bar{0})=\left(0,\bar{0}) . Now, we show that X 0 = ( Λ μ ) {X}_{0}=\left(\frac{\Lambda }{\mu }) is globally asymptotically stable for the subsystem (2.29) d S d t = F ( S , 0 ¯ ) = ω − μ S . \frac{{\rm{d}}S}{{\rm{d}}t}=F\left(S,\bar{0})=\omega -\mu S.

The integrating factor is e μ t {e}^{\mu t} , and therefore, after performing integration on equation (2.29), we obtain S ( t ) e μ t = Λ e μ t μ + c . S\left(t){e}^{\mu t}=\frac{\Lambda {e}^{\mu t}}{\mu }+c. As t → ∞ t\to \infty , we obtain S ( t ) = ω μ , S\left(t)=\frac{\omega }{\mu }, which is independent of c. This independency implies that X 0 = Λ μ 1 {X}_{0}=\frac{\Lambda }{{\mu }_{1}} is globally asymptotically stable for the subsystem d S d t = Λ − μ S \frac{{\rm{d}}S}{{\rm{d}}t}=\Lambda -\mu S . So, the assumption A 1 {A}_{1} is satisfied.

Now, we will show that assumption A 2 {A}_{2} holds. First, we will find the matrix A A . As per the theorem, A = D Y G ( X , Y ) A={D}_{Y}G\left(X,Y) at X = X 0 X={X}_{0} and Y = 0 ¯ Y=\bar{0} . Now, D Y G ( X , Y ) = − ( μ + π + γ 1 ) β N h S π − ( μ + γ 2 + d N h ) . {D}_{Y}G\left(X,Y)=\left[\begin{array}{cc}-\left(\mu +\pi +{\gamma }_{1})& \beta {N}_{h}S\\ \pi & -\left(\mu +{\gamma }_{2}+d{N}_{h})\end{array}\right].

At X = X 0 X={X}_{0} and Y = 0 ¯ Y=\bar{0} , we obtain A = − ( μ + π + γ 1 ) β N h Λ μ π − ( μ + γ 2 + d N h ) . A=\left[\begin{array}{cc}-\left(\mu +\pi +{\gamma }_{1})& \frac{\beta {N}_{h}\Lambda }{\mu }\\ \pi & -\left(\mu +{\gamma }_{2}+d{N}_{h})\end{array}\right]. Clearly, matrix A A has non-negative off-diagonal elements. Hence, A A is a M-matrix. Using G ^ ( X , Y ) = A Y − G ( X , Y ) \widehat{G}\left(X,Y)=AY-G\left(X,Y) , we obtain G ^ ( X , Y ) = G 1 ^ ( X , Y ) G 2 ^ ( X , Y ) = β N h I Λ μ − S 0 \widehat{G}\left(X,Y)=\left[\begin{array}{c}\widehat{{G}_{1}}\left(X,Y)\\ \widehat{{G}_{2}}\left(X,Y)\end{array}\right]=\left[\begin{array}{c}\beta {N}_{h}I\left(\frac{\Lambda }{\mu }-S\right)\\ 0\end{array}\right] Hence, G 1 ^ ( X , Y ) = β N h I ( Λ μ − S ) ≥ 0 \widehat{{G}_{1}}\left(X,Y)=\beta {N}_{h}I\left(\frac{\Lambda }{\mu }-S)\ge 0 because S ( t ) ≤ Λ μ S\left(t)\le \frac{\Lambda }{\mu } and G 2 ^ ( X , Y ) = 0 \widehat{{G}_{2}}\left(X,Y)=0 .

Thus, both Assumptions A 1 {A}_{1} and A 2 {A}_{2} are satisfied, and therefore, infection-free equilibrium point E 0 {E}_{0} is globally asymptotically stable provided R 0 < 1 {R}_{0}\lt 1 .□

A unique infected equilibrium point E 1 {E}_{1} exists for Systems (2.18)–(2.20) and remains locally asymptotically stable provided that the following conditions are satisfied:

• (i)

  R 0 > 1 {R}_{0}\gt 1 .

• (ii)

  ( A 1 B 1 − C 1 ) > 0 \left({A}_{1}{B}_{1}-{C}_{1})\gt 0 .