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 druglimited 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 nonhomogeneous 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 weaklyproliferative 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. (2018Greene et al. ( , 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, druginduced 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 nongenetic 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 . 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: uptake by the cells uptake by the cells 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: 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 mat i . 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 θ: 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: 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 dam i (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 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). rate p γ (three times faster that the damage rate (Gevertz et al., 2015;Perez-Velazquez et al., 2016): 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 tol i (t) (at a slow rate of tol = 0.01% of the baseline tolerance value; Gevertz et al., 2015;Perez-Velazquez et al., 2016): where the amount of accumulated drug C γ i (t) depends on its continuous absorption (at a constant rate ρ γ ) and decay (at a rate d γ ): 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 dam i (t) exceeds the tolerance level C tol i (t). 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).

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 topright 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 1 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 2 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. 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 topright 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 nearlyhypoxic 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 individualcell 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.