|Home | About | Journals | Submit | Contact Us | Français|
The chemotherapy late intensity schedule is revised to account for tumor growth instability, where a small tumor cell fraction emerges that exhibits a higher proliferation rate than the parent strain. Modeling this instability as simplified two-population dynamics we find that: 1) if this instability precedes the onset of treatment, the slope of the linear increase of the drug concentration for the standard “Norton-Simon late intensity schedule” changes and the initial value of the dose strongly depends on the balance of the two tumor cell populations and on their distinct growth rates; and, 2) if the instability trails the initial treatment, the effective chemotherapeutic drug concentration changes as well. Both cases point towards necessary improvements of the “Norton-Simon late intensity” schedule.
Cancer therapy is aimed at targeting specific tumor (re)growth patterns to inhibit, or at least control, local tumor proliferation and metastatic dissemination. Mathematical modeling provides a promising tool by which to shed light onto these growth dynamics from a theoretical perspective, using both experimental and clinical data as input. Such works range from investigating various forms of exponential growth kinetics, e.g.,(1, 2) to that of describing a Universal growth law (3), with intriguing implications of the latter to radiotherapy (4). For cytotoxic chemotherapy, a well documented case in which theoretical reasoning had an impact on the therapeutic approach is that of breast cancer where, following the “Norton-Simon” hypothesis (5), standard treatment has been adjusted to a “dose dense” regimen; that is, the agent is given at a greater dose rate (i.e., multiple pulses of chemotherapy) to optimize efficacy, and to keep toxic side effects comparatively low. The underlying hypothesis is that, in comparison to the host organ, cancerous tissue exhibits a denser structure with higher fractal mass dimension, and is composed of smaller, relatively faster growing sub-seeds, all of which follow the Gompertz law. Differing from the traditional “log-kill” concept (6), Norton and Simon then argued, based on a body of experimental and clinical studies on solid cancers, that “therapy results in a rate of regression proportional to the growth rate of an unperturbed neoplasm of that size” (7, 8).
This hypothesis has been tested in a phase III clinical study led by the Cancer and Leukemia Group B (CALGB 97-41) [http://www.calgb.org/](9) This large prospective randomized trial demonstrated that shortening the time interval between each chemotherapy cycle while maintaining the same dose size resulted in significant improvements in disease-free and overall survival in patients with node-positive breast carcinoma. More recently, also a phase I study in patients with metastatic breast cancer confirmed the therapeutic results predicted by the Norton-Simon model (10). However, the issue of drug resistance is only considered at the qualitative level within the framework of the Norton-Simon model, contrary to some leading work addressing this issue performed by Goldie and Coldman (11, 12). In a more recent theoretical work, (13) Coldman and Murray investigated a stochastic model of non-cell cycle phase specific cancer chemotherapy for exponentially growing tumours including the development of drug resistance and the effect of chemotherapy on normal cells. The authors concluded that early intensification, a feature shared with the Norton Simon model, is a common aspect of successful regimes where drug resistance is likely. We note that the role of intra-tumor heterogeneity and its implications for the evolution of drug resistance and cell kinetics has already been considered by Gardner (14), who proposed a kinetically tailored treatment (KITT model). Finally, Monro & Gaffney (15) recently proposed a model aiming at a unifying description based on well founded approximations, i.e. Gompertzian tumour growth, log cell kill by chemotherapy and Luria-Delbruck mutation to resistance (this assumption was based on the Nobel Prize winning work of Luria and Delbruck in 1943 (16), who showed that bacterial cultures developed resistance to bacteriophages at random). They were able to show that palliative continuous chemotherapy achieves optimal results for intermediate (rather than high) dosage levels with later (rather than earlier) intervention, due to the effects of competition between resistant and sensitive cells. However, their approach focused only on the different tumor response to therapy and not on the different growth rate (both before and after therapy) that represents intrinsic tumor heterogeneity.
Here, we present evidence that cancer's inherent clonal diversity results in non-homogenous growth patterns which in turn - even for its simplest form of two-population dynamics - will require a modification of the “late intensification schedule” that has been put forward by Norton-Simon (17). The next section introduces the underlying growth instability model and its relevance for the current standard schedule for both situations, i.e. if the instability arises prior or after the regimen. We also propose an approach more closely related to clinical therapy and compare the results with data on survival probability.
Let us recall that for a homogeneous system the general growth law can be written as
where N is the number of individual cells at time t and α is their specific growth rate. The growth laws can be classified according to the order of power expansion in α of the generating function (18)
For instance, exponential growth, the Gompertz and logistic laws as well as various power laws proposed to describe tumor growth can be obtained by using
If N is a population composed of n sub-populations, each proliferating with a different specific growth rate αi, i = 1,….,n,
then eqn.  can be restored, provided α is assumed to be the average of αi, i.e.
For simplicity, let us consider a cancer system where, at time tδ following a progression event, a small cancer cell subpopulation emerges which grows at a distinctively different replication rate. At any time one can write
Where N1 (t) and N2 (t) are the number of cells belonging to the two subpopulations, evolving according to eqns.:
where θ (t - tδ) is 1 for t ≥ tδ and 0 for t < tδ.
The initial value is given by N(0) = N1(0), at t = tδ the cell population N2 (t) starts to grow and, since the instability emerges when a small portion of tumor cells is more rapidly growing, let us consider N1 (tδ) >> N2 (tδ).
To clarify the implications of the previous equations, let us assume that the two populations both follow the Gompertz law (GL), i.e.
This translates into the idea that the same type of cancer, in its sub-clones N1 (t) and in N2 (t), follows the same functional form of the growth law but different values of the specific growth rate. In the simplest model of this instability, with K = K1g = K2g, the difference in the specific growth is related to the saturation value of the cell number N (t → ∞) = N∞ given by
According to previous eqns. [8,9], different specific growth rates are obtained for N1∞ ≠ N2∞. In order to evaluate the effects of the choices of the parameter values, various configurations will be considered.
In fig. 3a we show, for the same set of parameters used in fig. 2a, the dramatic impact of the unstable system when N2 (tδ) is 10-6 smaller than N1 (tδ), for tδ = 500 days. The effect of increasing the delay time is shown in Fig. 3b by considering the same set of values of the parameter as in Fig. 3a and tδ = 1000 days. Let us note that a similar analysis can be performed easily with other growth laws, for example the Universal law proposed by West et al (19, 20). The results are expected to be qualitatively the same as those obtained using the Gompertz law, which is appropriately used here to relate the effects of tumor instability to chemotherapy regimen.
To investigate the potential effects of such tumor instability on treatment, let us consider the `Norton-Simon' hypothesis with regards to its `late intensity' schedule for chemotherapy. Considering a monoclonal cell population that grows according to eqn., and assuming that the dose-response per cell to a non-specific chemotherapeutic cycle of treatment can be described by a function C(t), the time evolution of the system is given by the following equation:
The tumor stops growing, and eventually regresses, when the therapeutic dose, delivered at time t*, is larger than a critical value given by
Moreover, assuming that, due to therapeutic impact, the cancer cell number decreases exponentially with a rate ω, then
which implies that, substituting in eqn. ,
i.e. the time-dependent dose, according to the `Norton-Simon' hypothesis, depends on the tumor growth rate α(t) of a single Gompertzian growth. By substituting for N the expression given in eqn.  and collecting the constant and the time dependent terms, we obtain that C(t), in the case of Gompertzian growth rate, is a linearly increasing function of time:
where A = ω + Kg ln and B = ωKg and
Therefore a linear increase of the chemotherapeutic drug dose with time can reduce the number of tumor cells exponentially. However, chemotherapeutic treatments based on the previous considerations assume single population dynamics.
The cancer populations N1 and N2 follow Gompertzian growth kinetics with the corresponding specific rates and drug concentrations, that is
and therefore for the time evolution of the total number of tumor cells, N(t), after onset of treatment, i.e t > t* one arrives at the following equation:
We will now discuss which modification of the `late intensity schedule' can be relevant for an unstable, i.e. heterogeneous tumor exemplified by two-population dynamics, according to the previous equations.
Let us first assume that tumor instability starts after the beginning of the therapy treatment, i.e. tδ > t*. In this case, for t* < t < tδ, the time evolution is due to population N1 (t) only and the corresponding drug concentration for a late intensity schedule is
Note that, in this time interval, the growth rate of population N1 (t) is that of the `one-population Gompertzian law' with saturation value N∞, because N1 (t) = N(t).
By recalling that
therefore, by assuming that for t > tδ the exponential decrease of the total tumor cell number elicited by the therapy is still effective, i.e.
it turns out
On the other hand, for a `one-population Gompertzian growth', with saturation value N∞, one has
Let us now assume that t* > tδ, i.e. the treatment (delivered at time t* follows the onset of instability (at time tδ). Accordingly, the tumor grows following the `one-population Gompertz law' up until time tδ (N (t) = N1 (t) for t < tδ) when the population N2 (t) emerges with a different specific proliferation rate.
Through the aforementioned eq.(22), the drug concentration for t < t* turns out to be
The role of tumor instability is crucial for the time behavior of the drug concentration because it determines different populations, N1 (t*) and N2 (t*) at time t*, and different slopes in the linear time dependence of . In fact, during the interval tδ < t < t* the two populations, N1 and N2, follow a Gompertzian growth with different specific rates. As we shall see in the next section, the presence of instabilities, both before and after the time at which the therapy has been administrated, postulates further improvements of the “Norton-Simon late intensity” schedule based on an uniform rate of regression, ω, of the number of cells according to eqn.. The previous eqns.[27,30,31] show that, for the two population dynamics following a Gompertzian growth, the average drug concentration required to obtain such a result is different from the drug concentration for a single tumor Gompertzian population. This general result is obtained by mimicking the differential response of the two populations to the treatment by different drug concentrations C1 and C2 and the same ω (see eqns.[17,18,19]). However, one can consider another approach to the problem, that could, to some extent, be closer to the real setting, where a “cocktail” of drugs is given which are expected to be effective against a particular tumor on the basis of previous empirical data. Indeed, in the therapy modeled according to the previous treatment, one expects that the drug concentration is uniformly delivered among cells and that C1 = C2 = C but the “differential” effectiveness against a cancerous clones, i.e. the reduction (for unit time) of cancer cells in the two populations, is translated in a different rate of regression for populations N1 and N2, that is ω1 ≠ ω2.
To be more precise, let us consider the growth equations
and let us assume that, after the time t* at which the treatment started, the therapeutic result is
In this discussion the delay time, tδ, is neglected but it can be easily included. From the previous equations for population N1 it turns out that (one assumes for simplicity that K1g = K2g = K)
which quantifies the linearly increasing drug concentration required for population N1 to exhibit a decrease in time with constant rate ω1. Analogously, by assuming that ω2 is time independent, for population N2 the required drug concentration turns out to be
In other terms, by requiring that the same dose C(t) can reduce the two populations with different constant rates, one obtains that, at any time, the following condition should be fulfilled:
Apart from the trivial solution ω1 = ω2 and
eqn.  cannot be satisfied except when relaxing the requirement that some reduction rate is time independent. This implies that eqn. for the effect of the drug concentration on population N2, for example, must be modified with a time dependent ω2, i.e.
obtaining for ω2 the linear differential equation:
where the last term has been obtained by imposing that there is a decrease of population N1 with constant rate ω1. By defining the parameter
the solution of the previous equation can be written as
with the boundary condition ω2 (t*) = ω1 - β.
This result allows us to compare the predicted evolution of the two different populations with experimental data.
Using the model described above, we now evaluate its effects on chemotherapeutic treatment and its predictive power by using medical data on breast cancer survival probability.
The ratio in eqns.[27,30] describes the effects in the drug therapy with respect to the single Gompertzian growth when tumor instability starts after the beginning of the therapy, i.e. tδ > t*. For t* <t <tδ this ratio is 1 and therefore one can set t* = 0. In Fig.4 DII is plotted for tδ = 60 days and ω = 0.1 per month, and different values of are considered. For a one-population Gompertzian growth, of course, DII = 1. With this set of parameters DII changes less than 10%. Notice that, when , the value of DII can be less than one for some time intervals. This is due to the fact that at the onset of instability the total tumor mass splits in two subpopulations with different specific proliferation rates but also with different cell numbers at tδ. Therefore a depletion of the large population N1 (tδ) occurs and it requires some time for the growing instability, N2, to compensate this initial negative step. Then ratio in eqns.[30,31] (which is, by definition, equal to 1 for t < tδ) describes the effects for onset of tumor instability prior to the beginning of chemotherapy. D strongly depends on the specific proliferation rates of the two populations and on their relative fraction of the total cancer cell population at the onset of instability and of treatment. In Fig.5a the ratio D is plotted as a function of time for t* = 365 days and tδ =100 days for different values of .
In Fig.5b, D is plotted for the same tδ and t* values, and for the same ratios , but changing the relative population at the onset of instability and at saturation (i.e. when t →∞). Depending on the value of , the ratio D increases from 20% to 70% (see figure captions for details).
We note that the values of the ratios D and DII depend on many parameters, which should be properly taken into account in a clinical context. On the other hand the presence of instabilities, both before and after the time at which the therapy has been administrated, postulates further improvements of the “Norton-Simon late intensity” schedule. For this reason, in the previous section, a more general approach with different rates of regression for populations N1 and N2 (ω1 ≠ ω2) has been also considered. It can be applied to more realistic clinical cases. For example, it is well known that primary breast cancer cells can express a variety of surface receptors which correlates with clinical outcome. In particular, during their proliferative activity, breast cancer cell populations may start expressing receptors for hormones, mainly for estrogens and progesterone (21, 22). Let us denote, in a breast cancer, the original cell population N1, that fails to express estrogen receptor, as ER-, and the population N2, that starts to express such receptors, as ER+. To be able to compare our model with experimental data, one has to relate the survival probability of a patient at time t, P(t), after the end of chemotherapy which occurred at time te, assuming a survival tumor cell fraction at time te computed according to eqns.[36,37], in which t* = 0 : therefore
where eqn. has been used. By a simple calculation
At least for short time intervals after the end of chemotherapy, the survival patient probability should be proportional to the regrowth rate at the end of the therapy: larger regrowth rate corresponds to smaller survival probability. Therefore one can write
For not too large t - te, the value of the parameter α1 can be fitted by approximating the curve ER- (ER/PgR absent) in fig. 1 of Colleoni et al. (23). This givesα1 ≈ 0.071 per year. To evaluate the survival probability with ER+ we use the corresponding formula with α2 (te) (note that a different te is simply a rescaling of β). By comparison with the corresponding clinical data, the parameter β can be estimated by trial and error methods and turns out to be β~ 0.08 per year for times in the range of 1-4 years. On the basis of such values, the survival probabilities can be estimated and the comparison with the corresponding clinical data is given in Table 1. These results are in good agreement with the data reported in (23) for the curve ER+ (ER/PgR present) and suggest a generalization of eqn. in exponential form, that is
which gives the results reported in Table 1 by using β~ 0.13 per year. One should note that the previous comparisons have been done by considering ω1 time independent and without any optimization of the parameter β. This preliminary analysis does not show the saturation of the patient survival probability observed at long time intervals, however it suggests a simple relationship between the drug dose required to obtain a given regression rate and the survival probability of the patient which could be generalized in analogy to the so called “linear-quadratic model” used in radiotherapy. We shall discuss this point in detail in a forthcoming paper.
In here we propose a model aimed at gaining more insights into how tumors respond to therapy, taking into account some features which are normally disregarded, and in particular:
We specifically investigated the case of an initially monoclonal tumor in which, at time tδ, a secondary, faster replicating strain emerges from the parent tumor cell population. The resulting “two population dynamics” is described in a way that is somewhat different from the usual interspecific competition of the two populations for survival and access to resources (prey-predator model, see, for example (24)) since the details about the dynamical evolution of the two distinct populations in the tumor required by such a more “traditional” approach would be unavailable, at least until the tumor is surgically approached and a biopsy taken for histological evaluation. In our case, therefore, the competition between the two cell populations is encoded in the ratio and in the boundary condition for long intervals: formally the difference in specific growth is related to the saturation value of the cell number, N (t→δ) = N∞, given by eqn. .
Since the maximum number of individuals N∞, the ratio and the delay in the onset of growth between the two subpopulations results in different specific growth rates, a variety of distinct realistic situations can be described without necessity of other “ad-hoc” assumptions. For instance, the occurrence of “dormancy” phases during tumor development, which has been simulated by Speer et al.(25) assuming that some parameters describing the overall Gompertzian growth may experience stepwise variations (responsible for the appearance of “plateaus” in the tumor growth curve), maybe generated by proper choices of the delay time tδ (see Fig 3b). We also investigated how the conventional “late-intensity” chemotherapy schedule proposed by Norton and Simon, based on the underlying assumption of single-population dynamics (with the parameters of a unique homogeneous Gompertzian growth), has to be modified when tumor time evolution is non uniform. We showed that in both case, i.e. if the instability starts before or after the onset of treatment the chemotherapeutic dose strongly depends on the balance of the two populations and on their specific growth rates.
Finally, by introducing a drug concentration non linear in time, the model was evaluated with some available clinical data, comparing its predictions for the subpopulation ER+ that emerges from the original breast cancer cells. Comparing the predicted survival probability following therapy, we conclude that the model is capable to prefigure the different responses on the basis of a small number of parameters, whose values can be easily fitted from clinical data. Most importantly, our simplified analysis indicates already that a global, average chemotherapeutic approach will not be successful if growth instabilities were to occur. It further argues that a better phenomenological and theoretical understanding of such therapeutic effects requires a more accurate description of these instabilities.
CG gratefully acknowledges the support from Regione Piemonte (PSF 2007, 2008) and University of Torino (60% 2007, 2008). This work has been supported by NIH grant CA 113004 (TSD) and by the Harvard-MIT (HST) Athinoula A. Martinos Center for Biomedical Imaging and the Department of Radiology at Massachusetts General Hospital.