We formulate a tumor-immune interaction model with a constant delay to capture the behavior following application of a dendritic cell therapy. The model is validated using experimental data from melanoma-induced mice. Through theoretical and numerical analyses, the model is shown to produce rich dynamics, such as a Hopf bifurcation and bistability. We provide thresholds for tumor existence and, in a special case, tumor elimination. Our work indicates a sensitivity in model outcomes to the immune response time. We provide a stability analysis for the high tumor equilibrium. For small delays in response, the tumor and immune system coexist at a low level. Large delays give rise to fatally high tumor levels. Our computational and theoretical work reveals that there exists an intermediate region of delay that generates complex oscillatory, even chaotic, behavior. The model then reflects uncertainty in treatment outcomes for varying initial tumor burdens, as well as tumor dormancy followed by uncontrolled growth to a lethal size, a phenomenon seen *in vivo*. Theoretical and computational analyses suggest efficacious treatments to use in conjunction with the dendritic cell vaccine. Additional analysis of a highly aggressive tumor additionally confirms the importance of representation with a time delay, as periodic solutions are strictly able to be generated when a delay is present.

Delay differential equations (DDEs) are routinely used to reflect intricate behaviors caused by naturally occurring delays. We propose a set of differential equations to examine the tumor-immune dynamics following application of a dendritic cell therapy. Delay is used to reflect the immune system response time. This model provides a good fit to clinical data and exhibits naturally observed complexities, such as oscillatory behavior and bistability. Sensitivity in model outcomes to initial conditions is reflected and dependent on the immune system response time. The model additionally reveals regions of complex, chaotic behavior and establishes thresholds for tumor elimination and existence.

## I. INTRODUCTION

Since their landmark discovery in 1973 by Steinman and Cohn,^{1} dendritic cells (DCs) have been hailed as the most potent antigen-presenting cells with the singularly efficient ability to initiate and coordinate immune responses. Their critical role in linking the adaptive and innate immune responses has made DCs an attractive candidate in cancer treatments. The first clinical trial testing a DC-based immunotherapy was published in 1996, in which four patients with follicular B-cell lymphoma were infused with antigen-pulsed DCs.^{2} As a result of patients responding in measurably favorable manners, a substantial number of clinical trials studying DC vaccines have been and continued to be conducted. In this form of immunotherapy, DCs are extracted from the patient and sensitized with tumor-associated antigens. The activated DCs are injected into the patient where they interact with and activate naïve and memory T cells in the lymphoid organs. Following exposure to the tumor-associated antigen, the T cells become activated and differentiate into cytotoxic T lymphocytes (CTLs), also known as effector cells, which migrate to the tumor and mount a fight.

DC vaccination studies have recently resulted in the approval of new government-approved cancer treatments. In 2010, sipuleucel-T (Provenge®) became the first DC-based therapy approved by the FDA.^{3} Sipuleucel-T was shown to yield an increase in overall survival of 4.1 months and an improvement of 8.7% in the 3-year survival probability for men with metastatic castration-resistant prostate cancer. Seven years later, the Indian government agency (CDSCO—Central Drugs Standard Control Organization) approved the DC-based vaccine APCEDEN® for treatment of prostate, ovarian, non-small cell lung carcinoma, and colorectal cancer.^{4} The recent federal approval of these DC vaccines rekindled interest in DC research, though excitement waned as the clinical results proved unsatisfactory. Much of the ineffective behavior has been attributed to the immunosuppressive tumor microenvironment.^{5} Interest in enhancing the immune response of DC vaccines has given rise to many clinical trials and mathematical modeling efforts.^{6–11}

Ordinary differential equations (ODEs) are extensively used in mathematical modeling to better understand various life processes, such as cell interactions,^{12,13} disease transmission,^{14,15} and predator–prey relationships.^{16,17} Implicit in ODE models is the assumption that reaction times are instantaneous. As such, these models are often unable to capture much of nature’s true complexity, as they are approximations of reality. Delay differential equations (DDEs) reflect naturally occurring delays, such as the binding time required for cell activation, and thus exhibit more interesting and realistic dynamics.

In Dickman *et al*.,^{18} we formulated and analyzed a system of four ODEs describing the interactions between DCs in the blood, DCs in the tumor, effector cells, and tumor cells as DC therapy was administered. This model system was motivated and based on the earlier work of de Pillis *et al*.^{19} The simplified system in Dickman *et al*.^{18} is given by

where $Db,Dt,E$, and $T$, respectively, denote DCs in the blood, DCs in the tumor, effector cells, and tumor cells. The $\mu $ parameters reflect transport between compartments, while $vb$ and $vt$ reflect intravenous and intratumoral dosing. The meanings of the remainder of the parameters are reflected in Table I.

Para. . | Description . | Value . | Reference . |
---|---|---|---|

D_{i} | Activation rate of immature DCs by tumors | 3.4946 × 10^{−4}/day | Fitted |

δ_{D} | Natural death rate of DCs | 0.34/day; range: [.116,.5] | 18 |

s_{E} | Source of activated CTLs | 2.83 × 10^{−3}cells/day | Fitted |

c | Activation/proliferation rate of CTLs | 0.016 001/day | Fitted |

α | Proportion of DCs in blood | 0.2778 | Fitted |

r_{am} | Inactivation rate of activated CTLs | 0.002/day; range: [4 × 10^{−4} − 1.2] | 18 |

δ_{E} | Natural death rate of activated CTLs | 0.1155/day | 18 |

τ | Time from activation to killing tumor cells | 0.744 912 days | Fitted |

c_{e} | Max inactivation rate of activated CTLs by tumor cells | 5.539 × 10^{−14}/(cells×day); range: [9.42 × 10^{−14} − 10^{−3}] | 18 |

r | Tumor cell growth rate | 0.39/day; range: [0.17 − 0.69] | 18 |

k | Tumor cell carrying capacity | 10^{9} cells | 18 |

c_{t} | Max rate activated CTLs kill tumor cells | 3.5 × 10^{−5}/(cells×day); range: [0 − 1] | 18 |

Para. . | Description . | Value . | Reference . |
---|---|---|---|

D_{i} | Activation rate of immature DCs by tumors | 3.4946 × 10^{−4}/day | Fitted |

δ_{D} | Natural death rate of DCs | 0.34/day; range: [.116,.5] | 18 |

s_{E} | Source of activated CTLs | 2.83 × 10^{−3}cells/day | Fitted |

c | Activation/proliferation rate of CTLs | 0.016 001/day | Fitted |

α | Proportion of DCs in blood | 0.2778 | Fitted |

r_{am} | Inactivation rate of activated CTLs | 0.002/day; range: [4 × 10^{−4} − 1.2] | 18 |

δ_{E} | Natural death rate of activated CTLs | 0.1155/day | 18 |

τ | Time from activation to killing tumor cells | 0.744 912 days | Fitted |

c_{e} | Max inactivation rate of activated CTLs by tumor cells | 5.539 × 10^{−14}/(cells×day); range: [9.42 × 10^{−14} − 10^{−3}] | 18 |

r | Tumor cell growth rate | 0.39/day; range: [0.17 − 0.69] | 18 |

k | Tumor cell carrying capacity | 10^{9} cells | 18 |

c_{t} | Max rate activated CTLs kill tumor cells | 3.5 × 10^{−5}/(cells×day); range: [0 − 1] | 18 |

Though strictly consisting of ODEs, our model was able to produce complex dynamics, including oscillatory behavior and bistability. The existence of bistability suggested that more aggressive treatments were needed than would otherwise be expected, and the mathematical analysis revealed thresholds that guaranteed tumor elimination or existence. Motivated by this work, we seek to gain additional insights by further reducing the system and introducing delay.

The primary aim of this article is to investigate how incorporating the more realistic delay in interactions influences the dynamics of the system. In doing so, the Dickman *et al*.^{18} model is first reduced in two stages. By assuming a constant proportion of DCs in the blood, we collapse the DCs in the blood and tumor into a single equation. Through the introduction of delay, we account for the time it takes for activated CTLs to become effective in killing tumor cells. Numerical experiments, including a justification using clinical data and bifurcation diagrams, reveal that the model remains capable of representing complex tumor-immune behavior. A quasi-steady state assumption is employed to further reduce the model to a system of two DDEs. The model is again justified using clinical data. Following non-dimensionalization, mathematical analysis reveals the existence of periodic solutions and stability switches. Numerical experiments confirm the results of the mathematical analysis and demonstrate the sensitivity to the immune system response time. The limiting case when the immunosuppressive tumor environment is neutralized is analytically explored. Furthermore, the limiting case in which the death rate of effector cells by tumor cell interactions far exceeds the natural death rate is analyzed both numerically and analytically.

This article is organized in the following manner. A description of the intermediate model formulation along with supporting numerical experiments and preliminary analysis is given in Sec. II. The formulation and preliminary analysis of the reduced model are provided in Sec. III. A thorough analysis of the interior equilibria stability and corresponding numerical experiments is provided in Sec. IV. The stability of the boundary equilibria is provided in Sec. V. The special case when the tumor is highly aggressive is explored in Sec. VI. The main results are outlined and discussed in Sec. VIII.

## II. INTERMEDIATE MODEL

In order to better study the effect of a delayed immune response, we make several simplifying assumptions to reduce system (1) and allow for a more mathematically tractable model. Since our focus is no longer on different dosing techniques, there is no longer a need to represent the intratumoral and intravenous dosing as separate terms. We collapse the DCs into a single equation, assuming the proportion of DCs in the blood, and $DbDb+Dt$ is approximated by the constant $\alpha $. Additionally, we assume that the DC therapy takes $\tau $ units of time to become effective; i.e., the effector cells take $\tau $ units of time following activation to kill tumor cells. Thus, the activation of naïve and memory effector cells is assumed to occur at time $t\u2212\tau $. Since activated effector cells die at a rate of $\delta ~E$, the probability of survival for the effector cells is given by $e\u2212\delta ~E\tau $ during this period of delay. The assumptions give rise to the following model:

where $c~=c\alpha $, and we assume that the initial values are

The biological assumptions are justified by the fit to clinical data displayed in Fig. 1, as system (2) is capable of capturing the behavior of the melanoma-induced mice data from Lee *et al*.^{20} with a single set of biologically reasonable parameters given in Table I. Data from Lee *et al*.^{20} reflect the average tumor volume of melanoma-induced mice treated with $0,1\xd7105,7\xd7105,$ or $21\xd7105$ DCs on days 6, 8, and 10. The fit to the data remains comparable to that of Dickman *et al*.^{18} Numerical experiments seen in Fig. 2 indicate that the reduced model continues to give rise to interesting dynamics, such as oscillatory behavior and the presence of a singular Hopf bifurcation, as in the original ODE model given by (1). Singular Hopf bifurcations occur in systems with multiple time scales. In these systems, an equilibrium may cross a fold of the critical manifold. By the definition outlined in Guckenheimer,^{21} a singular Hopf occurs at a distance that is $O(\epsilon )$ from a fold point. This type of bifurcation is correlated with a canard explosion, which can be computationally observed in the form of a dramatic increase in the amplitude of the periodic orbits.

In this section, we confirm that the solutions to system (2) are biologically meaningful through establishing positivity and boundedness given appropriate initial conditions.

## III. PRELIMINARY ANALYSIS: A SIMPLIFIED MODEL

The rapid turnover of DCs plays a critical role in regulating the magnitude of an immune response, with studies finding defects in DC apoptosis giving rise to autoimmune disorders.^{22} Activated effector cells are similarly known to be short-lived, though some have increased lifespans upon differentiating to become memory cells. Recent studies have shown that myeloid DCs, a DC subset largely utilized in DC vaccines, have lower levels of anti-apoptotic BCL-2/BCL-xL and higher levels of pro-apoptotic molecules when compared to T cells.^{23} In particular, the anti-apoptotic BCL-xL is elevated in activated effector cells,^{24} thereby promoting cell survival. As such, the turnover for dendritic cells is quicker than that of the effector cells,^{9,25} and a quasi-steady state assumption can be applied to $D$. The system (2) then reduces to the simpler system given by

where $v~=v\delta D,D~i=Di\delta D$. The simplifications can again be justified through the fit to the Lee *et al*.^{20} data, as shown in Fig. 3, where the fitting is comparable to all previous model formulations.

Non-dimensionalization, taking $t\xaf=rt,E\xaf=Ek,T\xaf=Tk$, mathematically reduces the number of parameters. In the remainder of the paper, we will then consider the dimensionless system

where $\beta =sErk,\gamma =c~v~rk,\delta =\delta ~Er,\rho =c~D~ir,\eta e=cekr,\nu =(1\u2212\alpha )v~k,\mu =(1\u2212\alpha )D~i,\eta t=ctkr,and\tau \xaf=r\tau $ are the dimensionless, positive parameters (with the exception of $\gamma =\nu =0$ when no DC therapy is administered), and we drop the bars from variables and parameters. Values of the dimensionless parameters computed from the parameters of system (4) are given in Table II.

Parameter . | Value . |
---|---|

β | 7.39 × 10^{−12} |

δ | 0.3277 |

ρ | 1.45 × 10^{−5} |

η_{e} | 2.46 × 10^{−4} |

μ | 6.46 × 10^{−4} |

η_{t} | 9138.38 |

τ | 0.382 |

Parameter . | Value . |
---|---|

β | 7.39 × 10^{−12} |

δ | 0.3277 |

ρ | 1.45 × 10^{−5} |

η_{e} | 2.46 × 10^{−4} |

μ | 6.46 × 10^{−4} |

η_{t} | 9138.38 |

τ | 0.382 |

Positivity and boundedness of the dimensionless system (11) can similarly be proven as before. The tumor-free equilibrium is given by $(E0,T0)=((\beta +\gamma e\u2212\delta \tau )/\delta ,0)$. Through the next generation matrix approach,^{26,27} the basic reproduction number can be calculated as

viewing the tumor cells as an infectious disease, as in Dickman *et al*.^{18} Biologically, the basic reproduction can be interpreted as the ratio of the proliferation potential of a tumor cell to the strength of the immune response and crowding effects. We will see that the system (11) admits two positive equilibria when $R0<1$, thus remaining consistent with the existence of a backward bifurcation observed in the Dickman *et al*.^{18} model.

If $R0\u22651,$ then there exists a unique positive equilibrium, $E1=(E\u2217(T1\u2217),T1\u2217)$. For $\nu $ sufficiently small, there exist constants $C3$ and $Rcrit$ such that if $\eta e>C3$ and $Rcrit<R0<1,$ then there exists a positive equilibrium $E2=(E\u2217(T2\u2217),T2\u2217)$ in addition to $E1$.

*et al*.,

^{18}$\eta e$ can be chosen sufficiently large, taking

Figure 4 confirms the analytical findings since both conditions necessary for the existence of two positive equilibria hold. With $\eta e$ chosen sufficiently large, then $A1<0$, and $A12\u22124A0A2>0$ holds when $R0\u2208(Rcrit,1)$.

Figure 5 clearly illustrates how the interior equilibria both depend on $\tau $. For lower $\tau $, $E0$ exists, but neither interior equilibria are feasible. Increasing $\tau $ shifts the E nullcline down and allows for the existence of the interior equilibria $E1$ and $E2$. A further increase to the delay causes the intermediate equilibrium $E2$ to coincide with $E0$ at a finite $\tau $, corresponding to $R0=1$. The condition $Rcrit<R0<1$ required for two equilibria is thus shown to only be satisfied for finite $\tau $. Though $E2$ is not feasible for a higher value of $\tau $, the interior equilibrium $E1$ remains feasible as $\tau \u2192+\u221e$.

Thus, for $\eta e$ sufficiently large such that $\eta e>C3$, there always exists at least one positive equilibrium root when $R0\u2265Rcrit$. Additionally, the second interior equilibrium only exists when the delay $\tau $ and other parameters are chosen such that $Rcrit<R0<1$. When conducting analysis, it will, therefore, be critical to keep track of how the feasibility of interior equilibria changes as $\tau $ or other parameters vary.

## IV. STABILITY OF INTERIOR EQUILIBRIA

### A. No delay

To understand the effects of delay on system (5), we first establish the stability of the equilibria when there is no delay. Without delay, the system is given as follows:

where $\gamma \beta =\gamma +\beta $. Figure 6 depicts the nullclines and equilibria when the conditions $Rcrit<R0<1$ and $\eta e>C3$ are satisfied for the existence of two interior equilibria.

We determine the local stability of the interior equilibria through a geometric approach.

If $Rcrit<R0<1$ and $\eta e>C3$ are satisfied for the existence of $E1$ and $E2$ (the intermediate tumor), the interior equilibria of (11), then $E1$ is stable, and $E2$ is unstable.

More aggressive treatment than standard is then required to decrease $R0$ sufficiently outside the region of bistability. Reducing the strength of the immunosuppressive microenvironment would increase the efficacy of the DC vaccines, as $\eta e$ is derived from the inactivation rate of CTLs by the tumor ($ce$), and sufficient reduction of $\eta e$ eliminates the region of bistability. Agents to inhibit immunosuppressive factors, such as anti-interleukin-10 (anti-IL-10) and anti-transforming growth factor-$\beta $ (anti-TGF-$\beta $), would then be suitable candidates for combination therapies.^{28}

### B. Linear stability of interior equilibria

In this section, we study the linear stability of the interior equilibria, $E1$, $E2$, assuming that the conditions $\eta e>C3$ and $Rcrit<R0<1$ hold such that both interior equilibria are feasible. However, we again note that the condition $Rcrit<R0<1$ can only be satisfied up to a finite value of $\tau $, as $R0$ monotonically increases in $\tau $, thereby causing the intermediate equilibrium $E2$ to lose its feasibility, as depicted in Fig. 5. Hence, it is necessary to keep track of the feasibility of the equilibria during any analysis.

Let us first linearize (5) at the interior equilibrium $(E\u2217,T\u2217)$. Setting $x=E\u2212E\u2217,y=T\u2212T\u2217$, where $x,y$ are small, gives

If we take solutions to be of the form $(x,y)=(c1,c2)e\lambda t$, then non-trivial solutions exist if and only if the characteristic equation $D(\lambda ,\tau )=0$, where

The characteristic equation (13) can alternatively be expressed in the following form:

where $P(\lambda ,\tau )=\lambda 2+A(\tau )\lambda +C(\tau )$ and $Q(\lambda ,\tau )=B(\tau )\lambda +D(\tau )$ with

Note that the characteristic equation (14) involves delay not just in the $e\u2212\lambda \tau $ term, but in several other places as well, since $T\u2217$ and $E\u2217$ depend on $\tau $. Since our model (5) has delay dependent parameters, there is increased complexity in analyzing the system. Analytical results for systems with delay independent parameters have been thoroughly studied^{29} and result in clean, explicit calculations, such as the exact values of $\tau $ where stability switches occur. We apply the approach from Beretta and Kuang^{30} for studying these challenging characteristic equations, which analyzes the stability switches via computational means, as exact values of $\tau $ where these switches occur cannot be analytically found.

In order to apply the technique of Beretta and Kuang,^{30} we must verify that the following properties hold for all $\tau \u2208R+$:

$P(0,\tau )+Q(0,\tau )\u22600,\u2200\tau \u2208R+0$;

If $\lambda =i\omega $, $\omega \u2208R$, then $P(i\omega ,\tau )+Q(i\omega ,\tau )\u22600,\tau \u2208R$;

$lim\u2006sup{|Q|(\lambda ,\tau )/P(\lambda ,\tau )|:|\lambda |\u2192\u221e,Re\lambda \u22650}<1$ for any $\tau $;

$F(\omega ,\tau ):=|P(i\omega ,\tau )|2\u2212|Q(i\omega ,\tau )|2$ for each $\tau $ has a finite number of zeros; and

Each positive root $\omega (\tau )$ of $F(\omega ,\tau )=0$ is continuous and differentiable in $\tau $ whenever it exists.

We first have $P(0,\tau )+Q(0,\tau )=C+D$. By assumption, property (i) is established, implying that $\lambda =0$ is not a root of (13). Now, $P(i\omega ,\tau )+Q(i\omega ,\tau )=\u2212\omega 2+i\omega (A+B)+C+D\u22600$, satisfying property (ii). Clearly, property (iii) holds, as $lim\u2006sup|\lambda |\u2192\u221e|Q(\lambda ,\tau )P(\lambda ,\tau )|=0$. Additionally, we can clearly see that $F(\omega ,\tau ):=|P(i\omega ,\tau )|2\u2212|Q(i\omega ,\tau )|2$ has a maximum of four roots since

thus satisfying property (iv). Let $(\omega 0,\tau 0)$ be an arbitrary point in its domain such that $F(\omega 0,\tau 0)=0$. Then, we see that $F\omega $ and $F\tau $ exist and are continuous in a neighborhood of $(\omega 0,\tau 0)$. Additionally, $F\tau (\omega 0,\tau 0)\u22600$. Then, by the implicit function theorem, we have property (v) holds. Since the five properties hold, we can then apply the results from Beretta and Kuang.^{30}

As the roots of the characteristic equations are functions of delays, a stability switch may occur, where the stability of the equilibrium changes as the delay increases. If a stability switch occurs, we know that (4) has a pair of conjugate pure imaginary roots $\lambda (\tau )=\xb1i\omega (\tau ),\omega (\tau )\u2208R+$, which crosses the imaginary axis at some $\tau =\tau \u2217$, a necessary condition for a Hopf bifurcation to occur.

To find the value of $\tau \u2217$ where a stability switch may occur, without loss of generality, we let $\lambda (\tau )=i\omega (\tau )$. Substitution of $\lambda =i\omega $ into (13) yields the following:

By using Euler’s formula $ei\theta =cos\u2061\theta +isin\u2061\theta $, we simplify and solve for the real and imaginary parts to obtain

Then, $\omega $ must be a solution of

with its roots as

where $\Delta =(2C\u2212A2)2\u22124(C2\u2212D2).$ We know that $\omega \u22122$ is always negative, as

If $\omega +2<0$, a real $\omega $ does not exist, and there are no stability switches. We can see that there exists one real solution $\omega >0$ if and only if

By Theorem 1, $C+D>0$ at $E1$ and $C+D<0$ at $E2$. Since $D\u22650$, we can deduce $C<0$ at $E2$, making $C\u2212D<0$ at $E2$. Thus, a real $\omega $ does not exist, and there are no stability switches for $E2$. It then follows that $E2$ is unstable for all $\tau \u22650$. The sign of $C\u2212D$ proves more challenging to evaluate for $E1$. Therefore, for a stability switch to exist for $E1$, the following additional condition is needed:

Since a real $\omega $ exists for $E1$ when $C\u2212D<0$, the characteristic equation (13) has a pair of simple and conjugate pure imaginary roots $\lambda =\xb1i\omega $. By the application of Theorem 4.1 of Beretta and Kuang,^{30} we prove that this complex conjugate pair of eigenvalues crosses the imaginary axis, ensuring a stability switch for $E1$ when $C\u2212D<0$.

Assume that $I\u2286R+0$ is the set where $\omega (\tau )>0$ is a root of (19). For any $\tau \u2208I$, we use (17) to define the angle $\theta (\tau )\u2208[0,2\pi ]$ as follows:

For $\tau \u2208I$, we have the relationship

From (23), we see

If the conditions hold for an imaginary solution to exist, making $\omega $ feasible for $\tau \u2208I$, then we have the continuous and differentiable sequence of functions $I\u2192R$,

The characteristic equation (13) has a pair of simple conjugate pure imaginary roots $\lambda =\xb1i\omega (\tau \u2217)$ at $\tau \u2217\u2208I$ when $Sn(\tau \u2217)=\tau \u2217\u2212\tau n(\tau \u2217)=0$ for some $n\u2208N0$. By application of Theorem 4.1 in Beretta and Kuang,^{30} the following theorem then holds.

We determine the direction of the roots crossing the imaginary axis by evaluating

We compute

By the characteristic equation (14), we have

Evaluating at $\lambda =i\omega (\tau \u2217)$, we have

As the sign is complex to determine analytically, we computationally determine whether the roots cross the imaginary axis and a stability switch for $E1$ exists as $\tau $ increases, recalling the relationship

by Theorem 2.

We consider two scenarios to demonstrate the impacts of the delay $\tau $ and $\rho $ on $E1$. Recall that $\rho $ is the only parameter in system (5) that is dependent on $c$, the activation/proliferation rate of CTLs, when there is not a constant DC dose applied ($v=0$). The stability of the system is demonstrated when the CTL activation/proliferation rate is low, intermediate, and high. Figures 7, 10, and 12 depict the functions $S0(\tau ),S1(\tau ),$ and $S2(\tau )$ as defined above. For both scenarios, $C\u2212D<0$ up to a finite $\tau $ such that real $\omega $ exists and Theorem 2 can be applied.

Figure 7 illustrates the case of the low activation/proliferation rate of the CTLs ($\rho $ small). The function $S0(\tau )$ is shown to have roots at $\tau 01$ and $\tau 02$, with $\tau 01<\tau 02$. For the choice of parameters, no real roots exist for $Sn(\tau )$ when $n>0$. The behavior of $S0$ indicates that $E1$ is asymptotically stable, loses its stability for $\tau \u2208(\tau 01,\tau 02)$, and then regains its stability as the delay increases. The switches in stability suggested by Fig. 7 are confirmed through bifurcation diagrams in Fig. 8 and phase portraits depicted in Fig. 9. A Hopf bifurcation occurs when $S0$ initially crosses the $\tau $ axis and the characteristic equation about $E1$ admits two roots with a positive real part.

Figure 10 depicts the case of an intermediate activation/proliferation rate of CTLs ($\rho $ intermediate). Similar to Fig. 7, $E1$ loses stability and regains it as $\tau $ increases. However, there additionally exists a region in Fig. 10 such that an additional pair of unstable eigenvalues appears. Thus, when $\tau \u2208(\tau 11,\tau 12)$, there are two pairs of characteristic roots of the characteristic equation about $E1$ with positive real parts, giving a region of increased instability prior to $E1$ eventually regaining its stability. The number of eigenvalues with positive real parts jumps by two as the $Sn$ curves increase across the $\tau $ axis, making $E1$ more unstable, and decreases by two as the $Sn$ curves decrease across the $\tau $ axis, lessening the complexity. Similar to computational results outlined in Gourley and Kuang,^{31} when $\tau $ increases such that $\tau \u2209I$ and real $\omega $ is no longer feasible, the eigenvalues of the characteristic equation become real and negative. Thus, the interior equilibrium $E1$ continues to exist and be asymptotically stable when real $\omega $ is no longer feasible. $E1$ is then asymptotically stable for $\tau \u2208[0,\tau 01)\u222a(\tau 02,\u221e)$ and unstable for $\tau \u2208(\tau 01,\tau 02)$, with added instability for $\tau \u2208(\tau 11,\tau 12)$, as confirmed by the phase plots in Fig. 11.

Figure 12 considers the case of a high activation/proliferation rate of CTLs ($\rho $ large). For $(\tau 21,\tau 22)$, chaotic behavior is produced (in a sense that the characteristic equation about $E1$ admits three pairs of eigenvalues with positive real parts), as displayed in the time series plots in Fig. 13. The regions of chaos and oscillatory behavior reflect how different initial tumor burdens can experience drastically distinct outcomes, as seen in practice.

Eliminating the assumption $C+D\u22600$, there remains one critical case to consider.

#### 1. Critical case: $C+D=0$

In the critical case of $C+D=0$, Theorem 2 is no longer valid, and we must analyze the stability in a different way. When $C+D=0$, it follows that

and $E2$ is only feasible when $\tau =0$ and $E2$ and $E1$ coalesce into a single equilibrium. The characteristic equation (14) reduces to $\lambda 2+A(\tau )\lambda +C(\tau )+D(\tau )=0$ when $\tau =0$. With the roots given by $\lambda =\u2212A(\tau )<0$ and $\lambda =0$, the interior equilibrium $E1=E2$ is stable. When $\tau >0$, application of results from Kuang (Ref. 29, pp. 79, 80) prove that $E1$ loses its stability for $\tau >A(\tau )D(\tau )$.

## V. STABILITY OF THE BOUNDARY EQUILIBRIUM

In the subsequent analysis, we seek to understand the necessary conditions for tumor elimination; thus, we examine the local and global stability of the boundary equilibrium.

The boundary equilibrium $E0$ is locally asymptotically stable when $R0<1$ and unstable when $R0>1$.

The next result establishes the thresholds for the global stability of the tumor-free equilibrium. We only consider the case when the immunosuppressive environment of the tumor is neutralized, corresponding to $\eta e=0$.

Let $\eta e=0$. If $R0<1,$ then $limt\u2192\u221e(E(t),T(t))=((\beta +\gamma e\u2212\delta \tau )/\delta ,0)$.

The immunosuppressive nature of the tumor microenvironment has been attributed to unmet expectations in DC vaccine behavior.^{32} Inhibiting the immunosuppressive factors through a combination treatment allows for less DC vaccines to be necessary for tumor elimination.

## VI. SPECIAL CASE: *δ* = 0

In highly aggressive tumors, the death of the effector cells is dominated by the interaction with the tumor cells. We analyze the equivalent scenario in Secs. VI A and VI B, where $\delta =0$. The system (5) becomes

where $\gamma \beta =\gamma +\beta $.

We first establish stability when $\tau =0$, as it is necessary for understanding the effects of incorporating delay. In Sec. VI B, we extend these results to consider when $\tau >0$.

### A. No delay

With $\delta =\tau =0$, system (30) reduces to

A boundary equilibrium is no longer feasible. Hence, the interior equilibria, $E1$ and $E2$, are the only possible equilibria, where $E2$ corresponds to the intermediate tumor given when $T\u2217=\u2212A1\u2212A12\u22124A0A22A0$. Here, $A0,A1,$ and $A2$ are given by (8) with $\delta =\tau =0$.

The characteristic equation $\lambda 2+a1\lambda +b1=0$ with

admits the eigenvalues

Hence, the stability of the equilibria depends on $\Delta =a12\u22124b1$. The stability and existence of the equilibria are examined in the following five cases.

#### 1. Case 1: *b*_{1} = 0

We first must examine whether $b1=0$ is feasible for both real, positive equilibria. By (32b),

when $b1=0$. Hence, $b1=0$ iff $A12=4A0A2$, and $E1$ and $E2$ coalesce into a single equilibrium. Since $\lambda 1=0$ and $\lambda 2<0$ when $b1=0$, the equilibrium $E1=E2$ is stable.

#### 2. Case 2: *b*_{1} < 0

Note that $b1<0$ when

Thus, $b1<0$ is only feasible for $E2$. By (33), $\lambda 1>0$ and $\lambda 2<0$, and $E2$ is unstable. Therefore, positive $E1$ is not feasible, and $E2$ is unstable when $b1<0$.

#### 3. Case 3: $0<b1<a124$

From $b1=(1+\eta t)(2A0T\u2217+A1)T\u2217>0$, we have

a contradiction for $E2$. Thus, positive $E2$ is not feasible when $b1>0$. To examine the existence of $E1$, we see the inequality

holds if and only if

where

Now, it must be that

else $(\eta e\u2212(1+\mu ))2<0,$ a contradiction. Since $A1<0$ and $A0,A2>0$, then the inequality

holds, and the existence of $E1$ is satisfied when $0<b1<a124$. Then, $E1$ is asymptotically stable, as $\lambda 1,\lambda 2<0$. Thus, positive $E2$ is not feasible, and $E1$ is asymptotically stable when $0<b1<a124$.

#### 4. Case 4: $b1=a124$

As shown in the previous case, positive $E2$ is not feasible since $b1>0$. Now, $b1=a124$ when

with constants $c1,c2,c11$ given by (34). Then,

and $T\u2217>0$ only when $|c2|>|c1|$ since $A0>0$ and $A1<0$. By definition, $c1>0$. Suppose $c2>0$. Then, $|c2|>|c1|$ yields $\u221234+c11>14+c11,$ a contradiction. Now suppose $c2<0$. Then $|c2|>|c1|$ gives $14>c11,$ a contradiction, as shown in Case 3. Thus, neither positive $E1$ nor positive $E2$ is feasible when $b1=a124$.

#### 5. Case 5: $b1>a124$

By Case 3, positive $E2$ is not feasible since $b1>0$. Now, $b1>a124$ when

with constants $c1,c2,c11$ given by (34). As shown in Case 3, we must have $c11\u226514$. Thus, we have

a contradiction. Thus, neither $E1$ nor $E2$ is feasible when $b1>a124$.

We can summarize the preceding results in the following theorem.

For the system (31), the following statements are true, with $a1,b1$ given by (32).

$E1$ is only feasible when $0\u2264b1<a124$. $E2$ is only feasible when $b1\u22640$.

If $b1=0,$ then $E1$ and $E2$ coalesce into a single stable equilibrium.

If $b1<0,$ then $E2$ is unstable.

If $0<b1<a124,$ then $E1$ is asymptotically stable.

### B. With delay

The previous results lend themselves to an extension to determine the impact on the stability when $\tau >0$. We first linearize (30) at the interior equilibrium $(E\u2217,T\u2217)$. Setting $x=E\u2212E\u2217,y=T\u2212T\u2217$, where $x,y$ are small, gives

Thus, a non-trivial solution of the form $(x,y)=(c1,c2)e\lambda \tau $ exists if and only if $G(\lambda ,\tau )=0$, with

The characteristic equation (38) can be written in a simpler form as

with

Recall, purely imaginary roots $\lambda =i\omega $ of the characteristic equation indicate possibilities of stability switches. Suppose $\lambda =i\omega $ is a root of the characteristic equation (39). This assumption yields

where the real and imaginary parts are, respectively, given by

Thus,

Its roots are

By (40),

Since $a2>2c$, $\omega \u22122$ will always be negative. Additionally, it always holds that $a,d>0$, but sgn($c$) is unknown and dependent on the parameter choice. We examine several cases via application of Theorem 3.1 in Kuang (Ref. 29, p. 77).

#### 1. Case 1: |*c*| > *d*

Suppose $|c|>d$. Then, the right-hand side of (42) is always negative. By Theorem 3.1 of Kuang (Ref. 29, p. 77), the stability of $(E\u2217,T\u2217)$ does not change for any $\tau \u22650$. As the stability is unchanging, we determine the stability of $(E\u2217,T\u2217)$ when $\tau =0$.

If $c<0$, then $c+d<0$. By Theorem 4, $E2$ is unstable and the only positive equilibrium for the system (30) when $\tau =0$. If $c>0$, then $c+d>0$. It follows from the application of Theorem 4 that $E2$ is not feasible, and $E1$ is asymptotically stable for all $\tau \u22650$, as depicted by Fig. 14.

#### 2. Case 2: |*c*| < *d*

Suppose $|c|<d$. Evaluating (42) yields

Hence, there exists one imaginary root $\lambda =i\omega $ with a positive imaginary part for the characteristic equation ((39)). Since $\omega \u2212$ is always negative, we never have two imaginary roots with a positive imaginary part. When $\tau =0$, by Theorem 4, $E1$ is asymptotically stable and $E2$ is not feasible since $c+d>0$.

Define

where $0\u2264\theta 1\u22642\pi $. Theorem 3.1 of Kuang (Ref. 29, p. 77) gives that $E1$ is uniformly asymptotically stable for $\tau <\tau 0,1$ and unstable for $\tau >\tau 0,1$. By (41),

Thus, the bifurcation point occurs when

Figure 15 illustrates the loss of stability as $\tau $ increases.

#### 3. Case 3: |*c*| = *d*

Suppose $|c|=d$. Evaluating (42) gives

Thus, $\omega \u22122<0$ and $\omega +2=0$. Suppose $c>0$. Then, $c+d>0$, and Theorem 4 gives that $E1$ is the only positive equilibrium and is asymptotically stable when $\tau =0$. Now, if the stability of $E1$ changes at some $\tau >0$, the characteristic equation (39) admits a root $\lambda =u+iv$ with $u>0$ for some $\tau >0$. Suppose (39) admits a complex root $\lambda $ with a positive real part. Substituting $\lambda =u+iv$ in the characteristic equation (39) yields

The real and imaginary parts are, respectively, given by

Squaring and adding (44) yields

Since $u>0$ and $\tau \u22650$ by assumption, it follows that $d2e\u22122u\tau \u2264d2$. Therefore,

With $c,u>0,a2>2c,$ and $|c|=d$, expanding (46) gives

a contradiction. Thus, it follows that all roots of (38) have nonpositive real parts. As $\lambda (\tau )=0$ is clearly never a root of (39) when $c>0$ and $|c|=d$, then $E1$ is asymptotically stable for all $\tau \u22650$, as shown in Fig. 16.

We now suppose $c<0$. Then, $|c|=d$ implies that $c+d=0$. By Theorem 4, $E1=E2$, and the equilibrium is stable. To determine if the stability remains the same for all $\tau \u22650$, follow the assumption in Kuang (Ref. 29, pp. 79, 80) and let $\tau >ad>0$. Consider the following function:

Clearly, $F(0,\tau )=0$ and $lim\lambda \u2192\u221eF(\lambda ,\tau )=\u221e$. Additionally, there exists $N>0$ such that if $\lambda \u2265N,F(\lambda ,\tau )\u22650$. We also have

Then, $\u2202F(0,\tau )\u2202\lambda =a\u2212d\tau <0$ since $\tau >ad$ by assumption. Then, when $\tau >ad$, there exists $\delta (\tau )>0$ such that when $0<\lambda \u2264\delta (\tau ),F(\lambda ,\tau )<0$. Thus, there must exist $\lambda \u2217$ with $\delta (\tau )<\lambda \u2217\u2264N$ such that $F(\lambda \u2217,\tau )=0$, making (38) have a positive root. Hence, when $c<0$ and $|c|=d$, the equilibrium $E1=E2$ is stable for $\tau <ad$ and unstable for $\tau >ad$, as displayed in Fig. 17.

The preceding results are summarized in the following theorem.

For the system (31), the following statements are true, with $a,b,c,d$ given by (40).

If $|c|>d,$ $E2$ is unstable and the only feasible equilibrium for $\tau \u22650$ when $c<0$. For $c>0$, $E1$ is asymptotically stable and the only feasible equilibrium for $\tau \u22650$.

If $|c|<d,$ $E1$ is the only feasible equilibrium and is uniformly asymptotically stable for $\tau <\tau 0,1$ and unstable for $\tau >\tau 0,1,$ where $\tau 0,1$ is calculated by (43).

If $|c|=d,$ $E1$ is asymptotically stable and the only feasible equilibrium for $\tau \u22650$ when $c>0$. For $c<0,$ $E1=E2,$ and the equilibrium is stable when $\tau <ad$ and unstable when $\tau >ad$.

When the tumor is highly aggressive, a less rapid immune response (represented by delay) can lead to a loss of stability for a tumorous equilibrium, indicating a loss of tumor control by the immune system. The instability initially reflects a period of tumor growth, followed by uncertainty in tumorous outcomes.

## VII. HEALTH-CARE IMPLICATIONS

The reduction to a system of two equations yields a model simple enough to allow for clinical use, yet complex enough to produce rich dynamics. Three values, the activation and proliferation rate of CTLs, the time for activated CTLs to kill tumor cells, and the inactivation rate of CTLs by the tumor, are key in driving the complexity, and the importance of taking their measurements should, therefore, be emphasized to health-care professionals. In evaluating these parameters, as many of the other parameters are commonly known from past experiments, health-care workers would then be able to determine which parameter space the tumor falls in and what would accordingly be the best course of treatment. Additionally, the model results suggest to clinicians which combination treatments would be best-suited to pair with DC vaccines. For any patient, treatment to increase the speed at which activated CTLs leave the spleen to kill tumor cells would be a wise combination, as a quick response time results in tumor control, where the tumor is maintained at a low level. When CTLs are activated at a lower rate, the immune response has more flexibility in response time as opposed to when CTLs are activated at a rapid rate, where only a small window of response time results in tumor control. Additionally, treatment to lessen the immunosuppressive environment of the tumor allows for the use of less DC therapy to eliminate the tumor, with the model providing a threshold to guarantee tumor elimination.

## VIII. DISCUSSION

The underwhelming results of DC-based clinical trials have led to greater testing of DC vaccines in combination treatments to enhance efficacy. In order to understand when the vaccine will perform most effectively, it is necessary to understand the interactions of tumor and immune cells under different conditions. Mathematical models allow for exploration of these interactions and can aid clinicians in designing better suited monotherapy and combination treatments.

We propose a simple mathematical model that is capable of exhibiting complex dynamics. We incorporate a constant time delay to represent the time for the immune system to respond to the tumor. We have analytically and numerically proven the existence of a Hopf bifurcation and provided a threshold to ensure tumor existence ($R0>1$). Our reduced model is capable of exhibiting bistability in the region $Rcrit<R0<1$. Numerical experiments suggest that the threshold for tumor elimination is $R0<Rcrit<1$. In the special case when the immunosuppressive tumor microenvironment is neutralized ($\eta e=0$), less DC treatment is necessary, with the model (5) guaranteeing tumor elimination for $R0<1$. Combination treatments pairing DC vaccines with agents to block or neutralize immunosuppressive factors, such as anti-IL-10 and anti-TGF-$\beta $, would improve the efficacy of the DC vaccine response, lessening the amount of treatment needed for tumor elimination.

The model outcomes are shown to be sensitive to the time delay. While larger delays are commonly known to destabilize a system through Hopf bifurcation, our model exhibits richer dynamics than what is commonly observed. Instead of remaining unstable following Hopf bifurcation, our model regains its stability for a sufficiently large delay. Additionally, our model may produce regions of increased complexity, where the equilibrium becomes increasingly unstable and can even lead to chaotic outcomes (three pairs of complex conjugate eigenvalues with a positive real part).

Analytical work and numerical experiments additionally reveal the activation/proliferation rate of CTLs (captured by $\rho $) to be critical in the dynamics of the system (5). When the activation/proliferation of CTLs is low ($\rho $ small), the immune system can control the tumor cells for small $\tau $, allowing for the coexistence of both populations at a low level. As the delay increases, a Hopf bifurcation introduces oscillatory dynamics for a window of $\tau $. Further increase of $\tau $ leads to a return to a stable positive steady state, where a high level of tumor coexists with a low level of CTLs. Taking $\rho $ to be larger, the oscillatory dynamics persist for a wider region of delay and gain added complexity before regaining stability. These complex, aperiodic regions reflect the clinically relevant variability in treatment outcomes dependent on initial tumor burden. Furthermore, the oscillatory behavior represents the phenomenon of tumor dormancy for long periods of time, suddenly followed by tumor reappearance and growth to a lethal size for unknown reasons, which has been demonstrated *in vivo*^{33} as well as in a variety of mathematical models.^{12,34} In all cases, when the time for the immune system to respond is too long, the effector cells are unable to control the tumor cells at a low level, and the tumor cells increase to a fatally high level for all initial conditions, aligning with biological intuition.

An analysis of the system when effector cell death is dominated by the interaction with the tumor cells ($\delta =0$) reveals that an increased delay is strictly shown to have a destabilizing effect on the high tumor burden steady state when conditions allow for a change in the stability. The system with $\delta >0$ exhibits large delays stabilizing and intermediate delays destabilizing the high tumor burden steady state. Biologically this makes sense, as when the immune system takes too long to respond, then there is less resistance to tumor growth, and the tumor will persist. The results highlight the importance of incorporating delay in mathematical models to capture these dynamics and clinically measuring the delay to evaluate the best course of treatment.

While we represented a constant delay in the response of the immune system, other modes of incorporating delay would be worth exploring in future work. With innate delays present in additional cellular dynamics, such as the binding time required for activation, we could examine whether incorporating a second discrete delay would be significant in the system behavior. Though less frequent, mathematical models of tumors with two discrete delays have been studied.^{35–37} While the Beretta and Kuang^{30} method employed in this article has allowed for the evaluation of one-delay systems with delay dependent coefficients, an efficient method for evaluating a two-delay system with delay dependent coefficients remains an open problem.^{36}

## ACKNOWLEDGMENTS

Research of Y.K. was partially supported by the National Science Foundation (NSF) grants (Nos. DMS-1615879 and DEB-1930728) and an National Institutes of Health (NIH) grant (No. 5R01GM131405-02).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## REFERENCES

*et al*., “

*BIOMAT 2012*, edited by R. P. Mondaini (World Scientific, 2013), pp. 70–85.

*Delay Differential Equations with Applications in Population Dynamics*, Mathematics in Science and Engineering Vol. 191 (Academic Press, 1993).