Figure 1.
Tumor eradication dynamics derived at $\gamma = 1.26$.
Citation: Robert P. Gilbert, Philippe Guyenne, Ying Liu. Modeling of the kinetics of vitamin D$_3$ in osteoblastic cells[J]. Mathematical Biosciences and Engineering, 2013, 10(2): 319-344. doi: 10.3934/mbe.2013.10.319
| [1] | Yoshi Nishitani, Chie Hosokawa, Yuko Mizuno-Matsumoto, Tomomitsu Miyoshi, Shinichi Tamura . Effect of correlating adjacent neurons for identifying communications: Feasibility experiment in a cultured neuronal network. AIMS Neuroscience, 2018, 5(1): 18-31. doi: 10.3934/Neuroscience.2018.1.18 |
| [2] | Yoshi Nishitani, Chie Hosokawa, Yuko Mizuno-Matsumoto, Tomomitsu Miyoshi, Shinichi Tamura . Classification of Spike Wave Propagations in a Cultured Neuronal Network: Investigating a Brain Communication Mechanism. AIMS Neuroscience, 2017, 4(1): 1-13. doi: 10.3934/Neuroscience.2017.1.1 |
| [3] | Shinichi Tamura, Yoshi Nishitani, Chie Hosokawa . Feasibility of Multiplex Communication in a 2D Mesh Asynchronous Neural Network with Fluctuations. AIMS Neuroscience, 2016, 3(4): 385-397. doi: 10.3934/Neuroscience.2016.4.385 |
| [4] | Shun Sakuma, Yuko Mizuno-Matsumoto, Yoshi Nishitani, Shinichi Tamura . Learning Times Required to Identify the Stimulated Position and Shortening of Propagation Path by Hebb’s Rule in Neural Network. AIMS Neuroscience, 2017, 4(4): 238-253. doi: 10.3934/Neuroscience.2017.4.238 |
| [5] | Shun Sakuma, Yuko Mizuno-Matsumoto, Yoshi Nishitani, Shinichi Tamura . Simulation of Spike Wave Propagation and Two-to-one Communication with Dynamic Time Warping. AIMS Neuroscience, 2016, 3(4): 474-486. doi: 10.3934/Neuroscience.2016.4.474 |
| [6] | Anastasios A. Mirisis, Anamaria Alexandrescu, Thomas J. Carew, Ashley M. Kopec . The Contribution of Spatial and Temporal Molecular Networks in the Induction of Long-term Memory and Its Underlying Synaptic Plasticity. AIMS Neuroscience, 2016, 3(3): 356-384. doi: 10.3934/Neuroscience.2016.3.356 |
| [7] | T. Virmani, F. J. Urbano, V. Bisagno, E. Garcia-Rill . The pedunculopontine nucleus: From posture and locomotion to neuroepigenetics. AIMS Neuroscience, 2019, 6(4): 219-230. doi: 10.3934/Neuroscience.2019.4.219 |
| [8] | Dirk Roosterman, Graeme S. Cottrell . Astrocytes and neurons communicate via a monocarboxylic acid shuttle. AIMS Neuroscience, 2020, 7(2): 94-106. doi: 10.3934/Neuroscience.2020007 |
| [9] | Eduardo Mercado III . Learning-Related Synaptic Reconfiguration in Hippocampal Networks: Memory Storage or Waveguide Tuning?. AIMS Neuroscience, 2015, 2(1): 28-34. doi: 10.3934/Neuroscience.2015.1.28 |
| [10] | Gigi Tevzadze, Tamar Barbakadze, Elisabed Kvergelidze, Elene Zhuravliova, Lali Shanshiashvili, David Mikeladze . Gut neurotoxin p-cresol induces brain-derived neurotrophic factor secretion and increases the expression of neurofilament subunits in PC-12 cells. AIMS Neuroscience, 2022, 9(1): 12-23. doi: 10.3934/Neuroscience.2022002 |
The elaboration of mathematical models that describe interactions between cancer cells and the immune system is an actively developing field in mathematical medicine. One of the important goals of this modeling is to analyze the ultimate dynamics which is necessary for monitoring the long-term tumor dynamics and, in the most optimistic case, for the creation of tumor eradication algorithms as well. Oncolytic virotherapy of cancer attracts attention of many researchers, see e.g. papers [1,2,3,4,5,6,8,14] and so on. Eftimie et al. in [5] described and studied the seven-dimensional model characterizing behavior of the cancer population in the case when the dual immunization protocol, by vaccine viruses and oncolytic viruses, is applied. They obtained local stability conditions for equilibrium points of their model and numerically investigated multi-stability phenomenon which may exist in the case of four types of equilibrium points. Later, in [6], Eftimie et al. presented the five-dimensional system which is obtained by a simplification of the model from [5]. This system was derived from an experimental protocol of the paper [3] and can be written as
|
$
⋅x1=rx1(1−x1+x2κ)−a1x1x5b1+x1−a2x1x4b2+x4⋅x2=a1x1x5b1+x1−δx2−a2x2x4b2+x4⋅x3=p1x3x5b3+x5(1−x3m)⋅x4=p2x3x5b3+x5−a3x4−a4x1x4+σ⋅x5=δβx2−γx5.
$
|
(1.1) |
Below the system (1.1) is considered on the nonnegative orthant ${\bf R}_{+, 0}^{5}$ and it is denoted by $x_{1}$/$x_{2}$/$x_{3}$ /$x_{4}$/ $x_{5}$ the density of uninfected tumor cells/ infected tumor cells/ memory cells/ effector immune cells/ oncolytic (Vesicular Stomatitis) virus particles (OVPs) correspondingly. Parameters are supposed to be positive excepting $\sigma $, $\sigma \geq 0\text{, }\; $ and are defined as follows: $r$ is proliferation rate for uninfected tumor cells; $\kappa $ is carrying capacity for tumor cells; $a_{1}$ is infection rate of tumor cells with the oncolytic virus; $b_{1}$ is half-saturation constant for the tumor cells infected with the oncolytic virus; $\delta $ is rate at which the oncolytic virus kills the tumor cells; $a_{2}$ is lysis rate of tumor cells (infected and uninfected) by the immune cells; $b_{2}$ is half-saturation constant for the effector cells that support half the maximum killing rate; $p_{1}$ is proliferation rate of memory cells following secondary encounter with tumor antigens carried by virus particles; $b_{3}$ is half-saturation constant of viral antigens that induce half the maximum proliferation rate of immune cells; $m$ is carrying capacity for memory cells; $p_{2}$ is rate at which memory cells become effector cells following secondary encounter with tumor antigens carried by virus particles; effector cells die at rate $a_{3}$; $a_{4}$ is inactivation rate of effector immune cells by tumor cells; $\beta $ is number of OVPs released from an infected cell, capable of forming plaques; $\gamma $ is decay rate for the concentration of OVPs in the blood; $\sigma $ is the influx rate into the lymphoid tissue at which the naive cells become effector cells, see [5]. It is pointed out that σ has sufficiently small values and may be neglected, [6]. In this paper it is assumed that the rate σ reflects the presence of the influx of immune cells from the external source and the dependence of global behavior of (1.1) is explored with a change in $\sigma $.
Authors of [6] studied only the case $\sigma = 0$. They found equilibrium points and performed local stability analysis. As well, authors of [6] fulfilled numerical studies transient and asymptotic behavior of the system which are based on using bifurcation diagrams. In particular, they showed the existence of periodic orbits and a chaotic attractor for some values of parameters. Nevertheless, the rigorous exploration of long time-dynamics is of interest for better understanding potentialities of cancer virotherapy described by equations (1.1).
The purpose of this paper and its novelty consist in studies of the case of convergence dynamics of (1.1), with $\sigma \geq 0$. In particular, conditions are found under which every trajectory in the nonnegative orthant tends to one of equilibrium points. It is worth to remind that systems with convergence dynamics have been investigated in articles of Hirsch; R.A. Smith; Li and Muldowney, see [7,15,16], in the broader perspective of studies of competitive or cooperative systems; higher dimensional generalizations of Bendixon's criterion and Dulac's criterion. Our results are based on the LaSalle theorem and the localization method of compact invariant sets (LMCIS), see [10,11]. It is worth to mention that the LMCIS has been efficiently utilized in papers [12,13,17,18,19,20] for analysis of ultimate dynamics of various cancer models.
The dynamics of (1.1) has some distinctive features: the unbounded set of equilibrium points, positive Lagrange stability of trajectories and several invariant planes that are useful in this analysis. In the case of a sufficiently large $\sigma >0$ it is found that (1.1) acquires a new global qualitative property consisting in the fact that all trajectories are attracted to the tumor-free plane $x_{1} = x_{2} = 0$. This property may be characterized quantitatively as well. Namely, it is computed a bound $\sigma _{\min }^{(1)}$ called the first eradication bound such that if $\sigma \geq \sigma _{\min }^{(1)}$ then each $\omega $-limit set is located in $x_{1} = x_{2} = 0$ plane. This signifies that the population of tumor cells is under control after sufficiently large time of observations. From the biological point of view, this case may correspond to the immunotherapy injection of some agent like IL-2 which is combined with virotherapy. Although the tumor may be cleared for $\sigma \geq \sigma _{\min }^{(1)}$, but the harmful side effects from an over-activated immune system may eliminate the benefits of tumor clearance. We propose an approach for reducing this risk. We show that there is the sequence of eradication bounds $\{\sigma _{\min }^{(k)}\}, k = 1, 2, .., $with the property that for $\sigma \geq \sigma _{\min }^{(k)}$ each $\omega $-limit set is located in $x_{1} = x_{2} = 0$ plane. The value $\sigma _{\min }^{(k)}$ appears as the minimal positive root of some algebraic equation of degree $k.$ The existence of $\sigma _{\min }^{(k)}$ is established in Theorem 5.1. The formula for $\sigma _{\min }^{(2)}$ is provided; here we have that $\sigma _{\min }^{(1)}>\sigma _{\min }^{(2)}$. It is demonstrated numerically that $\sigma _{\min }^{(2)}>\sigma _{\min }^{(3)}.$ We expect that $\{\sigma _{\min }^{(k)}\}$ is a decreasing sequence. Figs. 1-3 show the tumor eradication dynamics with $\sigma = \sigma _{\min }^{(i)}, i = 1, 2, 3.$
In the present paper ultimate upper bounds are computed for variables of (1.1). The domain defined by these bounds contains all compact invariant sets including a chaotic attractor found in [6] for some values of parameters. Further, it is noted that in the case $\sigma = 0$ tumor-free TF-equilibrium points are always locally unstable [6]. Main results of this paper are
1) Theorem 4.2 containing conditions under which the $\omega $-limit set of any positive half trajectory in ${\bf R}_{+}^{5}$ is either one of TF-equilibrium points or one of tumor-only TO-equilibrium points and
2) Propositions 3;4 and Theorem 5.1 in which we give minimal bounds for $\sigma $ guaranteeing convergence dynamics of any positive half trajectory in ${\bf R}_{+}^{5}$ to one of tumor-only TO-equilibrium points.
Therefore conditions of this paper may be utilized for a prediction of the possible health scenario for a sufficiently long time of observations.
The results of this work are stated in terms of inequalities imposed on parameters of the model. This may be important in applications, because in medical models, as a rule, biological parameters are known in a certain range.
The present paper is organized as follows. Helpful auxiliary results are provided in Section 2 and in Appendix. In Section 3 ultimate dynamics of the system (1.1) is studied via computing ultimate bounds in case $\sigma \geq 0$. These bounds form a domain in which the chaotic attractor is detected for some values of parameters given in [6]. Further, Section 4 contains conditions under which (1.1) possesses convergent dynamics to one of TO- or TF-equilibrium points in case $\sigma = 0$. In Section 5 we describe the case when the immunotherapy is utilized $\sigma >0$. Here we show that the tumor clearance phenomenon may be achieved with a sufficiently large value of $\sigma $. The tumor eradication bounds for $\sigma $ are derived. Section 6 contains concluding remarks. Appendix contains the derivation of the cubic equation for finding $\sigma _{\min }^{(3)}$.
With the goal of describing helpful results the following objects are introduced:
1) a nonlinear system
| $ \dot{x} = v (x) $ | (2.1) |
where $v$ is a $C^{1} -$differentiable vector field; $x \in {\bf R}^{n}$ is the state vector;
2) a $C^{1} -$differentiable function $h (x)$ such that $h$ is not the first integral of (2.1);
3) the set $S (h): = \{x \in {\bf R}^{n} \mid L_{v} h (x) = 0\}$ where $L_{v} h$ is the Lie derivative;
4) values $h_{\inf }: = \inf \{h (x) \mid x \in S (h)\}$; $h_{\sup }: = \sup \{h (x) \mid x \in S (h)\}\text{.}$
If it is set up the localization problem of compact invariant sets contained in the set $M \subset {\bf R}^{n}$ then the sets
|
$
S(h;M):=S(h)∩M={x∈M∣Lvh(x)=0};hinf(M):=inf{h(x)∣x∈S(h;M)};hsup(M):=sup{h(x)∣x∈S(h;M)}.
$
|
are defined.
Assertion 1 [9,10,11] For any $h(x)\in C^{1}({\bf R}^{n})$ all compact invariant sets of the system (2.1) located in $M$ are contained in the localization set $K(h; M)$ defined by the formula $K(h; M) = \{x\in M\mid h_{\inf }(M)\leq h(x)\leq h_{\sup }(M)\}$ as well. If $M\cap S(h) = \varnothing $ then there are no compact invariant sets located in $M$.
The function $h$ mentioned in this assertion is called localizing.
Below the following notations are used: if the system is examined on the positive orthant $M = {\bf R}_{+, 0}^{n}$ then $K(h): = K(h; M)\text{.}\; $ If the system is considered on the domain $M: = {\bf R}_{+, 0}^{n}\cap \{x_{5}>0\}$ then this $M$ is indicated.
Further, it is utilized
Assertion 2 [12]. Let $h$ be a differential function which is bounded below on the set $D\subset {\bf R}^{n}$ and $L_{v}h(x)\mid _{D}\leq -\varepsilon $ for some $\varepsilon >0$. Then the solution $\varphi (x, t)$ of the system (2.1), with initial condition $\varphi (x, 0) = x\in D\text{, }$ leaves the set $D$ in a finite time.
To formulate the following useful result, the nonautonomous system
| $ \dot{y} = G (y , t) , y \in {\bf R}^{n} , $ | (2.2) |
is introduced; here $G$ is a $C^{1}-$differentiable vector function. The system (2.2) is called asymptotically autonomous with the limit system (2.1) if $G(x, t)\rightarrow v(x)$, with $t\rightarrow\infty$, for any compact subset in ${\bf R}^{n}$, see e.g. [21].
Theorem 2.1. (Markus), see e.g. [21]. Let us consider systems (2.1) and (2.2) in case $n = 2.$ Now if $\Omega $ is the $\omega $ -limit set of a forward bounded solution $x$ of (2.2) then $\Omega $ either contains equilibria of (2.1) or $\Omega $ is the union of periodic orbits of (2.1).
The Markus theorem is applicable to the nonnegative orthant ${\bf R}_{+, 0}^{2}$ in case when ${\bf R}_{+, 0}^{2}$ is a positively invariant domain.
The planes $x_{1} = 0$; $x_{1} = x_{2} = 0$; $x_{3} = 0$; $x_{3} = m$; $x_{2} = x_{5} = 0$ are invariant.
Further, if $\sigma = 0$ then there are four types of equilibrium points, [6]:
a) TF-equilibrium points are given by ${\bf E}_{1}: = (0, 0, e_{1 3}, 0, 0)^{T}; e_{1 3} \in {\bf R}_{ +, 0}^{1};$
b) TO-equilibrium points are given by ${\bf E}_{2}: = (\kappa, 0, e_{2 3}, 0, 0)^{T}; e_{2 3} \in {\bf R}_{ +, 0}^{1};$
c) a tumor with virus but no immune response TV-equilibrium point is given by ${\bf E}_{3}: = (e_{3 1}, e_{3 2}, 0, 0, e_{3 5})^{T}$, with
|
$
e31=b1γa1β−γ;e32=rδβ(κ−e31)(b1+e31)a1δβκ+b1rγ+e31rγ;e35=e35γβδ;
$
|
d) an internal tumor virus ITV-equilibrium point is given by ${\bf E}_{4}: = (e_{41}, e_{42}, m, e_{44}, e_{45})^{T}$, with components $e_{4j}, j = 1, 2, 4, 5\text{, }\; $defined by very cumbersome formulas and omitted here.
Since TF-equilibrium points and TO-equilibrium points form an unbounded set the system (1.1) does not possess the bounded positively invariant domain.
It is noted that if $\sigma = 0$ then each TF-equilibrium point $(0, 0, e_{1 3}, 0, 0)^{T}$ is connected with TO-equilibrium point $(\kappa, 0, e_{1 3}, 0, 0)^{T}$ by a heteroclinic orbit which is located in the invariant plane $x_{i} = 0, i = 2, 4, 5;x_{3} = e_{13}$ and is an arc directed from $(0, 0, e_{13}, 0, 0)^{T}$ to $(\kappa, 0, e_{13}, 0, 0)^{T}$. So heteroclinic orbits cover half strip $x_{3} \geq 0, x_{2} = x_{4} = x_{5} = 0.$
Next, upper ultimate bounds are got for all state variables of the system (1.1) with help of computing the localization domain containing all compact invariant sets.
1. Bound $x_{1\max }$. The localizing function $h_{1} = x_{1}$ is applied. Then
|
$
S(h1)∩{x1>0}={r(1−x1+x2κ)=a1x5b1+x1+a2x4b2+x4≥0}⊂{r(1−x1+x2κ)≥0}.
$
|
Therefore the localization set $K_{1\max }: = \left \{0 \leq x_{1} \leq \kappa : = x_{1\max }\right \}$ is derived.
2) Bound $x_{2\max }\; \text{.}\; $ Here the function $h_{2} = x_{1} +x_{2}$ is utilized. Then
| $ L_{f} h_{2} = r x_{1} \left (1 -\frac{x_{1} +x_{2}}{\kappa }\right ) -a_{2} \frac{(x_{1} +x_{2}) x_{4}}{b_{2} +x_{4}} -\delta x_{2} $ |
Therefore the inequality
| $ \delta (x_{1} +x_{2}) \leq (\delta +r) x_{1} -\frac{r}{\kappa } x_{1}^{2} $ |
holds on $S (h_{2})\text{.}$ Considering the last inequality within $K_{1\max }$ the following inequality holds
|
$
\delta (x_{1} +x_{2}) \leq \zeta : = \left \{δκ,δ>rκ(δ+r)24r,δ≤r \right \}\;\text{.}\;
$
|
Thus the localization set $K_{21}: = \left \{h_{2} \leq h_{2\max }: = \zeta \delta ^{ -1}\right \}$ is derived. So the upper bound $x_{2\max }$ given in
| $ K_{2\max }: = \left \{x_{2} \leq \zeta \delta ^{ -1}: = x_{2\max }\right \} $ | (3.2) |
is obtained.
3) Bound $x_{5\max }$. Here using $h_{3} = x_{5}$ it is noted that
| $ S \left (h_{3}\right ) \cap K_{2\max } \subset \left \{x_{5} \leq \frac{\delta \beta }{\gamma } x_{2\max }\right \} $ |
and so the localization set with the bound $x_{5\max }$ is obtained:
| $ K_{5\max }: = \left \{x_{5} \leq x_{5\max }: = \frac{\delta \beta }{\gamma } x_{2\max }\right \} . $ | (3.3) |
4) The bound $x_{3\max }$ is computed with respect to the domain ${\bf R}_{ +, 0}^{5} \cap \{x_{5} >0\}\; \text{.}\; $ The function $h_{4} = x_{3}$ provides the localization set:
| $ K_{3\max } (h_{4} , {\bf R}_{ + , 0}^{5} \cap \{x_{5} \gt 0\}): = \{0 \leq x_{3} \leq m ;x_{5} \gt 0\}\;\text{.} $ |
5) The bound $x_{4\max }$ is computed with respect to the domain ${\bf R}_{ +, 0}^{5} \cap \{x_{5} >0\}\; \text{.}\; $ Here the function $h_{5} = x_{4}$ is applied. Then
| $ L_{f} h_{5} = p_{2} \frac{x_{3} x_{5}}{b_{3} +x_{5}} -a_{3} x_{4} -a_{4} x_{1} x_{4} +\sigma $ |
and therefore the set $S (h_{5}) \cap K_{5\max } \cap \{x_{5} >0\}$ is contained in the set defined by the inequality
| $ x_{4} \leq \frac{p_{2} x_{3} x_{5}}{a_{3} (b_{3} +x_{5})} +\frac{\sigma }{a_{3}} \leq x_{4\max }: = \frac{p_{2} m x_{5\max }}{a_{3} (b_{3} +x_{5\max })} +\frac{\sigma }{a_{3}}\;\text{.}\; $ |
Therefore the localization set:
| $ K_{4\max }(h_{5} , {\bf R}_{ + , 0}^{5} \cap \{x_{5} \gt 0\}): = \left \{0 \leq x_{4} \leq x_{4\max } ;x_{5} \gt 0\right\} $ |
is got.
As a result, it is established
Proposition 1. All compact invariant sets in${\bf R}_{ +, 0}^{5} \cap \{x_{5} >0\}$ are located in the domain
| $ \Pi : = \cap _{i = 1 , 2 , 5}K_{i\max } \cap K_{3\max } (h_{4} , {\bf R}_{ + , 0}^{5} \cap \{x_{5} \gt 0\}) \cap K_{4\max } (h_{5} , {\bf R}_{ + , 0}^{5} \cap \{x_{5} \gt 0\})\;\text{.}\; $ |
Let $\varphi (x, t)$ be a solution of (1.1), with $\varphi (x, 0) = x$ for all $x$. Further, the following fact is used in this paper:
Lemma 3.1. Let $\sigma \geq 0$. The system (1.1) is positively Lagrange stable in ${\bf R}_{+, 0}^{5}$, i.e. for any point $x\in {\bf R}_{+, 0}^{5}$ the time functions $\varphi _{i}(x, t), t\geq 0$, $i = 1, ..., 5\text{, }\; $ are upper bounded. So there are no escaping to infinity trajectories in positive time.
Proof. 1. Taking the function $h_{1} = x_{1}$ and the domain $D_{1} = \{x_{1}\geq \kappa +\varepsilon _{1}\}\cap {\bf R}_{+, 0}^{5}$ for $\varepsilon _{1}>0$ it is easy to see that
| $ L_{f}h_{1}\mid _{D_{1}}\leq -r\varepsilon _{1}-\frac{r\varepsilon _{1}^{2}}{\kappa } \lt 0\text{.}\; $ |
So applying Assertion 2 it is established that each trajectory leaves $D_{1} $ in a finite time. So $\varphi _{1}(x, t)$ is upper bounded as a time function.
2. Since $x_{3} = m$ is an invariant plane any trajectory with an initial condition satisfying the inequality $x_{3}\leq m$ has the property that $\varphi _{3}(x, t)$ is upper bounded as a time function. From the other side, if a trajectory has an initial condition satisfying $x_{3}>m$ then since $\dot{x}_{3}(t)\leq 0$ on the domain $x_{3}>m$ then $\varphi _{3}(x, t)$ is upper bounded as a time function.
3. We have in ${\bf R}_{+, 0}^{5}$ that
| $ \dot{x}_{4}\leq -a_{3}x_{4}+p_{2}\varphi _{3}(x, t)+\sigma . $ |
Integrating this inequality we obtain that the function $\varphi _{4}(x, t)$ is upper bounded.
4. Taking the function $h_{2} = x_{1}+x_{2}$ and the domain
| $ D_{2} = \{x_{1}+x_{2}\geq h_{2\max }+\varepsilon _{2};x_{1} \lt \kappa +\varepsilon _{1}\}\cap {\bf R}_{+, 0}^{5}, $ |
for $\varepsilon _{2}>0$ it is easy to see that
| $ L_{f}h_{2}\mid _{D_{2}}\leq -r\varepsilon _{1}-\frac{r\varepsilon _{1}^{2}}{\kappa }-a_{2}(h_{2\max }+\varepsilon _{2})\inf\limits_{t \gt 0}\frac{\varphi _{4}(x, t)}{a_{2}+\varphi _{4}(x, t)} \lt 0\text{.}\; $ |
So applying Assertion 2 it is established that each trajectory leaves $D_{2}$ in a finite time. So both of $\varphi _{i}(x, t), i = 1, 2\; \text{, }\; $ are upper bounded time functions.
5. At last, it follows from integration of $\dot{x}_{5}$ with the upper bounded time function $\varphi _{2}(x, t)$ that $\varphi _{5}(x, t)$ is the upper bounded time function as well.
The main purpose of this section is to present attraction conditions to either TF- equilibrium points or TO- equilibrium points in case $\sigma = 0$. All statements are established in the form of simple inequalities imposed on the decay rate $\gamma $ of VSV particles in the blood. Here $\gamma $ can be considered as a control parameter. Firstly, it is introduced a system obtained from (1.1) by a restriction on the invariant plane $x_{3} = x_{4} = 0$:
|
$
⋅x1=rx1(1−x1+x2κ)−a1x1x5b1+x1;⋅x2=a1x1x5b1+x1−δx2;⋅x5=δβx2−γx5.
$
|
(4.1) |
By $g$ it is denoted its vector field. This system has the following equilibrium points
|
$
e1=(0,0,0)T;e2=(κ,0,0)T;e3=(e31,e32,e33)T,
$
|
with
|
$
e31=b1γa1β−γ;e32=rb1γ(κ−e31)(e31+γ)rb1γ(e31+γ)+e31κa1βδ;e33=δβγe32.
$
|
Before considering the system (1.1) it is formulated
Lemma 4.1. The system (4.1) has the localization domain$\Pi _{1}$ in ${\bf R}_{+, 0}^{3}$ defined by
| $ x_{1}\leq x_{1\max }: = \kappa ;x_{2}\leq x_{2\max };x_{5}\leq x_{5\max }\text{, }\; $ |
where $x_{2\max }$ is defined by (3.2). This polytope contains all compact invariant sets and the attracting set of (4.1).
These bounds are derived as in case of Proposition 1.
Proposition 2. Assume that
| $ \gamma \gt \gamma_{min}: = \frac{a_{1}\beta \kappa }{\kappa +b_{1}}. $ | (4.2) |
Then each trajectory of the system (4.1) in ${\bf R}_{+, 0}^{3}$ tends to ${\bf e}_{1}$ or ${\bf e}_{2}\text{.}\; $
Proof. Let us apply the function $h_{6} = x_{2}+\beta ^{-1}x_{5}$. Then
| $ L_{g}h_{6} = (\frac{a_{1}x_{1}}{b_{1}+x_{1}}-\frac{\gamma }{\beta })x_{5}. $ | (4.3) |
It is noted that all $\omega $-limit sets of the system (4.1) are located in $\Pi _{1}$. So one may examine the equality (4.3) on $\Pi _{1}$ and proceed as follows
| $ L_{g}h_{6}\mid _{\Pi _{1}}\leq (\frac{a_{1}\kappa }{b_{1}+\kappa }-\frac{\gamma }{\beta })x_{5}\mid _{\Pi _{1}}\leq 0. $ |
By using Lemma 4.1 and the LaSalle theorem it is established that the $\omega $-limit set of any point of the system (4.1) is not empty and belongs to the set $\{x_{5} = 0\}\text{.}\; $Therefore, in order to analyze the ultimate dynamics of (4.1) it should be examined the ultimate dynamics of the system
|
$
⋅x1=rx1(1−x1+x2κ)−a1x1φ5(x,t)b1+x1⋅x2=a1x1φ5(x,t)b1+x1−δx2,
$
|
(4.4) |
with the input $\varphi _{5}(x, t)$ satisfying the condition
| $ \underset{t\rightarrow \infty }{\lim }\varphi _{5}(x, t) = 0. $ |
The system (4.4) is asymptotically autonomous with the limit system
|
$
⋅x1=rx1(1−x1+x2κ)⋅x2=−δx2,
$
|
(4.5) |
All solutions of the system (4.5) in ${\bf R}_{+, 0}^{2}$ are convergent to the equilibrium point $(\kappa, 0)^{T}$. In virtue of the Markus theorem, all solutions of the system (4.4) are convergent to the equilibrium point $(\kappa, 0)^{T}\text{.}\; $ Therefore any solution of the system (4.1) tends to ${\bf e}_{1}$ or ${\bf e}_{2}.$
Remark 1. If the inequality
| $ 1 \lt \frac{a_{1}\beta }{\gamma } \lt 1+\frac{b_{1}}{\kappa } $ |
then there is the third equilibrium ${\bf e}_{3}$ which is not contained in ${\bf R}_{ +, 0}^{3}$ because $\kappa -e_{3 1} < 0$ and, therefore, $e_{3 2} < 0$. It follows from formulas for equilibrium points.
Now it is established
Theorem 4.2. Suppose that (4.2) is fulfilled. Then the$\omega $ -limit set in ${\bf R}_{+, 0}^{5}$ is either one of TF-equilibrium points or one of TO-equilibrium points.
It should be noted that the condition (4.2) is the condition of local stability of TO-equilibrium points, see [6].
Proof. Let us apply the function $h_{6} = x_{2}+\beta ^{-1}x_{5}$. In this case
| $ L_{f}h_{6} = a_{1}\frac{x_{1}x_{5}}{b_{1}+x_{1}}-\frac{\gamma }{\beta }x_{5}-a_{2}\frac{x_{2}x_{4}}{b_{2}+x_{4}}. $ | (4.6) |
Since all $\omega $-limit sets of the system (1.1) in ${\bf R}_{+, 0}^{5}\cap \{x_{5}>0\}$ are nonempty and are located in $K_{1\max }$ one may restrict the equality (4.6) on $K_{1\max }$ and proceed as follows
| $ L_{f}h_{6}\mid _{K_{1\max }}\leq \{(\frac{a_{1}\kappa }{b_{1}+\kappa }-\frac{\gamma }{\beta })x_{5}-a_{2}\frac{x_{2}x_{4}}{b_{2}+x_{4}}\}\mid _{K_{1\max }}\leq 0. $ |
Using the condition (4.2) and the LaSalle theorem it is possible to conclude that any $\omega $-limit set belongs to the set $\{x_{2} = x_{5} = 0\}\cup \{x_{4} = x_{5} = 0\}$. Hence, any $\omega $-limit set is one of TF-equilibrium points or TO-equilibrium points.
Below everwhere we take the case $a_{2}>r$. We shall utilize the notation $x_{1\max }^{(1)}: = x_{1\max }$.
Proposition 3. Then all compact invariant sets are located in the domain $K_{4\min }^{(1)}\cap K_{1\max }^{(2)}$, with
|
$
K(1)4min:={x4≥x(1)4min:=σa3+a4κ};K(2)1max:={x1≤x(2)1max:=κ−a2κx(1)4minr(b2+x(1)4min)}.{}
$
|
Here these formulas are valid under condition
| $ \sigma \lt \sigma _{\min }^{(1)}: = \frac{b_{2}r(a_{3}+a_{4}\kappa )}{a_{2}-r}. $ | (5.1) |
2) Further, if
| $ \sigma \geq \sigma _{\min }^{(1)} $ | (5.2) |
then all $\omega $ -limit sets are located in the plane $x_{1} = 0.$ Moreover, all $\omega $ -limit sets are located in the plane $x_{1} = 0;x_{2} = 0$ as well.
The value of $\sigma _{\min }^{(1)}$ is called the first eradication bound.
Proof. 1). Using the function $h_{5}$ one may obtain inequalities
| $ \sigma \leq p_{2}\frac{x_{3}x_{5}}{b_{3}+x_{5}}+\sigma = a_{3}x_{4}+a_{4}x_{1}x_{4}\leq (a_{3}+a_{4}\kappa )x_{4} $ |
within the set $K_{1\max }\cap S(h_{5})$ and, as a result, the localization set $K_{4\min }^{(1)}$ is got. Next,
| $ \dot{x}_{1}\mid _{x_{1} \gt 0}\leq (r-\frac{r}{\kappa }x_{1}-\frac{a_{2}x_{4}}{b_{2}+x_{4}})\mid _{x_{1} \gt 0} $ | (5.3) |
and the inequality (5.3) on the domain $K_{4\min }^{(1)}$ leads to the inequality
| $ 0 \lt \frac{r}{\kappa }x_{1}\leq r-\frac{a_{2}x_{4\min }^{(1)}}{b_{2}+x_{4\min }^{(1)}}\;\text{.}\; $ |
This implies the desirable assertion.
2) It follows from (5.3), the LaSalle theorem and the fact that the second equation of (1.1) with $x_{1}\equiv 0$ is globally asymptotically stable.
Now we discuss how we can decrease the value for $\sigma $ for which the tumor clearance can be reached.
One may iterate this argument and obtain sequences of localization sets
|
$
K(1)4min⊂K(2)4min⊂…K(2)1max⊃K(3)1max⊃…,
$
|
with
| $ x_{1\max }^{(k)} (\sigma ): = \kappa -\frac{a_{2} \kappa x_{4\min }^{(k -1)} (\sigma )}{r (b_{2} +x_{4\min }^{(k -1)}) (\sigma )} , k = 3 , 4 , \ldots $ |
and
| $ x_{1\max }^{(k)}(\sigma ) \gt x_{1\max }^{(k+1)}(\sigma ), k = 1, 2, ... $ |
So similarly as above, we establish
Proposition 4. If for given value of the influx rate $\sigma $ there is $k$ such that
| $ x_{1\max }^{(1)}(\sigma ) \gt x_{1\max }^{(2)}(\sigma ) \gt ... \gt x_{1\max }^{(k-1)}(\sigma ) \gt 0 $ |
but $x_{1\max }^{(k)}(\sigma)\leq 0$ then we have the tumor clearance for this value of $\sigma $.
The minimal positive root of the equation (5.4)
| $ x_{1\max }^{(k)}(\sigma ) = 0 $ | (5.4) |
is called the $k-$th eradication bound which will be denoted by $\sigma _{\min }^{(k)}$. In what follows, we examine the issue concerning the existence of a positive root of (5.4) for any $k$. Finding bounds for the value $\sigma _{\min }^{(k)}$ encounters computational difficulties and is not the subject of this study. We notice that though the function $x_{1\max }^{(k)}(\sigma)$ may possess more than one positive root, this fact does not contradict the biological meaning of the use of the treatment $\sigma $. This is because in Propositions 3 and 4 only sufficient conditions of the tumor clearance are given.
Below we shall utilize notations
| $ A_{k}(\sigma ) = \frac{P_{k}(\sigma )}{Q_{k}(\sigma )} = x_{1\max }^{(k)}(\sigma );B_{k}(\sigma ) = x_{4\min }^{(k)}(\sigma ) $ | (5.5) |
and for simplicity we shall omit $\sigma $ for functions introduced in (5.5). Further, our goal is to establish
Theorem 5.1. The equation (5.4) for $k\geq 2$ has at least one positive root.
Proof. Firstly, we prove
Lemma 5.2. We have for any $k\geq 2$ that 1) $\deg P_{k} = \deg Q_{k} = k-1;$ 2) coefficients of $P_{k}$ and $Q_{k}$ satisfy the following condition:
|
$
coefσ0(Pk)>0;coefσ0(Qk)>0;coefσk−1(Pk)<0;coefσk−1(Qk)>0.
$
|
Proof. We shall prove this statement by induction on $k.$ For $k = 2$ this result is valid since
|
$
A2=P2Q2=κ−κa2B1r(b2+B1)=rb2κ(a3+a4κ)+σκ(r−a2)r[b2(a3+a4κ)+σ].
$
|
Suppose that this assertion is proven for $k.$ Let us consider this assertion for $k+1.$ Then
| $ A_{k+1} = \frac{P_{k+1}}{Q_{k+1}} = \kappa -\frac{\kappa a_{2}B_{k}}{r(b_{2}+B_{k})} = \frac{\kappa rb_{2}+B_{k}\kappa (r-a_{2})}{r(b_{2}+B_{k})}. $ | (5.6) |
Using the formula
| $ B_{k} = \frac{\sigma }{a_{3}+a_{4}A_{k}} = \frac{\sigma Q_{k}}{a_{3}Q_{k}+a_{4}P_{k}} $ |
we proceed (5.6) as follows:
| $ A_{k+1} = \frac{\kappa rb_{2}(a_{3}Q_{k}+a_{4}P_{k})+\sigma Q_{k}\kappa (r-a_{2})}{r[b_{2}(a_{3}Q_{k}+a_{4}P_{k})+\sigma Q_{k}]}. $ | (5.7) |
Now we get from the inductive hypothesis and the formula (5.7) that the first assertion of this lemma is valid and, besides, we notice that
|
$
coefσ0(Pk+1)=κrb2[a3coefσ0(Qk)+a4coefσ0(Pk)]>0;coefσ0(Qk+1)=rb2[a3coefσ0(Qk)+a4coefσ0(Pk)]>0;coefσk(Pk+1)=κ(r−a2)coefσk−1(Qk)<0;coefσk(Qk+1)=rcoefσk−1(Qk)>0.
$
|
Therefore Lemma 5.2 is proven.
As a result, in virtue of Lemma 5.2 and the formula (5.7) we obtain the desirable assertion.
Now we compute the value $\sigma _{\min }^{(2)}$. Then we have that
| $ B_{2} = \frac{\sigma r(b_{2}+B_{1})}{r(a_{3}+a_{4}\kappa )(b_{2}+B_{1})-a_{2}a_{4}\kappa B_{1}}. $ |
Examining the equation $A_{3}(\sigma) = 0$ we get the equation
| $ rb_{2} = B_{2}(a_{2}-r). $ | (5.8) |
For our convenience we introduce notations
|
$
α1=b2a2−r;α2=rb2(a3+a4κ)2;α3=r(a3+a4κ)−a2a4κ
$
|
in (5.8) and after this we come to the equation
| $ \sigma ^{2}+\sigma (b_{2}(a_{3}+a_{4}\kappa )-\alpha _{1}\alpha _{3})-\alpha _{1}\alpha _{2} = 0. $ |
Its unique positive root has the form
| $ \sigma _{+} = \frac{\alpha _{1}\alpha _{3}-b_{2}(a_{3}+a_{4}\kappa )}{2}+\sqrt{\frac{(\alpha _{1}\alpha _{3}-b_{2}(a_{3}+a_{4}\kappa ))^{2}}{4}+\alpha _{1}\alpha _{2}} $ |
We notice that $\sigma _{\min }^{(2)} = \sigma _{+}.$ Finally, by the routine computations omitted here it is easy to establish that
| $ \sigma _{\min }^{(1)}+\frac{b_{2}(a_{3}+a_{4}\kappa )-\alpha _{1}\alpha _{3}}{2} \gt 0 $ |
and, its square exceeds
| $ \frac{(\alpha _{1}\alpha _{3}-b_{2}(a_{3}+a_{4}\kappa ))^{2}}{4}+\alpha _{1}\alpha _{2}. $ |
So, we have that $\sigma _{\min }^{(1)}> \sigma _{\min }^{(2)}$ and in case of using $\sigma \geq \sigma _{\min }^{(2)}$ we still achieve the tumor eradication.
Now taking parameters from Table 2 [6], we demonstrate by calculations that the tumor eradication bounds $\sigma _{\min }^{(1)}; \sigma _{\min }^{(2)}; \sigma _{\min }^{(3)}$ are decreasing: $\sigma _{\min }^{(1)} = 1118.3;\sigma _{\min }^{(2)} = 530.3;\sigma _{\min }^{(3)} = 419.9$.
On Figures 1; 2; 3 we show the tumor eradication dynamics corresponding these values of $\sigma$ at $\gamma = 1.26; 1.46$ and $1.66:$
It should be said that one can get improved bounds for $x_{2}$ and $x_{5}$ with help of using formulas (3.2)-(3.3). In this case the function $h_{2} = x_{1}+x_{2}$ it should be employed in the same way as above in the formula (3.1). So considering $S(h_{2})$ within $K_{1\max }^{(1)}\cap K_{4\min }^{(1)}$ it is easy to see that this set is located in the set defined by the inequality
| $ (\delta +\frac{a_{2} x_{4\min }}{b_{2} +x_{4\min }}) (x_{1} +x_{2}) \leq \zeta \;\text{.}\; $ |
Thus the localization set $K_{2 2}: = \left \{h_{2} \leq h_{2\max }: = \zeta \delta _{1}^{ -1}\right \}$ is obtained, with
| $ \delta _{1} = \delta +\frac{a_{2} x_{4\min }}{b_{2} +x_{4\min }}\;\text{.}\; $ |
As a result, one can get upper bounds indicated in the formula
|
$
K(1)2max:={x2≤ζδ−11:=x(1)2max}, K(1)5max:={x5≤x(1)5max:=σβγx(1)2max}.
$
|
Hence one may conclude:
Proposition 5. All compact invariant sets are located in the domain $K_{1\max }^{(1)}\cap K_{2\max }^{(1)}\cap K_{5\max }^{(1)}\cap K_{4\min }^{(1)}$.
Convergence conditions for trajectories of the model (1.1) are presented in this paper. Such kind of ultimate dynamics may be important for predicting the patient's health scenario, its monitoring, correction and a proper application of treatments. For this purpose one may vary the parameters $\gamma $ (the decay rate for the concentration of OVPs in the blood) and $\sigma $ (the rate of influx of immune cells from the external source).
1) Ultimate upper bounds for $x_{1}$-, $x_{2}$-cells populations and for concentration of OVPs ($x_{5}$) are obtained. Upper bounds for $x_{3}$-, $x_{4}$-cells populations are given provided we examine compact invariant sets outside the plane $x_{5} = 0.$ As a result, it is shown that ultimate dynamics of the system (1.1) outside the plane $x_{5} = 0$ is contained in the domain defined by these bounds. Then it is established that in case $\sigma >0$ upper bounds for $x_{1}$-, $x_{2}$-, $x_{5}$- variables and the lower bound for $x_{4}$- variable can be improved by some number of iterations.
2) In the case $\sigma = 0$ we examine one specific type of the convergence dynamics (1.1), which excludes the existence of TV-, ITV-equilibrium points, periodic orbits and the chaotic attractor. All results are expressed in terms of $\gamma $. It is computed the lower bound $\gamma _{\min }$ such that for $\gamma >\gamma _{\min }$ any trajectory is attracted to either one of TO-equilibrium points or one of TF-equilibrium points. This problem is solved in Theorem 4.2.
3) In Proposition 3 conditions for the global tumor elimination are presented for sufficiently large values of the influx rate $\sigma $ of immune cells from the external source $\sigma, \sigma \geq \sigma _{\min }^{(1)}$. The formula (5.2) describes how the eradication bound $\sigma _{\min }^{(1)}$ depends on parameters of the model. It is shown that the bound $\sigma _{\min }^{(1)}$ may be decreased if we apply the iterational procedure. This leads to a reduced risk of over-immunization. As a result, we get the sequence of bounds $\sigma _{\min }^{(k)}$ which are characterized in Theorem 5.1. It is demonstrated that $\sigma _{\min }^{(2)} < \sigma _{\min }^{(1)}$. The decrease of the sequence of eradication bounds $\{\sigma _{\min }^{(k)}\}$ seems clear from the point of view of common sense and is confirmed by calculations. However, the additional analysis is still required and deserves a separate research. On Figures 1; 2; 3 it is shown numerically that there is the tumor eradication dynamics at different values of $\gamma $ and $\sigma = \sigma _{\min }^{(1)}; \sigma _{\min }^{(2)}; \sigma _{\min }^{(3)}\text{.}\; $Our statements about the tumor eradication are presented in a convenient form for applications and may be helpful in administering immunotherapy.
Further, let us derive an equation for finding $\sigma _{\min }^{(3)}$. With this goal we notice that $_{.}$
| $ A_{2} = \frac{\theta _{1}+\sigma \theta_{2}}{\theta _{3}+\sigma \theta_{4}}, $ | (6.1) |
with
|
$
θ1=κrb2(a3+a4κ);abcdθ2=κ(r−a2);θ3=rb2(a3+a4κ);abcdeθ4=r.
$
|
Utilizing these formulas for $B_{2}$ we have that
| $ B_{2} = \frac{\sigma }{a_{3}+a_{4}A_{2}} = \frac{\theta _{3}\sigma +\theta _{4}\sigma ^{2}}{\theta _{5}+\theta _{6}\sigma }, $ |
with
| $ \theta _{5} = a_{3}\theta _{3}+a_{4}\theta _{1}; \phantom{abcd}\theta _{6} = a_{3}\theta _{4}+a_{4}\theta _{2}. $ |
Next, we compute
| $ A_{3} = \kappa -\frac{\kappa a_{2}B_{2}}{rb_{2}+rB_{2}} = \kappa \frac{\theta _{7}+\sigma \theta_{8}+\sigma ^{2}\theta_{9}}{\theta _{10}+\sigma \theta _{11}+\sigma ^{2}\theta _{12}}, $ |
with
|
$
θ7=rb2θ5;abcdθ8=rb2θ6+(r−a2)θ3;abcdθ9=(r−a2)θ4;abcdθ10=rb2θ6+rθ3;abcdθ11=rθ4.
$
|
Further,
| $ B_{3} = \frac{\sigma }{a_{3}+a_{4}A_{3}} = \frac{\sigma \theta _{7}+\sigma ^{2}\theta _{10}+\sigma ^{3}\theta _{11}}{\theta_{12}+\sigma \theta _{13}+\sigma ^{2}\theta _{14}}, $ |
with
| $ \theta_{12} = (a_{3}+a_{4}\kappa )\theta _{7};\theta _{13} = a_{3}\theta _{10}+a_{4}\kappa \theta _{8};\theta _{14} = a_{3}\theta _{11}+a_{4}\theta _{9} $ |
Finally, we notice that $\sigma _{\min }^{(3)}$ is a minimal positive root of the cubic equation $rb_{2}+B_{3}(r-a_{2}) = 0$ which can be written in the form
| $ \sigma ^{3}\theta _{11}(r-a_{2})+\sigma ^{2}(rb_{2}\theta _{14}+(r-a_{2})\theta _{10})+\sigma (rb_{2}\theta _{13}+(r-a_{2})\theta _{7})+rb_{2}\theta _{12} = 0. $ |
All authors declare no conflicts of interest in this paper.
| [1] | Biophys. J., 65 (1993), 1727-1739. |
| [2] | Comp. Biochem. Physiol., 115B (1996), 313-317. |
| [3] | Submitted to J. Theor. Biol., (2012). |
| [4] | Endocrinology, 120 (1987), 1173-1178. |
| [5] | J. Pharm. Exp. Therap., 307 (2003), 839-845. |
| [6] | Math. Mehods Appl. Sciences, 33 (2010), 2206-2214. |
| [7] | Biomechan. Model. Mechanobiol., 9 (2010), 87-102. |
| [8] | Springer, Berlin 1994. |
| [9] | Springer, Berlin, 1998. |
| [10] | Bone, 33 (2003), 206-215. |
| [11] | Cell, 93 (1998), 165-176. |
| [12] | Oxford University Press, New York, 1996. |
| [13] | J. Theor. Biol., 229 (2004), 293-309. |
| [14] | Bone Miner. Res., 14 (1999), 1543-1549. |
| [15] | J. Biol. Chem., 285 (2010), 31859-31866. |
| [16] | J. Cell Biochem., 88 (2003), 438-445. |
| [17] | Academic Press, New York, 1997. |
| [18] | Biophys. J., 97 (2009), 758-767. |
| [19] | J. Chem. Phys., 134 (2011), 154104. |
| [20] | Academic Press, New York, 1997, 105-125. |
| [21] | Cell, 89 (1997), 309-319. |
| [22] | J. Cell Biol., 179 (2007), 635-641. |
| [23] | Academic Press, New York, 1997, 219-261. |
Robert P. Gilbert, Philippe Guyenne, Ying Liu. Modeling of the kinetics of vitamin D$_3$ in osteoblastic cells[J]. Mathematical Biosciences and Engineering, 2013, 10(2): 319-344. doi: 10.3934/mbe.2013.10.319
DownLoad:
