Stability of a retrovirus dymanic model

TL;DR

This paper introduces a retrovirus dynamic model focusing on low pathogenicity and infected cell reproduction, analyzing the stability of a chronic disease state using Lyapunov functions.

Contribution

It presents a new mathematical model for retrovirus dynamics and applies Lyapunov methods to study stability of the chronic infection equilibrium.

Findings

Inner equilibrium is stable under certain conditions.

The model provides insights into chronic retrovirus infections.

Lyapunov functions effectively analyze stability in this context.

Abstract

A retrovirus dynamic model is proposed. We pay attention to the case when viral pathogenicity is low and the infected cells are able to reproduce. Using Lyapunov function method we study stability properties of an inner equilibrium of the model. The equilibrium represents a chronic disease steady state.

Equations

$$\left\{\begin{array}[]{l}\frac{dC(t)}{dt}=\lambda-dC(t)-\beta C(t)V(t),\\ \frac{dI(t)}{dt}=\beta C(t)V(t)-aI(t),\\ \frac{dV(t)}{dt}=aN\,I(t)-kV(t),\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\frac{dC(t)}{dt}=aC(t)\left(1-b_{11}C(t)-b_{12}I(t)\right)-\alpha\,C(t)V(t),\\ \frac{dI(t)}{dt}=a_{I}I(t)\left(1-b_{21}C(t)-b_{22}I(t)\right)+\alpha\,C(t)V(t)-mI(t),\\ \frac{dV(t)}{dt}=kmI(t)-\sigma V(t),\\ \end{array}\right.$$

$$C(0)=C_0\ge 0,I(0)=I_0\ge 0,V(0)=V_0\ge 0,$$

$$\left\{\begin{array}[]{l}0=aC(t)\left(1-b_{11}C(t)-b_{12}I(t)\right)-\alpha\,C(t)V(t),\\ 0=a_{I}I(t)\left(1-b_{21}C(t)-b_{22}I(t)\right)+\alpha\,C(t)V(t)-mI(t),\\ 0=kmI(t)-\sigma V(t).\\ \end{array}\right.$$

$$\frac{1}{\hat{I}}\left[a_I{b}_{22}-\frac{1}{\hat{I}}{a}_I\left(1-b_{21}\hat{C}-b_{22}\hat{I}\right)+m\right]\cdot \frac{1}{\hat{C}}\left[a{b}_{11}-\frac{1}{\hat{C}}a\left(1-b_{11}\hat{C}-b_{12}\hat{I}\right)\right]>$$

$$>\frac{1}{4}{\left(\frac{a{b}_{12}}{\hat{C}}+\frac{b_{21}}{\hat{I}}-{\hat{V}}^2\right)}^2.$$

$$W(C,I,V)\equiv A\cdot v\left(\frac{C}{\hat{C}}\right)+B\cdot v\left(\frac{I}{\hat{I}}\right)+D\cdot v\left(\frac{V}{\hat{V}}\right),$$

$${\dot{W}}_{(2)}(t)=A\left(1-\frac{\hat{C}}{C(t)}\right)\frac{1}{\hat{C}}\left(aC(t)\left(1-b_{11}C(t)-b_{12}I(t)\right)-\alpha C(t)V(t)\right)$$

$$+B\left(1-\frac{\hat{I}}{I(t)}\right)\frac{1}{\hat{I}}\left(a_{I}I(t)\left(1-b_{21}C(t)-b_{22}I(t)\right)+\alpha C(t)V(t)-mI(t)\right)$$

$$+D\left(1-\frac{\hat{V}}{V(t)}\right)\frac{1}{\hat{V}}\left(kmI(t)-\sigma V(t)\right).$$

$$S_1=A\frac{C-\hat{C}}{C\hat{C}}\left(aC\left(1-b_{11}C-b_{12}I\right)-\alpha CV-a\hat{C}\left(1-b_{11}\hat{C}-b_{12}\hat{I}\right)+\alpha \hat{C}\hat{V}\right).$$

$$S_1=-(C-\hat{C}{)}^2\cdot \frac{A}{\hat{C}}\left[a{b}_{11}-\frac{1}{C}a\left(1-b_{11}\hat{C}-b_{12}\hat{I}\right)+\frac{1}{C}\alpha \hat{V}\right]$$

$$-(C-\hat{C})(I-\hat{I})\cdot \frac{Aa{b}_{12}}{\hat{C}}-(C-\hat{C})(V-\hat{V})\cdot \frac{A\alpha }{\hat{C}}.$$

$$S_2=B\frac{I-\hat{I}}{I\cdot \hat{I}}\left[a_{I}I(1-b_{21}C-b_{22}I)+\alpha CV-mI\underbrace{-a_I\hat{I}(1-b_{21}\hat{C}-b_{22}\hat{I})-\alpha \hat{C}\hat{V}+m\hat{I}}\right]$$

$$=B\frac{I-\hat{I}}{I\cdot \hat{I}}[a_{I}I\left(b_{21}(\hat{C}-C)+b_{22}(\hat{I}-I)\right)+(I-\hat{I})a_I(1-b_{21}\hat{C}-b_{22}\hat{I})-m(I-\hat{I})$$

$$-(\hat{V}-V)C\alpha +(C-\hat{C})\alpha \hat{V}]$$

$$=-(I-\hat{I}{)}^2\cdot \frac{B}{\hat{I}}\left[a_I{b}_{22}-\frac{1}{I}{a}_I\left(1-b_{21}\hat{C}-b_{22}\hat{I}\right)+m\right]$$

$$=-(C-\hat{C})(I-\hat{I})\cdot \left(\frac{B{b}_{21}}{\hat{I}}-\alpha \hat{V}\right)+(V-\hat{V})(I-\hat{I})\cdot \frac{C}{I}\frac{B\alpha }{\hat{I}}.$$

$$S_3=-D\frac{V-\hat{V}}{V\hat{V}}\cdot \left[kmI-\sigma V\underbrace{-km\hat{I}+\sigma \hat{V}}\right]$$

$$=D\frac{V-\hat{V}}{V\hat{V}}\cdot \left[km(I-\hat{I})-\sigma (V-\hat{V})\right]=-(V-\hat{V}{)}^2\cdot \frac{1}{V}\frac{D\sigma }{\hat{V}}+(V-\hat{V})(I-\hat{I})\cdot \frac{1}{V}\frac{Dkm}{\hat{V}}.$$

$$-{\dot{W}}_{(2)}=(V-\hat{V}{)}^2\cdot {\omega }_{11}+(I-\hat{I}{)}^2\cdot {\omega }_{22}+(C-\hat{C}{)}^2\cdot {\omega }_{33}$$

$$+(V-\hat{V})(I-\hat{I})\cdot 2{\omega }_{12}+(V-\hat{V})(C-\hat{C})\cdot 2{\omega }_{13}+(C-\hat{C})(I-\hat{I})\cdot 2{\omega }_{21},$$

$$\overline{W_{(2)}}=(V-\hat{V}{)}^2\cdot \overline{{\omega }_{11}}+(I-\hat{I}{)}^2\cdot \overline{{\omega }_{22}}+(C-\hat{C}{)}^2\cdot \overline{{\omega }_{33}}$$

$$+(V-\hat{V})(I-\hat{I})\cdot 2\overline{{\omega }_{12}}+(V-\hat{V})(C-\hat{C})\cdot 2\overline{{\omega }_{13}}+(C-\hat{C})(I-\hat{I})\cdot 2\overline{{\omega }_{23}}.$$

$${\Delta }_2=\frac{1}{\overline{V}}\frac{D\sigma }{\hat{V}}\cdot \frac{B}{\hat{I}}\left[a_I{b}_{22}-\frac{1}{\overline{I}}{a}_I\left(1-b_{21}\hat{C}-b_{22}\hat{I}\right)+m\right]-\frac{1}{4}{\left(\frac{1}{\overline{V}}\cdot \frac{Dkm}{\hat{V}}+\frac{\overline{C}}{\overline{I}}\cdot \frac{B\alpha }{\hat{I}}\right)}^2.$$

$${\Delta }_2(\hat{C},\hat{I},\hat{V})=\frac{1}{\hat{V}}\frac{D\sigma }{\hat{V}}\cdot \frac{B}{\hat{I}}\left[a_I{b}_{22}-\frac{1}{\hat{I}}{a}_I\left(1-b_{21}\hat{C}-b_{22}\hat{I}\right)+m\right]$$

$$-\frac{1}{4}{\left(\frac{1}{\hat{V}}\cdot \frac{Dkm}{\hat{V}}+\frac{\hat{C}}{\hat{I}}\cdot \frac{B\alpha }{\hat{I}}\right)}^2$$

$${\Delta }_3=\overline{{\omega }_{11}}\cdot (\overline{{\omega }_{22}}\cdot \overline{{\omega }_{33}}-{\overline{{\omega }_{23}}}^2)-\overline{{\omega }_{12}}\cdot (\overline{{\omega }_{12}}\cdot \overline{{\omega }_{33}}-\overline{{\omega }_{23}}\cdot \overline{{\omega }_{13}})+\overline{{\omega }_{13}}\cdot (\overline{{\omega }_{12}}\cdot \overline{{\omega }_{23}}-\overline{{\omega }_{22}}\cdot \overline{{\omega }_{13}}).$$

Full text

††margin:

2018.12.29.AR

2010 Mathematics Subject Classification: 34D20; 93D05; 92B05.

Keywords: Evolution equations; Lyapunov stability; virus infection model;

anticancer virotherapy.

1 Introduction

There are two major modes of viral replication: lytic and non-lytic. The dynamics of lytic virus is considered in many papers (see, e.g. [9, 10]). Non-lutic is almost out of consideration. In the nonlytic reproduction virus stays dormant in an infected cell, that reproduce. Virus stays dormant until the cell stars to exhibit signs of exhaustion. When the fist such signs appear, virus starts fast replication in the cell, kills it and brakes the cell membrane releasing virus particles.

The idea of virotherapy consist in using viruses in delivery their genetic material (a piece of RNA or DNA) into host cells that need to be treated. There are several branches of virotherapy such as viral immunotherapy, anti-cancer oncolytic viruses and viral vectors for gene therapy. Anticancer virotherapy is a new and promising method of anticancer therapy.

However, the majority of viruses that are used for virotherapy are of highly virulent and highly pathogenic type, and, therefore, the infected cells usually do not reproduce (producing the virus instead). Typical mathematical models reflect this fact.

The following classical model of virus dynamics was proposed in [9, 10]

[TABLE]

where $C(t),I(t),V(t)$ represent the concentration (or the total number) of non-infected host cells, infected cells and free virions at time $t$, respectively. All the constants $\lambda,d,\beta,a,N,k$ in (1) are positive. The non-infected cells are produced at rate $\lambda$, die at rate $d$ and become infected at rate $\beta$. Infected cells die at rate $a$. Free virus is produced by infected cells at rate $aN$ and die at rate $k$. This model was extended in many directions including introduction more general nonlinear terms [4, 5], time delays [3, 16, 15, 11, 12, 13, 14], additional equations describing immune responses [15, 16, 20, 21, 22, 23, 24, 11, 12] and inhomogenuous in space terms (which lead to partial differential equations models) [17, 18, 19, 13, 14] (see also references therein for more information).

In this paper we, in contrast, consider a situation when viral pathogenicity is low and the infected cells are able to reproduce. Such a situation arise, for instance, when retrovirus are used for the therapy. This assumption leads to a principally different mathematical model and principally different outcomes. We propose the following virus dynamics model

[TABLE]

where unknowns $C(t),I(t),V(t)$ are as in the model (1) above.

We are interested in the stability properties of an equilibrium (a stationary solution) of the model (2). For the general Lyapunov stability theory see the original work [7].

2 Basic properties

First, one sees that for any non-negative initial data

[TABLE]

the system (2) has a unique global (defined for all $t\geq 0$) solution. Each coordinate is non-negative for all $t\geq 0$, which is a biologically important property of the model. It follows from the standard property $\frac{dC(t)}{dt}|_{t=\tau}\geq 0$ provided $C(\tau)=0.$ Similar properties are valid for $I(t)$ and $V(t)$. In the similar standard way one shows that any solution is bounded. Moreover there is a bounded invariant region in $\mathbb{R}^{3}$.

The next step of our study, is to look for possible equilibria of the model (2).

2.1 Sationary solutions

We are interested in stationary solutions of the system (2). As a constant in-time solutions, they satisfy

[TABLE]

It is easy to see that there is a unique stationary solution of (3) such that all the coordinates are positive (inner equilibrium). In this note we are interested in this inner equilibrium and do not discuss boundary stationary solutions (when at least one coordinate is zero). Let us denote this unique inner solution of (3) by $(\widehat{C},\widehat{I},\widehat{V})$.

Remark 2.1

We believe this stationary solution is the most important one from the biological point of view. The equilibrium represents a chronic disease steady state.

2.2 Lyapunov stability

Assume the parameters of the system (2) satisfy

[TABLE]

[TABLE]

Now we formulate the main result.

Theorem 2.2

Let the condition (4) be satisfied. Then there exist $\alpha_{0}>0,k_{0}>0$ such that for any $\alpha\in(0,\alpha_{0}),$ and $k\in(0,k_{0})$ the inner equilibrium $(\widehat{C},\widehat{I},\widehat{V})$ of the system (2) is locally asymptotically stable.

Proof of theorem 2.2. We use the Volterra function $v(s)\equiv s-\ln\,s-1$ to construct the following Lyapunov functional

[TABLE]

where $A,B,D$ are positive constants to be chosen below.

It is easy to check that $W(C,I,V)>0$ for all $(C,I,V)\neq(\widehat{C},\widehat{I},\widehat{V})$ and $W(\widehat{C},\widehat{I},\widehat{V})=0$.

We denote by $\dot{W}_{(2)}$ the derivative of $W$ along a solution of the system (2), which is $\dot{W}_{(2)}(t)=\frac{d}{dt}W(C(t),I(t),V(t))$. As usual, it is computed using the right-hand side of the system (2) and the property $\frac{d}{dt}v(s)=1-s^{-1}.$

We have

[TABLE]

[TABLE]

[TABLE]

We split the above sum (6) on there parts $\dot{W}_{(2)}(t)=S_{1}(t)+S_{2}(t)+S_{3}(t)$ to estimate them separately. We omit the time argument for short.

Using $-a\widehat{C}\left(1-b_{11}\widehat{C}-b_{12}\widehat{I}\right)+\alpha\,\widehat{C}\widehat{V}=0$, one can check that

[TABLE]

Some calculations give

[TABLE]

[TABLE]

In a similar way, using $a_{I}\widehat{I}\left(1-b_{21}\widehat{C}-b_{22}\widehat{I}\right)+\alpha\,\widehat{C}\widehat{V}-m\widehat{I}=0$ (see (3)), calculations give

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The third step of calculations give (remind that $\sigma\widehat{V}=km{\widehat{I}}$ )

[TABLE]

[TABLE]

Hence, combining (7), (8), (9), we arrive to the form of $\dot{W}_{(2)}(t)=S_{1}(t)+S_{2}(t)+S_{3}(t)$. Now our goal is to find sufficient local conditions for $\dot{W}_{(2)}(t)<0$ for all points (except the inner equilibrium $(\widehat{C},\widehat{I},\widehat{V}),$ where $\dot{W}_{(2)}(t)\equiv 0$). Our main idea is to compare the expression for $\dot{W}_{(2)}$ with auxiliary quadratic forms. Hence, knowing the conditions for the quadratic form to be positive (negative) defined, we find sufficient local conditions for $\dot{W}_{(2)}$ to be negative defined as a functional in a neighbourhood of $(\widehat{C},\widehat{I},\widehat{V})$.

Formulas (7), (8), (9) are prepared to write $\dot{W}_{(2)}$ in the following form (notice the sign)

[TABLE]

[TABLE]

where all functions $\omega_{ij}$ depend on coordinates $(C,I,V)$, parameters of the system, coordinates of the equilibrium $(\widehat{C},\widehat{I},\widehat{V})$ and coefficients $A,B,D$. Because of the dependence of $\omega_{ij}$ on coordinates $(C,I,V)$ the expression (10) is not a quadratic form. We proceed as follows. Let us choose an arbitrary point $(\overline{C},\overline{I},\overline{V})$ in a small neighbourhood of $(\widehat{C},\widehat{I},\widehat{V})$. We denote the values of $\omega_{ij}$ at this fixed point as $\overline{\omega_{ij}}\equiv\omega_{ij}(\overline{C},\overline{I},\overline{V})$ and arrive to the quadratic form

[TABLE]

[TABLE]

This auxiliary quadratic form $\overline{W_{(2)}}$ is designed in such a way that $\overline{W_{(2)}}(\overline{C},\overline{I},\overline{V})=-\dot{W}_{(2)}(\overline{C},\overline{I},\overline{V})$. Hence, a condition to be positive defined for $\overline{W_{(2)}}$ implies the desired property $\dot{W}_{(2)}(\overline{C},\overline{I},\overline{V})<0$. Since point $(\overline{C},\overline{I},\overline{V})$ is arbitrary chosen, we arrive to the local stability result.

We apply the classical Sylvester’s criterion to the quadratic form $\overline{W_{(2)}}$.

Let us outline the main technical idea to show that the form is positive defined about $(\widehat{C},\widehat{I},\widehat{V})$. The first leading principal minor of the quadratic form (10) reads $\Delta_{1}\equiv\overline{\omega_{11}}=\frac{1}{\overline{V}}\frac{D\sigma}{\widehat{V}}>0.$ Let us consider the second leading principal minor $\Delta_{2}=\overline{\omega_{11}}\,\cdot\overline{\omega_{22}}-\overline{\omega_{12}}^{2}.$ It reads

[TABLE]

We remind that we consider the inner equilibrium, so $\Delta_{2}\equiv\Delta_{2}(\overline{C},\overline{I},\overline{V})$ is continuous in a small enough neighbourhood of $(\widehat{C},\widehat{I},\widehat{V})$. Hence, the property $\Delta_{2}(\widehat{C},\widehat{I},\widehat{V})>0$ will guarantee the property $\Delta_{2}(\overline{C},\overline{I},\overline{V})>0$ in a small enough neighbourhood. We write down the expression of

[TABLE]

[TABLE]

and see that the property $\Delta_{2}(\widehat{C},\widehat{I},\widehat{V})>0$ can be reached by choosing small enough $k$ and $\alpha$.

The same line of arguments is used to show that the third leading principal minor $\Delta_{3}$ is positive. We write it in a short form as follows

[TABLE]

The detailed analysis of these three terms of $\Delta_{3}$ shows that the first term $\overline{\omega_{11}}\,\cdot(\overline{\omega_{22}}\cdot\overline{\omega_{33}}-\overline{\omega_{23}}^{2})$ can be made positive provided condition (4) is satisfied. Moreover, $\Delta_{3}$ is positive provided $k$ and $\alpha$ are small enough. We omit the calculation here.

It completes the proof of theorem 2.2.

Remark 2.3

We notice that the smallness of parameters $k$ and $\alpha$ is not the only possible way to reach the stability. We choose this case just because of the clear biological meaning.

Acknowledgement. This work was supported in part by GA CR under project 16-06678S.