Immuno-kinetics of immunotherapy: dosing with DCs

Therapeutic vaccines play a large role in the cast of immunotherapies that are now an essential component in most cancer treatment regimes. The complexity of the immune response and the ability of the tumour to mount a counter-offensive to this response have made it difficult to predict who will respond to what treatments, and for clinicians to optimise treatment strategies for individual patients. In this paper, we present a mathematical model that captures the dynamics of the adaptive response to an autologous whole-cell cancer vaccine, without some of the complexities of previous models that incorporate delays. Model simulations are compared to published experimental and clinical data, and used to discuss possible improvements to vaccine design.


Introduction
Immunotherapy was heralded as the 'Breakthrough' of the year in 2013, and laboratory and clinical trials show a variety of treatments that have great promise (Couzin-Frankel, 2013). Cancer vaccines are a particular type of immunotherapy, where the goal is to boost the patient's own immune response to tumour cells. Unlike preventative vaccines which are developed for a population and target a specific virus, a cancer vaccine is therapeutic: it is developed for a particular patient for their specific tumour. The complexity of the immune response, the heterogeneity of the many diseases we call 'cancer' and the variability between patients create challenges in designing an effective cancer vaccine. Mathematical models can help oncologists meet these challenges by (1) increasing the understanding of the mechanism underlying the observed immuno-kinetics; (2) identifying predictors of disease progression and response to therapies; (3) providing an in silico test bed for treatment scenarios and virtual clinical trials; (4) suggesting combinations of therapies and dosage timings to improve the effectiveness of the treatments. There are multiple kinds of immunotherapy used in fighting cancer Institute (2015), which include monoclonal antibodies, adoptive cell transfer and therapeutic vaccines. The greatest advantage of immunotherapy over chemotherapy and radiotherapy is its low toxicity and potential long-term benefits. In this paper, we focus on modelling the adaptive immune response to a dendritic cell vaccine. Dendritic cells (DCs) are a type of antigenpresenting cell. The role of these immune cells is to process foreign substances, called 'antigen', found in the body, and present the antigen to immune cells that can be activated to engage in an immune response specific to that antigen. Immature DCs circulate through the blood stream, or hang out in tissues that are close to the outside of the body, such as the skin, nose and lungs. Once they encounter antigen and process it, they become 'activated', and they migrate to the lymph organs, such as the spleen, where they present the antigen to immature T-cells. (See Figure 1).
A dendritic cell vaccine is created from the patient's own immune cells. Immature DCs are isolated from the patient's blood, and allowed to mature and proliferate in contact with the patient's tumour cells. Activated DCs are then injected back into the patient, with the goal of initiating a specific immune response by activating cytotoxic T lymphocytes (CTLs), a type of effector cell. These activated CTLs then go and fight the tumour, become memory cells to be available for a future tumour challenge, or die off (Porfyris & Kalomoiris, 2013).
The first DC vaccine to be approved by the FDA was for prostate cancer, and a large clinical trial started in 2010. Prostate cancer is generally slow-growing, which may allow the immune system to mount an effective anti-tumour response after stimulation. While the vaccine prolonged survival for several months, the results were not as positive as was hoped (Small et al., 2006). In this paper, we develop a mathematical model to explore some of the mechanisms behind the body's ability to mount an effective response.
In §2, we develop a mathematical model of the CTL activation in the spleen. This model improves on earlier models by explicitly tracking the DC-CTL conjugates during the activation phase. This avoids the need for delay equations and makes the model more descriptive of what is actually happening. In §3.1, we calibrate the model to vaccination experiments with mice and humans. In §3.2, a sensitivity analysis shows which model parameters most effect outcomes of interest, such as final tumour size, final CTL levels and memory cell levels. In §4, we apply the model to several situations that have been explored experimentally, and use model simulations to inform phenomena such as variability in response and immuno-exhaustion.

The mathematical model
Previous mathematical models DePillis, Gallegos, and Radunskaya (2013), Ludewig et al. (2004), and Radunskaya and Hook (2012) describe the number of dendritic and effector cells in the spleen, blood and/or tumour. Ludewig, et al. describe mathematical models for the predator-prey type of interaction between DC and CTLs, and transfer dynamics of the DC and activated, naive and memory CTLs in the spleen, blood and the liver with a nonlinear compartmental system (Ludewig et al., 2004). The model includes a DC-dependent retention of CTLs in the spleen, and a 'trapping effect' which controls the export rate of CTLs into the blood depending on the presence of DCs in the spleen. The trapping effect has been incorporated in order to comply with the observed CTL dynamics. DePillis, et al. extend this model by incorporating the tumour compartment, which in turn introduces the tumour population, and its interactions with the immune cells. A contribution of this model is the incorporation of the effects of DC-based vaccines on tumour growth. The tumour growth is modelled by logistic growth, as it was a better fit to the data, and the elimination of tumour cells by the activated CTLs is described by a ratio-dependent kill term.
When the activated DCs encounter the CTL cells, they bind together via receptors on the cell surface. There is a time lapse, the 'synaptic connection time', before the CTLs are completely activated, and can begin proliferating. In DePillis et al. (2013) and Ludewig et al. (2004) delay differential equations are used to describe the synaptic connection time. The delay equations account for the time lag, but they introduce unwanted complexities into the model. In particular, the presence of the delay means that a continuum of states is required to describe the system at any point in time: the model is infinite dimensional. In this work, we remove the delay by introducing two new state variables which represent the DCs bound to CTLs (W) and the rapidly proliferating population (P). In addition to reducing the model to a finite-dimensional system, we can also directly model this important part of the immune activation cascade.
We consider three locations for the cells: the blood, spleen and tumour. Within each location, we model the dynamics of the immune and tumour cells. The events that occur within each location and the connection between them are described in Figure 2. The blood primarily serves as the transporter of dendritic and effector cells between the spleen and tumour. In the spleen, the activation and proliferation of effector T-cells occurs. During the activation of the effector T-cells, some proportion of DCs enter a new state of bonded cells,  Notes: The parameters indicate the rate at which a cell moves, dies or changes; see Appendix A for details. The 'no entry' symbol denotes a death term. The dotted line indicates an interaction term. referred to as the synaptic connection state. After the synaptic connection time, effector cells begin proliferating and produce active effector cells. The effector cells either move to their target, turn into memory cells that remain in circulation until the next challenge by the same antigen, or enter an apoptotic or rapidly dying phase. This proliferation of cells better equips the immune system to fend off tumour cells.
The spleen is populated by dendritic (D spleen ) and effector cells. The effector cells are classified as naive, active E a spleen or memory E m spleen . There is an abundance of naive cells so we do not incorporate these into the model. The DCs bind with up to 10 naive effector cells to form a conglomerate cell W (Giese & Marx, 2014). When the conglomerate breaks apart (after the required activation time), the effector cells then enter a proliferating stage P where new active effector cells E a spleen are created. Memory effector cells can also bind with DCs and initiate a proliferating stage. Figure 3 shows the different stages and the transfer rate between those stages. The parameter μ represents the flow rate for each type of cell in and out of the spleen, a represents the death rates, and λ, β and γ are the conversions between types of cells (see Appendix A). The DCs of interest are those that recognize the tumour antigen and move to the tumour. Once a DC has arrived at the tumour, it has reached its destination thus we do not track it leaving the tumour. The rate at which effector cells are released into the blood (μ SB ) is lower when DCs are present in the spleen, describing the 'trapping' effect noted by some experimentalists (Ludewig et al., 2004).
The resulting model for the interaction between the DCs and the effector cells in the spleen is given in system (1).
The blood transports the dendritic and effector cells throughout the body. We assume there is no proliferation or cell death in the blood compartment, resulting in system (2). One treatment option for reducing the size of the tumour is to introduce additional DCs into the blood. This vaccine is composed of immature DCs harvested from the patient, and activated ex-vivo. The mature DCs are then injected into the bloodstream, described by v blood (t) in Equation (2a). We investigate the dosage and the frequency of the vaccine in Section 4. The DCs emigrate from the blood at a rate μ B , and enter the blood from the tumour at a rate of μ TB . Both the active and memory effector cells are released from the spleen dependent on the DC density, as in the Equation (1c). The function μ * SB + μ normal −μ *

1+
D spleen θ shut describes the effect that activated CTLs are retained in the spleen when DCs are present, this term is called the 'trapping term' in DePillis et al. (2013), Ludewig et al. (2004). These equations modelling the transport of cells in the blood are the same as those presented for the blood in DePillis et al. (2013).
The tumour in the model is composed of DCs, active effector cells and tumour cells. The number of effector cells and DCs travelling to the tumour is dependent on the tumour size, with maximal rates of μ BB and m, respectively. We assume a logistic growth for the tumour. Every encounter between an effector cell and tumour cell results in the annihilation or inactivation of the effector cell with probability c. The corresponding elimination rate for the tumour is more complicated due to the specific immune strength of the patient. We incorporate a ratio-dependent kill term, in Equation (4), based on agreement with previous experimental results, (see de Pillis, Radunskaya, & Wiseman, 2005). The effector cells do not leave the tumour compartment, unlike the DCs. In this work, we investigate the tumour as an alternative application site for the dendritic cell vaccine, hence the appearance of the term v tumor (t). We explore the tumour growth for various dosages and frequencies of the vaccine in Section 4.
The equations modelling the interaction of cells within the tumour are the same as those in DePillis et al. (2013).
We use the non-dimensional version of the model (1-3) when simulating. Details of the non-dimensionalization are given in Appendix B.

Parameters
In this section, we describe the algorithm used to fit the model to data. We also determine the sensitivity of the model to the parameters. All parameter descriptions can be found in Appendix A.

Estimating parameters
The parameters in the model are fit using the Nelder-Mead least squares algorithm by implementing the MATLAB function fminbdd. The data used are from an experimental study on mice Lee, Cho, and Lee (2007) who received an injection of PBS (control), 1 × 10 5 , 7 × 10 5 , or 21 × 10 5 DCs on days 6, 8, and 10. In these experiments, the vaccine was injected directly into the tumour. In the numerical fitting, we consider experimental day 6 to be simulation day 0. The initial size of the tumour is 5 × 10 5 cells, assuming the average side length of a square tumour cell is .01 mm. The resulting values for the parameters fit are given in Table A1 and the fit is shown in Figure 4. The model is able to mimic tumour growth under all four sets of experimental conditions with one set of parameter values. The exception is the tumour size on Day 14 using the highest dose of vaccine: while both the experimental and simulated tumours continue to grow over time, the mathematical model predicts a slightly smaller tumour size than that observed experimentally. We conjecture that this is due to the phenomenon known as 'immunoexhaustion', where increasing doses of vaccine have a decreasing ability to promote an effective immune response. This is discussed further in Section 4.4, where we suggest model refinements that include immuno-exhaustion.

Sensitivity analysis
In order to examine the sensitivity of the system to the choice of parameters, we implement the Latin Hypercube Sampling technique described in Blower and Dowlatabadi (1994) and calculate partially ranked correlation coefficients (PRCC). The parameters that have a monotonic relationship to the outcomes are a EaT , d, k, , μ B , μ BS , m, r, s. We consider the outcomes on day 15 after receiving three doses of vaccine consisting of 10 7 DCs on days 1, 3, and 5. Assuming each of the parameters are normally distributed, we randomly generate 1,00,000 sets of parameters selected from a normal distribution truncated at two standard deviations. The maximum and minimum values for each distribution are given in Table 1. Note: Parameter values are given in Table A1.
The sensitivity is quantified by calculating the PRCC. The results are shown in Figure 5. The parameters that have the greatest influence on the tumour size include d (negative), , r and s (positive). The parameters with the largest effect on the CTL and memory levels are μ B (negative), μ BS and m (positive). The p-value for all of these parameters is less than .001.
In order to further investigate the influence of these parameters on the size of the tumour, we vary several parameters that can potentially be regulated by clinical intervention. The parameters determined when fitting to the data in Lee et al. (2007) are used with a vaccination dosage of 7 × 10 5 DCs given on days 1, 2, 4 for 30 min. As d, r, s individually vary, Figure 6 shows the change of the outcome of the tumour. In all cases, there are parameters that cause the tumour to decrease and eventually disappear. These numerical experiments indicate bifurcation values for each of these parameter values. The clinical implication is that treatments need to move the effective parameter values across these bifurcation values.
We also explore the relationship between the death rate of the tumour cells by the CTLs (d) and the tumour growth rate (r). Figure 7 shows the needed r value for a given d value in order for the tumour to be eliminated when the medium dosing schedule (10 7 DCs given on days 1, 3, and 5) is implemented. The allowable tumour growth rate increases at a slower rate than the death rate of the tumour cells. Thus an increase in the growth rate of the tumour requires a larger increase in the death rate of the tumor cells by CTLs.
We do not further investigate the transfer rates of the DCs to which the memory cells and CTLs are most sensitive because we do not experience drastically different behaviour due to the low number of memory cells within the system. One would expect that the carrying capacity (k) would greatly influence the growth of the tumour. However, at the early time of 15 days the tumour size is not near carrying capacity which explains k not being one of the most sensitive parameters.

Results
The mathematical model developed in the previous sections gives a simple description of the kinetics of the adaptive immune response to a dendritic cell vaccine. Its potential usefulness is twofold: (1) it can be used to explore hypothetical mechanisms governing the immune response and the effectiveness; (2) it can be used to predict responses to hypothetical vaccination strategies. We calibrated the model to four different published experiments in order to test its explanatory and predictive power. These include three murine experiments, and one study on patients receiving a dendritic cell vaccine for prostate cancer. We note also that the experiments were done using different tumour lines, and using a variety of treatment protocols. The results reported in this section illustrate both the flexibility and the limitations of the current model.

Experiment 1: a murine model of a DC vaccine for ovarian cancer
In Shubina, Donenko, Akhmatova, and Kiselevskii (2014), an experiment was conducted on 8 − 12 week old male mice. The DC vaccine was created by removing DCs from mice  Note: Parameter combinations that lie above the curve result in the tumour growing to carrying capacity while those below the curve cause the tumour to be eliminated. bone marrow and exposing them to a type of ovarian cancer (CaO-1). The size of the tumour was converted to number of tumour cells by assuming the cells are circular with an average diameter of 17.6 microns (LLC, 2003). In the experiment, each mouse has 5, 00, 000 ovarian cancer tumour cells put intraperitoneally into the right side and 10 6 DCs into the left side. Due to this change in sides, we introduce the vaccine into the blood on days 1, 3, 5, 7, and 10. Figure 8 shows the agreement between the data and the simulation. Notice that the tumour continues to grow in both groups: the growth rate is slower when a DC vaccine is administered. On day 20, the tumour is 54% the size of the control and on day 30 the tumour is 58% the size of the control (Shubina et al., 2014). Using a least-squares fit as described in Section 3.1, we fit the model to the data above by adjusting the parameters given in Table 2. The resulting simulations are show in Figure 8.

Experiment 2: a study of a DC vaccine for patients with prostate cancer
The effects of dendritic cell therapies for humans with prostate cancer were studied in Profile (2006). We utilise the linear function found in Portz, Kuang, and Nagy (2012) to convert the given PSA numbers to tumour cells. The value for r, the intrinsic growth rate of the tumour cells, given in Portz et al. (2012) is .025 − .045 day, so we set r = .031, in the middle of this range. We fit our model to the results in Profile (2006) by adjusting the other parameters as given in Table 3. Figure 9 shows the tumour doubling time for patients receiving the control and patients receiving the vaccine. The time required for the tumour to double is longer when DC vaccines are employed. The model successfully captures the growth rate of the tumour.  Figure 9. Data from Profile (2006) shown with model simulations. In Profile (2006), it is reported that the doubling time increased from 6.7 months with no vaccine to 12.7 months with vaccine.
Notes: Patients with prostate cancer were administered a 30 min IV-infusion of .3 − 1.2 × 10 9 cells of the vaccine on weeks 0, 2, and 4. We utilize .8 × 10 9 cells of DC vaccine injected into the blood over a 30 min period.
We further investigated the variable, d, an immune strength parameter that represents the maximum per-cell kill rate of tumour cells by effector cells. Figure 10 shows that when d = .03 the tumour grows very slowly; above this value, the tumour grows and below this value, the tumour dies off. This suggests that this dendritic cell vaccine for prostate cancer could be made more effective by increasing the maximum per-cell kill rate of the tumour cells by the effector cells. One possibility in achieving this are the anti-PDL1 antibodies currently being developed, (e.g. Sagiv- Barfi et al., 2015). See also Section 5 for further discussion of this type of treatment.

Experiment 3: a murine model of a vaccine for prostate cancer
The sipuleucel-T vaccine modelled in Section 4.2 consisted of mature DCs, as well as some proteins, or cytokines, to promote the activation of cancer-specific immune cells. These cytokines were fused with prostatic acid phosphatase (PAP), a good target for vaccines against prostate-specific disease since it is found in most prostate cancers. In the hopes of improving the efficacy of the sipuleucel-T vaccine, researchers in Fujio et al. (2015) added additional cytokines to the treatment. They tested the new vaccine on mice who were injected with PAP-RM9 cells, i.e. tumour cells that should present a sensitive target for these PAP-fused cytokines.
The in vivo experiments in Fujio et al. (2015) were designed to test the effect of the DC vaccine on the very early stages of tumour growth. Therefore, they gave three vaccine injections, spaced two days apart (on days 0, 2, and 4), before injecting the tumour cells (on day 7). A control group was also injected with tumour cells, but received no vaccine. All  Table 3. mice were analysed and examined for tumour formation and tumour size on the twentyfirst day.
Antitumour effects of the treatment were evidenced in two ways: only 40% of the mice in the vaccinated group formed tumours at all (2 out of the 5 mice), and the tumours that did occur in the vaccinated group were 20% of the size of the tumours in the group with no vaccination. We tested our model against data from the control group and data from the vaccinated group. In order to capture the heterogeneity in the experimental groups, we sampled from a range of prostate tumour growth rates reported in the literature (Portz et al., 2012). We hypothesised that the additional cytokines, which included human interleukin-2 (IL2), IL4 and IL7, should increase the lytic efficacy of the effector cells. By adjusting the corresponding parameter, d, in the model, we were able to fit the model to the experimental results (see Figure 11). As hypothesised, the value of d increased to .2, an order-of-magnitude increase from the fit in the previous numerical experiment reported in Section 4.2.
This numerical experiment underlines another use for mathematical modelling. If a range of parameter values is known from previous experiments, in silico 'trials' can be performed by sampling from this parameter distribution. Thus, we can use the model to not only predict average behaviour, but also to show the range of behaviours that might occur when a new treatment is applied. For a more extensive example of this type of modelling, (see de Pillis, Radunskaya, & Savage, 2014).

A murine model of a DC vaccine to study immuno-exhaustion.
When used in isolation, therapeutic vaccines are not a consistently effective treatment. One obstacle to designing an effective vaccine is that the immune system's response drops off after continued treatment. Two causes of this response may be that the DCs cease to Results are consistent with experimental results reported in Fujio et al. (2015): tumours grew in all mice in the un-vaccinated group, but only two of the five mice in the vaccinated group showed any tumour growth after fifty days. The size of the tumours in these two mice was only 40% of the average tumour size in the unvaccinated group. proliferate into active effector cells but instead become ineffective or the memory cells cease to begin proliferating (see Jacobs et al., 2014 andRicupito et al., 2013). For example, a decrease in the total number of effector cells in observed in the experiments performed in Francesco Pappalardo, Ricupito, Topputo, and Bellone (2014). We adapt the model (1-3) to incorporate a drop off in active effector cells (E a spleen ). The relevant parameters are the reversion rate of the activated CTL to memory CTL (r am ), the rate at which proliferating cells become ineffective (they do not produce active CTL) (β d ), and the rate at which memory cells begin proliferating (λ m ). First we dramatically increase r am to get a higher per cent of memory effector cells (E m spleen ). Second we let β d and λ m change with time. We would like the rate at which proliferating cells become ineffective to dramatically increase at first and then level off in time (see Figure 12). At the same time, we want the probability that memory cells are activated and begin proliferating to decrease over time. While a more sudden decrease is more likely to be biologically accurate, we start with a linear decrease as a simple first model (see Figure 12). Finally, since there is no tumour present, we set the relevant transfer rates, μ TB and μ BB , to zero. With these adjustments to the parameters as seen in Table 4, we are able to mimic the behaviour observed in Francesco Pappalardo et al. (2014).
In Francesco Pappalardo et al. (2014), the authors test a dendritic vaccine in mice without the presence of a tumour. They found that after 154 days and four injections, the number of effector cells decreased, indicating immuno-exhaustion where the vaccine is becoming less effective over time. Despite the overall decrease, the per cent of the effector  cells which are memory cells increases, from about 20 to about 40% from day 63 to 126. After adjusting the parameters (see Table 4), we simulate the vaccine schedule in Francesco Pappalardo et al. (2014) and our model gives similar results as seen in Figure 13. In our simulation, the vaccine is administered into the blood in 4 doses of 5 × 10 5 cells given over a half hour period at days 0, 28, 70, and 112 with the experiment terminating at day 154. There is no tumour present so the initial conditions of all state variables are set to 0. As in Francesco Pappalardo et al. (2014), we see that both memory (E m spleen ) and active effector cells (E a spleen and P) decrease dramatically after the initial dose yet proportionally the memory cells increase over time. While our modelling of immuno-exhaustion here is phenomenological (we assume that the parameters β d and λ m are fixed functions of time), we are able to mimic what is observed experimentally. A next step in the modelling process would be to develop a more sophisticated, physiologically based model of immunoexhaustion.

Conclusions
We have developed a mathematical model of the activation of the adaptive immune response by a dendritic cell vaccine for cancer. By comparing the model to five experimental studies, we have shown that the model is flexible, and that it can capture the immune kinetics in both mice and humans, and for different cancer types.
We have identified several model parameters that most affect tumour size in the initial stages of the disease. These include the intrinsic growth rate of the tumour and the maximum lysis rate of tumour cells by effector cells. Some immunological treatments are under investigation that affect these parameters, and the model can be used to determine threshold values for these treatments. For example, we see in Figure 6 that there are critical values of the parameters d and r that determine whether the vaccine is effective in controlling the tumour and Figure 7 shows the relationship between these parameters and the tumour behaviour. These threshold values can be determined analytically and then used to determine patient specific treatment for a patient-specific set of parameters. The administration of anti-PD1 antibody could increase the parameter d, while a drug such as ibrutinib, which inhibits the survival signalling pathway could decrease r, (Sagiv-Barfi et al., 2015). In future work, we plan to present a stability analysis and a thorough bifurcation analysis of these two parameters that would inform an optimal prescription of this type of adjuvant immunotherapy.
Treatment scheduling is another important consideration in the design of vaccine protocols. Results from clinical trials of DC vaccines have been somewhat disappointing: often, the vaccines delay tumour growth, but do not ultimately control it, as seen in Profile (2006). A mathematical model can be calibrated to a particular type of patient using in vitro assays that determine tumour lysis by effector cells, and by estimating the intrinsic growth rate of the tumour from patient data. This model can then be used to optimize vaccine doses and timings for that particular patient. Both analytical and heuristic optimization techniques have been used in the context of cancer treatments (de Pillis, Fister, & Gu, 2008;de Pillis & Radunskaya, 2001, 2003Ledzewicz, Mosalman, and Schättler, 2013;Swierniak & Duda, 1994;Schättler & Ledzewicz, 2015;Swierniak, Ledzewicz, & Schattler, 2003;Villasana & Ochoa, 2004); in future work we plan to apply these tools to this model of DC vaccines in combination with other immunotherapies.
The success of cancer vaccines will depend on developing combination therapies and vaccination protocols that address immuno-exhaustion and immuno-suppression. In future work, we will test proposed mechanisms that contribute to immuno-exhaustion by incorporating the accumulation of antigen in the lymph organs. Using these more refined models, we will be able to suggest optimal dosing and delivery protocols. Advances in targeted delivery, triggered release via thermo-or sono-sensitive nanoparticles, and sustained release vehicles will enable the implementation of sophisticated schedules of combination therapies. Simulations of these delivery schedules using a validated mathematical model will significantly speed up the time needed to take these new treatment protocols from bench to clinic.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by NSF DMS [1239280].