Drug-Induced Resistance in Micrometastases: Analysis of Spatio-Temporal Cell Lineages

Resistance to anti-cancer drugs is a major cause of treatment failure. While several intracellular mechanisms of resistance have been postulated, the role of extrinsic factors in the development of resistance in individual tumor cells is still not fully understood. Here we used a hybrid agent-based model to investigate how sensitive tumor cells develop drug resistance in the heterogeneous tumor microenvironment. We characterized the spatio-temporal evolution of lineages of the resistant cells and examined how resistance at the single-cell level contributes to the overall tumor resistance. We also developed new methods to track tumor cell adaptation, to trace cell viability trajectories and to examine the three-dimensional spatio-temporal lineage trees. Our findings indicate that drug-induced resistance can result from cells adaptation to the changes in drug distribution. Two modes of cell adaptation were identified that coincide with microenvironmental niches—areas sheltered by cell micro-communities (protectorates) or regions with limited drug penetration (refuga or sanctuaries). We also recognized that certain cells gave rise to lineages of resistant cells (precursors of resistance) and pinpointed three temporal periods and spatial locations at which such cells emerged. This supports the hypothesis that tumor micrometastases do not need to harbor cell populations with pre-existing resistance, but that individual tumor cells can adapt and develop resistance induced by the drug during the treatment.

Introduction

Drug resistance is one of the main impediments in effective anti-cancer therapy. While tumors may first respond well to chemotherapeutic agents, they often start growing back during or after the treatment period and become tolerant to the treatment. Several different intrinsic mechanisms of drug resistance have been postulated (Holohan et al., 2013; Housman et al., 2014; Cree and Charlton, 2017), including alteration of drug targets, changes in the expression of efflux pumps, increased ability to repair DNA damage, reduced apoptosis, elevated cell death inhibition, and altered proliferation. Some extrinsic factors contributing to drug resistance have also been postulated. A pivotal role can be played by the tumor microenvironment (Correia and Bissell, 2012; Sun, 2016) due to reciprocal communication between tumor cells and the surrounding stromal components. This includes interactions with fibroblasts and emergence of cancer associated fibroblasts, cross-talk with immune cells, and sensing cues from the extracellular matrix (ECM), as well as ECM remodeling. Tumor cells can modify their microenvironment by creating specific niches, including pre-cancerous, pre-metastatic or stem cell niches (Barcellos-Hoff et al., 2013; Hambardzumyan and Bergers, 2015; Huch and Rawlins, 2017). Additionally, the changes in tumor vasculature and interstitial fluid pressure may lead to creation of regions that are poorly penetrated by therapeutics, forming drug-limited pharmacologic sanctuaries or refugia (Cory et al., 2013; Puhalla et al., 2015) that influence tumor response to therapeutics. However these extrinsic factors are still not well understood.

Of particular interest is the heterogeneous and dynamically changing tumor microenvironment. As a result, spatially and temporarily variable gradients of drugs can be formed in the stroma, and tumor cells can be, therefore, exposed to different drug levels during the treatment period. It has been shown experimentally by Wu et al. (2013) that aggressive breast tumor cells can respond to drug gradients by migrating toward the regions of high concentration of doxorubicin and low cell population, and are able to adapt to high drug levels. Subsequently, they become tolerant to the drug and repopulate the region despite the high drug concentration. Fu et al. (2015) used mathematical modeling to investigate how the heterogeneity in drug penetration through the microenvironment effects tumor response to treatment. They showed that resistance arises first in cells located in regions with poor drug penetration, named pharmacological sanctuaries, and then populate areas with higher drug levels. Our own research showed that the non-homogeneous drug distribution within the tumor tissue that results in emergence of tissue regions with poor drug penetration but with normal oxygenation levels may lead to the emergence of acquired resistance (Gevertz et al., 2015; Perez-Velazquez et al., 2016). Similar results were previously generated using different mathematical models. Chisholm et al. (2015) investigated transient emergence of a drug tolerant population of cells using models of reversible phenotypic evolution, and concluded that a combination of non-genetic instability, stress-induced adaptation and selection are responsible for the emergence of weakly-proliferative and drug-tolerant tumor cells. Cho and Levy (2017) used a continuous model to show that cancer cells of different resistance levels can coexist in spatially-different areas in tumor tissue. Feizabadi (2017) used mathematical modeling to show that certain chemotherapy strategies are highly unsuccessful, and even damaging to the patient, under the assumption that the drug can induce resistance during the treatment period. Greene et al. (2018, 2019) developed mathematical approach to differentiate between spontaneous and induced resistance to drugs and proposed in vitro experiments that can determine whether treatment can induce resistance. The authors also designed optimized treatment protocols that can prolong the time before resistance develops.

Several experimental studies considered scenarios in which resistance is acquired by the tumor cells as a result of their exposure to the drug, either through epigenetic alteration, drug-induced genetic changes or non-genetic phenotype switching. Pisco et al. (2013) and Pisco and Huang (2015) used a combination of laboratory experiments and mathematical modeling to show that the emergence of multi-drug resistance in leukemic cells can be induced by the lasting stress response to the drug. In this case, the tumor cells exploited their phenotypic plasticity by modifying efflux capacity in a non-genetic but inheritable way. Goldman et al. (2015) and Goldman (2016) showed that exposure of breast tumor cells to high concentration of taxanes can induce phenotypic transitions toward chemotherapy-tolerant stem-like state that can confer drug resistance. Moreover, the authors demonstrated that this adaptive resistance process can be halted by carefully designed order of administered drug combinations. Other examples of drug-induced resistance pointed to modifications in chromatin configuration in lung cancer cells (Dannenberg and Berns, 2010; Sharma et al., 2010), changes in expression of stress adaptation-related proteins in prostate cancer cells (Ferrari et al., 2017), or switching to mesenchymal phenotype in melanoma cells (Su et al., 2017) as the mechanisms of increased cell tolerance to the drug. In all these studies, the exposure of tumor cells to chemotherapy caused non-genetic changes that allowed the tumor cells to tolerate drug treatment and evade drug-induced death.

In this paper, we used mathematical modeling to examine how individual tumor cells can adapt to alterations in drug distribution within the tumor microenvironment in order to acquire resistance to the drug. By tracking cells individually and reconstructing their behavioral history, we were able to provide insights into the complex spatio-temporal changes that occur in cell microcommunities and to explain how they avoid drug-induced death leading to therapy failure. In particular, we developed a concept of 3D spatio-temporal lineage trees that trace both genealogy and spatial locations of cells that survived the simulated treatment. This is an extension of classical lineage trees used to depict tumor clonal expansion in a form of a flat graph with an initiating cell connected to its children cells, that are connected to their descendants until the terminal nodes are reached (Navin and Hicks, 2010; Davis et al., 2017). The 3D spatio-temporal lineage trees allow us to identify the cells that drive a resistant phenotype in the sense that all their successors are resistant to the drug. The existence of such “special” cells has been reported previously under various names: drivers (Hutchinson, 2016; Nikbakht et al., 2016), superstars (Cheeseman et al., 2014a, b), or starter cells (Perez-Velazquez et al., 2015). We refer to these cells as precursors of drug resistance. The current study focuses on analyzing the behavior of individual resistant cells which is an extension of our previous work at the population level. This approach allowed us to develop novel evaluation methods, such as the 3D lineage trees, and also to identify the third microenvironmental niche prone to the emergence of resistant cells. Overall, this paper contributes to a better understanding of drug-induced resistance.

Materials and Methods

We used a hybrid multi-cell lattice-free model (MultiCell-LF) that combines the off-lattice individual tumor cells with the continuous description of oxygen and a cytotoxic drug. The cells can physically interact with one another, and respond to levels of oxygen and cytotoxic drug absorbed from cell’s vicinity. Low levels of oxygen (hypoxia) result in cell quiescence (Qiu et al., 2017). Exposure to the drug leads to cell damage – while we model this as a generic process, one can consider a more specific processes, such as the DNA damage (genotoxicity; Swift and Golsteyn, 2014) or cell membrane damage (lysis; Collins and Kao, 1989). Moreover, cells can become more tolerant to the drug they are exposed to, as shown in Sharma et al. (2010) and Pisco and Huang (2015). This cell’s response to different levels of oxygen and drug is a mechanism of cell adaptation to the microenvironment.

Drug and Oxygen Kinetics

The model is defined on a small patch of the tumor tissue with four irregularly positioned stationary vessels (Figures 1A,B). Both drug γ and oxygen ξ are intravenously supplied, diffuse through the tissue, are absorbed by the cells, and the drug is subject to decay. We model a small drug molecule of diffusivity comparable to oxygen diffusion (Schmidt and Wittrup, 2009) and with the same supply rate from the vessels. However, the drug is absorbed by the cells twice faster than oxygen (Schmidt and Wittrup, 2009). Drug γ and oxygen ξ kinetics are given by the following equations:

FIGURE 1

Figure 1. Components of the MultiCell-LF model. (A) A snapshot showing an irregular drug gradient (high level-white, low level-black) and individual tumor cells color-coded according to their viability (low viability-dark green, high viability-light green). (B) The same time snapshot showing oxygen gradient (high level-white, low level-black) and tumor clones marked by a unique symbol assigned to their initial ancestor cell (65 different symbols). Red circles in both panels represent four non-symmetrically located vessels. (C) A flowchart showing the relationship between cell behavior and (from left to right) the tumor microenvironment; oxygen levels that regulate cell proliferation or quiescence; drug levels that regulate cell survival, adaptation or death; upon cell division daughter cells inherit from the mother cell: the damage, tolerance level and half of the accumulated drug. Panel (C) adopted from Shah et al. (2016).

$${\displaystyle \frac{\partial \mathrm{\gamma }(\mathbf{x},t)}{\partial t}}=\underset{diffusion}{\underbrace{D_{\mathrm{\gamma }}\cdot \mathrm{\Delta }\mathrm{\gamma }(\mathbf{x},t)}}-\underset{decay}{\underbrace{d_{\mathrm{\gamma }}\cdot \mathrm{\gamma }(\mathbf{x},t)}}$$

$$-\underset{uptakebythecells}{\underbrace{{\mathrm{\rho }}_{\mathrm{\gamma }}{\displaystyle {\sum }_i}{\mathrm{\chi }}_R({\mathbf{\text{X}}}_i,\mathbf{x})}}+\underset{supply}{\underbrace{S_{\mathrm{\gamma }}{\displaystyle {\sum }_j}{\mathrm{\chi }}_{Rv}(V_j,\mathbf{x})}}(1)$$

$${\displaystyle \frac{\partial \xi (\mathbf{x},t)}{\partial t}}=\underset{diffusion}{\underbrace{D_{\xi }\cdot \mathrm{\Delta }\xi (\mathbf{x},t)}}$$

$$-\underset{uptakebythecells}{\underbrace{{\mathrm{\rho }}_{\xi }{\displaystyle {\sum }_i}{\mathrm{\chi }}_R({\mathbf{\text{X}}}_i,\mathbf{x})}}+\underset{supply}{\underbrace{S_{\xi }{\displaystyle {\sum }_j}{\mathrm{\chi }}_{Rv}(V_j,\mathbf{x})}}.(2)$$

where, D_γ and D_ξ are the drug and oxygen diffusion coefficients, ρ_γ and ρ_ξ are the drug and oxygen uptake rates, S_γ and S_ξ are the drug and oxygen supply rates, and d_γ is the drug decay rate. In numerical implementation, we take the smaller of the cellular demand ρ_γ /ρ_ξ and the current drug/oxygen level to assure that both concentrations do not fall below zero. Here, x represents the Cartesian coordinate system, X_i are the coordinates of discrete cells, V_j are the coordinates of discrete vessels, and χ is the indicator function of the local neighborhood of radius R around the cells X_i or of radius Rv around the vessels V_j, respectively:

$${\mathrm{\chi }}_R(Y,x)=\{\begin{array}{cc}\hfill 1\mathit{\hspace{1em}}\hfill & \hfill if\mathit{\hspace{1em}}||x-Y||<R\hfill \\ \hfill 0\mathit{\hspace{1em}}\hfill & \hfill otherwise\hfill \end{array}(3)$$

The initial condition consists of a stable oxygen gradient and no drug. The sink-like boundary conditions are imposed to implicitly represent the lymphatic system.

Individual Cell Dynamics

Each cell C_i(t) is defined by its position X_i(t), a fixed radius R, cell current age A_i(t) and cell maturation age $A_i^{mat}$. Cell division takes place upon reaching maturation age (30 h with 5% fluctuations between the cells to avoid synchronized cell division (Mehrara et al., 2007; Hafner et al., 2016), provided that the host cell is not overcrowded by other cells (14 cells within 2 cell diameters), and it is not located in the hypoxic areas (Qiu et al., 2017). If the level of oxygen in a cell’s neighborhood falls below the hypoxia level (5% of vascular supply), the cell becomes quiescent and will not proliferate (flowchart in Figure 1C). Upon division of cell C_i(t), two daughter cells C_i1(t) and C_i2(t) are created instantaneously. One daughter cell takes the coordinates of the mother cell, whereas the second daughter cell is placed near the mother cell at a random angle θ:

$$C_{i2}\left(t\right)=C_i\left(t\right)+R(\cos \mathrm{\theta },\sin \mathrm{\theta }).(4)$$

The current ages of both cells are initialized to zero, and their cell maturation ages are inherited from their mother cell with a small noise term. To preserve cell volume, the repulsive forces are introduced between overlapping cells (Gevertz et al., 2015; Perez-Velazquez et al., 2016). Since both daughter cells are placed in a distance equal to one cell radius, the repulsive forces are exerted to push the cells apart until they reach the distance equal to cell diameter. The repulsive forces are defined as overdamped springs:

$${\displaystyle \frac{d{\mathbf{\text{X}}}_i}{dt}}={\displaystyle \frac{1}{\mathrm{\nu }}}{\displaystyle \sum _{j=1\mathrm{\dots }M}}{\mathbf{\text{f}}}_{i,j}^{rep},where$$

$$\mathrm{\hspace{1em}\hspace{1em}}{\mathbf{\text{f}}}_{i,j}^{rep}=\{\begin{array}{cc}\hfill F^{rep}\left(2R-||{\mathbf{\text{X}}}_i-{\mathbf{\text{X}}}_j||\right){\displaystyle \frac{{\mathbf{\text{X}}}_i-{\mathbf{\text{X}}}_j}{||{\mathbf{\text{X}}}_i-{\mathbf{\text{X}}}_j||}},\hfill & \hfill if||{\mathbf{\text{X}}}_i-{\mathbf{\text{X}}}_j||<2R\hfill \\ \hfill 0,\hfill & \hfill otherwise,\hfill \end{array}(5)$$

here, ν is the damping coefficient, F^rep is the repulsive spring stiffness, and 2R is the spring resting length; M denotes the number of cells that overlap with X_i.

Upon division, both daughter cells inherit mothers’ damage level and tolerance level, while the drug absorbed by the mother cell is split into half between both daughter cells (Schmidt and Wittrup, 2009; Greene et al., 2019).

Cell Resistance Mechanism

Cell’s resistance mechanism is modeled as a competition between the level of drug-induced damage accumulated by the cell and the level of damage that the cell can withstand (tolerance) without committing to death. However, the cell can also adapt by increasing its tolerance level if it is exposed to the drug for a certain time (flowchart in Figure 1C; Gevertz et al., 2015; Perez-Velazquez et al., 2016).

Cell damage $C_i^{dam}(t)$ is increased proportionally to the newly absorbed amount of drug (we assume that the drug absorbed in the past has already exerted its damage effect). The rate of internal drug decay is taken to be the same as in the extracellular microenvironment. This drug-induced damage can be counterbalanced by cell natural ability for damage repair at rate p_γ (three times faster that the damage rate (Gevertz et al., 2015; Perez-Velazquez et al., 2016):

$$\frac{d{C}_i^{dam}}{dt}=\underset{newlyabsorbeddrug}{\underbrace{{\mathrm{\rho }}_{\mathrm{\gamma }}\sum _x{\mathrm{\chi }}_R({\mathbf{\text{X}}}_{\mathbf{i}},\mathbf{x})\underset{drugdecay}{\underbrace{(1-d_{\mathrm{\gamma }})}}}}-\underset{repair}{\underbrace{p_{\mathrm{\gamma }}{C}_i^{dam}(t)}}(6)$$

Cell exposure to high drug concentrations γ_exp (at least 1% of the vascular supply) for a long enough time t_exp (at least 2% of the cell cycle) results in cell adaptation and in increased cell tolerance to the drug $C_i^{tol}\left(t\right)$ (at a slow rate of Δ_tol = 0.01% of the baseline tolerance value; Gevertz et al., 2015; Perez-Velazquez et al., 2016):

$$\frac{d{C}_i^{tol}}{dt}=\{\begin{array}{c}\hfill \begin{array}{cc}\hfill {\mathrm{\Delta }}_{tol}\hfill & \hfill if{C}_i^{\mathrm{\gamma }}\left(t\right)\ge {\mathrm{\gamma }}_{exp}\mathit{\hspace{1em}\hspace{1em}}foratime\ge t_{exp}\hfill \end{array}\hfill \\ \hfill \begin{array}{cc}\hfill 0\hfill & \hfill otherwise\hfill \end{array}\hfill \end{array}(7)$$

where the amount of accumulated drug $C_i^{\mathrm{\gamma }}\left(t\right)$ depends on its continuous absorption (at a constant rate ρ_γ) and decay (at a rate d_γ):

$$\frac{d{C}_i^{\mathrm{\gamma }}}{dt}=\underset{uptake}{\underbrace{{\mathrm{\rho }}_{\mathrm{\gamma }}\sum _{\mathbf{x}}{\mathrm{\chi }}_R({\mathbf{\text{X}}}_{\mathbf{i}},\mathbf{x})}}-\underset{decay}{\underbrace{d_{\mathrm{\gamma }}{C}_i^{\mathrm{\gamma }}\left(t\right)}}(8)$$

Similarly as for Equation (1), in numerical implementation we take into account that cell demand for the drug may exceed the amount available in cell microenvironment, thus we take the smaller of the cellular demand ρ_γ and the current drug level to assure that drug concentration is non-negative.

Cell death depends on whether cell damage $C_i^{dam}\left(t\right)$ exceeds the tolerance level $C_i^{tol}\left(t\right)$. The dead cells are removed from the system. Thus, cell resistance to the drug depends on competition between the level of its damage and the level of damage the cell can withstand without committing to death.

Initially, there is a small micrometastasis consisting of 65 cells with the same baseline tolerance level, no damage, and identical proliferative properties. Each cell responds to the environmental cues, such as the level of sensed oxygen (that regulates cell quiescence or proliferation) and the amount of absorbed drug (that induces cell damage and modulates cell adaptability). The levels of drug and oxygen that the cell is exposed to during its lifetime can vary because the cell can move from one part of the tissue to another, and because drug gradient can change if the overall number of tumor cells changes. The full flowchart of cell behavior is shown in Figure 1C. During the simulation, we trace location and viability (the difference between tolerance and actual damage) of each individual cell. When the values of cell damage and tolerance steadily diverge over time, the cell is considered resistant to the drug.

Viability Trajectories of Individual Cells

Cell viability is defined as a difference between the level of cell tolerance to drug-induced damage and actually accumulated damage. The larger the viability value, the more non-responsive to the absorbed drug the cell is. The viability trajectory shows how the viability value is changing in time for a given cell and all its predecessors. The viability trajectory is generated backward starting from the last iteration at which the cell was alive, and going back the cells’ lifespan, the lifespan of that cell’s mother, the mother’s mother, and so on until the one of the initial 65 cells is reached (compare Figures 2B, 3A–D, 4B, 5).

FIGURE 2

Figure 2. Distribution of non-resistant cells. (A) A temporal histogram of dying cells. (B) Individual viability trajectories for each dying cell. (C) A density map of final locations of cells before they were annihilated by the drug-induced damage. Colors correspond to the number of cells killed at the given location during the whole simulation.

FIGURE 3

Figure 3. Individual cell adaptation to drug exposure. (A) Temporal evolution of viability trajectories of all tumor cells that survived the treatment (gray lines) and the average viability value (yellow dashed line). (B–D) Viability trajectories for individual cells from three distinct populations of cells with (B) superlinear, (C) linear, and (D) intermediate viability patterns; the insets show drug uptake by each cell from a given subpopulation during the final period of treatment corresponding to 29 cell cycles of the simulated time: (B) significant decrease in drug uptake, (C) constant drug uptake, and (D) small or short-time decrease in drug uptake. (E) Spatial configurations of tumor cells at six different time points; colors correspond to cell subpopulations from (B–D); black dots indicate cells that did not survive to the end of simulation. The larger red circles represent four vessels. All time points are shown in terms of the number of cell cycles (cc).

FIGURE 4

Figure 4. 3D Spatio-temporal routes of the representative resistant cells. Twelve resistant cells with a common predecessor were traced in space and time. All line colors correspond to cell viability categories: concave (green), linear (red), and intermediate (blue). (A) Three snapshots of routes traversed by the selected cells within a dynamically changing drug gradient in the background (high to low: yellow-red-blue; four large circles indicate the vessels). (Ai) initial overlap in cell routes (time corresponding to 33 cell cycles, 3 cells); (Aii) separation of individual cell routes and significant distances traveled by the cells (about 67 cell cycles, 9 cells); (Aiii) final map of cell routes showing very small changes in cell locations over an extended time (about 196 cell cycles, 12 cells). (B) Viability trajectories for the selected cells (i–iii correspond to the time points from A). (C) Projections on the xz plane of the individual routes for all 12 cells. A black star indicates the initial position of the common predecessor cell.

FIGURE 5

Figure 5. Lineage trees of survivors and cellular precursors of drug resistance. (Ai) shows all successors of cell #6 that survived the treatment (lineage tree of survivors). The cells which have no dead successors (the precursor cells of resistance) are shown by red arrows. The corresponding spatio-temporal routes of all surviving cells that arose from cell #6 are shown in (Aii) and indicate changes in cells position during the treatment, and final cell locations within the tissue with a drug gradient. The corresponding viability trajectories are shown in (Aiii); line colors correspond to those in Figure 4. (Bi–iii) shows the lineage tree of survivors, the precursor of resistance and, the 3D spatio-temporal routs, and the viability curves for all daughter cells generated by cell #26. (Ci) shows distributions of all precursor cells along the full set of viability trajectories, with spatial localizations of all precursor cells within the tumor tissue shown in (Cii); colors represent time of cell birth: early (green), middle (magenta), and late (cyan); red circles show vessel locations. All time points indicate time corresponding to the number of cell cycles.

Classification of Viability Trajectories and Cell Adaptation Process

To classify how a given cell adapts to the drug exposure, we took into account both its viability trajectory and its recorded drug uptake over the last 20 cell cycles. Visually, there were two significantly distinguishable patterns: cells with constant drug uptake, and cells with rapidly increasing viability trajectories (concave shape). Therefore, the following classification criteria were chosen: (i) the amount of absorbed drug is constant; (ii) the viability curve is monotonically increasing over at least 95% of the considered time interval (numerical second derivative of the viability curve is negative); (iii) the remaining cases. As a result, we identified: (i) linear adaptation pattern (constant drug uptake and almost linear viability trajectories); (ii) a superlinear adaptation pattern (concave viability trajectories with diminished drug uptake); (iii) intermediate pattern where cellular uptake was diminishing but the viability trajectory was fluctuating for the majority of time (compare Figures 3B–D and the figure insets).

Cells 3D Spatio-Temporal Routes/3D Lineage Trees

A 3D cell route shows a spatio-temporal evolutionary history of a given tumor cell; that is, it shows all recorded locations of that cell and all cell’s predecessors within the tissue patch. The 3D route is created backwards in time by linking positions (in the x-z plane, locations within a tissue at a given time) of a given cell taken at consecutive time points (y-axis) until the cell’s birth time is reached, and then repeating this procedure for all cell’s predecessors until the beginning of the simulation. The 3D spatio-temporal routes can be traced for multiple cells of the same predecessor forming a 3D spatio-temporal lineage tree (compare Figures 4A, 5A,Bii). These 3D lineage trees are an extension of classical lineage trees used to depict tumor clonal expansion. These figures synthesize information regarding cell locations, cell lineage and time in one single image.

Lineage Trees of Survivors

For each of the initial 65 cells, the whole classical binary lineage tree can be constructed that contains all successors of this cell. A lineage tree of survivors is a subtree of the whole lineage tree and contains only these tree branches that lead to cells that survived the whole treatment (compare Figures 5A,Bi and Supplementary Figures S4–S18).

Precursors of Resistance

The precursors of resistance are these cells for which all successors survived the treatment at the end of simulation. The precursor of resistance is identified by inspecting the lineage tree generated by that cell; if the lineage tree does not contain any dead cells, its initiating cell is considered to be a precursor. We treat the cells that left the computational domain as alive, thus allow them to be successors of the precursor cells (compare Figures 5A,Bi).

Results

We previously analyzed a parameter space of this model and identified regimes for which the whole tumor developed resistance (Gevertz et al., 2015; Perez-Velazquez et al., 2016). Here, we summarize these results briefly. A small colony of 65 sensitive tumor cells was exposed to a drug diffusing from four irregularly placed vessels for the period of about 200 cell cycles. Initially, the tumor increased in size until some cells started responding to drug-induced damage and dying; but the remaining cells finally adapted their tolerance. After about 84 cell cycles, the tumor reached a stable population. The average cell viability showed also a steady increase that confirmed the emergence of a resistant tumor. The evolution of tumor resistance on the population level is presented in Supplementary Figure S1. This showed that a small homogeneous cell colony exposed to a drug gradient can acquire resistance. The final tumor contained the offsprings of 15 initial cells only; the successors of the remaining 50 initial cells went extinct. The rest of the paper is devoted to analysis of resistance at the individual cell level, whether spatial structure of the tumor and tumor microenvironment play a role in the emergence of resistance, and which cell lineages drive resistance of the whole tumor.

Temporal Distributions of Dying Cells Confirm the Drug-Induced Resistance

The fate of each cell depends on both the accumulated damage and the level of damage that the cell can withstand without committing to death. To determine how the resistance is acquired in individual cells, we need to understand the conditions leading to cell death. Initially, each cell has some baseline tolerance level and no accumulated damage. With time, the absorbed drug induces damage to the cell, while the cell can also adapt to the surrounding extracellular conditions that leads to increase in its tolerance. The cell dies when the level of cell damage exceeds the level of cell tolerance to damage. During the simulated treatment, the initially sensitive cells either develop resistance or respond to the treatment and die. In fact, about 75% of the initial 65 cells did not produce offsprings that were able to survive to the end of the treatment period. Since some cells located near the domain boundaries might have been pushed outside of the observed tissue region by the pressure from their growing neighbors, these cells are assumed to move to the other tissue areas and are removed from our system. Here, we only consider cells that remained inside the tissue domain until they were annihilated by the drug. The summary of spatial and temporal analysis of dying cells is shown in Figure 2.

During the initial period of treatment, no cells were dying (no cell counts in histogram in Figure 2A) since they must accumulate the drug-induced damage to overcome the baseline tolerance level. However, the viability trajectories for numerous cells decreased during this time (Figure 2B) confirming that these cells were accumulating damage. Each curve in this graph corresponds to one cell and traces in time the viability values of this cell and all its predecessors, back to one of the initial 65 cells. This period corresponds to a steady tumor growth shown in Supplementary Figures S1i,ii. The first peak in cell death histogram and a time interval when cell viability trajectories reached zero match the significant reduction in the overall tumor size (Supplementary Figure S1iii). The second peak in the death histogram is much smaller since a large number of cells have already developed resistance and only a small subpopulation of cells remained still sensitive to the drug (compare to steady tumor growth in Supplementary Figures S1iv,v). After the time corresponding to about 84 simulated cell cycles, no more cells have died. Similarly, all viability trajectories for these dying cells reached the zero value at or before this time (Figure 2B). This confirms that all remaining tumor cells in the observable tissue patch have developed a drug-induced resistance. Spatially, the tissue regions that are most prone to cell death are situated either near the single vessel in the top-right corner or in the region near the tissue center between the remaining three vessels (Figure 2C). It is worth noting, that in our previous work (Gevertz et al., 2015), we identified the model parameter regimes for which the tumors got extinct, thus the development of drug-induced resistance is not an intrinsic property of our model.

Cell Adaptation to Drug Exposure Can Progress in Three Distinct Ways

To determine how individual cells contributed to the overall tumor resistance, we analyzed the viability trajectories of each cell that survived the treatment (Figure 3A). These graphs confirm our previous observations of several phases in the evolution of resistance in the individual tumor cells: from initial identical viability values, to viability decrease due to the damage being accumulated, to a transient increase in viability when the mechanism of tolerance became activated (initial max), to prolonged reduction in viability values due to accumulated damage approaching the individual cell tolerance level leading to cells adaptation to the drug (prolonged min), to a continuous increase in cell viability when the tolerance mechanism gains a lead. Despite the fact that all surviving cells originated from identical predecessor cells and that they shared very similar viability trajectories for the first 55 cell cycles, we identified three patterns of cell adaptation that resulted in drug-induced resistance (Figures 3B,D).

The first cell subpopulation is characterized by rapid increase in viability values that form concave curves of distinct durations (Figure 3B). In all these cases, there is also a reduced absorption of the drug for a significant length of time (at least 29 cell cycles, inset in Figure 3B). The diminished drug uptake is a result of drug concentration being below the cell’s demand. A closer analysis of cell spatial distributions over time shows that this subpopulation occupied tissue regions distant from the vessels and, more importantly, was surrounded by other cells (green circles in Figures 3Ei–vi). This was a combined effect of cells’ proliferation and their passive relocation due to physical pressure from other growing cells. Since the cells remained in the areas poorly penetrated by the drug for a prolonged time, it resulted in rapidly increasing cell viability that is manifested by the concave shape of the viability curves. The second subpopulation consists of cells with nearly linear increase in viability values (Figure 3C) and with constantly high drug absorption (inset in Figure 3C). The early predecessors were located in between the three central vessels (red circles in Figures 3Eii,iii), and thus were exposed to moderate drug concentrations. This resulted in faster increase of drug-induced tolerance than drug-induced damage and in the steady increase in cell viability values. The fluctuations in almost linear viability patterns arose from competition between gained tolerance, acquired damage and damage repair. This led to repopulation of the space between the blood vessels (Figure 3Eiv). Subsequently, the cells were able to survive in the areas well penetrated by the drug, even in the vicinity of the blood vessels Figures 3Eiv–vi). The remaining resistant cells manifest an intermediate behavior with regards to drug absorption, as it decreases over a very short time near the end of the treatment period (Figure 3D) but not as pronounced as in the subpopulations with concave viability curves. This subpopulation also acquired a quite distinct spatial pattern on a border between two other subpopulations (Figures 3Eiii–vi). These cells are transient in the sense that their characteristics may change during the treatment. For example, a cell with linear viability values may become transient if it gets surrounded by other cells, and becomes protected from drug exposure (cell C₁ in Figure 4Aiii). Similarly, a transient cell can give birth to a cell that falls into the category of concave viability if it moves to the poorly penetrated area (cell C₂ in Figure 4Aiii).

The 3D Cellular Routes Delineate Spatio-Temporal Dynamics of Cell Adaptation

To more closely examine how cells from all three categories can adapt to the treatment, we selected 12 cells (four from each category) with one common predecessor (Figure 4) and traced their locations within the tissue during the whole simulation. The 3D spatio-temporal routes traversed by each cell are shown in Figure 4A at three different time points together with the drug profile at that time. Here, the xz-plane represents cell positions within the tissue, and y-axis represents the time. Note, that the drug distribution profile at each time point is different despite the continuous drug influx from the vessels because drug absorption depends on the total number of cells in the tissue, and this cell number varies in time. The presented exemplary cells were chosen intentionally to show a variety of spatial and temporal dynamics that may lead to cell survival, adaptation and acquired resistance. The corresponding 12 viability trajectories are presented in Figure 4B to confirm characteristics of each cell. Since these cells have a common predecessor, there is a period of time when both the viability trajectories and the 3D routes overlap and thus the number of observable curves is limited (Figures 4A,Bi). However, these curves eventually split up in both figures into eight separate lines (Figures 4A,Bii). Furthermore, individual cells were able to move at significant distances from the position of their common predecessor (Figures 4Aii,C). This was due to the pressure imposed by other growing cells. From this point on, the viability trajectories steadily increased (Figure 4Biii), but cells’ routes deviated only insignificantly forming almost horizontal lines (Figure 4Aiii). This was due to cell overcrowding by numerous neighbors that resulted in cells’ prolonged dormancy, without division. This ultimately contributed toward cell survival and steady increase in cell viability. We intentionally selected a case in which the initial predecessor cell was able to give rise to successors from each of the three categories. However, out of 15 initial cells which successors survived the whole treatment, four generated cells in all three categories, three produced cells in two categories and eight gave rise to cells in a single category.

Lineage Tree Analysis Identifies the Cells That Drive Resistance

Less than a quarter of cells that formed the initial micrometastasis (15 out of 65) produced successors that survived the whole chemotherapeutic treatment. Here we examined the lineages of each subpopulation in order to identify how drug resistance developed for each of them. We inspected the full lineage trees for each of the survived subpopulation and identified subtrees containing only those branches that led from the initial cell to cells that survived the whole treatment. The branches leading to dead cells were omitted. If one of the daughter cells left the domain, but the other survived, its symbol was indicated along the vertical line connecting that cell with its mother cell. These structures represent the lineage trees of survivors. Two representative examples generated by the initial cells with indices #6 and #26 are shown in Figures 5A,Bi. The corresponding spatio-temporal routes traversed by these cells are shown in Figures 5A,Bii. Additional snapshots of spatio-temporal routes at different time points are shown in Supplementary Figures S2, S3. The cell viability trajectories are shown in Figures 5A,Biii. These examples illustrate different cases of cellular adaptation observable among all survived subpopulations. The cells for which viability increases linearly are located in well-penetrated areas. These cells were able to survive the drug insult for a prolonged time since they were surrounded by other cells that absorbed the drug creating a protective niche (Figures 5Aii,iii). Cells with concave viability trajectories are located in poorly penetrated areas, often equidistant from the vessels, where damage induced by the drug is lower that the ability of the cell to repair damage (Figures 5Bii,iii). For some lineage trees of survivors, their spatio-temporal routes may have multiple spatially separated branches due to the proliferation and pressure from neighboring cells. In other cases, the routes do not deviate significantly in space and form horizontal lines. This is due to overcrowding that limits cell proliferation and migration (other 3D routs are discussed in Supplementary Figures S4–S18).

For each lineage tree of survivors, we identified the subtrees that do not contain any dead cells; that is, all branches of these subtrees point either to cells that survived the whole treatment or to cells that left the domain (these cells have positive viability values, so they are alive). The roots of such subtrees are considered to be the precursors of drug resistance, since none of their successors underwent drug-induced death. The precursor cells are indicated by black rectangles and red arrows in the trees shown in Figures 5A,Bi (for clarity, the branches leading to the cells that left the domain are omitted from the graphs). In total, there were 224 precursor cells emerging from all 15 lineage trees of survivors. They all are pictured in Figure 5Ci along the viability trajectories to indicate the time at which they emerge. In Figure 5Cii, these cells are cumulatively projected on the tissue space to show the initial locations of the precursor cells. The cell colors correspond to a time period at which they first appeared. The very first precursor cells have arisen in the area poorly penetrated by the drug between the single vessel in the top-right corner, and the three other vessels (cells shown in green in Figure 5C). Such areas are known as drug sanctuaries or refugia. The next cohort of precursor cells emerged in the center of the tissue between the three blood vessels (indicated by magenta dots in Figure 5C). While, in principle, these areas can be better penetrated by the drug, they actually form protective niches (protectorates) in which the precursor cells may be shielded from the exposure to the drug by the surrounding cells. The final cohort of precursor cells (indicated by cyan dots in Figure 5C) was emerging over a longer period of time and mostly in the areas located closer to the tissue boundaries in the hypoxic or nearly-hypoxic niches. Interestingly, none of the precursor cells were located directly at the concave viability trajectories. This indicates that all precursor cells emerged as a result of a direct competition between drug-induced cell damage and acquired tolerance, and that the increase in cell viability was amplified (in fast superlinear fashion) in cells that have already developed resistance.

Discussion

We presented here a study analyzing how resistant cell lineages arise in micrometastases exposed to a systemic chemotherapeutic treatment. This research is an extension of our previous work (Gevertz et al., 2015; Perez-Velazquez et al., 2016) that focused on the emergence of drug-induced resistance on a cell population level. While we followed the previous mathematical model setup and considered a small tumor growing in a heterogeneous microenvironment, the individual-cell perspective and novel evaluation methods allowed us to identify a new microenvironmental niche prone to the emergence of resistant cells. In addition to previously reported refugia characterized by low drug penetration due to their distance from the vasculature and the hypoxic or near-hypoxic niches in which cells were able to thrive and repair the drug-induced damage, we also located areas in which cells were not exposed to lethal drug concentrations because they were shielded by other cells absorbing the excess of the drug—the protectorates. We also recognized that certain cells gave rise to lineages of resistant cells (precursors of resistance) and correlated three temporal periods with three different spatial locations at which such cells emerged. This supports the hypothesis that tumor micrometastases do not need to harbor cell populations with pre-existing resistance, but that individual tumor cells can adapt and develop resistance induced by the drug during the treatment.

The novel analysis and visualization methods developed here, such as the lineage trees of survivors, the method to identify the precursors of resistance and the 3D sptatio-temporal routes and 3D lineage trees can enhance the library of tools used with other hybrid mathematical models (Kim et al., 2013; Karolak and Rejniak, 2019; Chamseddine and Rejniak, 2020) to analyze tumor evolution and clonality.

Moreover, we showed that once the cells have developed resistance, they were able to elevate their viability either in a fast superlinear manner or in a slower, linear fashion, depending whether they moved toward the refugia areas or not; a small population of transient cells that could transfer from the linear to superlinear populations was also observable. This is in line with the theory of mixed models of tumor evolution (Davis et al., 2017), in which different evolution forms can occur in parallel or can shift from one form to another as a result of changes in tumor size or due to microenvironmental selection forces.

Our results can be also placed within a context of tumor ecology (Kenny et al., 2006; Korolev et al., 2014), such as the ecological concepts of microenvironmental niche partitioning and niche construction (Scott and Marusyk, 2017). In the former case, different cell subpopulations are driven into distinct tissue compartments by the microenvironmental selection forces – we observed that certain cell subpopulations were harbored within the refugia areas or within the hypoxic niches. In the latter, the cells are able to modify their own surroundings to create a favorable microenvironment – we observed the formation of cellular protectorates characterized by microenvironmental conditions distinct from the surrounding areas. This spatial heterogeneity in tumor microenvironments is often referred as ecological habitats (Chang et al., 2017; Sala et al., 2017) that can lead to unique fitness landscapes and selection for different cell phenotypes and genotypes, even under the same extrinsic pressure such as anti-cancer therapy. Our simulations showed that individual cell viability was changing over time that encourages revisiting the idea of a static fitness landscape, and supports the view that cell fitness is not a constant value, but a function of the environmental context (Rozhok and DeGregori, 2015; Scott and Marusyk, 2017). While we did not explicitly model any genetic mutations, the observable changes in tumor cell viability could be potentially linked to changes in cell gene expression.

Ultimately, the link between ecological changes within the tumor microenvironment and tumor evolutionary changes will reflect on patients’ clinical outcome. While the systemic chemotherapy is often used in the clinical protocol in order to minimize the tumor metastatic spread, it should be taken into account that such therapy may stimulate progression of the nearly-killed cells toward resistance. Therefore, the approaches targeting the resistance-inducing strategies may prove more effective than targeting the tumor cells directly. This is similar in concept to eco-evo drugs from the field of microbial antibiotic resistance (Baquero et al., 2011). Some such preconditioning mechanisms have been tested in cancer cells and already showed promise (McDunn and Cobb, 2005; Pisco et al., 2013; Huang, 2014); however, more research in this area is needed.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Author Contributions

The computational model and analysis methods were conceived by KR and JP-V. All computer simulations were performed by KR and JP-V. KR and JP-V contributed to the manuscript writing.

Funding

This work was supported in part by the U01-CA202229 Physical Sciences Oncology Project (PS-OP) grant from the US National Institutes of Health (to KR). JP-V thanks the “Global Challenges for Women in Math Sciences” program of the Mathematics Faculty of the Technical University of Munich for the Entrepreneurial Award to further this collaborative project.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2020.00319/full#supplementary-material

References

Baquero, F., Coque, T. M., and de la Cruz, F. (2011). Ecology and evolution as targets: the need for novel eco-evo drugs and strategies to fight antibiotic resistance. Antimicrob. Agents Chemother. 55, 3649–3660. doi: 10.1128/AAC.00013-11

PubMed Abstract | CrossRef Full Text | Google Scholar

Barcellos-Hoff, M. H., Lyden, D., and Wang, T. C. (2013). The evolution of the cancer niche during multistage carcinogenesis. Nat. Rev. Cancer 13, 511–518. doi: 10.1038/nrc3536

PubMed Abstract | CrossRef Full Text | Google Scholar

Chamseddine, I. M., and Rejniak, K. A. (2020). Hybrid modeling frameworks of tumor development and treatment. Wiley Interdiscip. Rev. Syst. Biol. Med. 12:e1461. doi: 10.1002/wsbm.1461

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, Y. C., Ackerstaff, E., Tschudi, Y., Jimenez, B., Foltz, W., Fisher, C., et al. (2017). Delineation of Tumor Habitats based on Dynamic Contrast Enhanced MRI. Sci. Rep. 7:9746. doi: 10.1038/s41598-017-09932-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheeseman, B. L., Newgreen, D. F., and Landman, K. A. (2014a). Spatial and temporal dynamics of cell generations within an invasion wave: a link to cell lineage tracing. J. Theor. Biol. 363, 344–356. doi: 10.1016/j.jtbi.2014.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheeseman, B. L., Zhang, D., Binder, B. J., Newgreen, D. F., and Landman, K. A. (2014b). Cell lineage tracing in the developing enteric nervous system: superstars revealed by experiment and simulation. J. R. Soc. Interf. 11:20130815. doi: 10.1098/rsif.2013.0815

PubMed Abstract | CrossRef Full Text | Google Scholar

Chisholm, R. H., Lorenzi, T., Lorz, A., Larsen, A. K., de Almeida, L. N., Escargueil, A., et al. (2015). Emergence of drug tolerance in cancer cell populations: an evolutionary outcome of selection, nongenetic instability, and stress-induced adaptation. Cancer Res. 75, 930–939. doi: 10.1158/0008-5472.CAN-14-2103

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, H., and Levy, D. (2017). Modeling the dynamics of heterogeneity of solid tumors in response to chemotherapy. Bull. Math. Biol. 79, 2986–3012. doi: 10.1007/s11538-017-0359-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Collins, J. L., and Kao, M. S. (1989). The anticancer drug, cisplatin, increases the naturally occurring cell-mediated lysis of tumor cells. Cancer Immunol. Immunother. 29, 17–22.

PubMed Abstract | Google Scholar

Correia, A. L., and Bissell, M. J. (2012). The tumor microenvironment is a dominant force in multidrug resistance. Drug Resist. Update 15, 39–49. doi: 10.1016/j.drup.2012.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Cory, T. J., Schacker, T. W., Stevenson, M., and Fletcher, C. V. (2013). Overcoming pharmacologic sanctuaries. Curr. Opin. HIV AIDS 8, 190–195. doi: 10.1097/COH.0b013e32835fc68a

PubMed Abstract | CrossRef Full Text | Google Scholar

Cree, I. A., and Charlton, P. (2017). Molecular chess? Hallmarks of anti-cancer drug resistance. BMC Cancer 17:10. doi: 10.1186/s12885-016-2999-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Dannenberg, J. H., and Berns, A. (2010). Drugging drug resistance. Cell 141, 18–20. doi: 10.1016/j.cell.2010.03.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, A., Gao, R., and Navin, N. (2017). Tumor evolution: linear, branching, neutral or punctuated? Biochim. Biophys. Acta 1867, 151–161. doi: 10.1016/j.bbcan.2017.01.003

CrossRef Full Text | Google Scholar

Feizabadi, M. S. (2017). Modeling multi-mutation and drug resistance: analysis of some case studies. Theor. Biol. Med. Model. 14:6. doi: 10.1186/s12976-017-0052-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferrari, N., Granata, I., Capaia, M., Piccirillo, M., Guarracino, M. R., Vene, R., et al. (2017). Adaptive phenotype drives resistance to androgen deprivation therapy in prostate cancer. Cell Commun. Signal. 15:51. doi: 10.1186/s12964-017-0206-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, F., Nowak, M. A., and Bonhoeffer, S. (2015). Spatial heterogeneity in drug concentrations can facilitate the emergence of resistance to cancer therapy. PLoS Comput. Biol. 11:e1004142. doi: 10.1371/journal.pcbi.1004142

PubMed Abstract | CrossRef Full Text | Google Scholar

Gevertz, J. L., Aminzare, Z., Norton, K. A., Perez-Velazquez, J., Volkening, A., and Rejniak, K. A. (2015). “Emergence of anti-cancer drug resistance exploring the importance of the microenvironmental niche via a spatial model,” in Applications of Dynamial Systems in Biology and Medicine, Vol. 158, eds A. Radunskaya and T. Jackson (Berlin: Springer), 1–34. doi: 10.1007/978-1-4939-2782-1_1

CrossRef Full Text | Google Scholar

Goldman, A. (2016). Tailoring combinatorial cancer therapies to target the origins of adaptive resistance. Mol. Cell Oncol. 3:e1030534. doi: 10.1080/23723556.2015.1030534

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldman, A., Majumder, B., Dhawan, A., Ravi, S., Goldman, D., Kohandel, M., et al. (2015). Temporally sequenced anticancer drugs overcome adaptive resistance by targeting a vulnerable chemotherapy-induced phenotypic transition. Nat. Commun. 6:6139. doi: 10.1038/ncomms7139

PubMed Abstract | CrossRef Full Text | Google Scholar

Greene, J. M., Gevertz, J. L., and Sontag, E. D. (2019). Mathematical approach to differentiate spontaneous and induced evolution to drug resistance during cancer treatment. JCO Clin. Cancer Inform. 3, 1–20. doi: 10.1200/CCI.18.00087

PubMed Abstract | CrossRef Full Text | Google Scholar

Greene, J. M., Sanchez-Tapia, C., and Sontag, E. D. (2018). “Control structures of drug resistance in cancer chemotherapy,” in Proceedings of the IEEE Conference on Decision and Control, Tampa, FL.

Google Scholar

Hafner, M., Niepel, M., Chung, M., and Sorger, P. K. (2016). Growth rate inhibition metrics correct for confounders in measuring sensitivity to cancer drugs. Nat. Methods 13, 521–527. doi: 10.1038/nmeth.3853

PubMed Abstract | CrossRef Full Text | Google Scholar

Hambardzumyan, D., and Bergers, G. (2015). Glioblastoma: defining tumor niches. Trends Cancer 1, 252–265. doi: 10.1016/j.trecan.2015.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Holohan, C., Van Schaeybroeck, S., Longley, D. B., and Johnston, P. G. (2013). Cancer drug resistance: an evolving paradigm. Nat. Rev. Cancer 13, 714–726. doi: 10.1038/nrc3599

CrossRef Full Text | Google Scholar

Housman, G., Byler, S., Heerboth, S., Lapinska, K., Longacre, M., Snyder, N., et al. (2014). Drug resistance in cancer: an overview. Cancers 6, 1769–1792. doi: 10.3390/cancers6031769

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, S. (2014). The war on cancer: lessons from the war on terror. Front. Oncol. 4:293. doi: 10.3389/fonc.2014.00293

CrossRef Full Text | Google Scholar

Huch, M., and Rawlins, E. L. (2017). Cancer: tumours build their niche. Nature 545, 292–293. doi: 10.1038/nature22494

CrossRef Full Text | Google Scholar

Hutchinson, L. (2016). Genetics: defining driver mutations in the genomic landscape of breast cancer. Nat. Rev. Clin. Oncol. 13:327. doi: 10.1038/nrclinonc.2016.75

CrossRef Full Text | Google Scholar

Karolak, A., and Rejniak, K. A. (2019). Micropharmacology: an in silico approach for assessing drug efficacy within a tumor tissue. Bull. Math. Biol. 81, 3623–3641. doi: 10.1007/s11538-018-0402-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kenny, P. A., Nelson, C. M., and Bissell, M. J. (2006). The Ecology of Tumors: by perturbing the microenvironment, wounds and infection may be key to tumor development. Scientist 20:30.

PubMed Abstract | Google Scholar

Kim, M., Gillies, R. J., and Rejniak, K. A. (2013). Current advances in mathematical modeling of anti-cancer drug penetration into tumor tissues. Front. Oncol. 3:278. doi: 10.3389/fonc.2013.00278

PubMed Abstract | CrossRef Full Text | Google Scholar

Korolev, K., Xavier, J., and Gore, J. (2014). Turning ecology and evolution against cancer. Nat. Rev. Cancer 14, 371–380. doi: 10.1038/nrc3712

PubMed Abstract | CrossRef Full Text | Google Scholar

McDunn, J. E., and Cobb, J. P. (2005). That which does not kill you makes you stronger: a molecular mechanism for preconditioning. Sci STKE. 2005:e34.

PubMed Abstract | Google Scholar

Mehrara, E., Forssell-Aronsson, E., Ahlman, H., and Bernhardt, P. (2007). Specific growth rate versus doubling time for quantitative characterization of tumor growth rate. Cancer Res. 67, 3970–3975. doi: 10.1158/0008-5472.can-06-3822

PubMed Abstract | CrossRef Full Text | Google Scholar

Navin, N. E., and Hicks, J. (2010). Tracing the tumor lineage. Mol. Oncol. 4, 267–283. doi: 10.1016/j.molonc.2010.04.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikbakht, H., Panditharatna, E., Mikael, L. G., Li, R., Gayden, T., Osmond, M., et al. (2016). Spatial and temporal homogeneity of driver mutations in diffuse intrinsic pontine glioma. Nat. Commun. 7:11185. doi: 10.1038/ncomms11185

PubMed Abstract | CrossRef Full Text | Google Scholar

Perez-Velazquez, J., Gevertz, J. L., Karolak, A., and Rejniak, K. A. (2016). Microenvironmental niches and sanctuaries: a route to acquired resistance. Adv. Exp. Med. Biol. 936, 149–164. doi: 10.1007/978-3-319-42023-3_8

PubMed Abstract | CrossRef Full Text | Google Scholar

Perez-Velazquez, J., Quinones, B., Hense, B. A., and Kuttler, C. (2015). A mathematical model to investigate quorum sensing regulation and its heterogeneity in pseudomonas syringae on leaves. Ecol. Compl. 21, 128–141. doi: 10.1016/j.ecocom.2014.12.003

CrossRef Full Text | Google Scholar

Pisco, A. O., Brock, A., Zhou, J., Moor, A., Mojtahedi, M., Jackson, D., et al. (2013). Non-Darwinian dynamics in therapy-induced cancer drug resistance. Nat. Commun. 4:2467. doi: 10.1038/ncomms3467

PubMed Abstract | CrossRef Full Text | Google Scholar

Pisco, A. O., and Huang, S. (2015). Non-genetic cancer cell plasticity and therapy-induced stemness in tumour relapse: ‘What does not kill me strengthens me’. Br. J. Cancer 112, 1725–1732. doi: 10.1038/bjc.2015.146

PubMed Abstract | CrossRef Full Text | Google Scholar

Puhalla, S., Elmquist, W., Freyer, D., Kleinberg, L., Adkins, C., Lockman, P., et al. (2015). Unsanctifying the sanctuary: challenges and opportunities with brain metastases. Neuro Oncol. 17, 639–651. doi: 10.1093/neuonc/nov023

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, G. Z., Jin, M. Z., Dai, J. X., Sun, W., Feng, J. H., and Jin, W. L. (2017). Reprogramming of the tumor in the hypoxic niche: the emerging concept and associated therapeutic strategies. Trends Pharmacol. Sci. 38, 669–686. doi: 10.1016/j.tips.2017.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Rozhok, A. I., and DeGregori, J. (2015). Toward an evolutionary model of cancer: considering the mechanisms that govern the fate of somatic mutations. Proc. Natl. Acad. Sci. U.S.A. 112, 8914–8921. doi: 10.1073/pnas.1501713112

PubMed Abstract | CrossRef Full Text | Google Scholar

Sala, E., Mema, E., Himoto, Y., Veeraraghavan, H., Brenton, J. D., Snyder, A., et al. (2017). Unravelling tumour heterogeneity using next-generation imaging: radiomics, radiogenomics, and habitat imaging. Clin. Radiol. 72, 3–10. doi: 10.1016/j.crad.2016.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmidt, M. M., and Wittrup, K. D. (2009). A modeling analysis of the effects of molecular size and binding affinity on tumor targeting. Mol. Cancer Ther. 8, 2861–2871. doi: 10.1158/1535-7163.MCT-09-0195

PubMed Abstract | CrossRef Full Text | Google Scholar

Scott, J., and Marusyk, A. (2017). Somatic clonal evolution: a selection-centric perspective. Biochim. Biophys. Acta Rev. Cancer 1867, 139–150. doi: 10.1016/j.bbcan.2017.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, A. B., Rejniak, K. A., and Gevertz, J. L. (2016). Limiting the development of anti-cancer drug resistance in a spatial model of micrometastases. Math. Biosci. Eng. 13, 1185–1206. doi: 10.3934/mbe.2016038

CrossRef Full Text | Google Scholar

Sharma, S. V., Lee, D. Y., Li, B., Quinlan, M. P., Takahashi, F., Maheswaran, S., et al. (2010). A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell 141, 69–80.

Google Scholar

Su, Y., Wei, W., Robert, L., Xue, M., Tsoi, J., Garcia-Diaz, A., et al. (2017). Single-cell analysis resolves the cell state transition and signaling dynamics associated with melanoma drug-induced resistance. Proc. Natl. Acad. Sci. U.S.A. 114, 13679–13684. doi: 10.1073/pnas.1712064115

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y. (2016). Tumor microenvironment and cancer therapy resistance. Cancer Lett. 380, 205–215. doi: 10.1016/j.canlet.2015.07.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Swift, L. H., and Golsteyn, R. M. (2014). Genotoxic anti-cancer agents and their relationship to DNA damage, mitosis, and checkpoint adaptation in proliferating cancer cells. Int. J. Mol. Sci. 15, 3403–3431. doi: 10.3390/ijms15033403

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, A., Loutherback, K., Lambert, G., Estevez-Salmeron, L., Tlsty, T. D., Austin, R. H., et al. (2013). Cell motility and drug gradients in the emergence of resistance to chemotherapy. Proc. Natl. Acad. Sci. U.S.A. 110, 16103–16108. doi: 10.1073/pnas.1314385110

PubMed Abstract | CrossRef Full Text | Google Scholar