Control Theory and Cancer Chemotherapy: How They Interact

Control theory arises in most modern real-life applications, not least in biological and medical applications. In particular, in biological and medical contexts, the role of control theory began to take shape in the early 1980s when the first works appeared on the application of control theory in models of pharmacokinetics and pharmacodynamics for antitumor therapies. Forty years after those first works, the theory of control continues to be considered a mathematical analysis tool of extreme importance and usefulness, but the challenges it must overcome in order to manage the complexity of biological processes are in fact not yet overcome. In this article, we introduce the reader to the basic ideas of control theory, its aims and its mathematical formalization, and we review its use in cell phase-specific models for cancer chemotherapy. We discuss strengths and limitations of the control theory approach to the analysis pharmacokinetics and pharmacodynamics models, and we will see that most of them are strongly related to data availability and mathematical form of the model. We propose some future research directions that could prove useful in overcoming the these limitations and we indicate the crucial steps preliminary to a useful and informative application of control theory to cancer chemotherapy modeling.

1. Introduction

Control theory is concerned with establishing whether the evolution of a system is controllable, i.e., whether the evolution can be influenced by some external agent, called control. Optimal control theory deals with finding a control for the system over a period of time such that a performance criterion is optimized. Control theory (and its dual, i.e., observability theory) was applied to the analysis of deterministic models of biological and biochemical processes (Lecca and Re, 2019; Thomas et al., 2019; Wu et al., 2019), and in particular to models of cancer chemotherapy. Mathematically, cancer chemotherapy can be formulated as an optimal control problem. This approach is more relevant than ever today, as testified by recent works (Schättler and Ledzewicz, 2015a,b; Wang and Schättler, 2016; Angaroni et al., 2020; Jarrett et al., 2020; Sweilam et al., 2020). In literature, applications of optimal control theory to mathematical models of cancer biology and role of chemotherapy began to appear in the 1980s and have appeared with regularity in the following years to the present day. The reader can refer to some of these works, such as for example (Swan, 1980, 1988, 1990; Martin, 1992; Swierniak et al., 1996; Kimmel and Swierniak, 2005; Pillis et al., 2008; Collins et al., 2010; Batmani and Khaloozadeh, 2011; Ledzewicz et al., 2013; Ghaffari et al., 2014; Michor and Beal, 2015; Wang and Schättler, 2016; Carrère, 2017; Irurzun-Arana et al., 2018). Indeed, a tumor undergoing pharmacological treatment can be viewed as a control system with the state of the system, given by the number of cancer cells at time t, N(t), and the control input at time t, u(t). Usually in these models u(t) is the drug dosage or the effect the drugs have on normal and cancer cells. Since chemotherapeutic agents affect both healthy cells and cancer cells, the objective of the control problem is to minimize the number of cancer cells keeping the toxicity to the normal tissue at a safe level. The greater usefulness of the optimal control model of cancer chemotherapy lies in the capacity of determining drug schedules that most effectively reduce the size of a tumor after a fixed period of treatment has elapsed (Martin, 1992; Ledzewicz and Schaettler, 2004).

Mathematical models for cancer chemotherapy are necessarily cell-cycle specific, because most drugs are active in a specific phase of the cell-cycle as reported in Ledzewicz and Schaettler (2004), Schättler and Ledzewicz (2015a). The starting point of cell cycle is a growth phase G1 after which the cell enters a phase S where DNA synthesis occurs. Then a second growth phase G2 takes place in which the cell prepares for mitosis (phase M) during which cell division occurs. Each of the two daughter cells can either re-enter phase G1 or for some time may simply lie dormant (phase G0) until re-entering G1, thus re-starting the entire process. The molecules of the anticancer drug treatments can be cell-killing agents, but also blocking and recruitment agents. Blocking agents slow down the transitions of the cells through the cell cycle, i.e., they slow down the tumor growth; recruiting agents make cancer cells leave the dormant stage G0 where they are not susceptible to any chemotherapy (Ledzewicz and Schaettler, 2004; Schättler and Ledzewicz, 2015a,b). Usually, in these models the phases of the cell cycle are clustered into compartments with the state representing the average number of cells in each compartment and the control representing the dosages or effects of the various drugs. The number of compartments depends on the number and types of chemotherapeutic agents considered, since as we said a chemotherapeutic agents may acts on one or more phases. The in- and out-flows between the compartments in the presence of the control, determine the system dynamics. We consider extremely important as a future line of progress in this field, the construction of a pharmacokinetic/pharmacodynamic model that is less abstract than those usually used and which mainly considers two variables, such as drug dosage and the number of cells exposed to them. In particular, we are thinking of a model that describes the interactions at the molecular level of the drug with the biological networks that inside the cell are responsible for its processing and with the molecular targets for which it is designed.

2. The Optimal Control Problem

The formulation of an optimal control problem consists of (i) a mathematical description of the problem to be controlled; (ii) a statement of the physical constraints, and (iii) the specification of a performance criterion (Kirk, 1970). We restrict our focus to models of ordinary differential equations. Let x₁(t), x₂(t), …, x_n(t) be the state variables of the system at time t, and let u₁(t), u₂(t), …, n_m(t) the inputs to the systems at time t. The system can be described by n differential equations

$\begin{array}{l}\frac{d{x}_1}{dt}=f_1(x_1(t),x_2(t),\dots ,x_n(t),u_1(t),u_2(t),\dots ,u_m(t))\\ \frac{d{x}_2}{dt}=f_2(x_1(t),x_2(t),\dots ,x_n(t),u_1(t),u_2(t),\dots ,u_m(t))\\ ⋮\\ \frac{d{x}_n}{dt}=f_n(x_1(t),x_2(t),\dots ,x_n(t),u_1(t),u_2(t),\dots ,u_m(t))\\ \end{array}$

where f_i : ℝ → ℝ (i = 1, 2, …, n). In matrix form, we write

$\begin{array}{ll}\frac{d\text{u}}{dt}=\text{f}(\text{x}(t),\text{u}(t),t)& (1)\end{array}$

where

$\begin{array}{ll}\text{x}(t)=\left[\begin{array}{c}{x}_1(t)\\ x_2(t)\\ ⋮\\ x_n(t)\end{array}\right]\text{u}(t)=\left[\begin{array}{c}{u}_1(t)\\ u_2(t)\\ ⋮\\ u_m(t)\end{array}\right]& \text{f}(t)=\left[\begin{array}{c}{f}_1(\text{x}(t),\text{u}(t),t)\\ f_2(\text{x}(t),\text{u}(t),t)\\ ⋮\\ f_n(\text{x}(t),\text{u}(t),t)\end{array}\right].\end{array}$

The explicit presence of t in the right-hand side of Equation (1) indicates that the system is time-variant. For a time-invariant system instead we will write

$\begin{array}{ll}\frac{d\text{u}}{dt}=\text{f}(\text{x}(t),\text{u}(t)).& (2)\end{array}$

If a system is time-varying and linear its state equations are

$\begin{array}{ll}\frac{d\text{x}}{dt}=\text{A}(t)\text{x}(t)+\text{B}(t)\text{u}(t)& (3)\end{array}$

where A(t) and B(t) are, respectively, n × n and n × m matrices with time-dependent elements. State equations for linear time-invariant systems are instead of the form

$\begin{array}{ll}\frac{d\text{x}}{dt}=\text{A}\text{x}(t)+\text{B}\text{u}(t)& (4)\end{array}$

where A and B are time-independent matrices.

A history of the control input u(t) in the interval [t₀, t_f] is called a control history. A history of the state variables x(t) in the interval [t₀, t_f] is called a state trajectory. The physical constraints are defined by the boundary conditions, i.e., the values of the state variables and the control input at the initial time t₀ and t_f. A state trajectory that satisfies the state variable constraints during the entire time interval [t₀, t_f] is called an admissible trajectory. A control history that satisfies the control constraints during the entire time interval [t₀, t_f] is called an admissible control. The admissibility is a useful property because it reduces the number of values that can be assumed by the state variable and the controls. Thanks to this property, instead of investigating all controls and trajectories histories to identify which is the optimal ones, we can analyse only those controls and trajectories that are admissible (Kirk, 1970).

The state variables that can be measured are called outputs. Let y₁(t), y₂(t), …, y_p(t) denote the outputs, then

$\begin{array}{ll}\text{y}=\text{g}(\text{x}(t),\text{u}(t),t).& (5)\end{array}$

If the outputs are related to the states and controls by a time-varying relationship then

$\begin{array}{ll}\text{y}=\text{C}(t)\text{x}(t)+\text{D}(t)\text{u}(t)& (6)\end{array}$

where C and D are p × n and p × m time-dependent matrices. If the systems is time-invariant then C and D are time-independent matrices.

Once the model and the physical constraints are specified, last step in the formulation of the optimal control problem is the definition of a performance measure useful to quantitatively assess the performance of the system. An optimal control is defined as one that minimizes (or maximizes) the performance measure. The optimal control problem consists in finding an admissible control u_a which causes the system described by Equation (1) to follow an admissible trajectory x_a that minimizes the performance measure

$\begin{array}{ll}J=h_1(\text{x}(t_f),t_f)+{\displaystyle {\int }_{t_0}^{t_f}}{h}_2(\text{x}(t),\text{u}(t),t)dt& (7)\end{array}$

where h₁ and h₂ are scalar functions. Starting from the initial state x(t₀) = x₀ and applying a control u(t) for t ∈ [t₀, t_f] causes the systems to follow some trajectory. The performance measure J assign a real number to each trajectory. We refer the reader to Moore (2018), Lopes et al. (2009), Sargent (2000), and Alekseev et al. (1987) for a review about necessary optimality conditions.

Two important notes are duty here: (i) before attempting to determine the controls we have to make sure that the system is controllable, since we may not know in advance if a system is controllable and if an optimal control exists, and (ii) even if an optimal control exists, it may not be unique. In the first case we can resort to controllability theorems. In the second case, non-unique optimal controls may complicate the computational procedures, but they give the possibility of choosing among different control configurations (Kirk, 1970).

The general solution of Equation (4) is

$\begin{array}{ll}\text{x}(t)=e^{\text{A}t}[\text{x}(t_0)+{\displaystyle {\int }_{t_0}^t}{e}^{-\text{A}\tau }\text{B}\text{u}(\tau )d\tau ].& (8)\end{array}$

Kalman (1960) has shown that a linear time-invariant system is controllable if and only if the n × nm matrix

$\begin{array}{ll}\text{E}\equiv [\text{B}|\text{A}\text{B}|{\text{A}}^2\text{B}|\cdots |{\text{A}}^{n-1}\text{B}]& (9)\end{array}$

has rank n. This condition guarantees that there are a finite time t₁ ≥ t₀, and a control u(t) (t ∈ [t₀, t_f]) which moves the state x₀ = x(t₀)) to time t₁ (Kirk, 1970). In the following section, we specialize these concepts to a general dynamic model of cell-cycle specific antitumor treatments.

3. Cell-Cycle Specific Dynamic Model

We consider a multi-compartment and multi-drug model, that according to Ledzewicz and Schättler (Schättler and Ledzewicz, 2015a) is described by a bilinear system of the form

$\begin{array}{ll}\frac{dN}{dt}=N(A+{\displaystyle \sum _{i=1}^m}{u}_i{B}_i),N_0=N(0)& (10)\end{array}$

where N = N(t) is the number of cancer cell at time t, m is the number of drugs, A and B are n × n matrices, n is the number of compartments, and u_i represent the dosages of various drugs. The drugs performs a cytostatic action which kills both cancer and healthy cells.

A model for the objective function is given by the formalization of the aim therapeutic treatment, which is to kill the cancer cells while keeping the toxicity to the health cells acceptable. In literature there are many non-equivalent ways of modeling the objective function (Swan, 1990; Swierniak et al., 1996; Schättler and Ledzewicz, 2015b; Wang and Schättler, 2016; Angaroni et al., 2020). A model that is largely adopted was proposed by Świerniak (1995), Świerniak et al. (2016) with the objective to minimize the number of cancer cells at the end of a fixed therapy interval. In the context of this modeling, Ledzewiczs and Scättler (Schättler and Ledzewicz, 2015b) proposed a linear L1-type objective function of the form

$\begin{array}{ll}J=pN(T)+{\displaystyle {\int }_0^T}(qN(t)+bu(t))dt& (11)\end{array}$

where p = (p₁, p₂, …, p_n), q = (q₁, q₂, …, q_m), and b = (b₁, b₂, …, b_m). are non-zero row vectors of non-negative coefficients. The components of b corresponding to killing chemotherapeutic agents are positive. The term pN(T) is a weighted average of the total number of cancer cell at the end of a fixed therapy interval [0, T]. The term qN(t) prevents that the number of cancer cells rise to unacceptably high values at intermediate times. In this scheme, the toxicity of the drugs is only modeled indirectly through the term bu(t), i.e., the adverse effect on the healthy cells are represented only indirectly by minimizing the drug dosage u(t) in the objective function.

Usually the pharmacokinetics equation for the body/plasma drug concentration c(t) if a first-order differential equation as the following

$\begin{array}{ll}\frac{dc}{dt}=-(k_1+k_{2}u(t))c(t)+hu(t),c(0)=0& (12)\end{array}$

where k₁ and h are two positive constants, but k₂ is arbitrary.

Finally, the equations of the dynamics is

$\begin{array}{l}\frac{dN}{dt}=(A+cB)N,N(0)=N_0\\ \frac{dc}{dt}=-[k_1+k_{2}u(t)]c(t)+hu(t),c(0)=0.\end{array}$

This model is still incomplete, and thus unrealistic, since it lacks of the pharmacodynamics that takes into account the effect the drug concentration c(t) has on the cancer cells.

Generally, the effect can be modeled by a function

$s:[0,\infty )\to \left[0,1\right].$

Different models arise depending on the choice of the function s. Two commonly used models in practice are the so-called E_max model of the form

$\begin{array}{ll}{s}_1(c)=E_0+\frac{E_{\text{max}}c}{E{C}_{50}+c}& (13)\end{array}$

and the sigmoidal function

$\begin{array}{ll}{s}_2(c)=E_0+\frac{E_{\text{max}}{c}^k}{E{C}_{50}^k+c^k}& (14)\end{array}$

where E₀ is the baseline of the drug efficacy, E_max if the maximum efficacy, EC₅₀ is the potency and k is a positive integer greater than 1 (Macdougall, 2006; Choe and Lee, 2017). In the E_max model it is assumed that the drug becomes almost immediately effective, but then saturates at high concentrations. In the sigmoidal models the effectiveness at both lower and higher concentrations are more accurately approximated.

Ultimately, we obtain the following equations for dynamics

$\begin{array}{ll}\frac{dN}{dt}=[A+s(c)B]N,N(0)=N_0& (15)\end{array}$

$\begin{array}{ll}\frac{dc}{dt}=-[k_1+k_{2}u(t)]c(t)+hu(t),c(0)=0& (16)\end{array}$

and the objective function J remains the same as in Equation (11).

4. Model Simulation

We simulated the model given in Equations (15)–(16) considering one compartment only (A and B are 1 × 1 matrices, i.e., scalar variables) and by setting the following values of the parameters expressed in arbitrary units: E₀ = 0, E_max = 100, k₁ = 0.005, k₂ = 0.0004, EC₅₀ = 15, h = 0.001, k = 1, A = 0.001, B = 0.001, and a time periodic function for the drug dosage u(t) = 10 cos(0.1t) + 10, to simulate a periodic administration of the drug. The values of E_max, and EC₅₀ have been taken as in Felmlee et al. (2012). The parameters q, b and p in the Equation (11) for J have been set equal to 1. The R code implementing the simulation of the models and the parameter sensitivity analysis is reported in Tables 1 and 2, respectively for the sake of results' reproducibility. The numerical solution of the equation systems if showed in Figure 1. The concentration curve shows synchronous fluctuations to the boulders and the minima of the u(t) curve and after a certain time interval it assumes a constant trend as long as the drug continues to be administered. In this model, if you stop administering the drug, the concentration curve assumes a decreasing trend after a certain time interval.

TABLE 1

Table 1. R script implementing the model described by the Equations (15)–(16).

TABLE 2

Table 2. R functions implementing the sensitivity analysis of the model of Equations (15)–(16).

FIGURE 1

Figure 1. Numerical simulation of the system defined by Equations (15) and (16) with E₀ = 0, E_max = 100, k₁ = 0.005, k₂ = 0.0004, EC₅₀ = 15, h = 0.001, k = 1, A = 0.001, B = 0.001. As shown later, numerical simulations for times greater than 500 time units, shows that after the first 500 time units the oscillatory behavior disappears to give way to a plateau.

For the parameters A, B, k₁, k₂, and k we considered ranges of values (see Table 3) and carried out global sensitivity analysis of the model using the functions of the R library FME (Soetaert and Petzoldt, 2010). Figure 2 shows the sensitivity ranges of the concentration curve c(t) to the parameters listed in Table 3. The analysis shows that the model is sensitive to k₁, k₂, and h, and not sensitive to k, A and B. In Figure 3 sensitivity ranges of the concentration c(t) are reported for the combination of parameters k₁, k₂ and h.

TABLE 3

Table 3. Parameters' ranges used for the global sensitivity analysis of the model defined by Equations (15) and (16).

FIGURE 2

Figure 2. Sensitivity range of the concentration curve to the parameters listed in Table 3 (see text for R-code).

FIGURE 3

Figure 3. Sensitivity ranges of the concentration c(t) for 100 combination of parameters k₁, k₂ and h. We select the Latin hypercube sampling algorithm (see the R code in Table 2).

These ranges have been obtained with the procedure reported in Hearne (1985), where we selected the latin hypercube sampling algorithm (see the R code in Table 2).

We focus our analysis on k₂ (similar conclusions can be drawn for the parameter h since it is as k₂ a multiplicative factor of u(t)). The higher the positive k₂ value, the sooner the curve c(t) reaches the saturation point and the lower the concentration saturation value is in Figure 4, we see that how k₂ affects the value of J: a small increment of k₂ induces a steep decrement of J.

FIGURE 4

Figure 4. The value of objective function j rapidly decreases as k₂ increases.

An interesting case arises for negative values of k₂. For example, if we consider k₂ = −0.0004 (instead of k₂ = +0.0004), the curve of c(t) looks like in Figure 5: after a time interval in which it shows slight fluctuations corresponding to the fluctuations of u(t), its growth accelerates and the fluctuations become almost nil. In case of negative values of k₂ the control and the optimal control are not guaranteed. In this regard, Ledzewicz et al. (2013) comment saying for the bilinear model the sign k₂ matters and for k₂ > 0 singular controls satisfy the necessary conditions for optimality. There is an intuitive explanation for this: once the drug's concentration is built up, in this case the injection of smaller time-varying doses can be used to maintain a high effectiveness of the drug which by itself slowly decays.

FIGURE 5

Figure 5. Numerical simulation of the system defined by Equations (15) and (16) with E₀ = 0, E_max = 100, k₁ = 0.005, k₂ = −0.0004, EC₅₀ = 15, h = 0.001, k = 1, A = 0.001, B = 0.001. A negative value of k₂ causes the concentration curve to assume an increasing monotone trend, while a positive value of k₂ causes the curve to reach a plateau after an initial oscillatory behavior.

We finally note that the scarce influence of parameters A and B on the model reveals the scarce influence of the compartment of the model in which the drug-cell interactions take place. This fact is a limitation of these models, since the processes through which the drug performs its functions are different in different phases of the cell cycle.

5. Potential Future Developments in the Field

Considering the model described in the previous section, it may look like we always know the explicit formulas for the control(s), but, unfortunately, this is not the case in most of the applications of control theory to the experimental sciences and especially to biological sciences. Usually, a guess of the analytical form of u(t) comes from the first experimental regime designed for drug administration scheduling. However, the example shown in the previous section shows that the knowledge of the mathematical form for u(t) is not enough if an accurate knowledge of its parameters is missing. We indeed observed that the controllability and also the optimality of the control strongly depend on the parameters of u(t) (k₂ and h are two multiplicative factors of u(t)).

As Moore notes in Moore (2018), a frequent assumption in pharmacokinetic/pharmacodynamic modeling is that the structure of the equations is the same in animals and humans, and only the parameter values differ. In the case in which parameter values in a model are obtained by fitting the model to animal preclinical data. Allometric scaling is used to predict corresponding parameter values for a human population (Knibbe et al., 2005; Mahmood, 2009; Huang and Riviere, 2014). In spite of the fact that allometric scaling is a simple potential tool and rational option for the estimate of pharmacokinetic parameters in species for which there are no data available, is has no the same predictive power of the use of human data samples. if data from human cohorts are available, the best way to obtain the parameter values is by directly fitting the model to human clinical data. In addition to tying models closely to data, evaluation of a predicted optimal regimen with experimental data is very essential. Study outcomes from regimens that are predicted to be optimal should be compared with those from standard regimens (Moore, 2018). Such evaluation tests whether assumptions in the model, the objective functional, and the parameter estimates, are plausible and lead to desired outcomes (i.e., the model is optimally controllable by u(t)).

The experimental determination of a possible model for u(t) and the estimation of its parameters are not the last steps to be taken before applying the control theory to optimize or discard the proposed regime for u(t). Performing sensitivity analysis on the model is an equally important step that we suggest to perform before any attempt to fit the model to the available data. As we have seen in the simple example of the previous section, sensitivity analyses is crucial to identify parameters that the outcome is largely insensitive. The parameters that do not affect the model can be fixed, while the other can be estimated by fitting the model to data or from the literature. Finally, however, the presence of parameters to whose variations the model is insensitive deserve out attention, and has not to be welcome as the possibility to reduce the dimensions of the free parameter space. It could indicate the inadequacy of the model used to describe the kinetics and dynamics of the drug under investigation. This is the case for example of similar models to the one presented in this study, which in fact is representative of a populated class of similar models. The specificity of the kinetics and dynamics at the cell cycle phase that the chemotherapist demonstrates in vivo is not adequately described by these models. In such models, the controllability analysis does not depend on the actual values of these parameters, but it is the type of pharmacokinetics model which determines the class of optimal controls, as noted also by Ledzewicz et al. (2013). In this case the issue about the proper choice of pharmacokinetics in the cancer chemotherapy problems arises and become an essential item in the modeling of cancer chemotherapy. In fact, it must always be remembered that when we talk about controllability we always refer to the controllability of a model. If the model turns out to be very close to reality then, the terms model and system can be used interchangeably and it can be said whether the system (and not only the model) is controllable or not.

6. Conclusions

The techniques and examples in this article are intended to present to practitioners and modelers of pharmacokinetics/pharmacodynamics the basic concepts of the optimal control theory and its application to optimize drug regimens.

Considering a simple but widely used example we show that the prediction about optimal controllability depends on the formal characteristics of the model equation and the knowledge or the possibility to infer their parameters. Controllability is primarily a property of the model, which becomes a property of the system if the model is accurately realistic. By integrating clinical data with data concerning the molecular processes of drug metabolism, mathematical models that are able to accurately capture the drug kinetics and dynamics as well as the dynamics of the disease states, there are advantageous opportunities to take advantage of the application of optimal control theory. The final goal is to plan therapeutic regimens for preclinical or clinical use to predict optimal therapeutic regimens.

Until a few years ago the clinical application of mathematical models of tumor growth were limited because they require input data difficult to obtain with sufficient spatial resolution in patients even at a single time point—for example, extent of vascularization, immune infiltrate, ratio of tumor-to-normal cells, or extracellular matrix status (Yankeelov et al., 2013). Nowadays, however, the availability of data from advanced quantitative tumor imaging methods rekindles the interest in the development of mathematical models and paves the way to an effective application of control theory to predict optimal therapy delivery parameters and ultimately to control the response to treatment and the time to progression. Indeed, the number of publications on mathematical modeling of cancer is growing at an exponential rate (Brady and Enderling, 2019), and for this reason it is urgent to consider the application of control theory to the analysis of the available models to evaluate their effectiveness as tools to aid in the programming of an optimal dosage and drug administration scheduling.

In this direction, there is also another particularly important aspect that control theory must incorporate, namely the development of the chemo-resistance. In a recent study concerning the abiraterone treatment for metastatic castrate-resistant prostate cancer, Gatenby et al. (Zhang et al., 2017a) noted that with standard dosing, evolution of resistance with treatment failure (testified by radiographic progression) occurs at a median of about 16.5 months.

The authors observed that the conventional treatment strategy, which administers cytotoxic drugs at maximum tolerated dose until progression, can be evolutionarily unwise since it strongly selects for resistant phenotypes and eliminates potential competitors. These Darwinian dynamics can lead to rapid proliferation of resistant populations. On the basis of this idea, the authors developed an evolutionary game theory model of Lotka–Volterra equations with three competing cancer “species”: androgen dependent, androgen producing, and androgen independent. As expected, the model simulations with standard abiraterone dosing demonstrated strong selection for androgen-independent cells and rapid treatment failure. Gatenby et al. then hypothesized that the time to progression could be increased by integrating evolutionary dynamics into therapy scheduling models. In order to implement this integration, they framed the population dynamics before and during abiraterone therapy using a game theoretic model built on evolutionary first principles with parameter estimates inferred from clinical data. The results of Gatenby et al. model (Zhang et al., 2017a) demonstrated how a therapy treatment strategy synchronized with intra-tumoral evolutionary dynamics may prolong time to progression and significantly decrease total drug dose. The evolution-based strategy suggested by the model was then applied in pilot clinical trials.

Cancer is a complex multiscale process (Zhang et al., 2008; Deisboeck et al., 2011), whose modeling requires the integration of different type of data (Deisboeck et al., 2011). this raises the problem of developing control theory approaches to multiscale systems. Building a model for a multiscale system is generally challenging. In particular, a multiscale model of tumoral growth in response to a therapy, can include (i) equations spanning a wide range of time and length scales, (ii) stochasticity inherent to the intracellular and intercellular biochemical interactions, and, consequently, non-continuum and dynamically coupled continuum-non-continuum processes description. To all this is added a comprehensive knowledge of the physicochemical mechanisms as well as the values of thermodynamic and kinetic parameters of the biochemical reactions involving the drug and its target molecules. Each of these points poses a challenge to the control theory applied to tumoral growth models. Overcoming each of these challenges in a short time may be possible if in the near future we will be able to develop a control theory suitable for the analysis of controllability of multiscale hybrid stochastic systems.

From the theoretical point of view, in the last decade works have been published that expose the mathematical foundations of control theory or that do discuss the perspectives of its application to the case of multiscale systems (see for example (Christofides and Armaou, 2006; Christofides et al., 2009; Hartmann et al., 2014; Mei et al., 2016), and to stochastic systems (see for example (Braatz et al., 2006; Nisio, 2015; Zhang et al., 2017b). What is currently lacking and which could prove to be of great use for the design of dosages and administration regimens of an anticancer drug is the construction of a control theory of control for hybrid stochastic multiscale systems. The modeling and control of stochasticity could also provide further tools for controlling the evolution of drug resistance (as an extension of evolutionary dynamics models) and implementing appropriate countermeasures.

Author Contributions

PL has carried out the bibliographic research necessary to review the state of the art of the application of control theories to pharmacokinetic/pharmacodynamic models and is the author of the design and R codes that implement the analysis of these models.

Funding

This study has been supported by the project COMPANET (RTD 2019), funded by the Free University of Bozen-Bolzano, Italy.

Conflict of Interest

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

References

Alekseev, V. M., Tikhomirov, V. M., and Fomin, S. V. (1987). Optimal Control. Part of the Contemporary Soviet Mathematics Book Series. New York, NY: Springer Science+Business Media. doi: 10.1007/978-1-4615-7551-1

CrossRef Full Text | Google Scholar

Angaroni, F., Graudenzi, A., Rossignolo, M., Maspero, D., Calarco, T., Piazza, R., et al. (2020). An optimal control framework for the automated design of personalized cancer treatments. Front. Bioeng. Biotechnol. 8:523. doi: 10.3389/fbioe.2020.00523

PubMed Abstract | CrossRef Full Text | Google Scholar

Batmani, Y., and Khaloozadeh, H. (2011). “Multi objective optimization of drug regimens in cancer chemotherapy using a PK-PD model,” in 2011 19th Iranian Conference on Electrical Engineering (Tehran), 1.

Google Scholar

Braatz, R., Alkire, R., Seebauer, E., Rusli, E., Gunawan, R., Drews, T., et al. (2006). Perspectives on the design and control of multiscale systems. J. Process Control 16, 193–204. doi: 10.1016/j.jprocont.2005.06.001

CrossRef Full Text | Google Scholar

Brady, R., and Enderling, H. (2019). Mathematical models of cancer: When to predict novel therapies, and when not to. Bull. Math. Biol. 81, 3722–3731. doi: 10.1007/s11538-019-00640-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Carrère, C. (2017). Optimization of an in vitro chemotherapy to avoid resistant tumours. J. Theoret. Biol. 413, 24–33. doi: 10.1016/j.jtbi.2016.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Choe, S., and Lee, D. (2017). Parameter estimation for sigmoid emax models in exposure-response relationship. Transl. Clin. Pharmacol. 25:74. doi: 10.12793/tcp.2017.25.2.74

PubMed Abstract | CrossRef Full Text

Christofides, P. D., and Armaou, A. (2006). Control and optimization of multiscale process systems. Comput. Chem. Eng. 30, 1670–1686. doi: 10.1016/j.compchemeng.2006.05.025

CrossRef Full Text | Google Scholar

Christofides, P. D., Armaou, A., Lou, Y., and Varshney, A. (2009). Control and Optimization of Multiscale Process Systems. Boston, MA: Birkhäuser. doi: 10.1007/978-0-8176-4793-3

CrossRef Full Text | Google Scholar

Collins, C., Fister, K., and Williams, M. (2010). Optimal control of a cancer cell model with delay. Math. Modell. Nat. Phenomena 5, 63–75. doi: 10.1051/mmnp/20105305

CrossRef Full Text | Google Scholar

Deisboeck, T. S., Wang, Z., Macklin, P., and Cristini, V. (2011). Multiscale cancer modeling. Annu. Rev. Biomed. Eng. 13, 127–155. doi: 10.1146/annurev-bioeng-071910-124729

CrossRef Full Text | Google Scholar

Felmlee, M. A., Morris, M. E., and Mager, D. E. (2012). “Mechanism-based pharmacodynamic modeling,” in Methods in Molecular Biology (Totowa, NJ: Humana Press), 583–600. doi: 10.1007/978-1-62703-050-2_21

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghaffari, A., Nazari, M., and Arab, F. (2014). Optimal finite cancer treatment duration by using mixed vaccine therapy and chemotherapy: state dependent Riccati equation control. J. Appl. Math. 2014, 1–9. doi: 10.1155/2014/363109

CrossRef Full Text | Google Scholar

Hartmann, C., Latorre, J. C., Zhang, W., and Pavliotis, G. A. (2014). Optimal control of multiscale systems using reduced-order models. J. Comput. Dyn. 1, 279–306. doi: 10.3934/jcd.2014.1.279

CrossRef Full Text | Google Scholar

Hearne, J. (1985). Sensitivity analysis of parameter combinations. Appl. Math. Modell. 9, 106–108. doi: 10.1016/0307-904X(85)90121-0

CrossRef Full Text | Google Scholar

Huang, Q., and Riviere, J. E. (2014). The application of allometric scaling principles to predict pharmacokinetic parameters across species. Expert Opin. Drug Metab. Toxicol. 10, 1241–1253. doi: 10.1517/17425255.2014.934671

PubMed Abstract | CrossRef Full Text | Google Scholar

Irurzun-Arana, I., Janda, A., Ardanza-Trevijano, S., and Trocóniz, I. F. (2018). Optimal dynamic control approach in a multi-objective therapeutic scenario: application to drug delivery in the treatment of prostate cancer. PLoS Comput. Biol. 14:e1006087. doi: 10.1371/journal.pcbi.1006087

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalman, R. (1960). On the general theory of control systems. IFAC Proc. Vol. 1, 491–502. doi: 10.1016/S1474-6670(17)70094-8

CrossRef Full Text | Google Scholar

Kimmel, M., and Swierniak, A. (2005). “Control theory approach to cancer chemotherapy: benefiting from phase dependence and overcoming drug resistance,” in Tutorials in Mathematical Biosciences III Lecture Notes in Mathematics, Vol. 1872, ed A. Friedman (Berlin; Heidelberg: Springer), 185–221. doi: 10.1007/11561606_5

CrossRef Full Text | Google Scholar

Kirk, D. (1970). Optimal Control Theory; An Introduction. Englewood Cliffs, NJ: Prentice-Hall.

Google Scholar

Knibbe, C. A. J., Zuideveld, K. P., Aarts, L. P. H. J., Kuks, P. F. M., and Danhof, M. (2005). Allometric relationships between the pharmacokinetics of propofol in rats, children and adults. Brit. J. Clin. Pharmacol. 59, 705–711. doi: 10.1111/j.1365-2125.2005.02239.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lecca, P., and Re, A. (2019). Identifying necessary and sufficient conditions for the observability of models of biochemical processes. Biophys. Chem. 254:106257. doi: 10.1016/j.bpc.2019.106257

PubMed Abstract | CrossRef Full Text | Google Scholar

Ledzewicz, U., Mosalman, M. S. F., and Schättler, H. (2013). Optimal controls for a mathematical model of tumor-immune interactions under targeted chemotherapy with immune boost. Discrete Contin. Dyn. Syst. B 18, 1031–1051. doi: 10.3934/dcdsb.2013.18.1031

CrossRef Full Text | Google Scholar

Ledzewicz, U., and Schaettler, H. (2004). “Application of control theory in modelling cancer chemotherapy,” in Proceedings of the International Conference on Control, Automation and Systems (ICCAS) (Bangkok: The Shangri-La Hotel), 330–335.

Google Scholar

Lopes, S. O., Fontes, F. A. C. C., Simos, T. E., Psihoyios, G., and Tsitouras, C. (2009). “Necessary conditions of optimality for calculus of variations problems with inequality constraints,” in AIP Conference Proceedings, Vol. 1168 (Crete: Rethymno). doi: 10.1063/1.3241344

CrossRef Full Text | Google Scholar

Macdougall, J. (2006). “Analysis of dose-response studies-emax model,” in Statistics for Biology and Health (New York, NY: Springer), 127–145. doi: 10.1007/0-387-33706-7_9

CrossRef Full Text | Google Scholar

Mahmood, I. (2009). Pharmacokinetic allometric scaling of antibodies: application to the first-in-human dose estimation. J. Pharmaceut. Sci. 98, 3850–3861. doi: 10.1002/jps.21682

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, R. (1992). Optimal control drug scheduling of cancer chemotherapy. Automatica 28, 1113–1123. doi: 10.1016/0005-1098(92)90054-J

CrossRef Full Text | Google Scholar

Mei, Q., Xu, F., Chen, H., Li, Z., and Hu, Y. (2016). “Fast model predictive control based on multiscale system theory,” in 2016 12th World Congress on Intelligent Control and Automation (WCICA) (Guilin), 2517–2522. doi: 10.1109/WCICA.2016.7578600

CrossRef Full Text | Google Scholar

Michor, F., and Beal, K. (2015). Improving cancer treatment via mathematical modeling: surmounting the challenges is worth the effort. Cell 163, 1059–1063. doi: 10.1016/j.cell.2015.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, H. (2018). How to mathematically optimize drug regimens using optimal control. J. Pharmacokinet. Pharmacodyn. 45, 127–137. doi: 10.1007/s10928-018-9568-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Nisio, M. (2015). Stochastic Control Theory. Part of the book Series Probability Theory and Stochastic Modelling. Vol. 72. Tokyo: Springer. doi: 10.1007/978-4-431-55123-2

CrossRef Full Text

Pillis, L. G. D., Fister, K. R., Gu, W., Head, T., Maples, K., Neal, T., et al. (2008). Optimal control of mixed immunotherapy and chemotherapy of tumors. J. Biol. Syst. 16, 51–80. doi: 10.1142/S0218339008002435

CrossRef Full Text | Google Scholar

Sargent, R. (2000). Optimal control. J. Comput. Appl. Math. 124, 361–371. doi: 10.1016/S0377-0427(00)00418-0

CrossRef Full Text | Google Scholar

Schättler, H., and Ledzewicz, U. (2015a). “Cell cycle specific cancer chemotherapy for homogeneous tumors,” in Interdisciplinary Applied Mathematics (New York, NY: Springer), 41–114. doi: 10.1007/978-1-4939-2972-6_2

CrossRef Full Text | Google Scholar

Schättler, H., and Ledzewicz, U. (2015b). Optimal Control for Mathematical Models of Cancer Therapies. (New York, NY: Springer). doi: 10.1007/978-1-4939-2972-6_5

CrossRef Full Text | Google Scholar

Soetaert, K., and Petzoldt, T. (2010). Inverse modelling, sensitivity and Monte Carlo analysis in R Using Package FME. J. Stat. Softw. 33, 1–28. doi: 10.18637/jss.v033.i03

CrossRef Full Text | Google Scholar

Swan, G. W. (1980). Optimal control in some cancer chemotherapy problems. Int. J. Syst. Sci. 11, 223–237. doi: 10.1080/00207728008967009

CrossRef Full Text | Google Scholar

Swan, G. W. (1988). General applications of optimal control theory in cancer chemotherapy. Math. Med. Biol. 5, 303–316. doi: 10.1093/imammb/5.4.303

PubMed Abstract | CrossRef Full Text | Google Scholar

Swan, G. W. (1990). Role of optimal control theory in cancer chemotherapy. Math. Biosci. 101, 237–284. doi: 10.1016/0025-5564(90)90021-P

PubMed Abstract | CrossRef Full Text | Google Scholar

Sweilam, N. H., Al-Mekhlafi, S. M., Assiri, T., and Atangana, A. (2020). Optimal control for cancer treatment mathematical model using Atangana-Baleanu-Caputo fractional derivative. Adv. Diff. Equat. 334:2020. doi: 10.1186/s13662-020-02793-9

CrossRef Full Text | Google Scholar

Świerniak, A. (1995). Cell cycle as an object of control. J. Biol. Syst. 3, 41–54. doi: 10.1142/S0218339095000058

CrossRef Full Text | Google Scholar

Świerniak, A., Kimmel, M., Smieja, J., Puszynski, K., and Psiuk-Maksymowicz, K. (2016). “Cell cycle as an object of control,” in System Engineering Approach to Planning Anticancer Therapies (Cham: Springer International Publishing), 9–54. doi: 10.1007/978-3-319-28095-0_2

CrossRef Full Text | Google Scholar

Swierniak, A., Polanski, A., and Kimmel, M. (1996). Optimal control problems arising in cell-cycle-specific cancer chemotherapy. Cell Proliferat. 29, 117–139. doi: 10.1046/j.1365-2184.1996.00995.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, P. J., Olufsen, M., Sepulchre, R., Iglesias, P. A., Ijspeert, A., and Srinivasan, M. (2019). Control theory in biology and medicine. Biol. Cybernet. 113, 1–6. doi: 10.1007/s00422-018-00791-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., and Schättler, H. (2016). Optimal control of a mathematical model for cancer chemotherapy under tumor heterogeneity. Math. Biosci. Eng. 13, 1223–1240. doi: 10.3934/mbe.2016040

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, L., Li, M., Wang, J.-X., and Wu, F.-X. (2019). Controllability and its applications to biological networks. J. Comput. Sci. Technol. 34, 16–34. doi: 10.1007/s11390-019-1896-x

CrossRef Full Text | Google Scholar

Yankeelov, T. E., Atuegwu, N., Hormuth, D., Weis, J. A., Barnes, S. L., Miga, M. I., et al. (2013). Clinically relevant modeling of tumor growth and treatment response. Sci. Transl. Med. 5, 187ps9–187ps9. doi: 10.1126/scitranslmed.3005686

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Cunningham, J. J., Brown, J. S., and Gatenby, R. A. (2017a). Integrating evolutionary dynamics into treatment of metastatic castrate-resistant prostate cancer. Nat. Commun. 8:1816. doi: 10.1038/s41467-017-01968-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Wang, Z., Sagotsky, J. A., and Deisboeck, T. S. (2008). Multiscale agent-based cancer modeling. J. Math. Biol. 58, 545–559. doi: 10.1007/s00285-008-0211-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, W., Xu, H., Wang, H., and Lin, Z. (2017b). Stochastic systems and control: theory and applications. Math. Probl. Eng. 2017, 1–4. doi: 10.1155/2017/4063015

CrossRef Full Text | Google Scholar