The role of dose-density in combination cancer chemotherapy

TL;DR

This paper uses a mathematical model to evaluate dose-dense chemotherapy protocols, finding that low-dose maintenance between cycles may be more effective than simply shortening treatment intervals.

Contribution

It introduces a multicompartment model to compare dose-dense protocols and suggests a novel strategy of low-dose maintenance to prevent tumor relapse.

Findings

Low-dose maintenance can arrest tumor relapse.

Shortening intervals alone may not be optimal.

Model-based evidence supports alternative scheduling strategies.

Abstract

A multicompartment mathematical model is presented with the goal of studying the role of dose-dense protocols in the context of combination cancer chemotherapy. Dose-dense protocols aim at reducing the period between courses of chemotherapy from three to two weeks or less, in order to avoid the regrowth of the tumor during the meantime and achieve maximum cell kill at the end of the treatment. Inspired by clinical trials, we carry out a randomized computational study to systematically compare a variety of protocols using two drugs of different specificity. Our results suggest that cycle specific drugs can be administered at low doses between courses of treatment to arrest the relapse of the tumor. This might be a better strategy than reducing the period between cycles.

Equations

$$\frac{d{P}_j}{dt}={\delta }_{1j}g(\vec{P})P_j-{\delta }_{5j}n(\vec{P})P_j+\sum _{i=1}^5{w}_{ji}(\vec{P})P_i-k_j(\vec{P},{\vec{C}}_s,{\vec{C}}_n)P_j,$$

$$\begin{array}[]{lr}\dfrac{dM}{dt}=(r-v)M+uS-w(P)M+pQ\vskip 12.0pt plus 4.0pt minus 4.0pt\\ \dfrac{dS}{dt}=-uS+vM\vskip 12.0pt plus 4.0pt minus 4.0pt\\ \dfrac{dQ}{dt}=w(P)M-(p+n)Q,\end{array}$$

$$\begin{array}[]{lr}\dfrac{dP}{dt}=rM-nQ.\\ \end{array}$$

$$\wp (t)=\frac{M(t)+S(t)}{P(t)},$$

$$q=\frac{r}{K}\left(1+\frac{p}{n}\right)$$

$$M^{\ast }=\frac{K}{1+v/u+r/n},S^{\ast }=\frac{Kv/u}{1+v/u+r/n}\text{and}{Q}^{\ast }=\frac{Kr/n}{1+v/u+r/n}.$$

$${\wp }^{\ast }=\frac{1+v/u}{1+v/u+r/n},$$

$$n=r\left(1-f\right)\left(\frac{{\wp }^{\ast }}{1-{\wp }^{\ast }}\right).$$

$$k(P,C)=\frac{\mu }{1+\sigma P/K}\left(1-e^{-\rho C}\right),$$

$$\begin{array}[]{lr}\dfrac{dM}{dt}=F_{M}(M,S,Q)-\dfrac{\mu_{n}}{1+\sigma_{n}P/K}\left(1-e^{-\rho_{n}C_{n}}\right)M\vskip 12.0pt plus 4.0pt minus 4.0pt\\ \dfrac{dS}{dt}=F_{S}(M,S,Q)-\dfrac{\mu_{n}}{1+\sigma_{n}P/K}\left(1-e^{-\rho_{n}C_{n}}\right)S-\dfrac{\mu_{s}}{1+\sigma_{s}P/K}\left(1-e^{-\rho_{s}C_{s}}\right)S\vskip 12.0pt plus 4.0pt minus 4.0pt\\ \dfrac{dQ}{dt}=F_{Q}(M,S,Q)-\epsilon\dfrac{\mu_{n}}{1+\sigma_{n}P/K}\left(1-e^{-\rho_{n}C_{n}}\right)Q,\end{array}$$

$$\frac{dC}{dt}=I(t)-k_{e}C,$$

$$I(t)=\{\begin{array}{rl}\frac{D}{{\tau }_a}& \text{for}t(\text{mod}T)\in [0,{\tau }_a)\\ 0& \text{for}t(\text{mod}T)\in [{\tau }_a,T)\end{array}.$$

$$C(t)=\{\begin{array}{rl}\frac{D}{k_e{\tau }_a}(1-e^{-k_{e}t})+a_N(t)& \text{for}t(\text{mod}T)\in [0,{\tau }_a)\\ \frac{D}{k_e{\tau }_a}(e^{k_e{\tau }_a}-1)e^{-k_{e}t}+a_N(t)& \text{for}t(\text{mod}T)\in [{\tau }_a,T)\end{array},$$

$$a_N(t)=\frac{D}{k_e{\tau }_a}(e^{k_e{\tau }_a}-1)e^{-k_{e}t}\sum _{n=1}^N{e}^{-n{k}_{e}T}$$

$$C_{\min }^{\ast }=\frac{D}{k_e{\tau }_a}\frac{e^{k_e{\tau }_a}-1}{e^{k_{e}T}-1}$$

Full text

The role of dose-density in combination cancer chemotherapy

Álvaro G. López1*, Kelly C. Iarosz2, Antonio M. Batista3, Jesús M. Seoane1, Ricardo L. Viana4, Miguel A. F. Sanjuán1

1 Departamento de Física, Universidad Rey Juan Carlos, Móstoles, Madrid, Spain

2 Institute of Physics, University of São Paulo, São Paulo, SP, Brazil

3 Department of Mathematics and Statistics, State University of Ponta Grossa, Ponta Grossa, PR, Brazil

4 Department of Physics, Federal University of Paraná, Curitiba, PR, Brazil

*alvaro.lopez@urjc.es

Abstract

A multicompartment mathematical model is presented with the goal of studying the role of dose-dense protocols in the context of combination cancer chemotherapy. Dose-dense protocols aim at reducing the period between courses of chemotherapy from three to two weeks or less, in order to avoid the regrowth of the tumor during the meantime and achieve maximum cell kill at the end of the treatment. Inspired by clinical trials, we carry out a randomized computational study to systematically compare a variety of protocols using two drugs of different specificity. Our results suggest that cycle specific drugs can be administered at low doses between courses of treatment to arrest the relapse of the tumor. This might be a better strategy than reducing the period between cycles.

Author summary

Chemotherapeutic drugs for the treatment of cancer are currently delivered in periodic cycles to allow healthy tissues to recover from their cytotoxic effects. Unfortunately, this recovery window allows the tumor to regrow between cycles as well. It has been observed that a reduction of the periodicity from three to two weeks can be beneficial without introducing higher toxicities. These new chemotherapeutic protocols are known under the term dose-dense chemotherapy. We perform a computational study to assess the limitations of the dose-dense approach and suggest possible alternatives.

Introduction

Overcoming the toxic side-effects of chemotherapeutic drugs on healthy tissues in the short term, as well as the development of drug resistance of cancer cells to such medications in the long run, are the two key cornerstones in the improvement of cytotoxic drug therapy [1, 2, 3]. These two limiting factors put before, the rest of our research efforts must concentrate on the particular arrangement of protocols. Randomized trials carried out along the last forty years have demonstrated that the way in which protocols are organized is of great relevance as well [4, 5, 6, 7]. For instance, breast cancer trialists have revealed that, in addition to the drug dose-intensity, the periodicity of the cycles of chemotherapy affects the patients’ outcome as well [8]. Accordingly, a sequential administration of two drugs can be more beneficial than the same drugs given alternately [4]. In particular, a reduction of the period between cycles of chemotherapy from three to two weeks has introduced moderate but statistically significant benefits in disease free survival at five years and overall survival [8]. Fortunately, this can be achieved without entailing higher toxicities with the assistance of granulocyte colony-stimulating factor.

The increase in the frequency of drug administration has been termed dose-dense chemotherapy, because when the dose is graphically represented against time, these protocols look more tight [9]. Of course, from a physical point of view, the temporal density of a certain magnitude can also be raised by increasing the amount of such physical magnitude while keeping the same distribution over time. On the other hand, dose-intensity is commonly defined in the field of oncology as an average by dividing the total dose of a drug administrated during a treatment by the duration of such treatment [8]. Thus, as currently defined and utilized, dose-density and dose-intensity are two mudding and redundantly interrelated concepts. To avoid any sort of confusion we unwind these two ideas by adopting the following criteria [10]. The term dose-density is here reserved to denote a variation in the periodicity of the cycles of chemotherapy when the dose per cycle is kept constant, while the term dose-intensity of a drug is hereinafter referred to an increase in the amount of dose administered per cycle.

From a theoretical point of view, the ground on which dose-dense chemotherapy settles is known as the Norton-Simon hypothesis. This hypothesis states that the rate of destruction by chemotherapeutic drugs is proportional to the rate of growth of an identical unperturbed tumor [11]. This statement is founded on the notion of mitotoxicity, which proclaims that chemotherapy acts more efficiently on the cells that are rapidly proliferating within a tissue [1]. Consequently, since only a fraction of the cells that shape a solid tumor are actively dividing [12, 13], and bearing in mind that more cells in proportion turn into the quiescent state so as to render the growth of the tumor sigmoid [14], the impact of chemotherapy should be stronger for smaller tumors. This effect would hinder the efficiency of chemotherapy because, at some point of the growth curve, the regrowth between cycles would prevail over the destruction caused by the drugs during that same cycle. Note that tumors tend to regrow more rapidly after a cycle of chemotherapy, becoming less sensitive at the end of such period of time, when the subsequent cycle arrives.

In summary, the progressive doubling time dilation of a growing solid tumor defies the log-kill paradigm that has dominated traditional chemotherapy over decades [15], which is being curtailed to non-solid malignancies nowadays. These supporting arguments notwithstanding, it has been recently demonstrated that dose-density might be pertinent even in the absence of any hypothesis affecting the particular nature of the tumor’s growth [10]. However, enthusiasm must be prevented at the same time, since the benefits of dose-dense chemotherapy are still being the subject of intense experimental research and scientific debate [16]. We recall that the biological mechanism by which a certain drug inflicts its damage in practice might not correspond exactly to the one for which it was originally conceived [17].

Bearing all these facts in mind, the appearance of dose-dense protocols for the treatment of solid tumors has been at all events an important breakthrough in breast and ovarian cancer. Furthermore, it is also posing crucial questions concerning the nature of solid tumor growth and its dynamical response driven by complex combination therapies [9]. Always bounded to experimental verification, needless to say, many of these questions can be better addressed within the framework of mathematical and in silico modeling [18]. In the present work we propose a multidimensional ordinary differential equation (ODE) model for cancer chemotherapy to investigate under which conditions dose-dense protocols are more or less advantageous. As far as we are concerned, the present work represents the first mathematical and computational approach to the study of the dose-dense principle in a cytotoxic combination therapy [19, 20, 21, 22, 23]. We disregard toxicity models in this study, since it is not our purpose to develop specific protocols for a certain class of tumors, which could sometimes be regarded with skepticism [24]. We would rather find general guiding principles that may aid clinicians to bias their decisions in their proposal of better forthcoming chemotherapeutic protocols.

1 Methods

1.1 Model description

A general multicompartment model for cancer chemotherapy can be built from preceding modeling attempts [25, 26, 27, 28, 29, 30]. Since the cell cycle comprises four phases (mitosis, Gap 1, synthesis and Gap 2) and cells can also be quiescent (Gap 0), we consider a model with five different compartments representing the cell populations in each of them (see Fig. 1(a)).

The use of different compartments is relevant since some drugs destroy only cells in a specific phase of the cell cycle. We assign to each of these compartments a number $j$ between one and five, respectively. Therefore, the five cell populations can be represented by a vector $\vec{P}(t)$, whose components can be written in tensor form as $P_{j}(t)$. We allocate the first compartment $(j=1)$ to mitotic cells, while quiescent cells are inserted in the last compartment $(j=5)$. Following previous works [28], we assume that only quiescent cells die from necrosis. On the contrary, the apoptotic pathway might be activated in any compartment [32]. However, tumors frequently evolve towards states of low apoptotic index quickly [28]. Consequently, and for computational simplicity, we sweep apoptosis under the rug of the proliferation term, by considering that this pathway simply reduces the net growth rate of tumor cells. Finally, every equation is equipped with two more terms. The first of them represents the rate of flow between adjacent compartments, while the second one models the action of cytotoxic drugs on each compartment. Broadly speaking, cytotoxic drugs can be classified as cycle specific (CS) or cycle non-specific (CNS). The former inflict their damage at a specific point of the cycle of the cell, while the second can destroy cells independently of their kinetic status. Therefore, the two vector variables $\vec{C}_{s}(t)$ and $\vec{C}_{n}(t)$ symbolize the different types of drugs that have to be delivered to the patient, respectively. The dimension of each vector is given by the number of drugs scheduled in the treatment. In summary, the differential equation governing the dynamics of each population $P_{j}(t)$ in the system can be written as

[TABLE]

where $\delta_{ij}$ is the Kronecker delta and $j=1,...,5$. As thoroughly described ahead, the functions $g(\vec{P})$ and $n(\vec{P})$ are nonlinear functions representing the fractional cell growth of the tumor cells in the mitotic compartment and the fractional cell death due to necrosis in the quiescent compartment. The functions $w_{ji}(\vec{P})$ and $k_{j}(\vec{P},\vec{C}_{s},\vec{C}_{n})$ are nonlinear functions representing the fractional cell flow between compartments and the fractional cell kill of chemotherapeutic drugs on each compartment, respectively. If the cell cycle is regarded as a Markov chain, the terms $w_{ji}(\vec{P})$ are related to the transition matrix. On the other hand, the detailed functions $\vec{C}_{s}$ and $\vec{C}_{n}$ depend on the pharmacokinetical model used and the way in which drugs are administered. Note that this inhomogeneous five-dimensional model is constituted by two well differentiated parts. The three first terms on the right hand side of equation (1) define the cell kinetics and include the turnover of cells from one stage of the cell cycle to another, as well as the birth of tumor cells and their death from other causes different from cytotoxic treatment. The remaining fourth term on the right hand side of the same equation represents precisely the destruction of tumor cells by the cytotoxic agents. In the next lines we particularize the present model to meet the needs of the present work.

1.2 Kinetic model

In practice, one or more compartments of the previous general model can be merged into a single compartment, reducing the dimensionality of the problem (see Fig. 1(b)). For example, in the present work, there is only one CS drug and another CNS drug. Therefore, only three kinetic compartments are needed. If we assume that the cycle specific drug only targets cells in the phase of synthesis, as it frequently occurs, two compartments are required to represent cells that are actively dividing. One of these compartments contains cells that are in the synthesis phase at a particular instant of time $t$. This cell population is represented henceforth by the variable $S(t)$, while the remaining cells comprise all those cells that are in the complementary part of the cell cycle (mitosis, Gap 1 and Gap 2). We just gather these three compartments into a single one and refer to the cell population contained in it as mitotic, representing it mathematically by the letter $M(t)$. To conclude, one more cell population $Q(t)$ is included to model those resting cells that make up the quiescent compartment (see Fig. 2). Therefore, in this first study we shall utilize the three-dimensional model given by the following system of differential equations

[TABLE]

where $r$ stands for the net growth constant rate resulting from proliferation and apoptosis of the cells in the mitotic pocket and $n$ describes the rate at which quiescent cells die from necrosis. These terms correspond to the functions $g(\vec{P})$ and $n(\vec{P})$ in equation (1), which are chosen to be constants. The parameters $u,v$ are the flow constant rates of cells from the s-phase to the mitotic population and vice versa. The term $w(P)$ and $p$ are the rates describing the speed at which cells flow from the mitotic compartment to the $G_{0}$ phase and vice versa. Now, these functions correspond to the functions $w_{ji}(\vec{P})$ in the original model. The variable $P(t)$ represents the whole tumor cell population resulting from adding the cells in all the compartments, i. e., the relation $P(t)=M(t)+S(t)+Q(t)$ holds. Note that the total tumor cell population is governed by the net growth rate balanced by the rate of necrosis

[TABLE]

Most of the fractional cell functions are considered constant for simplicity, following previous works [28] and traditional quasispecies models [34]. However, the flow rate between mitotic cells and quiescent cells is more generally assumed to be governed by the linear term $w(P)M$. The reason why this is so is that we want the model to represent limited growth following a sigmoid function up to a certain maximum number of tumor cells $P^{*}=K$. This stable fixed value is frequently known as the carrying capacity of the tumor and we shall follow this terminology here. More specifically, we assume logistic growth [35, 37, 36] for mathematical simplicity. Nevertheless, other types of sigmoid functions, such as Gompertzian growth or a Richards’ type of growth can be used instead, since they are all generated by diffusive long range cell interactions [38]. In short, we consider that the flow rate of mitotic cells to the resting pocket obeys a linear function $w(P)=qP$, which gives rise to a logistic type of growth.

We now describe how to compute the model parameters describing the rates of flow $u,v,q$ between compartments in terms of the other parameters, which are more amenable to the experimental practice and for which we have several values in the literature [28]. These better known parameters are the constant rate of mitosis $r$, the turnover of mitotic to quiescent cells $p$, the proliferative index, here denoted as $\wp$ and defined as

[TABLE]

and the fraction of time $f$ that cells spend in the phase of synthesis.

The parameters $u$ and $v$ can be expressed in terms of $r$ by considering the time spent in each phase of the cycle by a human cell. These constant rates can be estimated by means of the following argument. If we name the total cell cycle time by $\tau$ and the cell cycle of the synthesis phase is denoted as $\tau_{S}$, while the remaining time of the cycle is written as $\tau_{M}$, we have $\tau=\tau_{S}+\tau_{M}$. In particular $\tau_{S}=f\tau$ and $\tau_{M}=(1-f)\tau$ are fulfilled. Now, for a geometric growth, the doubling time is inversely proportional to $r$, which allows us to write $u=r/f$ and $v=r/(1-f)$. On the side of the variables we have that at equilibrium the rates are zero $dM/dt=dS/dt=dQ/dt=0$ and, consequently, also the relation $dP/dt=0$ holds. Now, we name the values of the cell populations at equilibrium as $M^{*}$, $S^{*}$ and $Q^{*}$ and assume that they are different from zero. Then, bearing in mind that $P^{*}=M^{*}+S^{*}+Q^{*}=K$, it is straightforward to obtain from the second identity appearing in equation (2) that $S^{*}=M^{*}v/u$, and from equation (3) that $Q^{*}=M^{*}r/n$. Using these two results and the equations obtained from the conditions $dM/dt=dS/dt=0$, the constraint

[TABLE]

is implied. Furthermore, making use of equation (2) we can also derive the size of each cell population at equilibrium, which reads

[TABLE]

From these three values it is easy to check that the proliferative index at equilibrium is written as

[TABLE]

from which we can derive the following parameter relation for the constant rate of necrosis $n$ as

[TABLE]

In the numerical study conducted in the section containing our results we present all the parameter values of the model, which are derived from the elementary parameters $r,f,\wp^{*},p$ by means of the preceding equations.

1.3 The fractional cell kill of cytotoxic drugs

Once the kinetic model and its parameters have been established, we proceed to introduce the fractional cell kill of cytotoxic drugs. Just as a reminder, and because sometimes the term rate is loosely used to refer a constant rate [39, 40] or the fractional cell kill itself [40], we clarify this concept. It describes the rate at which a cell population is reduced to a fraction of itself after a fixed period of time. If the cell population presents a constant fractional killing, it means that it decays exponentially in time. The fractional cell kill should be then identified with the constant rate of the decay, being both in direct relation. Therefore, the fractional cell kill must be rigorously defined from a mathematical point of view as the rate at which the logarithm of a cell population $P$ decreases as a consequence of some biological process ($d\ln P/dt$), as it is done in previous works [40].

In the present case, such attrition is caused by cytotoxic agents. To work out the mathematical nature of the function $k(P,C)$ that describes the cell kill by a chemotherapeutic drug we appeal to past works on this subject [41]. There, the following function is proposed

[TABLE]

where $C(t)$ is the concentration of the drug at the tumor site and $\mu$ is the maximum fractional cell kill achievable by the drug, when it is ideally given at an infinite dose to sufficiently small (exponentially growing) tumors. The term $1-e^{-\rho C}$ is based on the Exponential Kill Model [33]. This mechanistic model was developed in conformity with in vitro data and has also been tested against in vivo experiments with mice [33, 36]. Contrary to traditional log-kill models, it exhibits saturation of the cell kill as the drug concentration is raised. The steepness of the dose-response curve is related to the sensitivity $\rho$ of a tumor cell line to a particular drug. Simply put, for a nascent tumor the parameter $\mu$ would be the maximum constant rate of destruction achievable at high doses. On the other hand, the parameter $\rho$ tells us how the response to the drug is modified as the cells become increasingly resistant. This reduction in sensitivity to the drugs is reflected by a smaller value of the parameter $\rho$. Finally, in another investigation [41] we derive and explain in great detail the Michaelis-Menten dependence of the cell kill rate on the total tumor size population $P(t)$. It has been introduced in the model in order to account for the Norton-Simon hypothesis, whose intensity is condensed in the parameter $\sigma$. As the tumor gets bigger and approaches its carrying capacity, the maximum fractional cell kill reduces from $\mu$ to the value $\mu/(1+\sigma)$, with $\sigma>0$.

Since we only deliver one cycle specific drug $C_{s}$ and another cycle non-specific drug $C_{n}$ in the foregoing investigation, the whole set of differential equations together with chemotherapy can be written as

[TABLE]

where the functions $F_{J}(M,S,Q)$ represent the kinetic equations for each compartment respectively, with $J=\{M,S,Q\}$. Note that a new parameter $\epsilon$ has been introduced in the third equation. It is linked to the extent to which non-specific drugs destroy cells in the quiescent phase. It is frequently assumed that, because of their mechanism of action, CNS drugs destroy quiescent cells to the same extent as they damage cycling cells [28]. From a logical point of view, this would be certainly so, even only by definition. However, although the mechanism of action can be clearly identified sometimes for a drug, on many occasions the death of the cell is only made effective once it reenters the cell cycle [17]. The notion of a cell kill that is totally independent of the kinetic status of the cell clearly challenges the paradigm of mitotoxicity. Even so, the fact that the proliferation rates in many chemosensitive human cancers are low clash with such paradigm, suggesting to some extent that the mechanisms of destruction might be different as those conventionally claimed [1, 17, 24]. Therefore, a value $\epsilon$ in the interval $[0,1]$ is deserved to modulate the effect of CNS drugs on resting cells.

1.4 Pharmacokinetics and protocols of chemotherapy

It remains to be introduced the pharmacokinetics of the present model, which determines the dynamics of the concentration of each drug $C(t)$ at the tumor site. As in previous works on chemotherapy [36, 37, 27], we assume a one-compartment model and first order pharmacokinetics, which in some situations can be considered as a good approximation to the clinical practice [42]. Hence, the differential equation governing the concentration of the drug simply reads

[TABLE]

where $I(t)$ is the function representing the rate of drug flow into the body (i. e. the instantaneous dose-intensity) while $k_{e}$ is the rate of elimination of the drug from the bloodstream. The half-life can be computed from this constant rate by taking $\ln 2/k_{e}$.

Concerning the drug delivery, we assume that every drug is administered intravenously at constant speed during a time $\tau_{a}$ in cycles of period $T$. Therefore, a typical protocol of chemotherapy can be schematically represented as in Fig. 3. If the total drug dose given in a course of chemotherapy is $D$, then the instantaneous dose-intensity $I(t)$ is mathematically expressed as the square function

[TABLE]

The solution to equation (11) with a drug input given by equation (12), when a number of $N+1$ chemotherapeutic cycles have been delivered, yields

[TABLE]

where $a_{N}(t)$ represents the decaying accumulated drug concentration at the $N+1$-th cycle of chemotherapy, which equals zero for the first cycle $(a_{0}(t)=0)$ and it is equal to

[TABLE]

for the subsequent cycles. As shown in Fig. 4, when the number of cycles tends to infinity ($N\rightarrow\infty$), the concentration reaches a steady state periodic oscillation $C^{*}(t)$ between a minimum value

[TABLE]

and a maximum value $C_{\max}^{*}=C_{\min}^{*}e^{k_{e}(T-\tau_{a})}$. Since the period between cycles of chemotherapy commonly spans a few weeks, whereas the drug is eliminated rather quickly in comparison (around a day), the accumulation of drug is not usually very important, unless $\tau_{a}$ is comparable to $T$.

2 Results

The most important part of the results presented here concerns the randomized computational study along which we compare several types of protocols in the six-dimensional parameter space spanned by the parameters $\mu_{s},\mu_{n},\rho_{s},\rho_{n},\sigma_{s}$ and $\sigma_{n}$. However, to acquire some intuition on the differences between cycle specific and cycle non-specific drugs, we thoroughly investigate the role of these drugs separately in the first place. These prior computations will enable a discussion on the role of the parameter $\epsilon$, which describes the relative effect of CNS drugs on resting cells, compared to their cycling counterparts. Moreover, they will permit us to uncover the evolution of the proliferative index along the treatment. As shown promptly, the proliferative fraction acts as a limiting factor to dose-densification when the quiescent compartment is barely chemosensitive.

As far as the kinetic parameters of the model are concerned, we give their values in this preamble. Such values shall be used for all the simulations, unless otherwise specified. We consider a tumor whose maximum rate of growth is $r=0.9~{}\text{weeks}^{-1}$ and with a value of the carrying capacity of $K=10^{9}$ cells, which are derived from experimental values in references [39, 36]. An extremely fast growing tumor would be one in which all cells were ceaselessly dividing through mitosis in the absence of apoptosis. Since the cell cycle of a human cell lasts approximately one day, this corresponds to an exponentially growing tumor with a constant rate value of $r=4.85~{}\text{weeks}^{-1}$. Therefore, we are considering a quite aggressive tumor for which the rate of apoptosis is rather low. The carrying capacity corresponds to a detectable premetastatic tumor mass of approximately one gram. We recall that limiting values of deadly tumor masses round numbers of $K=10^{12}$ cells [1]. The original size of the tumor is considered to be $P_{0}=9\times 10^{8}$ cells, divided among the three compartments as given by their equilibrium distribution, represented by equation $\eqref{eq:6}$. We take a typical value of $k_{e}=4.85~{}\text{weeks}^{-1}$ for the CNS drug and $k_{e}=9.55~{}\text{weeks}^{-1}$ for the CS one. These two values correspond to a half-life of approximately one day and less than a half a day, respectively [41]. The proliferative index is considered to be $\wp^{*}=0.35$, which is a reasonable intermediate value within the range appearing in previous works based on empirical data [28]. The value of $p$ is chosen to be $0.35~{}\text{weeks}^{-1}$, which is exactly the same value assumed in the last cited work, and which is also based on experimental evidence [43, 44, 45]. This value settles the constant rate of necrosis to approximately $n=0.3~{}\text{weeks}^{-1}$ by means of equation (8). Again, this amount is also of the order of other values appearing in the literature [28], which results in a necrotic fraction compatible with the experimental data gathered from in vivo experiments [46]. The value of $q$ is derived from equation (5) after an elementary computation, and results to be $q=1.95\times 10^{-9}~{}\text{weeks}^{-1}\text{cell}^{-1}$. To conclude, since the time spent in the s-phase of the cycle of a typical human cell is two fifths of the whole cycle, we have $f=0.4$. As we argued before $u=r/f$ and $v=r/(1-f)$, thus we obtain that $u=2.25~{}\text{weeks}^{-1}$ and $v=1.5~{}\text{weeks}^{-1}$.

2.1 Cycle specific and cycle non-specific drugs separately

We begin by studying the role of the administration of a CS drug alone. If desired, we consider $\mu_{n}=0.0~{}\text{weeks}^{-1}$ in the model equations (10). The value of the sensitivity is set to $\rho_{s}=1.0~{}\text{mg}^{-1}$, which is within values appearing in previous works for cells sensitive to methotrexate [33], an antifolate drug that acts during the s-phase of the cycle. For reasons explained somewhere else [10, 41], the maximum fractional cell kill is set to $\mu_{s}=4.0~{}\text{weeks}^{-1}$. As shown in previous works [41], this value corresponds to a moderate value of a drug that is capable of tumor shrinking. Finally, the parameter $\sigma_{s}=1.0$ is set, corresponding to a value for which the Norton-Simon hypothesis applies rather weakly [41]. It means that the rate of kill is reduced to half of its maximum value as the tumor gets close enough to the carrying capacity. A dose of a CS drug $D_{s}=90~{}\text{mg}$ is administered intravenously ($\tau_{s}=1~{}\text{min}$) every two weeks ($T_{s}=2~{}\text{weeks}$). The values of the dose and the drug periodicity are within values frequently appearing in standard protocols of chemotherapy111Information about standard protocols of chemotherapy has been drawn from $<http://www.bccancer.bc.ca/health-professionals/clinical-resources/chemotherapy-protocols>$. For example, the CS drug methotrexate can be administered at high doses of $D={50~{}\text{mg}~{}\text{m}^{-2}}$ for an adult male [47]. Since the average surface body area of an adult male is around $1.8~{}\text{m}^{2}$, a value of $D_{s}=90~{}\text{mg}$ can be estimated. Nevertheless, we recall that a great variability in the parameters $\mu$, $\rho$ and $\sigma$ is expected in practice [33]. For this reason, we design the randomized numerical trial presented in the next section.

As shown in Figs. 5(a) and (b), when a CS drug destroys tumor cells in its specific compartment, this cell population is naturally reduced. As cells in the s-phase are destroyed, less cells flow from the synthesis phase to the mitotic one, so that $M(t)$ mimics $S(t)$, although the variations are less pronounced, more smooth and present a delay of a couple of days. Once the drug is eliminated from the organism, the tumor cells in the mitotic pocket regrow during the remaining part of the cycle. As the dividing tumor cells are wiped out, the resting cells abandon their state of quiescence and reenter the cell cycle. This effect produces a continuous reduction of the quiescent compartment. As can be seen in Fig. 5(c), the regrowth of tumor cells between cycles of chemotherapy and the draining of resting cells reinforce each other contributing to an average increase of the proliferative index. However, this increase and the periodic oscillations around it are very modest for CS drugs.

Now we examine the effect of administering CNS drugs on the model. For this purpose, we consider a dose of drug $D_{n}=60~{}\text{mg}$ administered intravenously ($\tau_{s}=1~{}\text{min}$) every three weeks ($T_{n}=3~{}\text{weeks}$), which are typical for the drug doxorubicine, used in chemotherapy for locally advanced breast cancer. The parameter values used in these simulations are $\mu_{n}=4.0~{}\text{weeks}^{-1}$ for the maximum fractional cell kill, $\rho=1.0~{}\text{mg}^{-1}$ for the sensitivity and $\sigma=1.0$ for the parameter simulating the Norton-Simon effect. In these simulations we assume that CNS drugs do not affect the quiescent compartment, thus $\epsilon=0$. As shown in Fig. 6(a) and (b), CNS drugs destroy tumor cells in the whole dividing compartment, producing a greater tumor cell decay. However, the regrowth between cycles is also higher now, and since the quiescent cells are resuming their cell cycle, these effects synergize to produce oscillations that grow in amplitude shooting up the proliferative index. In Figs. 6(d) and (e) the same results are investigated, by considering a certain destruction of quiescent cells by the CNS drug, as is commonly accepted [28, 17]. However, following a principle of prudence, we consider a small value $\epsilon=0.2$, which is respectful with the paradigm of mitotoxicity as well [1]. This phenomenon produces an even stronger boost of the proliferative index. Now its oscillations are smaller, because the transfer of cells from the quiescent compartment reduces slightly. We note that a possible and simplified way of modeling the most strict version of the Norton-Simon hypothesis, would be to consider the parameters $\epsilon$ and $\sigma$ both equal to zero.

These results clearly suggest that, contrary to the philosophy of dose-density, allowing some time between cycles has the benefit of maintaining a high proliferative fraction. In turn, it is then expected that if we push the dose-density too much, protocols might become less effective, unless quiescent cells are as chemosensitive as cancer cells are. This idea is corroborated when we compute the same protocol in Figs. 7(a) and (b) with the same drug but at an increased frequency of one week. The resulting plot for the proliferative fraction (see Fig. 7(c)) indicates that some time is required for the quiescent cells to reenter the cell cycle and regrow, posing a limit on the frequency reduction craved by dose-dense protocols. Furthermore, and as a completely novel result, we show in Fig. 8 that when comparing the survival fraction at the nadir $s_{f}=P(t_{\text{nad}})/P(0)$ of two treatments for $\epsilon=0$, one being three times more dense that the other, there exist parameter values for which dose-density is not recommendable. We consider this statement an important revelation of the present model compared to previous investigations on this topic, where the dose-dense principle was unequivocally advocated [10, 41].

2.2 Randomized computational trial

Now we accomplish a randomized comparative study between several types of chemotherapeutic protocols. Inspired by clinical trials, we represent a vast number of different virtual patients (so to speak) by letting the parameters $\mu_{s}$, $\mu_{n}$ fluctuate in the range, $[0.1,4.5]$ with a uniform distribution. The parameters $\sigma_{s}$ and $\sigma_{n}$ are selected from the interval $[0.1,5.0]$ following the same recipe. Finally, the values of $\rho_{s}$, $\rho_{n}$ are picked at random in the range $[0.1,1.0]$ and $\epsilon$ is considered to be equal to $0.2$. This set of values is enough to model both chemosensitive and resistant cell clones [41], a considerable heterogeneity of drug responses to the drugs and a considerable degree of mitotoxicity. As far as we have investigated, doubling the size of the parameter intervals does not produce a substantial variation of our results. We use a total count of $400.000$ parameter values for each protocol, which has proven to be sufficient to give a convergent probability distribution at different scales and allow fast computations at the same time. Finally, for each point in the six-dimensional parameter space and every protocol, we compute the survival fraction at the nadir $s_{f}=P(t_{\text{nad}})/P(0)$, using the same initial conditions as in the simulations of the previous section. This effectiveness criterion has been adopted in previous works [28, 41]. Consequently, we recall that, since in the present study we are not explicitly modeling toxicity, the benefits or the effectiveness of a protocol should not be read from a clinical point of view, but from a theoretical one instead.

We now briefly discuss the design of the protocols. As seen in Fig. 9 seven different protocols are computed. A first reference protocol for which the two drugs (CS and CNS) are delivered concurrently every three weeks at doses appearing in the previous sections. The duration of the infusions also remains unchanged ($\tau=1~{}\text{min}$) with respect to the previous simulations. Then, two more protocols are computed to test the benefits of dose-densification in relation to the specificity of the drug. The former protocol is twice denser in the CS drug, while the latter is more tight for the CNS drug. The frequency is doubled in order to test protocols that are limiting in the time scale frequently used in the clinical practice (order of weeks). Now, the two next protocols are conceived to test an alternative way of avoiding tumor growth between cycles, compared to dose-densification. For this purpose, we use the CS-drug administered at constant total doses, and increase progressively the frequency of drug administration. Therefore, the fourth protocol for which the CS drugs is given at half the original dose ($D_{s}=45~{}\text{mg}$), but delivered every one and a half weeks. Then, the same total dose of $D_{s}=360~{}\text{mg}$ is delivered continuously. In this last case the drug concentration is very small to destroy the tumor, but enough to halt its regrowth between cycles of the CNS drug. Finally we introduce two additional protocols to compare the relative benefits of sequential against alternating protocols. To some extent, this evaluation allows us to ascertain the good performance of our model, since they are similar to real life trials carried out in the recent past [4, 9].

We begin by discussing the differences between the sequential and the alternating arms. As can be seen in Fig. 10, the model echoes the conclusions discussed in previous studies [9]. The number of events with low survival fractions are much smaller for the alternating protocol than for the sequential arm. In particular, we find that the latter is superior by more than a $10\%$ of events with one-log kill. This means that sequential administration is a better strategy, insofar as it is more dense in each drug. However, and as discussed in [11], this is especially true in this model because one of the drugs of the protocol is more effective than the other. Clearly put, our model does not tell the difference between the particular biochemical mechanisms of each drug, and therefore sequential and alternating protocols can not be distinguished by any means, as long as such drugs cause similar tumor cell destruction.

Concerning dose-dense protocols, the results appearing in Fig. 11 leave no doubt about the benefits of dose-densification for CNS drugs, where we obtain around a $13\%$ of events with one more log kills compared to the reference protocol. Unexpectedly, an increase in the frequency of drug administration does not seem to be a recommendable strategy if the CS drug is considered, since it barely augments the number of events with small survival fractions. In contrast, the number of events with survival fractions closer to one is noticeably higher. Probably, the reason behind this phenomenon is that when the cells in the synthesis phase are erased, it is useless to keep on trying to destroy this compartment, since it is already empty. Thus, this effect would be the dose-dense analog to the variation of the dose-response curves commonly observed for CS drugs, which limits the benefits of their dose intensification.

Regarding the possible benefits of protocols where the CS drug is given more frequently but at a constant equal total dose, we have found the following results. When the frequency is doubled and the dose per cycle is reduced to half of its control value, we see a modest increase in the protocol. However, when the CS drug is delivered continuously, despite being at quite low doses, we obtain a tremendous increase in overall tumor reduction. More specifically, a $20\%$ increase in cases bellow one-log kill is recorded. As we are about to show by going to smaller scales of the survival fraction, the impact is even more striking for this protocol. Accordingly, and as a colophon to this randomized numerical trial, we compare the relative benefits obtained for the third and the fifth protocols, which have been the two most outstanding settings.

These results are contrasted in Figs. 13(a). Not only the continuous CS protocol is about $6\%$ more destructive under one-log kill than the CNS-dense protocol. When the results are plotted by zooming inside the region of survival fractions between [math] and $0.1$, we find that more than $10\%$ of the cases for the CS-continuous protocol are under two-log kills. Moreover, by pushing further our scrutiny, we find that a non negligible $6\%$ of the cases lie under three-log kills. We believe that these results are very relevant from a theoretical point of view and, as discussed in the next section, have the virtue of suggesting a possible change of paradigm concerning the future of dose-dense chemotherapy.

3 Discussion and Conclusions

In the present paper we have introduced a quasispecies model of tumor growth for cancer chemotherapy that allows to represent tumor masses formed by cells with different kinetic status. Even though the dimensionality has been kept low in this initial study, the extension of the model to reproduce heterogeneous tumors is straightforward and can serve to explore the evolution of resistance in maintenance or palliative chemotherapy, in addition to adjuvant therapy. Moreover, the incorporation of other players to the dynamical scenario, as for example the cell-mediated immune response or other types of cells from mesenchymal origin, can enable a more accurate modeling of mixed therapies [37] and permit the simulation of toxicity as well.

The kinetic variability of our mathematical ODE model has allowed us to unveil limitations on the principle of dose-densification that were unobserved in previous works [10, 41]. These conclusions have been drawn on two different fronts. Firstly, our study sheds light into the role of specific drugs in dose-dense chemotherapy. As strongly suggested by our simulations, dose-densification of cycle-specific drugs is envisaged to be counterproductive. On the other hand, even though CNS drugs and dose-densification are expected to be generally beneficial, too dense protocols might not allow the quiescent compartment to flow into the replicating one, rendering dose-density less effective than a priori foreseen. In fact, we wonder if this effect is behind the failure of dose-density that has been reported recently in trials of dense cancer chemotherapy [16]. Evidently, this effect arises as long as the destruction in the quiescent compartment is not too high compared to the mitotic one. Otherwise, as $\epsilon$ is made close to $1.0$, dose-density for CNS drugs would be unconditionally recommended.

Secondly, our study suggests that other possible strategies can be devised to arrest the regrowth of the tumor between cycles. In the present work, we have proposed that the use of CS drugs at low doses administered continuously in addition to CNS drugs periodically delivered can be a much better schedule than dose-densification in the CNS drug alone. The authors are aware that continuous protocols can lead to tachyphilaxis and might also be too toxic [48], even for CS drugs. Therefore, it is not our desire to suggest this specific strategy as a viable optimal choice, but rather to set off a debate concerning alternatives to dose-densification, capable of producing the same effect in a more forceful and safe way.

As an example, in previous works [10] we reasoned that the use of targeted therapies using cytostatic agents to arrest the tumor revival between cycles might serve to complement cytotoxic chemotherapeutic drugs. Certainly, this way of tackling the problem presents other difficulties [10], but we note that a change in the point of view of dose-dense chemotherapy is entailed by this perspective. In lieu of trying to achieve maximum fractional cell kill by simply reducing the time between cycles, we would rather put the tumors to sleep during the meantime. In this way, stasis and destruction would cooperate to produce an enhanced overall destruction of the tumor. Among other potential benefits, this approach could allow to extend recovery times between cycles of chemotherapy. We should not forget that in order to damp the increased toxicity of dose-dense protocols, growth factors are currently delivered to the patients.

Thus, to conclude, we retake once last Paul Ehrlich’s old adage and suggest that the periodically thrown cannonballs of cytotoxic chemotherapy might be complemented with the magic bullets of some cytostatic targeted therapy administer in between [49]. Then, the following hesitation immediately arises. In order to produce the aforementioned effect without halting the rebound of the collaterally damaged healthy tissues, can cytostatic drugs be readily designed to impact more selectively the tumors than cytotoxic ones?

Aknowledgments

This work has been supported by the Spanish Ministry of Economy and Competitiveness and by the Spanish State Research Agency (AEI) and the European Regional Development Fund (FEDER) under Project No. FIS2016-76883-P. K.C.I. acknowledges FAPESP (2015/07311-7 and 2018/03211-6).

Additional information

Competing financial interests: The authors declare to have no competing financial interests.