The Goldilocks Window of Personalized Chemotherapy: Getting the Immune Response Just Right

Immune tolerance occurs when the immune system fails to respond to a target despite its potential to induce an immune response. In cancer, this failure leads to immune evasion and tumor growth. CD8⁺ effector T cells, also known as cytotoxic T lymphocytes (CTL), are an essential component of the adaptive immune system capable of responding to tumor antigens and inducing cell death. Immunologically inert tumors induce T-cell tolerance through multiple direct mechanisms such as inhibition of programmed death ligand 1 (PD-L1), activation of the T-cell regulatory protein CTLA4, and production of regulatory cytokines and metabolites (1), as well as indirect methods such as recruitment of regulatory T cells (Treg), myeloid-derived suppressor cells (MDSC), and tolerogenic dendritic cells (DC; ref. 2) Tregs inhibit CTL cytotoxic activity via cell-cell contact (3, 4) and through secreted factors such as TGFβ (5, 6). They have posed challenges for cancer immunotherapies as well as preventing the activation of the immune system during more traditional therapy approaches (4, 7).

Breaking tolerance requires removal of multiple suppressive factors and activation of cytotoxic immune cells. Chemotherapy, while toxic to CTLs, also has paradoxical and important immunostimulatory effects through dysregulation of the immunosuppressive tumor microenvironment by reducing regulatory cytokine levels, changes in oxygen levels, and reduced metabolites. Several chemotherapies, including cyclophosphamide, paclitaxel, gemcitabine, and 5-fluorouracil, can selectively target MDSCs and Tregs (8, 9).

In addition, highly cytotoxic chemotherapies with lymphodepleting effects create immunologic space (10, 11). During homeostasis, the body maintains T-cell pools at consistent levels. When these pools are depleted, T cells refill this space through antigen-independent proliferation, termed homeostatic repopulation, which favors memory T cells (12). This homeostatic proliferation breaks tolerance, temporarily restoring immune response to previously tolerated antigens (13). This was first characterized in the posttransplant setting where memory T cells lose peripheral tolerance during homeostatic repopulation, leading to graft rejection (12).

Chemotherapy-induced tolerance breaking is dynamic and transient, often requiring treatment breaks to achieve full effect. Various studies report that regulatory cells return 5–10 days post-treatment (8). Homeostatic repopulation following moderate lymphopenia can fully restore the lymphocyte pool as early as 4 days following therapy in murine models (14). Even in the case of nearly complete lymphodepletion using alemtuzumab in nonhuman primate transplant models, the T-cell pool is completely restored in 8 weeks, consisting of 96% memory T cells (15). An obvious question then arises: is there an optimal chemotherapy schedule that could maximize tumor kill and also enhance immune response?

To investigate this question, we created a mathematical model of the complex tumor-immune dynamics that occur during multiple cycles of chemotherapy. In particular, we investigated three, clinically relevant, therapeutic dynamics: immunodepletion, immunostimulation via vaccination, and immunosupportive prophylactics. We identified significant immune trade-offs during chemotherapy as well as the relevant patient metrics that determine the magnitude and severity of these compromises. Furthermore, by exploring the impact of clinically-established therapy, as well as more experimental treatment decisions, we illustrate a more complex interplay between chemotherapy and patient immune dynamics than has been previously investigated. Our results indicate that optimal chemotherapy requires identification of a “Goldilocks Window” in which treatment can both induce cytotoxic effects in the tumor and enhance the immune response to tumor antigens. Therefore, instead of the one-size-fits-all paradigm of fixed therapy regimens, patient immune biology should be a key consideration when developing personalized chemotherapy strategies.

A central assumption of this work is that a clinically detectable tumor has induced a tolerant state in which the immune system can no longer respond to tumor antigens. Chemotherapy temporarily removes this tolerance through lymphodepletion, which eliminates Tregs and allows a burst of immune response. However, the lymphodepletion itself also kills CTLs and therefore reduces the potential cytotoxic efficacy. This double-edged response to chemotherapy implies that there is an optimal therapeutic strategy.

We develop a mathematical model that includes five major populations of cells: tumor cells (T), CTLs (E), Tregs (R), memory T cells (M), and naïve T cells (N). Immune function is separated into two distinct temporal stages, relative to the time of application of each chemotherapy cycle: (i) a period of CTL expansion immediately following the application of chemotherapy (Fig. 1A); and (ii) CTL contraction as tolerance returns (Fig. 1B). The transition time between these phases remains poorly characterized, but empirically occurs 5–10 days after the expansion starts (16). This range has been observed in murine models and is dramatic, involving over a 90% decrease in population size (17). A central assumption of this work is that a clinically detectable tumor has induced a tolerant state in which the immune system can no longer respond to tumor antigens. Systemic lymphodepletion, including that caused by chemotherapy, can help break this tolerance. This can have drastically different effects depending on the type and strength of lymphodepletion (18, 19). First, chemotherapy can selectively reduce Tregs (20, 21, 22) helping to break peripheral tolerance. Second, strong lymphodepletion will cause homeostatic proliferation in the lymphoid compartment, further reducing tolerance.

However, dead immune cells cannot elicit cytotoxic effects or engage in homeostatic proliferation. This implies that there is an optimal therapeutic strategy. If the dose is too high, then the few remaining immune cells will not be able to take advantage of the tolerance breaking; if the dose is too low, then the lymphodepleting effects will be insufficient to break tolerance. In addition to these immune effects, the chemotherapy itself can induce cancer cell death affecting both the tumor size directly and releasing tumor antigens, adding another layer of complexity to the tumor-immune dynamics.

While the full course of lymphocyte recoveries are not observed in the treatment course, measurements of lymphocyte populations over time have shown that a stable equilibrium is reached between chemotherapeutic depletion and population sizes (23). Therefore, in the model, there is a window of 5 days immediately following each chemotherapy cycle in which the immune system is sensitive, and outside of these periods, it is tolerant. We explore the length of this window more thoroughly below.

Our efforts to use mathematical modeling to inform chemotherapy build upon previous immune and personalized medicine works. Mathematical models of tumor-immune activity are numerous, given the complexity of the mechanisms involved (see refs. 24–28) for examples relevant to the current work). Explained more fully below, we extend the modeling work of Robertson-Tessi and colleagues (29) to a more clinically oriented setting by simplifying the immunosuppressive dynamics while maintaining Treg recruitment and function. There have been efforts to study explicit spatial dynamics of the growing cancer cell population in the context of healthy tissue (30, 31). Here, we implement an implicit spatial limitation on growth (see explanation of f(T) below); our model may be extended in the future to incorporate explicit spatial dynamics. Our initial framework choices have been to incorporate patient immune parameters to build toward a model for personalized oncology (32).

During the phase in which the immune system is sensitive to the tumor, a few key processes occur. CTLs, which target and kill the tumor, are recruited from a memory cell population due to response to tumor antigens (16). Recent studies indicate that memory T cells make up the majority of T cells engaged in homeostatic repopulation (15, 33). These memory cells are constantly undergoing a low level of replenishing proliferation, but they only convert to CTLs during the sensitive expansion phase following lymphodepletion. During this phase, there is also tumor-mediated recruitment of Tregs. This eventually causes a significant shift in immune dynamics, leading to a contraction of the CTL compartment during the tolerized phase. Under tolerance, there is no longer a significant recruitment of CTLs from the memory cell compartment. Instead, while the existing CTLs do carry out some tumor-killing function, the Tregs decrease the CTL number.

Our immune tolerance model assumes that the growth of tumor cells (T) can be checked by CTLs (E). However, CTLs are themselves inhibited by Tregs (R) that are recruited at a rate σ by tumor antigens. This leads to CTL-mediated tumor cell death being moderated by the quantity of Tregs |( {\bi \frac{{\bm{R}}}{{{\bi{R}} \bi+ \bi{\bm{E}}}}} )$|⁠. CTLs exhibit different behaviors during immune expansion and immune contraction. This switching behavior is modeled with the Heaviside function |( {{\bi{H}}( {{{\bi{t}}_{{\bi{off}}}} - {\bi{t}}} )} )$|⁠. During the immune expansion phase, CTLs are recruited from the memory pool based on both available memory cells (M) and the tumor burden |( {\frac{{{\bi{TM}}}}{{{\bi{T}} + {\bi{M}}}}} )$|⁠. During immune expansion, the antigenicity of the tumor |( {\bi{\bialpha }} )$| induces differentiation to CTLs |( {\frac{{{\bi{TM}}}}{{{\bi{T}} + {\bi{M}}}}} )$|⁠. However, as immune tolerance sets in, there is a contraction in the CTL population, caused by degradation of CTLs by Tregs |( {\bi{b}} )$|⁠. During immune contraction, CTLs can convert back to memory T cells |({\bi{\biomega E}},{\bi{\ \biomega }}\, \lt \,{\bf 1})$|⁠. Finally, the total remaining lymphocyte population that is not sensitive to the tumor (N) replicates in a logistic growth model |{\bi r _{\bi{N}}}{\bi{N}}( {\bf 1} -{ \frac{{{\bi{M}} + {\bi{N}}}}{{{{\bi{K}}_{{\bi{max}}}}}}} )$|⁠.

Tumor growth dynamics (term 1) are approximated via a combination of exponential growth for smaller tumors and power law growth for larger tumors. This growth model includes a few key assumptions about the limitations that a growing tumor faces before clinical detection. In the absence of effector cells attacking the tumor population, tumor cells first grow exponentially but then transition to power-law growth. This growth dynamic is typical of early-stage, preclinical malignant growths and is based on mathematical modeling as well as experimental observation (29). Furthermore, there are also practical limitations to the biological validity of the tumor population sizes that the model can approximate. While the model can simulate unbounded tumor growth, this is obviously clinically impossible due to the resulting morbidity and eventual patient mortality. Here, we restrict the analysis to the range of tumor sizes that are typical for clinically detectable masses, namely |{\bmi{T}} \,\lt {\bf 10}^{\bf 10}}$| cells. The transition between exponential and power law growth dynamics is governed by f(T) as defined in Eq. B.

The function f(T) employs the method of modeling tumor growth in ref. 29. Beyond a certain size (T_trans), small tumors are not able to sustain their early exponential growth due to physical and nutrient limitations, and therefore transition to power law growth at larger tumor sizes. The smoothness of this transition is governed by the exponent P. The parameter r_T represents the tumor growth coefficient.

Term 2 of Eq. A represents the tumor loss due to killing by CTLs. Parameter k₀ represents the CTL cytotoxic efficacy, with the actual tumor kill rate dependent upon the relative numbers of tumor and CTLs |( {\frac{{{\bi{TE}}}}{{{\bi{T}} + {\bi{E}}}}} )$|⁠. An estimate of this efficacy was initially set at 1 day⁻¹ based on the potency of CTLs in preventing tumor growth when stimulated by multiple types of tumor antigen (34). In vivo killing capacities of CTLs have also been measured in the 1–10 day⁻¹ range by real-time imaging in viral systems, although there is significant heterogeneity (35). However, this rate is mitigated by the presence of Tregs, with b representing their inhibition efficacy. As Tregs increase in density, the CTL-mediated tumor death rate decreases.

CTL dynamics are modeled in two phases, expansion (terms 1–3) and contraction into tolerance (terms 4–6). Terms 1 and 4 switch between these phases via the Heaviside function, with t_off being the length of the expansion phase (5 days, unless noted) immediately following each round of chemotherapy. Terms 2 and 3 govern the growth of CTLs during immune sensitivity to the tumor. CTLs are generated based upon the antigenicity of the tumor (α) as well as the number of tumor and memory cells. Modulating this is an amplification rate, γ, since one memory cell can yield multiple CTLs. Term 2 accounts for the maximum number of lymphocytes that can be supported by the cytokine pool. This paradigm of CTL function being limited by cytokine availability is supported by lymphodepletion studies showing increased CTL activity when IL7 and IL15 cytokine-responsive cells were removed (36). When the immune compartment is full and in homeostasis, this term will be near zero, effectively shutting down CTL recruitment; however, immediately after a dose of chemotherapy, memory and naïve T cells are depleted, which promotes CTL expansion.

Term 5 represents the contraction of the CTL compartment that occurs due to immune tolerance. The death rate of CTLs during contraction, ρ, is due to decreases in the level of supportive cytokines. This rate is increased by the relative fraction of Tregs that are present, |\frac{{\bi{R}}}{{{\bi{R}} + {\bi{E}}}}$|⁠. The modifying constant c represents the sensitivity of CTL suppression to Tregs through a variety of mechanisms (37). Finally, term 6 represents the rate of conversion of CTLs back into memory cells, an active mechanism during immune contraction (38, 39).

Memory cells continually replenish themselves through homeostatic growth in term 1. Parameter r_M is the maximum memory cell growth rate, subject to a carrying capacity, K_max. During the immune expansion phase (terms 2–4), memory cells convert to CTLs, governed by the relative abundances of tumor and memory cells, |\frac{{{\bi{TM}}}}{{{\bi{T}} + {\bi{M}}}}$|⁠, as well as the antigenicity (α). As described in Eq. C, the rate of recruitment is moderated by the homeostatic fraction of the overall immune system (term 3). During the contraction phase (terms 5 and 6), memory cells are replenished from the CTL compartment. A fraction (ω) of the CTLs is successfully converted back to memory cells (38). Because of some loss and inefficiency in conversion, ω <1 (40).

Tregs are recruited by tumor cells with rate σ, and they decay with a rate δ_R (41, 42).

Naïve T-cell dynamics follow homeostatic proliferation with rate r_N, up to a common carrying capacity of K_max, which is the maximum number of memory and naïve T cells in the immune system (43).

The model was parameterized on the basis of literature sources when possible, as shown in Table 1. For many cases, there was evidence of variation in parameters, and some cannot be easily measured. Where possible, we have tried to make a biologically reasonable order-of-magnitude approximation. To address this parameter uncertainty, we explicitly consider the impact of parameter variation on model results.

To establish tolerance in the system and allow transients from initial conditions to dampen before applying therapy, the simulation was initialized with a tumor size of 10⁷ cells. Chemotherapy was started when the tumor reached 10⁸ cells and was simulated as periodic doses of cytotoxic therapy at 14-day intervals. In total, 10 cycles of chemotherapy were applied. At the time of each treatment cycle, all cell populations (immune and tumor) were instantaneously reduced by a fraction C₀ representing the cytotoxic effect of chemotherapy on that population.

This instantaneous death fraction can be understood as lethal dose (LD) values with, for example, C₀ = 0.5 representing LD₅₀. The choice for an instantaneous decrease is simplifying, allowing us to omit pharmacodynamics; however, this approach reflects the general potency of many therapy agents. For example, in vitro studies have shown that cellular uptake and incorporation into RNA for 5-fluorouracil occurs as soon as 3 hours after exposure (44). Uptake levels were directly shown to correlate with cytotoxicity. For doxorubicin, cytotoxicity studies have found that just 1 hour of exposure is enough to induce a 90% decrease in viable, colony-forming cells (45).

Immune cells were reduced by the same fraction (C₀) on each chemotherapy cycle. However, to account for tumor resistance to therapy, the fractional tumor reduction for cycle i (C_i) was linearly reduced with each cycle, such that the cytotoxic fraction on the last cycle was 75% of the initial fraction C₀. Approximating the impact of chemoresistance on drug efficacy is challenging since values vary for different classes of drugs. Furthermore, Hao and colleagues (46) noted dose-dependent differences of up to 400% between resistant and resensitized prostate cancer cell populations to docetaxel. Here, the value of 75% chemotherapy efficacy at the last cycle represents a 33% advantage of survivorship for a resistant population versus a susceptible population. It is a conservative estimate of the impact of resistance, but we believe it is reasonable given that tumor populations are unlikely to be entirely resistant. Varying this range is a relevant question for future research. For our purposes, C_i is given by:

The final tumor size after 10 cycles of chemotherapy was compared with the tumor size at the start of treatment (10⁸ cells) and evaluated according to RECIST categories. Specifically, a total loss of tumor (≤99% change in size) is a complete response (CR). A change between −30% and −99% is considered a partial response (PR). Tumor changes between −30% and +20% are classified as stable disease (SD) and increases of greater than +20% are seen as progressive disease (PD; ref. 47). While there are many different methods of measuring therapy efficacy impact on disease, RECIST categories were chosen here because they have correlated well with overall survival in patients across a variety of cancers.

The model was programmed in the Python language (ver. 2.7.11). The open-source packages Scipy (ver. 0.17.0), Numpy (ver. 1.10.4), and Matplotlib (ver. 1.5.1) were used for simulation of the ODEs as well as visualization of the results. The platform for the program was both an Intel Core i7-6820 HQ processor as well as the high performance computing cluster at H. Lee Moffitt Cancer Center and Research Institute (Tampa, FL). The source code is available at the github repository for the Integrated Mathematical Oncology department at github.com/MathOnco/Goldilocks.

Memory cell population sizes are variable among patients (48). Arstila and colleagues have estimated 10⁶–10⁷ T-cell clones in the human body with approximately 10⁵ memory T cells per antigen (48, 49). However, due to antigen responses being polyclonal, this suggests multiple orders of magnitude of potential variation in memory T-cell numbers. Therefore, varying doses of chemotherapy were simulated for a range of memory cell population sizes (Fig. 2A and B). Results from the model show that patient memory cell numbers significantly influence the optimum chemotherapy dose. Generally, there is a minimum memory cell population size that is necessary for any given strength of chemotherapy to be successful. Above this threshold, the more memory cells there are, the better the improvement with stronger doses of therapy. Conversely, this means that when memory cell populations are close to the minimum threshold, chemotherapy should be similarly weak for a more favorable outcome. Furthermore, if memory cells are below the minimum threshold, then the optimal strategy is to use strong chemotherapy (Fig. 2A and B), because the immune system will not contribute to tumor regression.

The double-edged nature of chemotherapy on the immune system can be better understood through the transient dynamics during therapy (Fig. 2C and D). In cases with stronger chemotherapy dosing, there is an early decrease in tumor population levels due to the cytotoxic strength of the therapy. However, we observe a trend that these strong therapies tend to lead to failure and larger final tumor sizes than if treated with a “weaker” chemotherapy regimen, which provides lower cytotoxicity on the tumor but maintain tumor size reduction for the duration of therapy.

This counterintuitive result stems from the fact that cytotoxicity alone is insufficient for suppressing tumor growth, especially due to the accumulating chemoresistance. Rather, it is the synergistic effect of cytotoxicity as well as the breaking of immune tolerance and consequent recruitment of CTLs that keeps tumor populations in check. Our in silico treatments consistently show that there is an inherent disadvantage to high-dose chemotherapy. There is a gradual decrease in the CTL population over multiple rounds of treatment due to the net loss that stronger dosing causes in memory T-cell populations (Fig. 2D). It is these memory cells that are affected the most by chemotherapy because they can only recover relatively slowly. If the cytotoxic pressure on memory cells is greater than the recovery rate of that compartment, then even with a resensitized immune system, expansion will lead to fewer total CTLs and ultimate treatment failure. In contrast, if the immunodepleting side effects of chemotherapy can be balanced with immune recovery, then more sustainable treatment responses are possible. In short, there is a tradeoff between having chemotherapy strong enough to sufficiently break tolerance, but mild enough to leave sufficient memory T cells for adequate CTL expansion. Akin to the story of Goldilocks and the three bears, the balancing of these two immunologic goals leads to an intermediary chemotherapy strength that is “just right.” In silico simulation shows that this “Goldilocks Window” is highly dependent upon patient-specific, preexisting memory T-cell populations.

While we identified this Goldilocks Window of optimal, submaximal chemotherapy dosing, we also sought to explore it in the broader context of patient-specific disease and immune variation. For tumor growth rates, we found that successful treatment outcomes are more sensitive to chemotherapy dosing for faster growing tumors and less sensitive for slower growing tumors. Experimental and model analyses have shown that selection pressures on growing tumors can lead to significant heterogeneity in metabolism and growth rates (49) In our framework, the tumor growth rate parameter (r_T) was set to the maximum speed for doubling during the exponential growth phase (1,000 cells^0.5 day⁻¹, representing a doubling time of 1 day), but we also explored faster and slower growth rates (Fig. 3A and B).

In slower growing tumors (r_T < 1,000 cells^0.5 day⁻¹), chemotherapy's cytotoxic effects are sufficient for tumor control. After the partial tumor clearance due to each cycle, there is regrowth of the cancer cell population (Fig. 2A and B). For slower growing tumors, there is less intercycle regrowth and therefore cancer cell populations can be controlled by chemotherapy alone without the need for CTL killing. The result is that, for slower growing tumors, there is no Goldilocks Window.

However, for faster growing tumors (r_T > 1,000 cells^0.5 day⁻¹), it becomes necessary to maintain chemotherapeutic strength within the Goldilocks Window to achieve optimal outcomes. For these tumors, regrowth between chemotherapy doses is significant and demands the addition of CTL-mediated tumor killing for disease control. Chemotherapy that is stronger than the Goldilocks Window hamstrings the patient's immune activation.

Importantly, for the most aggressively growing tumors, there is actually a “worst case scenario” of intermediary chemotherapy strength (Fig. 3B). Here, the worst option for chemotherapy is not the strongest possible dose, but a “mid-range” strength in treatment instead. At this chemotherapeutic strength, the drug alone is insufficient to cause a reduction in tumor size. However, the dose is still strong enough to lead to severe memory cell depletion, undermining any immune efforts at constraining tumor growth.

Separate from tumor parameters, patient immune characteristics can also impact the sensitivity of treatment outcomes to chemotherapy dosing. One important parameter we sought to explore was the rate of CTLs in killing tumor cells (k₀, Fig. 3C and D). Without changing the initial patient memory cell populations (M₀ = 10⁶ cells), or the tumor growth rate (r_T = 1,000 cells^0.5 day⁻¹), the CTL-mediated cytotoxicity rate was varied around the biologically realistic parameter of k₀ = 0.9 day⁻¹ (34). CTL efficacy was found to dramatically impact sensitivity chemotherapy dosing and the Goldilocks Window. Lower rates of CTL-mediated tumor cell death lead to greater sensitivity of treatment outcomes on chemotherapy dosing (Fig. 2C and D). With a lower value of k₀, more CTLs are necessary to exert the same degree of immune control over the tumor. However, strong chemotherapy on a patient with lower k₀ values prevents sufficient CTL expansion by rapidly diminishing the memory cell populations. Higher CTL killing rates, although, removed the restriction of this Goldilocks Window and made successful treatment outcomes less sensitive to dosing. While higher chemotherapy doses may lead to larger immune depletion, more efficient CTLs mean that these smaller immune populations still lead to successful treatment outcomes.

In addition, we examined the impact of changing the window duration for immune expansion immediately following each chemotherapy dose. Current literature indicates that immune contraction can begin to occur anywhere from 4 to 8 days after treatmen (8, 14, 17) When these extremes were explored (see Supplementary Fig. S1), there was no significant qualitative difference to our observation of a submaximal optimal dosing range when compared between 4 days (Supplementary Fig. S1A) and 8 days (Supplementary Fig. S1B). While a longer window of immune expansion (Supplementary Fig. S1B) leads to more favorable outcomes for more rapidly growing tumors when treated in the optimal dosing range, the actual presence of this submaximal dosing range does not change. Furthermore, there is almost no difference in the outcomes of patients who are overtreated. This also implicitly addresses our mathematical implementation of a switch via a Heaviside Function. Specifically, while there might be any number of less abrupt and more gradual transitions between immune expansion and immune contraction, our exploration of the dynamics at the extremes of this transition range would give an idea of what the intermediate behaviors due to a smoother transition might cause. That is, our qualitative results would not significantly change with a smoother function.

In a broader exploration of the model's immune parameters, a general trend was observed that a more robust immune response would improve the outcome (Supplementary Figs. S2–S7). For certain model parameters that were more difficult to accurately estimate from the literature, we explored their variation for the default tumor growth rate and a chemotherapy strength of 25%. If the patient had a stronger immune system characterized by lower CTL death rate (δ_E), lower sensitivity to Tregs (c), greater memory cell expansion (γ), regrowth (r_M), and back conversion (ω), the final tumor population was smaller. Furthermore, more robust antitumor immune responses led to greater maximum possible reductions over the range of chemotherapy. In addition, these changes led to expansions of the Goldilocks Window in terms of chemotherapy doses that could achieve tumor reductions.

In short, patient-specific disease and immune biology determines the sensitivity of treatment outcomes to chemotherapy dosing. For rapidly growing tumors, chemotherapy must be maintained in a submaximal Goldilocks Window to optimize drug and immune synergies. However, patient immune biology matters as well, with weaker immune characteristics also leading to a greater necessity to stay within the Goldilocks Window. Importantly, this presents potentially counterintuitive guidance because an initial motivation may suggest that, in a situation where a patient has a weaker immune system, chemotherapy strength should be increased to compensate. However, our model suggests that the lymphodepleting impact of heavy chemotherapy on an already weakened immune system will only worsen outcomes. When confronted with weaker patient immune systems, chemotherapy needs to be maintained within the Goldilocks Window for successful outcomes.

Patient-specific vaccines have become a recent hallmark in personalized cancer therapy. One of the first to acquire FDA approval was Sipuleucel-T, for treating metastatic castrate resistant prostate cancer (50). Each vaccine is tailored to a specific patient by culturing dendritic cells from patients using a specific tumor antigen. Reinjection into the patient would potentially stimulate a T-cell–mediated antitumor immune response. Three doses were administered in 2-week intervals with significant clinical responses being observed. Vaccination led to a 22% reduction in the relative risk of death, although there was no noticeable decrease in the rate of progression of disease (50). The specific effect on T cells was quantified by looking at T-cell receptor changes in response to vaccination. Certain receptor sequences were enriched, while others were significantly decreased (51), suggesting that the vaccine promoted an antigen-specific immune response against the tumor.

To study the effects and potential synergy of chemotherapy with this method of T-cell stimulation, we simulated a vaccine regime similar to that used for Sipuleucel-T (3 doses, spaced 14 days apart), with different vaccine strengths. Mathematically, this was modeled by modifying the ODEs that govern CTL expansion, without explicitly representing the complex DC-to-T-cell cascade that the vaccine induces. Other models have examined the DC cascade in more detail. For example, the explicit migration of dendritic cells between blood, spleen, and tumor have been modeled via delay-differential equations to better characterize the specific dose timing–dependent responses to therapy (28) For simplicity, here we focus solely on the net effect of the vaccine on T-cell numbers by changing the antigenicity parameter of the tumor, α, from a constant coefficient to a variable, time-dependent function, α_v(t):

Total antigenicity is modeled as the result of both the constant, baseline antigenicity of the tumor, α, and an exponentially decaying vaccine-augmented component, v, which decays with a half-life, t_half = 3 days, a biologically realistic timespan (52). This model of dynamic antigenicity can be expanded for multiple vaccinations, as used in the clinical protocol (Eq. I).

Here, H(t) is again the Heaviside function. The constant n_vac represents the total number of vaccine injections and t_n represents the time of the n^th vaccination. The ODEs used for the simulation of immune and tumor cell populations are then dependent on the instantaneous current value of α_v(t) throughout the course of simulated therapy.

Here, we explored a range of antigenic increases due to potential patient-to-patient variation in responses to immunostimulatory vaccines. While dendritic cell vaccines like Sipuleucel-T administer all of the available dendritic cells, responses in individual patients vary in how much the antigenicity is changed. In our range of exploration, though, we found some commonalities in vaccine interaction for chemotherapy.

Results show that vaccine therapy can improve outcomes, but only within a specific range of chemotherapy strengths (Fig. 4). Treatment outcomes improve when a vaccine is used with moderate chemotherapy (Fig. 4A), but for very high chemotherapy doses, the beneficial effects of a vaccine are diminished. As before, the underlying cause for decreasing efficacy is the persistent lymphodepletion due to the chemotherapy. Antigenicity augmentation due to vaccine stimulation is offset by reduced CTL expansion. However, very low-dose chemotherapy poses its own challenges, because with insufficient lymphodepletion, tolerogenic mechanisms and greater Treg recruitment inhibit any CTL response augmented by the vaccine. The immune system remains closer to tumor-tolerized homeostasis, and as a result, vaccine stimulation is mitigated because the immune system is already suppressed. Therefore, the width of the optimal window is not significantly affected by the vaccine because the vaccine has no effect on a highly depleted or tolerized immune system.

Therefore, even with immunostimulatory vaccines, there exists an optimal “Goldilocks” Window. Quantitatively, we define this window to be the region in which a therapy dose can offer at least a 20% reduction in tumor size because this is the necessary amount for disease to become classified as a partial response. In order for there to be this maximized benefit from vaccine application, the chemotherapy regimen must be “just right.” Chemotherapy must have sufficient lymphodepletion to resensitize the immune system, but must leave enough immune cells such that vaccine stimulation leads to a large CTL response. Similar to the results of chemotherapy without the vaccine, the specific range of this Goldilocks Window depends upon the initial patient memory cell (M₀) numbers.

We note that the small oscillations observed in the plots (Fig. 4B) are a result of the use of dual growth laws for the tumor. Essentially, giving the vaccine causes the tumor to dip into the faster exponential growth phase at an earlier chemo cycle than when chemo is given alone. Because the chemo cycles are discrete and instantaneous, this generates an effective step function to the response with increasing chemo dose, superimposed on the single-peaked optimal curve; this step function is further rounded by both the smoothing exponent P between the growth laws and the nonlinear interactions between tumor growth and immune response at small tumor sizes.

Chemotherapeutic lymphodepletion in the clinical setting can pose a serious threat to the safety of the patient through neutropenia (53), which commonly leads to dose reductions and disruptions to the standard schedule of therapy for patients. Consequently, multiple tools have been developed to help mitigate the effects of chemotherapy on the immune system. For example, it was recognized that dexamethasone treatment before carboplatin and gemcitabine could not only increase chemotherapy efficacy, but also reduce the lymphodepleting effects by preventing uptake in the spleen and bone marrow (54). In contrast, other aspects of cancer therapy can potentially hamper CTL responses to tumor insults. For example, G-CSF application has been shown to reduce CTL activation and could conceivably impede the impact of lymphodepletion as a break from immune tolerance (55, 56). More generally, however, the broader impact of immune system augmentation or suppression during therapy remains unexamined.

To examine the effect of attenuated or augmented lymphodepletion on therapy outcome, we allowed for variable chemotherapeutic toxicity to immune populations, as compared with the tumor population. Mathematically, this simply means modifying the chemotherapy dose by a scaling factor h. The effect of chemotherapy on immune cell populations at a given treatment time is:

where I₁ is the immunologic population size after application of chemotherapy, I₀ is the population size before therapy, and 0 < C < 1 is the dose strength. The specific numerical range in which h falls represents either attenuated or augmented chemotherapeutic toxicity. For values of 0 < h < 1, this represents an attenuated toxicity on the immune system relative to the toxicity on the tumor. In contrast, values of h > 1 represent higher toxicity on patient immune populations than on the tumor. This could be due to patient-dependent increased sensitivity to chemotherapy. The maximum possible reduction of cells by chemotherapy when modified by immune support is 100%. This leads to the resulting condition that hC < 1. For our in silico therapies, h was varied across the allowable ranges for three different strengths of chemotherapy. Values of C were chosen to represent lower (C = 0.25), middle (C = 0.6), and higher (C = 0.9) dose chemotherapy (Fig. 5A).

Interestingly, the results suggest that immune-supporting combination therapy has essentially no benefit when given with low-dose chemotherapy. As shown in Fig. 5, similar tumor reduction occurred for a wide range of values of h around h = 1. Furthermore, outcomes were worse when h was very low or very high. In situations where it was very low, final tumor sizes were large because a lack of lymphodepletion did not sufficiently break immune tolerance. In contrast, for larger h values, there was overdepletion that prevented an effective immune response despite significant tolerance breaking.

In contrast, high-dose chemotherapy saw treatment failure or success highly dependent upon the amount of immune support. Similar to low-dose therapy, a small value of h that mitigated the depleting effects of chemotherapy led to the best possible outcomes in terms of tumor shrinkage. Final tumor sizes were, in fact, multiple orders of magnitude lower than was possible with low-dose chemotherapy. As h increased (representing less toxicity mitigation) treatment outcomes rapidly worsened. The transition value h*, where the clinical outcome rapidly shifts, indicates a threshold effect with regard to immune support. For high chemotherapy doses, immune support treatments must have a significantly large mitigation (h < h*) of immunodepletion in order for successful treatment responses to occur. The position of this inflection point is influenced by the strength of the patient immune system to begin with. In expanded parameter analyses, the strength or weakness of the simulated patient's immune system led to changes in the upper bound of the Goldilocks Window (Supplementary Figs. S2–S7).

The moderate strength chemotherapy regimen yielded only partial benefits of either extreme. The greatest tumor reduction possible, with immune support, yielded tumors that were smaller than those achievable with low-dose chemotherapy. However, these tumors were still multiple orders of magnitude larger than those achievable with high-dose chemotherapy. For treatment failure at lower immune support (h > h*), tumor sizes were actually larger than when high dose chemotherapy failed.

Clinically, the results suggest that chemotherapy dose strength can be used to mitigate uncertainty regarding the amount of immune support a certain treatment will give to a specific patient. Low-dose therapy offers a wide range of potential immune support in which treatment can successfully reduce tumor sizes. The disadvantage is that the maximum tumor size reduction still leaves larger tumors than are possible using higher doses of chemotherapy. While our model has not analyzed this, a potential impact is that larger tumor sizes could lead to more heterogeneous populations and thus lead to a higher likelihood of resistant or metastatic populations. However, higher doses have a narrower range of immune support in which they are successful. Chemotherapy can be balanced, then, against how certain the clinician is of the benefit that G-CSF (or other immune supporting drug) will give. For patients where there is high certainty of a significant benefit due to the drug, high-dose therapy is optimal. In contrast, lower dosing should be used when the drug may have lower or variable efficacy.

Finally, we sought to investigate how variation in the effectiveness of these immune adjuvants might impact treatment outcomes in a group of patients. Chemotherapy treatment leads to a wide range of responses, both successful and unsuccessful, across multiple types of cancer (47). This variation has been attributed to disease variation, patient variation, and interactions between the two. However, less attention has been given to variable patient responses to secondary drugs, such as G-CSF, and how they impact therapy. Patient responses to these secondary drugs are currently poorly measured and could have significant implications for therapy outcomes.

To better explore the effect of variable patient responses to immune support drugs, cohorts of 500 patients were randomly generated from a normal distribution with a mean immune support response value of h = 0.8 and variance of 0.2. These values were chosen to center the distribution around the model-derived threshold value h* = 0.8. While not directly describing patient responses to immune support drugs, a normal distribution for selection was chosen due to the fact that immune cell counts have been found to be normally distributed in population cohorts (57).

Similar to our previous investigations, cohorts were then subjected to regimens of low (C = 0.4) and high (C = 0.8) chemotherapy strengths (Fig. 5B). Percent changes in tumor size after therapy were displayed for each individual patient in the cohort to generate a waterfall plot. In doing so, we used our model to simulate cohort responses as is commonly measured in aggregated studies of patient data (47). The waterfall plots (Fig. 5) illustrate that chemotherapy strength can significantly change the proportion of successfully responding patients in a population with variable responses to immune prophylactics. This is significant because the proportion of successful responses is often an important criterion for judging therapeutic efficacy. The simulated waterfall plots show how clinical outcomes could not only be the result of therapy, but also due to inherent immune variation within the cohort.

A major barrier to success for immunotherapy in cancer is tolerogenic mechanisms that reduce the immune response to tumor antigen (4, 7, 58). A potential solution has come from observations that lymphodepletion stimulates homeostatic proliferation in the immune system that can transiently restore an immune response. This has led to increasing efforts to selectively apply chemotherapy to improve outcomes from immunotherapy (59).

To better understand this potential synergy, we constructed a mathematical model to frame these complex dynamics and identify critical parameters that govern the clinical outcomes. Our studies focused on three clinically observed dynamics of immunodepletion, immunostimulatory vaccination, and immunosupportive prophylactics. With regard to immunodepletion, we demonstrated that chemotherapy results in a trade-off. At very high doses, chemotherapy has a maximal cytotoxic effect on the tumor but also maximally depletes T cells such that no effective CTL response can be mounted despite the transient loss of tolerance during reexpansion of the immune cells after completion of chemotherapy. Similarly, low doses of chemotherapy are insufficient to produce the post-treatment immune cell expansion that is necessary for reversal of immune tolerance.

Importantly, however, we find there is a Goldilocks Window of chemotherapy doses in which lymphodepletion causes adequate immune resensitization, but does not impose an overly large recovery burden. This window is governed by the patient-specific quantity of memory T cells so that larger pretreatment T-cell populations allow more favorable outcomes with higher doses of chemotherapy. In contrast, fewer pretreatment CTLs can limit the immune response even in the Goldilocks window of chemotherapy. Thus, there is a necessary “minimum efficacy” of CTLs for successful stimulation of immune response by chemotherapy. Below this threshold of immune activity, the benefit of chemotherapy is almost solely dependent on its inherent cytotoxicity (Fig. 6).

Our model also provides insight into the potential effects of variation in the tumor growth rate. In slower growing tumors, chemotherapy alone can be sufficient to achieve optimal treatment response. Treatment of faster growing tumors, however, is best when the chemotherapy is administered to enhance the immune response. Unfortunately, if the pretreatment population of CTLs is small, we find chemotherapy for rapidly growing tumors will be ineffective if it is both highly lymphodepleting and insufficiently cytotoxic to significantly reduce tumor growth. Assessing the clinical importance of this question is challenging because it remains unclear from the literature as to the actual size of the population of tumor-specific T cells that are present during treatment. In spite of these difficulties, the impact and existence of antitumor immunity has been bolstered by recent immunotherapies that act to remove inhibitions to T-cell action (60).

Chemotherapy is increasingly being used in concert with vaccines to help stimulate the patient immune system. We investigated the interactions between vaccines and lymphodepletion and found that, as before, there is a window of chemotherapy ranges in which vaccines can improve outcomes versus chemotherapy alone. At very high doses, however, the resulting lymphodepletion substantially reduces benefits of immune stimulation by vaccination. More broadly, other novel immunotherapies could also potentially be hampered by overdepletion of the immune system.

To further investigate the potential impact of this interaction, we modeled the effect of differential responses to immune prophylactics. G-CSF and other drugs have become common recourses in chemotherapy for mitigating the immunodepletion effects on patients (61). However, recent studies have suggested that T cell response is hampered by G-CSF administration (55). While G-CSF may help prevent neutropenia and cytopenia for patients, it may impede the ability of retolerized T cells to mount an antitumor response. In addition, responses to prophylactics are not constant but the significance of this variation remains relatively uninvestigated. Our model suggests that interpatient variation in prophylactic response can lead to drastically different outcomes for the same dosing of chemotherapy. Across larger samples, this variation can further interact with chemotherapy to be a significant determinant of whether the chemotherapy dose leads to more success or failure across a range of patients.

In the clinical literature, our model results cautioning about balancing chemotherapy and immunogenic effects has been echoed in multiple situations. Previous studies have explored the mechanisms of action in mAb-based treatments including targeting of HER2 (62, 63). When quantifying the impact of antibody-dependent cytotoxicity-mediated by CTLs, it was noted that addition of paclitaxel reduced the lasting impact of the immune response generated against the tumor. While in the short term, higher doses of chemotherapeutic agents could induce larger tumor reductions, mice that were given both antigen and large chemotherapy doses were more susceptible to tumor rechallenge. Similarly, in radiotherapy it has been found that CTL priming occurs due to antigen-dependent cell death (64). However, the addition of even a small amount of paclitaxel was found to induce a significant reduction in CTL numbers. Adjuvant chemotherapy regimens were found to significantly abrogate the immunogenic benefits of radiotherapy-induced immune responses, while immunotherapies increased the efficacies. This result is also significant because it implicitly addresses whether our results might hold when antigen increase, due to cell death, is accounted for. In this mouse model, even with tumor-cell death–mediated antibodies, the loss of T cells leads to a worse overall outcome (64). This presents a natural extension of our framework to be applied to a specific disease and chemotherapy dosing setting. While we created a general model of chemotherapy, there may be interesting dynamics unique to individual cancers that could be explored. It would also allow the employment of more complex pharmacodynamics modeling for specific treatment regimens.

In conclusion, our results suggest opportunities to increase the efficacy of immunotherapy with precise application of chemotherapy. Our model affirms the importance of CTL and memory T-cell expansion following chemotherapy to reduce immune tolerance to tumor antigens. However, we demonstrate that optimal chemotherapy requires identification of a Goldilocks Window in which treatment can both induce cytotoxic effects in the tumor and enhance the immune response to tumor antigens. Identifying optimal strategies for chemotherapy in each patient will likely benefit from the application of mathematical models that are parameterized by patient data pretreatment to generate an optimal treatment strategy for that patient. Importantly, these predicted strategies would most likely need to change as patient responses diverge from those predicted, leading to an iterative loop of “predict-apply-refine.” With the growing drive toward precision medicine, we believe that mathematical models are critical for the future of truly personalized therapy, where no two patients will receive the same therapeutic regimen, and where treatments adapt a change based on patient responses. The model presented here is a step toward describing the complex landscape of treatment decisions regarding dosing and combination of different therapies, and we have shown how these decisions can be sensitive to patient-specific parameters and guide clinical intuition.

No potential conflicts of interest were disclosed.

Conception and design: D.S. Park, M. Robertson-Tessi, K.A. Luddy, P.K. Maini, A.R.A. Anderson

Development of methodology: D.S. Park, M. Robertson-Tessi, K.A. Luddy, P.K. Maini, A.R.A. Anderson

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D.S. Park

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D.S. Park, M. Robertson-Tessi, K.A. Luddy, P.K. Maini, M.B. Bonsall

Writing, review, and/or revision of the manuscript: D.S. Park, M. Robertson-Tessi, K.A. Luddy, P.K. Maini, M.B. Bonsall, R.A. Gatenby, A.R.A. Anderson

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D.S. Park

Study supervision: M. Robertson-Tessi, P.K. Maini, A.R.A. Anderson

Other (source code maintenance): D.S. Park

D. Park was supported by a 2014 Marshall Scholarship from the Marshall Aid Commemoration Commission of Great Britain. A.R.A. Anderson, R.A. Gatenby, M. Robertson-Tessi, and K. A. Luddy were supported by a Physical Sciences-Oncology Center grant from the National Cancer Institute of the United States of America (grant number 1U54CA193489) and the Center of Excellence for Evolutionary Therapy at H. Lee Moffitt Cancer Center and Research Institute.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.