Abstract

The development of resistance is commonly attributed to genetic heterogeneity. However, heterogeneity of drug penetration of the tumor microenvironment both on the microscopic level within solid tumors as well as on the macroscopic level across metastases may also contribute to acquired drug resistance. Here we use mathematical models to investigate the effect of drug heterogeneity on the probability of escape from treatment and the time to resistance. Specifically we address scenarios with sufficiently potent therapies that suppress growth of all preexisting genetic variants in the compartment with the highest possible drug concentration. To study the joint effect of drug heterogeneity, growth rate, and evolution of resistance, we analyze a multi-type stochastic branching process describing growth of cancer cells in multiple compartments with different drug concentrations and limited migration between compartments. We show that resistance is likely to arise first in the sanctuary compartment with poor drug penetrations and from there populate non-sanctuary compartments with high drug concentrations. Moreover, we show that only below a threshold rate of cell migration does spatial heterogeneity accelerate resistance evolution, othenNise deterring drug resistance with excessively high migration rates. Our results provide new insights into understanding why cancers tend to quickly become resistant, and that cell migration and the presence of sanctuary sites with little drug exposure are essential to this end.

Author Summary

Here, we take into consideration important aspects of the tumor microenvironment, including spatial drug gradients and differential rates of cell proliferation. We show that the dependence of fitness on space together with cell migration facilitates the emergence of acquired resistance. Our analysis indicates that resistant cells that are selected for in compartments with high concentrations are likely to disseminate from sanctuary sites Where they first acquire resistance preceding migration. The results suggest that it would be helpful to improve clinical outcomes by combining targeted therapy With anti-metastatic treatment aimed at constraining cell motility as well as by enhancing drug transportation and distribution throughout all metastatic compartments.

Introduction

Targeted cancer therapy is currently an area of active research [6] , and is under rapid development [7—10]. Targeted agents can send cancer into remission, but the response is often short-lived [7—11]. The recurrence of cancer, if treated with single agents, is almost certain due to acquired drug resistance [12, 13]. Among efforts to understand the rapid acquisition of resistance by cancer cells, particular attention has been paid to the preexisting resistance arising prior to treatment [14—17].

Important aspects of the tumor microenvironment include spatial drug gradients and differential rates of cell proliferation [22, 23]. The heterogeneity of this kind can be found both on the microscopic level within solid tumors as well as on the macroscopic level across metastases [20, 22, 23]. It is not uncommon that incomplete drug distributions within and across metastatic lesions can compromise the efficacy of treatment [22, 23], which may be in part due to the differences in the ability of the drug to penetrate different tissues [24]. As shown in recent experimental and theoretical studies [25—28], drug gradients can help accelerate the evolution of antibiotic resistance. Moreover, a most recent overview of clinical and pharmacological data concerning distribution of many anticancer drugs in human solid tumors highlights the likely importance of insufficient and/ or heterogeneous exposure of cancer cells to effective drug levels in tumor resistance [29]. Thus, in order to improve efficacy of cancer therapy, it is highly relevant to investigate how heterogeneous levels of drug distribution in different parts of the tumor or across metastases affect the emergence of resistance.

By activating tissue invasion and metastases [36] , cancer cells are able to escape from the primary site and disseminate to distant parts of the body, causing life-threatening health problems [37, 38]. At the time of diagnosis and treatment (which generally occur late in the course of disease), a high proportion of common cancer patients have already had tumor cells disseminated to distant sites for years prior to presentation [30, 39]. What is more, clinical outcomes are often complicated by the presence of overt or occult micrometastases in patients [30, 39—42]. It is, therefore, of primary interest to understand the emergence of resistance, particularly in the setting of disseminated cancer.

This observation suggests that disseminated cancer, rather than only the primary tumor in situ, is actually under stress when treatment is started [30, 39, 41]. Moreover, a few studies observed that metastatic cells tend to quickly become chemoresistant [21, 44—49] , suggesting a positive relationship between metastatic phenotype and the rapid acquisition of drug resistance. It is possible that, with the migration and seeding dynamics [31, 34, 39, 50—53], these metastatic cells located at distant metastatic compartments can reseed each other, particularly in the presence of tumor sanctuary sites with very little drug exposure. Therefore, it is necessary to explicitly account for the role that the metapopulation structure of metastatic disease plays in the rapid emergence of resistance.

This important experimental result begs theoretical questions aimed at revealing the relevant pathways of evolving resistance. Do tumor cells first migrate from low concentration areas and then adapt to high concentration areas? Or alternatively, do they first acquire resistance in low concentration areas and then migrate to and populate areas of exceedingly high concentrations? How do cell motility and drug gradient, or more generally how does the spatial heterogeneity in drug concentrations affect the emergence of resistance? The goal of the present work is to answer these questions and provide qualitative insights by a simple conceptual model.

Spatial compartments mean different locations that harbor metastatic deposits (i.e., target organs, such as bone marrow, liver, brain and lung) [34]; migration means dissemination and seeding of cells between metastatic compartments [35]. We focus on the common population dynamics that govern the evolution of resistance for various cancers differing in their capacity to metastasize. Metastases of solid tumors (such as breast cancer [42] and melanoma [10]) tend to have well-defined spatial compartments because of low dissemination rate, whereas the compartment structure of liquid cancer (e.g. blood tumors [7]) is diminished by exceedingly high fluidity. By adjusting the migration rate, our model can be suited to study the specific kind of cancer in question.

In line with this, we use a stochastic, multi-type branching model to account for the fate of individual cells, particularly these drug-resistant mutations in establishing surviving lineages. Mathematical models of this kind have provided particularly useful insights into understanding evolutionary dynamics of cancer in response to treatment [11, 14—17, 63—66] (see a review in Ref. [67] ). A large set of previous models are focused on preexisting resistance in the primary tumor, arising from neutral evolution prior to treatment [15, 17]. Built on these prior studies, the current work incorporates the compartment structure of metastatic disease and quantifies the role spatial heterogeneity in drug concentrations plays in the evolution of resistance by metastatic cells during treatment.

Spatial heterogeneity in drug concentrations means that there exist sanctuary sites that are not or only partially penetrated by drugs. Therefore, reproductive fitnesses of cells depend not only on their cell types but also on their spatial locations, leading to a rugged fitness landscape (since drug concentrations in different spatial compartments are not necessarily continuous but discrete in space). As shown in previous studies [11, 14, 15] , if the fitness of resistant cells as compared to sensitive cells is neutral or even slightly advantageous in the absence of drugs, the acquisition of resistance not only becomes highly likely but also is accelerated, since there is no selection pressure against resistance. However, it is less clear about the most likely pathway to select for resistance, if resistance mutations incur a fitness cost in the absence of drugs while conferring an advantage over sensitive cells in the presence of drugs. Recent mathematical modeling with laboratory test using mice suggests that resistance carries a fitness cost [68]. In View of this, the present study is focused on the latter scenario with fitness cost of resistance, although our approach works for any fitness landscape. We show that resistant mutants are unlikely to emerge in situ in compartments of high drug concentrations, but arise through the mutation-migration pathway; namely, metastatic cells acquire costly resistance in the sanctuary sites preceding migrating to and thriving in harsh compartments containing high levels of drugs.

Results

We consider two spatial compartments with dichotomic distributions of drugs: compartment ‘0’ can hardly be penetrated by the cancer drug, thereby representing a perfect drug sanctuary site; compartment ‘1’ is distributed with an adequate amount of drugs that is able to completely wipe out any wild type cells. We assume that one point mutation is sufficient to confer high levels of resistance to the maximum possible concentration of drugs administered during therapy. We denote the genotypes of cells by the number of acquired point mutations: the wild type ‘0’ and the resistant type ‘1’ (see Materials & Methods for a detailed description of the model). This minimal model, albeit overly simplified, offers intuitive insights into understanding how spatial drug heterogeneity facilitates the acquisition of de novo resistance during cancer therapies. Later on, we will extend this simple model to more realistic cases with multiple cell types and with multiple compartments, where multiple point mutations can successively accrue to confer resistance to exceedingly high drug concentrations [54].

Without teasing out the effect of population size, we cannot clearly pinpoint the role of spatial drug heterogeneity in tumor resistance. Thus, let us start with our theoretical analysis with tracking the lineages derived from a single sensitive cell that can be initially placed in either compartment. Nevertheless, using the multiplicative properties of branching processes, we can calculate the probability of escape and the (conditional) average time to resistance for any given initial conditions of tumor size or metastases (see derivation details in 81 Text).

In most circumstances, it is the sanctuary compartment that persistently seeds the drug-containing compartment, which would have almost certainly become void otherwise due to effective levels of potent drugs. For this reason, in the presence of a sanctuary compartment, there exist two competing pathways to lead to the outgrowth of resistance in the compartment of high drug concentrations, as depicted in Fig. 1B and 1C. One is the “migration-mutation” pathway: sensitive cells first emigrate from the sanctuary compartment and then adapt in situ to the non-sanctuary compartment with high drug concentrations. The other is the “mutation-migration” pathway: sensitive cells first mutate and acquire resistance in the sanctuary compartment, and then migrate to and populate the compartment with high drug concentrations. To understand how the presence of tumor sanctuary sites and cell migration together affect the resistance evolution, we need to determine which pathway provides the more likely path to resistance. Because cells of the same types have different fitness in the two compartments, we regard migration as a sort of status change in spatial locations. In this way, the two competing pathways can be seen as three-type branching processes, respectively, with different fitness

The condition above can be greatly simplified when the two compartments provide exactly the same condition for population outgrowth in the absence of drugs. The rugged fitness landscape due to spatial drug heterogeneity, therefore, can be reflected solely by differences in proliferation rates (i.e., doo = dlo 2 dm 2 du). Moreover, the fitness cost of resistance, 5, for resistant cells located in the sanctuary compartment 0 can be parameterized as 1910 = (1—5)b00, while the fitness cost of sensitivity, 6, for wild type cells in the drug-containing compartment 1, 1901 = (1—6)b00. Substituting these parameterizations into (1), we arrive at the much simplified condition in terms of s and 6: We immediately observe that this condition is expected to be fulfilled in most cases since potent therapy should be characterized by 6 > 5.

2A). We also calculate and compare the average time to resistance (relapse time conditional on non-extinction) following each pathway, as shown in Fig. 2B. The result demonstrates that the simple condition as given in Eq. (2) works well in the limit of low migration rate and remains a good approximation for intermediate migration rates. Moreover, if resistant mutation is neutral or even advantageous in the sanctuary, then the mutation-migration pathway is always much faster than the migration-mutation pathway to result in resistance in the drug-containing compartment (Fig. 2C).

We thus use the parameter 0<gb g 1 to rescale the proliferation and death rates of cells in the sanctuary compartment 0 relative to that in the drug-containing compartment 1: d00/¢ = d10/ gb 2 dm and 1901 = (1—6) (900/ gb (which implies that in the absence of treatment, cells grow 1/ ([3 times as fast in compartment 1 as when located in compartment 0). Although one should refer to the general inequality (1) as the exact condition for the mutation-migration pathway to be predominant, simple algebra shows that the simplified inequality (2) is still a necessary condition in this case. In fact, reducing the rate of cell proliferation in the sanctuary equally affects the two pathways, delaying the time to resistance (Fig. 2D). Taken togethers, these results demonstrate that, under a wide variety of conditions (including the ranges of parameter values relevant to cancer), prevailing resistant cells in compartments with high levels of drugs are likely to originate from sanctuary sites, where they acquired resistance preceding migration.

We develop a numerical method to show the spatio-temporal dynamics of emerging drug resistance across spatial compartments (see Materials & Methods). Note that different from the constrained pathways analyzed in Fig. 2, in this case, migration is allowed to be bidirectional, and both pathways can be at work at the same time. Doing this enables us to study the emergence of resistance in a more natural and realistic setting. Fig. 3 shows the joint probability distribution of the numbers of resistant cells in the two compartments with respect to time. The skewed distribution in Fig. 3A suggests that the sanctuary compartment provides much more favorable condition to evolve resistance than the drug-containing compartment, and thus escaping from the drug-containing compartment to the sanctuary compartment is crucial to this end. Constantly seeding the drug-containing compartment with evolved resistant cells tends to make the distribution more balanced (cf Fig. 3B and 3C). As a result, resistance soon gets established in the drug-containing compartment, and its growth outpaces that in the sanctuary compartment (Fig. 3D), although the chance of having resistance increases steadily with time in both compartments. Similar results are obtained using different initial conditions, despite that resistance evolution is more likely (and sooner) to occur When the sensitive cell is initially placed in the sanctuary compartment than in the drug-containing compartment (cf. 81 Fig. and Fig. 3). Our results demonstrate that the sanctuary compartment, even though cells are slowly replicating therein, serves as an escape hatch, and indeed is most likely to be the breeding ground of resistance. Therefore, in the presence of sanctuary sites, the overwhelming outgrowth of resistance in the drug-contain-ing compartment is an inevitable outcome due to the mutation-migration pathway.

1—3), we now turn to predict outcomes of hypothetical treatments to eradicate (two) metastases in cancer patients (Fig. 4). We assume the total mass of metastasis is relatively small so that treatment starts without any preexisting resistance. In particular, we assume that the two metastatic lesions differ in size as well as in the level of drug penetration, because of different microenvironments. Specifically, cells in lesion 0 grow much more slowly than in lesion 1 in the absence of drugs, yet the drugs have better penetration of lesion 1 than of lesion 0. To this end, we continuously vary the difference of drug concentrations in the two compartments, AD 2 D1—D0, while keeping the total sum of concentrations constant. Let us now specifically incorporate into the simple model a Hill function that describes concentration-dependent killing efficacy of drugs (see Materials & Methods).

However, a worrying situation arises in the presence of a sanctuary with sufficiently low drug concentrations. With increasing difference of drug concentration between the two metastatic compartments, treatment responses exhibit a sharp transition from successful eradication to failed treatment due to acquired resistance (Fig. 4A). Our mathematical framework allows us to calculate the relapse curve, 1—p5(t), where ps(t) is the probability of no resistance by time 1‘ following treatment, and the (conditional) average time to patient relapse due to acquired resistance (see Materials & Methods). Fig. 4B shows that the relapse is accelerated by the spatial heterogeneity; the larger AD, the faster the relapse. In this simulated hypothetical patient with small lesions, the worst-case scenario is when the drug cannot penetrate lesion 0 at all, and thus the relapse occurs on average approximately 103 days. The relapse could have occurred within weeks if the treatment starts with much bigger lesions (82 Fig. ).

Noteworthy, cell lineages originating from lesion 1 are faster to evolve resistance than these originating from lesion 0 (Fig. 4B). This result is mainly due to the initial condition used: lesion 1 is much larger than lesion 0 that the influx of escaping cells to the sanctuary exceeds the number of cells in situ. Since we are considering a branching process, larger population size is more likely to generate resistant mutation. Indeed, as shown in 82 Fig., if the size of lesion 1 is smaller, cells lineages originating from lesion 1 are actually slower to evolve resistance than from lesion 0.

Under different assumptions of fitness effects of mutations as described in Eq. (3), however, multiple point mutations might be required to confer sufficient level of resistance to increasingly high concentrations in compartment 1, as AD increases. This variation does not change the general picture about how the presence of sanctuary sites impairs the effectiveness of cancer therapies (83A Fig. ), but the time to the sufficient levels of resistance may well depend on how many point mutations are needed to this end and exhibits an abrupt increase when AD is increased beyond a critical threshold value (the vertical line in 83B Fig.). For very large AD, only two-point mutants are able to survive in compartment 1, wherein the abundances of one-point mutants and sensitive cells are maintained by the ‘immigration-death’ dynamics of the branching process. As time passes by, two-point mutants eventually pop up, most likely in the sanctuary compartment 0, and subsequently immigrate to and populate the compartment 1.

Furthermore, only for migration rates below a certain critical threshold does the spatial heterogeneity in drug concentrations speed up the emergence of resistance. Excessively high migration rates actually slow down resistance emergence (S4 F ig. ), as the role of compartment structure is diminished by frequent migrations and consequently cells are exposed to the non-sanctuary compartment more often. It is worth noting that there exists an optimum migration rate that leads to the fastest emergence of resistance and that excessively high migration rates actually deter and delay the evolution of resistance (S4 Fig.). The results suggest that it may be helpful to improve clinical outcomes by combining targeted therapy with anti-metastatic treatment aimed at inhibiting cell motility as well as by enhancing drug transportation and distribution throughout all metastatic compartments. We also extend the simple model to more general cases with multiple cell types and with multiple compartments (see detailed mathematical description in 81 Text). In this extended model, multiple point mutations can consecutively accrue to confer resistance to increasingly high drug concentrations. In particular, we study and compare the impact of the two different schemes of cell migration, local versus global migration, on the evolution of resistance (S5 and S6 Fig.). The numerical results confirm that our conclusions derived from the simple model above remain qualitatively unchanged. Additionally, we observe that strong resistance evolves much faster for local migration than global migration, particularly for sensitive cells located at the sanctuary sites with low drug concentrations (S5 F ig.). In other words, sequential local migration over spatial gradient of drug concentration following newly accrued point mutations (i. e., the mutation-migration pathway) plays an important part in leading to rapid selection for high-level resistance. This may be of particular relevance for the evolution of resistance within a solid tumor [22, 54]. Furthermore, apart from oversimplified migration schemes addressed in the current study, this full model can be readily extendedcenter to integrate with a realistic vascular network that regulates the metastatic routes of circulating tumor cells among target organs [53].

Specifically, we consider a solid tumor consisting of 1011 cells which is about 4.6 cm in diameter [57]. Fig. 5A shows the schematic representation of tumor microenvironment of cells surrounding a blood vessel located in the center: the rate of proliferation of tumor cells decreases with increasing distance from the central blood vessel in the absence of treatment. Similarly, the delivery of cancer drugs is also compromised in the presence of treatment, thereby resulting in the spatial drug gradient as illustrated in Fig. 5B. The level of spatial drug heterogeneity is represented by 1/ TD: the larger of this value, the more poorly the drug penetrates distal tumor issues away from the nearest blood vessel. To make progress in our calculations, we artificially divide the tumor into M = 30 spatial compartments with consecutive concentric circles with equal interval in between from center to surface. We confirm that dividing more compartments leads to almost the same results as shown here in Fig. 5C and 5D. In this simulated example, two point mutations are required to confer full resistance to the maximum drug concentration in the center. As shown in Fig. 5C, the tumor can be eradicated under perfect drug penetration where the drug is almost homogeneously distributed throughout the entire tumor population (extremely small 1/ TD). In contrast, inadequate penetration of the tumor gives rise to sanctuary sites, these outer compartments that are most distant away from the blood vessel and thus exposed with the least amount of drugs. Therefore, cancer therapy fails certainly with large values of 1/ TD. In line with Fig. 4, relapse occurs sooner with increasing spatial drug heterogeneity, 1/ TD (Fig. 5D). Moreover, distal cells, although slowly proliferating, are more likely to generate resistance than these proximal cells that are affected most by the drug. These results quantitatively demonstrate that tumor microenvironment mediates cell sensitivity to drugs and thus plays an important role in drug resistance acquired during treatments [22, 29, 54].

Discussion

Cancer cells can migrate from one spatial compartment to another. As compartments may contain different levels of drugs, the cells experience distinct selection pressure for resistance in different compartments. We calculate the probability of resistance and the average time to resistance for sensitive cells originally located in different compartments. We show that the presence of sanctuary sites with poor drug penetration can speed up the emergence of acquired resistance to cancer therapy. Moreover, we show that resistance is unlikely to arise in situ within high concentration compartments, but rather that resistance first emerges in sanctuary sites and then spreads to and populates other compartments with high concentrations that would be able to completely inhibit growth of sensitive cells. This result is in line with a prior analysis of competing pathways to resistance, based on an ecological source-sink model with logistic population growth [76]. In spite of that, our study is specifically aimed at understanding the role of tumor micro-environment in the evolution of resistance during treatment.

Although well suited for studying experimental microbial evolution [26, 27], these models are not directly applicable to the context of cancer; especially solid tumours with intrinsic heterogeneity in their microenvironments warrant a thorough separate investigation [18—20]. Furthermore, in these prior models [26, 27], compartments are placed in an order with increasing drug concentration, sequential migration occurs only between the two nearest neighbor compartments, and evolution begins within the sanctuary while other compartments are initially void. In contrast, the present mathematical framework takes into account realistic concentration-dependent response for any initial population of cancer cells distributed over multiple metastatic compartments, and allows us to calculate the risk of acquiring resistance as well as to ascertain the timing of relapse. We note that it is promising for future studies to improve mathematical tractability of this problem addressed in this work, for example, by using various analytical techniques as detailed in Refs. [69, 70]. Moreover, our model can be readily extended to incorporate specific migration schemes of cancer cells, e. g., along with a more realistically connected vasculature [53]. In parallel, it is worth mentioning that the approach of using partial differential equations to describe spatio-temporal dynamics of cancer evolution [77—79] sheds a different yet useful light on the selection of resistance under cancer therapies [80, 81]. Taken together, our theoretical results improve our understanding of how metastatic cells can acquire resistance during treatment, especially in the presence of sanctuary sites.

It has been shown that migratory cells, though with lower growth potential than nonmigratory cells, can be selected for during therapy [82]. Extending this prior result, our results show that cell motility and the presence of sanctuary sites with little drug exposure are essential for the rapid acquisition of resistance by metastatic cells. Moreover, only for low migration rates below a certain threshold does spatial heterogeneity in drug concentration speed up resistance evolution. Arguably, this finding may help to explain qualitatively the differences of clinical successes in treating liquid cancer (such as chronic myeloid leukemia [7]) and solid tumors (such as melanoma [10]). Because of limited cell motility, metastases of solid tumor tend to have well-defined spatial compartments (namely, distal lesions), and thus are more prone to drug penetration problems. For this reason, additional attention should be paid to eliminate the sanctuary sites for cancer therapy [21].

It is promising for future work to study multi-drug resistance requiring multiple sets of resistant mutations, given that combination therapy is increasingly used in clinical setting [83, 84]. (We refer to Ref. [85] for a most recent development in this matter with a focus on multi-drug resistance to antiviral combination treatments.) In addition to the spatial heterogeneity in drug concentrations addressed here, we think that epi-static interactions of resistant mutations to each drug may also be important and deserve further investigation [86]. In general, our work provides a mathematical and computational framework for studying how various aspects of tumor microenvironment, and spatial drug heterogeneity in particular, influence cell sensitivity to drugs and thus mediate the evolution of resistance during treatment. More important, our results may allow us to better understand the recent in vitro experiment that shows the rapid emergence of high-level resistance can be facilitated by the presence of spatial drug gradients and the motility of metastatic cancer cells along the gradient [54]. In this work, we demonstrate that the mutation-migration pathway is most relevant for establishing sufficient levels of resistance in areas of high drug concentrations that would have wiped out all initial sensitive cells. This theoretical result can stimulate further empirical research that aims at better understanding the spatio-temporal evolutionary dynamics of resistance (e. g., the competing pathways of resistance as shown in Fig. 1B and 1C) using single-cell tracking and sequencing technologies.

This theoretical result can be extended to derive potential therapeutic strategies for increasing the efficacy of cancer therapy. Specifically, our results suggest that combining targeted therapy with anti-metastatic treatment might help improve clinical outcomes [87—89] , especially when treating disseminated cancer. As demonstrated in the present work, inhibition of cell migration between compartments not only suppresses the escape route of sensitive metastatic cells to sanctuary sites and also prevents the dissemination of evolved metastatic cells from sanctuary sites to these compartments with high concentrations, where resistance is strongly selected for. In addition to targeting cell dissemination, it is desirable to enhance drug transportation and distribution throughout all metastatic compartments, in order to deter the rise of resistance [29]. The current model is minimalistic, but allows proof of principle. We leave out many important issues, such as cancer stem cells [90—93] , cellular quiescence and cancer dormancy [94, 95] , and the inefficiency of metastatic processes [96—98]. In particular, we have considered the dynamics of the evolution of drug resistance in a situation in which cells are able to move freely from one compartment (tumor) to another. This first approximation, while enlightening, may overestimate some of the dynamics that would occur in a more realistically connected vascula-ture, especially when the inefficiencies of metastasis, due to filtration and dissemination in the vascular network, are considered [52, 53]. With the increasing understanding of the molecular biology of metastasis as well as clinical advances in treating metastasis, we believe that it will become feasible to obtain accurate estimations of key parameters regarding metastatic burden (location and size) and metastatic rates. Then a calibrated model of this sort as introduced here can be used to simulate patient responses in silica and predict outcomes of treatments to eradicate the disseminated cancer [11], as well as to derive more efficacious treatment strategies, particularly for overcoming the problem of imperfect drug penetrations [29].

Materials and Methods

Minimal model

In our model, we explicitly account for the compartment structure of tumor microenVironment as well as the spatial heterogeneity in drug concentrations across compartments. For proof of concept, we focus on the simplest possible case with only two compartments and two types of cells in the main text. Without loss of generality, we assume that drugs have better access to compartment 1 than to compartment 0. In contrast to conventional chemotherapy agents that have cytotoxic effects, most molecularly targeted cancer therapies have cytostatic effects on cancer cells [75]. Moreover, it is commonly found that the efficacy of drugs is concentration-dependent in pharmacological kinetics studies ranging from antimicrobial treatment to cancer therapy [99—102]. Therefore, it is plausible for us to specifically consider concentration-dependent inhibition of cell replication in response to cancer therapy.

Upon division, one of the daughter cells can mutate with probability u to become resistant. Denote by i the genotype of a cancer cell if it has acquired 1' point mutations (i = 0, 1). The fitness of a cell depends on its type and spatial location. Specifically, the replication rate of a cancer cell, 19 type i and spatial location j = 0, 1) as follows, ij, is determined by its geno-Here we use a Hill function for the drug response curve [103]. fij (05], respectively) is the division rate (the death rate, respectively) of a sensitive cell located in compartment j in the absence of drugs, 5 is the cost of resistance per point mutation in the absence of drugs, ICSO is the drug concentration that is needed to inhibit cell growth by one half of its original rate, p is the fold increase in ICSO per mutation, and m determines the steepness of the Hill function. Previous studies have fitted empirical data to similar Hill functions as given above [Eq. (3)], to indicate antibiotic resistance and antiviral resistance to various treatment regimens [99, 103]. Although a full characterization of cancer drug resistance in this way has yet to be done, our general results are not dependent on specific parameter choices of s, p, m and IC50. The death rate of a cell with genotype 1' and in compartment j is unaffected by the presence of drugs and equal to that of sensitive cells irrespective of their genotypes, dij = aj. The net growth rate is denoted by rij = b,-j—d,-j. Cells can migrate between the two compartments with rate v. The unit of all rates is per cell per day.

Generating function approach

Denote by Fly(X; t) the probability generating function for the lineages at time t initiated by a single ij-type cell, Where X 2 [x00, x01, x10, x11] T denotes the vector of dummy variables with elements xij representing each ij-type of cells. The backward equations for this branching process are (see 81 Text for how to derive them)

Probability of acquired resistance

Similar methods have been used in the literature [104, 105]. For example, the probability density of the number of resistant cells, pmn(t), in both compartments at time t, as shown in Fig. 3, is given by Using the Cauchy’s integral formula, we obtain Applying the trapezoid rule to approximate the double integral, we arrive at:

To reduce the aliasing effect arising from spectral methods, we set N = 1000(>> 20) for the results presented in Fig. 3. Conditional time to resistance We use the conditional probability, p5(t), of having no resistant cells in both compartments to determine the average time to resistance starting With sensitive cells in either compartment. Specifically, starting with a single sensitive cell in compartment 1, the conditional probability Then average time to resistance Tr (relapse time) can be calculated as follows,

Full model

Drugs are distributed over M spatial compartments according to given levels of spatial concentration heterogeneity. In particular, we consider two different schemes of migration: local migration versus global migration. Local migration means that compartments are situated on a “ring” where a cancer cell can only migrate to the two nearest neighbor compartments with equal probability v/ 2. In contrast, global migration means compartments are fully connected where a cancer cell is allowed to migrate from one compartment to any other one with equal probability v/ (M —1). This extended model allows us to study how metastatic cancer cells acquire increasing levels of drug resistance as they migrate along a spatial gradient of drug concentration and thrive in areas of excessively high drug concentrations, as shown in a recent in vitro experiment [54].

Supporting Information

Panels (A) - (D) plot the joint probability density distribution of the numbers of resistant cells in both compartments at time points 1‘: 10, 50, 100, 200, respectively, starting with a single sensitive cell placed in the sanctuary compartment. Parameters are the same as in Fig. 3.

Relapse occurs sooner with bigger tumor size at the start of therapy. Relapse is destined to happen in the presence of drug sanctuary and cell motility. The situation is even worse for really big tumor sizes at the start of therapy; relapse can happen within weeks.

(A) shows the escape probability as a function of the difference in drug concentrations between the two metastatic compartments, AD. (B) shows how the average time to resistance changes with increasing AD. The vertical line marks the critical value of AD above which two point mutations are required to confer sufficient levels of resistance to increasingly high concentrations in compartment 1 while compartment 0 is the sanctuary containing lower level of drugs. Parameters: p = 3.5, v = 10—4, and other parameters are the same as in Fig. 4.

Shown are the escape probabilities and the average conditional time to resistance for a single sensitive cell initially placed in each compartment, with increasing migration rate, v. The spatial heterogeneity in drug concentrations is realized by using a rescaled Normal distribution with a peak in the central compartment, and the level of heterogeneity is denoted by the standard deviation, 0, of the concentration distribution over compartments. Parameters: M = 20, n = 5,1(350 = 100, m = 2, p = 1.1, D = 50, s = 0.01, 190 = 0.2, do = 0.1, [,1 =10‘4,a= 24.7. (TIF)

Shown are the escape probabilities and the average conditional time to resistance for a single sensitive cell initially placed in each compartment, with increasing spatial heterogeneity, a. The spatial heterogeneity in drug concentrations is realized by using a rescaled Normal distribution with a peak in the central compartment, and the level of heterogeneity is denoted by the standard deviation, 0, of the concentration distribution over compartments. Parameters: M = 20, n = 5, ICSO = 100, m = 2, p = 1.1, D = 50, s = 0.01, 190 = 0.2, do = 0.1, [,4 =10‘4,v = 0.01. 81 Text. Supplementary Information for “Spatial heterogeneity in drug concentrations can facilitate the emergence of resistance to cancer therapy”. (PDF)

Acknowledgments

We thank Jacob G. Scott and Philip Gerlee for careful reading of and helpful comments on an earlier version of this manuscript.

Author Contributions

Performed the experiments: FF MAN SB. Analyzed the data: FF MAN SB. Contributed reagents/materials/analysis tools: FF MAN SB. Wrote the paper: FF MAN SB.