Optimal dosage protocols for mathematical models of synergy of chemo- and immunotherapy

The release of tumor antigens during traditional cancer treatments such as radio- or chemotherapy leads to a stimulation of the immune response which provides synergistic effects these treatments have when combined with immunotherapies. A low-dimensional mathematical model is formulated which, depending on the values of its parameters, encompasses the 3 E’s (elimination, equilibrium, escape) of tumor immune system interactions. For the escape situation, optimal control problems are formulated which aim to revert the process to the equilibrium scenario. Some numerical results are included.

1 Introduction

The release of tumor antigen during traditional cancer treatments, such as radio- (1) or chemotherapy (2), can lead to a stimulation of the innate immune response which, in some cases, is able to trigger protective antitumor immunity with possibly long-lasting effects (3, 4). For example, a well-documented abscopal effect of radiation therapy (i.e., the reduction of tumor metastases in areas well outside the field of radiation) is hypothesized of being immune mediated (5–8). These stimulatory effects are the basis for an observed synergy some of these therapies have with immunotherapy, e.g., check-point blockades.

Mathematical models of tumor growth and treatment have a long history going back to the 1980s with research unabatedly continuing (e.g., see (9)). Probably the earliest works on tumor-immune interactions are Stepanova’s paper (10) and (11) by Kuznetsov et al. while mathematical models including immunotherapies are more recent [e.g (12–20)]. Capturing all aspects of tumor-immune interactions in a mathematical model is difficult as the competitive interactions between tumor cells and the immune system are complex, to say the least, and still are the topic of intense medical research. While large-scale, agent-based, PDE, or hybrid models are more precise, they suffer from the inabilities to determine a large number of parameters.

In this paper, we consider a qualitative, low-dimensional mathematical model (21). Rather than analyzing the dynamical system for a specific set of parameters, the aim of a qualitative analysis is to understand the totality of all the possibilities the model allows for. Especially for nonlinear models this is an important aspect in pointing out both mathematical limitations and complexities of the underlying dynamics. Motivated by the papers (10, 11, 22, 23) here we consider a model for tumor immune system interactions which in addition to tumor volume and immuno-competent cell densities includes as a third separate variable tumor antigen. This model retains the main aspects of tumor immune system interactions called the three E’s of cancer immunoediting (24, 25). These are (i) elimination: the immune system is able to completely eradicate the tumor; (ii) equilibrium: adaptive immunity is able to maintain cancer in a benign state (tumor dormancy) (26), and (iii) tumor escape: tumor growth overcomes or evades the actions of the immune system (27). From a practical (therapeutical) point of view, the first case will never be seen while therapy will not be able to save the patient in the last. Only when it is possible to influence the dynamics (that is, tumor growth) through treatment which will lead to positive changes on a permanent basis, i.e., even after treatment has been stopped, a cure is possible. This requires that the system can be reverted back to the equilibrium situation.

Optimal control problems are optimization problems in time in which the full range of possibilities to influence the dynamics of a system, i.e., without a priori restrictions on these structures, is considered. This is very different from a “best in class” argument sometimes wrongly also called optimal control in many publications where only a limited number of possibilities is considered, often by mere simulations. Typically in optimal control problems the aim is to transfer the state of a dynamical system from a given initial point into a desired set of terminal states. This is accomplished by minimizing some suitably chosen criterion subject to the dynamics of the system and other constraints that need to be satisfied. Solutions then are functions of time which describe the best actions relative to the chosen criterion. Historically, there has been great success in applying optimal control to engineering problems (moon landing, autopilots on airplanes) and economics (portfolio optimization) while medical applications with its uncertainties in the dynamics—these are generally based on ad-hoc modeling premises, not like in physics on first principles—and usually a large parameter uncertainty in the model lack similar success stories. Nevertheless, the scheduling of therapeutic agents over time has all the characteristics of an optimal control problem (28) and there is an increasing bulk of literature in which optimal control is applied to medical problems, e.g., see (19, 29–33). The aim is to minimize some objective related to tumor burden and quality of life of the patient while the underlying system follows the processes of tumor development and treatment interactions. While direct clinical applications are a mere possibility for the future in our opinion, currently the use of optimal control techniques lies more in understanding the dynamics of mathematical models proposed to study medical processes (which contributes to a validation of these models) while solutions to optimal control problems can be helpful in identifying realistic therapy protocols that possibly can be explored in medical trials and practice (28, 31). The contributions of our paper lie in this direction. We formulate an optimal control problem for a qualitative mathematical model of tumor immune system interactions which considers the transfer of the state of the system from a malignant initial condition (corresponding to a tumor escape situation) into a benign state (corresponding to the equilibrium scenario) and show how geometric properties of the dynamics help in formulating and understanding the proper goal of treatment. We discuss the complexities of obtaining optimal controls in this case and highlight some numerical results.

2 Methods

A low-dimensional, qualitative model for tumor immune system interactions under chemo- and immunotherapy is formulated as a dynamical system and analyzed mathematically. For the medically relevant scenario of tumor escape (malignant), optimal control problems are formulated whose solutions would revert the system to the equilibrium case (benign).

2.1 Mathematical model

We consider the following dynamics (21) for tumor immune system interactions based on classical papers by Stepanova (10) and Kuznetsov et al. (11):

$\begin{array}{ll}\dot{x}=\xi x\left(1-\frac{x}{x_{\infty }}\right)-\theta xy-\alpha xu,& (1)\end{array}$

$\begin{array}{ll}\dot{y}=a(1-bx)yz+\gamma -\delta y-\kappa yu+\nu yv,& (2)\end{array}$

$\begin{array}{ll}\dot{z}=\sigma x+\psi xu-\mu z,& (3)\end{array}$

State variables are the tumor volume x, the immunocompetent cell density y, and tumor antigen z. The variable y is a non-dimensional, order of magnitude quantity which is related to various types of T-cells activated during the immune reaction and summarily represents the actions of the immune system. The variables u and v represent time-varying dose rates u = u(t) at which chemotherapy is given and a time-varying immune boost v = v(t). For simplicity, drug dose rates and concentrations are identified. (It is well-known how to deal with the required changes if standard pharmacokinetic models are included (34). All Greek letters and a and b are parameters which for the time under consideration are assumed constant. The meaning of variables and parameters is given in Table 1.

Table 1

Table 1 Variables and parameters.

Most of the terms in the equations are standard. Log-linear terms of the Skipper model (35) are used to formulate the damage done to the tumor through the concentrations of the agents and a logistic growth model is used for the tumor volume. This is merely for sake of specificity and analogous results hold qualitatively, for example, for a Gompertzian growth function. Equations 1, 2 follow the classical papers by Stepanova (10) and Kuznetsov (11) and have been taken over with only small changes. Equation 3 extends these earlier 2-dimensional models to include a direct link between tumor antigen z and the immuno-competent cell density y. This has led to the modified term a(1 − bx)yz used in Equation 2. In (10) this term instead was taken of the form a(1 − bx)y² with the justification that tumor antigen would be proportional to the tumor volume thus generating the factor y² in the interaction term as tumor antigen is not considered separately in that model. In (2) we have therefore replaced one the factors y with z restoring a direct link between these two variables. Equation 3 is based on a similar equation in (36) and models the evolution of tumor antigen. It is assumed that the tumor produces antigen z at rate σ which results in an intrinsic (i.e., not therapy induced) immunogenicity of the tumor. Antigen is cleared by the immune system at rate µ which leads to the creation of immune effector cells which generate a stimulating effect onto the proliferation of lymphocytes and thus a positive influx into the compartment determining the immunocompetent cell density y. This effect is represented by the term ayz in Equation 2. The term ψxu models the immuno-stimulatory aspect of therapy assuming that the tumor produces antigen at a dose dependent rate ψu with ψ modeling the therapy induced immunogenicity of the tumor.

2.2 The three E’s of immuno-editing

Depending on the value of the parameters, the dynamical Equations 1–3 properly replicate the full variety of medically realistic scenarios. In order not to be confusing with the medical notion of equilibrium, we use the terminology stationary point for the states which are obtained as solutions when the derivatives in Equations 1–3 are set to zero. There always exists a tumor-free stationary point given by $w_0=(0,\frac{\gamma }{\delta },0)$. It is stable if ξδ < θγ and unstable if ξδ > θγ. Intuitively, stability means that solutions of the dynamics which start near w₀ converge to w₀ in time. The relevant term is a difference between products of tumor stimulating parameters (the tumor growth rate ξ and natural death rate δ of immune cells) and tumor inhibiting parameters (the influx γ stimulating the immune system and the effectiveness θ of the immune system fighting the tumor). Stationary points with positive tumor volumes x_∗ are zeros of a cubic polynomial Q = Q(x) computed by eliminating y_∗ and z_∗ from the equations $\dot{y}$ = 0 and ż = 0. Given the logistic growth model used in Equation 1, only zeros in the range 0< x_∗ < x_∞ are viable solutions for the tumor volume.

The three E’s of immuno-editing correspond to the following scenarios:

2.2.1 Elimination

This situation arises if the tumor-free stationary point is stable and no stationary points with positive tumor volumes exist. All solutions of the dynamics converge to the tumor-free stationary point, i.e., the actions of the immune system are eliminating the tumor. While this is not a relevant scenario medically—in fact, it will never be observed—it nevertheless is part of the complete picture of tumor immune system interactions.

2.2.2 Equilibrium

This situation arises once the tumor-free stationary point becomes unstable and there exists a stationary point with small tumor volume x_∗ (and generally up-regulated y_∗) which is stable. We call this stationary point ‘benign’. There are two different scenarios mathematically which correspond to the medical notion of equilibrium: In the simpler one, the benign stationary point is the only stationary point with positive tumor volume and all trajectories converge to it. In the second case, called the bi-stable scenario, there exist three stationary points with positive tumor volumes labelled 0 < x_∗,b < x_∗,u < x_∗,m < x_∞. We call the stationary point with lowest tumor volumes, x_∗,b, benign and the one with highest tumor volume, x_∗,m, malignant. This is merely terminology, but it is somewhat justified by the fact that typically the tumor volume x_∗,b is small with high y_∗,b while x_∗,m is high (close to carrying capacity) with low y_∗,m. Both the benign and malignant stationary points are stable and we call their regions of attraction (i.e., the set of all initial conditions (x₀,y₀,z₀) from which the solution of the dynamics converges to the respective equilibrium point) the benign, respectively malignant regions. The third stationary point x_∗,u is unstable and there exists a surface (its 2-dimensional stable manifold) that passes through it which separates the benign from the malignant region, the so-called stability boundary. (For 2-dimensional systems the terminology separatrix is common (37)). Depending on where the initial condition for the system lies, as time evolves, for the system without any outside interventions the state will converge either to the benign stationary point— and this also corresponds to the medical equilibrium scenario—or it will converge to the malignant stationary point in which case tumor escape occurs. In a rather precise mathematical sense, the bifurcations (changes in stability) which arise as the values for parameters change characterise the transitions from the medical state of equilibrium to the one of tumor escape.

2.2.3 Escape

In addition to the situation just described, it is also possible that there exists just one viable equilibrium point which, however, has high tumor volume x_∗ and low y_∗, i.e., is malignant. In this case, unless somehow by means outside of the modeling done here a change in the values of the parameters can be achieved, after the termination of any treatment the state of the system will always converge to the malignant equilibrium point and it is not possible to revert to the medical condition of equilibrium. In this case, a cure is elusive.

We illustrate the role of the stability boundary in the bi-stable scenario through a 2-dimensional representation in Figure 1. Formally, we have dropped Equation 3 and replaced z in Equation 2 by the equilibrium relation µz = σx. This thus is not directly related to Equations 1–3, but is merely intended to give an illustration of the underlying geometric scenario. Mathematically it will look the same, but not quite as clearly visible, in a higher-dimensional setting. Figure 1 faithfully represents the dynamics in the bi-stable case when both a stable benign and malignant stationary point exist. Obviously, whether or not this is the case depends on the parameter values, but it will hold true for an open set, i.e., for a whole range of values. We note that the same feature is present in the original models by Stepanova (10) and Kuznetsov (11), but also in a more recent in spirit similar 2-dimensional model by Bekker et al. (37) (where a slightly simplified dynamics has been used). In that paper the effects of various immunotherapies on shifting the stability boundary are considered for the model. All these low-dimensional mathematical models clearly point out the stability boundary as the defining structure for tumor immune system dynamics. Such a stability boundary only exists in the bistable scenario and for therapy is the only relevant case.

Figure 1

Figure 1 A 2-dimensional illustration of the bi-stable scenario showing the benign stationary point (x_∗,b,y_∗,b) in green, the saddle point (x_∗,u,y_∗,u) in black and the malignant stationary point (x_∗,m,y_∗,m) in red. The dashed red curve is the stability boundary which separates trajectories which converge to the benign stationary point (above the red curve) from those which converge to the malignant stationary point (below the red curve).

2.3 Formulation of treatment as an optimal control problem in the bi-stable scenario

From a practical point of view, only if the current state of the system (initial condition) is malignant the question of treatment arises. Treatment then should aim to move the state into the benign region, possibly in an efficient way or, in other words, one simply wants to minimize the use of agents to limit side-effects. If ‘tumor escape’ can be reversed to the ‘equilibrium’ condition through therapy by moving the state of the system into the benign region, then—but this assumes that parameter values will not change—after therapy is halted convergence to the benign equilibrium point will occur.

Formulating the problem as an optimal control problem can help with singling out reasonable therapy protocols obtained by minimizing some criterion. For the problem considered here, the initial condition is the present state of the system, desired terminal states are in principle all points in the benign region, the dynamics is given by Equations 1–3, and other constraints that need to be considered are related to the side effects of therapy. Optimal controls then give the protocols on how to administer the therapeutic agents in time which are ‘best’ relative to the chosen criterion. Formulating this criterion is a relevant step in this process. If the objective function does not properly represent the overall goal of therapy results may simply not give beneficial suggestions for therapy protocols. When formulating this criterion to be minimized, the following aspects thus must be taken into account:

(1) Minimizing the objective must induce the system to move into the benign region. For this aim, the geometric shape of the boundary between the benign and malignant regions matters, but there exist many ways to realize such an objective. Here we use a penalty term of the form V_xx(T) + V_yy(T) + V_zz(T) evaluated at the terminal point w(T) where $\overrightarrow{V}=(V_x,V_y,V_z)$ is a suitable vector oriented to point from the benign into the malignant region (as we shall minimize the objective). Generally, one wants to minimize the tumor volume, but it also is the aim to up-regulate the immunocompetent cell density. Thus V_x should be positive while V_y may be allowed to be negative.

(2) Side effects of the therapies have not been included in the modeling. Hence these must now be incorporated indirectly by including penalty terms into the objective which limit the overall amounts of drugs given. The total amount of drugs given are measured by the so-called AUC (‘area under the curve’) in pharmacology. This quantity is given by the integral over the dose rate of the drugs: ${\displaystyle {\int }_0^{T}u(t)dt}$ and ${\displaystyle {\int }_0^{T}v(t)dt}$. Alternatively, a priori constraints on these amounts could be fixed and then the question would be how to best administer these amounts in time. In the literature often quadratic terms are used for the controls for mathematical expediency, but they have no pharmacological meaning.

(3) Mathematically, the existence of a solution needs to guaranteed.

All these considerations led us to formulate the following objective:

$\begin{array}{ll}J=J(u,v)=V_{x}x(T)+V_{y}y(T)+V_{z}z(T)+{\displaystyle {\int }_0^T(Au(t)+Bv(t)+C)dt.}& (4)\end{array}$

The objective Equation 4 is a weighted average of ‘good’ and ‘bad’ terms with the components of the vector $\overrightarrow{V}$ and the coefficients A, B and C weights. These are variables of choice which need to be chosen to strike a balance between the benefit at the terminal time T and the overall side effects. As it is standard in engineering approaches, these coefficients should be calibrated to fine-tune the response of the system.

3 Results

Understanding the geometric properties of the stability boundary gives relevant insights into the possible behavior of ad-hoc chosen therapy protocols. Given a malignant initial condition, it is a more than reasonable strategy to apply chemo- and/or immunotherapy for some time τ, probably chosen by medical guidelines. The time of administration, however, may be crucial. The example in Figure 2 highlights this importance as it shows that the timing of therapy can make a crucial difference and that what might overall constitute a ‘good’ strategy is anything from obvious. The graphs in the figure show the course of two trajectories in (x,y)-space (not their evolution in time) when both chemo- and immunotherapy are applied for time τ and then therapy is stopped. The specific numerical values are irrelevant as we merely want to illustrate a general phenomenon which always exists in the bistable scenario, simply caused by the presence of a stability boundary. The initial segment under therapy is shown as the magenta curve, the subsequent trajectories without treatment are shown in red and green, respectively. For the red portion, τ = 0.72 and this was not sufficient to move the state of the system into the benign region so that convergence to the malignant stationary point occurs. On the other hand, increasing the time just to τ = 0.73 the benign region is reached and subsequently the trajectory converges to the benign stationary point, i.e., the medical condition called equilibrium has been achieved. At the time when therapy is stopped, the states lie close to the stability boundary, but on opposite sides. For a while both trajectories still trace the stability boundary and in each case the tumor volume increases for some time and so does the immunocompetent cell density making it rather impossible to decide whether the course of action was successful or not. Separation of the trajectories only occurs when the state gets close to the unstable saddle point w_∗,u where the instability of the saddle becomes dominant and forces the trajectories to converge to one of the stable stationary points. Only at that time the separation between benign or malignant behavior becomes noticeable. If the state is in the benign region, eventually the reaction of the immune system will be strong enough to control the tumor. This clearly points to the importance of transferring the state of the system well into the benign region. Such an objective is easily incorporated in an optimal control framework, but it requires some knowledge about the geometry of the stability boundary.

Figure 2

Figure 2 Projections of trajectories which administer both chemo and immunotherapy (u ≡ 1 and v ≡ 1) from the initial condition (x₀,y₀,z₀) = (600,0.40,400) (shown as a blue dot) for time τ and then turn off therapy. The initial segment under therapy is shown by a magenta curve, the subsequent trajectories are shown in red (τ = 0.72), respectively green (τ = 0.73). After therapy is stopped both trajectories still show virtually the same behavior for some time (that is, tumor volume and immuno competent sell density evolve the same way in time) before a separation occurs. Only then the red trajectory converges to the malignant stationary point while the green one converges to the benign one. As this simulation demonstrates, a small change in the administration time of the agents can make a big difference and this may not be recognizable for quite some time from the time history of the states.

Considering treatment as an optimal control problem avoids such fallacies as the solutions can be forced to transfer the state well into the benign region. For the problem formulated here a theoretical analysis based on necessary conditions for optimality [e.g., see (38–41)] singles out the following types of therapy protocols as optimal (42, 43): (a) drug administration schedules which alternate between maximum dose-rates and rest-periods, so-called bang-bang controls in optimal control (43) and (b) specific intermittent administrations with particular time-varying dosage protocols at reduced rates. These rates are determined by mathematical formulas for what in optimal control theory are called singular controls. Such protocols can be ruled as being not optimal for immunotherapy based on a mathematical analysis of the problem with tools of geometric optimal control theory (Legendre-Clesch condition (38, 39, 41).

Figure 3 shows a numerically computed candidate bang-bang solution for the problem with a fixed terminal time T. The panel shows the controls u and v and a projection of the corresponding trajectory into (x,y)-space. Segments of the curves which correspond to immunotherapy only (u = 0 and v = 1) are shown as a magenta curve and segments of the curves which correspond to chemotherapy only (u = 1 and v = 0) are shown as a blue curve. The segment of the curve where both chemo- and immunotherapy are at full dose (u = 1 and v = 1) is shown as a brown curve and the segment of the curve where none of chemo- or immunotherapy is active (u = 0 and v = 0) is shown as a black curve. As before, this is merely an episodal illustration.

Figure 3

Figure 3 Example of a numerically computed candidate bang-bang trajectory: The top row shows the controls u (left, blue) and v (right, red) as functions of time while the bottom panel shows the projection of the corresponding trajectory into (x,y)-space. The initial point is marked with a red and the terminal point with a green dot. Black dots show the states on the trajectory where a change in the controls occurs. Initially only immunotherapy is given (magenta segment), then both chemo- and immunotherapy are active at the same time (brown segment) while immunotherapy is stopped and only chemotherapy is given along the blue segments. No drugs are administered along the black segment.

Unfortunately, and contrary to what seems to be claimed in some publications, there are no “fool-proof” numerical algorithms (not to mention software) which safely compute optimal controls for a mathematical problem of the type considered here. The Hamiltonian function for the control problem is not convex in the control and this precludes the use of simple two-point boundary algorithms. Discretization methods are notoriously unreliable when it comes to locating optimal singular controls often arbitrarily declaring that optimal solutions have been found simply when the value of the objective changes by little. The extremal shown in Figure 3 was painstakingly computed solving the very sensitive two-point boundary for bang-bang controls using our own code verifying that the result was an extremal. In line with the theoretical results, immunotherapy follows a bang-bang control and here is active at the beginning to be terminated at the second switching time t₂ = 0.630. The more potent chemotherapy is activated almost immediately at the first switching time t₁ = 0.055 and both therapies are on at full dose until the second switching time. The bulk of the reduction in tumor volume occurs during this time-interval. Chemotherapy is terminated at the third switching time t₃ = 1.195 and then only reactivated briefly at the fourth switching time t₄ = 2.4580 close to the end of the therapy horizon. During the long no treatment phase [t₃, t₄] the tumor volume first increases slightly, but then, as the immunocompetent density increases, starts to decrease again. This behavior is typical for the evolution of a trajectory in a benign region. Therapy then concludes with a brief segment of maximum dose chemotherapy after the prolonged no treatment phase [t₃, t₄]. While there are no claims made that the parameter values underlying this calculation are medically realistic (they were simply used from various sources to illustrate the structure of possible solutions) and the switching times are merely given to convey some sense of the timing for this particular example, both the controls and trajectories follow reasonable patterns.

The graph in Figure 3 represents a situation when chemotherapy is significantly more effective than immunotherapy—this was reflected in the numerical values chosen for this particular computation—and is only meant to be representative for such a scenario. Depending on the severity of side effects (represented in the weights of the objective) and the efficacy of the particular agents, different distributions of the administration of the agents will arise. Various examples of locally optimal bang-bang controls for a related optimal control problem are given in (42). Based on our numerical computations, bang-bang controls (administrations of the agents at full dose with rest periods) are the commonly observed optimal protocols for these types of problems.

4 Discussion

We formulated a 3-dimensional model for tumor-immune system interactions which in an attempt to more closely model the synergies traditional treatments might have with immunotherapies includes a separate equation for tumor antigen. For a mathematical model of tumor-immune interactions to be credible, it is in our opinion a necessary condition that it encompasses, within its range of parameters, the full spectrum of medically known tumor immune system interactions: elimination, equilibrium and tumor escape. This holds for Equations 1–3. Including tumor antigen as a separate state variable is an attempt to model synergistic effects which traditional treatments (chemotherapy considered here) have when combined with immunotherapies.

We reiterate that the model is qualitatively, not quantitative. The aim is to obtain information about the behavior of the dynamics overall, not about some particular situation. For the latter, a rather precise knowledge of the parameter values is required (which simply lies beyond our possibilities). Our emphasis here was on how understanding the behavior of the dynamics can help in the search for optimal dosage protocols. The optimal control problem allows to explore possible therapy protocols in silico suggesting what could be ‘good’ administration protocols relative to some chosen mathematical criterion. Unfortunately, there does not exist off-the-shelf software to solve such optimal control problems reliably and computing optimal controls. As there is great freedom in formulating this objective, there exists the danger of using inadequate criteria for mathematical expediency leading to not only not beneficial but possibly harmful outcomes. Thus we emphasize the need for an a posteriori analysis of the feasibility of the computed therapy protocols. In particular, as side effects of the therapies are only included indirectly while minimizing the objective, the feasibility of the computed strategies from this point of view needs to be checked. It is believed that solutions to these optimal control problems can aid in formulating realistic therapy protocols that can be explored in medical trials and practice.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding author.

Author contributions

UL: Conceptualization, Supervision, Writing – original draft, Writing – review & editing. HS: Conceptualization, Software, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Apetoh L, Ghiringhelli F, Tesniere A, Obeid M, Ortiz C, Criollo A, et al. Toll-like receptor 4-dependent contribution of the immune system to anticancer chemotherapy and radiotherapy. Nat Med (2007) 13):1050–9. doi: 10.1038/nm1622

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Wielkie K, Hahnfeldt P. Modeling the dichotomy of the immune response to cancer:cytotoxic effects and tumor-promoting inflammation. Bull Math Biol (2017) 79:1436–48. doi: 10.1007/s11538-017-0291-4

CrossRef Full Text | Google Scholar

3. Hagemann T, Balkwill F, Lawrence T. Inflammation and cancer: a double-edged sword. Can Cell (2007) 12:300–1. doi: 10.1016/j.ccr.2007.10.005

CrossRef Full Text | Google Scholar

4. Tang C, Wang X, Soh H, Seyedin S, Cortez MA, Krishnan S, et al. Combining radiation and immunotherapy: a new systemic therapy for solid tumors? Cancer Immunol Res (2014) 2:831–8. doi: 10.1158/23266066.CIR-14-0069

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Filatenkov A, Baker J, Mueller A, Kenkel J, Ahn G, Dutt S, et al. Ablative tumor radiation can change the tumor immune cell microenvironment to induce durable complete remissions. Clin Can Res (2015) 21:3727–39. doi: 10.1158/1078-0432.CCR-14-2824

CrossRef Full Text | Google Scholar

6. Golden EB, Demaria S, Schiff P, Chachoua A, Formenti S. An abscopal response to radiation and ipilimumab in a patient with metastatic nonsmall cell lung cancer. Can Immunol Res (2013) 1:365–72. doi: 10.1158/2326-6066.CIR-13-0115

CrossRef Full Text | Google Scholar

7. Poleszczuk JT, Enderling H. The optimal radiation dose to induce robust systemic anti-tumor immunity. Int J Mol Sci (2018) 19(11):3377. doi: 10.3390/ijms19113377

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Poleszczuk JT, Luddy KA, Prokopiou S, Robertson-Tessi M, Moros E, Fishman M, et al. Abscopal benefits of localized radiotherapy depend on activated T-cell trafficking and distribution between metastatic lesions. Cancer Res (2016) 76(5):1009–18. doi: 10.1158/0008-5472.CAN-15-1423

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Enderling H, Chaplain MAJ. Mathematical modeling of tumor growth and treatment. Curr Pharm Design (2014) 20(30):4934–40. doi: 10.2174/1381612819666131125150434

CrossRef Full Text | Google Scholar

10. Stepanova NV. Course of the immune reaction during the development of a Malignant tumour. Biophys (1980) 24:917–23.

Google Scholar

11. Kuznetsov VA, Makalkin IA, Taylor MA, Perelson AS. Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. Bul Math Biol (1994) 56:295–321. doi: 10.1007/BF02460644

CrossRef Full Text | Google Scholar

12. Bellomo N, Preziosi L. Modelling and mathematical problems related to tumor evolution and its interaction with the immune system. Math Comp Model (2000) 32:413–52. doi: 10.1016/S0895-7177(00)00143-6

CrossRef Full Text | Google Scholar

13. Bodnar M, Forys U. A model of immune system with time-dependent immune reactivity. Nonlin Anal: Theory Methods Appl (2009) 70(2):1049–58. doi: 10.1016/j.na.2008.01.031

CrossRef Full Text | Google Scholar

14. Kogan Y, Fory U, Shukron O, Kronik N, Agur Z. Cellular immunotherapy for high grade gliomas: mathematical analysis deriving efficacious infusion rates based on patient requirements. SIAM J Appl Math (2010) 70(6):1953–76. doi: 10.1137/08073740X

CrossRef Full Text | Google Scholar

15. Lai X, Friedman A. Combination therapy of cancer with cancer vaccine and immune checkpoint inhibitors: A mathematical model. PloS One (2017) 12(5):e0178479. doi: 10.1371/journal.pone.0178479

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Lai X, Stiff A, Duggan M, Wesolowski R, Carson W.E.III, Friedman A. Modeling combination therapy for breast cancer with BET and immune checkpoint inhibitors. Proc Natl Acad Sci USA (2018) 115(21):5534–9. doi: 10.1073/pnas.1721559115

PubMed Abstract | CrossRef Full Text | Google Scholar

17. d’Onofrio A. A general framework for modelling tumor-immune system competition and immunotherapy: Mathematical analysis and biomedial inferences. Physica D (2005) 208:202–35. doi: 10.1016/j.physd.2005.06.032

CrossRef Full Text | Google Scholar

18. d’Onofrio A. Metamodeling tumor-immune system interaction, tumor evasion and immunotherapy. Math Comp Model (2008) 47:614–37. doi: 10.1016/j.mcm.2007.02.032

CrossRef Full Text | Google Scholar

19. de Pillis LG, Radunskaya A. A mathematical tumor model with immune resistance and drug therapy: an optimal control approach. J Theo Med (2001) 3:79–100. doi: 10.1080/10273660108833067

CrossRef Full Text | Google Scholar

20. de Pillis LG, Radunskaya A, Wiseman CL. A validated mathematical model of cellmediated immune response to tumor growth. Canc Res (2005) 65:7950–8. doi: 10.1158/0008-5472.CAN-05-0564

CrossRef Full Text | Google Scholar

21. Ledzewicz U, Schättler H. On modeling the synergy of Cancer immunotherapy with radiotherapy. Commun Nonlin Sci Num Simul (2023) 118. doi: 10.1016/j.cnsns.2022.106987

CrossRef Full Text | Google Scholar

22. Ledzewicz U, Naghnaeian M, Schättler H. Optimal response to chemotherapy for a mathematical model of tumor-immune dynamics. J Math Biol (2012) 64:557–77. doi: 10.1007/s00285-011-0424-6

PubMed Abstract | CrossRef Full Text | Google Scholar

23. de Vladar HP, González JA. Dynamic response of cancer under the influence of immunological activity and therapy. J Theor Biol (2004) 227:335–48. doi: 10.1016/j.jtbi.2003.11.012

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Dunn GP, Old LJ, Schreiber RD. The three E’s of cancer immunoediting. Ann Rev Immunol (2004) 22:329–60. doi: 10.1146/annurev.immunol.22.012703.104803

CrossRef Full Text | Google Scholar

25. Dunn GP, Bruce AT, Ikeda H, Old LJ, Schreiber RD. Cancer immunoediting: From immunosurveillance to tumor escape. Nat Immunol (2002) 3:991–8. doi: 10.1038/ni1102-991

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Koebel CM, Vermi W, Swann J, Zerafa N, Rodig S, Old L, et al. Adaptive immunity maintains occult cancer in an equilibrium state. Nature (2007) 450:903–7. doi: 10.1038/nature06309

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Swann JB, Smyth MJ. Immune surveillance of tumors. J Clin Invest (2007) 117:1137–46. doi: 10.1172/JCI31405

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Schättler H, Ledzewicz U. Optimal control for mathematical models of cancer therapies. New York: Springer (2015).

Google Scholar

29. Castiglione F, Piccoli B. Cancer immunotherapy, mathematical modeling and optimal control. J Theo Biol (2007) 247:723–32. doi: 10.1016/j.jtbi.2007.04.003

CrossRef Full Text | Google Scholar

30. Ledzewicz U, Faraji Mosalman MS, Schättler H. Optimal controls for a mathematical model of tumor-immune interactions under targeted chemotherapy with immune boost. Discrete Continuous Dynamical Systems Ser B (2013) 18:1031–51. doi: 10.3934/dcdsb.2013.18.1031

CrossRef Full Text | Google Scholar

31. Jarrett AM, Faghihi D, Hormuth II D, Lima E, Virostko J, Biros G, et al. Optimal control theory for personalized therapeutic regimens in oncology: Background, history, challenges, and opportunities. J Clin Med (2020) 9(5):1314. doi: 10.3390/jcm9051314

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Ledzewicz U, Schättler H. Optimal and suboptimal protocols for a class of mathematical models of tumor anti-angiogenesis. J Theo Biol (2008) 252:295–312. doi: 10.1016/j.jtbi.2008.02.014

CrossRef Full Text | Google Scholar

33. Sharp JA, Browning A, Mapder T, Baker C, Burrage K, Simpson M. Designing combination therapies using multiple optimal controls. J Theor Biol (2020) 497:110277. doi: 10.1016/j.jtbi.2020.110277

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ledzewicz U, Schättler H. On the role of pharmacometrics in mathematical models for cancer treatments. Disc Cont Dyn Syst B (2021) 26:483–99. doi: 10.3934/dcdsb.2020213

CrossRef Full Text | Google Scholar

35. Skipper HE. On mathematical modeling of critical variables in cancer treatment (goals: better understanding of the past and better planning in the future). Bul Math Biol (1986) 48:253–78. doi: 10.1016/S0092-8240(86)90027-3

CrossRef Full Text | Google Scholar

36. Serre R, Benzekry S, Padovani L, Meille C, Andr´e N, Ciccolini J, et al. Mathematical modeling of cancer immunotherapy and its synergy with radiotherapy. Cancer Res (2016) 76:4931–40. doi: 10.1158/0008-5472.CAN-15-3567

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Bekker RA, Zahid MU, Binning JM, Spring B, Hwu P, Pilon-Thomas S, et al. Rethinking the immunotherapy numbers game. J ImmunoTher Cancer (2022) 10(7):2022. doi: 10.1136/jitc-2022-005107

CrossRef Full Text | Google Scholar

38. Bonnard B, Chyba M. Singular trajectories and their role in control theory. In: Mathématiques & Applications, vol. 40. Paris: Springer-Verlag (2003).

Google Scholar

39. Bressan A, Piccoli B. Introduction to the mathematical theory of control. Springfield, Mo, USA: AIMS (2007).

Google Scholar

40. Pontryagin LS, Boltyanskii VG, Gamkrelidze RV, Mishchenko EF. The mathematical theory of optimal processes. New York: Macmillan (1964).

Google Scholar

41. Schättler H, Ledzewicz U. Geometric optimal control. New York: Springer (2012).

Google Scholar

42. Ledzewicz U, Maurer H, Schättler H. Bang-bang optimal controls for a mathematical model of chemo- and immunotherapy in cancer. Discr Cont Dyn Syst B (2023). doi: 10.3934/dcdsb.2023141

CrossRef Full Text | Google Scholar

43. Ledzewicz U, Schättler H. On the optimal control problem for a model of the synergy of chemo- and immunotherapy. Opt Cont Meth Appl (2023). doi: 10.1002/oca.3016

CrossRef Full Text | Google Scholar