Article Contents
Article Contents

On a mathematical model of tumor-immune system interactions with an oncolytic virus therapy

• * Corresponding author: Sophia R-J Jang
• We investigate a mathematical model of tumor–immune system interactions with oncolytic virus therapy (OVT). Susceptible tumor cells may become infected by viruses that are engineered specifically to kill cancer cells but not healthy cells. Once the infected cancer cells are destroyed by oncolysis, they release new infectious virus particles to help kill surrounding tumor cells. The immune system constructed includes innate and adaptive immunities while the adaptive immunity is further separated into anti-viral or anti-tumor immune cells. The model is first analyzed by studying boundary equilibria and their stability. Numerical bifurcation analysis is performed to investigate the outcomes of the oncolytic virus therapy. The model has a unique tumor remission equilibrium, which is unlikely to be stable based on the parameter values given in the literature. Multiple stable positive equilibria with tumor sizes close to the carrying capacity coexist in the system if the tumor is less antigenic. However, as the viral infection rate increases, the OVT becomes more effective in the sense that the tumor can be dormant for a longer period of time even when the tumor is weakly antigenic.

Mathematics Subject Classification: Primary: 92D25, 92B05.

 Citation:

• Figure 1.  A schematic diagram of model (1), depicting the interactions between different populations is presented

Figure 3.  One solution of the immune $ZY_TY_V$ subsystem (8) is plotted to illustrate stabilization of the model, Proposition 3.2(b). The parameter values are given in text

Figure 2.  The parameter values are $h_T = 2.7\times 10^4, \ \delta_{YT} = 3.75\times 10^{-4}, \ r_T = 0.0192, \ a_{AT} = 0.0016$ and $k_{TA} = 1/24$ taken from the baseline values in [38] to illustrate oscillations of the $T_SY_T$ subsystem (5), where the time unit is an hour as given in [38]

Figure 4.  The plots demonstrate the existence of positive equilibria for system (13), which corresponds to the positive intersections of equations (11) and (12). Dimensionless parameter values in (a) are $h_T = 0.5, \ \delta_{YT} = 1, \ \hat Z = 0.5, \ a_{TZ} = 0.3$ and $a_{AT} = 1.6$ with $T_{S_c} = 0.833$ so that $g(0)<f(0)$ and $\delta_{YT}<a_{AT}$. The black vertical line segment is the line $T_S = T_{S_c}$. In (b), $h_T = 0.1$, $\delta_{YT} = 0.5, \ \hat Z = 0.6, \ a_{TZ} = 0.1$ and $a_{AT} = 0.15$ so that $g(0)> f(0)$. The curves have two positive intersections

Figure 5.  (a) Bifurcation diagram using $\beta$ as the bifurcation parameter, where all other parameter values are the same as in (23). The dashed lines at $\log(T_S) = -1$ (blue), $\log(T_S) = 8.7124$ (green), and $\log(T_S) = 3.9172$ (cyan) represent the unstable equilibria $E_0$, $E_1$, and $E_2$, respectively. The black dashed line between $E_2$ and $E_0$ represents an unstable positive equilibrium. (b) and (c) are the time series of cell populations for $\beta = 6\times 10^{-8}$. The other parameter values in (b) are same as those in (a) while $a_{TZ} = 2.4$, $s_{ZR} = 4.8$, $\delta_{ZR} = 1.658$, and $a_{ZZ} = 4.8$ in (c). The tick label -1 on the vertical axis represents a population level less than 0.1

Figure 6.  (a) Bifurcation diagram using $\beta$ and $a_{AT}$ as the bifurcation parameters, where all other parameter values are the same as in (23). (b) A closer look at the bifurcation curves for $\log(\beta)\in [-10.12,-9.97]$

Figure 7.  (a) Bifurcation diagram using $a_{AT} = 0.005$ and $\beta$ as the bifurcation parameters, where all other parameter values are the same as in (23). The blue dashed curve represents unstable positive equilibria and the red curve represents stable equilibria. The light blue and green curves represent stable limit cycles. (b) A closer look at the rectangle shown in (a)

Figure 8.  Time series of cell populations for $a_{AT} = 0.005$ and (a) $\beta = 9.4\times 10^{-9}$, (b) $\beta = 9.4\times 10^{-9}$, and (c) $\beta = 9.5\times 10^{-9}$. The tick label -1 on the vertical axis represents a population level less than 0.1

Figure 9.  (a) Bifurcation diagram using $a_{AT} = 0.008$, $a_{AI} = 0.1$, and $\beta$ and $a_{AI}$ as the bifurcation parameters, where all other parameter values are the same as in (23). (b) Bifurcation diagram using $a_{AT} = 0.008$, $a_{AI} = 0.1$, $b_{T} = 50$, and $\beta$ as the bifurcation parameter. The blue dashed curve represents unstable positive equilibria and the red curve represents stable equilibria. The green curve represents stable limit cycles. (c) Bifurcation diagram using $\beta$ and $r_T$ as the bifurcation parameters, where all other parameter values are the same as in (23)

Table 2.  Equilibria of $T_SZY_TY_V$ subsystem (10) and their biological interpretations

 Equilibrium Interpretation $E_{30}=(0,0,0,0)$ Extinction of all cell populations $E_{31}=(1, 0,0,0)$ Susceptible tumor only $E_{32}=(\bar T_S, 0, \bar Y_S,0)$ Coexistence of susceptible tumor and anti-tumor immune cells $E_{33}=(0,\hat Z,\hat Y_T,\hat Y_V)$ Immune cells only $E_{34*}=(\tilde T_S, \hat Z, \tilde T_T, \hat Y_V)$ Coexistence of susceptible tumor and immune cells

Table 3.  Existence and stability of boundary equilibria of system (4)

 Equilibrium Existence Asymptotic stability $E_{0}$ always unstable $E_{1}$ always $s_{ZR}a_{ZZ}<\delta_Z\delta_{ZR}$, $a_{AT}<\delta_{YT}(h_T+1)$, $b_T<\omega$ $E_{2}$ $a_{AT}>\delta_{YT}(h_T+1)$ $s_{ZR}a_{ZZ}<\delta_Z\delta_{ZR}$, $\bar T_S>(1-h_T)/2$, $\omega(\delta_T+\bar Y_T/h_I)>h_T\delta_T\bar T_S$ $E_{3}$ $s_{ZR}a_{ZZ}>\delta_Z\delta_{ZR}$ $\hat Y_T>h_T$ $E_{4*}$ $s_{ZR}a_{ZZ}>\delta_Z\delta_{ZR}$, $h_T\geq 1$, $\hat Y_T\delta_Z\delta_{ZR}$, $h_T<1$, $\hat Y_T\delta_{YT}$

Table 1.  Existence and stability of equilibria of $T_SZY_TY_V$ subsystem (10)

 Equilibrium Existence Asymptotic stability $E_{30}$ Always Unstable $E_{31}$ Always $s_{ZR}a_{ZZ}<\delta_{ZR}\delta_Z$ and $a_{AT}<\delta_{YT}(h_T+1)$ $E_{32}$ $a_{AT}>\delta_{YT}(h_T+1)$ $s_{ZR}a_{ZZ}<\delta_{ZR}\delta_Z$ and $\bar T_S>(1-h_T)/2$ $E_{33}$ $s_{ZR}a_{ZZ}>\delta_{ZR}\delta_Z$ $\hat Y_T>h_T$ (i.e., $a_{TZ}\hat Z>\delta_{YT}h_T$) $E_{34*}$ $s_{ZR}a_{ZZ}>\delta_Z\delta_{ZR}$, $h_T\geq 1$, $\hat Y_T\delta_Z\delta_{ZR}$, $h_T<1$, $\hat Y_T\delta_{YT}$ * This is a case for which the positive equilibrium is unique. System (10) may have more than one positive equilibrium.

Table 4.  Biological interpretations of boundary equilibria of full system (4)

 Equilibrium Biological meaning $E_0=(0,0, 0, 0,0,0)$ Extinction of tumor, virus and immune cells $E_1=(1, 0,0,0,0,0)$ Susceptible tumor only $E_2=(\bar T_S, 0,0,0,\bar Y_T, 0)$ Coexistence of susceptible tumor and anti-tumor immune cells $E_{3}=(0,0,0,\hat Z, \hat Y_T, \hat Y_V)$ Immune cells only $E_{4}=(\tilde T_S, 0,0,\hat Z, \tilde Y_T, \hat Y_V)$ Coexistence of susceptible tumor and immune cells
•  [1] J. Aguirre-Ghiso, Models, mechanisms and clinical evidence for cancer dormancy, Nat. Rev. Cancer, 7 (2007), 834-846.  doi: 10.1038/nrc2256. [2] B. K. Al-Ramadi, M. J. Fernandez-Cabezudo, H. El-Hasasna, S. Al-Salam, S. Attoub, D. Xu and S. Chouaib, Attenuated bacteria as effectors in cancer immunotherapy, N.Y. Acad. Sci, 1138 (2008), 351-357.  doi: 10.1196/annals.1414.036. [3] L. Allen, An Introduction to Mathematical Biology, Prentice-Hall, New Jersey, 2006. [4] K. S. Cheng, S. B. Hsu and S. Lin, Some results on global stability of a predator-prey system, J. Math. Biol., 12 (1981), 115-126.  doi: 10.1007/BF00275207. [5] M.-H. Chou, H.-C. Wei and Y.-T. Lin, Oregonator-based simulation of the Belousov-Zhabotinskii reaction, Int. J. Bifurc. Chaos Appl. Sci. Eng., 17 (2007), 4337-4353.  doi: 10.1142/S0218127407019998. [6] B. S. Choudhury and B. Nasipuri, Efficient virotherapy of cancer in the presence of immune response, Int. J. Dynam. Control, 2 (2014), 314-325.  doi: 10.1007/s40435-013-0035-8. [7] A. L. de Matos, L. S. Lranco and G. McFadden, Oncolytic viruses and the immune system: The dynamic duo, Mol. Ther. Methods Clin. Dev., 17 (2020), 349-358.  doi: 10.1016/j.omtm.2020.01.001. [8] L. G. de Pillis, A. E. Radunskaya and C. L. Wiseman, A validated mathematical model of cell-mediated immune response to tumor growth, Cancer Res., 65 (2005), 7950-7958.  doi: 10.1158/0008-5472.CAN-05-0564. [9] R. Eftimie, J. Dushoff, B. W. Bridle, J. L. Bramson and D. J. D. Earn, Multi-stability and multi-instability phenomena in a mathematical model of tumor-immune-virus interactions, Bull. Math. Biol., 73 (2011), 2932-2961.  doi: 10.1007/s11538-011-9653-5. [10] R. Eftimie and G. Eftimie, Tumour-associated macrophages and oncolytic virotherapies: A mathematical investigation into a complex dynamics, Lett. Biomath., 5 (2018), 6-35.  doi: 10.30707/LiB5.2Eftimiea. [11] H. Fukuhara, Y. Ino and T. Todo, Oncolytic virus therapy: A new era of cancer treatment at dawn, Cancer Sci., 107 (2016), 1373-1379.  doi: 10.1111/cas.13027. [12] H. Fukuhara and T. Todo, Oncolytic herpes simplex virus type 1 and host immune responses, Curr. Cancer Drug Targets, 7 (2007), 149-155.  doi: 10.2174/156800907780058907. [13] T. F. Gajewski, H. Schreiber and Y.-X. Fu, Innate and adaptive immune cells in the tumor microenvironment, Nat. Immunol., 14 (2013), 1014-1022.  doi: 10.1038/ni.2703. [14] S. Gujar, J. G. Pol, Y. Kim, P. W. Lee and G. Kroemer, Antitumor benefits of antiviral immunity: An underappreciated aspect of oncolytic virotherapies, Trends Immunol., 39 (2018), 209-221.  doi: 10.1016/j.it.2017.11.006. [15] Y. Guo, B. Niu and J. Tian, Backward Hopf bifurcation in a mathematical model for oncolytic virotherapy with the infection delay and innate immune effects, J. Biol. Dyn., 13 (2019), 733-748.  doi: 10.1080/17513758.2019.1667443. [16] D. Haddad, Genetically engineered vaccinia viruses as agents for cancer, treatment, imaging, and transgene delivery, Front. Oncol., 7 (2017), 1-12.  doi: 10.3389/fonc.2017.00096. [17] B. Ingalls,  Mathematical Modeling in Systems Biology: An Introduction, The MIT Press, Cambridge, 2013. [18] S. R.-J. Jang and H.-C. Wei, Deterministic predator-prey models with disease in the prey population, J. Biol. Syst., 28 (2020), 751-784.  doi: 10.1142/S0218339020500151. [19] A. L. Jenner, C.-O. Yun, P. S. Kim and A. C. F. Coster, Mathematical modelling of the interaction between cancer cells and an oncolytic virus: Insights into the effects of treatment protocols, Bull. Math. Biol., 80 (2018), 1615-1629.  doi: 10.1007/s11538-018-0424-4. [20] J. M. Jeschke, M. Kopp and R. Tollrian, Consumer-food systems: why type I functional responses are exclusive to filter feeders, Biol. Rev., 79 (2004), 337-349.  doi: 10.1017/S1464793103006286. [21] P.-H. Kim, J.-H. Sohn, J.-W. Choi, Y. Jung, S. W. Kim, S. Haam and C.-O. Yun, Active targeting and safety profile of PEG-modified adenovirus conjugated with herceptin, Biomaterials, 32 (2011), 2314-2326.  doi: 10.1016/j.biomaterials.2010.10.031. [22] Y. Kitajima and K. Miyazaki, The critical impact of HIF-1a on gastric cancer biology, Cancers, 5 (2013), 15-26.  doi: 10.3390/cancers5010015. [23] N. L. Komarova and D. Wodarz, ODE models for oncolytic virus dynamics, J. Theor. Biol., 263 (2010), 530-543.  doi: 10.1016/j.jtbi.2010.01.009. [24] Y. Kuang and H. I. Freedman, Uniqueness of limit cycles in Gause-type models of predator-prey systems, Math. Biosci., 88 (1988), 67-84.  doi: 10.1016/0025-5564(88)90049-1. [25] H. -Z. Li, X. -D. Liu, R. Yan and C. Liu, Hopf bifurcation analysis of a tumor virotherapy model with two time delays, Physica A, 553 (2020), 124266. doi: 10.1016/j. physa. 2020.124266. [26] J. Li, J.-N. Chen, T.-T. Zeng, F. He, S.-P. Chen, S. Ma, J. Bi, X.-F. Zhu and X.-Y. Guan, CD133+ liver cancer stem cells resist interferon-gamma-induced autophagy, BMC Cancer, 16 (2016), 1-11.  doi: 10.1186/s12885-016-2050-6. [27] X. Li and J.-X. Xu, A mathematical prognosis model for pancreatic cancer patients receiving immunotherapy, J. Theor. Biol., 406 (2016), 42-51.  doi: 10.1016/j.jtbi.2016.06.021. [28] Y. Louzoun, C. Xue, G. B. Lesinski and A. Friedman, A mathematical model for pancreatic cancer growth and treatments, J. Theor. Biol., 351 (2014), 74-82.  doi: 10.1016/j.jtbi.2014.02.028. [29] A. Magen, J. Nie, T. Ciucci and et al., Single-cell profiling defines transcriptomic signatures specific to tumor-reactive versus virus-responsive CD4$^+$ T cells, Cell Reports, 29 (2019), 3019-3032.  doi: 10.1016/j.celrep.2019.10.131. [30] K. J. Mahasa, A. Eladdadi, L. de Pillis and R. Ouifki, Oncolytic potency and reduced virus tumorspecificity in oncolytic virotherapy. A mathematical modelling approach, PLoS ONE, 12 (2017), e0184347, 1–25. doi: 10.1371/journal. pone. 0184347. [31] K. J. Mahasa, R. Ouifki, A. Eladdadi and L. de Pillis, Mathematical model of tumor-immune surveillance, J. Theor. Biol., 404 (2016), 312-330.  doi: 10.1016/j.jtbi.2016.06.012. [32] G. Marelli, A. Howells, N. R. Lemoine and Y. Wang, Oncolytic viral therapy and the immune system: A double-edged sword against cancer, Front. Immunol., 9 (2018), 1-9.  doi: 10.3389/fimmu.2018.00866. [33] D. McDonald and O. Levy, Innate immunity, in Clinical Immunology, 5th Edn., (eds R. Rich, T. Fleisher, W. Shearer, H. Schroeder), A. Frew, and C. Weyand (London: Elsevier), 2019, 39–53. [34] K. W. Okamoto, P. Amarasekare and I. T. D. Petty, Modeling oncolytic virotherapy: Is complete tumor-tropism too much of a good thing?, J. Theor. Biol., 358 (2014), 166-178.  doi: 10.1016/j.jtbi.2014.04.030. [35] A. Reynolds, J. Rubina, G. Clermont, J. Day, Y. Vodovotz and G. B. Ermentrout, A reduced mathematical model of the acute inflammatory response: I. Derivation of model and analysis of anti-inflammation, J. Theor. Biol., 242 (2006), 220-236.  doi: 10.1016/j.jtbi.2006.02.016. [36] E. Sadurska, Current views on anthracycline cardiotoxicity in childhood cancer survivors, Pediatr. Cardiol., 36 (2015), 1112-1119.  doi: 10.1007/s00246-015-1176-7. [37] L. A. Segel and M. Slemrod, The quasi-steady-state assumption: A case study in perturbation, SIAM Rev., 31 (1989), 446-477.  doi: 10.1137/1031091. [38] K. M. Storey, S. E. Lawler and T. L. Jackson, Modeling oncolytic viral therapy, immune checkpoint inhibition, and the complex dynamics of innate and adaptive immunity in glioblastoma treatment, Front. Physiol., 11 (2020), 1-18.  doi: 10.3389/fphys.2020.00151. [39] Y. Tao and Q. Guo, The competitive dynamics between tumor cells, a replication-competent virus and an immune response, J. Math. Biol., 51 (2005), 37-74.  doi: 10.1007/s00285-004-0310-6. [40] H. R. Thieme, Convergence results and a Poincare-Bendixson trichotomy for asymptotically autonomous differential equations, J. Math. Biol., 30 (1992), 755-763.  doi: 10.1007/BF00173267. [41] A. Timalsina, J. P. Tian and J. Wang, Mathematical and computational modeling for tumor virotherapy with meated immunity, Bull. Math. Biol., 79 (2017), 1736-1758.  doi: 10.1007/s11538-017-0304-3. [42] H.-C. Wei, Numerical revisit to a class of one-predator, two-prey models, Int. J. Bifurc. Chaos Appl. Sci. Eng., 20 (2010), 2521-2536.  doi: 10.1142/S0218127410027143. [43] H.-C. Wei, A modified numerical method for bifurcations of fixed points of ODE systems with periodically pulsed inputs, Appl. Math. Comput., 236 (2014), 373-383.  doi: 10.1016/j.amc.2014.03.054. [44] H.-C. Wei, A mathematical model of intraguild predation with prey switching, Math. Comput. Simul., 165 (2019), 107-118.  doi: 10.1016/j.matcom.2019.03.004. [45] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer, New York, 2003. [46] D. Wodarz, Viruses as antitumor weapons, Cancer Res., 61 (2001), 3501-3507. [47] D. Wodarz and N. Komarova, Towards predictive computational models of oncolytic virus therapy: Basis for experimental validation and model selection, PLoS ONE, 4 (2009), e4217. doi: 10.1371/journal. pone. 0004271. [48] K. H. Wong, A. Lu, X. Chen and Z. Yang, Natural ingredient-based polymeric nanoparticles for cancer treatment, Molecules, 25 (2020), 3620. doi: 10.3390/molecules25163620. [49] J. T. Wu, H. M. Byrne, D. H. Kirn and L. M. Wein, Modeling and analysis of a virus that replicates selectively in tumor cells, Bull. Math. Biol., 63 (2001), 731-768.  doi: 10.1006/bulm.2001.0245.

Figures(9)

Tables(4)