Optimal chemotherapy counteracts cancer adaptive resistance in a cell-based, spatially-extended, evolutionary model

Orlando et al [10] proposed a deterministic PB model to describe cancer growth and adaptive response to chemotherapy. The main equation is a simple logistic growth:

$$\begin{align}\hfill \frac{\mathrm{d}}{\mathrm{d}t}N(t)& =N(t)\left[r\left( 1 - \frac{N(t)}{K} \right) - {d}_{1}{y}_{1}(t) - {d}_{2}{y}_{2}(t)\right.\hfill \\ \hfill & \quad \left.-\beta {y}_{1}(t){y}_{2}(t)\right],\hfill \end{align} \tag{ 1 }$$

in which N(t) is the cancer size at time t (proportional to the number of cells) and y₁(t) and y₂(t) are the concentrations of the two administered drugs in the tumor body. In the absence of therapy, the tumor grows to a maximum size K (the so-called carrying capacity of the logistic growth), with exponential rate r at low size. The two drugs induce pharmacologic mortality of the tumor cell at per-capita rates d₁ and d₂ per-unit of drug concentration, respectively. Their joint action can cause an extra mortality (synergistic drugs, β > 0) or reduce their effectiveness (antagonistic drugs, β < 0); however, in the present work we only consider the case of non-interacting drugs (β = 0).

Since there is no spatial structure, the drug concentrations are assumed uniform throughout the tumor. Their clearance is modeled with first order pharmacokinetics

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{y}_{1}(t)={w}_{i}(t)-{z}_{i}{y}_{i}(t),\quad i=1,2,\end{equation} \tag{ 2 }$$

where z₁ and z₂ are per-unit drug clearance rates and w₁(t) and w₂(t) are the rates of drugs administration, i.e., the control variables at the disposal of the oncologist (drug mass per unit of time, assuming a unit dilution volume).

Cancer adaptive response to chemotherapy is described by a two-dimensional cell phenotype: one component, u ≥ 0, measuring the cell overall investment in drug resistance and one, u₁ ∈ [0, 1], controlling the allocation of the acquired resistance to drugs 1 and 2. The larger is u, the larger is the fraction of the cell energy spent in the development and functioning of drug resistant mechanisms (u = 0 and u → ∞ corresponding to 0% and 100% of the energy, respectively); the smaller is consequently the cell energy spent in growth and duplication. The phenotype u therefore characterizes an energy allocation trade-off for the tumor between drug resistance and growth (see equation (3) below). Drug resistance can be fully specialized to one drug (u₁ = 1 to drug 1; u₁ = 0 to drug 2) or show any intermediate level of generalization toward both drugs. The resulting resistances to drugs 1 and 2 are defined by the products uu₁ and uu₂, where the allocation to drug 2, u₂, is bound to u₁ by one of the following drug resistance allocation trade-offs:

• (a)  
  Neutral, implemented by the linear bound u₂ = 1 − u₁. It is the case of noninterfering drug resistances, in which the total resistance u(u₁ + u₂) equals the energy investment u.
• (b)  
  Pro-specializing, implemented by the concave bound ${u}_{2}=1-\sqrt{1-{({u}_{1}-1)}^{2}}$. It is the case of negatively interfering drug resistances, in which the total resistance is lower than the energy investment u (lowest if the cell maximally generalizes on both drugs at ${u}_{1}={u}_{2}=1-\sqrt{2}/2$, see figure 1) and equals u only if the cell specializes on a single drug.
• (c)  
  Pro-generalizing, implemented by the convex bound ${u}_{2}=\sqrt{1-{u}_{1}^{2}}$. It is the case of positively interfering drug resistance, in which the total resistance is higher than the energy investment u (highest if the cell maximally generalizes on both drugs at ${u}_{1}={u}_{2}=\sqrt{2}/2$, see figure 1) and equals u only if the cell specializes on a single drug.

The trade-off names refer to the fact that the interference between the allocated drug resistances favors, in terms of the total acquired resistance, cells that are specialized on a single drug or generalist. Since there is no individual structure, (u, u₁) represents the mean phenotype of the population.

The tumor maximum size K and the mortality rates d₁ and d₂ are controlled by the phenotype in the following way. The maximum size is assumed to drop with the energy investment in drug resistance as a half-Gaussian bell

$$\begin{equation}K(u)={K}_{\mathrm{max}}\enspace \mathrm{exp}\left( - \frac{{u}^{2}}{2{\sigma }_{k}^{2}}\right),\end{equation} \tag{ 3 }$$

while the mortality rates decrease along with the specific drug resistance

$$\begin{equation}{d}_{i}(u,{u}_{i})=\frac{1}{{k}_{i}+u{u}_{i}},\quad i=1,2,\end{equation} \tag{ 4 }$$

where k_i is the baseline level of resistance to drug i. These phenotypic-dependences define the energy allocation trade-off for the tumor between drug resistance and growth. The cell energy allocated to the biological processes conferring drug resistance is invested at the expense of energy otherwise allocated to cell growth and duplication. This reflects the real behaviour of many cancer cells that under chemotherapy slow down their cell cycle and develop drug resistance [45, 46].

The phenotype is assumed to be fully genetically determined (no epigenetic nor plastic responses are considered), with the two components controlled by non-interacting sets of genes. This allows to consider mutations acting independently on u and u₁. The fitness of a mutation in evolutionary theory is the exponential rate of growth of the mutant population during its invasion transient [47, 48]. Denoting the mutant phenotype by (v, v₁), the fitness is given by the per-capita growth rate of the mutant (v, v₁)-population at low (technically zero) size, while invading the resident (u, u₁)-population at size N:

$$\begin{align}\hfill G(v,{v}_{1},N,{y}_{1},{y}_{2})& =r\left( 1 - \frac{N}{K(v)} \right) - {d}_{1}(v,{v}_{1}){y}_{1}\hfill \\ \hfill & \quad -{d}_{2}(v,{v}_{2}){y}_{2} - \beta {y}_{1}{y}_{2}.\hfill \end{align} \tag{ 5 }$$

When the resident population is at equilibrium, its growth rate is vanishing, i.e., G(u, u₁, N, y₁, y₂) = 0, so that the sign of the fitness, positive/negative, discriminates invading mutations from disadvantageous ones. For example, if treatment is interrupted after the development of some level of drug resistance u in the resident population, mutant cells with lower drug resistance have positive fitness (because of the higher maximum size K(v) > K(u), see equation (3), at y₁ = y₂ = 0), so that natural selection drives evolution toward the loss of the acquired drug resistance. In a deterministic setting, an invading mutation replaces the former resident type [49], so that the fitness gradient with respect to the mutant phenotype (v, v₁) gives the direction of natural selection. Taking into account the stochastic nature of mutations, but looking at the limit in which they are extremely rare and each entailing an extremely small phenotypic change, the expected evolutionary dynamics of the mean phenotype is described by the following two equations:

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}u(t)=s {\left.\frac{\partial }{\partial v}G(v,{u}_{1}(t),N(t),{y}_{1}(t),{y}_{2}(t))\right\vert }_{v=u(t)} ,\end{equation} \tag{ 6a }$$

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{u}_{1}(t)=s {\left.\frac{\partial }{\partial {v}_{1}}G(u(t),{v}_{1},N(t),{y}_{1}(t),{y}_{2}(t))\right\vert }_{{v}_{1}={u}_{1}(t)}\hspace{-14.22636pt}\quad \end{equation} \tag{ 6b }$$

(known as the canonical equations of adaptive dynamics [47, 48]), in which s is a time-scale parameter that goes to zero in the above limit. To describe real situations, s is taken proportional to the frequency of mutations and to the variance of the entailed phenotypic changes. The larger s, the lower is however the accuracy of equations (1)–(6) in describing the real demographic and evolutionary dynamics. When the frequency of mutations is of the order of one every few cell duplications, the paradigm underneath equation (6) of each invading mutant substituting the former resident type gets lost, as new mutations appears while the substitution transient is not yet over. In these situations a stochastic cell-based approach is more accurate.

Orlando et al's model of cancer growth and evolution under chemotherapy is composed of equations (1)–(6). The optimal control problem formulated in [10] finds the treatment (w₁(t), w₂(t)) in a given timespan [0, t_f] that minimizes the final tumor size N(t_f), starting from a fully grown untreated tumor

$$\begin{equation}N(0)={K}_{\mathrm{max}},\hspace{5.0pt}{y}_{1}(0)={y}_{2}(0)=0,\hspace{5.0pt}u(0)=0,\hspace{5.0pt}{u}_{1}(0)={u}_{2}(0),\end{equation} \tag{ 7 }$$

subject to the dynamics (1)–(6) and to the toxicological constraints

$$\begin{equation}{w}_{1}\leqslant 10,\hspace{5.0pt}{w}_{2}\leqslant 10,\hspace{5.0pt}{y}_{1}+{y}_{2}\leqslant 10.\end{equation} \tag{ 8 }$$

The optimal solutions under the three drug resistance allocation trade-offs (a)–(c) have a common beginning, during which both drugs are delivered at the maximum admissible rate to increase drug concentrations up to hit the toxicological constraints. Then, the solutions corresponding to the neutral and pro-generalizing trade-offs are constant and keep equal maximal drug concentrations for the rest of the treatment. Differently, under the pro-specializing trade-off, the optimal therapy is dynamic and alternates the two drugs. The therapy exploits the pro-specializing trade-off to select for tumor resistance for one drug and then hit the tumor with the other drug toward the end of the treatment.

Unlike PB models, which are defined in terms of averaged per-capita demographic rates, cell-based (CB) models are IB. The model assumptions are made at the level of the single individual, the cell, about its characteristic parameters and about the procedural rules of interaction with the neighboring cells and with the surrounding environment. The population-level behaviour emerges bottom-up, as a result of the model simulation [50].

Though, in principle, any physiological and genetic detail can be included in a CB setting, we limit our description to a 2D spatial extension of Orlando et al's [10] PB model (described in the previous section). As discussed by Orlando and coauthors, the challenge for CB models of cancer growth, evolution, and treatment is computational, especially in the optimization of treatment. Our work is therefore to be considered a step forward to make the challenge feasible.

We represent space by a square grid of 100 × 100 nodes, with connection degree four, i.e., each node communicates with its four neighbors. Each node is a cell, which can be healthy or tumoral, though only the dynamics of the latter will be considered in line with Orlando et al's [10] model.

Each node is characterized by a state profile that includes the following information:

• The type of cell, healthy or tumoral.
• For tumor cells, the phenotypic strategy (u, u₁): u ⩾ 0 is the overall energy investment in drug resistance; u₁ ∈ [0, 1] is the resistance allocation to drug 1.
• The position of the cell in the grid, row r, top-to-bottom, and column c, left-to-right.
• The drug concentrations y₁ and y₂ at the location (r, c).

Time proceeds in discrete steps, each corresponding to an hour of real life. From a time step t to the next, the state of all nodes is updated according to the procedural rules described in the remaining part of this section. Some of the rules are stochastic, to better describe the intrinsic uncertainties of cellular processes, such as duplication, mutation, and death, and the regulating probabilities are scaled to the assumed time step. (An alternative, that we do not follow, is to stochastically draw the timing of the various events according to their rates of occurrence [51].)

2.2.1. Duplication

The cell cycle—cell growth up to cell division—of a cancer cell not investing energy in drug resistance (phenotype u = 0) is assumed to correspond to the time step. Although cancer cells often show cycles longer than 1 h [52], we consider a fast-growing tumor.

The energy allocation trade-off between growth and drug resistance is introduced by reducing the probability to complete the cell cycle in one time step along with the energy investment in resistance (phenotype u) [45, 46]. Specifically, we set the probability b(u) that the cell cycle is completed at the end of each time step (independently from one step to the next), by means of the following exponential relation:

$$\begin{equation}b(u)=\mathrm{exp}(-{b}_{1}u),\end{equation} \tag{ 9 }$$

where parameter b₁ is calibrated (in section 2.3) to obtain an unviable tumor at a value of u at which the tumor maximum size K(u) of Orlando et al's model is negligible.

Note that the slope of function (9) at u = 0 is negative, in contrast to the zero slope of the maximum size K(u) in equation (3). This is a minor change, with respect to Orlando et al's model, that is however justified by energetic principles and will be quantified in section 2.3.

To introduce the inter-cellular competition for energy and space—described by the logistic growth in Orlando et al's model—the duplication of tumor cells (cell division) is assumed to occur at the completion of the cell cycle only in the presence of neighboring healthy cells. That is, we assume that the tumor grows by displacing healthy tissues.

Specifically, a neighboring location is randomly selected and duplication occurs only if the location is not occupied by another cancer cell. To minimize the spatial bias, cancer cells are considered for duplication in a new random order at each time step.

2.2.2. Mutation

Like Orlando et al's PB model, our CB extension describes the genetic of cancer at phenotypic level, and the two considered phenotypes are assumed to be fully genetically controlled by non-interacting sets of genes. Genetic mutations are therefore implemented as independent phenotypic variations at cell duplication, where only the daughter cell occupying the new location potentially mutates because of chromosomes copying errors.

The two phenotypes independently mutate with a probability μ. The extent of the mutation is drawn from a normal distribution with zero mean and standard deviation σ for the overall energy investment in drug resistance u and σ₁ for the resistance allocation u₁ to drug 1. If the resulting mutant phenotypes are unfeasible (u < 0 or u₁ ∉ [0, 1]) they are forced to their extremal values. Parameters μ, σ, and σ₁ are calibrated in section 2.3 to reflect the desired cancer evolutionary response.

Note that, because of the nonlinearity of the concave and convex bounds imposed by the pro-specializing and pro-generalizing allocation trade-offs (see cases 2 and 3 in section 2.1), mutations in u₁ are not equivalent to opposite mutations in u₂ (see figure 1). By implementing mutations on u₁, the result is that symmetric scenarios, such as the administration of only one drug starting from zero resistance, do not yield symmetric model simulations. For example, after that cancer cells develop some resistance to drug 1 under the pro-specializing trade-off, a successful (positive) mutation of u₁ (mutation at point A in the figure) corresponds to a smaller step (AB) on the concave trade-off curve, than the one produced by the opposite u₁-mutation while administrating drug 2 at the symmetric point of the curve (step A'B').

Zoom In Zoom Out Reset image size

To make mutations in u₁ equivalent to opposite mutations in u₂, we redefine the phenotype u₁ as the fraction of the overall energy investment u in drug resistance allocated to drug 1, u₂ = 1 − u₁ being the fraction allocated to drug 2 (so that point (u₁, u₂) lies on the linear curve of figure 1, independently of the type of resistance allocation trade-off). Because of the interference between the two allocations, we assume that the amounts of energy actually available for effective drug resistant mechanisms against drugs 1 and 2 are given by the products ur₁ and ur₂, respectively, where r₁ and r₂ are the coordinates of the point of the trade-off curve sharing the same ratio u₂/u₁. That is, we look at the intersection of the trade-off curve with the ray passing through (u₁, u₂) (see figure 1). In formulas, we have

$$\begin{equation}{r}_{1}({u}_{1})=\frac{{u}_{1}}{1+\sqrt{2{u}_{1} (1-{u}_{1})}},\hspace{5.0pt}{r}_{2}({u}_{1})=\frac{1-{u}_{1}}{1+\sqrt{2{u}_{1} (1-{u}_{1})}}\end{equation} \tag{ 10a }$$

for the pro-specializing trade-off and

$$\begin{equation}{r}_{1}({u}_{1})=\frac{{u}_{1}}{\sqrt{1-2{u}_{1} (1-{u}_{1})}},\hspace{5.0pt}{r}_{2}({u}_{1})=\frac{1-{u}_{1}}{\sqrt{1-2{u}_{1} (1-{u}_{1})}}\end{equation} \tag{ 10b }$$

for the pro-generalizing one. Consequently, we redefine the pharmacologic mortality rates in (4) as

$$\begin{equation}{d}_{i}(u,{u}_{1})=\frac{1}{{k}_{i}+u\enspace {r}_{i}({u}_{1})},\quad i=1,2.\end{equation} \tag{ 11 }$$

2.2.3. Mortality

We assume that the energy allocation trade-off also affects the cell death process. We therefore add to the pharmacologic mortality rates in (11) a natural mortality (not drug-induced)

$$\begin{equation}m(u)={m}_{0}+{m}_{1}u,\end{equation} \tag{ 12 }$$

that increases with the energy investment in resistance. Parameters b₁ (in equation (9)) and m₁ (here in equation (12)) measure the trade-off strength and are calibrated in section 2.3. In particular, m₁ takes into account that the more energy is invested in drug resistance, the less remain for mechanisms resisting cellular death and avoiding immune destruction.

The total (per-cell) mortality rate during time step t is hence

$$\begin{equation}d(u,{u}_{1},t)=m(u)+{d}_{1}(u,{u}_{1}){y}_{1}(t)+{d}_{2}(u,{u}_{1}){y}_{2}(t),\end{equation} \tag{ 13 }$$

where y_i (t) is the concentration of drug i at the cell location. Because drug diffusion is simulated on a finer time scale with N_y sub-steps (see next section), the probability of cell death is evaluated at each sub-step k based on the drug concentrations y_i (t + k/N_y ), k = 1, ..., N_y . Assuming an independent Poisson process for each cell within each time sub-step, each cancer cell independently dies in sub-step (t, k) with probability

$$\begin{equation}1-\mathrm{exp}(-d(u,{u}_{1},t + k/{N}_{y})).\end{equation} \tag{ 14 }$$

2.2.4. Drugs diffusion

Being endowed with a 2D spatial structure, our model can describe the diffusion dynamics of the two considered drugs through the tissues represented by the cell-grid. We imagine that a capillary vessel is located in the grid top side, bringing oxygen and nutrients as well as the drugs. We assume that the drugs percolate top-to-bottom, following the pressure-gradient from the top artery region towards the bottom venous one, with no distinction between healthy and tumor cells.

Because the time scale of fluid diffusion through our body tissues is typically shorter than 1 h, we divide the time step into N_y equal sub-steps, t + k/N_y , k = 1, ..., N_y , and we update the drug concentrations y_i (t + k/N_y ) accordingly, i = 1, 2. We also assume that the oncologist can directly control the drug concentrations w_i (t) at the vessel in each time step (or, more typically, on a daily or weekly base), thus neglecting the vascular transport dynamics. Thanks to this simplification, the limitations imposed by the chemotherapeutic protocol on drug supply (see the constraints on w₁ and w₂ in (8)) imply that the constraints on the drug concentrations remain satisfied throughout the treatment.

The resulting diffusion dynamics, from sub-step (t, k) to sub-step (t, k + 1) is defined by the following equations:

$$\begin{equation}{y}_{i,1,c}\left.t +(k+1)/{N}_{y}\right.={w}_{i}(t)\end{equation} \tag{ 15a }$$

$$\begin{align}\hfill & {y}_{i,r,c}\left.t + (k+1)/{N}_{y}\right.=(1-z){y}_{i,r,c}\left.t +k/{N}_{y}\right.\hfill \\ \hfill & +z{y}_{i,r-1,c}\left.t +k/{N}_{y}\right.\hspace{17.07164pt}\hfill \end{align} \tag{ 15b }$$

r > 1, c ⩾ 1, where z is the fraction of the vascularizing blood that goes, every 60/N_y minutes, to the neighboring cell one row below. Left and right diffusion compensate, because we assume that the controlled concentrations w_i (t) are homogeneously set in the grid top row (see equation (15a) for any column c). Parameter z is calibrated in section 2.3 to match the drug clearance rates of Orlando et al's model.

In the following subsections we present a series of simulations that we run to calibrate the parameters of our CB model to reflect the dynamical characteristics of Orlando et al's PB model. The involved parameters are listed in the subsection title. All model parameters and their values are summarized in table 2. The model has been implemented and numerically analyzed in the MATLAB environment.

2.3.1. The demographic time scale (m₀, k₁, k₂)

2.3.2. The energy allocation trade-off (b₁, m₁)

To calibrate the energy allocation trade-off, we simulate the untreated growth of a non-adaptive tumor starting from four cells in the middle of our grid. We artificially test the unrealistic scenario in which a drug-free tumor is characterized by resistance u > 0 and we stop the simulation when the tumor is either extinct or reaches a stationary size. Extinction is unlikely at low u, whereas it certainly occurs if u is so high to excessively limit growth. Convergence to a stationary solution is obtained when the number of cancer cells, averaged over the last week, stays within a 1% variation for a week. The transition, for increasing u, from a positive tumoral regime to extinction corresponds to the exchange of stability between the two solutions (a bifurcation in the systems theory terminology).

The energy allocation trade-off between growth and drug resistance is controlled by parameters b₁ and m₁ in equations (9) and (12), respectively. We assume that the growth limitation equally splits between reduced cell duplication and increased mortality. We therefore set b₁ = m₁ and we search for the value that is consistent with the tumor maximum size of Orlando et al's model. Specifically, though the support of the Gaussian bell used in equation (3) is technically unbounded, most of the density occurs within three standard deviations, so that we look for the value of b₁ yielding the above described extinction transition at about u = 90 ≈ 3σ_K (σ_K = 30 in Orlando et al's simulations). The result is obtained by averaging the regime tumor size over 10 simulations for each of a sequence of values of u. Figure 2 shows the result for b₁ = 1/185.

Zoom In Zoom Out Reset image size

2.3.3. Mutations (μ, σ, σ₁)

2.3.4. Drugs diffusion (N_y ,z)

We test the calibration of our model by comparing the regime level of the energy investment in drug resistance u developed with constant maximum dosage of both drugs (w₁(t) = w₂(t) = 5 at all time steps) under the neutral resistance allocation trade-off (case 1 in section 2.1) against the level reached in Orlando et al's model.

With the parameter calibration of the previous sections (summarized in table 2, except for parameters b₁ and m₁ that were set at 1/185 in section 2.3.2), our model yields a regime value u ≈ 19, whereas u converges to about 22.5 in Orlando et al's model. To improve the matching, we slightly modify the parameters b₁ and m₁—respectively regulating the energy allocation trade-off between drug resistance and duplication and mortality—to the value reported in table 2, for which the tumor extinction shown in figure 2 occurs at about u = 100.

The result is shown in figure 3, averaged over 10 simulations starting from the same 10 fully-grown untreated tumors (initial phenotypes u = 0, u₁ = 0.5) used in calibration (see sections 2.3.2 and 3). Convergence is obtained in about five months, consistently with the calibration of the mutation process (section 2.3.3). The solid black line is the average u (averaged over the grids of the 10 simulations at each time step); the shaded area indicates the standard deviation (computed over the 10 grids); the dashed line is the solution of Orlando et al's equations (1)–(6). Note that the time scale of convergence in Orlando et al's model (reported on the top time axis) is much longer than ours (one time unit of equations (1)–(6) corresponds to 100 time steps, about 4 days, in our model, so that the top time axis is rescaled to simultaneously show convergence). This is due to the missing separation, in Orlando et al's model, between the demographic and evolutionary time-scales (as already commented in section 2.3.3). However, by suitably rescaling time, the solution of Orlando et al's model (dashed line) is essentially within a standard deviation from the one of our calibrated model, indicating the agreement between the two models under a static chemotherapeutic protocol.

Zoom In Zoom Out Reset image size

To look for the optimal chemotherapy schedule, we assume that the drug supplies w₁(t) and w₂(t) are changed only weekly over a five-month period. Moreover, with no other toxicological constraint, the optimal solution necessarily lies on the bound w₁(t) + w₂(t) = 10 of maximum admissible dosage (equation (8)). The control variables are hence given by the 5 × 4 = 20 weekly dosages of drug 1, w_1,k (w_2,k = 10 − w_1,k ), k = 1, ..., 20, that we pack in vector W. The control target is a functional of the cancer grow and evolutionary dynamics along the treatment, that we average over 10 realizations starting from the same fully-grown untreated tumors (initial phenotypes u = 0, u₁ = 0.5) used in calibration (see sections 2.3.2 and 3). For example, to implement Orlando et al's target on the final tumor size, we look at the minimum number of cancer cell in the last week of treatment. Notwithstanding the intrinsic stochasticity, that we smooth out by averaging over 10 realizations, we qualitatively consider the target as a function F(W) of the control variables.

To minimize the target, we implement an adaptive version of the classical gradient descent method. Starting from a given control W⁽⁰⁾ yielding target F(W⁽⁰⁾), the idea is to produce a sequence controls W⁽ⁱ⁾, i ≥ 1, along which the target is reduced to a (possibly local) minimum. To find W⁽ⁱ⁾ from W⁽ⁱ⁻¹⁾, we first perturb the control variables in W⁽ⁱ⁻¹⁾ one at a time by a fixed amount Δw (we use Δw = 0.5, i.e., 5% of the variable range; we use Δw = −0.5 if, when moving variable k, ${W}_{k}^{(i-1)}+{\Delta}w > 10$) and compute the kth component of the target gradient as

$$\begin{equation}\nabla {F}_{k}({W}^{(i-1)})=\frac{F({W}^{(i-1)}+{E}^{(k)}{\Delta}w)-F({W}^{(i-1)})}{{\Delta}w},\end{equation} \tag{ 16 }$$

where E^(k) is the kth vector of the canonical basis of the strategies space (${E}_{k}^{(k)}=1$; all other components being zero).

We then move the control variables in the direction

$$\begin{equation}{h}^{(i)}=-\nabla F({W}^{(i-1)})/{\Vert}\nabla F({W}^{(i-1)}){\Vert},\end{equation} \tag{ 17 }$$

opposite to the gradient, by an amount r_i that is the step size at the step i of the algorithm (variables that go out of range, w_1,k < 0 or w_1,k > 10 are truncated at 0 and 10, respectively). We first try the last step size, i.e., r_i = r_i−1 (r₀ = 10). If we find a better control, i.e. F(W⁽ⁱ⁻¹⁾ + r_i h⁽ⁱ⁾) < F(W⁽ⁱ⁻¹⁾), the step is accepted, we set W⁽ⁱ⁾ = W⁽ⁱ⁻¹⁾ + r_i h⁽ⁱ⁾, and r_i is increased by 50% for the next step; otherwise the step size r_i is halved and the target F(W⁽ⁱ⁻¹⁾ + r_i h⁽ⁱ⁾) recomputed, until a better solution is found or the step size goes below a minimum threshold that we set at 10⁻⁵. If the minimum step size is reached, the gradient is recomputed at W⁽ⁱ⁻¹⁾ with halved perturbation Δw = 0.25 and the search is repeated starting with step size r_i−1. If the minimum step size is again reached, the same is repeated with doubled Δw = 1, before the algorithm stops proposing W⁽ⁱ⁻¹⁾ as best solution.

To cope with the possible local minima, we always repeat the algorithm for each of three starting constant solutions W⁽⁰⁾, the two single-drug protocols (w_1,k = 10 and w_1,k = 0 for all k) and the one that equally administers both drugs (w_1,k = 5 for all k).

In line with Orlando and colleagues [10], we look for the five-month chemotherapy schedule that minimizes the cancer cells number at some point during the final week of treatment. We note that with the calibration of the model parameters performed in section 2.3 and the dosage constraint w₁(t) + w₂(t) ≤ 10, the therapy is not able to eradicate the tumor, which persists even with the constant dosages w₁(t) = w₂(t) = 10 violating the constraint. Minimizing the tumor size at the end of a finite time horizon, could therefore be the optimal strategy before a removal surgery.

We find the same optimal weekly schedule, independently of the drug resistance allocation trade-off (neutral, pro-specializing, or pro-generalizing). The optimal solution exploits the evolutionary double bind mechanism [10, 41], administering only one drug but the final week, in which the other drug is used. This solution has a clear meaning. Giving only one drug during a time interval, the oncologist selects the cancer cells for a specialized resistance. This is true under any of the three resistance allocation trade-offs, because the pharmacologic mortality due to a single drug decreases by specializing the resistance along all the trade-off curves (see figure 1). At the beginning of the final week, the oncologist switches to the other drug and, after a washout period (which is rather fast, few hours, in our setting), cancer cells find themselves maximally sensitive to the new drug, while wasting the actual energy investment in drug resistance at the expense of reduced growth and survival rates. This is the worst scenario for the tumor, so that many cancer cells get killed and the tumor size rapidly drops.

Because of the assumed rapid evolution of cancer cells, especially the rapid reallocation of drug-resistance (phenotype u₁) resulting from the model calibration, once the new demographic regime is reached, the remaining cell start to revert the allocation, thus creating a more favorable environment for tumor regrowth. This is clearly visible in figure 4 (though the target in section 3.2 is different), when the drug concentrations are suddenly changed (red lines in the top panels). See, e.g., the middle panels under the pro-specializing trade-off: after a quick drop in about two days (blue in the top panel), the tumor size increases as long as the cancer cells reallocate their resistance to the new drug (blue in the bottom panel).

Zoom In Zoom Out Reset image size

The minimum tumor size is reached at a specific time during the last week, to which the oncologist should carefully pay attention—possibly supported by tumour mass scanning—to optimally plan the removal surgery. The effective minimum, when it is achieved, and the obtained improvement with respect to the best and worst static protocols depend on the drug resistance allocation trade-off. Numerical results are reported in table 3.

We now imagine that the oncologist wants to keep the tumor mass under control along the treatment, e.g., with monthly monitoring mid-points. To this purpose, we minimize the averaged number of cancer cells during the last weeks of the five months of the protocol.

We discover that the optimal schedule is dynamic for all the investigated drug resistance allocation trade-offs (neutral, pro-specializing, or pro-generalizing). The three optimal solutions (numerically obtained and averaged over 10 realizations starting from the fully-grown untreated tumors used in calibration in sections 2.3.2 and 3) are displayed in figure 4 and share a common beginning: after inducing cancer cells to specialize by administering only one drug for three weeks (drug 1 in the figure), the other drug (drug 2), to which the tumor is at that point highly sensitive, is used during the first monitoring week to kill an effective number of cancer cells. Also note the significant initial drop of about 70% of the tumor size, from the untreated regime of about 8600 cells, due to the pharmacologic mortality. This drop is very fast (one-two days) also because of the quick pharmacokinetics (clearance in few hours).

Under the pro-specializing trade-off (figure 4 middle panels), cancer cells are not able to revert the allocation of resistance from drug 1 to drug 2 during the first monitoring week. Not even the generalized resistance at u₁ = 0.5 is reached. Thus, in the first week of the incoming month, the drug used in the monitoring week is still the most effective. Moreover, while steering the tumor to exchange its specialized resistance from drug 1 to drug 2, the two actual allocations r₁(t) and r₂(t) must follow the concave trade-off curve of figure 1, along which the total resistance is lower than a fully specialized one and minimal at the point of generalized resistance u₁ = 0.5. For both these reasons, after the monitoring week, the optimal schedule continues to administer the same drug also in the incoming month, up to the next monitoring week. At the beginning of the next monitoring week, the situation is the same as the previous month but with the allocation of resistance and drug concentrations exchanged. The oncologist can again exploit the evolutionary double bind mechanism [10, 41], hitting the cancer with the effective drug, and so on for the rest of the protocol.

Under the pro-generalizing trade-off, after the first three common weeks, to optimally exploit the evolutionary trade-off, the cancer cells must be hit with a mix of the two drugs (see figure 4, right, where w₁(t) is nonzero in the fourth week). This is essentially due to the different speeds of evolution of the resistance allocation phenotype u₁ toward specialized and generalized resistance. The first is smaller than the second, because of the shape of the convex trade-off curve (see figure 1, convex curve, where mutations close to specialized resistance to drug i yield minor changes in the actual allocation r_i (t), whereas large changes occur in both r₁(t) and r₂(t) close to generalized resistance); the situation being reversed under the pro-specializing trade-off. Thus, at the beginning of the first monitoring week, the tumor is not fully specialized to drug 1, and this makes it easier to reach and pass generalized resistance u₁ = 0.5 during the week. Moreover, evolution toward generalized resistance is sped up by the evolutionary trade-off, so that using only drug 2 would lead to a strong specialization to drug 2 and to a tumor regrowth that would increase our target function. The optimal drug mix turns out to allow the tumor to pass generalized resistance, so that it is convenient to switch back to drug 1 in the incoming second month. In the first week the tumor is indeed more sensitive to drug 1, though eventually the target requires to specialize the tumor to drug 1, to be again severely attacked by drug 2 during the second monitoring week. The control strategy is repeated in the following months, though the drug mix used in the monitoring weeks changes according to the evolution of the energy investment in resistance u (black line in figure 4, right-bottom), which modulates the speed of evolution of the resistance allocation u₁.

Under the neutral trade-off (figure 4, left column), the optimal schedule found by our algorithm follows the same logic of the solution found for the pro-generalizing trade-off (right column). That is, the logic of using one drug to induce tumor specialization during the three consecutive unmonitored weeks and of hitting the tumor with the other drug during the monitoring weeks. The difference is that only one drug is used in the monitoring weeks (drug 2 in the figure), instead of a mix as it occurs under the pro-generalizing trade-off. The reason is again the relative speeds of evolution of the resistance allocation phenotype u₁ toward specialized and generalized resistance. Indeed, under the neutral trade-off and compared to the pro-generalizing one, cancer cells evolve faster/slower toward specialized/generalized resistance, so that the tumor turns out to be more specialized at the beginning of a monitoring week and can be hit with the effective drug at maximum dosage (drug 2) without inducing excessive tumor specialization and regrowth. During the monitoring week, the allocation phenotype u₁ crosses generalized resistance u₁ = 0.5, so that the best strategy is to switch back to the first administered drug (drug 1) to induce tumor specialization in the new incoming month.

We note that, under the neutral trade-off, the optimal schedule we found for the pro-specializing trade-off is equally optimal (differences are within the variability due to the model stochasticity). Indeed, the neutral one is an intermediate case between the pro-specializing and pro-generalizing trade-offs. Table 4 reports the obtained target values (average tumor size in the monitoring weeks; results averaged over 10 realizations). Because the identification of the type of drug resistance allocation trade-off could be a critical issue for the oncologist, we apply the three identified best protocols to all cases. Interestingly, the target loss of using the best schedule for the wrong trade-off is limited and, in any case, smaller than the loss of using the best static protocol.

Also note the significantly increased computational time with respect to table 3. The reason is that we incrementally optimized one more month of treatment at a time, starting from the obtained solution for the already optimized months and from the three constant solutions (the two single-drug and the one that equally administers both drugs) for the new month.

If the oncologist's goal is to reduce the tumor mass as much as possible in a prefixed time horizon, for example to schedule a removal surgery, the optimal protocol is the same independently of the drug resistance allocation trade-off. As exposed in section 3.1, the optimal protocol is to steer cancer cells to a specialized resistance by administering only one drug and then hit the tumor mass with the other drug, to which cancer cells are maximally sensitive. During this second phase, cancer cells suffer the highest mortality, but the surviving ones reallocate their resistance to the new drug and the mass eventually starts to regrowth. The tumor mass hence reaches a minimum to which the oncologist must pay careful attention.

Note that, when only one drug is administered, cancer cells specialize for that drug, even under the pro-generalizing resistance allocation trade-off. Moving along the trade-off curve (see figure 1) toward specialization to drug i, the resistance to drug i increases toward 1, so that the pharmacologic mortality induced by drug i monotonically decreases, i = 1 or 2. Once full specialization is reached, cancer cells are in the most favorable condition for the given drug dosage. Administering a single drug at constant maximum dosage is therefore the worst protocol among the ones using maximum admissible drug dosage.

Conversely, the best constant protocol is the one that equally administers both drugs at maximum total dosage. Under the neutral resistance allocation trade-off this is easy to see. Imagine a sequence of different constant protocols obtained from the worst one administering only drug 1 by replacing an increasing amount of drug 1 with the same amount of drug 2. The tumor cells optimal allocation of drug resistance then moves from fully specialized to generalist, while the resulting mortality rate increases along the sequence (see equation (13)). Because of the symmetry between drugs 1 and 2, the maximum mortality that can be achieved with a static protocol is hence the one obtained by administering equal constant amounts. Interestingly, the 'bang-bang' control—solutions in which only one drug is administered at a time—optimally exploits the evolutionary double bind mechanism [10, 41] and outcompetes the best static protocol by 40%–50% (detailed results in table 3).

Regarding the target with monthly monitoring mid-points (section 3.2), depending on the drug resistance allocation trade-off, we found two different logics of control. One makes mainly use of one drug, during the first three unmonitored weeks of each month to induce cells specialization, and administers the other drug in the monitoring weeks to severely reduce the tumor mass; the other is a bang-bang solution that alternates the two drugs, with the switch at the beginning of each monitoring week. The latter is optimal under the pro-specializing trade-off, whereas a version of the former that slightly mix the two drugs in the monitoring weeks is optimal under the pro-generalizing trade-off. For the neutral trade-off, the two bang-bang versions of the two solutions are equally optimal.

What makes one of the two logics optimal has to do with the relative speeds of evolution of the resistance allocation phenotype toward specialized and generalized resistance. Under the pro-specializing trade-off, cancer cells evolve faster/slower toward specialized/generalized resistance, whereas the situation is reversed under the pro-generalizing trade-off. In the first case, cancer cells are not able to specialize to the new drug during the monitoring week, because to do so they have to first evolve to generalized resistance. Consequently, they are still more sensitive to the new drug at the beginning of the new month. Moreover, by inducing full specialization to the new drug during the first three weeks of the new month, the optimal treatment also exploits the higher evolutionary cost, with respect to the pro-generalizing trade-off, of evolving along the concave trade-off curve (see figure 1).

Conversely, the tumor crosses generalized resistance and specializes to the new drug during the monitoring weeks under the pro-generalizing trade-off, also because the starting point at the beginning of the week is less specialized than under the pro-specializing trade-off, because of the slower evolution toward specialization. Actually, using only the new drug would induce a too strong specialization and a consequent tumor regrowth, so that the new drug is optimally mixed with the previous one used to specialize the tumor.

In clinical situations, the type of drug resistance allocation trade-off is most of the times unclear, if not unknown. On the one hand, our model suggests a way to estimate the trade-off curve. For example, a long-term experiment could monitor the reallocation of resistance from drug 1 to drug 2. First only drug 1 is administered to induce a fully specialized resistance, then only drug 2 is given to drive the reallocation. During this transient, the pharmacological mortality to drug 2 is monitored, and the one to drug 1 could be estimated on short experiments in which some cells are removed and treated with drug 1. From the pharmacological mortalities, a measure of the relative drug resistance could be inferred by a model like equation (4).

On the other hand, our results show that the bang-bang dynamic solutions identified by the two logics of control both achieve significantly better targets than the best static solution (see table 4, where n-best is the best protocol for the neutral trade-off, i.e., the bang-bang solution that makes mainly use of one drug, while ps-best is the best protocol for the pro-specializing trade-off, i.e., the bang-bang solution that alternates the two drugs). Under the uncertainty on the type of trade-off, the best protocol for the pro-specializing trade-off (ps-best in table 4) is preferable, because it is equally optimal for the neutral trade-off and quasi-optimal for the pro-generalizing one.

Confronting our results relative to the target at the end of treatment (section 3.1) against those obtained by Orlando and colleagues [10], we find a qualitative agreement only under the pro-specializing trade-off of drug resistance allocation. Both optimal solutions—our own and the one of Orlando et al's PB model—are dynamic and exploit the tumor evolutionary double bind. Besides the fact that the control variables—the drug dosages—are continuous in Orlando et al's model and discrete on a weekly base in our model, the solution logic is qualitatively the same: induce tumor specialization under the effect of a single drug and finally hit it with the other one.

Quantitatively speaking, however, there are significant differences. In the PB model, tumor cells are all identical, so they all homogeneously specialize under the effect of a single drug and they all suffer the highest mortality when the oncologist switches to the other drug. In the CB setting, the stochasticity of mutations, as well as that of demographic processes, always leaves some unspecialized cell that are immediately drug resistant when facing the new drug. In a way, the more realistic setting confers the tumor a more reactive adaptive resistance. This makes our solution sharper (bang-bang solutions) than the continuous oscillations obtained in the PB model. Moreover, our optimal strategy makes use of a single drug throughout the treatment, but the last period (the last week), to induce the highest specialization of the tumor cells. Conversely, the fact that full specialization is more easily reached in the PB setting makes possible to exploit the evolutionary double bind several times during treatment, though this option seems an artefact of the unrealistic setting. Without considering the toxicological or clinical constraints of administering one drug for a relatively long period, the best strategy to target the tumor size at the end of treatment is to exploit the evolutionary double bind trap only in the last phase of the cure.

Moreover, Orlando and colleagues focus on a very short treatment. Because they do not clearly set the time-scale separation between the demographic and evolutionary processes (see the comment in section 2.3.3 on the calibration of the mutation process), their unit of evolutionary time turns out to be equivalent to about 4 days in our simulations, so that the eight-unit span of their simulations must be confronted to a month in ours. Moreover, by setting the evolutionary rates (mutation probability and standard deviation of the phenotypic mutational step) to reach a regime of the average cells investment in drug resistance in about five months, we end up with a more aggressive (fast-adaptive) tumor than Orlando and colleagues. Indeed, the cells investment in drug resistance (phenotype u) reaches values of about 10 in one month in the optimal solutions in figure 4, whereas it nearly reaches 1 in Orlando et al's optimal solution under the pro-specializing trade-off (see the middle-right panel in figure 3 in [10]). With such a small resistance to allocate, the gain offered by the evolutionary trap is very limited. The target gain (tumor size at the end of treatment) with respect to the best static protocol is negligible (and not published in [10]), whereas the gain against the worst static protocol that administer only one drug is small (7.81%).

Under the neutral and pro-generalizing resistance allocation trade-offs, Orlando and colleagues find the best static protocol to be optimal, with limited gain against the worst protocol (4.92% and 14.6%, respectively). This is again due to the little resistance developed on the short cure horizon. Static and dynamic solution are practically equally effective in that context, having the cancer little resistance to allocate. Conversely, the explicit spatial and genetic heterogeneities of our model show that it is possible to take significant advantage of the neutral and pro-generalizing resistance allocation trade-offs, as well, using the same dynamical solution designed for the pro-specializing trade-off.

Our relatively simple test of applying optimal control techniques to a CB model of cancer growth and evolution can be extended in many directions. The challenge is to increase the model realism, in line with the recent advancements in cancer modeling [37, 38, 40] possibly tailoring to specific cancers and to personalized medicine, while keeping feasible the computational burden of control.

For example, cell metabolic processes and oxygen/nutrients transport could be explicitly modeled, together with mechanisms of angiogenesis, invasion of new tissues, and metastasis [36, 37]. Tumor geometry could be described with heterogeneous 3D structures, taking explicit account for the tumor micro-environment and for the location of the tumor in the human body [56].

As for the mechanisms of drug resistance, phenotipic plasticity could take into account the resistance developed by the cell in response to changes in its environment, such as the administration of a new drug. Especially the fast reallocation of specific resistance to drugs could be better described as a plastic [57, 58] or epigenetic [59, 60] response, rather than by permanent DNA modifications. Multi-drug resistance could also be added, letting mechanisms that are specific to a drug to partially give resistance to other drugs as well [61].

More realism could be achieved also on the clinical side. Longer time horizon could be considered, with more realistic clinical constraints, such as concentrated drug administration in day-hospital, rather than distributed over the week, and break periods in case of adverse circumstances.

More realism or specificity could be also introduced on the target of control. Multiple targets could be mixed in the objective function to be optimized, or some could be expressed in terms of solution constraints. For example, following the recent modeling literature on cancer, enforcing a static tumor volume while controlling time to tumor regrowth or phenotypic heterogeneity seem to be important targets, rather than tumor size [40, 42, 55].

On the side of the control method, gradient descend could be jointly used together with gradient-free exploration methods, such as those based on genetic algorithms, to avoid the costly computation of gradients at each step and, at the same time, allow a more global exploration of the control space.

Envisaging the eventual application in-vivo of our results, especially in the context of personalized medicine, ideas from model predictive control could be incorporated to close the loop with the patient. The idea is to plan the optimal cure on a calibrated model initialized on the measured patient's state, apply the cure on a relatively short time horizon, and then replan it from the actual state reached by the patient.

Finally, modeling advancements should also be accompanied with experimental ones, and modeling could actually suggest experimental directions. For example, the two-drug resistance allocation trade-off introduced by Orlando and colleagues [10] and reviewed here (see section 2.2.2) seems to play a role in determining the optimal cure, so that its experimental identification is interesting (as discussed in section 4.1). In the context of cancer drug resistance, experiments targeted to evaluate multi-drug resistance are also very promising.