Dynamic analysis of an age structure model for oncolytic virus therapy

1.   Introduction

According to the latest statistics published in the American Cancer Society's the American Cancer Society Clinical Cancer Journal ^[1], cancer affects 18 million people and causes 9.6 million deaths worldwide each year. Cancer is one of the world's deadliest diseases. Unfortunately, the number of deaths from cancer is increasing rapidly. Despite the best efforts of scientists to research traditional cancer treatment methods, they still have problems with large side effects and low efficacy, and cancer treatment is still an underdeveloped science. Current studies and clinical trials have shown that oncolytic viruses have gradually shown their therapeutic advantages due to their high tumor-killing efficiency, good targeting, and few side effects ^[2]. The oncolytic viruses genetically engineered not only kill tumor cells but also activate the body's defenses by stimulating the immune system, attracting more immune cells to kill the remaining cancer cells without destroying the body's normal cells. Oncolytic therapy is a safe strategy for biological the immunotherapy of tumors ^[3,4,5].

Oncolytic viruses first appeared in the early 20th century. Ribacka et al.^[6] reported for the first time a case of tumor regression in a cervical cancer patient infected with the rabies virus. In 1991, literature ^[7] showed that transgenic HSV played a certain role in the treatment of glioblastoma. In the 21st century, at least three oncolytic therapeutics have been approved to enter the global market, namely, Regvir (ECHO-7 virus), Enco (Regubi, and T-VEC (herpes simplex virus) ^[8]. In recent decades, research on oncolytic therapy has made remarkable progress in the medical field ^[9,10,11].

In addition, mathematicians use clinical trial results and associated data analysis to establish mathematical models to analyze the dynamics of oncolytic virotherapy both in the short and long term ^[12,13,14], potentially providing key guidance for cancer treatment. Some studies incorporated viruses into ordinary differential models at a constant rate in an attempt to explore the dynamic properties of tumor cells and their associated immune responses but ignored studies of free viruses ^[15,16,17]. Some studies have described complex communication between tumor cells and viruses, but have ignored the immune response of tumor cells ^{[18,19,20,21]}. There were also oncolytic virus therapy models that demonstrate interactions between tumor cells, effector T cells, and virions ^[22,23]. In 2001, Wodarz et al.^[24] first proposed a mathematical model for oncolytic virus therapy, using two sets of ordinary differential equations to simulate virus interactions with tumor cells. Although virus populations were modeled implicitly, these mathematical models provide useful insights into the dynamic properties between oncolytic viruses and tumor cell populations. In 2008, Bajzer et al.^[25] the first time proposed a model that explicitly simulated virus population and improved the logistic growth of tumor cells, and conducted bifurcation analysis on the system, obtaining the stable oscillation of Hopf bifurcation and further developing the mathematical model of oncolytic therapy. In 2017, Titze et al.^[26] based on the literature ^[25] model, provided an adenovirus treatment containing three different soluble tumor glioblastoma cell growth under the ordinary differential models to the characterization of tumor cell growth and the relationship between the different soluble tumor viruses, by optimizing the model data, obtained the parameter estimate and predicted the long-term tumor recurrence. In 2021, Vithanage, G. et al.^[27] proposed a mathematical model of the interaction between the tumor-immune system and oncolytic virus therapy, studying the bistability of the immune system with immunosuppression in a large parameter space. The results showed that tumors of all sizes can be eliminated when the tumor killing rate of immune cells exceeds a critical value, and the critical killing rate of tumors depended on whether the immune system was immunosuppressive or proliferative. In 2022, Xiao et al.^[28] proposed a delay differential equation that specifies the activation time of cytotoxic T lymphocyte responses, and it was shown that this delay time is critical for oncolytic therapy for cancer research. Khaphetsi Joseph Mahasa et al.^[29] developed an ordinary differential model to investigate the main question of whether single or multiple sequential doses of CAR T cells during oncolytic virus therapy can have a synergistic effect on tumor reduction. The results showed that the simultaneous administration of oncolytic viruses and CAR T cells had a stronger synergistic effect on tumor cell reduction than the combination of first or later administration of CAR T cells. Therefore, it can be observed that previous models of oncolytic therapy studies only considered changes in cell density with time ^[27,29,30].

In recent years, the age structure model of viral infection has attracted extensive attention from researchers because the mortality of infected cells varies with their lifespan ^[31,32,33]. Data showed that the rate of production of the virus is initially low but increases as infected cells age ^[34]. So it is appropriate to assume that the productivity of the virus and the death rate of the infected cell are related to the age of the infected cell. In 2004, Nelson et al.^[35] developed and analyzed an age structure model for HIV-1 infection that correlated apoptosis rate and virus particle production in infected cells with infection age. Through numerical simulation, it was concluded that the time at which virus levels peak was not only determined by the initial conditions, it also depended on how quickly virus production reached its maximum. In 2017, considering the dynamic characteristics of HIV host cells, Wang et al.^[36] established a mixed model of age structure in which the rate of virus production and the rate of apoptosis of infected cells were related to the age of infection, and explored the effects of two modes of infection, intercellular and intracellular, on virus generation and transmission. In 2017, Y. Zhao, M. Li and S. Yuan ^[37] proposed an SEIR epidemic model using age grouping to study the impact of age on TB transmission in mainland China from 2005 to 2016. Through the estimation of model parameters and sensitivity analysis, it concluded that different age groups had different impacts on TB and proposes two effective measures that can help achieve the goal of the World Health Organization strategy to end TB. In conclusion, models of viral infection age structure have been extensively used in HIV, tuberculosis, hepatitis B, and other studies ^{[38,39,40,41]}.

However, to our knowledge, there are few models to consider that oncolytic viruses and effector T cell production rates depend on the age of tumor cell infection. In addition, we take into account the fact that the limited infective capacity of uninfected tumor cells. Therefore, based on the literature^[23], we propose an age-structured model of oncolytic virus therapy with Holling-Ⅱ functional response. It shows the interaction between uninfected tumor cells, infected tumor cells of different ages, oncolytic viruses, and immune response, which is beneficial for us to better understand the underlying phenomena in tumor therapy. Compared with ^[23], considering the limited ability of uninfected tumor cells to infect, and the age of the infected tumor cells has important implications for oncolytic therapy, the dynamics of the system include the existence and uniqueness of solutions, the stability of the infection-free state and the stability of the infection state were studied. The numerical results are more comprehensive. The results show that when $R_{0} < 1$ and the age of infected tumor cells is in the $\left [0,100 \right]$, the lysis rate of infected tumor cells is higher than 100 years of age. At this point, the density of uninfected and infected tumor cells is reduced to tiny, achieving the effect of tumor treatment. In addition, when $R_{0} > 1$ and the age of the infected tumor cells change significantly within $\left [0, 50 \right]$, the density of uninfected tumor cells first decreases and then continues to increase to a high value, and the density of effector T cells and virus particles remains stable. Therefore, in this case, oncolytic therapy needs to be combined with other cancer therapies to achieve effective treatment.

The rest of this paper is distributed below. In Section 2, we came up with the suggested model. In Section 3, the existence and uniqueness of the solution as well as the existence of infection-free and infection-free homeostasis are discussed. In Section 4, the local and global stabilities of both the infection-free and infection-steady states are proved. In addition, the uniform persistence of the system (2.2) is given. See Section 5 for numerical simulations. In Section 6, we provide some brief discussion.

2.   Model development

In 2020, Zachary Abernathy et al.^[23] studied a soluble tumor model in which viruses interact with tumor cells, based on the short-term dynamics of proposed logical growth instead of exponential growth to allow for long-term population dynamics, and simplified the immune response. Based on this model, the response of effector T cells to infected tumor cell populations and the long-term kinetics of oncolytic virus therapy on tumor cell populations were studied.

Here, $r$ represents the growth rate of uninfected tumor cells, $k$ represents the total carrying capacity of tumor cells, $\beta$ represents the rate of uninfected tumor cells becoming infected, $\gamma _{U}$ denotes the rate at which uninfected tumor cells decay through T cells, $\gamma _{I}$ denotes the rate at which infected tumor cells decay via T cells, $\delta _{I}$ denotes the decay rate of infected tumor cells, $c$ represents the rate at which T cells grow through infected tumor cells, $\delta _{E}$ denotes the decay rate of effector T cells, $\delta _{V}$ represents the rate of decay for virions, $\alpha$ represents the number of virions released via infected cell lysis, and $N$ represents the virotherapy dosage.

However, their model (2.1) ignored the limited infective capacity of uninfected tumor cells and differences in viral response between tumor cells of different ages. Therefore, we propose an age-structured model (2.2) of oncolytic virus therapy with a Holling-Ⅱ functional response that can describe the complex mechanisms of action between uninfected tumor cells $U$, infected tumor cells $I$, effector T cells $E$ and viruses $V$:

with initial and boundary conditions

where $U\left (t \right), E\left (t \right), V\left (t \right)$ represent the density of uninfected tumor cells, effector T cells and virions at time $t$, respectively. $I\left (t, a \right)$ represents the density of infected tumor cells with infection age $a$ at time $t$. $r$ represents the rate at which tumor cells grow. $k$ represents the tolerance of the internal environment to the tumor cells, $\beta$ represents the infection rate of the virus to the tumor cells, and $\gamma _{U}$ is the rate at which uninfected cells decay through T cells. Considering the predator-prey relationship between uninfected tumor cells and viruses, the Holling-Ⅱ functional response function relationship was applied. $D$ is a positive constant. $\delta _{I}\left (a \right)$ is the decay rate of infected cells at age $a$, and $\gamma _{I}\left (a \right)$ is the rate of decay for infected tumor cells of age $a$ via T cells. $c(a)$ is the rate of T cell growth via infected tumor cells of age $a$, $\delta _{E}$ is the rate of decay for effectors T cells. Infected tumor cells release virus particles at a rate of $\alpha\delta _{I}\left (a \right)$. The mortality rate of virus particles is $\delta _{V}$. $N$ is the Dirac function representing the dose of virus injected at time $t$. Let

here, $T(t)$ is the total cell density of the tumor, $P(t)$ is the expansion of the effector T cells at time $t$, and $C(t)$ is the density of the virus at time $t$. Let's assume all the parameters in (2.2) are positive.

3.   Major properties of solutions for system (2)

3.1.   Existence and uniqueness of solutions

To explore the existence and uniqueness of the solution of the system (2.2), make the following denote

and

Define operator $\mathbb{P}:\varpi \left (\mathbb{P} \right)\subset \mathbb{X}\rightarrow \mathbb{X}$ by

here, $\varpi \left (\mathbb{P} \right) = \left \{ 0 \right \}\times \mathbb{W}^{1, 1} \left (0, \infty \right)\times \Re \times \Re \times \Re$. If $\lambda \in \mathbb{C}$ and $Re\left (\lambda \right) > -\tau$, $\tau = \underset{0\leq a\leq \infty }{\min} \left \{ \delta _{I}\left (a \right) \right \}$, there is $\lambda \in \omega \left (\mathbb{P} \right)$, $\omega \left (\mathbb{P} \right)$ is the resolvent set of $\mathbb{P}$, it is given by the following formula: Figures

Rewriting the system (2.2) becomes the following Cauchy problem:

and

$F$ is a nonlinear mapping defined on a bounded set from $\overline{\varpi \left (\mathbb{P}\right)}$ to $\mathbb{X}$, with Lipschitz continuity. In light of this fact, we apply the results given in ^[42] and obtain the following theorem.

Theorem 1. For every $\xi\in \mathbb{X}_{0+}$ there is a definite semi-flow $\left \{ G\left (t \right) \right \}_{t\geq 0}$ on $\mathbb{X}_{0+}$ and a unique continuous map $G\in\mathbb{C}([0, \infty], \mathbb{X}_{0+})$, which represents an integrated solution to the abstract Cauchy problem (3.4); that is to say,

and

Furthermore, the semi-flow $\left \{ G\left (t \right) \right \}_{t\geq 0}$ is asymptotically smooth and bounded dissipative.

Let

where $\kappa = \min \left (\frac{r}{k}, \underline{\delta _{I}} \right)$, $\underline{\delta _{I}} = \underset{a\geq 0}{\inf}\delta _{I}\left (a \right)$, $\eta = \underset{a\geq 0}{\sup} \alpha\delta _{I}\left (a \right)$, $\overline{c} = \underset{a\geq 0}{\sup} c\left (a \right)$. We obtain the corresponding theorem as follows.

Theorem 2. $\Omega$ is a set of positive invariants under semi-flow; that is $G\left (t \right)\Omega \subset \Omega$. Moreover, $\Omega$ attracts all positive solutions of (3.4).

Proof. Let $\overline{\delta _{I}} = \underset{a\geq 0}{\sup}\delta _{I}\left (a \right)$, $\overline{\gamma _{I}} = \underset{a\geq 0}{\sup}\gamma _{I}\left (a \right)$, $\theta = \max \left (\frac{r}{k}, \bar{\delta _{I}} +\frac{\bar{\gamma _{I}}\bar{c}r }{\kappa \delta _{E}}\right)$. For any $\left (\left (0, I\right), U, E, V \right)\in \Omega$ and $t\geq 0$.

That's the integral of the (3.10) above

Therefore, $T\left (t \right)\leq \frac{r}{\kappa }$ if $\left (\left (0, I_{0} \right), U_{0}, E_{0}, V_{0} \right)\in \Omega$.

Similarly, we have

In addition,

We have,

Therefore, $G\left (\Omega \right)\subset \Omega$ is the set of positive invariants. In addition, it's clear that $\Omega$ attracts all the positive solution of (3.4) through (3.11), (3.12) and (3.13). The proof is complete.

3.2.   Existence of steady states

Theorem 3. For the system (2.2), there always exists an infection-free steady state $P^{0}$, and there is an infection steady state $P^{\ast }$ if $R_{0} > 1$.

Proof. System (2.2) always exists an infection-free steady state $P^{0} = \left (U^{0}, 0, 0, V^{0}\right)$, where $V^{0} = \frac{N\left (t \right)}{\delta _{V}}$. There may have an infectious steady state $P^{*} = \left (U^{*}, I^{* }\left (a \right), E^{*}, V^{*} \right)$ meeting the following equations:

When $U^{*}\neq 0$, by solving system (3.16), we get

where $\Gamma _{1}\left (a \right): = exp\left (-\int_{0}^{a} \delta _{I}\left (\sigma \right)d\sigma \right), \Gamma _{2}\left (a, f \right): = exp\left (-\int_{0}^{a}\gamma _{I}\left (\sigma \right)fd\sigma \right), f\left (\cdot \right)\in L_{+}^{1}\left (\Re _{+}, \Re \right)$.

Let

They're all positive.

According to the calculation system (3.16) obtained the quadratic function of $U^{*}$, make it as $g\left (U^{*} \right)$:

Because of $\beta H_{3}\left (E^{*} \right) > \delta _{V}$, the coefficient of $g\left (U^{*} \right)$ is positive. So if and only if $g\left (0\right)\leq 0$, there is a unique positive real root for $g\left (U^{*} \right)$.

Therefore, the basic regeneration number of the system (2.2) is defined as

System (2.2) exists with an infection steady state $P^{\ast }$ provided that $R_{0} > 1$. This completes the proof.

4.   Stability of steady state

4.1.   Local stability in infection-free state

To demonstrate the local stability of infection free status. We evaluate the Jacobian matrix at $P^{0}$ and obtain the following theorem.

Theorem 4. If $R_{0} < 1$, the infection-free steady state $P^{0}$ is locally stable.

Proof. Evaluate

The eigenvalues are $r-\frac{2rU^{0}}{k}-\frac{\beta V^{0}}{U^{0}+D}+\frac{\beta U^{0}V^{0}}{\left (U^{0}+D \right)^{2}}, -\delta _{I}, -\delta _{E}, -\delta _{V}$, respectively.

If the eigenvalues are all negative, we have $\lambda _{1} < 0$

So

The equation about $U^{0}$ is obtained from $P^{0}$ of system (2.2).

When $r\delta _{V}D > \beta N$, the equation (4.4) has positive real roots. Therefore, we get

It shows that all eigenvalues are negative, then $P^{0}$ is locally stable. The proof is complete.

4.2.   Global stability for infection-free state

Denote

and $W\left (t \right) = \mho \left (t \right)+F\left (t \right)+S\left (t \right)$. Using the formal integration of the system (2.2), we have

where $l_{U} = e^{-\left (r-\frac{r}{k} \right)a}, l_{E} = e^{-\delta _{E}a}, l_{V} = e^{-\delta _{V}a}$. Then, the second equation of (2.2) is integrated along the characteristic lines, we have

Therefore,

So for any particular $t_{0} > 0$, we define the operator $X\left (\mathfrak{f} \right), T\left (\mathfrak{f} \right), V\left (\mathfrak{f} \right), Z\left (\mathfrak{f} \right)$ for $\mathfrak{f}\left (\cdot \right)\in \mathbb{C}\left (\left [0, t_{0} \right]; \Re _{+}\right)$ as

where

Hence

where

Hence the theorem for the global stability of the infection-free steady state follows.

Theorem 5. If $R_{0} < 1$, the infection-free steady state $P^{0}$ is globally asymptotic stable.

Proof. Define

Under the assumption $R_{0} < 1$, it has $\varepsilon _{0}$ small enough, such that $\Psi _{\varepsilon _{0}}\left (a \right)da < 1$. For the above $\varepsilon _{0}$, there is $T > 0$, such that $U_{0}e^{\delta _{U}t} < \varepsilon _{0}$ for any $t > T$. So, for any $t > T$, we have

where

In the same way that $\underset{x\rightarrow +\infty}{\lim} H_{4}(t) = 0$, we get $\underset{x\rightarrow +\infty}{\lim} H_{5}(t) = 0$.

Let $\bar{\mho}$ is the solution to the renewal integral equation below:

Therefore, we get $\mho\left (t \right)\leqslant \bar{\mho}\left (t \right)$. In addition, $\underset{x\rightarrow +\infty}{\lim} H_{5}(t) = 0$ and $\int_{0}^{t}\Psi _{\varepsilon _{0}}\left (a \right)da < 1$ if $R_{0} < 1$, using the Paly-Wiener theorem from the literature ^[43,44], $\underset{x\rightarrow +\infty}{\lim} \bar{\mho}(t) = 0$ is obtained, which means $\underset{x\rightarrow +\infty}{\lim} {\mho}(t) = 0$. Hence, when $R_{0} < 1$, according to (4.11), $\underset{t\rightarrow +\infty}{\lim} I\left (t, a \right) = 0$ for any $a\in \left (0, +\infty \right)$, by(4.8), (4.9), (4.10), $\underset{t\rightarrow +\infty}{\lim} U\left (t \right) = U^{0}$, $\underset{t\rightarrow +\infty}{\lim} E\left (t \right) = 0$, $\underset{t\rightarrow +\infty}{\lim} V\left (t \right) = V^{0}$. The proof is complete.

4.3.   Uniform persistence

In this part, we will show the uniform persistence of the system (2.2) when $R_{0} > 1$.

Define

and $\partial M_{0} = \Omega \setminus M_{0}$.

Lemma 1. The subsets $M_{0}$ and $\partial M_{0}$ are positive invariants under the semi-flow $\left \{ G\left (t \right) \right \}_{t\geq 0}$; that is, $G\left (t \right)M_{0}\subset M_{0}$ and $G\left (t \right)\partial M_{0}\subset \partial M_{0}$. In addition, for every $\Im \in \partial M_{0}$, $G\left (t \right)\Im \rightarrow \zeta ^{0}$ as $t\rightarrow \infty$, where

is an infection-free steady state of $\left \{ G\left (t \right) \right \}_{t\geq 0}$.

Proof. Denote

For any $\left (\left (0, I \right), U, E, V\right)\in M_{0}$, we have

where $\theta _{1} = \left \{ \delta _{V}, \bar{\delta _{I}}+\frac{\bar{\gamma _{I}\bar{c}}b}{m\delta _{E}} \right \}$, $\mathbb{J}\left (t \right)\geq \mathbb{J}_{1}\left (0\right)e^{-\theta _{1}t} > 0$. So we know that $G\left (t \right)M_{0}\subset M_{0}$.

In addition, for any $\left (\left (0, I\right), U, E, V \right)\in M_{0}$, we have

which follows that $N-\delta _{V}V\left (t \right) = 0$ if $V_{0}\in \partial M_{0}$. So, $G\left (t \right)\partial M_{0}\subset\partial M_{0}$. That is, $V\left (t \right) = \frac{N}{\delta _{V}}$ means that $B\left (t \right) = 0$. As proof of Theorem 5, it follows that for every $\Im \in \partial M_{0}$, $G\left (t \right)\Im \rightarrow e^{0}$ as $t\rightarrow +\infty$. Then using the conclusion in ^[45,46], we get the following theorem.

Theorem 6. If $R_{0} > 1$, the semi-flow $\left \{ G\left (t \right) \right \}_{t\geq 0}$ is uniformly persistent for the pair $\left (\partial M_{0}, M_{0} \right)$; that is, there exists $\varepsilon > 0$ such that $\underset{t\rightarrow +\infty}{\lim}\inf d\left (G\left (t \right)\Im, \partial M_{0}\right) \geq \varepsilon$ for any $\Im \in M_{0}$.

Proof.

Application of Theorem 4 in ^[45], $\left \{ G\left (t \right) \right \}_{t\geq 0}$ is uniformly persisent only in a below case

where

The existence of $\Im \in \mathbb{W}^{s}\left (\left \{ \zeta^{0} \right \} \right)\cap M_{0}$ is assumed in a paradoxical manner, then there exists $t_{0} > 0$, such that $V\left (t_{0} \right)+\int_{0}^{\infty }I\left (t_{0}, a \right)da > 0$. Since $G\left (t \right) M_{0}\subset M_{0}$, we have $V\left (t \right)+\int_{0}^{\infty }I\left (t, a \right)da > 0$ for any $t\geq t_{0}$. Let

Then

Define

From $\Im \in \mathbb{W}^{s}\left (\left \{ \zeta^{0} \right \} \right)$, we get $\underset{x\rightarrow +\infty}{\lim}V\left (t \right) = V^{0}$, $\underset{x\rightarrow +\infty}{\lim}U\left (t \right) = 0$, and $\underset{x\rightarrow +\infty}{\lim}E\left (t \right) = 0$, hence, $\underset{x\rightarrow +\infty}{\lim}\Im \left (t, 0 \right) = H_{3}\left (0 \right)$. When $R_{0} > 1$, taking $\varepsilon _{0} > 0$ and $\varepsilon _{1} > 0$ are small enough such that $0 < \frac{\beta V^{0}\varepsilon _{1}\left (\varepsilon _{0}+D \right)+\beta H_{3}\left (0 \right)\varepsilon _{0}\left (\varepsilon _{1} +D\right)}{\left (\varepsilon _{1}+D \right)\left (\varepsilon _{0} +D\right)\delta _{V}} < R_{0}-1$. For the above $\varepsilon _{0}$ and $\varepsilon _{1}$, there is $t_{1} > 0$ such that $V\left (t \right) > V^{0}-\varepsilon _{0}$ and $\Im \left (t, 0 \right) > H_{3}\left (0 \right)-\varepsilon _{1}$ for all $t > t_{1}$. Furthermore, for $t > t_{1}$, we get

which means that $\mathbb{J}_{1}\left (t\right)$ is a nondecreasing function for $t\geq t_{1}$. Thus, $\mathbb{J}_{1}\left (t\right)\geq \mathbb{J}_{1}\left (t_{2}\right) > 0$ for all $t\geq t_{2}$ with $t\geq t_{2} = \max \left \{ t_{0}, t_{1} \right \}$, which prevents $\left (I\left (t, a \right), V\left (t \right) \right)$ from converging to $\left (0_{L^{1}}, 0 \right)$ at $t\rightarrow +\infty$. There is an inconsistency with $\Im \in \mathbb{W}^{s}\left (\left \{ \zeta^{0} \right \} \right)$. This section demonstrates the stability of the system infected steady state $P^{\ast } = \left (U^{\ast }, I^{\ast }\left (a \right), E^{\ast }, V^{\ast } \right)$ when $R_{0} > 1$. The proof is complete.

4.4.   Local stability of the infected state

In this section, we discuss the local stability with respect to the infected steady state.

Theorem 7. The infection stable state $P^{\ast }$ is locally stable when $R_{0} > 1$.

Proof. Let

and assume

We have the linearized equation for $P^{\ast }$, which is as follows:

Solve the system (4.21), we get

where

and

Substitute $I_{1}(a), V_{1}$ into the fourth equation of (4.21), we get

We can obtain that the left modulus is greater than the right modulo $\vartheta$ of all complex roots of (4.27) with nonnegative real parts. So all the roots of (4.21) have negative real parts. According to the linear stability principle where the evolution equation as system (2.2), the infected steady state $P^{\ast }$ is locally stable if $R_{0} > 1$. The evidence is complete.

4.5.   Global stability of the infection state

The purpose of this section is to discuss the global stability properties of the system infection steady state by constructing Lyapunov functions.

Define a positive function

Denote

and $\mu \left (0 \right) = N_{1}$, ${\mu \left (a \right)}' = -\alpha \left (\sigma \right)\delta _{I}\left (\sigma \right)+\left (\delta _{I}\left (a \right)+\gamma _{I}\left (a \right)E^{\ast } \right)\mu \left (a \right)$, and $\mu \left (a \right)$ is well defined and bounded.

Theorem 8. If $R_{0} > 1$ and $c\left (a \right) = \epsilon \gamma _{I}\left (a \right)\mu \left (a \right)$ for particular constant $\epsilon > 0$, the infected steady state $P^{\ast }$ of the system is globally asymptotically stable.

Proof. Define a function

$q$ must be positive and has a minimum and is unique at $1$, which satisfies $q\left (1 \right) = 0$.

Then, represent a Lyapunov function

where

Take the derivatives of $Q_{U}\left (t \right), Q_{I}\left (t \right), Q_{E}\left (t \right)$, and $Q_{V}\left (t \right)$ with respect to time.

From (3.17) to know

Subsequently,

Equality is true if and only if

Therefore, when $U^{\ast }\neq 0$, obviously to get $I\left (t, a \right) = I^{\ast }, E\left (t \right) = E^{\ast }, V\left (t \right) = V^{\ast }$. Thus, the maximal invariant set held by (4.38) includes only the infected steady state. The proof is complete.

5.   Numerical simulations

In this section, we confirm the theoretical results through numerical simulation. The solution of the system (2.2) tends to be infection-free homeostasis when the basic reproduction number $R_{0} < 1$. When $R_{0} > 1$, the system solution tends to the infected steady state. Numerical simulations also confirm the stability of the steady state.

The parameter values applied in the program, from the parameter values cited in the literature ^[23,47] and corrected some parameter values, are finally shown in Table 1 below. Furthermore, the maximum age of infected tumor cells is set to 100 years, then the expression of the three age-dependent functions is as follows:

In Figure 1, the steady-state stability of infection-free is demonstrated. It shows the changes of uninfected tumor cells, effector T cells, virions with $t$, and infected tumor cells with $t$ and $a$ of cell infection when $D = 1000$, $N = 8.6382\times 10^{11}$. Both the density of uninfected tumor cells (Figure 1(a)) and the density of infected tumor cells (Figure 1(d)) eventually tend to the smallest possible value (close to $0$) under the effect of effector T cells and virions. According to the research in literature ^[27], effector T cells can kill tumor cells and also affect the virus particles in treatment. Figure 1(b) and Figure 1(c) show the density changes of effector T cells and virions in the process of tumor cells being eliminated. It can be seen that when the lytic number of effector T cells decreases, the number of virions continues to increase and tends to plateau, which is consistent with the fact. Figure 1(d) characterizes that when $R_{0} = 0.0108 < 1$, the infection-free steady state of the system is globally stable. It shows that the lysis rate of infected tumor cells in the $\left [0,100 \right]$ is higher than that of those aged after $100$. At this point, the density of uninfected and infected tumor cells is reduced to tiny, achieving the effect of tumor treatment.

Since the value of $R_{0}$ is inversely proportional to the value of $D$, when $D$ decreases, the value of $R_{0}$ increases, even more than $1$. Using the example, we get $R_{0} = 134.7832 > 1$ when $D = 0.01$. Conclusions of Theorem $6$, Theorem $7$, and Theorem $8$ imply that the infection stable state $P^{*}$ is globally asymptotically stable. In Figure 2, when $D = 0.01$, $N = 1.07978\times 10^{11}$, the results demonstrate the global stability of the infected steady state. It shows the changes in uninfected tumor cells, effector T cells, virions with $t$, and infected tumor cells with $t$ and $a$ of cell infection. At this point, the density of uninfected tumor cells (Figure 2(a)) first decreases and then increases to a stable value (between $7\times 10^{8}$ and $8\times 10^{8}$). This indicates that the tumor cells cannot be eliminated, which also causes the density of effector T cells (Figure 2(b)) to decrease after a brief increase and eventually maintain at a low level. In addition, affected by the number of effector T cells, the density of virions (Figure 2(c)) could not reach the number required for oncolytic virus therapy. In this case, oncolytic virus therapy needs to be combined with other cancer treatments to achieve effective treatment, or the injection dose of oncolytic virus needs to be adjusted to achieve a better therapeutic effect. This needs to be further studied.

6.   Discussions

As we all know, cancer is one of the serious diseases threatening the health of people all over the world. In this paper, we consider a model of oncolytic virus therapy with a Holling-Ⅱ functional response function and age-dependent infection to investigate the theoretical significance of oncolytic therapy as a cancer treatment strategy. Firstly, by defining the semi-flow of the system (2.2), the existence and uniqueness of the solution is studied, and the global attraction of the solution is proved. Then we get the expression for the basic reproduction number $R_{0}$. The Jacobian matrix is used to confirm that the infection-free steady state of the system is locally stable. In addition, we obtain that the infection-free steady state of the system is globally stable by updating the integral equation. If $R_{0} < 1$, the infection-free steady state $P^{0}$ is globally stable. The local stability of the infection is established by linearization. Meanwhile, by constructing the Lyapunov function, a sufficient condition for global stability of infection steady state is obtained. If $R_{0} > 1$, the infection steady state $P^{*}$ is globally stable.

Biologically, numerical simulations can both confirm theoretical results and suggest how to change treatment regimens to achieve the most effective cancer treatment. Not only have we discussed the effect of age on infected tumor cells and Holling-Ⅱ response function to oncolytic viruses treating cancer, but our model also presents two cancer treatment states. The first possibility is a disease-free state. When $R_{0} < 1$, the infection-free steady state is globally stable. At this point, the number of tumor cells is extremely low and does not endanger health, but is never eliminated, and the number of effector T cells and virions remains stable to stimulate the body's immune response and prevent the growth and spread of cancer cells. The model predicts that viral therapy could serve as a new adjuvant therapy in preparation for surgery or radiation therapy. The second possibility is the persistence of cancer. Here $R_{0} > 1$, meaning that the infection's steady state is globally stable. It means that the tumor is no longer sufficient to treat with oncolytic viruses and will continue to develop or metastasize.

In contrast to previous studies ^[23], the highlights of this paper include the following (i) consideration of an age-structured model of oncolytic therapy involving Holling-Ⅱ functional response; (ii) the age parameter is included in the expression for $R_{0}$, which means whether the system (2.2) has stability is affected by the age factor; (iii) the steady state stability of the system (2.2) is discussed more extensively and its biological significance is explained by numerical simulations.

Our model shows that the infection age of tumor cells and the saturation of virions are key parameters for oncolytic therapy. It remains challenging to study age-structure models of oncolytic virus therapy for cancer if continuous changes in viral injection dose over time are taken into account.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (NSFC: 11271068, 11961024), the Joint Training Base Construction Project for Graduate Students in Chongqing (JDLHPYJD2021016), Group Building Scientific Innovation Project for universities in Chongqing (CXQT21021).

Conflict of interest

The authors declare there is no conflict of interest.