Bone remodeling platforms: Understanding the need for multicellular lab-on-a-chip systems and predictive agent-based models

1.   Introduction

Malignant tumors, recognized as one of the most lethal diseases globally, arise from uncontrolled cell growth ^[1,2]. Traditional treatment approaches encompass surgery, chemotherapy, and radiotherapy. An emerging strategy in the field is cancer immunotherapy, a therapeutic approach that entails stimulating and leveraging the immune system to combat the cancer ^[3].

Certain viruses possess the remarkable ability to precisely select and invade cancer cells while leaving normal cells unharmed. These specialized viruses are referred to as oncolytic viruses (OVs) and are pivotal in oncoviral therapy (OVT). OVT has demonstrated significant anti-tumor effects in diverse mouse models and clinical studies ^[3]. Apart from their oncolytic potential to direct towards cancer cells and destroy them, OVs can be engineered to incorporate specific genes or cytokines, enhancing their ability to activate T cell responses ^[3].

Cancer therapies encounter several formidable challenges, with drug resistance emerging as a major concern. Drug resistance is a frequent outcome in cancer treatment and is categorized as acquired or primary (inherent) resistance. Inherent resistance is intrinsic, whereas acquired resistance develops after exposure to drugs. In ^[4], Noll et al. presented significant confirmation of inherent resistance in the measles vaccine virus (MeV) against NCI-60 solid cancer cells. The authors concluded that inherent resistance in oncolytic viral therapy may be linked with epigenetic and genetic alterations following malignant transition.

While several computational and mathematical models have been devised to investigate chemotherapy resistance ^{[5,6,7,8,9,10]}, there has been relatively less exploration of models focused on other therapeutic resistances. Notably, studies such as ^[10,11] delve into acquired/primary resistance in targeted therapy, while ^[12] utilizes spatial probabilistic models to investigate resistance in oncolytic virotherapy. In their investigation, Bhatt et al. ^[12] developed probabilistic, cell-based models to explore the dynamics of resistance in OVT. The model encompasses four distinct populations: normal healthy cells, infected cells, infection-sensitive cancer cells, and infection-resistant cancer cells. The authors systematically explored the model's outcomes, analyzing their dependence on various parameter values and assumptions. Key aspects investigated include the impact of the ratio of the death rate of infected cells to the viral spread rate on cancer eradication, the influence of the timing of virotherapy, and the effects of parameters such as the production rate of resistant cancer cells, virus diffusion distance, and the probability of virus infection.

The study ^[12] utilized simulations based on three different types of cell-based models: 2D lattice, 3D lattice, and Voronoi, demonstrating that the results exhibit qualitative similarity with some variations. Furthermore, the authors examined the scenario of allowing normal healthy cells to become infected. They found that while this approach may increase the likelihood of complete cancer elimination, representing a successful outcome, it comes with the risk of elevating the population of resistant cancer cells, potentially leading to therapy failure. The study suggests that improved therapeutic efficacy can be achieved by sensitizing healthy normal cells to infection ^[12].

Despite being relatively unexplored, resistance to OVT significantly limits therapeutic effectiveness. The article ^[13] extensively explores a spectrum of mechanisms contributing to resistance in OVT, shedding light on the multifaceted challenges faced in utilizing viruses as therapeutic agents against cancer. For example, interferon-mediated resistance stands out as a formidable obstacle, as host cells activate interferon pathways to establish an anti-viral state. In addition, epigenetic modifications, encompassing changes in DNA methylation and histone modifications, have been implicated in altering the gene expression landscape, potentially limiting the permissiveness of host cells to oncolytic viruses. The hypoxic microenvironment within tumors presents another layer of complexity, with low oxygen concentrations inhibiting viral replication. Furthermore, virus-entry barriers pose challenges by hindering efficient entry into target cells, while spatiotemporal restrictions to viral spread highlight the intricate dynamics within tumors that may limit the uniform distribution of therapeutic viruses. These identified mechanisms underscore the need for a nuanced understanding to develop strategies that can effectively overcome OVT resistance and enhance the therapeutic potential of oncolytic viral therapies in the realm of cancer treatment.

The viral cycle in oncolytic viral therapy typically consists of several stages, including attachment and entry, uncoating and gene expression, and replication and assembly ^[14]. Once the virus attaches to specific receptors on the surface of cancer cells, it enters the cancer cell, either through the cell membrane or via receptor-mediated endocytosis. Inside the cancer cell, the viral genome is released from its protein coat. The virus then utilizes the cellular machinery to express its genes. The viral genome is replicated, and new viral components are synthesized within the host cell. The components are assembled to form new viral particles. These processes take time to complete ^[14].

To our knowledge, the study ^[15] represents a pioneering use of deterministic models to explore OVT resistance, although it adopts a generic approach without considering specific mechanisms of resistance. In ^[15], the authors developed systems of differential equations to study the effects of OVT resistance, with the time delay associated with the infection process. The research presented conditions for the existence of equilibria and their stability, deriving critical delays for Hopf bifurcations. The Hopf bifurcation marks the point at which a stable equilibrium in the model becomes unstable, leading to the emergence of periodic oscillations. This can correspond to the transition from steady tumor growth to oscillatory tumor dynamics, where tumor size fluctuates over time. Understanding the occurrence of Hopf bifurcation in cancer models can provide insights into how tumors respond to different treatment strategies.

The study ^[15] concluded that if tumor cells resistant to treatment cannot revert naturally to sensitive cells or through other mechanisms, every cancer cell will inevitably develop resistance to the OVT, even when a time delay in the viral infection is present. Additionally, the authors proved that the delay in infection process cannot destabilize those equilibria with no viruses. The numerical investigations conducted by the authors showed that the existence of resistant cancer cells, along with time delay in the virus infection, substantially amplifies the number of tumor cells, particularly during periods of instability in the population interaction, i.e., when the interaction is not at equilibrium.

The immune system, a complex network of organs, cells, and proteins, serves as the body's defense against infections ^[1]. Malignant tumors often express antigens that can trigger an immune response, creating an anti-tumor effect [1, 2]. According to the Tumor Immuno-Surveillance Hypothesis, the immune system can potentially inhibit the growth of small tumors, eliminating them before they become clinically evident [1, 2]. Building upon the work in ^[15], this study introduces effector cells into the tumor-virus interaction. These immune cells, activated by tumor cells ^[16,17,18], have the ability to eliminate both tumor cells and viruses in oncolytic viral therapy. The primary objective is to investigate the impact of immune cells on the dynamics of the tumor-virus interaction. In addition, as the immune activation process in response to the presence of the tumor is not instantaneous, we introduce a time delay in the activation of immune cells. This results in a model formulated as delay differential equations (DDEs). In contrast to ^[15], where delay is incorporated into the viral cycle, this study considers a delay in immune activation, as it takes days or even weeks for the adaptive immune response to become established. We derive critical delays for which Hopf bifurcations occur at certain boundary equilibria, where one or more populations become extinct.

Contrary to the findings in ^[15], the introduction of immune cells in the tumor-virus interaction prevents all tumor cells from becoming resistant when there is no conversion from resistance to sensitivity, provided that the proliferation rate of immune cells exceeds their death rate. However, if the immune cells have a smaller recruitment rate, then eventually all tumor cells will become resistant. In ^[15], the model exhibits a unique positive equilibrium, whereas the model with the addition of immune cells into the interaction may have three positive equilibria. This new phenomenon can lead to a more complex interaction. Additionally, the numerical simulations conducted in this study demonstrate that the inclusion of immune cells significantly reduces the total tumor count in both stable and unstable scenarios compared to the results presented in ^[15]. Furthermore, this study includes global sensitivity analysis and numerical bifurcation analysis. In particular, combined therapies of OVT and CAR T-cell therapy are considered in Section 5.

The subsequent sections unfold as follows: model development is presented in Section 2, followed by the analysis of ordinary differential equations and DDE models in Sections 3 and 4, respectively. Section 5 delves into the numerical investigation, and the final section concludes with a brief discussion. Appendixes A and B discuss the models of special cases, while Appendix C presents the proofs of some analytical results.

2.   Model formulation

Cancer cells are categorized as either infected or uninfected. The uninfected tumor cells are further separated into resistant or sensitive to OVT. Let $T_r$ and $T_s$ denote the numbers of uninfected resistant and sensitive cancer cells, respectively. The variable $T_i$ denotes the class of infected tumor cells, and $V$ represents the compartment of free viruses. The number of immune cells at the tumor site is denoted by $Z$. A non-sensitive cancer cell will produce resistant descendants except if it transforms into sensitivity before birth occurs. We assume that the resistant and sensitive tumor cells have logistic growth with intrinsic growth rates $r_r$ and $r_s$, respectively, and a carrying capacity denoted by $K$. While resistant and sensitive tumor cells may differ in their response to oncolytic viral therapy, they still share the same physical environment within the tumor. This environment provides the necessary resources for cell growth and imposes limitations on cell proliferation. Therefore, both types of tumor cells are subject to the same constraints imposed by factors such as nutrient availability and spatial limitations, and consequently, they have the same carrying capacity.

Several researchers propose that non-sensitive cancer cells exhibit a smaller growth rate in comparison to cancer cells that are not resistant due to the associated costs of resistance ^[5,12]. To incorporate this finding, we will assume $r_r < r_s$ in numerical simulations. The rate of mutation from sensitivity to resistance is represented by the parameter $a$. Because infected cancer cells have a short life duration and do not require substantial materials, their geographical occupancy, proliferation, and nutrient absorption are not considered. Refer to ^[19,20], in which the tumor's carrying capacity does not include infected cancer cells.

Immune cells are capable of killing sensitive, resistant, and infected tumor cells. Regarding the killing of viruses, it typically involves a separate immune response, such as the innate immune response mediated by phagocytes and the adaptive immune response mediated by antibodies and cytotoxic T cells ^[21]. The immune cells considered in the model include both innate and adaptive systems, where the immune cells are not activated by viruses directly but through infected tumor cells. This modeling aspect is consistent with the works ^[16,17,18]. Additionally, ^[22] distinguishes adaptive immune cells into those that are either tumor-specific or virus-specific, where virus-specific adaptive immune cells do kill viruses. Therefore, even the adaptive immune cells are capable of killing viruses according to Murphy et al. ^[21] and the studies by Vithanage et al., and Storey et al. ^[17,18,22]. The killings of sensitive, resistant, and infected tumor cells, as well as viruses, by immune cells are modeled using a simple mass action approach, with maximal rates denoted by $k_s$, $k_r$, $k_i$, and $k_v$, respectively.

The interaction between sensitive cancer cells and viruses obeys the simple law of mass action where $\beta$ denotes the maximum rate. The parameter $b$ denotes an extra mortality rate of infected cancer cells due to infection, which quantifies pathogenicity of the virus. When infected cancer cells die, new viral progeny are discharged. The number of viruses released per infected tumor cell is represented by $q$. The viral clearance rate is represented by the parameter $c$. Every parameter is positive except potentially for $d$ and $a$. Here, $d \geq 0$ denotes the rate of transition from resistance to sensitivity. Parameters $a$ and $d$ are identified as the Darwinian effect in ^[5] when studying acquired resistance in chemotherapy. With respect to resistance in OVT, the transition from resistance to sensitivity could be influenced by the solid tumor's geometric characteristics. After the demise of certain infected cancer cells, specific resistant cancer cells may become susceptible to infection and consequently transition to a sensitive state. Moreover, the drug ruxolitinib is employed to deal with pancreatic ductal adenocarcinoma (PDAC) cells that are not sensitive to oncolytic vesicular stomatitis virus (VSV), improving the susceptibility of resistant pancreatic cancer cells to VSV ^[23]. This treatment may result in an increase in $d$ because of resistance therapy.

The immune cells are recruited by tumor cells using a Michaelis-Menten function $\frac{\alpha(T_s+T_r+T_i)Z}{m+T_s+T_r+T_i}$, where $\alpha$ represents the maximal proliferation rate, and $m$ is the half-saturation constant. Immune cells can be eliminated or depleted through various processes, such as cell death (apoptosis), clearance by other immune cells, or inhibition of immune cell production ^[3]. We assume that this eliminate rate is constant and is denoted by the parameter $\gamma$. To incorporate the time delay between the encounter with tumor cells and the proliferation of immune cells, we introduce a discrete delay parameter, denoted as $\tau$, into the activation term.

From the above discussion, the interaction between various types of tumor cells, viruses, and immune cells can be captured by the following system:

The initial data on $[-\tau, 0]$ are given by

where $C([-\tau, 0], \mathbb{R}_+)$ represents the Banach space of all continuous functions mapping the interval $[-\tau, 0]$ into $\mathbb{R}_+ = \{x\in\mathbb R: x\geq 0\}$. A conceptual diagram of model (2.1) is given in Figure 1.

Let

and denote the right hand side of (2.1) as

In the following, we verify that model (2.1) is well-posed and ensure that the model is biologically feasible. Its proof is relegated to Appendix C.

Theorem 2.1. The solution to the systems (2.1) and (2.2) exists for $t > 0$, and is nonnegative and bounded.

We next proceed to study model (2.1) without considering the delay in the activation of immune cells.

3.   Analysis of ODE models

In this section, we shall investigate model (2.1) without delay, $\tau = 0$. Specifically, we study the ODE model without OVT in Section 3.1. The full system (2.1) with $\tau = 0$ is investigated in Section 3.2. The models of the special cases $d = 0$ and no resistant tumor cells are briefly discussed at the end of this section, and a more detailed discussion is presented in Appendix A.

3.1.   The model with no viral therapy

Suppose $d > 0$. When $\tau = 0$ and there is no viral therapy, $T_i(0) = 0 = V(0)$, $T_i(t) = 0 = V(t)$ for $t > 0$, and (2.1) is reduced to

In addition, $Z(0) = 0$ implies $Z(t) = 0$ for $t > 0$, and (3.1) reduces to the two-dimensional $T_sT_r$ system. The $T_sT_r$ subsystem of (3.1) has been studied in ^[15]. The subsystem has two equilibria: $(0, 0)$ and $(\hat T_s, \hat T_r)$, with

where $(\hat T_s, \hat T_r)$ is globally asymptotically stable as shown in ^[15]. Consequently, (3.1) always has two equilibria $(0, 0, 0)$ and $(\hat T_s, \hat T_r, 0)$. The Jacobian matrix of (3.1) evaluated at $(0, 0, 0)$ and $(\hat T_s, \hat T_r, 0)$ is given respectively by

It is evident that the stability of the equilibrium point $(0, 0, 0)$ hinges on the eigenvalues of the matrix

In the proof of Proposition 3.1 in ^[15], it is established that the eigenvalues of this matrix have positive real parts. Therefore, it follows that the equilibrium point $(0, 0, 0)$ is unstable. From (3.3), the stability of $(\hat T_s, \hat T_r, 0)$ is determined by the eigenvalues of the left upper $2\times 2$ submatrix, which is Eq (3.3) of ^[15], and the sign of

It was shown in [15, Proposition 3,1] that the eigenvalues of the submatrix have negative real parts. As a result, $(\hat T_s, \hat T_r, 0)$ is asymptotically stable if $\frac{\alpha K}{m+K} < \gamma$ and unstable if $\frac{\alpha K}{m+K} > \gamma$. Let

Based on the proof provided in Appendix C, we can conclude that $(\hat T_s, \hat T_r, 0)$ is globally asymptotically stable in $\Gamma$ if $\frac{\alpha K}{m+K} < \gamma$.

Proposition 3.1. Let $\tau = 0$, $d > 0$, and $\frac{\alpha K}{m+K} < \gamma$. The equilibrium $(\hat T_s, \hat T_r, 0)$ is globally asymptotically stable in $\Gamma$ for (3.1).

The fraction $\cfrac{\alpha K}{m+K}$ is the maximal proliferation rate of immune cells induced by cancer cells. If this maximal rate is lower than the rate of immune cell death $\gamma$, the effector cells cannot sustain themselves, and the immune-free equilibrium $(\hat T_s, \hat T_r, 0)$ becomes globally asymptotically stable.

A positive equilibrium is an equilibrium in which each of its components is positive. To study the existence of positive equilibria, we assume $\frac{\alpha K}{m+K} > \gamma$ by Proposition 3.1, since if the inequality is reversed, solutions all converge to $(\hat T_s, \hat T_r, 0)$, and a positive equilibrium cannot exist. It follows that $\alpha > \gamma$. Define

Setting $Z^{\prime} = 0$ and $Z > 0$, we obtain $T_s+T_r = \xi$. Next, $T^{\prime}_s = 0$ and $T^{\prime}_r = 0$, imply respectively,

Setting the above two expressions of $Z$ equal, leads to the following second degree polynomial equation in $T_s$:

where

In addition, $g(\xi) = -k_s\xi^2a < 0$, and $g(0) = {\mathbf{ a}_3} > 0$. Since $T_r = \xi-T_s$ holds for the two components of any positive equilibrium, we therefore require the root $T_s < \xi$ as well. Also, setting $0 = T^{\prime}_s+T^{\prime}_r$, we have

provided $\xi < K$. i.e., $\xi < K$ is a necessary condition for the existence of a positive equilibrium. Notice $\xi < K$ is equivalent to $\frac{\alpha K}{m+K} > \gamma$, under which $(\hat T_s, \hat T_r, 0)$ is unstable.

We separate the discussion of $g(T_s) = 0$ into the cases of ${\mathbf{ a}_1} < 0$, ${\mathbf{ a}_1} > 0$, and ${\mathbf{ a}_1} = 0$. If ${\mathbf{ a}_1} < 0$, then it is evident that (3.1) exhibits an interior equilibrium uniquely since $g(\xi) < 0$. If ${\mathbf{ a}_1} > 0$, then since $g(\xi) < 0$, $g(x) = 0$ has 2 positive roots $x_1$ and $x_2$ with $x_1 < \xi < x_2$. It follows that (3.1) also has a unique positive equilibrium. If ${\mathbf{ a}_1} = 0$, then ${\mathbf{ a}_2} = -\xi(ak_s+dk_r) < 0$, and thus there exists a unique positive equilibrium. We conclude that (3.1) exhibits a unique interior steady state $(\bar T_s, \bar T_r, \bar Z)$ if $\frac{\alpha K}{m+K} > \gamma$.

Let $\cfrac{\alpha K}{m+K} > \gamma$ and let $J = (\bar j_{ij})$ denote the Jacobian matrix of (3.1) at $(\bar T_s, \bar T_r, \bar Z)$, where

Observe that

and the characteristic polynomial is given by

where

Notice

and thus $b_3 > 0$.

By the Routh-Hurwitz condition ^[24], ${Re}(\lambda) < 0$ for all roots of $\bar P_2(\lambda) = 0$ if and only if $b_1 > 0$, $b_3 > 0$, and $b_1b_2 > b_3$. Since $b_1 > 0$ and $b_3 > 0$, $(\bar T_s, \bar T_r, \bar Z)$ is locally asymptotically stable for system (3.1) if the condition

is satisfied. We cannot obtain a simplified expression of (3.11) involving components of the equilibrium, so we leave (3.11) as it is. The above discussion is summarized as follows.

Proposition 3.2. Let $d > 0$ and $\cfrac{\alpha K}{m+K} > \gamma$. Then, (3.1) exhibits a unique interior steady state $(\bar T_s, \bar T_r, \bar Z)$, which is asymptotically stable if (3.11) is valid, and unstable if (3.11) is reversed.

The existence condition of a unique interior equilibrium implies that

and therefore the total tumor load in $(\bar T_s, \bar T_r, \bar Z)$ is smaller than the corresponding tumor burden in $(\hat T_s, \hat T_r, 0)$. The immune cells are able to reduce the total tumor load if their maximal proliferation rate is greater than their death rate.

3.2.   The full model

We assume $d > 0$ and investigate the full system (2.1) without delay by examining the existence and stability of equilibria. The Jacobian matrix has the following entries:

Denote the two equilibria by $E_0 = (0, 0, 0, 0, 0)$ and $E_1 = (\hat T_s, \hat T_r, 0, 0, 0)$, where $\hat{T_s}+\hat{T_r} = K$. At $E_0$, the Jacobian matrix becomes

Eigenvalues of $J(E_0)$ are the eigenvalues of

Therefore, $E_0$ is always unstable. At $E_1$,

The eigenvalues consist of the eigenvalues of the upper left $2\times 2$ and the lower right $3\times 3$ submatrices. It follows that $E_1$ is locally asymptotically stable if

The stability of $E_0$ and $E_1$ is summarized below.

Proposition 3.3. Let $\tau = 0$ and $d > 0$. System (2.1) always has two equilibria $E_0 = (0, 0, 0, 0, 0)$ and $E_1 = (\hat T_s, \hat T_r, 0, 0, 0)$, where $E_0$ is unstable. The equilibrium $E_1$ is asymptotically stable if $c > \beta q\hat T_s$ and $\cfrac{\alpha K}{m+K} < \gamma$, and unstable if either $c < \beta q\hat T_s$ or $\cfrac{\alpha K}{m+K} > \gamma$.

In the following, we show that there will be no viruses and infected tumor cells if the viral death rate $c$ is sufficiently large. Since the proof is similar to that given in ^[17], it is omitted.

Proposition 3.4. If $c > \beta qK$, then solutions of (2.1) with $\tau = 0$ satisfy $\lim_{t \to \infty}{T_i(t)} = 0 = \lim_{t \to \infty}{V(t)}$.

Based on Proposition 3.4, we provide a global asymptotic stability result of $E_1$ under some constraints. Its proof is presented in Appendix C. Let

Theorem 3.1. Let $\tau = 0$ and $d > 0$. If $c > \beta qK$ and $\frac{\alpha K}{m+K} < \gamma$, then $E_1 = (\hat T_s, \hat T_r, 0, 0, 0)$ is globally asymptotically stable in $D$.

The lump parameter $\beta q K$ can be viewed as the maximal virus production rate. If this rate is smaller than the viral clearance rate $c$, and the maximal proliferation rate of immune cells is not large enough, $\frac{\alpha K}{m+K} < \gamma$, then neither the effector cells nor the viruses can persist in the system.

Suppose now the maximal proliferation rate of effector cells exceeds their death rate, $\frac{\alpha K}{m+K} > \gamma$. Then (2.1) has an equilibrium of the form $E_2 = (\bar T_s, \bar T_r, 0, 0, \bar Z)$, where $(\bar T_s, \bar T_r, \bar Z)$ is the positive equilibrium of the virus-free subsystem (3.1). Its stability depends on the eigenvalues of

where $\bar J_{ik}$, $1\leq i, k\leq 2$, and $\bar j_{13}, \bar j_{23}, \bar j_{31}$ are defined in (3.8). By interchanging columns 3 and 5, and the corresponding rows, $J(E_2)$ is similar to

It follows that the characteristic polynomial of $J(E_2)$ is the product of $\bar P_2(\lambda)$ defined in (3.9) and

Since $j_{33} < 0$, $j_{44} < 0$, $j_{34} > 0$ and $j_{43} > 0$, in addition to the condition of (3.11), one would also require $j_{33}j_{44}-j_{34}j_{43} > 0$, i.e.,

for $E_2$ to be locally asymptotically stable.

Proposition 3.5. Let $\tau = 0$ and $\cfrac{\alpha K}{m+K} > \gamma$. Then, (2.1) has the equilibrium $E_2 = (\bar T_s, \bar T_r, 0, 0, \bar Z)$, where $E_2$ is asymptotically stable if (3.11) and (3.16) are satisfied, and $E_2$ is unstable if the inequality in (3.11) or (3.16) is reversed.

Observe that system (2.1) has an equilibrium of the form $E_3 = (T_s^*, T_r^*, T_i^*, V^*, 0)$, provided $c < \beta q\hat T_s$, where $(T_s^*, T_r^*, T_i^*, V^*)$ is the unique positive equilibrium in the tumor-virus subsystem without immune cells by a previous analysis ^[15]. In particular,

The equilibrium $E_3$ has the Jacobian matrix given by

where

The characteristic polynomial of $J(E_3)$ can be expressed by

where $a_i, \ 1\leq i\leq 4$, are given as

By the Routh-Hurwitz criterion ^[24], $E_3$ is asymptotically stable if

We now summarize the stability properties of equilibrium point $E_3$.

Proposition 3.6. Let $\tau = 0$ and $\beta q \hat T_s > c$. Then, (2.1) has an equilibrium of the form $E_3 = (T_s^*, T_r^*, T_i^*, V^*, 0)$, where $E_3$ is asymptotically stable if the inequalities in (3.20) are satisfied, and $E_3$ is unstable if $a_{55} > 0$ or $a_1a_2a_3 < a_3^2+a_4a_1^2$.

The conditions for the existence and local stability of the boundary equilibria of (2.1) with $\tau = 0$ are summarized in Table 1.

In the cases where resistant tumor cells cannot be converted to sensitivity, for example, if there are no drugs available to reverse the resistance with OVT for a certain tumor type, then the parameter $d$ is set to $0$. We can provide a similar analysis as in the previous discussion and obtain some biological conclusions. In particular, according to Theorem B.1, every sensitive tumor cell becomes resistant, and there are no immune cells involved in the interaction if the following two conditions are met: First, the death rate $c$ of viruses must exceed its maximum production rate $\beta q K$. Second, the immune cell activation rate $\frac{\alpha K}{m+K}$ must be smaller than its clearance rate $\gamma$. This result differs from the findings presented in ^[15], where it was shown that every tumor cell becomes resistant if $c > \beta q K$, as there were no immune cells to control the tumor.

For the situation of no resistant tumor cells, it represents the best-case scenario where the tumor cells are not resistant to OVT. From the analysis given in Appendix A.2, we conclude that even in the absence of resistance to OVT, therapy can fail if the virus death rate is high and the immune activation rate is low, as demonstrated in Theorem A.2.

4.   Analysis of DDE models

The main goal of this section is to study the effects of delay $\tau$ on the stability of equilibria of model (2.1). Our discussion parallels the ODE model discussed in Section 3. Specifically, the model of no OVT is studied in Section 4.1, and Section 4.2 investigates the full system (2.1). The models with $d = 0$ and with no resistant tumor cells are briefly summarized at the end of Section 4.2, with a more detailed presentation given in Appendix B. For an analysis of the local stability of equilibria in DDE models, we refer the reader to the works of Hale ^[25], Kuang ^[26], and Smith ^[27].

4.1.   The model of no OVT

Let $d > 0$. The three-dimensional model of delay differential equations without infected cancer cells is given as

The stability of an equilibrium $(\hat X_0, \hat X_0)$ of (4.1) depends on the stability of the linearized system

at the zero equilibrium. Here,

where $\hat f$ refers to the right-hand side of (4.1) with $\hat X$ and $\hat Y$ similarly defined as in (2.3). The characteristic equation for the linearized system is given as

At $\hat X_0 = (0, 0, 0)$, since

does not depend on $\tau$, $(0, 0, 0)$ remains unstable for $\tau > 0$. At $\hat X_0 = (\hat T_s, \hat T_r, 0)$,

where the upper left $2 \times 2$ submatrix is the Jacobian matrix of the two-dimensional $T_sT_r$ ODE system evaluated at $(\hat T_s, \hat T_r)$, which was analyzed in [15, Proposition 3.1] with eigenvalues consisting of negative real parts. Therefore, the stability of $(\hat T_s, \hat T_r, 0)$ depends on the zeros of

Assume $\frac{\alpha K}{m+K} < \gamma$, i.e., $(\hat T_s, \hat T_r, 0)$ is stable asymptotically for the corresponding non-delay system. Consider

It is known that $\lambda = \cfrac{\alpha K}{m+K}-\gamma < 0$ when $\tau = 0$. If $(\hat T_s, \hat T_r, 0)$ alters its stability for some $\tau > 0$, then (4.2) will possess a pair of purely imaginary roots $\pm i\omega$, $\omega > 0$, at some $\tau > 0$. Replacing $\lambda$ by $iw$ in (4.2), we obtain the following two equations:

Then, from the trigonometric identity $\cos^2(\theta)+\sin^2(\theta) = 1$ for $\theta\in\mathbb{R}$, it follows that $\omega > 0$ must satisfy

The quadratic equation has a unique root

by the assumption. We obtain a contradiction and conclude that all the roots $\lambda$ of (4.2) satisfy $Re(\lambda) < 0$. Therefore, $(\hat T_s, \hat T_r, 0)$ is asymptotically stable for all $\tau\geq 0$ if $\frac{\alpha K}{m+K} < \gamma$. The proof of the last part of Proposition 4.1(b) is presented in Appendix C.

Proposition 4.1. The following statements hold for system (4.1).

(a) $(0, 0, 0)$ is always unstable for $\tau \geq 0$.

(b) $(\hat T_s, \hat T_r, 0)$ is asymptotically stable for all $\tau\geq 0$ if $\frac{\alpha K}{m+K} < \gamma$, and it is is unstable for all $\tau\geq 0$ if $\frac{\alpha K}{m+K} > \gamma$.

Let $\frac{\alpha K}{m+K} > \gamma$ and the condition in (3.11) be satisfied. That is, (3.1) has the unique positive equilibrium $(\bar T_s, \bar T_r, \bar Z)$ by Proposition 3.2 and is moreover asymptotically stable for the ODE model (3.1). Since

the stability of $(\bar T_s, \bar T_r, \bar Z)$ depends on the eigenvalues of

where $\bar j_{ik}$, $1\leq i, k\leq 2$, $\bar j_{13}, \bar j_{23}$, and $\bar j_{31}$, defined in (3.8), are the entries of the Jacobian matrix of the ODE model (3.1) evaluated at $(\bar T_s, \bar T_r, \bar Z)$. The determinantal equation for (4.3) can be written as

where

If $(\bar T_s, \bar T_r, \bar Z)$ undergoes a stability change for some $\tau > 0$, then (4.4) will have a pair of purely imaginary eigenvalues $\lambda = \pm i \omega$, $\omega > 0$, at some $\tau > 0$. Replacing $\lambda$ by $i\omega$ in (4.4), we arrive at

By squaring both sides of (4.6) and (4.7), and adding the resulting two equations, we have

In terms of $x = \omega^2$, Eq (4.8) becomes a polynomial equation of degree three:

Further, Eqs (4.6) and (4.7) imply

If (4.9) has no positive roots, then (4.4) has no pure imaginary roots, and $(\bar T_s, \bar T_r, \bar Z)$ remains stable for all $\tau > 0$. Applying the Routh-Hurwitz criterion ^[24], all roots of (4.9) have negative real parts if and only if the following conditions are satisfied:

We conclude the following with respect to the stability of $(\bar T_s, \bar T_r, \bar Z)$.

Proposition 4.2. Let $\cfrac{\alpha K}{m+K} > \gamma$ and suppose (3.11) holds true. Then, $(\bar T_s, \bar T_r, \bar Z)$ is asymptotically stable for $\tau > 0$ if the inequalities in (4.12) are satisfied.

On the other hand, if $F(x) = 0$ has a simple positive root, the following result is valid. The proof is presented in Appendix C.

Theorem 4.1. Let $\frac{\alpha K}{m+K} > \gamma$ and assume $(\bar T_s, \bar T_r, \bar Z)$ is asymptotically stable for (3.1). If (4.9) has at least one positive simple root, then one can find $\tau_0 > 0$ for which $(\bar T_s, \bar T_r, \bar Z)$ is asymptotically stable for $\tau \in [0, \tau_0)$ and ${\rm Sign}(\frac{dRe(\lambda)}{ d\tau}|_{\tau = \tau_0})\neq 0$.

It is expected that a Hopf bifurcation occurs at $\tau = \tau_0$, causing the equilibrium $(\bar T_s, \bar T_r, \bar Z)$ to become unstable.

4.2.   The full model

We assume $d > 0$. To study the effect of delay $\tau$ on the equilibria for the full model (2.1), we examine the eigenvalues of the Jacobian matrix $A+Be^{-\lambda \tau}$, where $A = f_X(X_0, X_0)$, $B = f_Y(X_0, X_0)$, and $(X_0, X_0)$ is an arbitrary equilibrium. At $E_0 = (0, 0, 0, 0, 0)$,

does not depend on $\tau$, and therefore $E_0$ is always unstable for $\tau\geq 0$. At $E_1 = (\hat T_s, \hat T_r, 0, 0, 0)$,

has eigenvalues consisting of the eigenvalues of the upper left $2\times 2$ and the lower right $3\times 3$ submatrices. We assume $\beta q \hat T_s < c$ and $\frac{\alpha K}{m+K} < \gamma$ to ensure that $E_1$ is asymptotically stable for the ODE system (2.1) with $\tau = 0$. It follows that the stability of $E_1$ for the DDE model (2.1) then depends on Eq (4.2). By Proposition 4.1, we conclude that if $\beta q \hat T_s < c$ and $\frac{\alpha K}{m+K} < \gamma$, all of the eigenvalues of $A+Be^{-\lambda \tau}$ have negative real parts, and $E_1 = (\hat T_s, \hat T_r, 0, 0, 0)$ is asymptotically stable for all $\tau \geq 0$. If $\beta q \hat T_s > c$, then $E_1$ is clearly unstable. If $\frac{\alpha K}{m+K} > \gamma$, one can also show that $E_1$ is unstable by the proof of Proposition 3.5(b). The above discussion is summarized as follows.

Proposition 4.3. The following statements are valid for system (2.1).

(a) $E_0 = (0, 0, 0, 0, 0)$ is always unstable for $\tau\geq 0$.

(b) $E_1 = (\hat T_s, \hat T_r, 0, 0, 0)$ is asymptotically stable for $\tau\geq 0$ if $\beta q \hat T_s < c$ and $\frac{\alpha K}{m+K} < \gamma$, whereas $E_1$ is unstable for all $\tau\geq 0$ if $\beta q \hat T_s > c$ or $\frac{\alpha K}{m+K} > \gamma$.

Observe that Proposition 4.3 implies the delay in the immune activation has no effect on the tumor-free equilibrium $E_0$ as well as the immune-free equilibrium $E_1$.

Let $\frac{\alpha K}{m+K} > \gamma$. Then, $E_2 = (\bar T_s, \bar T_r, 0, 0, \bar Z)$ exists, and its stability depends on the eigenvalues of

with $\bar j_{ij}$ given in (3.8). The characteristic equation can be written as

where $p_i, \ 1\leq i\leq 3$, and $q_i, 1\leq i\leq 2$, are defined in (4.5).

Assume $E_2$ is asymptotically stable for the ODE model, i.e., both conditions (3.11) and (3.16) hold. Then the two roots of

have negative real parts. The stability of $E_2$ therefore depends on the roots of (4.4), and we have the following conclusions.

Theorem 4.2. Let $d > 0$ and let $\cfrac{\alpha K}{m+K} > \gamma$. Then, $E_2 = (\bar T_s, \bar T_r, 0, 0, \bar Z)$ exists. Suppose (3.11) and (3.16) are satisfied. In addition, if (4.12) is true, then $E_2$ is asymptotically stable for $\tau\geq 0$.

The delay may not affect the stability of the virus-free equilibrium $E_2$ in the presence of immune cells, provided that the parameters satisfy the conditions outlined in Theorem 4.2. On the other hand, if $F(x) = 0$ has positive root, then the following is true.

Theorem 4.3. Let $d > 0$ and let $\cfrac{\alpha K}{m+K} > \gamma$, and (3.11) and (3.16) be satisfied. Assume $F(x) = 0$ exhibits at least one simple positive root. Then, one can find $\tau_0 > 0$ such that $E_2 = (\bar T_s, \bar T_r, 0, 0, \bar Z)$ is asymptotically stable for $\tau\in[0, \tau_0)$ and $\cfrac{dRe(\lambda)}{d\tau}|_{\tau = \tau_0}\neq 0$.

It follows that the delay $\tau$ in the immune cell proliferation can destabilize the equilibrium $E_2$. As $\tau$ increases beyond the critical value $\tau_0$, the virus-free equilibrium $E_2$ becomes unstable, leading to oscillations in the tumor-virus-immune interaction as a result of the delay in immune recruitment.

Let $c < \beta q \hat T_s$. Then, $E_3 = (T_s^*, T_r^*, T_i^*, V^*, 0)$ exists, and the stability of $E_3$ depends on the corresponding Jacobian matrix, written as

where $a_{ij}$, $1\leq i, j\leq 5$, are defined in (3.18). The characteristic equation is given by

with $a_i$, $1\leq i\leq 4$, defined in (3.18). Assume $a_{55} < 0$ and $a_1a_2a_3 > a_3^2+a_1^2a_4$ such that $E_3$ is asymptotically stable for the ODE model. That is,

has only roots with negative real parts. Let

where $\tau > 0$ is fixed. If $g(\lambda)$ has a pair of pure imaginary roots $\lambda = \pm iw$, $w > 0$, then $w^2 = a_{55}(a_{55}+2\gamma)$, where $a_{55} < 0$ and $a_{55}+2\gamma > 0$. We obtain a contradiction and conclude that all roots of $g(\lambda)$ have negative real parts. On the other hand, if $a_{55} > 0$, then $g(0) = -a_{55} < 0$, $g(\infty) = \infty$, and $g^\prime(\lambda) > 0$ for $\lambda\geq 0$. Thus, $g(\lambda) = 0$ has one positive root. Therefore, $E_3$ is unstable. Further, $E_3$ is clearly unstable if $a_1a_2a_3 < a_3^2+a_1^2a_4$. Below is a summary of the discussion.

Theorem 4.4. Let $\beta q\hat T_s > c$. Then, the equilibrium $E_3 = (T_s^*, T_r^*, T_i^*, V^*, 0)$ exists, and the following statements hold true.

(a) $E_3$ is asymptotically stable for $\tau\geq 0$ if $a_{55} < 0$ and $a_1a_2a_3 > a_3^2+a_1^2a_4$.

(b) $E_3$ is unstable for $\tau\geq 0$ if either $a_{55} > 0$ or $a_1a_2a_3 < a_3^2+a_1^2a_4$.

Comparing Theorem 4.4 with Proposition 3.6, we observe that the delay $\tau$ in immune cell proliferation does not affect the stability of the immune-free tumorous equilibrium $E_3$.

For the scenario of no conversion from resistance to sensitivity, it is concluded from Appendix B that the delay in the immune activation has no effects on the stability of the two equilibria $E_0 = (0, 0, 0, 0, 0)$ and $E_1^0 = (0, K, 0, 0, 0)$. In addition, if $\cfrac{\alpha K}{m+K} > \gamma$, then $E^0_2 = (0, \xi, 0, 0, \tilde{Z^0})$ exists, where $\xi = \cfrac{m\gamma}{\alpha-\gamma}$ and $\tilde{Z^0} = \dfrac{r_r}{k_r}(1-\dfrac{\xi}{K})$. It follows from Appendix B.1 that the delay can destabilize $E_2^0$.

In the best scenario where there is no resistance to OVT, the four-dimensional $T_sT_iVZ$ system is further reduced to the two-dimensional $T_sZ$ system, as shown in (14.5), if no viruses and infected tumor cells are present. System (14.5) possesses an equilibrium of the form $(\xi, \tilde Z)$, provided that $\cfrac{\alpha K}{m+K} > \gamma$, where $\xi$ is defined in (3.5) and $\tilde Z$ is given in (14.6). Therefore, the four-dimensional system has an equilibriun of the form $(\xi, 0, 0, \tilde Z)$. This equilibrium is no longer asymptotically stable for $\tau > 0$. The time delay $\tau$ can destabilize this equilibrium. However, it can be verified that the delay has no destabilization effects on the equilibria $E_0 = (0, 0, 0, 0)$ and $\tilde E_1 = (K, 0, 0, 0)$. In addition, when the parameters satisfy the conditions given in Theorem C.2(b), the tumorous equilibrium $\tilde E_2 = (\xi, 0, 0, \tilde Z)$ in the presence of immune cells becomes unstable, as the delay $\tau$ passes beyond the critical value $\tau_0$. Consequently, the tumor-virus-immune interaction may exhibit oscillations due to the time delay in the immune activation. For the equilibrium $\tilde {E_3} = (\tilde T_s, \tilde T_i, \tilde V, 0)$, since there are no immune cells present in the equilibrium, it is trivial to verify that the delay cannot have any effect on its stability.

5.   Numerical simulation

In this section, we employ numerical tools to study system (2.1), addressing the effects of immune cells compared to the previous work in ^[15] in Section 5.1, incorporating global sensitivity analysis as discussed in Section 5.2, and covering numerical techniques for bifurcation analysis in Section 5.3. Section 5.4 explores the impact of parameters such as the virus transmission rate and half-saturation constant of immune activation on treatment success. Additionally, Section 5.5 considers a combination therapy of OVT and CAR T-cell therapy.

5.1.   Effects of immune cells

In our prior study ^[15], we explored the interactions between tumor cells and oncolytic viruses. This section now focuses on demonstrating the impact of immune cells, in comparison with our previous work. Initially, we compare the outcomes of models that do not incorporate delays in either viral infection or immune cell activation. For this comparison, we set the common parameter values in both the four-dimensional and five-dimensional systems to be identical to those in the parameter set given in [15, Eq (4.1)], as follows:

The initial conditions chosen are $(7\times 10^6, 10^2, 0, 5\times 10^6)$ for the model with no immune cells and $(7\times 10^6, 10^2, 0, 5\times 10^6, 2\times 10^3)$ for the five-dimensional system. In numerous mouse experiments, such as those documented in ^[28,29], researchers typically inoculate mice subcutaneously with either $10^6$, $2\times 10^6$, or $5\times 10^6$ cancer cells. At the time of injection, it is generally assumed that no immune cells are present in the mice. To allow the immune cells to proliferate to a quantity of $2\times 10^3$, we increase the tumor burden to $7\times 10^6$ cells, while maintaining a small number of resistant tumor cells, around $100$ in this case. The doses of oncolytic viruses (OVs) administered vary, ranging from $2\times 10^6$ in ^[29] to $10^8$ pfu in ^[28]. For this part of the numerical simulation, we select $5\times 10^6$, but slight variations in the initial conditions should not impede our biological conclusions.

Figure 2(a) depicts the outcome of the four-dimensional $T_sT_rT_iV$ model before incorporating immune cells, while Figure 2(b)–(d) depicts the outcome of the ODE five-dimensional model (2.1) with $\tau = 0$ for various parameter values. Specifically, the additional parameters for the five-dimensional $T_sT_rT_iVZ$ model are set as follows: (b) $\gamma = 0.03, \ k_s = 1.9\times 10^{-4}, \ k_r = 1.9\times 10^{-4},$ $k_i = 4, \ k_v = 0.025, \ m = 500$ with varied $\alpha = 0.03, \ 0.035, \ 0.04$; (c) $\alpha = 0.035, \ k_s = 1.9\times 10^{-4},$ $k_r = 1.9\times 10^{-4}, \ k_i = 4, \ k_v = 0.025, \ m = 500$ with varied $\gamma = 0.03, \ 0.035, \ 0.04$; and (d) $\alpha = 0.035, \ \gamma = 0.03, \ k_s = 1.9\times 10^{-4},$ $k_r = 1.9\times 10^{-4}, \ k_i = 4, \ k_v = 0.025$ with varied $m = 500, \ 5000, \ 10^{5}$.

In Figure 2(b), the ratio $\cfrac{K}{m+K} = 0.999$, indicating that when $\alpha = 0.03 = \gamma$, the immune system is weak, and the total tumor load is large. As $\alpha$ is increased, the total tumor burden is significantly reduced. This finding is further confirmed by the scenarios shown in plots (c) and (d). Obviously, the solutions of those in Figure 2(b) with $\alpha = 0.03$, and (c) with $\gamma = 0.05$ and $0.1$ converge to an equilibrium level. By utilizing a three-dimensional $T_sT_rZ$ plot, we can determine that the remaining solutions in Figure 2(b) with $\alpha = 0.035$ and $0.4$ converge to a periodic solution, and so does the solution in (c) with $\gamma = 0.03$. Moreover, each solution presented in (d) also converges to a periodic solution.

Through further adjustments to these parameters, we compared the outcomes and consistently observed that the immune-incorporated system maintained a lower total tumor count compared to the four-dimensional model without immune cells, especially when the recruitment rate of immune cells exceeds the death rate.

Next, we compare the results of the corresponding five-dimensional delay model $T_sT_rT_iVZ$ with the four-dimensional $T_sT_rT_iV$ ODE model, which does not include immune cells. Here, we consider the delay in immune cell activation with $\tau = 5$. In previous results without delay, as shown in Figure 2(b), an increase in the maximal immune cell proliferation rate $\alpha$ visibly lowered the total tumor burden. However, as depicted in Figure 3(a), with a small $\alpha$ value of $0.03$ and delay, the stability of the equilibrium remains unchanged. When $\alpha$ is increased to $0.035$, the amplitude of oscillation is larger in the delay model, and the total tumor load may exceed that for $\alpha = 0.03$. A similar phenomenon is observed as $\alpha$ is further increased to $0.04$. However, the tumor burden is clearly smaller than that of the model without immune cells.

Regarding the immune cell death rate $\gamma$, it is apparent from Figure 3(b) that immune cells are not effective in controlling the tumor if the death rate $\gamma = 0.1$ is large. As $\gamma$ is reduced to $0.03$, it is evident that the immune cells can lower the total tumor load even with a recruitment delay as shown in Figure 3(b). Further, the tumor load attains a larger maximum value with a delay compared to that without delay presented in Figure 2(c). A similar conclusion is also obtained as $m$ is varied. The time delay in immune recruitment can increase the amplitude of the tumor burden.

5.2.   Sensitivity analysis

The feasible parameter ranges of the model are compiled from the literature and presented in Table 2. The killing rate $k_r$ of resistant cancer cells by immune cells was estimated over a broader range of $10^{-4}$–$10^{-2}$ cell$^{-1}$ day$^{-1}$, while the killing rate $k_s$ of susceptible tumor cells by immune cells is set at $10^{-5}$–$10^{-3}$ cell$^{-1}$ day$^{-1}$. From the previous subsection, the initial conditions for Eq (2.1) with $\tau = 0$ are set as $(7\times 10^6, 10^2, 0, 5\times 10^6, 2\times 10^3)$. For Eq (2.1) with $\tau > 0$, the history function is set as

and

In OVT, OVs can be administered following single dose or multiple dose regimens ^[35]. In this paper, single dose regimens are used, and OVs are administered at time $t = 0$. Population dynamics are solved numerically using the built-in function ode23tb in Matlab for $\tau = 0$, and a history for population levels is created for solving Eq (2.1) with $\tau > 0$.

Sensitivity analysis is a common method used to identify critical inputs in a model. Global sensitivity analysis (GSA) is a method to analyze the impact of each input or parameter on the output uncertainty when all model inputs or parameters change randomly ^[44]. In this subsection, we conduct GSA utilizing the Partial Rank Correlation Coefficient (PRCC) to assess the relationship between each parameter and tumor size, $T_s+T_r$. Furthermore, GSA aids in selecting parameters for the bifurcation analysis to be performed in the next section. The parameter ranges are shown in Table 2. For the introduction and implementation of this method, we refer to the works by Marino et al. ^[44] and Vithanage et al. ^[17,18].

According to Liu et al. ^[45] and De Matos et al. ^[46], rapid immune-mediated clearance of oncolytic viruses can affect the efficacy of OVT. Following virus administration, oncolytic virus infection and replication within the tumor play a crucial role in the initial days ^[47,48]. Selecting $T = 7$ for a short time interval enables the identification of parameters potentially related to the early-stage anti-tumor or anti-viral effects. In contrast, choosing $T = 21$ for a longer time interval helps pinpoint parameters associated with the asymptotic tumor size and the treatment outcome. GSA will be conducted with two endpoints, $T = 7$ and $T = 21$.

In Figure 4, PRCC is computed for each parameter. The criteria for the strength of a correlation in ^[49] are used in this paper. A correlation is considered non-important if $|\mbox{PRCC}| < 0.2$, weak if $0.2 < |\mbox{PRCC}| < 0.5$, strong if $0.5 < |\mbox{PRCC}| < 0.7$, and very strong if $0.7 < |\mbox{PRCC}|$. Figure 4(a), (b) depicts the PRCC for each parameter in Eq (2.1) with $\tau = 0$. The maximum immune cell proliferation rate, $\alpha$, shows a strong negative correlation with the tumor size at $T = 7$. The virus transmission rate, $\beta$, exhibits a weak positive correlation with the tumor size when $\tau = 0$. All other parameters show unimportant correlations with the final tumor size.

It has been reported that the activation of the adaptive immune response takes at least several days ^[43], and Pulendran et al. and Sun et al. have noted that after infection, T lymphocytes go through an expansion phase, reaching their peak on day 8 in response to antigen stimulation ^[42,50]. To examine the effect of delay in immune cell proliferation, we choose a larger delay in its range (Table 2). In Figure 4(c), (d), we consider Eq (2.1) with $\tau = 7$. The maximum immune cell proliferation rate, $\alpha$, and the killing rate, $k_s$, have strong negative correlations with the tumor size at $T = 7$. All other parameters show weak or unimportant correlations with the final tumor size.

Figure 4 shows consistency of important parameters in all cases. Parameters $\alpha$ and $k_s$ related to immune cell proliferation and immune response have strong correlations with the tumor size. We explore the correlations by studying the tumor dynamics, $T_s(t)+T_r(t)$ for $0\leq t \leq 21$, with selected parameter values for $\alpha$, $k_s$, and $\beta$. All other parameter values are fixed as follows:

Figure 5(a), (d) shows the tumor dynamics $T_s+T_r$ for selected $\alpha$ values. Larger proliferation rates result in smaller tumor sizes. Equation (2.1) models a delay in immune cell proliferation. Figure 5(d) shows a delay in tumor cell elimination compared with Figure 5(a), which reflects a delay in immune cell proliferation. Figure 5(b), (e) shows larger tumor killing rates, $k_s$, resulting in smaller tumor sizes. Interestingly, Figure 5(b), (e) suggests that the killing rate, $k_s$, must exceed $5\times 10^{-4}$ to effectively destroy the tumor cells. Hale et al. have reported that the immune response is important for the success of treatment and preventing tumor recurrence ^[51].

Figure 5(c), (f) shows the same tumor dynamics when $\beta < 10^{-6}$, implying that the virus transmission rate is insufficient to induce an anti-tumoral effect if $\beta < 10^{-6}$. In contrast, the tumor population drops dramatically initially when $\beta = 10^{-4}$ and $\beta = 10^{-3}$, indicating that tumor cells are infected after viruses are administered. As the virus transmission rate increases, more tumor cells are infected, leading to a smaller population level of $T_s+T_r$. The tumor size for $\tau = 7$ is slightly smaller than that for $\tau = 0$, implying that a delay in immune cell proliferation may reduce anti-viral response. Based on the above numerical result, we assume $\beta\geq 10^{-6}$ from this point onward.

5.3.   Numerical methods for bifurcation analysis

Bifurcation analysis is the study of changes in the asymptotic behavior of a model under parameter variation and can provide a complete picture of asymptotic dynamics on a parameter domain. Bifurcation analysis usually begins with the identification of equilibria. Given a set of parameter values, the numerical method constructed in ^[18] can be directly applied to locate the equilibrium $E_2$ and all positive equilibria. The mathematical analysis given in Section 3 can be used to find all other boundary equilibria. The numerical method in ^[52,53] is then applied to compute the bifurcations of equilibria. To study the effect of OVT, the virus transmission rate $\beta$ is used as a bifurcation parameter. Note that Figure 5(c), (f) shows that the virus transmission rate $\beta$ must be greater than $10^{-6}$. We use the range $[10^{-6}, 1]$ for $\beta$. From GSA, the immune response is an important factor associated with the tumor size. The half-saturation constant, $m$, of the immune cell proliferation rate measures the steepness of the curve for the activation of immune cells in response to the presence of tumor cells. A small half saturation constant results in a sharp increase in immune cells with increasing tumor cells. We use $m$ as the other bifurcation parameter.

Mahasa et al. ^[16] have suggested a half saturation constant of 40 cells in their model. The term describing immune cell proliferation in their model is independent of immune cells. With such a small half-saturation constant in Eq (2.1), immune cell proliferation is almost constant in the presence of tumor cells. The process of activation of immune cells involves the interaction between tumor cells and APCs (antigen presenting cells) ^[54]. Once activated, CD8+ T cells undergo clonal expansion with interleukin-2 (IL-2) stimulation ^[55]. Based on the biological mechanism of the activation and expansion process, the immune cell proliferation term in the fifth equation of the ODE system is modeled as a function of both tumor cells and immune cells. To match a similar proliferation rate in the model by Mahasa et al. ^[16], the half-saturation constant $m$ may be a large number. Storey et al. ^[22] has suggested a range of $40-10^5$ cells for $m$ based on the work by Banerjee et al. ^[56]. The model considered in ^[56] does not consider a half saturation constant in the recruitment term for immune cells. Thus, the range of the half-saturation constant, $m$, has not been validated or justified in these works ^[16,22,56].

It has been shown that a small half-saturation constant may cause destabilization and produce sustained oscillations in population sizes ^[57,58,59]. Incorporating a delay in the ODE system can also lead to population oscillations in the solution. Using a range of small half-saturation constants may hinder the study of the destabilization effect caused by the time delay. Furthermore, as pointed out above, the half-saturation constant $m$ may be a large number, and it affects the tumor size in the equilibria. Therefore, we use a wide range for $m$. Let $[10^{-6}, 1]\times [10^5, 10^9]$ be the parameter domain for $\beta$ and $m$. Other parameters are fixed as in Eqs (5.3) and (5.4). An adaptive grid method in Section 2.1 in ^[52] is applied with a 10 by 10 grid in the parameter domain. The regions containing bifurcation curves are refined with 6 levels. The code is implemented in Matlab, where the Newton method is used to find the equilibrium that has no analytic solution, and the built-in function eig is used to solve the eigensystem of the Jacobian matrix. The bifurcation diagram is shown in Figure 6(a).

All six equilibria coexist in a region of the parameter domain. Proposition 3.3 has shown that the tumor-free equilibrium $E_0$ is unstable. From Eqs (3.2), (5.3) and (5.4), along with Proposition 3.3, the equilibrium $E_1$ is unstable in the parameter domain (Figure 6(a)). The numerical simulation agrees with Proposition 3.3. There are one saddle node bifurcation curve $SN$, one Hopf bifurcation curve $Hopf_0$, and two transcritical bifurcation curves $Tr_1$ and $Tr_2$. Both positive equilibria appear through saddle-node bifurcations on the curve $SN$. We use the notation $E_{41}$ for the positive equilibrium that might be stable. The equilibrium $E_{42}$ is never stable. The positive equilibrium $E_{41}$ changes its stability through Hopf bifurcations on the curve $Hopf_0$ when $\tau = 0$. Other bifurcations such as bifurcations of limit cycles may occur in the parameter domain. These bifurcations do not affect the stability property of the equilibria shown in Figure 6(a) and are out of the scope of this paper. Therefore, they will not be studied in this paper. The positive equilibria $E_{41}$ loses its stability through transcritical bifurcations, $Tr_2$, and the equilibrium $E_3$ becomes stable when $m$ increases. Note that the region above $Tr_2$ satisfies the condition in Eq (3.20), and the curve $Tr_2$ satisfies the condition $a_{55} = 0$ in Eq (3.20). The equilibrium $E_2$ is stable when $(\beta, m)$ lies below the transcritical bifurcation curve $Tr_1$ for $\tau = 0$, and it is unstable when $(\beta, m)$ lies above $Tr_1$. Our numerical computation agrees with Proposition 3.5. The conditions Eqs (3.11) and (3.16) are satisfied in the region where $E_2$ is stable, and the curve $Tr_1$ satisfies the left hand side of Eq (3.16) equal to zero. Our numerical computation agrees with Proposition 3.6. Bistability exists in a region below $Tr_1$.

Consider Eq (2.1) with $\tau > 0$. Recall that ODE and DDE systems have the same equilibria. The local stability of an equilibrium depends on the signs of the real parts of the solutions to $det(A-\lambda I +Be^{-\lambda \tau}) = 0$, where $A$ and $B$ are defined in Section 4.2. If $\lambda = 0$ is a solution to $det(A-\lambda I +Be^{-\lambda \tau}) = 0$, it is also a solution to $det(A+B-\lambda I) = 0$. System (2.1) goes through the same bifurcations as the ODE system does on the bifurcation curves $SN$, $Tr_1$, and $Tr_2$. It is known that incorporating a time delay $\tau$ into the ODE system may produce destabilization effects and cause fluctuations in population dynamics through Hopf bifurcations ^[60]. A numerical method is proposed to identify such Hopf bifurcations for system (2.1) in a parameter domain when $\tau$ is fixed.

Based on the bifurcation curves for $\tau = 0$ shown in Figure 6, we focus on the Hopf bifurcations for system (2.1) in the region where $E_{41}$ is stable. At a Hopf bifurcation of $E_{41}$ for $\tau > 0$, $det(A-\lambda I +Be^{-\lambda \tau}) = 0$ has a pair of pure imaginary roots. All other solutions have negative real parts. Starting with any Hopf bifurcation point on $Hopf_0$, say $(\beta_0, m_0)$, let $\pm i\omega_0$ be the pure imaginary eigenvalues for $A+B$. We fix $\beta_0$ and use $\omega_0$ and $m_0$ as an initial guess to find $\omega$ and $m$ for a Hopf bifurcation point of system (2.1) with a fixed $\tau > 0$. It may be assumed that all other solutions of $det(A-\lambda I +Be^{-\lambda \tau}) = 0$ at the parameter point $(\beta_0, m)$ have negative real parts. Otherwise, another bifurcation curve other than $Hopf_0$, $SN$, $Tr_1$, and $Tr_2$ would appear in Figure 6(a). It is worth noting that $\omega$ and $m$ appear in several entries of the matrix $A-i\omega I +Be^{-i\omega\tau}$. The Newton method is not practical since it involves taking derivatives with respect to $\omega$ and $m$. In this paper, the Nelder-Mead simplex method ^[61], which is a direct method, is employed for finding $\omega$ and $m$. The detailed algorithm can be found in ^[61], or the built-in function "fminsearch" in Matlab can be directly applied for this purpose.

Mahasa et al. ^[16] studied a mathematical model of OVT with a time delay in the stimulation of virus-specific immune cells. The time delay $\tau$ in their work was fixed at 7 hours, and their work showed that the time delay $\tau = 7$ hours does not affect the stability of the virus-free equilibrium. In this work, we use $\tau = 1$ day and $\tau = 7$ days, which are at both ends of the parameter range shown in Table 2. When $(\beta_0, m)$ is not close enough to $(\beta_0, m_0)$, the continuation technique with some smaller increments in $\tau$ helps with the success of identifying $m$ for $\tau = 1$. With each point on $Hopf_0$ as an initial guess for the corresponding point on $Hopf_1$, parallel computation helps accelerate the process. Hopf bifurcation curves $Hopf_1$ with $\tau = 1$, $Hopf_2$ with $\tau = 1$, and $Hopf_3$ with $\tau = 7$ are computed and shown in Figure 6(b). The equilibrium $E_{41}$ changes its stability on the curves $Hopf_1$ and $Hopf_3$ when $\tau = 1$ and $\tau = 7$, respectively. The curve $Hopf_3$ is above $Hopf_1$ which is above $Hopf_0$, indicating that the presence of time delay destabilizes the system. For Eq (2.1) with $\tau = 1$, the curve $Hopf_2$ shows Hopf bifurcations where $E_2$ changes its stability. The equilibrium $E_2$ is stable in the region bounded by the curves $Tr_1$ and $Hopf_2$. This shows again time delay exhibits a destabilization effect in Eq (2.1). Theorem 4.2 gives a condition for the stability of $E_2$ when $\tau > 0$. The stability condition is given in Eq (4.12). Equation (4.12) is independent of $\beta$. Therefore, we examine the condition with the parameter values used in Figure 6. According to Theorem 4.2, the equilibrium $E_2$ is stable for $\tau > 0$ when $m > 1.208 \times 10^7$. We then use the numerical method proposed in this section to compute the Hopf bifurcation for $E_2$ for $\tau < 1000$ and find that $E_2$ is stable when $m > 1.205 \times 10^7$ for $\tau < 1000$. Our numerical computation agrees with Eq (4.12). Since such large $\tau$ values are not within the plausible range of $\tau$, the figure is not shown in this paper.

5.4.   Virus transmission rate and the steepness of immune cell proliferation

In this subsection, we study the effect of virus transmission rate $\beta$ and the steepness of immune cell proliferation $m$ on the outcome of treatment. Let $\beta = 0.1$, which corresponds to a relatively large virus transmission rate for the range of $\beta$ in Table 2. Let $m = 2\times 10^6$. Note that this set of parameter values lies in the region where both $E_2$ and $E_{41}$ are stable, as shown in Figure 6(a). The tumor persists with a population size of $T_s+T_r+T_i = \gamma m/(\alpha-\gamma)$. Furthermore, the tumor-free equilibrium $E_0$ is unstable, meaning that either the treatment may fail or tumor relapse may occur. This simulation can help identify factors that may cause unsuccessful treatment or tumor recurrence and further help design strategies to improve OVT treatment.

Consider the initial condition $(7\times 10^6,100, 0, 5\times 10^6, 2\times 10^3)$ for treating a tumor of $7\times 10^6$ cells with a therapeutic dose of $5\times 10^6$ pfu. Figure 7(a) shows that the population dynamics approach the stable equilibrium $E_{41}$. The stability of $E_{41}$ implies that the virus is persistent. Virus persistence has been observed in experimental studies and proven as a means to slow down tumor growth ^[62]. The sensitive tumor cells (red curve) are infected quickly after the treatment, and thus the infected tumor cells (black curve) increase rapidly. The tumor cells respond to the treatment initially but develop resistance to the therapy causing an increase in resistant tumor cells (blue curve). Note that the elimination of resistant tumor cells relies on immune cells (green curve), and the immune system is not able to control the resistant tumor cells. Strategies to reduce resistant tumor cells or enhance the immune response to resistant tumor cells are required to overcome this challenge.

Let $\tau = 1$, and all parameter values remain the same as used for Figure 7(a). Figure 7(b) shows population oscillations. The equilibrium $E_{41}$ becomes unstable, demonstrating the destabilization effect caused by time delay. A slightly larger $m$ value is used with $\tau = 1$ and $\tau = 7$. Let $m = 3\times 10^6$. The equilibrium $E_{41}$ is stable when $\tau = 1$ and unstable when $\tau = 7$, where the population dynamics are similar to the population dynamics in Figure 7(a), (b), respectively. Equation (2.1) with $\tau > 0$ exhibits a destabilization effect. Treatment is not successful in either case, $\tau = 0$ or $\tau > 0$.

Figure 7(a), (b) shows that the OVT is effective initially when the transmission rate $\beta = 0.1$, where the susceptible tumor cells are infected and killed. The immune cells go through the contraction phase after eliminating the infected cells. However, part of the susceptible tumor cells develop resistance to treatment, as shown in Figure 7(c), which presents a closer look at Figure 7(a). The resistant tumor population grows progressively, and the treatment ultimately fails. Consider $m = 10^4$, which corresponds to a quicker immune cell proliferation in response to tumor growth. Figure 7(d) shows that a quicker immune response eliminates infected cells and viruses but not the other two populations of tumor cells. The anti-viral response, where viruses are cleared by the immune system, is a challenge in OVT ^[45,46]. Furthermore, a low half-saturation constant tends to cause an oscillation of populations ^[57,58]. Such persistent large scale oscillations are unlikely to be seen in tumor immune interactions ^[63]. Nevertheless, oscillations in tumor and immune cell populations may account for the transitional state between elimination and equilibrium phases in immunoediting ^[63].

Let $\beta = 10^{-6}$, which corresponds to a small virus transmission rate compared with the previous simulation. Let $m = 2\times 10^6$. Figure 8 shows that the OVT is ineffective when the virus transmission rate $\beta$ is small. The infected tumor cell population is small and eliminated soon, but the immune system is not able to eliminate the other two populations of tumor cells. This result agrees with the result shown in Figure 5(c), (f). Again, a quicker immune response by using a smaller half saturation constant $m$ or a time delay $\tau > 0$ does not improve the outcome of OVT but causes an oscillation of populations. The figures are not shown here as the population dynamics are similar to those in Figure 7(d).

It is known that the interactions between the immune system and oncolytic viruses, resulting in both anti-viral and anti-tumoral responses, can have both beneficial and detrimental effects on the outcome of OVT ^[3,45,64]. Several studies have suggested the use of combination therapy with other immunotherapeutic approaches to improve the efficacy of treatment ^[64,65,66]. Guedan and Alemany ^[67] have published a review article addressing the potential role of combination therapy using OVT and CAR T-cell therapy in fighting cancer. It is expected that combination therapy using OVT and CAR T-cell therapy can improve the outcome of the OVT studied in this subsection.

5.5.   Combination therapy using OVT and CAR T-cell therapy

Chimeric antigen receptor T-cell (CAR T-cell) immunotherapy is a novel revolutionary cancer treatment. The process of CAR T-cell therapy involves deriving T cells from the patient and engineering these T cells in vitro to promote the recognition of cancer cells and improve T-cell function. CAR T cells are then expanded and infused back into the patient to enhance the elimination of tumor cells ^[68,69,70]. Preclinical studies in combination therapies of OVT and CAR T-cell therapy have been conducted to investigate the synergistic effects ^[71], and the results have shown sustained anti-tumor immune response yielding promising treatment outcomes ^[72].

According to Figure 5 (c), (d), the killing rate must be at least $5\times 10^{-4}$ to effectively fight tumor cells. From Section 5.4, the development of resistance to OV is a challenge to the success of OVT. Assume that the engineered T cells have higher killing rates of susceptible and resistant cells than used in Section 5.4. Let $k_s = 5\times 10^{-4}$, $k_r = 5 \times 10^{-3}$. Consider $\beta = 0.1$, and $m = 2\times 10^6$. Gruber et al. studied the relationship between tumor cells and T lymphocytes in breast cancer patients, and the patients had CTL counts of 365 $\pm$ 194 cells/$\mu$L (Mean $\pm$ SD) ^[73]. Within two standard deviations of the mean, it is safe to say that CTL counts for most patients range from several tens to 750 cells/$\mu L$. The dose levels for CAR T-cell therapy are reported to range from 60 million to 600 million cells ^[74,75]. For a person with body weight of 70 kg, about 7.2 percent of body weight is blood ^[76]. That is about 5 liters of blood. The concentration of CAR T cells in blood is about 12–120 cells/$\mu L$. The goal of CAR T cell therapy is to generate a sufficient number of CAR T cells to go into the tumor site to effectively target and attack the cancer cells ^[77]. From the simulation in the previous section, the numbers of effector cells in the stable equilibria $E_2$ (Figure 8) and $E_{41}$ (Figure 7(a)) are 43,000 and 2000, respectively. In this example, we simulate the concurrent administration of both agents. Using the initial condition $(7\times 10^6,100, 0, 5\times 10^6, 8000)$, where the number of the initial immune cells is 8000, Figure 9(a) shows that combination therapy is able to eradicate tumor cells. After the elimination of tumor cells, the immune cells go through the contraction phase. Further simulation shows that monotherapy using either CAR T-cell immunotherapy ((Figure 9(b)) or OVT (Figure 7(a)) is not able to eliminate the tumor. It has been reported that CAR T cell therapy is ineffective for solid tumors ^[77]. Combination therapy produces a synergistic effect ^[72].

Consider a larger tumor with $10^9$ cells, which corresponds to a tumor of approximately 1 cm in diameter ^[34]. It is known that virus clearance mediated by innate and adaptive immune cells is a significant challenge in OVT ^[78]. The innate immune response acts rapidly ^[79,80]. We use time delay to simulate the strategy to circumvent this virus clearance in the initial days. Fu et al. ^[62] conducted experiments to show that genetically coating oncolytic virus with CD47 allows OVs to evade the immune response and promotes virus persistence at the tumor site. We use a lower $k_v$ value to simulate evasion of the immune response and to enhance virus persistence. Let $k_v = 0.03$ in this simulation. Figure 10(a) shows that the combination treatment is successful with a low dose of 5000 CAR T cells, compared with the simulation in Figure 9(a). The results of further simulation of other cases are summarized in Table 3. Combination therapy allows the use of lower doses of each agent and achieves better efficacy. According to Eldar-Boock et al. ^[81], combination therapy can reduce the toxicity and side effects due to lower doses.

Recall that $E_0$ is unstable. This implies a small perturbation may cause population dynamics to move away from $E_0$. In Figure 10(b) we simulate that one susceptible cell is formed at $t = 1000$. Note that 5–10% of T cells become memory T cells during the contraction phase^[82,83]. These memory T cells can quickly respond to previously encountered pathogens ^[82]. The number of immune cells increases in response to the presence of tumor cells. However, the immune cells cannot control the growth of tumor cells leading to tumor recurrence. It has been reported that tumor recurrence is a formidable challenge for cancer treatment ^[84].

6.   Conclusions

Drug resistance is a common phenomenon encountered in cancer therapy, posing a critical barrier that requires urgent overcoming. The issue of viral resistance is also encountered in oncolytic viral therapy. While various mathematical and computational models that have appeared in the literature focus on investigating resistance to chemotherapy, only one article ^[15] employs differential equations to address resistance in OVT.

In this study, we extended the work of ^[15] by incorporating immune cells into the tumor-virus interaction. These immune cells have the capability to eliminate both tumor cells and viruses. We established the existence and stability of boundary equilibria and presented global stability results. We proved that the delay in immune cell recruitment cannot destabilize boundary equilibria in which immune cells are not present. For those boundary equilibria with the presence of immune cells, we provided sufficient conditions based on model parameters for which the delay has no destabilization effects. Additionally, we derived a critical delay value, under which a boundary equilibrium with the presence of immune cells is destabilized when the delay passes through this critical value, leading to tumor oscillations. The introduction of immune cells results in an additional virus-free equilibrium when the recruitment rate of immune cells is sufficiently high. At this virus-free equilibrium, the total tumor burden is smaller compared to the virus-free equilibrium without immune cells. Consequently, we conclude that immune cells can effectively reduce the tumor load under conditions of sufficient immune strength. We also demonstrated numerically, using reasonable parameter values in Section 5.1, that the inclusion of immune cells can lower the tumor burden compared to the model where immune cells are not included.

GSA showed that the virus transmission rate $\beta$ and the parameters associated with immune response, such as immune cell proliferation rate $\alpha$ and tumor killing rate $k_s$, were important factors that affected the final tumor population (Figure 4). Strengthening the immune system is crucial for treatment success and preventing tumor recurrence ^[51]. The numerical simulation demonstrated that the virus transmission rate must exceed a certain threshold to induce an anti-tumoral effect (Figure 5(c), (f)). A numerical method of Hopf bifurcations was developed for system (2.1) with a time delay. Bifurcation analysis (Figure 6), using the virus transmission rate and the half-saturation constant of immune cell proliferation rate as bifurcation parameters, showed that the system (2.1) exhibited rich dynamics. The results of bifurcation analysis agreed with the mathematical analysis in Sections 3 and 4. The time delay $\tau$ in the activation of immune cells had a destabilization effect (Figures 6 and 7(a)–(d)).

Numerical simulation showed that the OVT was effective if the virus transmission rate was sufficiently high; however, the treatment may ultimately fail due to the development of resistance (Figure 7(c)). Resistance to oncolytic viruses has been reported to jeopardize the success of treatment ^[4]. Another simulation for a lower virus transmission rate showed that the viruses were cleared out before susceptible cells were infected (Figure 8), and thus the OVT was ineffective. Virus clearance is known to be a challenge in OVT ^[45,46]. Literature has suggested that OVT as monotherapy frequently fails to control tumor progression and sustain a reduction in tumor masses due to anti-viral effect or development of resistance ^[3,45,64,85]; combination therapy is a preferable solution to overcome such challenges ^{[65,67,86,87]}.

To enhance the success of treatment, combination therapy using OVT and CAR T-cell therapy was considered in this paper. Numerical simulation showed that the combination therapy produced synergistic effects and enhanced the success of treatment (Figure 9(a)). Monotherapy using either OVT or CAR T-cell therapy is ineffective (Table 3). Further simulation studied treatment strategies to improve the outcome of treatment (Table 3). Circumventing immune clearance of OVs in the initial days allowed the use of a lower dose of CAR T cells and concurrently promoted the efficacy of the combination therapy. This agrees with the viewpoint of experimental and clinical reports. Virus infection and replication during the initial days of OVT significantly impact its efficacy. ^[47,48]; one goal of combination therapy is to improve efficacy while reducing toxicity ^[81]. Additional simulation showed that Prolonged virus persistence enhanced the outcome of treatment, which agreed with the experimental report that prolonging virus persistence improved the efficacy of OVT ^[62].

Finally, our simulation showed that cancer recurrence occurred after a tumor cell was formed (Figures 10(b)). Evidence has suggested that cancer recurrence has been a major challenge after the initial treatment response ^[88]. Our simulation results agreed with clinical and experimental research studies. It is expected that these results can be applied to clinical practice. Further studies including clinical tests are needed to further verify these results.

The immune cell assumptions made in this work were simplified. In future research, we intend to introduce more complexity into the interaction. We propose segregating them into two groups: tumor-specific immune cells and virus-specific immune cells. Tumor-specific immune cells will be activated by both resistant and sensitive tumor cells, capable of targeting and eliminating both types of cancer cells. On the other hand, virus-specific immune cells will be recruited by infected tumor cells and viruses, focusing solely on eliminating infected tumor cells and viruses. This revised model will encompass six nonlinear ordinary differential equations. Furthermore, we may introduce time delays, either in the viral cycle or in immune cell recruitment, to further enhance the model's accuracy.

Use of AI tools declaration

The authors declare they have not used AI tools in the creation of this article.

Acknowledgments

We appreciate all four reviewers for their valuable suggestions, which have significantly improved the paper. The research of Hsiu-Chuan Wei was partially supported by the National Science and Technology Council of Taiwan under the grant NSTC 112-2115-M-035-006. The work of S. Jang was partially supported through Travel Support for Mathematicians from the Simons Foundation under MPS-TSM-00002532.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

A.   Special cases of ODE models

The ODE models of the special cases $d = 0$ and no resistant cancer cells are briefly discussed in this appendix.

A.1.   The case of $d=0$

Since $d = 0$, it follows that $(\hat T_s, \hat T_r, 0, 0, 0) = (0, K, 0, 0, 0)$. If the maximal immune cell proliferation rate $\cfrac{\alpha K}{m+K}$ is greater than its death rate $\gamma$, then the interaction is able to support another equilibrium, $E_2^0 = (0, \xi, 0, 0, \tilde Z^0)$, where

It can be verified that $E_2^0$ is asymptotically stable if

and unstable if the above inequality is reversed. Moreover, if

then the model (2.1) with $\tau = 0$ and $d = 0$ has an equilibrium of the form $E_3^0 = (\bar T_s^0, \bar T_r^0, 0, 0, \bar Z^0)$, where

Observe that the last inequality in (A.3) implies

Therefore, the existence of $E_3^0$ implies that $E_2^0$ is unstable. The Jacobian matrix at $E_3^0$ has the same entries as $J(E_2)$ defined in (3.14). We can conclude that $E_3^0$ is asymptotically stable if (3.11) and (3.16) are satisfied with $d$ being replaced by $0$. The existence and stability of boundary equilibria of (2.1) with $\tau = 0$ and $d = 0$ are summarized in Table A1.

Comparing Table A1 with Table 1, it is observed that $E_1^0 = E_1$ when $d = 0$, as in this instance $\hat T_s = 0$. Therefore, the condition $c > \beta q \hat T_s$ given in Table 1 is trivially true and the stability of $E_1^0$ in Table A1 reduces to $\cfrac{\alpha K}{m+K} < \gamma$ only. Further, the second inequality in (A.3) implies $\xi < K$, i.e., $\cfrac{\alpha K}{m+K} > \gamma$. Thus, the equilibrium $E_3^0$ given in Table A1 is the equilibrium $E_2$ given in Table 1, and therefore the stability of $E_3^0$ is based on $E_2$. The existence of $E_3$ in Table 1 requires that $c < \beta q \hat T_s$, which is impossible in the case $d = 0$, and there is no such corresponding equilibrium in Table A1.

Parallel to Theorem 3.1, we can also derive sufficient conditions for the global stability of $E_1^0$ as follows, without providing a proof.

Theorem A.1. Let $\tau = 0$ and $d = 0$. The equilibrium $E_1^0 = (0, K, 0, 0, 0)$ is globally asymptotically stable in $D$ for (2.1) if $\frac{\alpha K}{m+K} < \gamma$ and $c > \beta q K$.

A.2.   The scenario of no resistant tumor cells

Assume that there are no resistant tumor cells initially, $T_r(0) = 0$, and no mutation, $a = 0$. Then, $T_r(t) = 0$ for all $t > 0$, and model (2.1) with $\tau = 0$ becomes a four-dimensional system involving only $T_s, T_i, V$, and $Z$. If, in addition, there are no infected tumor cells and viruses initially, $T_i(0) = 0 = V(0)$, then $T_i(t) = 0 = V(t)$ for all $t > 0$, and we obtain the following two-dimensional system:

The system (A.5) has three equilibria, where $(0, 0)$ and $(K, 0)$ always exist, and $(\xi, \tilde Z)$ is biologically feasible only if $\cfrac{\alpha K}{m+K} > \gamma$, with $\xi = \dfrac{\gamma m }{\alpha-\gamma}$ defined in (3.5) and

It is straightforward to see that $(0, 0)$ is always unstable, while $(K, 0)$ is asymptotically stable if $\dfrac{\alpha K}{m+K} < \gamma$ and unstable if $\dfrac{\alpha K}{m+K} > \gamma$. We assume $\dfrac{\alpha K}{m+K} > \gamma$ so that $(\xi, \tilde Z)$ exists. It follows from the Dulac criterion ^[24] that (A.5) has no positive periodic solutions, which implies that $(\xi, \tilde Z)$ is globally asymptotically stable in $\mathring {\mathbb R}^2_+$ for the two-dimensional model (A.5) by the Poincar$\acute{\rm e}$-Bendixson Theorem ^[24]. As a result, we see that $\tilde E_2 = (\xi, 0, 0, \tilde Z)$ is asymptotically stable if

Let $c < \beta qK$ so that $\tilde E_3 = (\tilde T_s, \tilde T_i, \tilde V, 0)$ exists, where

See Section 3.3 of ^[15]. Applying the Routh-Hurwitz criterion ^[24], without delving into all the details, $\tilde E_3$ is locally asymptotically stable if it satisfies the following conditions

The above results are summarized in the following table.

Parallel to Theorem 3.1, we can also derive sufficient conditions for the global stability of $\tilde E_1$ as follows, without providing a proof.

Theorem A.2. Let $\tau = 0$, $a = 0$, and $T_r(0) = 0$. The equilibrium $\tilde E_1 = (K, 0, 0, 0)$ is globally asymptotically stable in $\{(T_s, T_i, V, Z)\in\mathbb{R}_+^4: T_s > 0\}$ if $\frac{\alpha K}{m+K} < \gamma$ and $c > \beta q K$.

Notice the stability of $\tilde E_2$ given in Table A2, (A.7), is the stability condition (3.16) for $E_2$ in Table 1, while the other condition (3.11) is trivially true for the case of no resistant tumor cells $a = 0$ and $T_r\equiv 0$. In addition, the equilibrium $\tilde E_1$ in Table A2 corresponds to the equilibrium point $E_1 = (\hat T_s, \hat T_r, 0, 0, 0)$ in Table 1 with $T_r\equiv 0$, and they have the same existence condition. In addition, the inequalities (3.20) are equivalent to those given in (A.9), and therefore both equilibria have the same stability conditions.

B.   Special cases of DDE models

B.1.   The case of $d=0$

The resulting model always has equilibria $E_0 = (0, 0, 0, 0, 0)$, and $E_1^0 = (0, K, 0, 0, 0)$. Similar to Section 4.2, it can be seen that $E_0$ is always unstable for $\tau\geq 0$ and $E_1$ is asymptotically stable for $\tau\geq 0$ if $\cfrac{\alpha K}{m+K} < \gamma$, and unstable for $\tau\geq 0$ if $\cfrac{\alpha K}{m+K} > \gamma$.

Let $\cfrac{\alpha K}{m+K} > \gamma$. Then, $E^0_2 = (0, \xi, 0, 0, \tilde{Z^0})$ exits, where $\xi = \cfrac{m\gamma}{\alpha-\gamma}$, and $\tilde{Z^0} = \dfrac{r_r}{k_r}(1-\dfrac{\xi}{K})$. It can be argued that $E_2^0$ is asymptotically stable for $\tau\geq 0$ if certain conditions on the parameters are met. On the other hand, under certain constraints on the parameters, there exists a critical delay $\tau_c$ beyond which $E_2^0$ becomes unstable as $\tau$ increases. Indeed, let

We summarize the following without proof as the proof is similar to that given previously.

Theorem B.1. Let $d = 0$, $\cfrac{\alpha K}{m+K} > \gamma$, $a_{11} < 0$, and $p_2-q_2 < 0$. Then, there exists $\tau_0 > 0$ such that $E^0_2$ is asymptotically stable for $\tau \in [0, \tau_0)$ and $\dfrac{d(Re\lambda)}{d\tau}\Big|_{\tau = \tau_0} > 0$.

Assume (A.3) so that $E_3^0 = (\bar T_s^0, \bar T_r^0, 0, 0, \bar Z^0)$ exists. Suppose that (3.11) and (3.16) are satisfied with $d = 0$, i.e., $E_3^0$ is asymptotically stable for the ODE model with $d = 0$. Then, one can also find a critical delay $\tau_c$ such that $E_3^0$ is asymptotically stable for $\tau\in[0, \tau_c)$, and the transversality condition for $\lambda$ can also be verified. We do not present the results here since the proof and conclusion are similar.

B.2.   The scenario of no resistant tumor cells

The proof of the following proposition is routine and is omitted.

Proposition 8.1. The following is true for system (2.1) with $a = 0$ and $T_r(0) = 0$.

(a) $E_0 = (0, 0, 0, 0)$ is unstable for $\tau\geq 0$.

(b) $\tilde{E_1} = (K, 0, 0, 0)$ is asymptotically stable for $\tau\geq0$ if $c > \beta qK$ and $\dfrac{\alpha K}{m+K} < \gamma$, and unstable if either $c < \beta q K$ or $\dfrac{\alpha K}{m+K} > \gamma$.

Let

The delay $\tau$ can affect the stability of $\tilde E_2$ as illustrated below.

Theorem B.2. Let $\cfrac{\alpha K}{m+K} > \gamma$. Then, $\tilde E_2 = (\zeta, 0, 0, \tilde Z)$ exists, and

(a) $\tilde{E_2}$ is asymptotically stable for $\tau\geq 0$ if $(k_i\tilde Z+b)(k_v\tilde Z+c) > qb \beta\zeta$ and $\tilde p_2-\tilde q_2 > 0$.

(b) If $(k_i\tilde Z+b)(k_v\tilde Z+c) < qb \beta\zeta$ or $\tilde p_2-\tilde q_2 < 0$, then there exists $\tau_0 > 0$ such that $\tilde E_2$ is asymptotically stable for $\tau\in[0, \tau_0)$ and $\dfrac{d(Re\lambda)}{d\tau}\Big|_{\tau = \tau_0}\neq 0$.

C.   Proofs of mathematical results

Proof of Theorem 2.1. Since $f$ and $\cfrac{\partial f}{\partial X}$ exist and are continuous on $\mathbb{R}_+^5\times\mathbb{R}_+^5$, (2.1) has a unique solution on $[-\tau, t_0)$ for some $t_0 > 0$. Moreover, if $X, Y \geq 0$ with $x_j = 0$, then $f_j(X, Y)\geq 0$, and thus solutions remain nonnegative on $[-\tau, t_0)$ by ^[27].

As

and

we have $T_s(t)+T_r(t)\leq K$ on $[0, t_0)$. Next,

on $[0, t_0)$ for some $C_1 > 0$ and $b_1 = \min\{r_s, r_r, b\} > 0$. It follows that $T_i(t)\leq \hat{C_1}$ for $0\leq t < t_0$ for some $\hat{C_1} > 0$. Then, $V^{\prime}(t)\leq qb\hat{C_1}-cV$, and we have $V(t)\leq C_2$ on $[0, t_0)$ for some $C_2 > 0$.

Notice

Let

Then, by the given assumption,

on $[0, t_0)$ for some $\hat{M} > 0$, and

Hence, $M(t)$ is bounded on $[0, t_0)$, and solutions cannot blow up as $t \uparrow t_0$. Therefore, solutions exist on $[0, \infty)$ and remain nonnegative. Similar arguments can be applied to show that solutions are bounded on $[0, \infty)$. ■

Proof of Proposition 3.1. Clearly, $(\hat T_s, \hat T_r, 0)$ is asymptotically stable, and

for all $t \geq 0$ implies $\lim_{t\to\infty} Z(t) = 0$ by the assumption. Thus, (3.1) is asymptotically autonomous ^[30] to the two-dimensional $T_sT_r$ subsystem, under which $(\hat T_s, \hat T_r)$ is globally asymptotically stable in the region for which $T_s(0)+T_r(0) > 0$ by ^[15]. Therefore, $(\hat T_s, \hat T_r, 0)$ is globally asymptotically stable in $\Gamma$. ■

Proof of Theorem 3.1. By Proposition 3.4, $\lim\limits_{t \to \infty}{T_i(t)} = 0 = \lim\limits_{t \to \infty}{V(t)}$, and thus system (2.1) with $\tau = 0$ is asymptotically autonomous ^[30] to the three-dimensional model (3.1), where $(\hat T_s, \hat T_r, 0)$ is globally asymptotically stable in $\Gamma$ by Proposition 3.1. Since $c > \beta qK > \beta q \hat T_s$ and $\frac{\alpha K}{m+K} < \gamma$, $E_1$ is asymptotically stable. It follows that $E_1$ is globally asymptotically stable in $D$. ■

Proof of Proposition 4.1. We only need to prove the second part of (b) for $\tau > 0$. Let

with $\tau > 0$. Then,

for $\lambda\geq 0$, and

Therefore, $h(\lambda) = 0$ has one positive root, and hence $(\hat T_s, \hat T_r, 0)$ is unstable for all $\tau > 0$. ■

Proof of Theorem 4.1 Suppose (4.9) has $j$ positive roots $x_n, \ 1\leq j\leq 3$, each of which is simple. Then, (4.4) has $j$ simple pure imaginary roots

Clearly, for each $n$, $1\leq n \leq j$, (4.10) and (4.11) have a unique solution $\omega_n \tau$ in $(0, 2\pi]$, $\omega_n = \sqrt{x_n}$. Let

if $\rho_s > 0$, and

if $\rho_s\leq 0$. Define

and set $\omega_0 = \omega_{n_0}$. We verify the transversality condition at $\tau = \tau_0$, i.e., $\frac{d Re(\lambda)}{ d\tau}|_{\tau = \tau_0}\neq 0$, by using

The transversality condition is a necessary condition for a Hopf bifurcation to occur ^[27].

Implicitly differentiate (4.4) with respect to $\tau$, and we obtain

It follows that

and by (4.8),

As a result,

The denominator of the above expression is clearly positive, and the numerator is $F^{\prime}(x_0)$, where $x_0 = \omega_0^2$ and $F$ is defined in (4.9). Since $x_0$ is a simple root, $F^{\prime}(x_0)\neq0$, and the transversality condition is satisfied. ■