Recent data have demonstrated that cancer drug resistance reflects complex biologic factors, including tumor heterogeneity, varying growth, differentiation, apoptosis pathways, and cell density. As a result, there is a need to find new ways to incorporate these complexities in the mathematical modeling of multidrug resistance. Here, we derive a novel structured population model that describes the behavior of cancer cells under selection with cytotoxic drugs. Our model is designed to estimate intratumoral heterogeneity as a function of the resistance level and time. This updated model of the multidrug resistance problem integrates both genetic and epigenetic changes, density dependence, and intratumoral heterogeneity. Our results suggest that treatment acts as a selection process, whereas genetic/epigenetic alteration rates act as a diffusion process. Application of our model to cancer treatment suggests that reducing alteration rates as a first step in treatment causes a reduction in tumor heterogeneity and may improve targeted therapy. The new insight provided by this model could help to dramatically change the ability of clinical oncologists to design new treatment protocols and analyze the response of patients to therapy. Cancer Res; 73(24); 7168–75. ©2013 AACR.

Major Findings

We suggest that chemotherapeutic treatment acts as a selection process in the effective drug concentrations range, whereas genetic/epigenetic alterations act as a diffusion process that results in trait spread based on different stress signals. Application of our model to cancer treatment suggests that reducing the alteration rate as a first step in treatment causes a reduction in tumor heterogeneity and may improve targeted therapy.

Quick Guide to Equations and Assumptions

A structured population model

Here, using a continuous deterministic population model, we aim to describe the problem of multidrug resistance (MDR) by including key mechanisms that control the dynamics of the cancer-related system. Such a model should be viewed as a simplification of the full, complex system. We consider all variables and parameters (e.g., death, mutation rates, etc.) to be functions of the resistance level (i.e., traits), in which a trait is assumed to be a continuous variable. On the basis of the approach of Lorz and colleagues (16) of using a system of integrodifferential equations parameterized by the variable of the resistance level (x), we provide a basis for structured population models designed to estimate the intratumoral heterogeneity as a function of the resistance level (x ∈ [0, 1]) and time (t), with density, ρ, dependence. The evolution over time of each subpopulation n(x, t) is based on the rate of cell division, r(x), the natural rate of death, d(x), and the initial system conditions, n(x, 0). In addition, there are contributions from all other subpopulations, depending on the genetic and epigenetic alteration function, M(x, y), that can be thought of as a spread or diffusive effect. Of note, θ(x) denotes the fraction of cells with trait x that can carry out new modifications, where $0\, \le \,\theta (x)\, \le \,1$⁠. These alterations, M, can be either short term or long term in relation to the daughter cells. If a drug is included, there is the additional rate of drug-induced death as a function of the dose, c(x, α), in which the drug is assumed to be applied uniformly in time. Of note, $f[\rho (t)]$ and $g[\rho (t)]$ are included to add density dependence to the division and death rates, where $\mathop {\lim}\limits_{\rho \to \infty} {{f(\rho)} \mathord{\left/{\vphantom {{f(\rho)} {g(\rho)}}} \right. \kern-\nulldelimiterspace} {g(\rho)}}\, = \,0^ +$⁠. Mathematically, the system is modeled by equation A. All variables are listed and defined in Supplementary Table S1. For further mathematical descriptions and assumptions of the system, see the Supplementary Materials.

formula

We demonstrate in detail how the model can qualitatively predict three common biologic observations based on density limitation, but with different assumptions about the cell division, death, and alteration rates of cancer cells. In the first case, we study the cell density response as a function of the cell division and death rates, in which the drug-induced death rate could be generalized as a function of the drug dose and the resistance level. We show how the relationship between these rates can create a heterogeneous tumor over time, and illustrate how different treatments could decrease the variation in the heterogeneity by selection, in the effective concentrations range. In the second case, we discuss the effects of genetic/epigenetic permanent alterations that occur in different fixed rates. In the third case, we present the contributions of additional temporal changes that can evolve from specific environmental causes. For full details about the parameters and numerical results, see the Supplementary Materials.

Resistance to chemotherapy remains a major cause of the failure of cancer treatment. This resistance results from numerous mechanisms. The “traditional” understanding of multidrug resistance (MDR) and its driving mechanisms oversimplifies the complexity of a perturbed cellular cancer network and focuses on several pathways/gene families. From that perspective, drug resistance is associated instead with the induction of drug efflux, activation of DNA repair, variations of target proteins, decreased drug uptake, altered metabolisms, sequestration, and changes in apoptotic pathways (1, 2).

Recently, intratumoral heterogeneity has also been found to be a major facilitator of drug resistance. Intratumoral heterogeneity refers to differences between cancer cells originating within the same tumor. Many primary human tumors have been found to contain genetically distinct cellular subpopulations reported to be mainly the result of stochastic processes and microenvironment signals (reviewed in ref. 3). In addition to the genetic differences within a tumor, therapeutic resistance can also be caused by several other nongenetic processes, such as epigenetic changes associated with chromatin modification or DNA methylation (4–7). One study of these processes was done by Kreso and colleagues (8, 9), who found, in a system with a single genetic clone, that there was functional variability among the tumor cells. Clearly, the integration of both genetic and nongenetic assumptions as well as heterogeneity should be included in the design of new experimental and computational models to have a better description of and ultimately, a solution to the problem of MDR.

To study intratumoral heterogeneity, we consider a mathematical modeling approach. As mentioned above, drug resistance is far from being a black-and-white, resistant or not resistant phenomenon. Accordingly, a continuous variable is a more appropriate way to describe, estimate and measure the resistance level. This estimator is a key variable in any mathematical representation of the MDR system and without it a comprehensive mathematical model cannot be developed.

Several direct and indirect approaches have been suggested to estimate the drug resistance level, depending on the type of data that is analyzed. For instance, in in vitro experiments, the dose-response assay (e.g., the MTT assay) can quantify the number of surviving cells after exposure to different drug concentrations for a certain time period and can be presented by “killing curves.” The 50% inhibitory concentration (IC50) values can be defined as the drug concentrations required to reduce cell viability to 50% of the untreated control population. Thus, for example, the resistance level can be described here by the IC50 value of each clone in the global population. A similar trend in killing curves would be expected, to some extent, for other drugs with similar features (targets, mechanisms, etc.). A linear generalization of such an approach would be the number of different drugs that can be separately applied to those cells and yet the cells still survive, in which the level of resistance can be calculated as a score of two variables: the number of drugs and the IC50 value of each drug. A nonlinear generalization would be the survival percentage of the treated population with drug combinations administered at the same time point. In all of these cases, the higher the score, the higher the resistance level.

Unfortunately, most clinical data do not include the IC50 values and the in vitro conclusions have not led to success in the clinic (10). Usually, clinical data include the physiologic properties that describe the progress, extent, or severity of a tumor (“staging”). All assignments of cancer stage are made at the time of diagnosis, before any treatment is given and thus cannot directly assess the resistance level. Combining clinical data with gene expression and survival data from the same patients can help to categorize them as “good” or “poor” responders, and a score for their resistance level can be calculated. Accordingly, any theoretical model should include subpopulations with resistance levels that can vary within the interval between “good” and “poor” scores.

The number of cellular mutations has also been proposed as a way to characterize resistance level. Because of the stochastic nature of the mutation process, there are mutations that do not necessarily contribute to cancer progression, and are not essential to the resistance level of a tumor. Yet such mutations certainly increase the intratumoral heterogeneity. Of course, once the number of passenger mutations accumulates to a certain level, they can be expected to have a global effect on tumor growth and sensitivity to certain drugs (11). The number of mutations does not necessarily go hand in hand with the resistance level, but rather the type of mutated pathways affects the evolution of MDR. For instance, mutations in the apoptosis pathway cause a decrease in the death rate (12), mutations in the RAS–RAF pathway cause increased cell proliferation and resistance to apoptosis (13), and mutator genes increase genetic alterations throughout the genome. Moreover, there are certain genes that promote genetic stability, including DNA repair genes, DNA damage sensor genes, and cell-cycle checkpoint genes. Changes in these stability genes affect the mutation rate and therefore this rate is not necessarily a constant parameter over time and space (14).

Many mathematical models related to MDR have been developed (reviewed in ref. 15). Most of these models quantify the resistance level as proportional to the number of mutations. Some models quantify the resistance level by considering the expression of a single gene. Our work can be viewed mainly as an extension of the mathematical model of Lorz and colleagues, who developed a physiologically structured population model that is structured according to a continuous variable that can represent the level of a phenotypic trait that has been selected as relevant to describe population heterogeneity (16). Lorz and colleagues assumed that the cell division rate is dependent on the cell density of normal cells, whereas cancer cells are assumed to grow exponentially with a low mutation rate. Inspired by their model, we propose an extended mathematical model for cancer cells with cell density dependence on cell division and death rates, and account for both genetic and epigenetic changes. Moreover, we describe the drug-induced death rate as a function of the dose and the trait. This work is important for optimizing the efficacy of treatment. It has obvious implications for understanding the complexity of cancer growth dynamics and for the analysis of drug resistance data.

Heterogeneity arising from treatment

The first case we study is a hypothetical scenario in which the system cannot have genetic or epigenetic changes, and its dynamic is solely driven by cell division and death rates (Fig. 1A). Because the mass of the tumor is limited by the density functions [Supplementary Materials, Eq. 7, mathematical proof can be found elsewhere (J. Greene, O. Lavi, M.M. Gottesman, and D. Levy, submitted)], the model predicts a stable distribution of cells with different traits after a finite time without any treatment (Fig. 1B, C0 lines for IC1 and IC2). The administration of chemotherapy to such a system leads to changes in the original distribution and a selection process that depends on the drug-induced death rate, $c(x,\alpha)$⁠. This rate of death is a function of the drug dose, α, and the resistance level, x.

Figure 1.

Heterogeneity arising from the administration of chemotherapeutic drugs, without genetic/epigenetic changes. A, rates of cell division, r(x); natural death, d(x); and drug-induced death for a given dose, c(x) (see Supplementary Materials for all functions and simulation details]. All rates are functions of the trait x. The two types of treatments are administered separately, in which drug 1 is a more effective treatment than drug 2. B, numerical predictions of the net growth [n(x, t)] after a certain time (t = 12.5), when no alterations (e.g., mutations) occur, with different initial conditions (solid and dashed lines, IC1 vs. IC2). In addition to the two treatments, the control (growth without treatment) is also plotted and labeled as C0. This plot shows the selection process resulting from treatment with respect to the trait, and demonstrates the time scale it would take to develop a fast-growing tumor. C, initial conditions, IC1: n(x, t = 0) = 1 for all x (solid curves); IC2: n(x ≤ 0.25, t = 0) = 0, or n(x > 0.25, t = 0) > 0 (dashed curves). D and E, cell density, ρ(t), for the two types of treatments: drugs 1–2, with the initial condition of IC1.

Figure 1.

Heterogeneity arising from the administration of chemotherapeutic drugs, without genetic/epigenetic changes. A, rates of cell division, r(x); natural death, d(x); and drug-induced death for a given dose, c(x) (see Supplementary Materials for all functions and simulation details]. All rates are functions of the trait x. The two types of treatments are administered separately, in which drug 1 is a more effective treatment than drug 2. B, numerical predictions of the net growth [n(x, t)] after a certain time (t = 12.5), when no alterations (e.g., mutations) occur, with different initial conditions (solid and dashed lines, IC1 vs. IC2). In addition to the two treatments, the control (growth without treatment) is also plotted and labeled as C0. This plot shows the selection process resulting from treatment with respect to the trait, and demonstrates the time scale it would take to develop a fast-growing tumor. C, initial conditions, IC1: n(x, t = 0) = 1 for all x (solid curves); IC2: n(x ≤ 0.25, t = 0) = 0, or n(x > 0.25, t = 0) > 0 (dashed curves). D and E, cell density, ρ(t), for the two types of treatments: drugs 1–2, with the initial condition of IC1.

Close modal

Drug-induced death rate based on a fixed dose.

If the dose is administered at a fixed value α, then $c(x,\alpha)$ can be replaced by $c(x)$⁠. The drug-induced death rates of two qualitative examples of different drugs are plotted in Fig. 1A as CDrug1 and CDrug2. The spikes in Fig. 1B labeled as CDrug1 and CDrug2 show tumor homogeneity, because most cells have similar traits. The CDrug1 curve has a spike at a higher resistance level (i.e., x value). Figure 1A shows that higher x values correspond to lower values of r(x), i.e., slower cell division rates, an assumption made when the model was developed. If only a few cells with different resistance levels survive the treatment, that would be sufficient to allow the distribution of intratumoral heterogeneity to return to its original stable profile (Fig. 1B, C0: IC1 curve), because the cell division and death rates have not changed. The time needed to reach that control state varies based on the treatment and its drug-induced death rate. Although if after the treatment, several subpopulations located in the nonzero range of the original stable profile (Fig. 1B, C0: IC1 curve) die out, the system would not restore this distribution (Fig. 1B, C0: IC2 curve). Kreso and colleagues provide evidence for a slowly proliferating cell population in primary human colorectal cancer cells that still retains potent tumor propagation potential, thus preferentially driving tumor growth after chemotherapy (8). A similar mechanism that reduces cell proliferation and induces a dormant nondividing state was reported by Lewis, in his review on MDR in microbiology (17). Kreso and colleagues noted that although cancer cells may have a uniform genetic lineage, individual tumor cells are functionally heterogeneous, with a wide variety in terms of growth dynamics and response to therapy. Dividing cells are most likely to be eradicated first, whereas the relatively slower growing/dormant cells are the major population during tumor reinitiation after chemotherapy. Of course, in the study of Kreso and colleagues there could be nongenetic alterations that could change one resistance level to another. In our theoretical case, we show that even if the cells are unable to mutate/change to another resistance level, it is possible to have the same dynamic, if only a few slowly proliferating cells that resist the drug remain. Once there is a density limitation, the distribution of cells as a function of the trait would be totally different from the control.

Drug-induced death rate as a function of the dose.

Drugs are usually categorized by their biochemical mechanisms of action, targets, molecular structure, gene expression, and in general by the cell response. However, from a mathematical perspective, there is another dimension to drug efficacy, which is the effect of the drug on the shape of the survival curve. Given two subpopulations, resistant and sensitive cells, the efficacy of a treatment could be estimated by the IC50 values, the slope of the survival curves, and the distance between the two curves of the two cell types (Fig. 2A and B). The slopes could be compared at their IC50 values (i.e., at their maximum values). The differences in treatment outcomes are not expressed solely by the percentage of cells that survive, but also by the effect on the remaining cells as a function of the dose and the resistance level. A better treatment would result in a lower number of surviving cells and a more homogeneous population. This can be expressed as survival curves with low IC50 values (solid vs. dashed lines, Fig. 2A), with steeper slopes (slope1 vs. slope2, Fig. 2B), and a smaller distance between the two curves (Δ1 vs. Δ2, Fig. 2A) with different levels of resistance. These three factors influence the drug efficacy.

Figure 2.

Death rate as a function of dose and resistance level. A and B, two common types of treatment (dashed and solid curves) for cells with different resistance levels (red curves for sensitive cells and blue for resistant cells). The cells may react differently to certain drugs. For example, a very toxic drug produces a lower IC50 value (solid vs. dashed lines), with a smaller distance between the two curves representing these cells (Δ1 vs. Δ2). Also, the slope can demonstrate the response to a treatment, in which the steeper slope is the desired one (slope1 vs. slope2). C and D, two examples of death rates $[c_1 (x,\alpha)\,vs.\,c_2 (x,\alpha)]$⁠, in which the heights and gradients of the functions vary according to the type of treatment. E and F, the cell density as a function of the dose and duration of a treatment for the cases described in C and D.

Figure 2.

Death rate as a function of dose and resistance level. A and B, two common types of treatment (dashed and solid curves) for cells with different resistance levels (red curves for sensitive cells and blue for resistant cells). The cells may react differently to certain drugs. For example, a very toxic drug produces a lower IC50 value (solid vs. dashed lines), with a smaller distance between the two curves representing these cells (Δ1 vs. Δ2). Also, the slope can demonstrate the response to a treatment, in which the steeper slope is the desired one (slope1 vs. slope2). C and D, two examples of death rates $[c_1 (x,\alpha)\,vs.\,c_2 (x,\alpha)]$⁠, in which the heights and gradients of the functions vary according to the type of treatment. E and F, the cell density as a function of the dose and duration of a treatment for the cases described in C and D.

Close modal

Data from survival curves, drug sensitivity assays, and apoptosis assays can be used to estimate the drug-induced death rate in at least two ways when incorporating the features that were mentioned above. Here, we describe two different functions [$c_1 (x,\alpha)$ and $c_2 (x,\alpha)$] to accommodate the two different biologic datasets that indirectly measure the drug-induced death rate as a function of the dose (α) or the trait (x). Each dataset supplies different pieces of information that are needed to estimate the death rate. Thus, different assumptions should be made when using each type of dataset. The first type involves data that measure the percentage of apoptotic cells for a range of resistance levels. Usually, these kinds of data are based on different staining of a tumor from a patient treated with a specific dose. Thus, the data include information on the death rate as a function of the trait, and can be used to infer the dependence of function $c(x,\alpha)$ on the dose. The second type of dataset is mainly from in vitro experiments. These types of data measure the survival curves for a range of drug concentrations for a few cell populations. The data include information on the death rate as a function of the dose, and can be used to infer the dependence of function $c(x,\alpha)$ on the trait. In general, both functions $c_1 (x,\alpha)$ and $c_2 (x,\alpha)$ would be expected to be nonincreasing functions of the dose, for a given resistance level. For high doses, the rate is at its maximum value, whereas for low doses, the rate is at its minimum value. Thus, there is a range of effective doses that this function of the death rate can meaningfully model and a sigmoid function is a very natural function to use in this case. The second variable of the death rate is the trait (i.e., resistance level). For a given dose in the effective range, the rate would be expected to decrease as a function of the trait. Theoretically, with these two types of data and the drug efficacy factors, the drug-induced death rate is illustrated in Fig. 2C and D. Figure 2C demonstrates the case in which this rate can be written as a product of the death rate for a typical dose and with the dose response as a decreasing sigmoid. The function in Fig. 2D, in contrast, demonstrates a case where the survival curves are used to describe the function $c(x,\alpha)$ as a complex, nonseparable dependence. All functions are listed in Supplementary Table S2.

Using our model, we estimated the cell density as a function of the dose and duration of a treatment (Fig. 2E and F). Our results show, in general, the desired sigmoid response to a drug as a function of the dose and describe two factors that are important to the efficacy of a treatment, the effective-concentration range and the critical-drug concentration (α*). The effective-concentration range is the range of drug concentrations in which the cell density is decreased. It is preferable to have as wide a range as possible. The critical-drug concentration is the value at which the cell density reaches a threshold for the first time (e.g., minimum tumor size that is detectable by the scanning instrument). Clearly, a low threshold is preferred. Furthermore, for each death rate (C1 or C2), we compared the cell density for two different drugs (Supplementary Fig. S1A and S1C). A more efficient treatment has a wider effective-concentration range, a lower critical-drug concentration, and a steeper cell-density slope, as occurs with drug 1 in both cases (Supplementary Fig. S1B and S1D).

In addition, during treatment the cell density is shown to be a more complex curve than a simple sigmoid curve. The success of a treatment is related to the shape of the corresponding density curve. For an example, see the density curves based on treatment of drug 1 versus drug 2 in both Supplementary Fig. S1B and S1D. The density response is a nonlinear function of the cell division rate, drug-induced death rate, and the density limitation in the effective concentrations range. For many drugs, the cell division rate also changes as a function of the dose, with higher doses usually leading to lower division rates. In practice, most drugs that cause a decrease in cell division rate also cause a decrease in the death rate. The strategy of combining multiple drugs, in which the first drug causes an increase in the death rate and the second drug causes a decrease in the cell division rate, may result in a better response, at lower drug concentrations.

Heterogeneity arising from both epigenetic and genetic changes

Spontaneous genetic mutations occur infrequently, but some cancer cells have increased rates of genetic change for various reasons. For example, in breast cancer, mutation of a single gene, Brca2, can lead to increased rates of mutations in other genes (18). The tumor microenvironment can also cause genetic changes. As tumors are exposed to repeated cycles of hypoxia and reoxygenation, DNA repair pathways are downregulated, thus leading to genetic instability (19, 20). Some of the mechanisms that lead to stable genetic changes are due to epigenetic alterations. Genome-scale genomic and epigenomic analyses have only recently revealed the complexity of the nonlinear relationship between these two types of changes and cancer. Although genetic mutations alter the sequence of DNA, epigenetic changes do not (reviewed in refs. 21, 22). Yet epigenetic changes can be inherited by causing mutations and silencing important genes by methylation (e.g., MGMT, CDKN2B, and RASSF1A), with a greater rate of epimutations than the rate of spontaneous genetic mutations (21). On the other hand, not all epigenetic changes are inherited, so they may lead to relatively rapid phenotype switching. In addition, these two types of changes can evolve from different causes. Genetic mutations can occur as a result of random processes or may be induced by stress (e.g., hypoxia, nutrient starvation, toxic molecules etc.), whereas epigenetic alterations are mainly linked to stress. It has also been shown that a central epigenetic control circuit can be disrupted by a genetic mutation (22). Clearly, the complexity of genetic and epigenetic networks has not yet been fully characterized, though it is clear that the outcome largely depends on external stresses. Thus, in the following cases we assume an indirect relationship between both networks that is based on an external stress.

Stable changes.

The second case considers the possibility of having epigenetic and/or genetic stable alterations when cells undergo division. Of note, M(y, x) denotes the probability that, given an alteration, a mother cell with trait y will yield a daughter cell with trait x. We considered M(y, x) to be a Gaussian distribution confined to (0, 1) with mean y and variance ${{\varepsilon ^2} \mathord{\left/{\vphantom {{\varepsilon ^2} 2}} \right. \kern-\nulldelimiterspace} 2}$⁠. The variance is one of the most important factors that must be considered because it reflects the impact of an alteration on the global system. Hence, the changes in variance can be written as a function of the external stress and time, $\varepsilon (stress, t)$⁠. As time progresses, different traits become advantageous or disadvantageous, and the overall dynamics are determined by various rates, mutation parameters, and the initial distribution of the cells. Figure 3 summarizes several examples of three variances $0\, = \,\varepsilon _0 \, < \,\varepsilon _2 \, < \,\varepsilon _1$ before (Fig. 3A) and after (Fig. 3B and C) treatment in a given effective dose. From these examples, it is obvious that different assumptions cause a divergence in predictions related to heterogeneity, even when using a single model. The general trend in these results indicates that mutations/epimutations (even at low rates) can serve as a diffusion process, over x, because the cells with a given trait may spread to an interval of x values that initially was set to zero. As expected, cases with high variation have high stable heterogeneity, whereas small variation slightly changes their dynamic from the nonaltered case. All other values and sequences of ϵ are constrained in that range, and affect the time to reach the control state $(\varepsilon _0 \, = \,0)$⁠.

Figure 3.

Heterogeneity arising from genetic/epigenetic changes. A, numerical predictions of the net growth [n(x, t)] after a certain time (t = 12.5), when no chemotherapeutic drugs were administered. The effects of genetic/epigenetic changes act as a diffusion process, where $\varepsilon _0 = 0,\,\,\varepsilon _1 = 0.01,\,\,\varepsilon _2 = 1$ with two initial conditions (solid and dashed lines, IC1 vs. IC2). For all simulation details, see Supplementary Materials. The initial conditions are shown in Fig. 1C. This graph shows the cases in which the diffusion potentially controls the tumor dynamic and its heterogeneity. B and C, net growth [n(x, t)] plotted at three time points: t1 = 2.5; t2 = 5; t3 = 12.5, with initial condition IC1. The administration of the drug takes place over a certain period of time (t1 < t < t2), for a given dose. These graphs demonstrate the effect of varying the alteration rates on the dynamic in the presence of anticancer drugs (⁠$\varepsilon _{\rm 1} \,{\rm =}\,{\rm 0}{\rm .01}$ in B; $\varepsilon _{2}\, =\, 1$ in C).

Figure 3.

Heterogeneity arising from genetic/epigenetic changes. A, numerical predictions of the net growth [n(x, t)] after a certain time (t = 12.5), when no chemotherapeutic drugs were administered. The effects of genetic/epigenetic changes act as a diffusion process, where $\varepsilon _0 = 0,\,\,\varepsilon _1 = 0.01,\,\,\varepsilon _2 = 1$ with two initial conditions (solid and dashed lines, IC1 vs. IC2). For all simulation details, see Supplementary Materials. The initial conditions are shown in Fig. 1C. This graph shows the cases in which the diffusion potentially controls the tumor dynamic and its heterogeneity. B and C, net growth [n(x, t)] plotted at three time points: t1 = 2.5; t2 = 5; t3 = 12.5, with initial condition IC1. The administration of the drug takes place over a certain period of time (t1 < t < t2), for a given dose. These graphs demonstrate the effect of varying the alteration rates on the dynamic in the presence of anticancer drugs (⁠$\varepsilon _{\rm 1} \,{\rm =}\,{\rm 0}{\rm .01}$ in B; $\varepsilon _{2}\, =\, 1$ in C).

Close modal

Temporal changes.

In the third case, the additional contribution of the M function is only valid for a specific period of time, when certain signals stress the tumor. We studied two possible mathematical modifications to system A. The first system assumes an accumulation of changes that occur with different time scales. Thus, we separate the alteration function into two functions. The first type of changes is related to hereditary alterations, and it is therefore applied for all time points, and the second type is applied only to a limited period of time because it is based on temporal epigenetic variations (marked as M1 and M2, respectively). As mentioned above, there is an indirect assumption of interaction between the M functions, which is qualitatively expressed by the activity timing of M2 and the ϵ of both functions. We studied the effect of ϵ, and its variation in each M function, on the dynamic with drug for a given dose, c(x). Supplementary Fig. S2 and Table S3 summarize several examples of different variances for both functions after stress.

These results highlight the tradeoff between changing the alteration rate and the drug selection process, which in turn affects the heterogeneity of the cancer cell population. If the alteration rate remains low during the entire process, the treatment would be more efficient, and the surviving cells would constitute a more homogenous population (Supplementary Fig. S2B). Once this selection is achieved, it is recommended that the second treatment should be more specifically targeted to the surviving population. However, if the first treatment is known to increase the alteration rate, a combination of drugs should be applied instead, in which one of these drugs (such as a drug that affects pattern of gene expression) aims to reduce this rate. On the other hand, if the initial alteration rate is too high and there is no effective treatment to reduce it, the model predicts (see Fig. 1A and red curve in Supplementary Fig. S2C) that increasing the alteration rate, at the time of treatment, should also result in a better response to the selective drug. In any event, at the end of the treatment, the rate of alteration must be decreased as much as possible.

The second modification to system A includes the temporal epigenetic changes by varying the parameters θ and ϵ based on the external stress $[\theta (stress, t),\,\varepsilon (stress, t)]$⁠, with a single M function (Supplementary Materials, system; ref. 12). During treatment, the results show that if θ increases, with a constant set of parameters, the density decreases but with a similar Gaussian shape of distribution. The cell density becomes significantly low when the ϵ reaches a certain level (Supplementary Table S4). A possible biologic interpretation for the parameter θ is the percentage of proliferating cancer cells based on stress.

Enormous progress has been made in understanding the molecular mechanisms leading to cancer, but our perspective on resistance to treatment is still in its infancy. A systems approach combines empirical, mathematical, and computational methods to gain an understanding of this complex phenomenon. Recent data show that the diverse intratumoral heterogeneity is mainly the result of stochastic processes and microenvironment signals. Thus, it is important to quantify the resistance levels in a nondiscrete manner. Resistance evolves parallel to the progression of cancer. However, not all abnormal changes contribute to the resistance level. Hence, the next experimental efforts should incorporate a characterization of the intratumoral heterogeneity based on the resistance level because this classification plays key roles in the development of MDR and in modeling such a system.

Here, we have presented a mathematical analysis of the MDR problem with the assumption that resistance is induced by adaptation to drug and/or environmental signals in a deterministic manner. We have used physiologically structured equations in which a continuous variable represents the resistance level. Such a variable may be generated by integrating a collection of physiologic measurements, as discussed in the introduction. Using our model, one can explain how the density plays an important part in the dynamic of a tumor. Adding this density dependence results in different distribution/heterogeneity.

Many observations of heterogeneous tumors have been reported (3), but not all of those tumors were sampled after drug treatments (23). Therefore, the drug was not necessarily the only component that affected the heterogeneity and the alteration rate. This rate varies widely across mammalian genomes (24), and this variation has important consequences for our understanding of the global evolutionary process, and in particular about cancer. We have sought to investigate such cases by examining heterogeneity as it relates to resistance level and time, in which the alteration rate varies over time with or without chemotherapeutic treatment and external stress. We assume that the alteration rate could be affected by either the treatment or the external stress, but the induced death rate varies only based on the treatment. Certainly, there is a possibility of having a mutation in the apoptosis pathway and therefore affecting the induced death rate, but there is no clear way of expressing this rate as a function of the trait.

Our results suggest that treatment acts as a selection process in the effective drug concentrations range, whereas genetic/epigenetic changes act as a diffusion process based on different stress signals. In a case with no changes and no treatment, the dynamic is determined solely by the cell density, cell division, and natural death rates as functions of the resistance level. We further examined the effects of two types of mathematical modifications. In this model, the resulting dynamic and response are expressed by only three main elements: density [$\rho (t)$], heterogeneity [$n(x, t)$], and alteration rate (ϵ).

Thus, the design of a new treatment protocol should also incorporate the individual status about these variables that should be estimated at the time of diagnosis. The common approach for selecting a protocol is to find and target genes that are significantly expressed in most cells in that tumor. Successful treatment is determined by the substantial reduction in the tumor size. The next obvious questions would be the following: What will happen with the remaining cells? Will these cells evolve into a more aggressive and resistant tumor or not? If these cells are low in number, how can we identify them? The model predicts that it is recommended to develop/bioengineer a treatment that targets the alteration rate of all cells to reduce this rate as much as possible before any other treatment, so the heterogeneity will not increase over time. In any event, it is predicted that it is better to first treat the small subset of cells with a high alteration rate than the bulk of the tumor using a specifically targeted therapy but with a low alteration rate. If this rate is very high and cannot be affected by external intervention, during treatment it is predicted that the increase in the alteration rate increases the spread to other traits that are not necessarily better in their survival; therefore, the density decreases and the entire system results in a better response. Increasing the alteration rate per se is an extreme approach; therefore, there is a need to bioengineer such a modification in a temporary and reversible manner. Another possibility is to increase the percentage of cells that could be exposed to alterations. By the end of the treatment, there is a need to reduce the alteration rate as much as possible.

Our findings imply that at certain levels of resistance, where the division rate is higher than the death rate (including the contribution of the alteration rate), this deterministic approach can be applied. It is possible that when the difference between those rates is low and there is additional stochastic noise, the stochastic results will differ from the deterministic prediction. Also, the study conducted in this article focused on the role of cell density and intratumoral heterogeneity in drug resistance. If several drugs are administered, a full study will be required to consider the multidimensional analog in which the resistance level is a vector in a space whose dimension is the number of drugs being administered. A complete study of this nature is beyond the scope of this article, as it will require detailed information about the cross-reactivity between the various drugs. An alternative approach would be to view the level of drug resistance as the magnitude of such a vector of resistance. In this case, changing the drug can be modeled by shifting the values of the distributions with respect to the trait. A full study addressing these two issues is left for future work.

No potential conflicts of interest were disclosed.

Conception and design: O. Lavi, J.M. Greene, D. Levy, M.M. Gottesman

Development of methodology: O. Lavi, J.M. Greene, D. Levy

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): O. Lavi, J.M. Greene, D. Levy

Writing, review, and/or revision of the manuscript: O. Lavi, J.M. Greene, D. Levy, M.M. Gottesman

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): O. Lavi, M.M. Gottesman

Study supervision: O. Lavi, M.M. Gottesman

The authors thank George Leiman for editorial assistance.

The work was supported by the Intramural Research Program of the NIH, Center for Cancer Research, National Cancer Institute and was supported in part by a seed grant from the UMD-NCI Partnership for Cancer Technology. The work of D. Levy was supported in part by the joint NSF/NIGMS program under grant number DMS-0758374 and by grant number R01CA130817 from the National Cancer Institute.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Fodale
V
,
Pierobon
M
,
Liotta
L
,
Petricoin
E
. 
Mechanism of cell adaptation: when and how do cancer cells develop chemoresistance?
Cancer J
2011
;
17
:
89
95
.
2.
Gillet
JP
,
Gottesman
MM
. 
Mechanisms of multidrug resistance in cancer
.
Methods Mol Biol
2010
;
596
:
47
76
.
3.
Saunders
NA
,
Simpson
F
,
Thompson
EW
,
Hill
MM
,
Endo-Munoz
L
,
Leggatt
G
, et al
Role of intratumoural heterogeneity in cancer drug resistance: molecular and clinical perspectives
.
EMBO Mol Med
2012
;
4
:
675
84
.
4.
Sanz-Moreno
V
,
Gadea
G
,
Ahn
J
,
Paterson
H
,
Marra
P
,
Pinner
S
, et al
Rac activation and inactivation control plasticity of tumor cell movement
.
Cell
2008
;
135
:
510
23
.
5.
Spencer
SL
,
Gaudet
S
,
Albeck
JG
,
Burke
JM
,
Sorger
PK
. 
Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis
.
Nature
2009
;
459
:
428
32
.
6.
Gupta
PB
,
Fillmore
CM
,
Jiang
G
,
Shapira
SD
,
Tao
K
,
Kuperwasser
C
, et al
Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells
.
Cell
2011
;
146
:
633
44
.
7.
Sharma
SV
,
Lee
DY
,
Li
B
,
Quinlan
MP
,
Takahashi
F
,
Maheswaran
S
, et al
A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations
.
Cell
2010
;
141
:
69
80
.
8.
Kreso
A
,
O'Brien
CA
,
van Galen
P
,
Gan
OI
,
Notta
F
,
Brown
AM
, et al
Variable clonal repopulation dynamics influence chemotherapy response in colorectal cancer
.
Science
2013
;
339
:
543
8
.
9.
Marusyk
A
,
Polyak
K
. 
Cancer. Cancer cell phenotypes, in fifty shades of grey
.
Science
2013
;
339
:
528
9
.
10.
Robey
RW
,
Massey
PR
,
Amiri-Kordestani
L
,
Bates
SE
. 
ABC transporters: unvalidated therapeutic targets in cancer and the CNS
.
Anticancer Agents Med Chem
2010
;
10
:
625
33
.
11.
McFarland
CD
,
Korolev
KS
,
Kryukov
GV
,
Sunyaev
SR
,
Mirny
LA
. 
Impact of deleterious passenger mutations on cancer progression
.
Proc Natl Acad Sci U S A
2013
;
110
:
2910
5
.
12.
Buchholz
TA
,
Davis
DW
,
McConkey
DJ
,
Symmans
WF
,
Valero
V
,
Jhingran
A
, et al
Chemotherapy-induced apoptosis and Bcl-2 levels correlate with breast cancer response to chemotherapy
.
Cancer J
2003
;
9
:
33
41
.
13.
Wong
KK
. 
Recent developments in anti-cancer agents targeting the Ras/Raf/MEK/ERK pathway
.
Recent Pat Anticancer Drug Discov
2009
;
4
:
28
35
.
14.
Koi
M
,
Boland
CR
. 
Tumor hypoxia and genetic alterations in sporadic cancers
.
J Obstet Gynaecol Res
2011
;
37
:
85
98
.
15.
Lavi
O
,
Gottesman
MM
,
Levy
D
. 
The dynamics of drug resistance: a mathematical perspective
.
Drug Resist Updat
2012
;
15
:
90
7
.
16.
Lorz
A
,
Lorenzi
T
,
Hochberg
ME
,
Clairambault
J
,
Perthame
B
. 
Populational adaptive evolution, chemotherapeutic resistance and multiple anti-cancer therapies
.
ESAIM-Math Model Num Analysis
2013
;
47
:
377
99
.
17.
Lewis
K
. 
Persister cells, dormancy, and infectious disease
.
Nat Rev Microbiol
2007
;
5
:
48
56
.
18.
Tutt
AN
,
van Oostrom
CT
,
Ross
GM
,
van Steeg
H
,
Ashworth
A
. 
Disruption of Brca2 increases the spontaneous mutation rate in vivo: synergism with ionizing radiation
.
EMBO Rep
2002
;
3
:
255
60
.
19.
Klein
TJ
,
Glazer
PM
. 
The tumor microenvironment and DNA repair
.
Semin Radiat Oncol
2010
;
20
:
282
7
.
20.
Reynolds
TY
,
Rockwell
S
,
Glazer
PM
. 
Genetic instability induced by the tumor microenvironment
.
Cancer Res
1996
;
56
:
5754
7
.
21.
Rando
OJ
,
Verstrepen
KJ
. 
Timescales of genetic and epigenetic inheritance
.
Cell
2007
;
128
:
655
68
.
22.
Shen
H
,
Laird
PW
. 
Interplay between the cancer genome and epigenome
.
Cell
2013
;
153
:
38
55
.
23.
Gerlinger
M
,
Rowan
AJ
,
Horswell
S
,
Larkin
J
,
Endesfelder
D
,
Gronroos
E
, et al
Intratumor heterogeneity and branched evolution revealed by multiregion sequencing
.
N Engl J Med
2012
;
366
:
883
92
.
24.
Hodgkinson
A
,
Eyre-Walker
A
. 
Variation in the mutation rate across mammalian genomes
.
Nat Rev Genet
2011
;
12
:
756
66
.