Natural Killer Cells Recruitment in Oncolytic Virotherapy: A Mathematical Model

Abstract

In this paper, we investigate how natural killer (NK) cell recruitment to the tumor microenvironment (TME) affects oncolytic virotherapy. NK cells play a major role against viral infections. They are, however, known to induce early viral clearance of oncolytic viruses, which hinders the overall efficacy of oncolytic virotherapy. Here, we formulate and analyze a simple mathematical model of the dynamics of the tumor, OV and NK cells using currently available preclinical information. The aim of this study is to characterize conditions under which the synergistic balance between OV-induced NK responses and required viral cytopathicity may or may not result in a successful treatment. In this study, we found that NK cell recruitment to the TME must take place neither too early nor too late in the course of OV infection so that treatment will be successful. NK cell responses are most influential at either early (partly because of rapid response of NK cells to viral infections or antigens) or later (partly because of antitumoral ability of NK cells) stages of oncolytic virotherapy. The model also predicts that: (a) an NK cell response augments oncolytic virotherapy only if viral cytopathicity is weak; (b) the recruitment of NK cells modulates tumor growth; and (c) the depletion of activated NK cells within the TME enhances the probability of tumor escape in oncolytic virotherapy. Taken together, our model results demonstrate that OV infection is crucial, not just to cytoreduce tumor burden, but also to induce the stronger NK cell response necessary to achieve complete or at least partial tumor remission. Furthermore, our modeling framework supports combination therapies involving NK cells and OV which are currently used in oncolytic immunovirotherapy to treat several cancer types.

Appendices

Well-Posedness

The proposed model in this study describes the temporal evolution of cells and virus populations, and therefore, the cell concentrations should remain nonnegative and bounded. Here, we establish the well-posedness of the immune-free submodel (Eqs. (3.1.1)–(3.1.3)) and its proofs. We should emphasize that the nonnegativity and boundedness of the immunocompetent model directly follow from the same type arguments presented herein, and hence, we shall omit their proofs. We also provide a brief outline for extending the solutions to our immunocompetent system (see Eqs. (2.1.1)–(2.1.4)). Note that when solving for one state variable, for simplicity, other state variables appearing in the equation are treated as constants.

Theorem A.1

Well-Posedness

• (i) (nonnegativity of solutions) Given that the nonnegative initial conditions $(T_{u0}> 0, T_{i0}> 0, V_0 > 0)$, the corresponding solutions $(T_{u}(t), T_{i}(t), V(t))$ will remain nonnegative for all $t \in [0, \infty )$.

• (ii) (boundedness of solutions and invariant region) The model system is bounded and the invariant region is given by ${\varvec{\Omega }}_{\mathbf{V}} = \{(T_{u}, T_{i}, V) \in \mathbb {R}_{+}^3 \mid \ 0 \hbox {\,\,\char 054\,\,}T_{u} \hbox {\,\,\char 054\,\,}\theta _{T},\ 0 \hbox {\,\,\char 054\,\,}T_{i} \hbox {\,\,\char 054\,\,}\frac{1}{\delta } \beta T_{u} V,\ 0 \hbox {\,\,\char 054\,\,}V \hbox {\,\,\char 054\,\,}\frac{1}{\gamma }(\mu _{v}(t) + b T_{i})\}$. Moreover, the domain ${\varvec{\Omega }}_{\mathbf{V}}$ is positively invariant for the model and therefore biologically meaningful for the cell concentrations and regarded as a “global” domain. The corresponding dimensionless system, Eqs. (3.2.6)–(3.2.8), is valid under the following positively invariant domain: $\varvec{\Omega }_{\mathbf{V1}} = \{(T_{u}, T_{i}, V) \in \mathbb {R}_{+}^3 \mid \ T_{u} \ge 0,\ T_{i} \ge 0,\ 0 \hbox {\,\,\char 054\,\,}T_{u} + T_{i} \hbox {\,\,\char 054\,\,}1,\ V \ge 0\}.$

• (iii) (existence and uniqueness) For any nonnegative initial values of the model state variables, a solution to the model exists and is unique in the positively invariant domain $\varvec{\Omega }_{\mathbf{V}}$ for all time $t > 0$.

1.1 (i) Nonnegativity of Solutions

Proof

Considering Eq. (3.1.1), we have a Bernoulli differential equation. We now rewrite Eq. (3.1.1) as $\frac{\mathrm{d}T_{u}}{\mathrm{d}t} = \alpha T_{u} - \frac{\alpha T_{u}^2}{\theta _{T}} - \frac{\alpha T_{u} T_{i}}{\theta _{T}} - \beta T_{u} V = \left( \alpha - \frac{\alpha T_{i}}{\theta _{T}} - \beta V\right) T_{u} - \frac{\alpha T_{u}^2}{\theta _{T}}.$ Thus, the Bernoulli standard form of Eq. (3.1.1) is $\frac{\mathrm{d}T_{u}}{\mathrm{d}t} + \left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) T_{u} = -\frac{\alpha T_{u}^2}{\theta _{T}}$.

Here, $n = 2$, $P(t) = \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha$, $Q(t) = -\frac{\alpha }{\theta _{T}}$. The integrating factor is $I(t) = e^{\int {(1 - 2)\left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) }dt} = e^{ - \left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) t}.$ The solution is therefore given by: $T_{u}^{-1} = e^{\left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) t}\left[ \int {(-1)\left( -\frac{\alpha }{\theta _{T}}\right) e^{ - \left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) t}}dt + c\right] = e^{\left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) t}\left[ \frac{\alpha }{\theta _{T}}\int {e^{ - \left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) t}}dt + c\right] = e^{nt}\left[ -\frac{\alpha }{\theta _{T} n} e^{-nt} + c\right] ,\quad \text {where} \quad n = \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha = -\frac{\alpha }{\theta _{T} n} + ce^{nt}.$

At $t = 0$, $T_{u} = T_{u0}$, solving for c we obtain $c = T_{u0} + \frac{\alpha }{\theta _{T} n}$. Now substituting c in the equation for $T_{u}^{-1}$ above, we get $T_{u}^{-1} = \frac{(\theta _{T} n T_{u0} + \alpha )e^{nt} - \alpha }{\theta _{T} n}$, and thus,

$$\begin{aligned} T_{u}&= \frac{\theta _{T} n}{\left( \theta _{T} n T_{u0} + \alpha \right) e^{nt} - \alpha }\\&= \frac{\theta _{T} \left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) }{\left( \theta _{T} T_{u0} \left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) + \alpha \right) e^{\left( \beta V + \frac{\alpha T_{i}}{\theta _{T}} - \alpha \right) t} - \alpha } >0. \end{aligned}$$

Considering Eq. (3.1.2), we have a first-order ODE which can easily be written as $\frac{\mathrm{d}T_{i}}{\mathrm{d}t} + \delta T_{i} = \beta T_{u} V$. Using the following integrating factor $I(t) = e^{\int {\delta }dt} = e^{\delta t},$ and solving we get, $T_{i} = \frac{1}{\delta } \beta T_{u} V + ce^{-\delta t}.$ At $t = 0$, $T_{i}(t) = 0$, solving for c we obtain $c = -\frac{1}{\delta } \beta T_{u} V.$ Substituting c in the equation of $T_{i}$ above, we get $T_{i} = \frac{1}{\delta } \beta T_{u} V - \frac{1}{\delta } \beta T_{u} V e^{-\delta t} = \frac{1}{\delta } \beta T_{u} V \left( 1 - e^{-\delta t}\right) > 0.$

Similarly, considering Eq. (3.1.3), we have a first-order ODE which can be rewritten as $\frac{\mathrm{d}V}{\mathrm{d}t} + \gamma V = \mu _{v}(t) + b T_{i}$. The integrating factor is given by $I(t) = e^{\int {\gamma }dt} = e^{\gamma t},$ which is used to solve for $V = \frac{1}{\gamma }\left( \mu _{v}(t) + b T_{i}\right) + ce^{-\gamma t}.$ At $t = 0$, $V(t) = 0$, solving for c we obtain $c = -\frac{1}{\gamma }(\mu _{v}(t) + b T_{i}).$ Substituting c in the equation of V above, we get $V = \frac{1}{\gamma }(\mu _{v}(t) + b T_{i})(1 - e^{-\gamma t}) > 0.$

1.2 (ii) Boundedness of Solutions and Invariant Region

In this section, we discuss the boundedness of the solutions of our model.

For the uninfected population, Eq. (3.1.1) can be rewritten as

$$\begin{aligned} \frac{\mathrm{d}T_{u}}{\mathrm{d}t} = \alpha T_{u} - \frac{\alpha T_{u}^2}{\theta _{T}} - \frac{\alpha T_{u} T_{i}}{\theta _{T}} - \beta T_{u} V. \end{aligned}$$

From which we have

$$\begin{aligned} \frac{\mathrm{d}T_{u}}{\mathrm{d}t}&\hbox {\,\,\char 054\,\,}\alpha T_{u} - \frac{\alpha T_{u}^2}{\theta _{T}}\nonumber \\&= \frac{\theta _{T} \alpha T_{u} - \alpha T_{u}^2}{\theta _{T}}\nonumber \\&\quad \implies \frac{\theta _{T} dT_{u}}{\theta _{T} \alpha T_{u} - \alpha T_{u}^2} \hbox {\,\,\char 054\,\,}dt\nonumber \\&\quad \frac{\theta _{T} dT_{u}}{\alpha T_{u}(\theta _{T} - T_{u})} \hbox {\,\,\char 054\,\,}dt \end{aligned}$$

(A.2.1)

Integrating the left-hand side of Eq. (A.2.1) using partial fractions, we have

$$\begin{aligned} \frac{1}{T_{u}(\theta _{T} - T_{u})}&= \frac{A}{T_{u}} + \frac{B}{\theta _{T} - T_{u}} \end{aligned}$$

Comparing coefficients and solving for A and B, we get $A = \frac{1}{\theta _{T}}$ and $B = \frac{1}{\theta _{T}}$. So

$$\begin{aligned} \frac{1}{T_{u}(\theta _{T} - T_{u})} = \frac{1}{\theta _{T} T_{u}} + \frac{1}{\theta _{T} (\theta _{T} - T_{u})}. \end{aligned}$$

Integrating the above yields

$$\begin{aligned} \frac{1}{\theta _{T}} (ln(T_{u}) - ln(\theta _{T} - T_{u})). \end{aligned}$$

So, the integral for Eq. (A.2.1) is

$$\begin{aligned} \frac{\theta _{T}}{\alpha }\left( \frac{ln(T_{u})}{\theta _{T}} - \frac{ln(\theta _{T} - T_{u})}{\theta _{T}}\right) \hbox {\,\,\char 054\,\,}t + c. \end{aligned}$$

At $t = 0$, $T_{u} = T_{u0}$

$$\begin{aligned} c = \frac{\theta _{T}}{\alpha }\left( \frac{ln(T_{u0})}{\theta _{T}} - \frac{ln(\theta _{T} - T_{u0})}{\theta _{T}}\right) , \end{aligned}$$

so

$$\begin{aligned} \frac{\theta _{T}}{\alpha }\left( \frac{ln(T_{u})}{\theta _{T}} - \frac{ln(\theta _{T} - T_{u})}{\theta _{T}}\right)&\hbox {\,\,\char 054\,\,}t + \frac{\theta _{T}}{\alpha }\left( \frac{ln(T_{u0})}{\theta _{T}} - \frac{ln(\theta _{T} - T_{u0})}{\theta _{T}}\right) \\ \frac{1}{\alpha }(ln(T_{u}) - ln(\theta _{T} - T_{u}))&\hbox {\,\,\char 054\,\,}t + \frac{1}{\alpha }(ln(T_{u0}) - ln(\theta _{T} - T_{u0}))\\ ln(T_{u}) - ln(\theta _{T} - T_{u})&\hbox {\,\,\char 054\,\,}\alpha t + ln(T_{u0}) - ln(\theta _{T} - T_{u0})\\ ln\left( \frac{T_{u}}{\theta _{T} - T_{u}}\right)&\hbox {\,\,\char 054\,\,}\alpha t + ln\left( \frac{T_{u0}}{\theta _{T} - T_{u0}}\right) \\ ln\left( \frac{T_{u}}{\theta _{T} - T_{u}}\right) - ln\left( \frac{T_{u0}}{\theta _{T} - T_{u0}}\right)&\hbox {\,\,\char 054\,\,}\alpha t\\ ln\left( \frac{T_{u}(\theta _{T} - T_{u0})}{T_{u0}(\theta _{T} - T_{u})}\right)&\hbox {\,\,\char 054\,\,}\alpha t\\ \frac{T_{u}(\theta _{T} - T_{u0})}{T_{u0}(\theta _{T} - T_{u})}&\hbox {\,\,\char 054\,\,}e^{\alpha t}\\ T_{u}(\theta _{T} - T_{u0})&\hbox {\,\,\char 054\,\,}e^{\alpha t} T_{u0}(\theta _{T} - T_{u})\\&= e^{\alpha t}T_{u0}\theta _{T} - e^{\alpha t}T_{u}T_{u0}\\ T_{u}(\theta _{T} - T_{u0}) + e^{\alpha t}T_{u}T_{u0}&\hbox {\,\,\char 054\,\,}e^{\alpha t}T_{u0}\theta _{T}\\ T_{u}((\theta _{T} - T_{u0}) + e^{\alpha t}T_{u0})&\hbox {\,\,\char 054\,\,}e^{\alpha t}T_{u0}\theta _{T}\\ T_{u}&\hbox {\,\,\char 054\,\,}\frac{e^{\alpha t}T_{u0}\theta _{T}}{\frac{e^{\alpha t}}{e^{\alpha t}}\left( \theta _{T} - T_{u0}\right) + e^{\alpha t}T_{u0}}\\&= \frac{T_{u0} \theta _{T}}{\frac{(\theta _{T} - T_{u0})}{e^{\alpha t}} + T_{u0}}. \end{aligned}$$

Taking the limit supremum on both sides, we have

$$\begin{aligned} \limsup \limits _{t\rightarrow \infty } T_{u}&\hbox {\,\,\char 054\,\,}\limsup \limits _{t\rightarrow \infty } \frac{T_{u0} \theta _{T}}{\frac{(\theta _{T} - T_{u0})}{e^{\alpha t}} + T_{u0}}\\&= \frac{T_{u0}\theta _{T}}{T_{u0}}\\&= \theta _{T} , \ \text {which is the upper bound for } T_{u}. \end{aligned}$$

Now, considering the infected tumor cell population, rewriting Eq. (3.1.2) leads to

$$\begin{aligned} \frac{\mathrm{d}T_{i}}{\mathrm{d}t} + \delta T_{i} = \beta T_{u} V, \end{aligned}$$

which is now in first order linear standard form and its solution is given by

$$\begin{aligned} T_{i} = \frac{1}{\delta } \beta T_{u} V \left( 1 - e^{-\delta t}\right) . \end{aligned}$$

Taking limits on both sides of the above equation,

$$\begin{aligned} \limsup \limits _{t\rightarrow \infty } T_{i}&= \lim _{t\rightarrow \infty } \frac{1}{\delta } \beta T_{u} V \left( 1 - e^{-\delta t}\right) \\&\hbox {\,\,\char 054\,\,}\lim _{t\rightarrow \infty } \frac{1}{\delta } \beta T_{u} V\\&= \frac{1}{\delta } \beta T_{u} V, \ \text {Which is the upper bound fot } T_{i}. \end{aligned}$$

For the virus population, Eq. (3.1.3) can be rewritten as

$$\begin{aligned} \frac{\mathrm{d}V}{\mathrm{d}t} + \gamma V = \mu _{v}(t) + b T_{i}, \end{aligned}$$

which gives rise to the solution

$$\begin{aligned} V = \frac{1}{\gamma }(\mu _{v}(t) + b T_{i})\left( 1 - e^{-\gamma t}\right) . \end{aligned}$$

Taking limits on both sides of the above equation,

$$\begin{aligned} \limsup \limits _{t\rightarrow \infty } V&= \lim _{t\rightarrow \infty } \frac{1}{\gamma }(\mu _{v}(t) + b T_{i})\left( 1 - e^{-\gamma t}\right) \\&\hbox {\,\,\char 054\,\,}\lim _{t\rightarrow \infty } \frac{1}{\gamma }(\mu _{v}(t) + b T_{i})\\&= \frac{1}{\gamma }(\mu _{v}(t) + b T_{i}), \ \text {which is the upper bound for } V. \end{aligned}$$

We can therefore conclude that the system above is bounded and the invariant region is given by

$$\begin{aligned} \varvec{\Omega }_{\mathbf {V}} \!=\! \left\{ (T_{u}, T_{i}, V) \in \mathbb {R}_{+}^3 \mid \ 0 \!\hbox {\,\,\char 054\,\,}\! T_{u} \!\hbox {\,\,\char 054\,\,}\! \theta _{T},\ 0 \!\hbox {\,\,\char 054\,\,}\! T_{i}\!\hbox {\,\,\char 054\,\,}\! \frac{1}{\delta } \beta T_{u} V,\ 0 \!\hbox {\,\,\char 054\,\,}\! V \!\hbox {\,\,\char 054\,\,}\! \frac{1}{\gamma }(\mu _{v}(t) \!+\! b T_{i})\right\} . \end{aligned}$$

Note that $\varvec{\Omega }_{\mathbf{V}}$ is a positively invariant domain for the system (Eqs. (3.1.1)–(3.1.3)). Moreover, notice also that $0 \le T_{u}(t) + T_{i}(t) \le \theta _{T}$ for $t > 0$. This shows that the total tumor burden $T_{u}(t) + T_{i}(t)$ cannot exceed the carrying capacity $\theta _{T}$. This proof is analogous to the proof in Appendix A of Dingli et al. (2006). $\varvec{\Omega }_{\mathbf{V}}$ is also a biologically feasible region for the state variables. Hence, we shall also regard this region as a “global” domain.

1.3 (iii) Existence and Uniqueness of Solutions

Proof

Since the right-hand side of system is $\mathbf {C}^{1}$ (class of continuously differentiable functions) satisfies the properties of locally Lipschitz functions, then the existence and uniqueness of solutions of the system are ascertained by the Cauchy–Lipschitz theorem (Diekmann et al. 1995; Schatzman 2002).

For the immunocompetent model (see Eqs. (2.1.1)–(2.1.4)), it is not difficult to establish that the solutions will remain nonnegative for $t\ge 0$ since $\frac{\mathrm{d}N}{\mathrm{d}t} \ge 0$ whenever $N(t) = 0$ and $N(t) \ge 0$, Also it can easily be verified that N is bounded by $0 \le N \le \frac{s_N}{d}$. Thus, we define the invariant region for the system (Eqs. (2.1.1)–(2.1.4)) as

$$\begin{aligned}&\varvec{\Omega }_{\mathbf {VN}} = \left\{ (T_{u},T_{i},V,N) \in \mathbb {R}_{+}^4 \mid \ 0 \hbox {\,\,\char 054\,\,}T_{u} \hbox {\,\,\char 054\,\,}\theta _{T},\ 0\hbox {\,\,\char 054\,\,}T_{i} \right. \\&\left. \hbox {\,\,\char 054\,\,}\frac{1}{\delta } \beta T_{u} V,\ 0 \hbox {\,\,\char 054\,\,}V \hbox {\,\,\char 054\,\,}\frac{1}{\gamma }(\mu _{v}(t) + b T_{i}),\ 0 \hbox {\,\,\char 054\,\,}N \hbox {\,\,\char 054\,\,}\frac{s_N}{d}\right\} . \end{aligned}$$

Hence, the existence and uniqueness of solutions theorem for our immunocompetent model are preserved. Note that the right-hand side of our immunocompetent system is $\mathbf {C}^{1}$ on $\mathbb {R}_{+}^4$.

Stability Theorems of the Steady State ${\mathbf {SS}}_{\mathbf{1}}$

1.1 Proof of Theorem 2 (Local Stability)

Proof

We note that $s_{1} = -r <0$ and $s_{2} = \frac{1}{2}(-\omega _{v} - 1 - \sqrt{\left( \omega _{v}-1\right) ^2 + 4\lambda _{v}\rho }) < 0$, for all nonnegative parameter values, and $s_{3} = \frac{1}{2}(-\omega _{v} - 1 + \sqrt{\left( \omega _{v}-1\right) ^2 + 4\lambda _{v}\rho })$ can either be negative, positive and zero. If $\sqrt{\left( \omega _{v}-1\right) ^2 + 4\lambda _{v}\rho } < 1 + \omega _{v}$, then $s_{3} < 0$. Since $1 + \omega _{v}$ is positive, then it follows that $(\omega _{v}-1)^2 + 4\lambda _{v}\rho < (1+\omega _{v})^2$, which is equivalent to $\rho < \frac{\omega _{v}}{\lambda _{v}}$. Thus, when $\rho < \frac{\omega _{v}}{\lambda _{v}}$, all three eigenvalues ($s_{1}, s_{2}$ and $s_{3}$) are negative. Hence, $\mathbf {SS_{1}}$ is locally asymptotically stable. Equivalently, since $P(s)= s^{3} + \left( r + \omega _{v} + 1\right) s^{2} + \left( - \lambda _{v} \rho + r \omega _{v} + r + \omega _{v}\right) s - r \lambda _{v} \rho + r \omega _{v} = s^{3}+a_{2}s^{2}+a_{1}s+a_{0}=0$, where $a_{2}=\omega _{v}+r+1, \ a_{1}=\omega _{v}r+r+\omega _{v}-\rho \lambda _{v}, \ a_{0}=\omega _{v}r-\lambda _{v}\rho r$, we notice that $a_{2}>0$ and that if $\mathcal {R}_{0}<1$ then $a_{1}>0,a_{0}>0$ and $a_{1}a_{2}-a_{0}>0$. Hence, by the Routh–Hurwitz criterion (Zi 2011; Harris 2002), we deduce that $\mathbf {SS_{1}}$ is locally asymptotically stable. Similarly, whenever $\rho > \frac{\omega _{v}}{\lambda _{v}}$, then $\sqrt{\left( \omega _{v}-1\right) ^2 + 4\lambda _{v}\rho } > 1+\omega _{v}$. This implies that $s_{3} > 0$. Hence, $\mathbf {SS_{1}}$ is unstable.

1.2 Proof of Theorem 3 (Global Stability)

Before we provide a detailed proof of this theorem, note that given the positively invariant domain $\varvec{\Omega }_{\mathbf{V1}}$, the steady-state solution $\mathbf {SS_{1}}$ should be globally asymptotically stable in the whole domain $\varvec{\Omega }_{\mathbf{V}}$. To prove this, we first construct a Lyapunov function based on a specified range of the parameter $\rho$. We simplify the model by the translating the state variables as follows: $T_{u} = 1 - \hat{T}_{u},\ T_{i} = \hat{T}_{i}$, and $V = \hat{V}$. For notation convenience, we drop the hats over the state variables and have the following system:

$$\begin{aligned} \frac{\mathrm{d}T_{u}}{\mathrm{d}t}&= -rT_{u} + rT_{i} + rT_{u}^2 + \lambda _{v}V - rT_{u}T_{i} - \lambda _{v}T_{u}V \end{aligned}$$

(B.2.1)

$$\begin{aligned} \frac{\mathrm{d}T_{i}}{\mathrm{d}t}&= \lambda _{v}V - \lambda _{v}T_{u}V - T_{i} \end{aligned}$$

(B.2.2)

$$\begin{aligned} \frac{\mathrm{d}V}{\mathrm{d}t}&= \rho T_{i} - \omega _{v}V, \end{aligned}$$

(B.2.3)

and the invariant domain $\varvec{\Omega }_{\mathbf{V1}}$ is translated to

$$\begin{aligned} \varvec{\Omega }_{\mathbf{V2}} = \left\{ (T_{u}, T_{i}, V) \in \mathbb {R}_{+}^3 \mid \ T_{u} \ge 0,\ T_{i} \ge 0,\ 0 \hbox {\,\,\char 054\,\,}T_{u} - T_{i} \hbox {\,\,\char 054\,\,}1,\ V \ge 0\right\} . \end{aligned}$$

(B.2.4)

Now we give the proof of Theorem 3:

Proof

Given the nonnegative initial conditions $(T_{u0},T_{i0},V_{0})$ in $\varvec{\Omega }_{\mathbf{V2}}$, then by Theorem A.1 and A.3, the corresponding solutions of the dimensionless state variable satisfy $0 \le T_{u}(t) \le 1,\ 0 \le T_{i}(t) \le 1,$ and $V(t) \ge 0$. It suffices to show that if $T_{i}(t)$ and V(t) approach zero, then $T_{u}(t)$ also approaches zero. When $0< \rho < 1$, we define a Lyapunov function

$$\begin{aligned} G\left( T_{u},T_{i},V\right) = \frac{1}{2}\lambda _{v}\rho \omega _{v}T_{i}^2 + \lambda _{v}^2\rho T_{i}V + \frac{1}{2}\lambda _{v}^2 V^2 \end{aligned}$$

Using Theorem A.1, it is easy to check that $G\left( T_{u},T_{i},V\right) > 0$. The orbital derivative is given by

$$\begin{aligned} \dot{G}\left( T_{u},T_{i},V\right)&= \lambda _{v}\rho \omega _{v} T_{i}\dot{T}_{i} + \lambda _{v}^2\rho \dot{T}_{i}V + \lambda _{v}^2\rho \dot{V}T_{i} + \lambda _{v}^2 V\dot{V} \\&= \lambda _{v}\rho \omega _{v} T_{i}\left( \lambda _{v}V - \lambda _{v}T_{u}V - T_{i}\right) + \lambda _{v}^2\rho \left( \lambda _{v}V - \lambda _{v}T_{u}V - T_{i}\right) V \\&\ \ \ + \lambda _{v}^2\rho \left( \rho T_{i} - \omega _{v}V\right) T_{i} + \lambda _{v}^2 V\left( \rho T_{i} - \omega _{v}V\right) \\&= \lambda _{v}\rho \left( \lambda _{v}\rho - \omega _{v}\right) T_{i}^2 - \lambda _{v}^2\rho \omega _{v} T_{i}T_{u}V + \lambda _{v}^2\left( \lambda _{v}\rho - \omega _{v}\right) V^2 - \lambda _{v}^3\rho T_{u}V^2 \\&\ \ \ + \lambda _{v}\omega _{v}\left( \rho - 1\right) T_{i}V. \end{aligned}$$

Since $\rho < \frac{\omega _{v}}{\lambda _{v}}$, which means, $\lambda _{v}\rho - \omega _{v} < 0$, and $0< \rho < 1$, implies that $\rho - 1 < 0$, then it follows that $\dot{G}\left( T_{u},T_{i},V\right) < 0$. Therefore, using this Lyapunov function, we have $T_{i}(t) \rightarrow 0$ and $V(t) \rightarrow 0$ as $t \rightarrow +\infty$ when $\rho < \frac{\omega _{v}}{\lambda _{v}}$.

Now, we show that $T_{u}(t) \rightarrow 0$ as $t \rightarrow +\infty$ too. From Eq. (B.2.1), we have

$$\begin{aligned} \frac{\mathrm{d}T_{u}}{\mathrm{d}t}&= -rT_{u} + rT_{i} + rT_{u}^2 + \lambda _{v}V - rT_{u}T_{i} - \lambda _{v}T_{u}V \\&= -rT_{u}(1-T_{u}) + rT_{i}(1-T_{u}) + \lambda _{v}V(1-T_{u}) \\&= (1-T_{u})\left[ rT_{i} + \lambda _{v}V - rT_{u}\right] \\&\le rT_{i} + \lambda _{v}V - rT_{u}. \end{aligned}$$

Solving for $T_{u}$, we obtain

$$\begin{aligned} 0 \le T_{u}(t) \le T_{u0}e^{-rt} + e^{-rt} \int _0^t \left( rT_{i}(s) + \lambda _{v}V(s)\right) e^{rs} ds. \end{aligned}$$

Taking the limit on both sides, we have

$$\begin{aligned} \lim _{t\rightarrow \infty } T_{u}(t)&\le T_{u0} e^{-rt} + \lim _{t\rightarrow \infty } e^{-rt} \int _{0}^{t} \left( rT_{i}(s) + \lambda _{v}V(s)\right) e^{rs} ds \\&\le 0 + \lim _{t\rightarrow \infty } \frac{\int _{0}^{t} \left( rT_{i}(s) + \lambda _{v}V(s)\right) e^{rs} ds}{e^{rt}} \\&\le \lim _{t\rightarrow \infty } \frac{1}{r} \frac{\left( rT_{i}(t)+\lambda _{v}V(t)\right) e^{rt}}{e^{rt}} \\&\le \frac{1}{r}\left[ \lim _{t\rightarrow \infty } rT_{i} + \lim _{t\rightarrow \infty } \lambda _{v}V(t)\right] \\&\le \frac{1}{r}(0+0)\ \ \ (\text {since both } T_{u}(t) \rightarrow 0 \text { and } V(t) \rightarrow 0 \text { as } t \rightarrow +\infty .) \\&\le 0. \end{aligned}$$

Thus, we notice that $T_{u}(t) \rightarrow 0$ as $t \rightarrow +\infty$. Hence, given an initial condition in the domain $\varvec{\Omega }_{\mathbf{V2}}$, then $T_{u}(t),\ T_{i}(t)$ and V(t) all approach the origin as $t \rightarrow +\infty$. Therefore, $\mathbf {SS_{1}}$ is a global attractor for the system (3.2.6-3.2.8).

1.3 Proof of Theorem 4 (Threshold Local Stability)

Proof

We use a center manifold theorem to reduce the system (B.2.1–B.2.3) to its local center manifold. First, we segregate the system into two parts, the part with zero eigenvalue and the one with negative eigenvalues. Considering the linear part of system (B.2.1–B.2.3), then the corresponding matrix is

$$\begin{aligned} L = \left( \begin{array}{rrr} -r &{} r &{} \lambda _{v} \\ 0 &{} -1 &{} \lambda _{v} \\ 0 &{} \rho &{} -\rho \lambda _{v} \end{array}\right) \end{aligned}$$

The eigenvalues and corresponding eigenvectors of L are

$$\begin{aligned} \lambda _{1}&= -r, \ \ \ \text {with} \ \ \ V_{1} = \left( 1,0,0\right) ^{\text {T}} \\ \lambda _{2}&= -\left( 1+\rho \lambda _{v}\right) , \ \ \ \text {with} \ \ \ V_{2} = \left( \rho \lambda _{v}-r,1+\rho \lambda _{v}-r,\rho (1+\rho \lambda _{v}-r)\right) ^{\text {T}} \\ \lambda _{3}&= 0, \ \ \ \text {with} \ \ \ V_{3} = \left( \lambda _{v}(r+1),r\lambda _{v},r\right) ^{\text {T}}. \end{aligned}$$

Let $T = \left( V_{1},V_{2},V_{3}\right)$ be the transformation matrix and let $Y = \left( T_{u},T_{i},V\right) ^{\text {T}}$, then we can write the system (B.2.1–B.2.3) into a reduced system

$$\begin{aligned} \frac{\mathrm{d}Y}{\mathrm{d}t} = LY + N, \end{aligned}$$

(B.3.1)

where $N = \left( rT_{u}^2-rT_{u}T_{i}-\lambda _{v}T_{u}V,- \lambda _{v}T_{u}V,0\right) ^{\text {T}}$. Now let $Y = TX$, then we have

$$\begin{aligned} \frac{\mathrm{d}X}{\mathrm{d}t} = T^{-1}LTX + T^{-1}N, \end{aligned}$$

(B.3.2)

where the diagonal matrix $T^{-1}LT$ is define as

$$\begin{aligned} T^{-1}LT = \left( \begin{array}{rrr} -r &{} 0 &{} 0 \\ 0 &{} -(1+r\lambda _{v}) &{} 0 \\ 0 &{} 0 &{} 0 \end{array}\right) , \end{aligned}$$

and $T_{u} = x_{1} + (\rho \lambda _{v}-r)x_{2} + \lambda _{v}(r+1)x_{3}$, $T_{i} = (1+\rho \lambda _{v}-r)x_{2} + r\lambda _{v}x_{3}$, and

$V = -\rho (1+\rho \lambda _{v}-r)x_{2} + rx_{3}$. Now considering the last term, $T^{-1}N$, in Eq. (B.3.2), and let $T^{-1}N = \left( n_{1},n_{2},n_{3}\right) ^{\text {T}}$, then expressing $n_{i},\ i=1,2,3$ in terms of $x_{i}$, we have

$$\begin{aligned} n_{1}&= rT_{u}^{2} - rT_{i} T_{u} - \lambda _{v}T_{u} V \\&= A_{11}x_{1}^{2} + A_{12}x_{1}x_{2} + A_{13}x_{1}x_{3} + A_{22}x_{2}^{2} + A_{23}x_{2}x_{3} + A_{33}x_{3}^{2}. \\ n_{2}&= \frac{1}{\lambda _{v}^{2} \rho ^{2} - \lambda _{v} r \rho + r - 1}\left( \left( \lambda _{v}T_{u} V + r T_{i} T_{u} - rT_{u}^{2} \right) \left( \left( \lambda _{v}^{2} r + \lambda _{v}^{2}\right) \rho ^{2} + r^{2}\right. \right. \\&\ \ \ \ \ - \left. \left. \left( \lambda _{v} r^{2} + \lambda _{v} r - \lambda _{v}\right) \rho \right) + r\lambda _{v} T_{u} V \right) \\&= B_{11}x_{1}^{2} + B_{12}x_{1}x_{2} + B_{13}x_{1}x_{3} + B_{22}x_{2}^{2} + B_{23}x_{2}x_{3} + B_{33}x_{3}^{2}. \\ n_{3}&= \frac{1}{\lambda _{v}^{2} \rho ^{2} - \lambda _{v} r \rho + r - 1}\left( {\left( \lambda _{v} T_{u} V + r T_{i} T_{u} - r T_{u}^{2}\right) } {\left( \lambda _{v}^{2} \rho + \lambda _{v}\right) } - r \lambda _{v}^{2} T_{u} V \right) \\&= C_{11}x_{1}^{2} + C_{12}x_{1}x_{2} + C_{13}x_{1}x_{3} + C_{22}x_{2}^{2} + C_{23}x_{2}x_{3} + C_{33}x_{3}^{2}. \end{aligned}$$

where the coefficients $A_{ij}, B_{ij}$ and $C_{ij},\ i,j = 1,2,3$ can be easily determined. Then, the transformed system can now be written as

$$\begin{aligned} \frac{\mathrm{d}Z}{\mathrm{d}t}&= BZ + \left( \begin{array}{r} n_{1} \\ n_{2} \end{array}\right) \end{aligned}$$

(B.3.3)

$$\begin{aligned} \frac{\mathrm{d}x_{3}}{\mathrm{d}t}&= Ax_{3} + n_{3}. \end{aligned}$$

(B.3.4)

where

$$\begin{aligned} B = \left( \begin{array}{rr} -r &{} 0 \\ 0 &{} -(1+r\lambda _{v}) \end{array}\right) ,\ \ \ \ \ \ \text {and} \ \ \ \ \ A = \left( 0\right) . \end{aligned}$$

It is easy to verify that the matrix B has negative eigenvalues and A has zero eigenvalue. It can also be easy checked that each $n_{i},\ i = 1, 2, 3$, is a $C^{2}$ differentiable function, $n_{i}(0,0,0) = 0$ and $Dn_{i}(0,0,0) = 0$, where $Dn_{i}$ is the variational matrix of the function $n_{i}$. Then, by the Center Manifold Theorem (Carr 1981), there exists a center manifold given by

$$\begin{aligned} Z = h(x_{3}) = \left( \begin{array}{r} h_{1}(x_{3}) \\ h_{2}(x_{3}) \end{array}\right) \end{aligned}$$

(B.3.5)

with $h(0) = 0$ and $Dh(0) = 0$, and it satisfies the equation

$$\begin{aligned} Bh(x_{3}) + \left( \begin{array}{r} n_{1}\left( h(x_{3}),x_{3}\right) \\ n_{2}\left( h(x_{3}),x_{3}\right) \end{array}\right)&= Dh(x_{3}) \cdot \left( Ax_{3} + n_{3}\left( h(x_{3}),x_{3}\right) \right) \end{aligned}$$

(B.3.6)

$$\begin{aligned}&= Dh(x_{3}) \cdot n_{3}\left( h(x_{3}),x_{3}\right) . \end{aligned}$$

(B.3.7)

Let $u = x_{3}$, then we can approximate h(u) as follows:

$$\begin{aligned} h(u)&= \left( \begin{array}{r} h_{1}(u) \\ h_{2}(u) \end{array}\right) \end{aligned}$$

(B.3.8)

$$\begin{aligned}&= \left( \begin{array}{r} c_{2}u^{2} + c_{3}u^{3} + c_{4}u^{4} + O\left( u^{5}\right) \\ d_{2}u^{2} + d_{3}u^{3} + d_{4}u^{4} + O\left( u^{5}\right) \end{array}\right) . \end{aligned}$$

(B.3.9)

For simplicity, we consider the order up to 5, and we will know later if it is enough. Then, we calculate the $n_{i},\ i = 1, 2, 3$, as follows:

$$\begin{aligned} n_{1}(h(u),u)&= n_{1}\left( h_{1}(u),h_{2}(u),u\right) \\&= A_{33}u^{2} + O\left( u^{4}\right) \\ n_{2}(h(u),u)&= n_{2}\left( h_{1}(u),h_{2}(u),u\right) \\&= B_{33}u^{2} + O\left( u^{4}\right) \\ n_{3}(h(u),u)&= n_{3}\left( h_{1}(u),h_{2}(u),u\right) \\&= C_{33}u^{2} + O\left( u^{4}\right) . \end{aligned}$$

We substituting $n_{i},\ i = 1, 2, 3$, into Eq. (B.3.7), and comparing the coefficients on both sides of the equation, we obtain $C_{33} = \frac{-r\lambda _{v}^2(r+1)^2(\lambda _{v}\rho + 1)}{\lambda _{v}^{2} \rho ^{2} - \lambda _{v} r \rho + r - 1} < 0$. Now reducing the system (B.2.1–B.2.3) to its local center manifold, which is a single equation, we have

$$\begin{aligned} \frac{\mathrm{d}x_{3}}{\mathrm{d}t}&= n_{3}(h(x_{3}),x_{3}) = C_{33}x_{3}^{2} + O\left( x_{3}^{4}\right) . \end{aligned}$$

(B.3.10)

Eq. (B.3.10) governs the stability of the zero solution of the system (B.3.3–B.3.4). Since $C_{33} < 0$, then the zero solution, $x_{3} = 0$, is locally asymptotically stable. Therefore, we conclude that the trivial solution of the system (B.2.1–B.2.3) is locally asymptotically stable when $\rho = \rho _{c} = \frac{\omega _{v}}{\lambda _{v}}$.

1.4 Proof of the Coexistence Steady State ($\mathbf {SS_{2}}$) Stability

Proof

If $\mathcal {R}_{0}>1$, then we know that there is an additional steady state (Coexistence Steady State ($\mathbf {SS_{2}}$)) given by

$$\begin{aligned} \mathbf {SS_{2}}=\left( {\frac{\omega _{v}}{\rho \lambda _{v}},}{\frac{\omega _{v}r\left( \lambda _{v}\rho -{\omega _{v}}\right) }{\lambda _{v}\rho \left( \omega _{v}r+{\lambda _{v}\rho }\right) },}{\frac{r\left( \lambda _{v}\rho -\omega _{v}\right) }{\lambda _{v}\left( \omega _{v}r+{\lambda _{v}}\rho \right) }}\right) . \end{aligned}$$

Furthermore, the characteristic equation at $\mathbf {SS_{2}}$ is

$$\begin{aligned} {z}^{3}+b_{2}{z}^{2}+b_{1}z+b_{0}=0 \end{aligned}$$

with

$$\begin{aligned} \begin{array}{l} b_{2}={\frac{\lambda _{v}\rho \omega _{v}+\omega _{v}r+\rho \lambda _{v}}{ \rho \lambda _{v}}} \\ b_{1}={\frac{\omega _{v}r\left( \lambda _{v}\rho r+{\omega _{v}} ^{2}r+\lambda _{v}\rho \omega _{v}+\rho \lambda _{v}\right) }{\rho \lambda _{v}\left( \omega _{v}r+\rho \lambda _{v}\right) }}{\frac{\omega _{v}r\left( \rho \lambda _{v}-\omega _{v}\right) }{\rho \lambda _{v}}} \\ b_{0}={\frac{\omega _{v}r\left( r\rho \lambda _{v}{\omega _{v}}^{3}-{\ \rho } ^{3}{\lambda _{v}}^{3}+{\rho }^{2}{\lambda _{v}}^{2}{\omega _{v}} ^{2}+\lambda _{v}\omega _{v}{r}^{2}\rho +{r}^{2}{\omega _{v}}^{3}+{\lambda _{v}}^{2}{\rho }^{2}r+3r\rho \lambda _{v}{\omega _{v}}^{2}+3{\rho }^{2}{ \lambda _{v}}^{2}\omega _{v}+\lambda _{v}\omega _{v}r\rho +{\ \lambda _{v}} ^{2}{\rho }^{2}\right) }{{\rho }^{2}{\lambda _{v}}^{2}\left( \omega _{v}r+\rho \lambda _{v}\right) }}. \end{array} \end{aligned}$$

Moreover, we have $b_{2}>0$ and if $\mathcal {R}_{0}>1$, then $b_{1}>0,b_{0}>0$ and

$$\begin{aligned} b_{1}b_{2}-b_{0}=&\Bigg (\frac{\omega _{v}r\left( r\rho \lambda _{v}{\omega _{v}} ^{3}-{\rho }^{3}{\lambda _{v}}^{3}+{\rho }^{2}{\lambda _{v}}^{2}{\omega _{v}} ^{2}+\lambda _{v}\omega _{v}{r}^{2}\rho +{r}^{2}{\omega _{v}}^{3}\right. }{{\rho }^{2}{\lambda _{v}}^{2}\left( \omega _{v}r+\rho \lambda _{v}\right) } \\&+ \frac{\left. {\lambda _{v}}^{2}{\rho }^{2}r+3r\rho \lambda _{v}{\omega _{v}}^{2}+3{\rho }^{2}{ \lambda _{v}}^{2}\omega _{v}+\lambda _{v}\omega _{v}r\rho +{\ \lambda _{v}} ^{2}{\rho }^{2}\right) }{{\rho }^{2}{\lambda _{v}}^{2}\left( \omega _{v}r+\rho \lambda _{v}\right) }\Bigg ) >0 \end{aligned}$$

Hence, by Routh–Hurwitz criterion (Zi 2011; Harris 2002), we deduce that $\mathbf {SS_{2}}$ is locally asymptotically stable.

1.5 Proof of Theorem 5

Proof

The equilibrium points of (3.5.5) are given by the solutions of

$$\begin{aligned} \left\{ \begin{array}{l} \left( r\left( 1-T_{u}\right) -c_{T}N\right) T_{u}=0 \\ f_{N}-\eta T_{u}N-\kappa _{N}N=0. \end{array} \right. \end{aligned}$$

(B.5.1)

From (B.5.1)$_{1}$ we obtain $T_{u}=0$ or $T_{u}={\frac{r-c_{T}N}{r}.}$

• 1. If $T_{u}=0$, then we obtain the equilibrium point, $E_{00}=\left( 0, \frac{f_{N}}{\kappa _{N}}\right) .$ The eigenvalues of the variational matrix at $E_{00}$ are $-\kappa _{N}$ and $r-\frac{f_{N}}{\kappa _{N}}c_{T}.$ Hence, $E_{00}$ is locally asymptotically stable if $r<{\dfrac{c_{T}f_{N}}{\kappa _{N}}}$.

• 2. If $T_{u}={\dfrac{r-c_{T}N}{r},}$ then $r>{\dfrac{c_{T}f_{N}}{\kappa _{N}}}$, then from (B.5.1)$_{2}$ we obtain

  $$\begin{aligned} \frac{\eta c_{T}}{r}N^{2}-\left( \eta +\kappa _{N}\right) N+f_{N}=0. \end{aligned}$$

  (B.5.2)

  The discriminant of (B.5.2) is $\Delta =\left( \eta +\kappa _{N}\right) ^{2}-4\dfrac{\eta c_{T}f_{N}}{r}>(\eta -\kappa _{N})^{2}$ (since $r>{\frac{ c_{T}f_{N}}{\kappa _{N}}}$). Hence, (B.5.2) has two real roots $N_{\pm }={\dfrac{r\left( \eta +\kappa _{N}\pm \sqrt{\Delta }\right) }{2\eta c_{T}}}$ which are positive (since their product and sum positive). Moreover, we have $T_{u_{\pm }}={\frac{r-c_{T}N_{\pm }}{r}=\frac{\eta -\kappa _{N}\mp \sqrt{ \Delta }}{2\eta }}$ which, by $\Delta >(\eta -\kappa _{N})^{2},$ implies that $T_{u_{+}}<0$ and $T_{u_{-}}>0.$ This leads to the following virus-free equilibrium point

  $$\begin{aligned} E_{0}=\left( T_{u_{-}},N_{_{-}}\right) . \end{aligned}$$

  The eigenvalues of the variational matrix at $E_{0}$ are given by the roots of the following characteristic equation

  $$\begin{aligned} z^{2}+\left( \left( r+\eta \right) T_{u_{-}}+\kappa _{N}\right) z+r\sqrt{ \Delta }T_{u_{-}}=0. \end{aligned}$$

  (B.5.3)

  Hence, $r>{\dfrac{c_{T}f_{N}}{\kappa _{N}}}$, then $T_{u_{-}}>0$ implying that (B.5.3) does not have roots with negative real parts.

1.6 The Basic Reproductive Number of the Virus-Free Sub-model

We use the next-generation method described in den Driessche and Watmough (2002) to calculate the basic reproductive number at $E_{0}$. The vector formed by the rates of new infections is given by

$$\begin{aligned} \mathcal {F}=\left[ \begin{array}{c} \lambda _{v}T_{u}V \\ 0 \end{array} \right] \end{aligned}$$

The vector formed by the other transfer rates is

$$\begin{aligned} \mathcal {W}=\left[ \begin{array}{c} T_{i}+c_{TV}T_{i}N \\ -\rho T_{i}+\omega _{v}V+c_{V}VN \end{array} \right] \end{aligned}$$

The next-generation matrix is given by $M=FV^{-1},$ where $F=D\mathcal {F}_{ E_{0}}=\left[ \begin{array}{cc} 0 &{} \lambda _{v}T_{u0} \\ 0 &{} 0 \end{array} \right]$ and

$W=D\mathcal {W}_{E_{0}}=\left[ \begin{array}{cc} 1+c_{TV}N_{0} &{} 0 \\ -\rho &{} \omega _{v}+c_{V}N_{0} \end{array} \right] .$ We obtain

$$\begin{aligned} M=\left[ \begin{array}{cc} \frac{\rho \lambda _{v}\left( r-c_{T}N_{0}\right) }{r\left( c_{TV}N_{0}+1\right) \left( \omega _{v}+c_{V}N_{0}\right) } &{} \frac{\lambda _{v}\left( r-c_{T}N_{0}\right) }{r\left( \omega _{v}+c_{V}N_{0}\right) } \\ 0 &{} 0 \end{array} \right] . \end{aligned}$$

Thus,

$$\begin{aligned} \mathcal {R}_{0}=\frac{\rho \lambda _{v}\left( r-c_{T}N_{0}\right) }{r\left( c_{TV}N_{0}+1\right) \left( \omega _{v}+c_{V}N_{0}\right) }. \end{aligned}$$

Model Parameters and Initial Conditions

Known model parameter values and their ranges were taken from available literature, while unknown parameter ranges were estimated using values that seemed biologically reasonable. The baseline parameters used in the model simulations are given in Table 1. Since various mechanisms might lead to NK cell activation and recruitment, we parametrize our model based on the fact that NK cell activation/recruitment is dependent on interactions between the NK cells and infected tumor cells (Leung et al. 2020) or an ICD of infected cells (Bommareddy et al. 2019; Somma et al. 2019). Thus, we adopt the activation/recruitment rate for NK cells as $\xi _{N} = 1\times 10^{-5}$ day$^{-1}$ (Guo et al. 2019) as our baseline value, but we shall vary this parameter in the numerical analysis to explore the potential effects of induced NK cell-mediated responses in oncolytic virotherapy. For the other parameters which govern the cytopathicity of OV, we also sweep parameters within their realistic biological ranges. Note that while we fix the values of most of the parameters used in our simulations, small variations in the values of OV or NK cell related parameters allow our model to capture different dynamics of the OV–tumor–NK cell interactions. Furthermore, in this study we do not consider a specific virus as our baseline state variable because we are interested in the OV-induced NK cell-mediated responses in general. We explore a plausible range of recruitment rates to capture the range of likely behavior of slow and fast replicating viruses. Notably, our model captures different aspects of slow replicating viruses, such the recombinant measles virus (MV-I98A-NIS) (Kemler et al. 2019; Ennis et al. 2010), which may lead to slow recruitment of activated NK cell, and the fast replicating viruses, such as the oncolytic herpes simplex virus (HSV) (Alvarez-Breckenridge et al. 2012b), which may lead to rapid recruitment of activated NK cells (Alvarez-Breckenridge et al. 2012b). It is also known, however, that the same oncolytic virus might replicate differently in different regions (Kemler et al. 2019). Hence, we do not consider specific virus kinetics in this study, but rather focus on the ability of OV infection in recruiting activated NK cells. For the sake of parameter estimation, we consider the virus kinetics of oncolytic vesicular stomatitis virus (VSV), which is also a fast replicating virus (Mahasa et al. 2017; Eftimie and Eftimie 2019), to estimate the virus lysis rate, $\delta$, and the viral clearance rate, $\gamma$. Tumor- and NK cell-related kinetic parameters are taken from the experimental–mathematical models available in the literature and their sources are provided in Table 1. Unless otherwise stated, we assume the following initial conditions: $T_{u} = 1\times 10^{6}$ cells, $T_{i} = 0$ cells and $V = 1\times 10^{6}$ PFU.