Analysis of tumor-immune dynamics in a delayed dendritic cell therapy model

Since their landmark discovery in 1973 by Steinman and Cohn,¹ 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.² 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.³ 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.⁴ 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.⁵ 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.,¹⁸ 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.¹⁹ The simplified system in Dickman et al.¹⁸ 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 $μ$ 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.

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.¹⁸ 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.

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 $α$⁠. Additionally, we assume that the DC therapy takes $τ$ units of time to become effective; i.e., the effector cells take $τ$ 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−τ$⁠. Since activated effector cells die at a rate of $δ~E$⁠, the probability of survival for the effector cells is given by $e−δ~Eτ$ during this period of delay. The assumptions give rise to the following model:

where $c~=cα$⁠, and we assume that the initial values are

Parameter values of system (2) are specified in Table I.

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.²⁰ with a single set of biologically reasonable parameters given in Table I. Data from Lee et al.²⁰ reflect the average tumor volume of melanoma-induced mice treated with $0,1×105,7×105,$ or $21×105$ DCs on days 6, 8, and 10. The fit to the data remains comparable to that of Dickman et al.¹⁸ 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,²¹ a singular Hopf occurs at a distance that is $O(ε)$ 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.

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.²² 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.²³ In particular, the anti-apoptotic BCL-xL is elevated in activated effector cells,²⁴ 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δD,D~i=DiδD$⁠. The simplifications can again be justified through the fit to the Lee et al.²⁰ data, as shown in Fig. 3, where the fitting is comparable to all previous model formulations.

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

where $β=sErk,γ=c~v~rk,δ=δ~Er,ρ=c~D~ir,ηe=cekr,ν=(1−α)v~k,μ=(1−α)D~i,ηt=ctkr,andτ¯=rτ$ are the dimensionless, positive parameters (with the exception of $γ=ν=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.

Positivity and boundedness of the dimensionless system (11) can similarly be proven as before. The tumor-free equilibrium is given by $(E0,T0)=((β+γe−δτ)/δ,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.¹⁸ 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.¹⁸ model.

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

Figure 5 clearly illustrates how the interior equilibria both depend on $τ$⁠. For lower $τ$⁠, $E0$ exists, but neither interior equilibria are feasible. Increasing $τ$ 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 $τ$⁠, corresponding to $R0=1$⁠. The condition $Rcrit<R0<1$ required for two equilibria is thus shown to only be satisfied for finite $τ$⁠. Though $E2$ is not feasible for a higher value of $τ$⁠, the interior equilibrium $E1$ remains feasible as $τ→+∞$⁠.

Thus, for $ηe$ sufficiently large such that $ηe>C3$⁠, there always exists at least one positive equilibrium root when $R0≥Rcrit$⁠. Additionally, the second interior equilibrium only exists when the delay $τ$ 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 $τ$ or other parameters vary.

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 $γβ=γ+β$⁠. Figure 6 depicts the nullclines and equilibria when the conditions $Rcrit<R0<1$ and $ηe>C3$ are satisfied for the existence of two interior equilibria.

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

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 $ηe$ is derived from the inactivation rate of CTLs by the tumor (⁠$ce$⁠), and sufficient reduction of $ηe$ eliminates the region of bistability. Agents to inhibit immunosuppressive factors, such as anti-interleukin-10 (anti-IL-10) and anti-transforming growth factor-$β$ (anti-TGF-$β$⁠), would then be suitable candidates for combination therapies.²⁸ 

In this section, we study the linear stability of the interior equilibria, $E1$⁠, $E2$⁠, assuming that the conditions $η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 $τ$⁠, as $R0$ monotonically increases in $τ$⁠, 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∗,T∗)$⁠. Setting $x=E−E∗,y=T−T∗$⁠, where $x,y$ are small, gives

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

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

where $P(λ,τ)=λ2+A(τ)λ+C(τ)$ and $Q(λ,τ)=B(τ)λ+D(τ)$ with

Note that the characteristic equation (14) involves delay not just in the $e−λτ$ term, but in several other places as well, since $T∗$ and $E∗$ depend on $τ$⁠. 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²⁹ and result in clean, explicit calculations, such as the exact values of $τ$ where stability switches occur. We apply the approach from Beretta and Kuang³⁰ for studying these challenging characteristic equations, which analyzes the stability switches via computational means, as exact values of $τ$ where these switches occur cannot be analytically found.

In order to apply the technique of Beretta and Kuang,³⁰ we must verify that the following properties hold for all $τ∈R+$⁠:

• $P(0,τ)+Q(0,τ)≠0,∀τ∈R+0$⁠;

• If $λ=iω$⁠, $ω∈R$⁠, then $P(iω,τ)+Q(iω,τ)≠0,τ∈R$⁠;

• $lim sup{|Q|(λ,τ)/P(λ,τ)|:|λ|→∞,Reλ≥0}<1$ for any $τ$⁠;

• $F(ω,τ):=|P(iω,τ)|2−|Q(iω,τ)|2$ for each $τ$ has a finite number of zeros; and

• Each positive root $ω(τ)$ of $F(ω,τ)=0$ is continuous and differentiable in $τ$ whenever it exists.

We first have $P(0,τ)+Q(0,τ)=C+D$⁠. By assumption, property (i) is established, implying that $λ=0$ is not a root of (13). Now, $P(iω,τ)+Q(iω,τ)=−ω2+iω(A+B)+C+D≠0$⁠, satisfying property (ii). Clearly, property (iii) holds, as $lim sup|λ|→∞|Q(λ,τ)P(λ,τ)|=0$⁠. Additionally, we can clearly see that $F(ω,τ):=|P(iω,τ)|2−|Q(iω,τ)|2$ has a maximum of four roots since

thus satisfying property (iv). Let $(ω0,τ0)$ be an arbitrary point in its domain such that $F(ω0,τ0)=0$⁠. Then, we see that $Fω$ and $Fτ$ exist and are continuous in a neighborhood of $(ω0,τ0)$⁠. Additionally, $Fτ(ω0,τ0)≠0$⁠. 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.³⁰ 

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 $λ(τ)=±iω(τ),ω(τ)∈R+$⁠, which crosses the imaginary axis at some $τ=τ∗$⁠, a necessary condition for a Hopf bifurcation to occur.

To find the value of $τ∗$ where a stability switch may occur, without loss of generality, we let $λ(τ)=iω(τ)$⁠. Substitution of $λ=iω$ into (13) yields the following:

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

Then, $ω$ must be a solution of

Using (15) to rewrite (18), we reach

with its roots as

where $Δ=(2C−A2)2−4(C2−D2).$ We know that $ω−2$ is always negative, as

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

By Theorem 1, $C+D>0$ at $E1$ and $C+D<0$ at $E2$⁠. Since $D≥0$⁠, we can deduce $C<0$ at $E2$⁠, making $C−D<0$ at $E2$⁠. Thus, a real $ω$ does not exist, and there are no stability switches for $E2$⁠. It then follows that $E2$ is unstable for all $τ≥0$⁠. The sign of $C−D$ 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 $ω$ exists for $E1$ when $C−D<0$⁠, the characteristic equation (13) has a pair of simple and conjugate pure imaginary roots $λ=±iω$⁠. By the application of Theorem 4.1 of Beretta and Kuang,³⁰ we prove that this complex conjugate pair of eigenvalues crosses the imaginary axis, ensuring a stability switch for $E1$ when $C−D<0$⁠.

Assume that $I⊆R+0$ is the set where $ω(τ)>0$ is a root of (19). For any $τ∈I$⁠, we use (17) to define the angle $θ(τ)∈[0,2π]$ as follows:

For $τ∈I$⁠, we have the relationship

From (23), we see

If the conditions hold for an imaginary solution to exist, making $ω$ feasible for $τ∈I$⁠, then we have the continuous and differentiable sequence of functions $I→R$⁠,

The characteristic equation (13) has a pair of simple conjugate pure imaginary roots $λ=±iω(τ∗)$ at $τ∗∈I$ when $Sn(τ∗)=τ∗−τn(τ∗)=0$ for some $n∈N0$⁠. By application of Theorem 4.1 in Beretta and Kuang,³⁰ 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

Substituting (28) into (27) yields

Evaluating at $λ=iω(τ∗)$⁠, 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 $τ$ increases, recalling the relationship

by Theorem 2.

We consider two scenarios to demonstrate the impacts of the delay $τ$ and $ρ$ on $E1$⁠. Recall that $ρ$ 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(τ),S1(τ),$ and $S2(τ)$ as defined above. For both scenarios, $C−D<0$ up to a finite $τ$ such that real $ω$ exists and Theorem 2 can be applied.

Figure 7 illustrates the case of the low activation/proliferation rate of the CTLs (⁠$ρ$ small). The function $S0(τ)$ is shown to have roots at $τ01$ and $τ02$⁠, with $τ01<τ02$⁠. For the choice of parameters, no real roots exist for $Sn(τ)$ when $n>0$⁠. The behavior of $S0$ indicates that $E1$ is asymptotically stable, loses its stability for $τ∈(τ01,τ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 $τ$ axis and the characteristic equation about $E1$ admits two roots with a positive real part.