The Role of Cell Density and Intratumoral Heterogeneity in Multidrug Resistance

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 (IC₅₀) 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 IC₅₀ 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 IC₅₀ 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 IC₅₀ 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.

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, C₀ lines for IC₁ and IC₂). 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.

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 C_Drug1 and C_Drug2. The spikes in Fig. 1B labeled as C_Drug1 and C_Drug2 show tumor homogeneity, because most cells have similar traits. The C_Drug1 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, C₀: IC₁ 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, C₀: IC₁ curve) die out, the system would not restore this distribution (Fig. 1B, C₀: IC₂ 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.

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 IC₅₀ 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 IC₅₀ 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 IC₅₀ values (solid vs. dashed lines, Fig. 2A), with steeper slopes (slope₁ vs. slope₂, Fig. 2B), and a smaller distance between the two curves (Δ₁ vs. Δ₂, Fig. 2A) with different levels of resistance. These three factors influence the drug efficacy.

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 (C₁ or C₂), 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.

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.

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)$|⁠.

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 M₁ and M₂, respectively). As mentioned above, there is an indirect assumption of interaction between the M functions, which is qualitatively expressed by the activity timing of M₂ 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.