A model of dendritic cell therapy for melanoma

Dendritic cells are a promising immunotherapy tool for boosting an individual’s antigen-specific immune response to cancer. We develop a mathematical model using differential and delay-differential equations to describe the interactions between dendritic cells, effector-immune cells, and tumor cells. We account for the trafficking of immune cells between lymph, blood, and tumor compartments. Our model reflects experimental results both for dendritic cell trafficking and for immune suppression of tumor growth in mice. In addition, in silico experiments suggest more effective immunotherapy treatment protocols can be achieved by modifying dose location and schedule. A sensitivity analysis of the model reveals which patient-specific parameters have the greatest impact on treatment efficacy.

1. Introduction

A promising immunotherapy approach to treating certain cancers involves the use of dendritic cells (DCs). DCs are part of the antigen-specific (adaptive) immune response and function as antigen-presenting cells. Immature DCs are derived in the bone marrow and reside in peripheral tissues. Upon encountering pathogen, DCs begin to mature, and travel to the lymphoid organs where they stimulate differentiation and maturation of cytotoxic T lymphocytes (CTLs). Some of these activated CTLs then travel to the infected tissue to form part of the adaptive immune response, while others become memory cells that are ready to mount a rapid response in case of a rechallenge by the pathogen.

Previous studies have established the efficacy of dendritic cell treatments for tumors in the murine system (DeMatos et al., 1998; Fields et al., 1998; Lee et al., 2007; Yamaguchi et al., 2007; Shinagawa et al., 2008). In these studies, DCs have been shown both to inhibit the growth of nascent tumors and to provide a memory response to previously encountered antigen. In the clinic, researchers have been able to extract immature dendritic cells from patients, culture them ex vivo, and load them with tumor antigens to create an individual-based vaccine that can boost a patient’s response against their own cancerous cells (Pilon-Thomas et al., 2004; Taquet et al., 2008). The success of clinical trials of DC vaccines has resulted in the recent FDA approval of the first cancer vaccine for prostate cancer (Cheever, 2011). Despite promising clinical responses in vaccine trials, it remains difficult to predict which patients will actually respond to these vaccines and why (Trefzer et al., 2005; Boon et al., 2006). Mathematical models of DC therapy can provide insight into the mechanisms driving the kinetics of the immune response that may lead to these disparate patient responses.

Cell trafficking is an important aspect of the DC-mediated immune response. DCs must travel from the tumor to the peripheral lymph organs via the blood, and activated CTLs must travel from the lymph organs back to the tumor. Ludewig et al. (2004) have developed a model describing DC and CTL trafficking in mice. The model includes activated and memory CTLs to capture both the immediate and long-term effect of DC injections. The DC trafficking model of Ludewig et al. was carefully calibrated using experimental data from murine studies.

In this paper we present an extension and modification of the model in Ludewig et al. (2004). Our extended model includes a tumor compartment to allow for analysis of various DC treatments and their effect on tumor growth, as well as the long-term behavior of the system. We find relevant model parameters using the data collected by Lee et al. (2007) describing tumor growth in response to varying levels of DC injections. We compare model simulations of various DC doses, injection sites, and dose times. We include a comparison to a prophylactic dosing schedule presented by Preynat-Seauve et al. (2007).

2. The Model

The compartment model proposed by Ludewig et al. (2004) includes dendritic cells, activated CTLs, and memory CTLs. Our extended model includes tumor cells in addition to these immune cell populations. Adding a tumor compartment requires the determination of tumor-immune system parameters such as immune cell trafficking rates to and from the tumor, effector cell deactivation rates by tumor cells, effector cell death rates, intrinsic tumor growth rates, and tumor cell kill rates by effector cells. We note that this compartment model does not account for the geometry of the system. In particular, it does not explicitly incorporate the distance between the spleen and the tumor. However, in murine models, the transit times between compartments are small relative to the tumor growth time scale, so this simplification is reasonable. In this section we describe the processes included in the mathematical model.

2.1. Model Definition

Our model consists of three compartments: the spleen, the blood, and the tumor. Dendritic cells and active effector cells can move between the blood and spleen compartments, and between the blood and tumor compartments. We assume that memory effector cells can move between the spleen and the blood compartments. The system is not conservative: all types of cells are cleared through the blood, immune cells are created in response to the presence of tumor, and tumor cells grow according to a logistic growth law.

The nine state variables in our model are:

• D_blood, the number of dendritic cells in the blood compartment;

• D_spleen, the number of dendritic cells in the spleen compartment;

• $E_{blood}^a,$ the number of activated CTLs in the blood compartment;

• $E_{spleen}^a,$ the number of activated CTLs in the spleen compartment;

• $E_{blood}^m,$ the number of memory CTLs in the blood compartment;

• $E_{spleen}^m,$ the number of memory CTLs in the spleen compartment;

• $E_{tumor}^a,$ the number of activated CTLs in the tumor compartment;

• T, the number of tumor cells;

• D_tumor, the number of dendritic cells in the tumor compartment, the tumor-infiltrating DCs.

We present the system of nine differential equations in groupings representing the blood, spleen, and tumor compartments. The model parameters are described in detail in Table A1 in Appendix.

2.1.1. Blood compartment

The equations describing DC and CTL flow in the blood are given by:

$\begin{array}{lll}\frac{d}{dt}{D}_{blood}& =-{\mu }_B{D}_{blood}+{\mu }_{TB}{D}_{tumor}+v_{blood}\left(t\right)& \end{array}$ (1)

$\begin{array}{lll}\frac{d}{dt}{E}_{blood}^a& ={\mu }_{SB}\left(D_{spleen}\right)E_{spleen}^a-{\mu }_{BB}{E}_{blood}^a& \end{array}$ (2)

$\begin{array}{lll}\frac{d}{dt}{E}_{blood}^m& ={\mu }_{SB}\left(D_{spleen}\right)E_{spleen}^m-{\mu }_{BB}{E}_{blood}^m& \end{array}$ (3)

where, as holds throughout the model, the μ parameters represent flow rates between compartments. We include the “trapping” term from Ludewig et al. (2004) which describes the observed phenomenon of activated CTLs being held back in the spleen in the presence of DCs:

$\begin{array}{rcll}{\mu }_{SB}\left(D_{spleen}\right)& ={\mu }_{SB}^{*}+\frac{\Delta \mu }{1+\frac{D_{Spleen}}{{\theta }_{shut}}},& & \\ \Delta \mu & ={\mu }_{SB}^{Normal}-{\mu }_{SB}^{*}.& & \end{array}$

The function v_blood(t) allows us to model injections of DCs into the blood. For example, two doses of 7 × 10⁵ each given on Day 0 and Day 7 could be described by the function:

$v_{blood}\left(t\right)=\left\{\begin{array}{cc}\frac{7\times {10}^5}{1∕48},\hfill & t\in \left(0,1∕48\right)\cup \left(7,7+1∕48\right),\hfill \\ 0\hfill & \text{otherwise}\text{.}\hfill \end{array}\right.$ (4)

2.1.2. Spleen compartment

The differential equations for the spleen compartment describe interaction, death, and recruitment of DCs and CTL. The equations include one delay which represents the synaptic connection time: the contact time required between DCs and effector cells in the spleen before proliferation can begin. The other interactions we account for in the system do not involve a required contact time and thus are modeled without delay. The dynamics of the populations in the spleen are described by:

$\begin{array}{llll}\frac{d}{dt}{D}_{spleen}& =MaxD\left(1-e^{\left(\frac{-{\mu }_{BS}{D}_{blood}}{MaxD}\right)}\right)& & \\ & -a_D{D}_{spleen}-b_{DE}{E}_{spleen}^a{D}_{spleen}& \end{array}$ (5)

$\begin{array}{llll}\frac{d}{dt}{E}_{spleen}^a& ={\mu }_{BSE}{E}_{blood}^a-{\mu }_{SB}\left(D_{spleen}\right)E_{spleen}^a+b_a{D}_{spleen}{E}_{spleen}^m& & \\ & +a_{E_{a}S}\left(D{C}_{on}{E}_{naive}-E_{spleen}^a\right)-r_{am}{E}_{spleen}^a& \\ & +b_p\frac{D_{spleen}\left(t-{\tau }_D\right)E_{spleen}^a\left(t-{\tau }_D\right)}{{\theta }_D+D_{spleen}\left(t-{\tau }_D\right)}& & \end{array}$ (6)

$\begin{array}{llll}\frac{d}{dt}{E}_{spleen}^m& =r_{am}{E}_{spleen}^a-\left(a_{E_m}+b_a{D}_{spleen}+{\mu }_{SB}\left(D_{spleen}\right)\right)E_{spleen}^m& & \\ & +{\mu }_{BSE}{E}_{blood}^m.& \end{array}$ (7)

Note that the term in equation (6),

$D{C}_{on}=\left\{\begin{array}{cc}0\hfill & \text{if }{D}_{spleen}\left(t\right)=0\hfill \\ 1\hfill & \text{if }{D}_{spleen}\left(t\right)>0.\hfill \end{array}\right.$

indicates that we do not allow for new CTLs in the absence of DCs. Thus, the populations we model only exist due to the presence of tumor and mature DCs.

The first term in equation (5) reflects our assumption that there is a maximum rate at which mature DCs can enter the spleen. This is in agreement with observations that DCs cannot enter the spleen at an unlimited rate. Based on a range of values for the maximum rate we have set MaxD to 400 (cells per hour), reflecting the parameter fit obtained with the data from Lee et al. (2007) and Preynat-Seauve et al. (2007). As noted above, in equation (6) the final term introduces a delay, τ, into the system that reflects the synaptic connection time. Mathematically, this delay introduces more complexity into the system, especially regarding the stability analysis of the equilibria (see Stability Analysis below).

2.1.3. Tumor compartment

The tumor compartment contains activated effector CTLs, DCs, and tumor cells. The interactions of these populations within the tumor are described by:

$\begin{array}{lll}\frac{d}{dt}{E}_{tumor}^a& ={\mu }_{BTE}\left(T\right)E_{blood}^a-a_{E_{a}T}{E}_{tumor}^a-c{E}_{tumor}^{a}T,& \end{array}$ (8)

$\begin{array}{lll}\frac{d}{dt}T& =rT\left(1-\frac{T}{k}\right)-\mathcal{D}T.& \end{array}$ (9)

$\begin{array}{lll}\frac{d}{dt}{D}_{tumor}& =\frac{mT}{q+T}-\left({\mu }_{TB}+a_D\right)D_{tumor}+v_{tumor}\left(t\right)& \end{array}$ (10)

where

${\mu }_{BTE}\left(T\right)={\mu }_{BB}\left(T∕\left(\alpha +T\right)\right),$

and

$\mathcal{D}=d\frac{{\left(\frac{E_{tumor}^a}{T}\right)}^l}{s+{\left(\frac{E_{tumor}^a}{T}\right)}^l.}$ (11)

The function v_tumor(t) is similar to v_blood(t) in the blood compartment, allowing us to inject DCs intratumorally in order to compare treatment protocols.

Note that in equation (9), tumor growth is fit to a logistic function as in previous models (de Pillis and Radunskaya, 2003; de Pillis et al., 2005, 2007, 2009; Cappuccio et al., 2006). The behavior of this particular model is robust to the choice of growth function, for example a Gompertz growth law gives similar results. However, we choose the logistic law since it provides a good fit to the experimental data we are using for model calibration (See Figure 1). Cytolysis of tumor cells by activated CTLs [equations (9) and (11)] is a ratio-dependent kill term introduced in de Pillis and Radunskaya (2003). Experimental results from Diefenbach et al. (2001) support ratio-dependent, antigen-specific killing, and the term has been employed to success in previous models (de Pillis and Radunskaya, 2003; de Pillis et al., 2005, 2007, 2009). The importance of tumor-infiltrating dendritic cells has been demonstrated in several studies. See, for example, Preynat-Seauve et al. (2007). We allow DCs in the tumor to increase as a saturation-limited function of the size of the tumor population.

FIGURE 1

Figure 1. Fits to data from Lee et al. (2007) and corresponding residuals. The vertical bars in the graph on the left are the error bars from the experimental data, and the solid line is the outcome of the model simulation using the parameters that minimize the distance to the median of those data. In the graph on the right, the thin vertical bars are the same error bars from the data, and the solid rectangles show the “residuals,” that is, the distances between the simulated outcomes and the data. The estimated parameter values are given in Table A1 in Appendix.

With this model we simulate a variety of treatment scenarios, including those investigated in Lee et al. (2007) and Preynat-Seauve et al. (2007). The model offers insight into how best to harness the tumor controlling potential of DCs.

3. Results

3.1. Parameter Determination

All unknown parameters were fit to data from Lee et al. (2007) using a Nelder–Mead least-squares algorithm. The data in Lee et al. (2007) were collected from C57BL/6 female mice who were subcutaneously inoculated with 5 × 10⁵ B16F10 melanoma cells. We take this as day 0 for the purpose of fitting unknown parameters so that we may use 5 × 10⁵ tumor cells as an initial condition. Injections of 1 × 10⁵, 7 × 10⁵, or 21 × 10⁵ DCs were given at days 6, 8, and 10, following inoculation with tumor cells. Additionally, a control group was injected with PBS according to the same schedule (Lee et al., 2007). See Figure 1.

3.2. Equilibria and Stability Analysis

In order to determine the long-term behavior of the system, we find the equilibria and determine their stability.

3.2.1. Determination of the equilibria

The system has multiple equilibrium values, determined by setting equations (1–10) to zero. One solution to this system is the zero, or disease-free, equilibrium. To find the remaining non-zero equilibria, we first write all the state variables at equilibrium as functions of T, then search for the values of T that solve all equations simultaneously. We use asterisks to denote the value of the variables at equilibrium. Therefore, if there exists a non-zero value T* that satisfies equation (9), we can use equation (8) to obtain

$E_{tumor}^{a*}={\left\{\frac{s{\left(T^{*}\right)}^l\left(k-T^{*}\right)}{T^{*}-k\left(1-d∕r\right)}\right\}}^{1∕l}.$ (12)

The value for $D_{tumor}^{*}$ can also be found in terms of T* using equation (10):

$D_{tumor}^{*}=\frac{m{T}^{*}}{q+T^{*}}\frac{1}{{\mu }_{TB}+a_D},$ (13)

where we have replaced DC_death with its assumed constant value a_D.

Given $E_{tumor}^{a*}$, we can use equation (8) to determine the equilibrium value of the active effector cells in the blood:

$E_{blood}^{a*}=\frac{\left(a_{E_{a}T}+c{T}^{*}\right)\left(\alpha +T^{*}\right)}{{\mu }_{BB}{T}^{*}}{E}_{tumor}^{a*}.$ (14)

We use equations (13) and (1) to obtain $D_{blood}^{*}$ in terms of T*:

$D_{blood}^{*}=\frac{{\mu }_{TB}{D}_{tumor}^{*}}{{\mu }_B}=\frac{{\mu }_{TB}}{{\mu }_B}\left(\frac{1}{{\mu }_{TB}+a_D}\right)\frac{m{T}^{*}}{q+T^{*}}.$

Equation (3) gives an expression for $E_{blood}^{m*}$ in terms of $E_{spleen}^{m*}$,

$E_{blood}^{m*}=\frac{{\mu }_{SB}\left(D_{spleen}\right)}{{\mu }_{BB}}{E}_{spleen}^{m*}.$

Turning to the spleen compartment, we have:

$E_{spleen}^{a*}=\frac{{\mu }_{BB}}{{\mu }_{SB}\left(D_{spleen}\right)}{E}_{blood}^{a*}.$

Using (14) and (12), this gives $E_{spleen}^{a*}$ in terms of T*. According to equation (5), knowing $D_{blood}^{*}$ allows determination of $D_{spleen}^{*}$. Using equation (6) results in the following quadratic equation for $D_{spleen}^{*}$:

$\begin{array}{ll}0& =-{\theta }_{shut}\left({\mu }_{SB}^{*}+\Delta \mu \right){\mathcal{D}}_{in}\left\{{\theta }_{shut}{\mu }_{SB}^{*}{a}_D+\Delta \mu {\theta }_{shut}{a}_D\right.\\ & \left.+{\theta }_{shut}{b}_{DE}{\mu }_{BB}{E}_{Blood}^{a*}-{\mu }_{SB}^{*}{\mathcal{D}}_{in}\right\}{D}_{spleen}^{*}\\ & \left\{{\mu }_{SB}^{*}{a}_D+b_{DE}{E}_{Blood}^{a*}{\mu }_{BB}\right\}{\left(D_{spleen}^{*}\right)}^2,\end{array}$ (15)

where

${\mathcal{D}}_{in}=MaxD\left\{1-exp\left(\frac{-{\mu }_{BS}{D}_{blood}}{MaxD}\right)\right\}.$

Solving this quadratic equation yields two different, relevant equilibrium values for $D_{spleen}^{*}$. From equation (7) we get a value for $E_{spleen}^{m*}$ for each value of $D_{spleen}^{*}$:

$E_{spleen}^{m*}=\frac{r_{am}{E}_{spleen}^{a*}}{a_{E_m}+b_a{D}_{spleen}^{*}+{\mu }_{SB}\left(D_{spleen}^{*}\right)\left\{1-{\mu }_{BSE}∕{\mu }_{BB}\right\}}.$

Finally, from equation (6), the roots of the following function, expressible in terms of one variable, T*, yield equilibrium values for T*.

$\begin{array}{l}Z(T^{*})={\mu }_{BSE}{E}_{blood}^{a*}-{\mu }_{SB}(D_{spleen}^{a*})E_{spleen}^{a*}+b_a{D}_{spleen}^{*}{E}_{spleen}^{m*}\\ +a_{E_{a}S}(D{C}_{on}{E}_{naive}-E_{spleen}^{a*})-r_{am}{E}_{spleen}^{a*}\\ +b_p\frac{D_{spleen}^{*}{E}_{spleen}^{a*}}{{\theta }_D+D_{spleen}^{*}}.\end{array}$ (16)

From equation (12), we see that as long as the values of T*, and thus the roots of (16), lie between k(1 − d/r) and k, a non-zero equilibrium state exists. Recall that d, r, and k are the parameters that represent the tumor cell kill, intrinsic growth rates, and tumor carrying capacity, respectively. The function Z(T) is plotted in Figure 2 for the parameter set given in Table A1 in Appendix.

FIGURE 2

Figure 2. Graph of the function Z given in equation (16). The two zeros, marked by dots, correspond to equilibrium values of T. The function becomes complex as T → 10⁹. The two curves correspond to the two roots of equation (15).

3.2.2. Stability of the equilibrium points

A stability analysis of the system of delay equations (1–10) can be carried out by analyzing the linear approximation to the system at an equilibrium point. Since the term $\mathcal{D}\left(T,E_{tumor}^a\right),$ given in equation (11), is not differentiable at (0, 0), the system of DEs is not differentiable and, hence, has no linear approximation at the origin. Although we cannot use the linearization in this case, we do have numerical simulations that indicate that the tumor free equilibrium is initially unstable, but gains stability as the value of d, the immune strength parameter, is increased (see Figure 9). There is ongoing investigation of the analytical nature of the stability of the disease-free equilibrium.

At other equilibria, the linearization is given by two matrices of partial derivatives, J₀ and J_τ. To simplify the notation, we denote the nine state variables as x₁ through x₉ and the delayed state variables as z_i(t) = x_i(t − τ). If the rate of change of x_i is denoted by:

$\frac{d{x}_i}{dt}=F_i\left(x_1,\dots x_9,z_1,\dots z_9\right)$

then the entries of the Jacobians are:

$J_0\left(i,j\right)=\frac{\partial F_i}{\partial x_j},J_{\tau }\left(i,j\right)=\frac{\partial F_i}{\partial z_j}$

The formulas for the entries of J₀ and J_τ are given in Appendix B.

The eigenvalues of the derivative matrices at an equilibrium, E can be determined by finding the roots of the characteristic polynomial P(λ, τ), where P is defined by:

$P\left(\lambda ,\tau \right)=det\left(J_0\left(\mathbf{\text{E}}\right)+e^{-\lambda \tau }{J}_{\tau }\left(\mathbf{\text{E}}\right)-\lambda I\right)$

For most systems of delay-differential equations, determining the roots of the characteristic polynomial is a non-trivial process. Given our equations and particular parameter set, there exist two positive equilibrium values of T* (see Figure 2). However, only one of these values result in positive (biologically relevant) equilibrium values for all other state variables (see Appendix B). For this biologically relevant equilibrium, in the non-delay (or τ = 0) case, all eigenvalues of the derivative matrix J₀ have negative real part (see Appendix B). Because this is equivalent to the non-delay case, we know that when τ = 0, the non-trivial equilibrium is stable. It is possible for an equilibrium to change stability as the delay increases from zero. In our case numerical simulations suggest that this equilibrium maintains its stability even for large values of the delay, τ (see Figure 3).

FIGURE 3

Figure 3. This figure presents a sampling of several simulations for various values of the delay, τ and several different initial conditions (IC). The solid black line (topmost curve) shows a simulation with initial value 9.82 × 10⁸ with τ = 16. The other curves (from top to bottom) are the results of simulations using (τ, IC) pairs: (2, 9.817 × 10⁸), (0.5, 9.835 × 10⁸), (0.1, 9.811 × 10⁸), (8, 9.8 × 10⁸). For each value of τ, we simulated several initial conditions ranging between 9.80 × 10⁸ and 9.82 × 10⁸ (other simulations not shown) and, in each case, the cell populations approached the same equilibrium value.

3.3. Calibration and Validation of DC Effect on Tumor Growth

In this section we discuss the validation of the model. Starting with parameter values estimated in Ludewig et al. (2004), we then calibrate DC and CTL dynamics against the data provided by the experiments in Lee et al. (2007). In the discussion section we present a number of numerical experiments in which we explore the difference between intravenous and intratumoral DC injections, as well as modifications of the dose timings with hypothetically improved treatment schedules.

Experiments carried out in Lee et al. (2007) give tumor growth data, both in the presence and absence of DC treatment. We fit our intrinsic tumor growth parameters to the PBS melanoma growth data provided by Lee et al. (2007). Since specific trafficking parameters to and from the tumor have not been measured, we used the tumor growth data provided to infer the parameter values needed for the tumor compartment DC and CTL dynamics.

In Lee et al. (2007), groups of three 6–8-week-old female C57BL/6 mice were challenged with 5 × 10⁵ B16F10 melanoma cells on day 1, then treated with DC injections starting on days 6, 8, and 10. In separate experiments, DC doses of size 1 × 10⁵, 7 × 10⁵, and 21 × 10⁵ were administered. Lee et al. point out that the largest DC dose is most effective at slowing tumor growth. In fact, according to Lee et al., the largest dose regimen of 21 × 10⁵ DCs injected three times provided up to 41% tumor growth suppression as compared to the control mice. Survival time for these mice was increased by approximately 60%. We note that, even with the most aggressive DC treatment attempted, tumor growth was not completely suppressed.

In Figure 4, left panel, we see the change in tumor volume over 20 days, and compare tumor growth with no DC treatment to growth with varying levels of DC treatment. Simulated DC doses are 1 × 10⁵, 7 × 10⁵, and 21 × 10⁵, reflecting the laboratory experiments of Lee et al. (2007). After a tumor challenge of 5 × 10⁵ B16F10 melanoma cells on day 1, DC injections of the specified doses were then given intratumorally on days 6, 8, and 10. The simulation results fall well within the data ranges provided by Lee et al. (2007).

FIGURE 4

Figure 4. Fractionated dosing comparison. Intratumoral injections. Compare original DC dosing schedule (left panel) to hypothetical fractionated dosing schedule (right panel).

4. Discussion

In Section 3 we validated and calibrated the model, and analyzed the long-term behavior of the system. We are now in a position to explore hypothetical treatment variations. In this section, we discuss the effects of varying treatment protocols, and possible implications for patients. In Section 4.1, we compare intratumoral DC injections to hypothetical intravenous DC injections. We will see that when injecting the smallest dose of DCs, hypothetical intravenous DC injections are more effective at suppressing tumor growth than are DCs injected directly into the tumor. However, intratumoral injections are more effective than intravenous injections when the highest DC dose is used. In Section 4.2, we explore the effect of modifying dose timings. We will see that fractionated doses that are administered intravenously delay tumor growth significantly. We will see that earlier treatment initiation also helps suppress tumor growth, but more so with intratumoral injections. Up to this point, we have found ways to slow tumor growth by varying dose timing and location, but have not been able to completely eliminate a tumor. In Section 4.3, we explore the effects of prophylactic DC dosing. We find that prophylactic DC dosing actually allows us to eliminate a tumor under the right circumstances. We will see that as long as the CTL immune response is sufficiently strong, as reflected by the immune strength parameter d, a tumor that is introduced after DCs are injected can be completely suppressed.

4.1. Intratumoral versus Intravenous Treatment

We first compare the effect of treatment at two different injection sites. In the work of Preynat-Seauve et al. (2007) DC trafficking resulting from different injection sites was compared. They observed that there is a “trapping effect” within the tumor: DCs injected intratumorally do not reach the lymph nodes in significant numbers, indicating that the DCs are “trapped” for a time within the tumor. Note that this trapping is a different phenomenon from the one described in Section 2.1 which referred to activated CTLs being held back in the spleen in the presence of DCs. Preynat-Seauve et al. observed that subcutaneous DC injections resulted in DCs getting to the lymph nodes in greater numbers. In our numerical experiments, we compare intratumoral DC injections (as was done in the Lee et al. (2007) experiments) to intravenous DC injections. We note that intravenous injections and subcutaneous injections are not equivalent, but both approaches do avoid the trapping effect of intratumoral injections.

In Figure 5, left panel, we see the effects on melanoma growth of hypothetical intravenous DC injections over a 20 day period. Simulated DC doses of size 1 × 10⁵, 7 × 10⁵, and 21 × 10⁵ are administered. The no-treatment tumor growth case is also included for comparison. After a tumor challenge of 5 × 10⁵ B16F10 melanoma cells on day 1, DC injections of the specified doses were then given on days 6, 8, and 10.

FIGURE 5

Figure 5. Fractionated dosing comparison. Intravenous injections. Compare original DC dosing schedule (left panel) to hypothetical fractionated dosing schedule (right panel). 5 × 10⁵ B16F10 melanoma cells on day 1, the original dosing calls for DC injections of 7 × 10⁵, and 21 × 10⁵ administered every other day on days 6, 8, 10 (pictured in graphs 0.25 × 10⁵, 1.75 × 10⁵, and 5.25 × 10⁵ twice a day on days 6, 7, 8, 9, 10, 11.

It is interesting to note that in these simulations, all three intravenous dose responses appear to be nearly equally effective. There is little difference between low-dose and high-dose intravenous injection outcomes, while there is a significant difference between low-dose and high-dose intratumoral injection outcomes. In addition, when we compare Figure 5, left panel, to Figure 4, left panel, we see that all three intravenous doses control the tumor growth about as effectively as the mid-sized 7 × 10⁵ intratumoral dose (that is, the intravenous doses are all more effective than the lowest tumor dose, but less effective than the highest tumor dose). The reason for this result can be explained mathematically by the presence of the MaxD term in the model, equation (5). This term limits the rate at which DCs can enter the spleen, which in turn limits how saturated with DCs the spleen can get. Consequently, this limited inflow rate works against any treatment that attempts to send DCs into the spleen too quickly. MaxD term limits DC inflow into the spleen, the total number of DCs in the spleen over MaxD limit. In the case of the low-dose injections, the DCs injected intravenously all enter the spleen, since their entry rate is not limited by MaxD, while the low intratumoral dose suffers some DC loss from the tumor. This can explain why the low-dose intravenous injection is more effective than the low-dose intratumoral injection. However, in the high-dose injection cases, the number of DCs entering the spleen from the high-dose intravenous treatment is being limited by the MaxD inflow ceiling. On the other hand, the high-dose intratumorally injected DCs enter the spleen more slowly, so even though there is still loss from the tumor, a greater total number of DCs remain in the spleen longer than in the case of the high-dose intravenous injection.

As discussed earlier, the model does not explicitly account for the distance between the spleen and the tumor. If we were to extend the model so that it could apply to human subjects, these distances could vary significantly between individuals. One possible model extension would incorporate the effect of transit times with a partial differential equation that includes a distance L along which cells diffuse and convect.

4.2. Modified Dosing

In this section we explore the hypothetical effect of modifying dose timings. We first investigate the effect of administering the same total medication over a 20 day time frame, but with more frequent injections; that is, a fractionated dosing schedule. The original dosing schedule starts on day 6, and administers a total of three doses spaced apart by 2 days. The hypothetical fractionated dosing schedule we explore also starts on day 6, but administers doses twice a day at 1/4 the original dose.

In Figure 4, we compare simulated melanoma growth in response to DCs administered according to the original protocol with a hypothetical fractionated dosing schedule. Injections are given intratumorally. The no-treatment tumor growth case is also included for comparison. After a tumor challenge of 5 × 10⁵ B16F10 melanoma cells on day 1, the original dosing calls for DC injections of 1 × 10⁵, 7 × 10⁵, and 21 × 10⁵ administered every other day on days 6, 8, 10. The hypothetical fractionated schedule administers doses of 0.25 × 10⁵, 1.75 × 10⁵, and 5.25 × 10⁵ twice a day on days 6, 7, 8, 9, 10, 11. The total DC treatment administered is the same in the original and fractionated dosing scenarios. It is clear that the fractionated schedule does not improve outcomes in the case of intratumoral injections.

In Figure 5, we again compare the original DC dosing schedule to a fractionated dosing schedule, but we now use intravenous injections. As before, after a tumor challenge of 5 × 10⁵ B16F10 melanoma cells on day 1, the original treatment schedule calls for DCs of doses 1 × 10⁵, 7 × 10⁵, and 21 × 10⁵ given every other day on days 6, 8, 10. The hypothetical fractionated schedule administers doses of 0.25 × 10⁵, 1.75 × 10⁵, and 5.25 × 10⁵ twice a day on days 6, 7, 8, 9, 10, 11. The total amount of DC administered is the same in both scenarios. The simulations highlight that although the fractionated schedule does not improve outcomes in the case of intratumoral injections, greater tumor control is observed when the fractionated treatment is administered intravenously. Although tumor growth is slowed with the intravenously dosed fractionated schedule, it is not completely controlled, and the tumor still eventually grows.

We next investigate the effect of starting the DC treatment regimen earlier than day 6. In this case, we compare fractionated DC doses both intratumorally and intravenously, but with treatment initiated on day 3 instead of day 6. The experimental outcomes are pictured in Figure 6. After a tumor challenge of 5 × 10⁵ B16F10 melanoma cells on day 1, DC injections of 1 × 10⁵, 7 × 10⁵, and 21 × 10⁵ are administered both intratumorally and intravenously on days 3, 5, 7 (pictured in graphs in the top row of Figure 6). The hypothetical fractionated schedule (pictured in the second row of Figure 6) administers doses of 0.25 × 10⁵, 1.75 × 10⁵, and 5.25 × 10⁵ twice a day on days 3, 4, 5, 6, 7, 8. The total DC treatment administered is the same in all scenarios. We see that in both the intratumoral and intravenous dosing cases, tumor growth is slowed when treatment starts on day 3. However, as before, fractionating the intratumoral doses does not have much effect (bottom left panel), but does slow tumor growth even further when administered intravenously (bottom right panel). Interestingly, the earlier start day has less effect when administered intravenously than it does when administered intratumorally, as can be seen in Figure 7. Here we compare non-fractionated intravenous dosing starting on day 6 (left panel) and on day 3 (right panel). Initial values and doses follow the original schedule. The result with fractionated dosing is similar, but is not pictured. We see that there is some improvement with the early start intravenous dosing, but the improvement is not as large as the improvement seen with the intratumoral doses started on day 3, as pictured in the left column of Figure 6.

FIGURE 6

Figure 6. Early treatment initiation: start day 3. Original dosing (top row) and fractionated dosing (bottom row) compared. 5 × 10⁵ B16F10 melanoma cells on day 1, DC injections of 7 × 10⁵, and 21 × 10⁵ are administered both intratumorally and intravenously 0.25 × 10⁵, 1.75 × 10⁵, and 5.25 × 10⁵ twice a day on days 3, 4, 5, 6, 7, 8. Treatment is given both intratumorally (left column) and intravenously (right column).

FIGURE 7

Figure 7. Compare original to earlier treatment initiation, intravenous, non-fractionated dosing. Treatment start day 6 (left) compared to start day 3 (right).

4.3. Prophylactic Vaccination

So far, we have been investigating the responses of our system to DC treatment after a tumor challenge. However, DC treatments have also been considered to have potential as prophylactic vaccines. For example, the work of Preynat-Seauve et al. (2007) details a variety of studies on tumor growth in mice inoculated with DC treatments prior to a tumor challenge. Although our model has not been constructed specifically to investigate preventative vaccination, we did see some interesting results when simulating such treatment. In one of the experiments from Preynat-Seauve et al. (2007), DC cells from tumors were cultivated and injected into B6C3F1 mice. Vaccination was performed twice, once weekly, with 10⁵ tumor-infiltrating dendritic cells. The authors state that this number corresponds to the total number of CD11c+ cells recovered from a single 1 cm-diameter tumor. Two weeks after the last injection, mice were challenged with 2 × 10⁵ melanoma cells (either K1735 or B16F10). According to the study, vaccinated mice were protected for 22 days, whereas naive mice succumbed to the tumor challenge.

Using the same parameter values we determined through fitting to the data in Lee et al. (2007), a simulation of pre-vaccination with mature DCs showed no particular benefit. In order to determine the sensitivity of the system to a change in parameters we used the Latin Hypercube sampling method described in Blower and Dowlatabadi (1994) to compare simulations with 50,000 randomly generated parameter sets. The effect of a change in parameter values on tumor size and CTL levels after 45 days was quantified by calculating the partially ranked correlation coefficients for each parameter that showed a monotonic relationship to the outcomes. The results are shown in Figure 8.

FIGURE 8

Figure 8. Partially Ranked Correlation Coefficients (PRCC) for two outcomes: final tumor size, and final CTL levels. Parameters with negative PRCCs relative to final tumor levels are negatively correlated with tumor growth. Increasing such parameters would be beneficial to the patient. An increase in parameters with positive PRCCs relative to final tumor levels, or with negative PRCCs relative to final CTL levels could be harmful to the patient.

The sensitivity analysis indicates that tumor levels are sensitive to the parameters d, m, and μ_TB. Since d, the fractional tumor kill rate by CTLs, has the potential to be manipulated through treatments (c.f. Chakraborty et al., 2003), we suggest that this parameter might play an important role in the vaccine’s effectiveness. In subsequent simulations, when d was increased we observed a protective effect of prophylactic vaccination with DCs. In Figure 9 we see tumor growth both without and with DC vaccination with d = 0.35, 0.85, 1.0, and 1.25. The results for d = 1 are what interest us. When d = 1.25, the immune system is sufficiently effective to suppress tumor growth without treatment intervention. In contrast, when d = 1 or less, then without DC vaccination, the tumor grows rapidly. However, if the tumor challenge has been preceded by two doses of a DC vaccine, then even when d = 1, tumor growth is suppressed. We also extended the simulations out 165 days after the tumor challenge (not pictured). For d ≥ 1, the tumor shrank to zero after vaccination and did not regrow. Numerically, this indicates that the zero tumor equilibrium is stable.

FIGURE 9

Figure 9. Prophylactic vaccination and the effect of varying immune strength parameter d. Top panel, no vaccine. Bottom panel, vaccinate with DC treatments, days 0 and 7, 1 × 10⁵ DCs per dose. Tumor challenge on day 21, with 2 × 10⁵ tumor cells. Dosing follows Preynat-Seauve et al. (2007) experiment. With d = 1, tumor is controlled as a result of vaccination.

5. Future Directions

In this paper, we presented a model of Dendritic Cell trafficking and interaction with a tumor cell population. With this model, we achieved simulation outcomes that quantitatively match published data from studies on mice (Lee et al., 2007) that first were challenged with tumor and subsequently treated with DC therapy. We then used the model to test a wider variety of hypothetical treatment scenarios. In addition, we examined the effects of prophylactic vaccination with DCs. The simulation results from the prophylactic vaccination scenarios that were discussed in the previous section are preliminary, but do show qualitative agreement with a different set of data from published laboratory experiments on mice (Preynat-Seauve et al., 2007).

In future work, we will investigate how to scale this model to reflect tumor growth and DC trafficking in humans. This will involve a careful examination of the effect of distances between the tumor site and the lymph organs. In addition, we will further investigate the effects of prophylactic vaccination. Our goal is to determine which parameters and model terms need tuning to achieve quantitative as well as qualitative outcomes that reflect the laboratory data.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Appendix

A. Parameter Values

The parameters are described in Table A1. Refer to Ludewig et al. (2004) for justification and confidence intervals for the parameters measured in that study.

TABLE A1

Table A1. Parameter values.

B. Stability Analysis: Jacobians

Recall that J₀ is the Jacobian for the non-delayed system, in other words, its entries are given by J₀(i, j) = dF_i/dx_j, and J_τ(i, j) = dF_i/dz_j is the Jacobian for the delayed system. For computational ease, we rename the state variables as follows:

$\begin{array}{c}{D}_{blood}=x_1,D_{spleen}=x_2,E_{blood}^a=x_3,\\ E_{spleen}^a=x_4,E_{blood}^m=x_5,E_{spleen}^m=x_6,\\ E_{tumor}^a=x_7,T=x_8,D_{tumor}=x_9.\end{array}$

There are only two delayed variables:

$D_{spleen}\left(t-\tau \right)=z_2,E_{spleen}^a\left(t-\tau \right)=z_4.$

We perform the stability analysis without any treatment, i.e., v_blood = v_tumor = 0. The non-zero elements of J₀ are calculated to be:

$\begin{array}{l}{J}_0(1,1)=-{\mu }_B\\ J_0(1,9)={\mu }_{TB}\\ J_0(2,1)=-{\mu }_{BS}{e}^{\left(\frac{-{\mu }_{BS}{x}_1}{MaxD}\right)}\\ J_0(2,2)=-a_D-b_{DE}{x}_4\\ J_0(2,4)=-b_{DE}{x}_2\\ J_0(3,2)=\frac{-\Delta \mu }{{\theta }_{shut}}\frac{x_4}{{(1+x_2/{\theta }_{shut})}^2}\\ J_0(3,3)=-{\mu }_{BB},\\ J_0(3,4)={\mu }_{SB}(D_{spleen}),\\ J_0(4,2)=\frac{-\Delta \mu }{{\theta }_{shut}}\frac{x_4}{{(1+x_2/{\theta }_{shut})}^2}+b_a{x}_6\\ J_0(4,3)={\mu }_{BSE}\\ J_0(4,4)={\mu }_{SB}(D_{spleen})-a_{E_{a}S}-r_{am}\\ J_0(4,6)=b_a{x}_2\\ J_0(5,2)=\frac{-\Delta \mu }{{\theta }_{shut}}\frac{x_6}{{(1+x_2/{\theta }_{shut})}^2}\\ J_0(5,5)=-{\mu }_{BB},J_0(5,6)={\mu }_{SB}(D_{spleen})\\ J_0(6,2)=-\left(b_a-\frac{\Delta \mu }{{\theta }_{shut}}\frac{1}{{(1+x_2/{\theta }_{shut})}^2}\right)x_6\\ J_0(6,4)=r_{am}\\ J_0(6,5)={\mu }_{BSE}\\ J_0(6,6)=-(a_{E_m}+b_a{x}_2+{\mu }_{SB}(D_{spleen}))\\ J_0(7,3)={\mu }_{BTE}(T)\\ J_0(7,7)=-a_{E_{a}T}-c{x}_8\\ J_0(7,8)=\frac{{\mu }_{BB}\alpha }{{(\alpha +x_8)}^2}{x}_3-c{x}_7\\ J_0(8,7)=\frac{dsl{(x_7/x_8)}^{l-1}}{{\left\{s+{(x_7/x_8)}^l\right\}}^2}\\ J_0(8,8)=r-\frac{2r{x}_8}{k}-\frac{dsl{(x_7/x_8)}^l}{{\left\{s+{(x_7/x_8)}^l\right\}}^2}+\frac{d{(x_7/x_8)}^l}{s+{(x_7/x_8)}^l}\\ J_0(9,8)=\frac{mq}{{(q+x_8)}^2}\\ J_0(9,9)=-({\mu }_{TB}+a_D)\end{array}$

where ${\mu }_{BTE}\left(T\right)={\mu }_{B}B\frac{T}{\alpha +T},and{\mu }_{SB}\left(D_{spleen}\right)={\mu }_{SB}^{*}+\frac{\Delta \mu }{1+D_{spleen}∕{\theta }_{shut}}$.

There are only two non-zero elements of J_τ:

$J_{\tau }\left(5,4\right)=b_p\frac{z_5{\theta }_D}{{\left({\theta }_D+z_4\right)}^2},J_{\tau }\left(5,5\right)=\frac{b_p{x}_4}{{\theta }_D+z_4}.$

As referenced in the text, there are two positive tumor equilibrium values, T*:

$T_1^{*}\approx 9.72\times 10^8;T_2^{*}\approx 9.81\times 10^8.$

However, the first one results in a non-biologically relevant equilibria as it corresponds to a negative equilibrium value for x₂, x₄, x₅. For the non-delay case (i.e., τ = 0,) there are no J_τ elements. Thus we consider the eigenvalues of $J_0$. At the positive (biologically relevant) equilibria, the eigenvalues of J_0, as computed via MatLab to four decimal places, are

$\begin{array}{l}\overrightarrow{e}=[-27.0720,-11.1021,-0.6765,-0.4719,-0.3918,\\ -0.3092,-0.2321,-5.7000,-5.7000].\end{array}$

Each eigenvalue clearly satisfies ℛ(λ) < 0, indicating that, in the non-delay case, the biologically relevant equilibrium is stable.

Acknowledgments

Ami Radunskaya is partially supported by NSF grant DMS-1016136.

References

Blower, S. M., and Dowlatabadi, H. (1994). Sensitivity and uncertainty analysis of complex models of disease transmission: an HIV model, as an example. Int. Stat. Rev. 62, 229–243.

CrossRef Full Text

Boon, T., Coulie, P. G., Van den Eynde, B. J., and van der Bruggen, P. (2006). Human T cell response against melanoma. Annu. Rev. Immunol. 24, 175–208.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cappuccio, A., Elishmereni, M., and Agur, Z. (2006). Cancer immunotherapy by interleukin-21: potential treatment strategies evaluated in a mathematical model. Cancer Res. 66, 7293–7300.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chakraborty, M., Abrams, S. I., Camphausen, K., Liu, K., Scott, T., Norman Coleman, C., et al. (2003). Irradiation of tumor cells up-regulates Fas and enhances CTL lytic activity and CTL adoptive immunotherapy. J. Immunol. 170, 6338–6347.

Pubmed Abstract | Pubmed Full Text

Cheever, M. A. (2011). PROVENGE (sipuleucel-T) in prostate cancer: the first FDA-approved therapeutic cancer vaccine. Clin. Cancer Res. 17, 3520–3526.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

de Pillis, L. G., Fister, K. R., Gu, W., Collins, C., Daub, M., Gross, D., et al. (2009). Mathematical model creation for cancer chemo-immunotherapy. Comput. Math. Methods Med. 10, 165–184.

CrossRef Full Text

de Pillis, L. G., Gu, W., and Radunskaya, A. E. (2006). Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations. J. Theor. Biol. 238, 841–862.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

de Pillis, L. G., and Radunskaya, A. (2003). “A mathematical model of immune response to tumor invasion,” in Computational Fluid and Solid Mechanics 2003, ed. K. J. Bathe (Elsevier Science Ltd), 1661–1668.

de Pillis, L. G., Radunskaya, A. E., and Wiseman, C. L. (2005). A validated mathematical model of cell-mediated immune response to tumor growth. Cancer Res. 65, 7950–7958.

Pubmed Abstract | Pubmed Full Text

de Pillis, L. G., Radunskaya, A. E., and Wiseman, C. L. (2007). Comment on: a validated mathematical model of cell-mediated immune response to tumor growth. Cancer Res. 67, 8420.

CrossRef Full Text

DeMatos, P., Abdel-Wahab, Z., Vervaert, C., Hester, D., and Seigler, H. (1998). Pulsing of dendritic cells with cell lysates from either B16 melanoma or MCA-106 fibrosarcoma yields equally effective vaccines against B16 tumors in mice. J. Surg. Oncol. 68, 79–91.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Diefenbach, A., Jensen, E. R., Jamieson, A. M., and Raulet, D. H. (2001). Rae1 and H60 ligands of the NKG2D receptor stimulate tumour immunity. Nature 413, 165–171.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fields, R. C., Shimizu, K., and Mule, J. J. (1998). Murine dendritic cells pulsed with whole tumor lysates mediate potent antitumor immune responses in vitro and in vivo. Proc. Natl. Acad. Sci. U.S.A. 95, 9482–9487.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lee, T.-H., Cho, Y.-H., and Lee, M.-G. (2007). Larger numbers of immature dendritic cells augment an anti-tumor effect against established murine melanoma cells. Biotechnol. Lett. 29, 351–357.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ludewig, B. B., Krebs, P., Junt, T., Metters, H., Ford, N. J., Anderson, R. M., et al. (2004). Determining control parameters for dendritic cell-cytotoxic t lymphocyte interaction. Eur. J. Immunol. 34, 2407–2418.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pilon-Thomas, S., Verhaegen, M. E., and Mulé, J. J. (2004). Dendritic cell-based therapeutics for breast cancer. Heart Dis. 20, 65–71.

Preynat-Seauve, O., Contassot, E., Schuler, P., French, L. E., and Huard, B. (2007). Melanoma-infiltrating dendritic cells induce protective antitumor responses mediated by t cells. Melanoma Res. 17, 169–176.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shinagawa, N., Yamazaki, K., Tamura, Y., Imai, A., Kikuchi, E., Yokouchi, H., et al. (2008). Immunotherapy with dendritic cells pulsed with tumor-derived gp96 against murine lung cancer is effective through immune response of CD8+ cytotoxic T lymphocytes and natural killer cells. Cancer Immunol. Immunother. 57, 165–174.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Taquet, N., Roberts, L. K., Burkeholder, S., Phipps, P., Finholt, J., Walters, L., et al. (2008). Dendritic cell vaccine production facility: from design to operation. Bioprocess. J. 7, 28–33.

Trefzer, U., Herberth, G., Wohlan, K., Milling, A., Thiemann, M., Sharav, T., et al. (2005). Tumour-dendritic hybrid cell vaccination for the treatment of patients with malignant melanoma: immunological effects and clinical results. Vaccine 23, 2367–2373.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yamaguchi, S., Tatsumi, T., Takehara, T., Sakamori, R., Uemura, A., Mizushima, T., et al. (2007). Immunotherapy of murine colon cancer using receptor tyrosine kinase epha2-derived peptide-pulsed dendritic cell vaccines. Cancer. 110, 1469–1477.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text