A novel dynamic model describing the spread of virus

1. Introduction

Mathematical modeling of biological processes is a field of research that facilitates our understanding of complex biological systems and variables by simulating and analyzing several biological scenarios through specific models based on mathematical theories [1,2,3]. It provides a robust approach to describe mechanisms behind the dynamics occurring in biological processes by melting mathematical theories and models in the same pot with biological simulations [1]. This approach is playing a crucial role in predicting the spread rates of various major public health problems caused by viral diseases such as Ebola [4], Influenza [5], Cancer [6], Zika [7], Usutu [8], and Covid-19 [9].

Simultaneously, the field of mathematical modeling and simulations has gained considerable attention in various areas of human physiology, e.g., the analysis of muscle structure [10] and the study of brain activity [11,12]. Furthermore, mathematical modeling is also actively used to optimize medication usage and treatment process [13,14,15]. In this study, we focused mainly on modeling and simulating the spread of viruses and the behavior of cells in the human immune system through the dynamical modeling approach. Dynamic models are widely recognized for their pivotal role in describing the interactions among uninfected cells, free viruses, and immune responses [16,17,18,19]. The latest study findings highlighted the ability of complex models to describe human biology. For instance, Nowak et al. proposed a three-dimensional dynamic model for viral infection, which utilized numerical methods from autonomous dynamical systems [17,18,19]. Giesl and Wendland characterized a Lyapunov function as a solution of a suitable linear first-order partial differential equation approximating it using radial basis functions [20]. Yang and Wang formulated a mathematical model employing non-constant transmission rates, which varied with environmental conditions and the epidemiological status and reflected the impact of the ongoing disease control measures [21]. Despite significant efforts in designing mathematical models for virus dynamics, characterizing their behavior remains challenging. Even though several models have been proposed over the past two decades, their results bring different outcomes. However, a deterministic compartmental (SEIR) model has been successfully devised [22]. In mathematical epidemiology, numerous models have been proposed to analyze and predict epidemic outbreaks [23,24,25]. One such approach is the stereographic Brownian diffusion epidemiology model (SBDiEM), which provides a novel spatial perspective for the dynamic publishing, prediction, and modeling of infectious diseases [26,27]. The SBDiEM model can be adapted to identify past outbreaks and track the spread of viruses, offering valuable insights for national health systems, international stakeholders, and policymakers, thereby contributing to innovative perspectives in the field [28].

Kahajji et al. conducted a study focusing on the transmission dynamics of viruses, developing a separate mathematical model to describe the spread of the virus among animals in different regions. Their study highlights the importance of implementing effective campaigns to prevent individuals from moving between regions, promoting participation in quarantine centers, utilizing awareness campaigns targeted at virus prevention, and implementing security measures and health protocols within the region [29]. They estimated outbreak dynamics through mathematical modeling and provided decision guidelines for successful outbreak control. Furthermore, their model provided a valuable tool for estimating vaccination effectiveness and quantifying the impact of relaxing political measures, such as total lockdowns, shelter-in-place orders, and travel restrictions, for both low-risk subgroups and the population as a whole [30]. Numerous studies have highlighted the importance of models and simulations in understanding complex biological phenomena and tackling emerging diseases. However, the diversity of biological factors affecting human health and the emergence of new diseases warrants further exploration.

Moreover, the mathematical approaches identified by the World Health Organization (WHO) can be essential in providing evidence-based information to healthcare decision-makers and policymakers [31].

In this study, we have theorized a nonlinear mathematical model of virus transmission based on the SEIR model. We consider the following mathematical model concerning the initial value problem for the following nonlinear systems:

$$\left\{\begin{array}{cc}\hfill \dot{T}\left(t\right)& =a-{\beta }_{1}V\left(t\right)T\left(t\right)-d_{1}T\left(t\right),\hfill \\ \hfill \dot{I}\left(t\right)& =qT\left(t\right)V\left(t\right)-{\beta }_{2}E\left(t\right)I\left(t\right)-d_{2}I\left(t\right),\hfill \\ \hfill \dot{E}\left(t\right)& ={\beta }_{3}I\left(t\right)E\left(t\right)-d_{3}E\left(t\right),\hfill \\ \hfill \dot{V}\left(t\right)& =bI\left(t\right)-cV\left(t\right)\hfill \end{array}\right.$$

(1)

$$T\left(t_0\right)=T_0,I\left(t_0\right)=I_0,E\left(t_0\right)=E_0,$$

(2)

$$V\left(t_0\right)=V_0,t_0\in \left[0,\right.\left.a\right),$$

where $T=T\left(t\right)$, $I=I\left(t\right)$, $E\left(t\right)$ and $V\left(t\right)$ denotes the concentration of uninfected cells, infected cells, effector immune cells, and free viruses at time $t\in (0,m)$, respectively.

Uninfected cells are supplied at a rate a, and uninfected hepatocytes (target cells, T) are infected by virus V at ${\beta }_1$. They are $d_i$ ($i=1,2,3$) that die naturally at the rate q is the rate constant characterizing infection of the infected cells. Effector cells mediate infection by eliminating productively infected cells at a rate of ${\beta }_2$. Effector immune cells E are supplied to the presence of tumor cells, stimulating the immune response. The virus activates the effector immune cells at the rate of ${\beta }_3$. The infected cells produce new viruses at the rate of b during their life. The constant $c>0$ is the rate at which the viruses are cleared [32,33,34]. They worked on the local and global dynamic model of cancer tumor growth [35]. Recently, studies have been carried out on dynamic models [6,35,36,37].

2. Boundedness and Dissipativity

In this section, we showed that the model is bounded by negative divergence, positively invariant with respect to a region in ${\mathbb{R}}_{+}^4$ and dissipative. As we are interested in biologically relevant solutions of the system, the next results show that the positive octant is invariant and that the upper limits of trajectories depend on the parameters.

We put

$$T\left(t\right)=x_1\left(t\right),I\left(t\right)=x_2\left(t\right),E\left(t\right)=x_3\left(t\right),V\left(t\right)=x_4\left(t\right).$$

Then the problem $\left(1\right)-\left(2\right)$ is reduced the following form:

$${\dot{x}}_1\left(t\right)=a-{\beta }_1{x}_4\left(t\right)x_1\left(t\right)-d_1{x}_1\left(t\right),$$

(3)

$${\dot{x}}_2\left(t\right)=q{x}_1\left(t\right)x_4\left(t\right)-{\beta }_2{x}_3\left(t\right)x_2\left(t\right)-d_2{x}_2\left(t\right),$$

$${\dot{x}}_3\left(t\right)={\beta }_3{x}_2\left(t\right)x_3\left(t\right)-d_3{x}_3\left(t\right),$$

$${\dot{x}}_4\left(t\right)=b{x}_2\left(t\right)-c{x}_4\left(t\right),$$

$$x_1\left(t_0\right)=x_{10},x_2\left(t_0\right)=x_{20},x_3\left(t_0\right)=x_{30},$$

(4)

$$x_4\left(t_0\right)=x_{40},t_0\in \left[0,\right.\left.a\right).$$

Let

$$x=x\left(t\right)=\left(x_1,x_2,x_3,x_4\right),x_j=x_j\left(t\right),j=1,2,3,4,$$

(5)

$$f_1\left(x\right)=a-{\beta }_1{x}_4\left(t\right)x_1\left(t\right)-d_1{x}_1\left(t\right)$$

$$f_2\left(x\right)=q{x}_1\left(t\right)x_4\left(t\right)-{\beta }_2{x}_3\left(t\right)x_2\left(t\right)-d_2{x}_2\left(t\right)$$

$$f_3\left(x\right)={\beta }_3{x}_2\left(t\right)x_3\left(t\right)-d_3{x}_3\left(t\right),f_4\left(x\right)=b{x}_2\left(t\right)-c{x}_4\left(t\right).$$

Here,

$${\mathbb{R}}_{+}^4=\left\{x=\left(x_1,x_2,x_3,x_4\right)\in {\mathbb{R}}^4,x_k>0\right\},$$

$$\Omega =\left\{x\in {\mathbb{R}}_{+}^4\text{:}\right.{\beta }_3{x}_2-{\beta }_1{x}_4-{\beta }_2{x}_3\le d_1+d_2+d_3+c\}$$

Consider the problem $\left(3\right)-\left(4\right)$ with $t_0=0.$

Condition 1. 

Assume the following assumption is satisfied

$${\beta }_1,{\beta }_2,d_1,d_2,d_3,c>0,{\beta }_3<0.$$

Theorem 1. 

Let the Condition 1 holds. Then the system $\left(3\right)$ is with the negative divergence and is dissipative in the domain $\Omega \subset {\mathbb{R}}_{+}^4$.

Proof. 

Indeed, from $\left(3\right)$ and $\left(4\right)$ we have

$$\frac{\partial f_1}{\partial x_1}+\frac{\partial f_2}{\partial x_2}+\frac{\partial f_3}{\partial x_3}+\frac{\partial f_4}{\partial x_4}=-\left({\beta }_1{x}_4+d_1\right)-$$

$${\beta }_2{x}_3\left(t\right)-d_2+{\beta }_3{x}_2\left(t\right)-d_3-c.$$

Hence, by Condition 1, the system $\left(3\right)$ is dissipative on the domain $\Omega .$ but there is no definition of Condition 1. □

3. The Local Stability of Equilibria Points

In this section, we derive the stability properties of equilibria points of the system $\left(1\right)$. Let

$${\mathbb{R}}_{+}^4=\left\{x\in {\mathbb{R}}^4\text{:}{x}_i\ge 0,i=1,2,3,4\right\},B_r\left(\overline{x}\right)=\left\{x\in {\mathbb{R}}^4,{∥x-\overline{x}∥}_{{\mathbb{R}}^3}<r\right\}.$$

Condition 2. 

Let

$$\frac{b{d}_3}{c}\ne d_1,b{d}_3\ne d_{1}c,\frac{b{a}_{33}{a}_{24}+a_{32}{a}_{23}c+a_{33}{a}_{21}{a}_{14}b}{a_{21}{a}_{14}b+b{a}_{33}{a}_{24}-a_{32}{a}_{23}}<0,d_3>{\beta }_3{\overline{x}}_2.$$

(6)

Theorem 2. 

Assume that the Condition 2 is satisfied. There is a point $P=\left(x_1,x_2,x_3,x\right)$ that is a equilibria points of the system $\left(1\right)$ in ${\mathbb{R}}_{+}^4$.

Proof. 

It is sufficient to find the solution of the following system of algebraic equation in $x_1$, $x_2$, $x_3$, $x_4$:

$$a-\left({\beta }_1{x}_4-d_1\right)x_1=0,q{x}_1{x}_4-{\beta }_2{x}_3{x}_2-d_2{x}_2=0,$$

$${\beta }_3{x}_2{x}_3-d_3{x}_3=0,b{x}_2-c{x}_4=0.$$

(7)

From first and second equations we have

$${\overline{x}}_1=\frac{a}{\left({\beta }_1{x}_4-d_1\right)},q{x}_1{x}_4-{\beta }_2{x}_3{x}_2-d_2{x}_2=0.$$

(8)

From third and fourth equations we get

$$\left({\beta }_3{x}_2-d_3\right)x_3=0,x_4=\frac{b}{c}{x}_2.$$

(9)

If $x_3\ne 0$ we get that ${\overline{x}}_2=\frac{d_3}{{\beta }_3}$. By $\left(9\right)$, then we dedused that ${\overline{x}}_4=\frac{b{d}_3}{c{\beta }_3}$. Hence, from $\left(8\right)$ we have

$${\overline{x}}_1=\frac{a}{\left(\frac{b{d}_3}{c}-d_1\right)},$$

$${\overline{x}}_3=\frac{1}{{\beta }_2{\overline{x}}_2}\left[q{\overline{x}}_1{\overline{x}}_4-d_2{\overline{x}}_2\right]=\frac{abq}{\left(b{d}_3-d_{1}c\right){\beta }_2}-\frac{d_2}{{\beta }_2}.$$

Thus we obtain that the system $\left(1\right)$ have a unique equalibria point $P\left({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3,{\overline{x}}_4\right)$, where

$${\overline{x}}_1=\frac{a}{\left(\frac{b{d}_3}{c}-d_1\right)},{\overline{x}}_2=\frac{d_3}{{\beta }_3},{\overline{x}}_3=\frac{abq}{\left(b{d}_3-d_{1}c\right){\beta }_2}-\frac{d_2}{{\beta }_2},{\overline{x}}_4=\frac{b{d}_3}{c{\beta }_3}.$$

(10)

□

Remark 1. 

For the point to have the biological meaning of stability point $P\left({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3,{\overline{x}}_4\right)$, it should be:

$$d_1<\frac{b{d}_3}{c},b{d}_3-d_{1}c>0,{\beta }_2{\overline{x}}_3+d_2>0.$$

(11)

We show here, the following results:

Theorem 3. 

Assume that the Condition 2 is satisfied. Suppose the estimate $\left(11\right)$ holds. Then the point $P\left({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3,{\overline{x}}_4\right)$ is locally stable point for the system of $\left(1\right).$

Proof. 

Consider the linearized matrix of $\left(1\right)$, i.e. the Jacobian matrıx according to system $\left(1\right)$ at point $P\left({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3,{\overline{x}}_4\right)$ is the following:

$$A=\frac{Df}{Dx}=\left[\begin{array}{cccc}\frac{\partial f_1}{\partial x_1}& \frac{\partial f_1}{\partial x_2}& \frac{\partial f_1}{\partial x_3}& \frac{\partial f_1}{\partial x_4}\\ \frac{\partial f_2}{\partial x_1}& \frac{\partial f_2}{\partial x_2}& \frac{\partial f_2}{\partial x_3}& \frac{\partial f_2}{\partial x_4}\\ \frac{\partial f_3}{\partial x_1}& \frac{\partial f_3}{\partial x_2}& \frac{\partial f_3}{\partial x_3}& \frac{\partial f_3}{\partial x_4}\\ \frac{\partial f_4}{\partial x_1}& \frac{\partial f_4}{\partial x_2}& \frac{\partial f_4}{\partial x_2}& \frac{\partial f_4}{\partial x_3}\end{array}\right]=\left[\begin{array}{cccc}{a}_{11}& 0& 0& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ 0& a_{32}& a_{33}& 0\\ 0& b& 0& -c\end{array}\right],$$

(12)

where ${\beta }_3{x}_2\left(t\right)x_3\left(t\right)-d_3{x}_3\left(t\right)$

$$a_{11}=-\left({\beta }_1{\overline{x}}_4+d_1\right),a_{14}=-{\beta }_1{\overline{x}}_1,a_{21}=q{\overline{x}}_4,a_{22}=-\left({\beta }_2{\overline{x}}_3+d_2\right),$$

$$a_{23}=-{\beta }_2{\overline{x}}_2,a_{24}=q{\overline{x}}_1,a_{32}={\beta }_3{\overline{x}}_3,a_{32}={\beta }_3{\overline{x}}_3,a_{33}={\beta }_3{\overline{x}}_2-d_3.$$

(13)

□

The eigenvalues of the matrix A can found as the solutions of the following equations

$$A-\lambda I=\left[\begin{array}{cccc}{a}_{11}-\lambda & 0& 0& a_{14}\\ a_{21}& a_{22}-\lambda & a_{23}& a_{24}\\ 0& a_{32}& a_{33}-\lambda & 0\\ 0& b& 0& -\left(c+\lambda \right)\end{array}\right]=$$

$$\left[a_{11}-\lambda \right]\left[\begin{array}{ccc}{a}_{22}-\lambda & a_{23}& a_{24}\\ a_{32}& a_{33}-\lambda & 0\\ b& 0& -\left(c+\lambda \right)\end{array}\right]-$$

$$a_{21}\left(i\right)\left[\begin{array}{ccc}0& 0& a_{14}\\ a_{32}& a_{33}-\lambda & 0\\ b& 0& -\left(c+\lambda \right)\end{array}\right]=$$

(14)

$$\left(a_{11}-\lambda \right)\left[\right.-\left(c+\lambda \right)\left(a_{22}-\lambda \right)\left(a_{33}-\lambda \right)-$$

$$b{a}_{24}\left(a_{33}-\lambda \right)+\left.a_{32}{a}_{23}\left(c+\lambda \right)\right]+a_{21}{a}_{14}b\left(a_{33}-\lambda \right)=0.$$

Let we put

$$\left(c+\lambda \right)\left(a_{22}-\lambda \right)\left(a_{33}-\lambda \right)=0.$$

i.e. ${\lambda }_1=-c$, ${\lambda }_2=a_{22}$ and ${\lambda }_3=a_{33}$ are the eigenvalues of A. Then other solutions of A can be obtained by solving the equation

$$b{a}_{24}\left(a_{33}-\lambda \right)+\left.a_{32}{a}_{23}\left(c+\lambda \right)\right]+a_{21}{a}_{14}b\left(a_{33}-\lambda \right)=$$

(15)

$$\left(a_{21}{a}_{14}b+b{a}_{33}{a}_{24}-a_{32}{a}_{23}\right)\lambda =b{a}_{33}{a}_{24}+a_{32}{a}_{23}c+a_{33}{a}_{21}{a}_{14}b=0$$

By solving the equation $\left(15\right)$ we get the fourth eigenvalue of the matrix A

$${\lambda }_4=\frac{b{a}_{33}{a}_{24}+a_{32}{a}_{23}c+a_{33}{a}_{21}{a}_{14}b}{a_{21}{a}_{14}b+b{a}_{33}{a}_{24}-a_{32}{a}_{23}}.$$

For local stability of the system $\left(1\right)$ it is sufficient to show that all eigenvalues of the matrix A are negative. Indeed, by $\left(10\right)$ and $\left(11\right)$, we have

$${\lambda }_1=-c<0,{\lambda }_2=a_{22}=-\left({\beta }_2{\overline{x}}_3+d_2\right)<0,$$

(16)

$${\lambda }_3=a_{33}={\beta }_3{\overline{x}}_2-d_3<0,{\lambda }_4=\frac{b{a}_{33}{a}_{24}+a_{32}{a}_{23}c+a_{33}{a}_{21}{a}_{14}b}{a_{21}{a}_{14}b+b{a}_{33}{a}_{24}-a_{32}{a}_{23}}<0.$$

By assumption $\left(11\right)$, and by $\left(10\right)$ we see that,

$${\overline{x}}_1\ge 0,{\overline{x}}_2\ge 0,{\overline{x}}_3\ge 0,{\overline{x}}_4\ge 0.$$

Hence by $\left(16\right)$ we get

$${\lambda }_1<0,{\lambda }_2<0,{\lambda }_3<0.$$

Moreover, it should be

$${\lambda }_4=\frac{b{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2-c{\beta }_2{\beta }_3{\overline{x}}_2{\overline{x}}_3-b{\beta }_1{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2{x}_4}{-b{\beta }_{1}q{\overline{x}}_1{\overline{x}}_4+b{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2+{\beta }_2{\beta }_3{\overline{x}}_2{\overline{x}}_3}<0.$$

(17)

The estimate $\left(17\right)$ satisfies if:

$$b{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2-c{\beta }_2{\beta }_3{\overline{x}}_2{\overline{x}}_3-b{\beta }_1{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2{\overline{x}}_4<0,$$

(18)

$$-b{\beta }_{1}q{\overline{x}}_1{\overline{x}}_4+b{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2+{\beta }_2{\beta }_3{\overline{x}}_2{\overline{x}}_3>0..$$

or

$$b{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2-c{\beta }_2{\beta }_3{\overline{x}}_2{\overline{x}}_3-b{\beta }_1{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2{\overline{x}}_4>0,$$

(19)

$$-b{\beta }_{1}q{\overline{x}}_1{\overline{x}}_4+b{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2+{\beta }_2{\beta }_3{\overline{x}}_2{x}_3<0;$$

Since $x_k\ge 0$ the second inequality in $\left(18\right)$ satisfied for all $x\in {\mathbb{R}}_{+}^4$, when

$$b{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2-c{\beta }_2{\beta }_3{\overline{x}}_2{\overline{x}}_3-b{\beta }_1{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2{\overline{x}}_4<0,$$

i.e. by $\left(10\right)$ if

$$\frac{bqa{d}_{3}c}{b{d}_3-d_{1}c}<\frac{c{d}_{3}abq+b{\beta }_{1}qa{d}_3^{2}b}{\left(b{d}_3-d_{1}c\right)}.$$

By assumption $\left(11\right)$, the above inequality is satisfied when

$$d_{3}c<c{d}_3+{\beta }_1{d}_3^{2}b,$$

that it is clear that holds for all $x\in {\mathbb{R}}_{+}^4$. Since $b{\beta }_{1}q{\overline{x}}_1{\overline{x}}_4+b{\beta }_{3}q{\overline{x}}_1{\overline{x}}_2+{\beta }_2{\beta }_3{\overline{x}}_2{\overline{x}}_3\ge 0$ for all $x\in {\mathbb{R}}_{+}^4$, the inequality $\left(19\right)$is not satisfied in ${\mathbb{R}}_{+}^4$ . Hence we obtained that all eigenvalues of the matrix are negative under our assumptions.

In Figure 1 we compare the number of viruses with the number of infected cells. Both the number of viruses and the number of infected cells decrease over time. The second Figure 2 compares the number of viruses with the number of uninfected cells. In this case, the rates of change between the two are in inverse proportion to each other. Finally, Figure 3 shows a comparison between the number of infected cells and the number of effector immune cells. It is noticeable that both the infected cells and the immune cells in this figure rapidly decrease over time.

4. Lyapunov Stability of Equilibria Points

Let $E\left(\overline{x}\right)$ is a equilibria point, where $\overline{x}=\left({\overline{x}}_1,{\overline{x}}_2,{\overline{x}}_3,{\overline{x}}_4\right)\in {\mathbb{R}}_{+}^4$ is defined by $\left(10\right)$. In this section, we show the following results: Let $A=A\left(\overline{x}\right)$ be the linearized matrix with respect to equilibria $E\left(\overline{x}\right)$ point defined by $\left(12\right)$, i.e.

$$A=\left[\begin{array}{cccc}{a}_{11}& 0& 0& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ 0& a_{32}& a_{33}& 0\\ 0& b& 0& -c\end{array}\right],$$

where $a_{ij}$ are defined by $\left(13\right)$. We consider the Lyapunov equation

$$BA+A^{T}B=-I,B=B\left(\overline{x}\right)=\left[\begin{array}{cccc}{b}_{11}& b_{12}& b_{13}& b_{14}\\ b_{21}& b_{22}& b_{23}& b_{24}\\ b_{31}& b_{32}& b_{33}& b_{34}\\ b_{41}& b_{42}& b_{43}& b_{44}\end{array}\right],b_{ij}=b_{ji}.$$

(20)

It is clear that

$$BA=\left[\begin{array}{cccc}{b}_{11}& b_{12}& b_{13}& b_{14}\\ b_{21}& b_{22}& b_{23}& b_{24}\\ b_{31}& b_{32}& b_{33}& b_{34}\\ b_{41}& b_{42}& b_{43}& b_{44}\end{array}\right]\left[\begin{array}{cccc}{a}_{11}& 0& 0& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ 0& a_{32}& a_{33}& 0\\ 0& b& 0& -c\end{array}\right]=$$

$$\left[\begin{array}{cccc}-a_{11}{b}_{11}+a_{21}{b}_{12}& a_{22}{b}_{12}+a_{32}{b}_{13}+b{b}_{14}& a_{23}{b}_{12}+a_{33}{b}_{13}& a_{14}{b}_{11}+a_{24}{b}_{12}-c{b}_{14}\\ a_{11}{b}_{21}+a_{21}{b}_{22}-& a_{22}{b}_{22}+a_{32}{b}_{23}+b{b}_{24}& a_{23}{b}_{22}+a_{33}{b}_{23}& a_{14}{b}_{21}+a_{24}{b}_{22}-c{b}_{24}\\ a_{11}{b}_{31}+a_{21}{b}_{32}-& a_{22}{b}_{32}+a_{32}{b}_{33}+b{b}_{34}& a_{23}{b}_{32}+a_{33}{b}_{33}& a_{14}{b}_{31}+a_{24}{b}_{32}-c{b}_{34}\\ a_{11}{b}_{41}+a_{21}{b}_{42}-& a_{22}{b}_{42}+a_{32}{b}_{43}+b{b}_{44}& a_{23}{b}_{42}+a_{33}{b}_{43}& a_{14}{b}_{41}+a_{24}{b}_{42}-c{b}_{44}\end{array}\right],$$

$$A^{T}B=\left[\begin{array}{cccc}-a_{11}& a_{21}& 0& 0\\ 0& a_{22}& a_{32}& b\\ 0& a_{23}& -a_{33}& 0\\ a_{14}& a_{24}& 0& -c\end{array}\right]\left[\begin{array}{cccc}{b}_{11}& b_{12}& b_{13}& b_{14}\\ b_{21}& b_{22}& b_{23}& b_{24}\\ b_{31}& b_{32}& b_{33}& b_{34}\\ b_{41}& b_{42}& b_{43}& b_{44}\end{array}\right]=$$

$$\left[\begin{array}{cccc}{a}_{11}{b}_{11}+a_{21}{b}_{21}& a_{11}{b}_{12}+a_{21}{b}_{22}& a_{11}{b}_{13}+a_{21}{b}_{23}& a_{11}{b}_{14}+a_{21}{b}_{24}\\ a_{22}{b}_{21}+a_{32}{b}_{31}+b{b}_{41}& a_{22}{b}_{22}+a_{32}{b}_{32}+b{b}_{42}& a_{22}{b}_{23}+a_{32}{b}_{33}+b{b}_{43}& a_{22}{b}_{24}+a_{32}{b}_{34}+b{b}_{44}\\ a_{23}{b}_{21}+a_{33}{b}_{31}& a_{23}{b}_{22}+a_{33}{b}_{32}& a_{23}{b}_{23}+a_{33}{b}_{33}& a_{23}{b}_{24}+a_{33}{b}_{34}\\ -a_{14}{b}_{11}+a_{24}{b}_{21}-c{b}_{41}& a_{14}{b}_{12}+a_{24}{b}_{22}-c{b}_{42}& a_{14}{b}_{13}+a_{24}{b}_{23}-c{b}_{43}& a_{14}{b}_{14}+a_{24}{b}_{24}-c{b}_{44}\end{array}\right].$$

$\left(20\right)$ reduced to the following equation

$$BA+A^{T}B=\left[\begin{array}{cccc}{g}_{11}& g_{12}& g_{13}& g_{14}\\ g_{21}& g_{22}& g_{23}& g_{24}\\ g_{31}& g_{32}& g_{33}& g_{34}\\ g_{41}& g_{42}& g_{43}& g_{44}\end{array}\right]=-I,$$

(21)

where

$$g_{11}=2a_{11}{b}_{11}+2a_{21}{b}_{12}=-1,$$

$$g_{12}=\left(a_{11}+a_{22}\right)b_{12}+a_{32}{b}_{13}+b{b}_{14}+a_{21}{b}_{22}=0,$$

$$g_{13}=a_{23}{b}_{12}+\left(a_{11}+a_{33}\right)b_{13}+a_{21}{b}_{23}=0,$$

$$g_{14}=a_{14}{b}_{11}+a_{24}{b}_{12}+\left(a_{11}-c\right)b_{14}+a_{21}{b}_{24}=0,$$

$$g_{22}=2a_{22}{b}_{22}+2a_{32}{b}_{23}+2b{b}_{24}=-1,$$

$$g_{23}=a_{23}{b}_{22}+\left(a_{22}+a_{33}\right)b_{23}+a_{32}{b}_{33}+b{b}_{34}=0,$$

(22)

$$g_{24}=a_{14}{b}_{12}+a_{24}{b}_{22}+\left(a_{22}-c\right)b_{24}+a_{23}{b}_{34}+b{b}_{44}=0,$$

$$g_{33}=2a_{23}{b}_{23}+2a_{33}{b}_{33}=-1,$$

$$g_{34}=a_{14}{b}_{13}+a_{24}{b}_{23}+\left(a_{33}-c\right)b_{34}+a_{23}{b}_{24}=0,$$

$$g_{44}=2a_{14}{b}_{14}+2a_{24}{b}_{24}-2c{b}_{44}=-1.$$

Main and associated determinants of the system $\left(22\right)$ in $b_{11},b_{12},$$b_{13},b_{14},b_{22},b_{23},b_{24},b_{33},b_{34},b_{44}$ is the following

$$\Delta =\left|\begin{array}{cccccccccc}2a_{11}& 2a_{21}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& a_{11}+a_{22}& a_{11}+a_{33}& a_{21}& 0& 0& 0& 0& 0& 0\\ 0& a_{23}& a_{11}+a_{33}& 0& 0& a_{21}& 0& 0& 0& 0\\ a_{14}& a_{24}& 0& \left(a_{11}-c\right)& 0& 0& a_{21}& 0& 0& 0\\ 0& 0& 0& 0& 2a_{22}& 2a_{32}& 2b& 0& 0& 0\\ 0& 0& 0& 0& a_{23}& a_{22}+a_{33}& 0& a_{32}& b& 0\\ 0& a_{14}& 0& 0& a_{24}& 0& \left(a_{22}-c\right)& 0& a_{23}& b\\ 0& 0& 0& 0& 0& 2a_{23}& 0& 2a_{33}& 0& 0\\ 0& 0& a_{14}& 0& 0& a_{24}& a_{23}& 0& \left(a_{33}-c\right)& 0\\ 0& 0& 0& 2a_{14}& 0& 0& 2a_{24}& 0& 0& -2cb\end{array}\right|,.$$

$${\Delta }_1=\left|\begin{array}{cccccccccc}-1& 2a_{21}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& a_{11}+a_{22}& a_{11}+a_{33}& a_{21}& 0& 0& 0& 0& 0& 0\\ 0& a_{23}& a_{11}+a_{33}& 0& 0& a_{21}& 0& 0& 0& 0\\ -0& a_{24}& 0& \left(a_{11}-c\right)& 0& 0& a_{21}& 0& 0& 0\\ -1& 0& 0& 0& 2a_{22}& 2a_{32}& 2b& 0& 0& 0\\ 0& 0& 0& 0& a_{23}& a_{22}+a_{33}& 0& a_{32}& b& 0\\ 0& a_{14}& 0& 0& a_{24}& 0& \left(a_{22}-c\right)& 0& a_{23}& b\\ -1& 0& 0& 0& 0& 2a_{23}& 0& 2a_{33}& 0& 0\\ 0& 0& a_{14}& 0& 0& a_{24}& a_{23}& 0& \left(a_{33}-c\right)& 0\\ -1& 0& 0& 2a_{14}& 0& 0& 2a_{24}& 0& 0& -2cb\end{array}\right|,$$

$${\Delta }_2=\left|\begin{array}{cccccccccc}2a_{11}& -1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& a_{11}+a_{33}& a_{21}& 0& 0& 0& 0& 0& 0\\ 0& 0& a_{11}+a_{33}& 0& 0& a_{21}& 0& 0& 0& 0\\ a_{14}& 0& 0& \left(a_{11}-c\right)& 0& 0& a_{21}& 0& 0& 0\\ 0& -1& 0& 0& 2a_{22}& 2a_{32}& 2b& 0& 0& 0\\ 0& 0& 0& 0& a_{23}& a_{22}+a_{33}& 0& a_{32}& b& 0\\ 0& 0& 0& 0& a_{24}& 0& \left(a_{22}-c\right)& 0& a_{23}& b\\ 0& -1& 0& 0& 0& 2a_{23}& 0& 2a_{33}& 0& 0\\ 0& 0& a_{14}& 0& 0& a_{24}& a_{23}& 0& \left(a_{33}-c\right)& 0\\ 0& -1& 0& 2a_{14}& 0& 0& 2a_{24}& 0& 0& -2cb\end{array}\right|,$$

...

$${\Delta }_{10}=\left|\begin{array}{cccccccccc}2a_{11}& 2a_{21}& 0& 0& 0& 0& 0& 0& 0& -1\\ 0& a_{11}+a_{22}& a_{11}+a_{33}& a_{21}& 0& 0& 0& 0& 0& 0\\ 0& a_{23}& a_{11}+a_{33}& 0& 0& a_{21}& 0& 0& 0& 0\\ a_{14}& a_{24}& 0& \left(a_{11}-c\right)& 0& 0& a_{21}& 0& 0& 0\\ 0& 0& 0& 0& 2a_{22}& 2a_{32}& 2b& 0& 0& -1\\ 0& 0& 0& 0& a_{23}& a_{22}+a_{33}& 0& a_{32}& b& 0\\ 0& a_{14}& 0& 0& a_{24}& 0& \left(a_{22}-c\right)& 0& a_{23}& 0\\ 0& 0& 0& 0& 0& 2a_{23}& 0& 2a_{33}& 0& -1\\ 0& 0& a_{14}& 0& 0& a_{24}& a_{23}& 0& \left(a_{33}-c\right)& 0\\ 0& 0& 0& 2a_{14}& 0& 0& 2a_{24}& 0& 0& -1\end{array}\right|..$$

We assume that $\Delta \ne 0$. Then by by solving $\left(22\right)$ with respect to $b_{ıj}$ by Kramer method, we obtain

$$b_{11}=\frac{{\Delta }_1}{\Delta },b_{12}=b_{21}=\frac{{\Delta }_2}{\Delta },...,b_{44}=\frac{{\Delta }_n}{\Delta }.$$

(23)

Theorem 4. 

Assume the Condition 2 holds, $\Delta \ne 0$. Suppose $a_{ij}$ such that $b_{ii}>0$, $i=1,2,3,4$, $b_{ij}\ge 0$ for $i,j=1,2,3,4$, when $i\ne j$. Then the system $\left(3\right)$ is asymptotically stable at the equilibria point $E\left(\overline{x}\right)$ in the sense of Lyapunov.

Proof. 

By assumptions the function $P_A\left(x\right)$ associated with tha matrix A defined by

$$P_B\left(x\right)=x^{T}Bx=\sum _{i,j=1}^4{b}_{ij}{x}_i{x}_j$$

is positive defined in ${\mathbb{R}}^4.$ Hence, all eigenvalues ${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4$ of the the matrix $B=B\left(\overline{x}\right)$ is positive in ${\mathbb{R}}^4$, i.e. $P_B\left(x\right)$ is a positive defined Lyapunov function candidate (see e.g. $\left[22,23\right]$). By $\left[12,\mathrm{Corollary}8.2\right]$. We need now to determine a domain $\Omega$ on which ${\dot{P}}_B\left(x\right)$ is negatively defined. By assuming $x_k\ge 0$, $k=1,2,3,4$ we will find the solution set of the following inequality

$$f_1\left(x\right)=a-{\beta }_1{x}_4\left(t\right)x_1\left(t\right)-d_1{x}_1\left(t\right)$$

$$f_2\left(x\right)=q{x}_1\left(t\right)x_4\left(t\right)-{\beta }_2{x}_3\left(t\right)x_2\left(t\right)-d_2{x}_2\left(t\right)$$

$$f_3\left(x\right)={\beta }_3{x}_2\left(t\right)x_3\left(t\right)-d_3{x}_3\left(t\right),f_4\left(x\right)=b{x}_2\left(t\right)-c{x}_4\left(t\right).$$

$${\dot{P}}_B\left(x\right)=\sum _{j=1}^4\frac{\partial {\dot{P}}_B}{\partial x_j}{f}_j\left(x\right)=$$

$$2\left(b_{11}{x}_1+b_{12}{x}_2+b_{13}{x}_3+b_{14}{x}_4\right)\left(a-{\beta }_1{x}_4{x}_1-d_1{x}_1\right)+$$

$$2\left(b_{21}{x}_1+b_{22}{x}_2+b_{23}{x}_3+b_{24}{x}_4\right)\left(q{x}_1{x}_4-{\beta }_2{x}_3{x}_2-d_2{x}_2\right)+$$

$$2\left(b_{31}{x}_1+b_{32}{x}_2+b_{33}{x}_3+b_{34}{x}_4\right)\left({\beta }_3{x}_2{x}_3-d_3{x}_3\right)+$$

(24)

$$2\left(b_{31}{x}_1+b_{32}{x}_2+b_{33}{x}_3+b_{34}{x}_4\right)\left(b{x}_2-c{x}_4\right)\le 0.$$

Hence, the system $\left(3\right)$ is asymptotically stabile at $E\left(\overline{x}\right)$ in the Lyapunov sense when,

$$a-{\beta }_1{x}_4{x}_1-d_1{x}_1\le 0,q{x}_1{x}_4-{\beta }_2{x}_3{x}_2-d_2{x}_2\le 0,$$

$${\beta }_3{x}_2{x}_3-d_3{x}_3\le 0,b{x}_2-c{x}_4,$$

i.e. the system $\left(3\right)$ is asymptotically stabile at $E\left(\overline{x}\right)$ in the Lyapunov sense in the following domain

$${\Omega }_1=\left\{x\in {\mathbb{R}}_{+}^4\text{:}\right.\left({\beta }_1{x}_4+d_1\right)x_1\ge a,\left({\beta }_2{x}_3+d_2\right)x_2\ge q{x}_1{x}_4,$$

$$\left.{\beta }_3{x}_2\le d_3\right\},x_4\ge \frac{b}{c}{x}_2.$$

(25)

□

Theorem 5. 

Assume the Condition 2 holds, $\Delta \ne 0$. Suppose $a_{ij}$ such that $b_{ii}>0$, $i=1,2,3,4$ and $b_{ij}\le 0$ for $i,j=1,2,3,4$ when $i\ne j$. Then the system $\left(3\right)$ is asymptotically stable at the equilibria point $E\left(\overline{x}\right)$ in the sense of Lyapunov.

Proof. 

$$P_B\left(x\right)=x^{T}Bx=\sum _{i,j=1}^4{b}_{ij}{x}_i{x}_j=$$

$$\frac{1}{4}{b}_{11}{\left(x_1+\frac{4b_{12}}{b_{11}}{x}_2\right)}^2+\left[\frac{1}{3}{b}_{22}-\frac{4b_{12}^2}{b_{11}}\right]x_2^2+\frac{1}{4}{b}_{11}{\left(x_1+\frac{b_{13}}{b_{11}}{x}_3\right)}^2+$$

$$\left[\frac{1}{3}{b}_{33}-\frac{4b_{13}^2}{b_{11}}\right]x_3^2+b_{11}{\left(x_1+\frac{b_{14}}{b_{11}}{x}_4\right)}^2+\left[\frac{1}{3}{b}_{44}-\frac{4b_{14}^2}{b_{11}}\right]x_4^2+$$

$$\frac{1}{3}{b}_{22}{\left(x_2+3\frac{b_{23}}{b_{22}}{x}_3\right)}^2+\left[\frac{1}{3}{b}_{33}-\frac{9b_{23}^2}{b_{22}}\right]x_3^2+\frac{1}{3}{b}_{22}{\left(x_2+\frac{3b_{24}}{b_{22}}{x}_4\right)}^2+$$

$$\left[\frac{1}{3}{b}_{44}-\frac{9b_{24}^2}{b_{22}}\right]x_4^2+\frac{1}{3}{b}_{33}{\left(x_3+\frac{3b_{34}}{b_{33}}{x}_4\right)}^2+\left[\frac{1}{3}{b}_{44}-\frac{9b_{34}^2}{b_{33}}\right]x_4^2\ge 0,$$

when

$$\frac{1}{3}{b}_{22}\ge \frac{4b_{12}^2}{b_{11}},\frac{1}{3}{b}_{33}\ge \frac{4b_{13}^2}{b_{11}},\frac{1}{3}{b}_{44}\ge \frac{4b_{14}^2}{b_{11}},\frac{1}{3}{b}_{33}\ge \frac{9b_{23}^2}{b_{22}},$$

$$\frac{1}{3}{b}_{44}\ge \frac{9b_{24}^2}{b_{22}},\frac{1}{3}{b}_{44}\ge \frac{9b_{34}^2}{b_{33}}.$$

Then by reasoning as in Theorem 4. we obtain the conclusion. □

Remark 2. 

Assume the Condition 2 holds, $\Delta \ne 0$. Suppose $a_{ij}$ such that $b_{ii}>0$, $i=1,2,3,4$, $b_{ij}\ge 0$ for $i,j=1,2$ and $b_{ij}\le 0$ for $i,j=3,4$ when $i\ne j$ or $b_{ij}\le 0$ for $i,j=1,2$ and $b_{ij}\ge 0$ for $i,j=3,4$ when $i\ne j$.

Similar way as in Theorem 4 we get that Then the system $\left(3\right)$ is asymptotically stable at the equilibria point $E\left(\overline{x}\right)$ in the sense of Lyapunov.

5. Discussion

The spread rate of the virus in the population and its interaction with immune system cells were examined through a SEIR model. Stability analysis was performed by finding the equilibrium points of the nonlinear model created with four variables. Specifically, these variables $T=T\left(t\right),I=I\left(t\right),E\left(t\right),$ and $V\left(t\right)$were the concentration of uninfected cells, infected cells, effector immune cells, and free viruses, respectively. An approximate solution has been pursued since theoretical nonlinear equations usually have no deterministic solution. Here, the coefficient of the variables is considered equal to $0.1.$ In addition, global stability analysis was discussed by constructing the Lyapunov function, which validates the robustness of the dynamic model in disease-free equilibrium or endemic equilibrium. This analysis provides an understanding of the long-term behavior of the model and its impact on the progression of viral infections.

In Figure 1 and Figure 2, the interaction of the virus between infected cells and normal uninfected cells was studied. Figure 1 shows a direct proportionality between cells’ dynamics, while Figure 2 shows an inversion of the cells’ dynamics. Finally, Figure 3 shows the interaction between immune and infected cells. This Figure shows a direct proportionality between immune and infected cell dynamics.

The nonlinear model aimed to capture the interaction among the virus cells in different conditions, i.e., when virus cells are in contact with infected cells, when the virus is in relation with uninfected cells, or when there are infected and uninfected cells alone. It is worth noting that these relationships were obtained within certain mathematical conditions, which can affect the model results. The real-life condition could be more complex than those represented in this study. However, this model could pose a remarkably well-structured nonlinear mathematical fashion for analyzing how virus cells interact with the other cells within a human body. Moreover, it provides interesting information on the process dynamics over time, i.e., from the past to the future.

Future studies need to focus and validate the results of the model through experimental studies and physical observations. Although the model provides important views of the interaction dynamics between virus cells and other cells in the human body, real-life conditions may be more complex than those represented in this study. Therefore, empirical studies would be crucial to validate the accuracy and applicability of the model in a more real-life scenario.

6. Conclusion

In this work, the interaction between viruses and immune cells was investigated. Cases of normal cells infected with the virus were also included in the model. A four-variable dynamic model of normal cells, infected cells, effector immune cells, and free viruses was constructed. Equilibrium points of this mathematical model were found, and Lyapunov stability analyses were performed. Relationships among the virus, the infected and uninfected cells were discussed, as reported in Figure 1, Figure 2 and Figure 3.

A directly proportional relationship between viruses and infected cells was found in the decrease and increase of the cells’ dynamics, as shown in Figure 3. In Figure 2, virus and uninfected cells were compared. Rates of change are inversely proportional to each other. Infected and immune cells are proportionally reduced rapidly in Figure 3.

This model, which could be validated by clinical trials, poses itself as a possible aid for decision-makers to structure their health policies according to evidence-based mathematical models.