Spatial Heterogeneity in Drug Concentrations Can Facilitate the Emergence of Resistance to Cancer Therapy

Acquired resistance is one of the major barriers to successful cancer therapy. 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, otherwise 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.


Full model
For completeness, here we present in details the model description and the analytical methods. In this paper, we use a simple spatial model to study the evolution of resistance to cancer therapy, particularly in the presence of sanctuary sites with low or no drug concentration. Our approach is based on branching process theory. The present study focuses only on (single-) drug resistance arising during treatment, rather than pre-existing resistance prior to treatment. When cancer patients with metastatic diseases are administered potent drugs, it is often the case that the drug is not uniformly distributed throughout the body compartments of the cancer patient due to various factors such as drug penetration problems. This leads to non-homogeneous drug concentrations across different body compartments, which may harbor small metastatic populations. These compartments with low drug concentrations, which are insufficient to fully inhibit the replication of cancer cells, become the sanctuaries from which drug resistant mutants can emerge. Through migration these resistant mutants can subsequently populate compartments with high drug concentration, where the sensitive type cannot survive. By such intuitive reasoning, we hypothesize that the spatial heterogeneity in drug concentrations can speed up the evolution of acquired resistance to cancer therapy. In the rest of this supplementary information, we shall derive mathematical results to investigate this hypothesis. Our theoretical results shed new light onto the emergence of drug resistance and treatment failure of cancer patients with metastatic diseases.
Let us consider a spatial compartment model of metastatic diseases. There are M compartments in total and the drug concentration in each compartment i is denoted by D i . We restrict our analysis on spatially heterogeneous drug concentrations, but without time-changing concentration fluctuations. Metastatic cancer cells have already acquired motility and are assumed to migrate between compartments with rate v. The effect of the drug can either be cytostatic inhibiting cell growth and division or be cytotoxic killing cells directly, or both. Without loss of generality, we assume here that the drug inhibits cell growth. Reproduction is subject to mutation. Upon division, one of the two daughter cells mutates with rate u. By acquiring step-wise point mutations, the mutant becomes more and more adapted to drug environments and can eventually survive in high-concentration compartments. Denote by i the genotype of a cancer cell if it has acquired i point mutations. The replication rate of a cancer cell, b i j , is determined by its genotype i and spatial location j as follows, (1) Here we use a Hill function for the drug response curve. β j is the division rate of a wild type cell in compartment j in the absence of drug, s is the cost of mutation in the absence of drug, IC 50 is the drug concentration that is needed to inhibit cell growth by one half of its original rate, ρ is the fold increase in IC 50 per mutation, and m determines the steepness of the Hill function. We note that drug mutations are deleterious in the absence of drugs. Thus the replication rate of a mutant with genotype i and in compartment j is (1 − is)β j in the absence of drugs, whereas its IC 50 is increased to ρ i IC 50 in the presence of drugs. Throughout this supplementary information, unless stated otherwise, the death rate of a cancer cell with genotype i within compartment j, d i j , is all the same across compartment and equal to that of the wild type in the absence of drugs, The migration pathway is specified by the underlying connections between compartments. For simplicity, 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). We restrict our subsequent analysis to these two foregoing migration patterns, though it is straightforward to extend it to more complicated migration patterns in our mathematical model.
In what follows, we shall use a continuous-time multi-type branching process to describe resistance evolution.

Analysis
As for multi-type branching process, we first need to construct the appropriate type space using the combination of genotype g i ∈ {0, 1, · · · , n − 1} and compartment location l j ∈ {0, 1, · · · , M − 1}. Thus a cell's type can be denoted by two-bit strings g i l j ∈ {0, 1, · · · , n − 1} × {0, 1, · · · , M − 1} and thus there are n × M different types. For example, a cell with "11"-type means that this cell has accumulated 1 point mutation and is located in compartment 1. Therefore, a cell's type, notated by its genotype and spatial location, can change as results of genetic mutation or migration.
Next let us introduce the probability generating functions. Denote by F i j (X; t) the probability generating function for the lineages at time t initiated by a single i j-type cell, where X = [x 00 , . . . , x n−1,M−1 ] denotes the vector of dummy variables with elements x i j representing each i j-type of cells. During an infinitesimal time interval ∆t, with probability b i j (1 − u)∆t an i j-type cell gives birth to two identical daughter cells, with probability b i j u∆t it gives rise to an identical offspring and a mutated offspring, with probability v∆t it migrates to another compartment, with probability d i j ∆t it dies, and with probability 1 − (b i j + d i j + v)∆t nothing changes. We can use the probability generating function approach to account for what can happen to an i j-type cell during an infinitesimal time interval ∆t, Using the multiplicative property of branching process (the independency of each individual lineage), the generating function F i j (X; t + ∆t) satisfies the recursive backward equation, where F is the vector of the probability generating functions with elements F i j (X; t). We obtain the following backward equation for global migration (0 ≤ i < n − 1, 0 ≤ j ≤ M − 1).
Since we do not consider more than i = n − 1 point mutations, for i = n − 1 we have The initial condition for these ODE's is F i j (X; 0) = x i j for i = 0, . . . , n − 1 and j = 0, . . . , M − 1. For any given profile of the initial numbers of cells, N = {N i j }, the corresponding probability generating function F N is a power compound of F i j Analogously with Eq. (4), we can derive the backward equation for local migration, In the above equation, the subscripts [( j − 1 + M)%M] and [( j + 1)%M] account for the two neighbor compartments the cell may migrate to. In the following unless stated otherwise, we have d i j = d 0 for all types of cells, while b i j depends on the drug concentration in compartment j as in Eq. (1). The fate of a i j-type cell can either go extinct in a relatively short time or successfully establish surviving lineages with expanding populations. The extinction probability, which is the opposite of escape probability, can be obtained by solving the fixed points of the ODE's given above. We can also calculate the probability that there is no evolved fully resistant strain in any compartment by time t, conditional on non-extinction.
For arbitrary number of genotypes and arbitrary number of compartments, we cannot find closed form probability generating functions for these nonlinear ODE's. To circumvent this difficulty, we rely on numerical ODE solvers to obtain the quantities of interest. For instance, if we solve the above nonlinear ODE's with the initial condition x i j = 0, we can obtain the extinction probability F i j (0; t) for each i j-type cell at time t. For t → ∞, F i j (0; t) converges to the steady equilibrium value (i.e., the smallest of all admissible fixed points). Using the initial condition X noR = [x i j ] where x i j = 1 for i < n − 1 and x i j = 0 for i = n − 1, the conditional probability P noR (t) that there is no evolved fully resistant strain in any compartment by time t, starting with an i j-type cell, can be calculated as, In general, for any given profile of the initial numbers of cells, N = {N i j }, the extinction probability P E (t) at time t is given by The conditional probability P noR (t) is given by

Simple model
For the simple model with 2 genotypes and 2 compartments, we are able to derive closed-form approximations for special cases. The backward equations are as follows, These are coupled nonlinear Riccati equations, for which we do not have closed-form explicit solutions. For simplicity, let us further assume dichotomous drug environments: compartment 0 is drug free, while compartment 1 is distributed with a drug that inhibits cell growth. We have d i j = d 0 .
In the drug-free compartment 0, both sensitive and resistant cells have supercritical replication potential in drug free compartment, but resistant cells have slightly lower replication rates than sensitive type due to the cost of the resistant mutation. We have b 00 > d 00 , b 10 > d 10 , and b 10 < b 00 . In the drug-containing compartment 1, sensitive cells have a subcritical replication potential while resistant cells still have a supercritical replication potential. We have b 01 < d 01 , b 11 > d 11 , and b 11 > b 01 .
For such a fitness landscape, the drug-free compartment provides far more favorable condition for breeding resistance than the drug-containing compartment. Therefore, the likely origin for treatment failure, these resistant cells that populate the drug-present compartment, should be attributed to the "mutation-migration" pathway than the "migration-mutation" pathway. In other words, it is less likely that resistance evolves in situ in the drug-containing compartment than that resistance mutants are originated from the drug-free compartment and then migrate to the drug-containing compartment.

Drug-environment-dependent escape
To verify these arguments above, let us first calculate and compare the likelihood of obtaining an escape mutation within one separate compartment (without any migration at all). For simplicity, the compartment is only seeded with a single sensitive cell. Denote by b 0 and b 1 the division rate of sensitive cells and resistant ones in that compartment. The death rate of the two types of cells is the same, d 0 = d 1 . We have b 0 > b 1 > d 0 for compartments with low concentrations, whereas b 0 < d 0 < b 1 for compartments with high concentrations. Let f 0 (x, y, t) denote the probability generating function for the lineages starting with a single sensitive cell, and let f 1 (x, y, t) denote the probability generating function for the lineages starting with a single resistant cell. We have the following differential equations: The initial condition is f 0 (x, y, 0) = x and f 1 (x, y, 0) = y.
For sufficiently small mutation rates u 1, we can derive a closed form approximation to the equations above. Let r i = b i − d i denote the net growth rate of type-i cells, i = 0, 1. We have r 0 > r 1 > 0. Neglecting the loss of sensitive cells due to rare mutations, we have This ODE above is explicitly soluble, and we obtain: Denote by X 0 (t) the number of sensitive cells at time t. Previous results show that X 0 (t)e −r 0 t is a non-negative martingale and X 0 (t)e −r 0 t → W 0 for t → ∞. The Laplace transform of W 0 is given by For t → ∞, using e −θe −r 0 t − 1 ≈ −θe −r 0 t , we obtain Accordingly, we obtain the limiting distribution of W 0 for t → ∞ as follows.
where δ 0 is a point mass at zero (corresponding to extinction), and Exponential( r 0 b 0 ) is an exponential distribution with the exponent r 0 b 0 . It is easy to see that, conditional on non-extinction, X 0 (t)e −r 0 t → V 0 = Exponential( r 0 b 0 ). Thus the Laplace transform of V 0 is The rate at which resistant mutants arise is ub 0 X 0 (t). Hence the conditional probability that the arrival time of the first mutant τ 1 > t is The median arrival time We can see that t1 2 is monotonically decreasing with increasing b 0 (u 1). That is, the less harsh the compartment is (b 0 ↑), the sooner resistant mutants arise (t1 2 ↓). In this sense, compartments with lower drug concentrations are more likely the breeding ground for resistance, since the replication of sensitive cells is less affected than in compartments with higher drug concentrations.
For small mutation rate u 1, the escape probability P 1 that a single sensitive cell establishes lineages without going extinction is approximately In this unfavorable case, sensitive cells do not have a chance to survive but resistant cells do. The approximation method used above does not work here, since the branching process of sensitive cells is subcritical. To establish surviving lineages, sensitive cells must evolve resistant mutation right before quickly dying out. We cannot find a closed form for the arrival time of the first mutant for this case, but it is easy to see that if the mutant does arise (though with a lower chance), it takes much longer time due to the smaller birth rate than in benign environments. The escape probability of a single sensitive cell is approximately which is much smaller than Eq.(21). Moreover, we should note that the escape probability P 1 is of the order √ u if sensitive cells are almost critical |r 0 | → 0:

Competing pathways to the emergence of resistance
Returning to the simple model introduced above, we are now primarily concerned with which pathway leads faster to resistant cells that populate the drug-present compartment. To this end, we compare two distinct pathways. The migration-mutation pathway: sensitive cells first migrate from the drug-free compartment to the drug-present compartment, and then acquire resistance in situ. The mutation-migration pathway: sensitive cells acquire resistant mutation in the drug-free compartment, followed by migration to the drug-present compartment. As shown above, we cannot obtain closed-from solutions for this four-type branching process. To simplify our mathematical derivations, let us consider unidirectional emigration flow only from the drug-free compartment to the drug-present compartment. To compare which pathway is faster, let us further artificially separate these two pathways, making only one pathway at work at a time, instead of considering them simultaneously. In this way, each pathway can seen as a three-type branching process: For the migration-mutation pathway, µ 0 = v and µ 1 = ub 01 ; for the mutation-migration pathway, µ 0 = ub 00 and µ 1 = v. The difference of the fitness landscape between these two pathways lies in the intermediate type 1, namely, the division rate b 01 of sensitive cells in the drug-present compartment versus the division rate b 10 of resistant cells in the drug-free compartment. It is therefore useful to first study a general three-type branching process and then obtain the results for each pathway by substituting the specific fitness landscape for each pathway.
For rare mutations (u 1) and low migration rates (v 1), we are able to derive explicit closed-form approximations for the conditional probability of no type 2 cells by time t. Let b i and d i be the birth and death rate of type i cells, i = 0, 1, 2. Let r i = b i − d i be the net growth rate of type i cells. We have r 0 > 0, r 2 > 0, and r 1 < r 0 . Conditional on non-extinction, the number of type 0 cells at time t is given by Because type 1 is less fit than type 0 (r 1 < r 0 ), we can approximate the number of type 1 cells at time t as where x(t) ∼ y(t) means x(t)/y(t) → 1 for large t.
The conditional probability that there are no surviving type 2 cells by time t is where the integral θ 2 (t) = µ 0 µ 1 r 2 r 0 −r 1 t 0 e r 0 s b 2 −d 2 e −r 2 (t−s) ds. Substituting the fitness landscape of each pathway, we find a simple condition for the mutation-migration pathway faster than the migration-mutation pathway in leading to surviving resistant cells in the drug-containing compartment if the follow inequality holds: .
Concerning cytostatic drugs that inhibits cell division, we have d 00 = d 10 = d 01 = d 11 = d 0 , b 00 = b 0 , b 10 = (1 − s)b 0 and b 01 = (1 − δ)b 0 . s denotes the fitness cost of resistant mutation in the absence of drugs, and δ denotes the fitness cost of sensitive cells in the presence of drugs.
Then the above condition simplifies to For cells in compartments with potent drugs and with sufficiently high drug concentrations, we have δ > s. More generally, the emergence of resistance is more contingent on the mutationmigration pathway than the migration-mutation pathway, if the difference of drug concentration between two neighboring compartments leads to the landscape satisfying the inequality (27).

Numerical results
In what follows, we report the numerical results on the full model. To avoid boundary problems, we consider generic, symmetric spatial distributions of drug concentrations over compartments, according to a rescaled Normal distribution such that the degree of heterogeneity can be tuned by the variance of the distribution while keeping the total sum of concentrations constant. We compare the extinction probabilities and the time to resistance for sensitive cells initially located in each compartment, with varying migration rate and the spatial heterogeneity in drug concentrations.
The results are consistent with those reported in the main text (Figs. S5 and S6). Sanctuary sites are compartments with low drug concentration, which are responsible for producing the resistance that populates compartments with high drug concentration. Selection for resistance is more contingent on the mutation-migration pathway than the migration-mutation pathway, particularly in the presence of sanctuary sites. Moreover, only for small migration rates below a certain threshold does the spatial heterogeneity in drug concentrations help speed up the resistance evolution. However, excessively high migration rates actually slow down resistance emergence, as the role of compartment structure is diminished by frequent migration.
Furthermore, we observe interesting results when comparing local migration to global migration. Both the extinction probability and the time to resistance depend more sensitively on migration rate for global migration than for local migration (Fig. S5). This is because global migration favors the dissemination of resistance all over the place, whereas local migration constrains the migration within neighboring compartments with reduced drug difference such that accruing one point mutation is sufficient to offset the increasing drug concentration. As a consequence, resistance evolves much faster for local migration than global migration, particularly for sensitive cells located at the sanctuary sites with low drug concentrations. On the contrary, the escape probability of sensitive cells located in the compartments with high concentrations is hampered more by local migration than by global migration. In this case, local migration limits its dispersal range, while global migration increases the chance of moving to sanctuary sites thus leading to greater escape probabilities. More specifically, compared to global migration, local migration speeds up the evolution of resistance for sensitive cells initially located at sanctuary sites, while delaying the evolution of resistance for sensitive cells initially located at harsh compartments with high levels of drugs. Regarding the latter case, moreover, there exists an optimal intermediate global migration rate with which the time to resistance is minimal.
Sanctuary sites for sensitive cells result when varying from homogeneous drug distributions to increasingly heterogeneous drug concentrations across spatial compartments. It is obvious that compartments with low concentrations become the shelter for sensitive cells and thus provide benign condition for breeding resistance that can be selected for in compartments with excessively high concentrations (Fig. S6). Local migration takes less time to evolve resistance than global migration, if these cells are initially located at the sanctuaries. Opposite results are obtained if these cells are initially located at harsh compartments. Under local migration, there exists an optimal level of heterogeneity in drug concentrations that most delays the emergence of resistance. In contrast, increasing the spatial heterogeneity in drug concentrations helps prolong the time to resistance under global migration.
Altogether, these results confirm that our conclusions derived from the simple model in the main text are robust with respect to model extensions and varying model parameters. In addition, the scheme of migration plays an important role in the evolution of acquired resistance. However, future studies are needed in order to characterize, and eventually manipulate, the migration routes for metastatic cells. We think that constraining escape route of metastatic cells to sanctuary sites and limiting dissemination of evolved metastatic cells from sanctuary sites are vital to eliminate disseminated cancer.