Analysis of a delayed and diffusive oncolytic M1 virotherapy model with immune response

https://doi.org/10.1016/j.nonrwa.2020.103116Get rights and content

Abstract

Oncolytic virotherapy (OVT) is a promising therapeutic approach that uses replication-competent viruses to target and kill tumor cells. Alphavirus M1 is a selective oncolytic virus which showed high efficacy against tumor cells. Wang et al. (2016) studied an ordinary differential equation (ODE) model to verify the potent efficacy of M1 virus. Our purpose is to extend their model to include the effect of time delays and anti-tumor immune response. Also, we assume that all elements of the extended model undergo diffusion in a bounded domain. We study the existence, non-negativity and boundedness of solutions in order to verify the well-posedness of the model. We calculate all possible equilibrium points and determine the threshold conditions required for their existence and stability. These points reflect three different fates for OVT: partial success, complete success, or complete failure. We prove the global asymptotic stability of all equilibrium points by constructing suitable Lyapunov functionals, and verify the corresponding instability conditions. We conduct some numerical simulations to confirm the analytical results and show the crucial role of time delays and immune response in the success of OVT.

Introduction

Cancer is a group of diseases caused by uncontrolled growth of abnormal cells. It can spread to other parts of the body and lead to death. In 2018, the total number of new cancer cases were estimated to be 18.1 million and the number of deaths 9.6 million [1]. The number of new cases is expected to rise to 19.3 million by 2025 [2]. The most common conventional therapies used to treat cancer include radiotherapy and chemotherapy [2], [3]. However, these drugs are not selective and can damage normal healthy cells and cause many side effects like hair loss, fatigue and mouth sores [2]. Oncolytic virotherapy (OVT) is an experimental cancer treatment which has occupied a large area in experimental and theoretical research over the last few years [4], [5]. It uses oncolytic viruses which selectively target tumor cells while avoiding normal cells. Oncolytic viruses replicate inside tumor cells, which eventually leads to cell death. Thereupon, the new released virus particles proceed to kill other tumor cells [6]. The ability of oncolytic viruses to reduce tumor size or remove the tumor depends on the efficacy of OVT.

Many oncolytic viruses have been used in clinical trials such as reovirus, measles virus, vesicular stomatitis virus, herpes virus and adenovirus [7], [8], and great results have been achieved [5], [8], [9], [10]. However, there are many obstacles to this treatment approach which can affect the ability of oncolytic viruses to eradicate or control the tumor [5], [6], [8], [11], [12]. One of the main obstacles is caused by the immune system [6], [7], [11], [12]. The anti-tumor immune response can destroy the tumor cell before the viruses get the opportunity to replicate. This may influence the size of virus population and thus weaken the efficacy of treatment. Therefore, oncolytic viruses should be engineered for rapid removal of tumor cells before the accumulation of immune responses [6], [8], [11]. Although the immune system can have a negative impact on OVT, it can enhance and facilitate the success of the therapy [3], [7], [11], [13]. In fact, the relation between OVT and the immune system is poorly understood and it is an active area of research. Another important factor that can affect the success of OVT is the time period of some processes during viral infection [5], [14], [15].

Mathematical models have been considered as an important tool to understand the complicated dynamics of OVT and to determine the best treatment strategies. Many of these models are similar to HIV and HBV viral infection models [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29]. Some of these models were formulated using ordinary differential equations (ODEs). For example, Ratajczyk et al. [30] studied a model for glioma treatment that combines OVT with TNF-α (tumor necrosis factor alpha) inhibitors. TNF-α, which is produced by macrophages, induces apoptosis of infected tumor cells and thus constrains the replication of oncolytic viruses. The inhibition of TNF-α can increase the efficiency of OVT. In [12], the same model was analyzed as an optimal control problem. Malinzi et al. [4] proposed a model to study the effect of combining oncolytic virotherapy with chemotherapy. The model describes the interaction between tumor cells, OVT, chemotherapy and immune response. They showed that OVT can enhance the ability of chemotherapy to eradicate the tumor if the correct dosage combination is used. Okamoto et al. [11] developed a model to show that reducing the selectivity of oncolytic viruses can facilitate tumor eradication before the start of adaptive immune response. Viral therapy models were also formulated using delay differential equations (DDEs). For instance, Wang et al. [5] extended the basic OVT model [31] and studied a system of nonlinear DDEs. They demonstrated that the virus burst size, which represents the number of new viruses formed inside the tumor cell, and the period of viral lytic cycle are two important factors that can influence the outcome of OVT. Hence, these two factors should be taken into account when developing new viruses for virotherapy. Kim et al. [14] considered a delayed model with cytotoxic T lymphocyte (CTL) immune response to study the impact of a time delay on OVT. They took the time delay as a bifurcation parameter and determined the existence conditions of Hopf bifurcation. Ashyani et al. [15] argued that cancer virotherapy fails in the second injection of the virus since CTLs suppress the infection before the virus removes the tumor. However, they showed using a delayed model that there is a narrow interval for the time delay where the therapy can be successful. Remarkably, the aforementioned models ignored the spatial variations in the distribution of cells and viruses. Nouni et al. [32] proposed a delayed model with saturated infection rate. They discussed the existence and local stability of all possible equilibria.

The diffusion of oncolytic viruses and cells within the tumor may play a crucial role in the delivery and success of cancer therapy [6], [13]. For this purpose, many mathematical models in the form of partial differential equations (PDEs) have been developed. For example, Tao and Guo [6] assumed that the immune response and viruses undergo diffusion and studied the spatiotemporal dynamics of uninfected tumor cells, infected tumor cells and immune response with a time delay. They found that the immune response may restrain the replication of oncolytic viruses and thus weaken the efficacy of OVT. Malinzi et al. [2] analyzed models that study the effect of the interaction between oncolytic viruses and CTL immune response on the outcome of OVT. In addition, they performed a stability analysis and obtained traveling wave solutions. Wang et al. [10] introduced a reaction–diffusion model with a time delay. They suggested different treatment strategies for OVT. In each strategy, they computed the optimal dosage required to completely eradicate the tumor depending on gene mutations in the cell.

In 2014, Lin et al. [33] identified a naturally occurring alphavirus M1 as a selective oncolytic virus which targets cancer cells that lack Zinc-finger antiviral protein (ZAP). The virus showed high efficacy in killing tumor cells without hurting normal cells. To get a better understanding of the role of M1 virus, Wang et al. [34] proposed an ODE model to study the effect of using oncolytic M1 virotherapy on the growth of normal and tumor cells. The normal and tumor cells compete on a limited nutrient source. The model takes the form dH(t)dt=κdH(t)β1H(t)N(t)β2H(t)Y(t),dN(t)dt=α1β1H(t)N(t)(d+η1)N(t),dY(t)dt=α2β2H(t)Y(t)β3Y(t)V(t)(d+η2)Y(t),dV(t)dt=μ+α3β3Y(t)V(t)(d+η3)V(t),where H(t), N(t), Y(t) and V(t) denote the concentrations of nutrient, normal cells, tumor cells and free M1 virus particles, respectively. The parameter κ represents the nutrient recruitment rate, while μ is the minimum effective dosage of oncolytic M1 virus. The nutrient uptake rates of normal cells and tumor cells are given by β1HN and β2HY, respectively. The growth rate of normal cells as a consequence of consuming nutrient is given by α1β1HN, while the growth rate of tumor cells is given by α2β2HY. M1 virus kills tumor cells at rate β3YV, while it is generated by tumor cells at rate α3β3YV. The parameter d is the washout rate constant of nutrient and bacteria. The parameters η1, η2 and η3 are the natural death rate constants of normal cells, tumor cells and M1 virus, respectively.

Our purpose in this paper is to extend model (1) by including three important effects that were not considered in [34]. First, the effect of anti-tumor immune response on the effectiveness of the therapy. Second, the effect of time lags during the process of consuming nutrient on tumor eradication. Third, the effect of diffusion on the spatio-temporal distributions of M1 virus and cells. The paper is organized as follows. In Section 2, we formulate the modified model. In Section 3, we prove the well-posedness of the model including the existence, non-negativity and boundedness of solutions. Also, we study all possible equilibrium points and determine the corresponding existence conditions. The last goal was not accomplished in [34] because they concentrated only on a tumor-free equilibrium point. In Section 4, we show that all equilibrium points are globally asymptotically stable and we confirm the local instability conditions. In Section 5, we perform some numerical simulations to support the theoretical results of the previous sections. Finally, Section 6 states the conclusion.

Section snippets

A delayed reaction–diffusion oncolytic M1 virotherapy model

In this section, we extend model (1) by adding two discrete time delays, diffusion and CTL immune response. All model’s components are assumed to diffuse freely in a connected and bounded domain Ω with a smooth boundary Ω. Thus, we study the following delayed partial differential equation model: H(x,t)t=DHΔH(x,t)+κdH(x,t)β1H(x,t)N(x,t)β2H(x,t)Y(x,t),N(x,t)t=DNΔN(x,t)+α1β1ea1τ1H(x,tτ1)N(x,tτ1)(d+η1)N(x,t),Y(x,t)t=DYΔY(x,t)+α2β2ea2τ2H(x,tτ2)Y(x,tτ2)β3Y(x,t)V(x,t)β4Y(x,t)Z(x,t)(

Basic properties

This section shows the existence, non-negativity and boundedness of solutions of model (2)–(4). Also, it determines all possible equilibrium points and the threshold conditions required for their existence.

Theorem 1

Assume that DH=DN=DY=DV=DZ=D. Then, for any given initial conditions satisfying (3) there exists a unique non-negative and bounded solution defined on Ω̄×[0,+).

Proof

Let X=CΩ̄,R5 be the Banach space of continuous functions from Ω̄ to R5 with the norm ψX=supxΩ̄|ψ(x)|, where || is the Euclidean

Global properties

In this section, we show the global asymptotic stability of all equilibrium points mentioned in Theorem 2 by constructing suitable Lyapunov functionals. The construction of these functionals follows the method used in [39], [40], [41], where we first assume a Lyapunov functional with general coefficients and then solve a system of equations to find the values of the coefficients. Moreover, we verify the instability conditions of the equilibrium points.

Theorem 3

The competition-free equilibrium E0 is

Numerical simulations

In this section, we conduct some numerical simulations in order to enhance the theoretical results presented in the previous sections. We take the spatial domain as Ω=[0,2] with space step size Δx=0.02 and time step size Δt=0.1. The parameters β1, β2, β3, β4, η1, η2, η3, η4 and α2 of model (2) are taken as free parameters. The rest of the parameters are fixed in Table 1. To verify the global stability of the six equilibria of model (2), we have the following cases:

  • (i)

    We take α2=0.8, β1=0.03, β2=0.

Conclusion

In this paper, we analyzed a delayed and diffusive oncolytic virotherapy model with immune response. We found that the model has six possible equilibrium points which are defined and stable under the following conditions:

  • (i)

    The competition-free equilibrium E0 is defined and globally asymptotically stable if R11 and R0R˜p. These threshold conditions determine when both the normal and tumor cell populations become extinct. Thus, this case does not show the effect of OVT or immune response on the

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors are grateful to the handling editor and unknown referees for many constructive suggestions, which helped to improve the presentation of the paper.

Funding

This project was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, Saudi Arabia , under Grant No. (D-255-130-1440). The authors, therefore, gratefully acknowledge the DSR technical and financial support.

References (42)

  • ZhangS. et al.

    Dynamic analysis and optimal control for a model of hepatitis C with treatment

    Commun. Nonlinear Sci. Numer. Simul.

    (2017)
  • PanS. et al.

    Threshold dynamics of HCV model with cell-to-cell transmission and a non-cytolytic cure in the presence of humoral immunity

    Commun. Nonlinear Sci. Numer. Simul.

    (2018)
  • XuR.

    Global dynamics of an HIV-1 infection model with distributed intracellular delays

    Comput. Math. Appl.

    (2011)
  • HattafK. et al.

    Global stability for reaction–diffusion equations in biology

    Comput. Math. Appl.

    (2013)
  • WangZ. et al.

    A mathematical model verifying potent oncolytic efficacy of M1 virus

    Math. Biosci.

    (2016)
  • International Agency for Research on Cancer

    The GLOBOCAN 2018 Database

  • MarelliG. et al.

    Oncolytic viral therapy and the immune system: a double-edged sword against cancer

    Front. Immunol.

    (2018)
  • MalinziJ. et al.

    Enhancement of chemotherapy using oncolytic virotherapy: mathematical and optimal control analysis

    Math. Biosci. Eng.

    (2018)
  • TaoY. et al.

    The competitive dynamics between tumor cells, a replication-competent virus and an immune response

    J. Math. Biol.

    (2005)
  • WeinL.M. et al.

    Validation and analysis of a mathematical model of a replication-competent oncolytic virus for cancer treatment: implications for virus design and delivery

    Cancer Res.

    (2003)
  • SuY. et al.

    Optimal control model of tumor treatment with oncolytic virus and MEK inhibitor

    Biomed. Res. Int.

    (2016)
  • Cited by (44)

    • Global stability of a diffusive HTLV-I infection model with mitosis and CTL immune response

      2023, Advances in Epidemiological Modeling and Control of Viruses
    • Stability of a secondary dengue viral infection model with multi-target cells

      2022, Alexandria Engineering Journal
      Citation Excerpt :

      The results of Theorems 3–6 were illustrated using numerical simulations. Our proposed model can be developed in several directions to take into account different biological effects such as time delay [37–41], reaction–diffusion [42–44]. For more generalizations of our work, we mention the readers to the recent papers [45,46].

    • A study on fractional tumour–immune–vitamins model for intervention of vitamins

      2022, Results in Physics
      Citation Excerpt :

      Mayer et al. [10] suggested a simple model with ODEs to illustrate the immune system’s reaction when pathogens attack the body, like tumour cells or viruses. Several mathematical models have been suggested for cancer study by researchers with ODEs, PDEs and delay differential equations [11–17] to interrogate the outcome of tumour evolution on other cells dynamics. Several mathematical models have been established to determine the primary causes of cancer, like hormones, obesity, and diet, but these components have remained unknown mostly.

    • A general non-local delay model on oncolytic virus therapy

      2022, Applied Mathematical Modelling
      Citation Excerpt :

      Furthermore, we can expand our model to include spatial terms. Many previous works [5,6,12,41] provided reasonable methods to understand how the diffusion of oncolytic viruses affects virus therapy. On the other hand, the stochastic term also can be considered in our model.

    • Analysis of a within-host HIV/HTLV-I co-infection model with immunity

      2021, Virus Research
      Citation Excerpt :

      We studied the effect of HIV infection on HTLV-I mono-infection dynamics and vice versa. Our proposed HIV/HTLV-I co-infection model can be generalized and extended to incorporate different biological phenomena such as reaction-diffusion (Bellomo et al., 2019; Elaiw and AlAgha, 2020a,b) and stochastic interactions (Gibelli et al., 2017). Ahmed Elaiw: Conceptualization, Methodology, Software.

    View all citing articles on Scopus
    View full text