{'fraction_non_alphanumeric': 0.03789102385593614, 'fraction_numerical': 0.009059686252668708, 'mean_word_length': 4.662952060555089, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Convection shapes the trade-off between antibiotic efficacy and the selection for resistance in spatial gradients

Since penicillin was discovered about 90 years ago, we have become used to using drugs to eradicate unwanted pathogenic cells. However, using drugs to kill bacteria, viruses or cancer cells has the serious side effect of selecting for mutant types that survive the drug attack. A key question therefore is how one could eradicate as many cells as possible for a given acceptable risk of drug resistance evolution. We address this general question in a model of drug resistance evolution in spatial drug gradients, which recent experiments and theories have suggested as key drivers of drug resistance. Importantly, our model takes into account the influence of convection, resulting for instance from blood flow. Using stochastic simulations, we study the fates of individual resistance mutations and quantify the trade-off between the killing of wild-type cells and the rise of resistance mutations: shallow gradients and convection into the antibiotic region promote wild-type death, at the cost of increasing the establishment probability of resistance mutations. We can explain these observed trends by modeling the adaptation process as a branching random walk. Our analysis reveals that the trade-off between death and adaptation depends on the relative length scales of the spatial drug gradient and random dispersal, and the strength of convection. Our results show that convection can have a momentous effect on the rate of establishment of new mutations, and may heavily impact the efficiency of antibiotic treatment.


Introduction
The emergence of drug resistance represents one of the major clinical challenges of the current century [1][2][3]. Microbial pathogens quickly acquire resistance to new antibiotics [4], while solid tumors often regrow after treatment because of resistance mutations that arise during tumor growth [5]. In addition to genomic studies examining the molecular causes of resistance [6,7], the dynamics of drug resistance evolution has recently attracted wide interest [8,9], with the dual goal of understanding the emergence of resistance and developing novel strategies to prevent or control its spread [10,11]. Next-generation sequencing and high-throughput experimental techniques enable the quantitative study of resistance evolution but require the development of new theories to appropriately interpret experimental results [12].
In many realistic systems, an evolving population interacts with its surroundings and exhibits a well-defined spatial structure (for instance, in tumors and biofilms [5]). It has recently been shown that this spatial structure can strongly influence the subclonal structure and the adaptation of spatially expanding populations, both from de novo and pre-existing mutations [13][14][15]. In addition, the presence of spatial drug gradients is well documented both in the outside environment [16,17] as well as within biofilms [18] and the human body [19][20][21][22]. The presence of spatial [12,[23][24][25][26][27] and temporal [10,[28][29][30] heterogeneities has been shown to facilitate the emergence of drug-resistant phenotypes and enable populations to reach a higher degree of resistance than in homogeneous drug concentrations. For instance, in a microfluidic experiment, a spatial gradient gave rise to a higher rate of adaptation of bacterial populations [26]. Similarly, microbes growing on soft agar plates with gradually increasing antibiotic concentrations were able to rapidly evolve resistance to high levels of antibiotics, while sudden jumps to unsustainably high concentrations dramatically slowed down adaptation [12]. Moreover, many realistic growth scenarios of bacterial populations may be subject to convection driving them up or down the gradient. Examples include the gut, arteries, and urethra in the human body [31][32][33][34], but also flows in aquatic environments, like ocean and river currents [35], or flow in pipes and catheters [36].
A number of recent theoretical studies have investigated how gradients speed up the evolution of drug resistance [37][38][39]. Greulich et al. [38] considered a population adapting to a smooth gradient, which gradually lowers the growth rate of susceptible individuals. Hermsen et al. [37] studied resistance evolution in a series of sharp step-like increases in concentration, where a novel resistance mutations was necessary for survival in the next step (the "staircase" model); Hermsen [39] later proposed a generalization of the staircase model to continuous gradients. These previous studies focused on the speed of adaptation, i.e., how quickly the population evolves to tolerate high concentrations of antibiotics. In the context of the emergence of drug resistance, however, this observable alone ignores a crucial reality of antibiotic treatment: efficient drug treatment first and foremost aims to kill as many bacteria as possible, while limiting the rise of resistance mutation [40,41]. How this apparent trade-off can be optimized for populations in  Figure 1. Sketch of our modeling setup. We assume that, initially, a purely wildtype population has reached a steady-state density profile (blue) in the presence of a steady-state antibiotic concentration gradient (orange). Resistance mutations occur spontaneously in randomly drawn individuals and disperse, proliferate and die until extinction or ultimate fixation. Convection can either drive the population towards the antibiotic (co-flow, v > 0) or away from it (counter-flow, v < 0). Our goal is to analyze the establishment probability of resistance mutation for a given rate of wild-type killing.
spatial gradients to prevent the evolution of drug resistance has so far been unexplored.
Here, we present simulations, rationalized by a comprehensive analytical framework, of populations evolving resistance in a variety of spatial antibiotic concentration gradients and under the influence of convection. We measure the establishment probability of individual resistant mutants arising in a region occupied by susceptible wild type and find that successful, "surfing", mutations arise in a localized population patch close to the population front, the size of which depends on the relative strength of bacterial diffusion, antibiotic gradient steepness, and convection. We find that shallow gradients and flow towards higher antibiotic concentrations promote wild-type death, at the cost of increasing the establishment probability of resistance mutations. Conversely, populations in steep gradients and subject to flow towards lower antibiotic concentrations give rise to fewer drug-induced wild-type deaths but also produce fewer resistance mutants. We introduce the notion of a treatment efficiency, which quantifies this inherent trade-off between adaptation and death, and find it to be strongly modulated by gradient steepness and convection.

Definition of the model
Consider a population in a fixed antibiotic concentration gradient. Where the antibiotic concentration is too high, the population cannot survive unless it evolves resistance to the antibiotic, and thus at steady-state, a population density profile c(x) will develop. For an individual mutation conferring antibiotic resistance that occurs at position x in the population we define the probability u(x) that this mutation will be successful, i.e., establish first locally and eventually colonize areas where the antibiotic concentration is too high to allow for growth of the wild type. In doing so, we implicitly assume a very low mutation rate and neglect clonal interference to focus on the fate of individual mutations (SI section 6). We come back to this question in the Discussion.
A drug gradient arises upon administering an antibiotic that is introduced at one location in the system (the source) and flows out of the system at another (the sink). The antibiotic may be influenced by convection, e.g., through blood flow or peristaltics in the gut, and be subject to degradation. At steady-state, a concentration gradient between source and sink is established, which may take a range of shapes, from a sharp, steplike gradient for strong convection, to a shallow gradient over the whole system size. Using a simple reaction-diffusion model, we show in the Supplementary Information (SI section 1) how the gradient depends on the distance between source and sink, the diffusivity of the antibiotic, and the speed of convection (see SI). To keep the discussion general, here we approximate the features of a typical antibiotic gradient by modeling it as a sigmoidal function that changes over a characteristic length scale λ, which sets the gradient steepness, without trying to associate a particular value of λ with a specific combination of real-life parameters.
The wild-type population growing in this gradient will in general also be subject to convection with flow speed v. However, it is important to note this flow speed may be very different from the convection that affects the antibiotic gradient. As an example, the flow speed inside a blood vessel depends on the distance from the wall, and thus, surface-bound populations such as biofilms may experience much lower flow speeds than antibiotic in the bulk fluid. More generally, it is to be expected that fluid flow differentially affects bacterial cells and antibiotic molecules due to their different sizes, potential porosity of the surroundings, adhesion effects, etc, such that convective flows are arguably the rule rather than the exception. Therefore, without loss of generality, we may ignore the details of how a particular steady-state antibiotic gradient is generated, and focus on how the combination of antibiotic gradient and convection shapes the emergence of resistance.
We model the effect of the antibiotic by a drug-induced death rate b(x) of susceptible wild-type individuals giving rise to a net growth rate s(x) of the wild type that ranges from the maximal growth rate a 0 (in antibiotic-free regions) to some negative net growth rate (where the antibiotic concentration is high). For simplicity, we model the death rate such that at the highest concentration wild-type individuals typically die within one generation. Given a net growth rate profile, we can compute the steady-state wild-type population density c(x), whence we obtain the number B of drug-induced wild-type deaths per generation, The total death rate B quantifies the efficacy of the antibiotic; it corresponds to the rate at which the population is hindered in its growth by the antibiotic. Note that since we only consider steady-state population densities and thus do not explicitly simulation the eradication of the wild-type population, our analysis relies on the immune system or slowly progressing antibiotic concentration gradients to completely eradicate the wild-type population [41].
To quantify the emergence of resistance, we measure the local mutant establishment probability u(x). Since the probability that a mutation occurs in the first place is proportional to the wild-type population density c(x), it follows that successful mutants can only arise where both the wild-type population density and the establishment probability are high (see Fig. 2). A measure for how readily new resistant mutants establish is thus given by the product of wild-type population density and the establishment probability [42], summed over the whole population [43]. We call this measure the rate of adaptation, R, The rate of adaptation, R, is proportional to the rate at which new resistance mutations arise (at a low mutation rate µ, see SI section 6) and establish in the population. Alternatively, the rate of adaptation, R, can be understood as a measure proportional to the mean establishment probability c(x)u(x)/ c(x), i.e., the probability that a mutation arising anywhere in the population establishes.
Finally, we define the treatment efficiency Q as the ratio of the effective reduction in growth and the rate of adaptation to the antibiotic: a treatment is deemed particularly efficient if it can reduce the growth of the wild-type population while at the same time hindering the emergence of resistant phenotypes as much as possible. In our model, this corresponds to defining treatment efficiency Q = B R . ( As we shall see below, the rate of adaptation, R, and the total death rate B typically follow the same trends, e.g., they are both larger in shallow gradients than in steep ones, as may perhaps be suspected intuitively. By contrast, the treatment efficiency Q will turn out to be less accessible to intuition and require a detailed understanding of the population density and establishment probability profiles. In the following, we first present simulation results and then turn to analytical theory to rationalize our findings.

Simulation results
We simulate a population of wild-type individuals on a lattice of L demes, where each individual can migrate into a neighboring deme, replicate, and die, using a Gillespie algorithm [44,45] (see Methods). Wild-type population growth is limited to a carrying capacity K by a logistic death term [46]. The population is subject to convective flow with speed v and set in a one-dimensional antibiotic gradient (see the sketch in Fig. 1), which interpolates between maximal and zero antibiotic concentration over a characteristic length scale λ. The antibiotic gradient induces a wild-type death rate b(x), giving rise to an effective growth rate s(x) for the wild-type. For simplicity, we choose the functional form such that wild-type growth is strongly inhibited for x 0 and proceeds unhindered with growth rate a 0 for x 0. Following the equilibration to the steady-state profile c(x) of the wild type, a single resistant mutant is inserted into the population at position x. The resistant mutant has the same birth rate as the wild type and is subject to the same carrying capacity as the wild-type, but it does not suffer from an increased death rate due to the antibiotics. We follow the mutant clone until it either goes extinct or reaches the far end of the simulation box, in which case we consider the mutant established. The establishment probability u(x) is then equal to the fraction of simulations in which a mutant introduced at x managed to establish. Fig. 2 shows the resulting population density c(x) (top row) and establishment probability u(x) profiles, for three different gradients (columns). Within each panel, different colors represent the profiles resulting under different flow speeds. Here, we adopt the usual convention of positive flow speed pointing to the right, which in our case corresponds to flow towards high antibiotic concentrations. We will refer to such flow as co-flow, and represent it graphically in hues of cyan. Conversely, negative flow speeds point to the right, towards lower antibiotic concentration; we term this counterflow and use red tones.  Figure 3.
Quantifying the trade-off between killing of wild-type cells and establishment of resistant mutants. (a) the total death rate B and the rate of adaptation, R, for different antibiotic gradients (shallow gradient, λ = 100, green; intermediate gradient, λ = 10, blue; step, purple). Counter-flow (v < 0) leads to drastic decrease of both B and R, while co-flow has a milder effect, in accordance with the individual population and establishment probability profiles in Fig. 2

. (b)
The treatment efficiency Q, defined as the number B of drug-induced deaths per generation divided by the rate of adaptation, R, (see eq. 3). Counter-flow can increase the treatment efficiency by an order of magnitude in shallow gradients because R is reduced by the effect of convection on both u(x) and c(x), while B only captures changes in c(x).
In the absence of convection, shown in black, the wild-type population density c(x) roughly follows the net growth rate s W T (x) (gray dotted line), and the establishment probability u(x) is generally high where c(x) is low, and vice-versa, since the mutants compete for resources with the wild-type.
Convection affects the wild-type population density profile as intuitively expected, by stretching or compressing the population spatial range for co-flow and counter-flow, respectively. Interestingly, while counter-flow significantly alters the profile, the effects of co-flow are hindered by the drug profile, which prevents the wild-type from growing too far into the antibiotic region.
The effect of convection on the establishment probability is more complex, due to the generation of two competing processes that either help or hinder the mutant success. On the one hand, the changes in the wild-type profile described above alter the competition with the mutants, reducing the establishment probability in co-flow conditions, and increasing it in counter-flow conditions. On the other hand, mutants are also transported towards the direction of flow: their success rate is thus increased by co-flow and reduced by counter-flow. The presence of these two opposing forces generates interesting cut-offs along the profile, which we characterize in the Theory section.
The relative strength of the two competing effects can be quantified by the product u(x)c(x), which represents the probability density of successful mutants (bottom row in Fig. 2). Here we see that co-flow widens the probability peak, while counter-flow shrinks it and moves it away from the antibiotic region. This suggests that transportation of the mutants in the direction of flow rather than competition with the wild-type is the dominant effect in the mutant success rate. Since the area under these curves is the rate of adaptation, R, defined in eq. 2, we find that both the rate of adaptationand the total death rate B are higher in co-flow conditions and shallow gradients and lower in steep gradient with counter-flow, as shown in Fig. 3a. Although R and B follow the same rough trends with flow speed v, the rate of adaptation, R, is generally more strongly affected by flow than the total death rate B. This is because flow alters both c(x) and u(x), but leaves the local death rates b(x) unchanged.
This has dramatic effects on the treatment efficiency Q = B/R, shown in Fig. 3b. Shallow gradients are generally conducive to a larger treatment efficiency Q than steep gradients. In particular, strong counter-flow conditions (v < 0) in shallow gradients give rise to very small R compared to B, such that the treatment efficiency can be more than an order of magnitude larger than in the no-flow case. In steep gradients, although both R and B change over two orders of magnitude for different flow speeds, the treatment efficiency Q remains roughly constant because both R and B are equally affected by flow (Fig. 3a, purple lines).
In summary, our simulations show that the population density profile c(x) as well as the local establishment probability u(x) of resistance mutants are both strongly influenced by environmental parameters, in particular, by the steepness of the antibiotic gradient and the strength of convection. In the following, we lay out a mathematical model based on the theory of branching random walks that reproduces our key findings and elucidates the relative roles of gradient steepness and flow speed; detailed calculation are mostly relegated to the Supplementary Information.

General framework
We assume that the wild-type population density c(x, t) (rescaled by the carrying capacity K) is described by the reaction-diffusion equation [47] where a W T (x) is the local wild type birth rate, b W t is the local antibiotic-induced death (b) Co-flow (ν > 0, cyan) gives rise to broader algebraic establishment probability profiles, while counter-flow (ν < 0) gives rise to an exponential decrease. Both cases are asymptotically captured by the analytical approximations (solid lines, eqs. S9, S15 and S16. is the local net growth rate of the wild type, and v is the external flow speed, representing, e.g., blood flow. This model is a straight-forward generalization for the standard Fisher model [48] to spatially inhomogeneous growth rates. Our model ensures that the steady-state local population density c SS depends explicitly on the local death rate b W T (x) when the death and birth rate profiles change sufficiently slowly in space, Our model, like the original Fisher equation, ignores the discrete nature of individuals, which has been shown to significantly alter the tip of the population front [49,50], especially when the carrying capacity is small [51] (see SI section 5 for a detailed discussion). Nevertheless, we expect good agreement between our model and simulations in terms of the total death rate B and the rate of adaptation, R, , whose values do not significantly depend on the population profile at the tip of the front. Given a single resistance mutation arises in the population at position x, its probability u(x, t) to survive for a time t can be derived by modeling the mutant lineage as a branching random walk; it obeys a nonlinear reaction-diffusion equation (see Ref. [43] and SI section 2), where a M T (x) is the local birth rate of the mutants and is their net growth rate. The establishment probability u(x), i.e., the ultimate survival probability for a mutation born at position x, as measured in simulations, is given by the steady-state solution of eq. 7.
To mimic the situation in our simulations, we assume that the birth rate of wildtype and mutant is identical and constant, a W T (x) = a M T (x) = a 0 , while the wild-type drug-induced death rate b W T (x) ranges from −a 0 to a 0 . This implies that effect of the antibiotic is to increase the death rate of the wild type, while the drug-induced death rate of the resistant mutant is zero. The effective growth rate of the mutants is thus determined purely through competition with the wild type, i.e., s MT To model random dispersal and external flow, we have included diffusion and convection terms in eqs. 5 and 7 (note the difference in sign between the convection terms). To get a feel for solution to the set of equations 5 and 7, we first study the case v = 0 before turning on convection (v = 0).

No flow
We begin by considering the simplest functional form for the antibiotic gradient -a step-like increase in concentration at x = 0 that gives rise a net growth rate of a 0 for x < 0 and −a 0 for x > 0, i.e., Such a sharp gradient could emerge, for instance, at the boundary of different tissues or organs with different affinities to store antibiotics [19,25]. Upon rescaling the spatial coordinate by the characteristic length scale = D/a 0 , which can be intuitively understood as the typical distance that a mutant individual travels through random dispersal before replicating, eq. 5 for the wild-type population density in this case becomes where ξ = x/ . Its solution, given analytically in SI section 3, transitions quickly from 1 to 0 over a distance ∼ 2 (Fig. 4a).
To find the establishment probability u(ξ) far from the transition, we solve which, given the exact solution c(ξ), can simply be integrated numerically (see Methods).
To make analytical progress, we approximate the wild-type population density with a step, i.e., c(ξ) ≈ Θ(−ξ) and ask for solutions to eq. 10 far from the transition region. For ξ > 0, u(ξ) also approaches 1 exponentially quickly over a distance , since the lack of competition with the wild-type in this region facilitates the establishment of resistance mutations. Instead, far inside the bulk of the population (ξ < 0), we find that the establishment probability exhibits a long tail, u(ξ) ∼ ξ −2 (Fig. 4a). The long tail implies that mutations arising deep inside the bulk of the population still have a relatively high chance of establishing, which explain the asymmetry of the simulation curves in Fig. 2c.
For gradients varying over length scales longer than , we model gradients decaying over a characteristic length scale λ with a net growth rate s(ξ) of the form eq. 4. Introducing the steepness parameter m = /λ, which quantifies the relative length scales associated with diffusion and drug gradient, we can write the steady-state equation for the wild-type population density as  Figure 5. In shallow gradients, counter-flow gives rise to sharp cut-offs in the population density (ξ * c ) and establishment probability (ξ * u ) profiles (see Fig. 2a). The cut-off positions are well-captured by the analytical approximations (solid lines, eqs. 15 and 16). Squares are cut-offs extracted from simulation, dots and dotted line from numerical solutions to eq. 5 and 7 without and with a growth rate cut-off (SI section 5).
While eq. 11 does not have an analytical solution, we can change variables to ξ = ξm to see that the diffusion term in eq. 11 can be neglected when m 1 and thus c(ξ) ≈ s W T (ξ). More generally, for mutants in shallow gradients with a mutant net growth rate s M T (ξ) (see eq. 7), we get the eventual establishment probability of the familiar form This so-called quasi-static approximation is a straight-forward extrapolation of the wellmixed result and has been used to model the establishment probability by Greulich et al. [38]. Once solutions for c(ξ) and u(ξ) are found, they can be used to compute the rate of adaptation, R, and total death rate B. We find asymptotically and which agrees well with the numerical result in the two limiting cases λ and λ (Fig. 6). Thus, the treatment efficiency Q is constrained to a relatively small range in the absence of convection.

Co-flow and counter-flow
We now investigate the influence of convection on the population density c(x), the establishment probability u(x), and finally, the treatment efficiency Q, which, in the absence of convection, is naturally constrained to a relatively small range, see eqs. 13 and 14. We observed in simulations (Fig. 2) that convection will incur either a depletion or enrichment of both mutant and wild types in the antibiotic region, depending on whether the convection points towards higher or lower antibiotic concentrations (coor counter-flow, respectively). To study the effects of convection analytically, we set v = 0 in eqs. 5 and 7, and introduce the dimensionless flow speed ν = v/ √ Da 0 . The characteristic flow speed v c is closely related to the Fisher wave speed v F , which corresponds to the expansion speed of a freely growing population. If convection is too strong, i.e., the flow speed is great than the Fisher wave speed, (|ν| > 2), no steady-state population density exists because the population is washed away (see SI section 5), and we are therefore constrained to |ν| < 2.
In step-like gradients (m 1), even in the presence of convection, the population density c(ξ) approaches its asymptotic values exponentially fast, as can be seen by expanding eq. 5 to first order around its fixed points and solving the resulting linear differential equations (see SI section 4). As in the no-flow case, we can therefore approximate c(ξ) with a step function to find approximate solutions for u(ξ) far from the step gradient. With co-flow (ν > 0), the diffusion term in eq. 7 can be neglected to first order and hence the establishment probability has a very broad tail, u(ξ) ∼ (1 + ξ/ν) −1 (Fig. 4b, SI eq. S15). With counter-flow (ν < 0), the diffusive term is essential because the convection term cannot balance the non-linearity since both are negative. This gives rise to a rapid exponential decay u(ξ) ∼ e −νξ (Fig. 4b).
In shallow gradients (m 1), co-flow has little effect on the population density (and correspondingly establishment probability) profiles because the convection term νc (ξ) in eq. 5 is roughly proportional to m 1. In counter-flow condition, our model reproduces the cut-off in the population density c(ξ) seen in our simulations, which traces the net growth rate s W T (ξ) until it drops sharply. The position of this cut-off can be estimated from the definition of the Fisher wave speed ν F , which depends on the wild-type net growth rate s W T (ξ) (see also SI section 5). In shallow gradients, s W T (ξ) changes little around the cut-off, such that the (rescaled) Fisher wave speed becomes position-dependent ν F (ξ) = 2 s W T (ξ)/a 0 . Hence, no wild-type population density can be maintained in regions where the flow speed ν becomes larger than the local Fisher wave speed ν F (ξ), and solving for ξ gives the position of the cut-off Similarly, we find two cut-offs in the establishment probability profile u(ξ): since in shallow gradients, u(ξ) ≈ 1 − c(ξ), u(ξ) approaches 1 sharply when c(ξ) drops to zero at ξ * c . In addition, since the mutants have a local net growth rate that depends on the wild-type density, their speed limit is by ν F (ξ) = 2 1 − c(ξ), which gives another cut-off at the position  When ν < − √ 2, the two cut-offs reduce to one, and both c(ξ) and u(ξ) exhibit a sharp transition at ξ * c (see Figs. 2a). Fig. 5 compares the cut-off positions in the establishment probability profiles from the numerical model and simulations. While there is excellent agreement between this theoretical predictions, eqs. 15 and 16, and the numerical evaluation of the model, as well as with the cut-off ξ * u of the simulated establishment probability u(ξ), the cut-off ξ * c of the simulated population density c(ξ) appears shifted towards higher ν. This effect can be traced back to number fluctuation at the front of the population where the population size is small, which leads to a significant shift in the cut-off depending on the carrying capacity K (SI section 5, Fig. 5, black dotted line).
By integrating the numerical solutions to eqs. 5 and 7, we can compute the total death rate B and the rate of adaptation, R, and compare with the simulation results. As shown in Fig. 7, our model reproduces the phenomenology of the simulation very well, except in strong counter-flow conditions, where the number fluctuations at the front give rise to mostly quantitative differences.

Discussion
We presented stochastic simulations and an analytical model to study the fates of individual resistance mutations in spatial antibiotic gradients, for the first time also including the effects of convection. Our analytical model allows us to identify two characteristic length scales, λ and , which, together with the growth rate a 0 , determine the fate of resistance mutations (see Fig. 8). The gradient length λ describes the characteristic length scale over which the antibiotic concentration varies. Individuals travel a characteristic distance = D/a 0 through random dispersal with migration rate D before replicating. The ratio m = /λ determines whether a gradient is shallow or steep. If the gradient is shallow (λ ), the antibiotic concentration does not  Figure 7. Comparison between simulations and numerical solutions. (a) Total death rate B and rate of adaptation, R, computed from simulations (symbols) and numerical solutions to eqs. 5 and 7 (lines), in three different gradients (as in Fig. 3). The agreement is very good, except for strong counter-flow, where number fluctuations at the front become important, as explained in the main text. (b) The treatment efficiency Q found in simulations is also well-captured by our model, except for strong counter-flow.
appreciably change over distances accessible to a single individual in one generation, such that the population is locally adapted to the antibiotic concentration and changes only on length scales of order λ. If the gradient is steep (λ ), an individual close to the gradient can travel between regions of high and low antibiotic concentration in its lifetime such that the population density is constant far from the gradient and changes relatively abruptly over a distance around the gradient.
From the migration length , we also identify the characteristic speed v c = √ Da 0 = a 0 , which is closely related to the expansion speed v F = 2v c of a freely expanding population [48]. Hence, if convection is too strong, i.e., if the convection speed |v| > v F , the population will not be able to grow against the flow and be simply washed away ( Fig. 8; this akin to the extinction transition of populations due to convection [52,53]. For the emergence of resistance, we are thus naturally restricted to the region |v| < v F . It is illuminating to estimate the typical range of and v c for microbial communities. A typical (non-motile) bacterial cell may have a diameter of 1µm, swimming in a medium of viscosity comparable to that of water (e.g., blood [54]), which gives a diffusivity of order 0.1-1µm 2 /s. The motion of motile bacteria is characterized by much larger diffusivities, up to tens or even hundreds of µm 2 /s [55]. Together with a typical growth rate of 0.5−2hr −1 , this gives a possible range for between 50µm and several millimeters. For comparison, in a microfluidic experiment by Zhang et al. [26] the length scale on which the drug gradient varied was ≈ 200µm such that λ/ ∼ 1. For another recent experimental study by Baym et al. [12] with a reported spreading velocity of 40mm/hr, we find ≈ 1 − 20mm. Thus, depending on the properties of the bacteria and the antibiotic gradient, both shallow (λ ) and steep (λ ) gradients can plausibly arise. Similarly, with our range of parameters for D and a 0 , we find typical values for v c ≈ 0.1 − 10µm/s, which is on the same order of magnitude as estimated flow speeds in the gut, and achieved in an artificial gut microfluidic system [34].
We have quantified the emergence of resistance by computing the rate of adaptation for a wide range of gradient steepnesses and convection speeds. As a general rule, our model predicts that the rate of adaptationin shallow gradients is roughly proportional to the gradient length λ, as long as λ is smaller than the system size. This is in agreement with previous results from Greulich et al. [38], who employed the quasistatic approximation, eq. 12. As we have shown, this approximation breaks down when λ ≈ . In the opposite limit, λ , i.e., in steep gradients, we find a finite rate of adaptation, in agreement with the so-called staircase model of Hermsen et al. [37]. This is because mutants born in low antibiotic concentration regions can migrate into regions of high concentration, where they enjoy a big selective advantage. Our model serves as a generalization of these two previous model by capturing the full crossover from the "mutate and migrate" mechanism governing adaptation in steep gradients to the "local" adaptation in shallow gradients [56].
Convection can have a momentous influence on the rate of adaptation in our model. Co-flow (convection towards higher antibiotic concentration) increases the rate of adaptation by extending the range of the wild-type population further into high concentration regions, where resistance mutants have a strong advantage. At the same time, co-flow transports resistant mutants born in low-concentration regions (where the wild-type population density is high) to the front. In the same manner, counter-flow (convection towards lower antibiotic concentration) can restrict the range of the wildtype population and prohibit resistance mutation from establishing when they arise far away from the antibiotic gradient. This can lead to a decrease of the rate of adaptation by several orders of magnitude for strong counter-flow, until eventually the population goes extinct when counter-flow becomes too strong (see Fig. 8).
In the context of antibiotic treatment of an infection, the eradication of the infecting bacterial population is paramount [40] and therefore, we have argued that focusing on the rate of adaptation alone may be misleading. To quantify the trade-off between wild-type killing and the emergence of resistance, we measured the reduction in pathogenic growth relative to the rate of adaptation. This measure, which we call the treatment efficiency, has an intuitive meaning: a high treatment efficiency implies a strong reduction in the growth of the bacteria before a resistance mutation arises and establishes. In our model, the treatment efficiency is strongly affected by counterflow, where it can be an order of magnitude higher than in a no-flow scenario. This is ultimately a consequence of the rate of adaptation being doubly affected by convection, because convection alters both the wild-type population density (thus changing the local competition that resistant mutant clones face) and the dynamics of individual mutants, which are carried away from the high antibiotic concentrations they require in order to establish. Despite its intuitive meaning, however, the treatment efficiency as defined here makes no predictions about optimal treatment regimens; it merely serves as a quantification of the inherent trade-off between wild-type eradication and selection for resistance.
Our model is restricted to scenarios where interference between multiple mutant clones can be neglected, since we follow the establishment of individual mutations. Therefore, our model, operates exclusively in the "mutation-limited" regime [39], where the rate of adaptation is dominated by the waiting time until the establishment of mutations. By contrast, both Hermsen [37,39] and Greulich [38] consider the establishment of many, potentially contemporaneous mutations. In very shallow gradients, this can lead to a "dispersion-limited" regime if the mutational supply is large; the rate of adaptation is then dominated by the speed with which established mutations invade previously uninhabitable territory. We derive upper limits for the mutation rate µ per generation for our model to be applicable (SI section 6) and find that in steep gradients, clonal interference is negligible as long as µK 1, where K is the typical population size in a population patch of size . This result is similar to what is expected in well-mixed populations [57] (see SI section 6). In shallow gradient, the corresponding condition is µK ( /λ) 2 . Thus, if the gradient is too shallow for a given mutation rate, mutations enter the population too fast and clonal interference effects are expected. For resistance mutations with small target sizes, mutation rates can be very small, typically less than 10 −6 [58,59], and thus our approximation may be accurate even for relatively large carrying capacities K. Even when the mutation rate is high, we expect convection and spatial gradients to have the same qualitative effects on the establishment of resistance mutations, and thus our results should remain qualitatively correct even beyond the single-mutation regime considered here.
Our analytical model extends previous models [37][38][39] by incorporating arbitrary gradients and convection, and it can serve as a general branching process framework to study adaptation in spatially heterogeneous environments. For instance, we can easily incorporate antibiotic sanctuaries which are predicted to facilitate the emergence of resistance [23,25]. Other, more complex gradients may be necessary to describe cases in which the death rate depends on the antibiotic concentration on a non-linear manner [37]. Our model can also easily be re-interpreted to apply to a broad range of different ecological scenarios, like heterogeneous nutrient concentration. To make our model more realistic, it would be useful to model the bacteria as having a finite size, such that the population front can advance through mere growth, even against strong counter-flow [60]. For strong co-flow, individuals may also de-adhere and be carried away from the bulk population, thus founding extant colonies that enjoy large growth rates in the absence of competition for resources. Such processes can be studied by generalizing the diffusion term in our model to a long-range dispersal term as used frequently to model epidemics [61]. Since long-range dispersal can allow mutants far from the population front to escape the bulk population, we expect it to increase the total establishment probability and thus the rate of adaptation relative to short-range dispersal as discussed here. Another interesting generalization of the model would be to extend the model to two-dimensional populations since real biofilms typically grow as two-dimensional communities, with complex spatial patterns. The establishment of beneficial mutations in microbial colonies has recently been discussed [15,62]. Due to the particular strength of genetic drift at the front of such populations, beneficial mutations first have to reach a threshold size (depending on the strength of the selective advantage) neutrally before they become established. Once the mutant clone reaches the threshold size, the selective advantage of the mutants can deterministically drive them to fixation in the population. During the initial phase, the mutant clone is contained between boundaries with characteristic stochastic properties that are not captured in our one-dimensional model [14]. However, if the threshold size is small, the boundary fluctuations will not have a large impact on the growth of mutant clones. In such cases, we expect our results to apply also to two-dimensional populations.
The emergence of drug resistance remains a topic of significant interest, both from a scientific and a public health point of view. Considerable effort is brought forward to create novel antibiotics [63] and new therapy strategies are developed that attempt to limit the emergence of resistance [64,65], but more research is needed to understand how resistance evolves in complex spatio-temporal settings like the spatial gradients discussed in this paper. In particular, as we have shown here, convection constitutes an important factor in shaping the adaptation to antibiotics in spatial concentration gradients and should receive more attention from both theorists and experimentalists.

Individual-based simulation
We perform individual-based, stepping stone simulations where both wild-type and mutants are modeled explicitly. The population is divided into demes with carrying capacity K = 100 on a one-dimensional lattice. Wild-type and mutants replicate at a rate a 0 and migrate at rate D independently on their position. Wild type in deme x die at a rate b(x) = 1 + tanh(mx). Since we assume throughout that the antibiotic leaves the mutants unaffected, mutants do not die in our simulations. Analogously to a Gillespie algorithm, in each simulation step, a birth, death or migration event is performed according to its relative rate [44,45], as follows.
• Birth. Birth events occur at a total rate equal to a 0 x c(x), where c(x) is the total number of individuals in deme x. For each birth event, a source individual is selected at random and replicated into a random target site between 1 and K within the same deme. Because the target site can either already be filled with an individual or be empty, this move effectively translates into logistic growth within the deme.
• Death. In our model, only the wild type can die, thus deaths have a total rate corresponding to i b(x)c wt (x), where c wt (x) represents the number of wild type individuals in deme i. To perform a death event, first, a deme x is picked proportionally to its relative death rate b(x)c wt (x), and then, a random wild type within the same deme selected to be removed.
• Migration. Migrations are performed at a rate D x c(x) by picking a random individual and swapping it with a randomly selected target site from one of the two neighboring demes. As in the case of birth events, the target site can either correspond to an individual, or to an empty site.
• Time step. Time is tracked by sampling a time interval δt from an exponential distribution with rate x [(a 0 + D)c(x) + b(x)c wt (x)], as in a standard Gillespie algorithm. The total elapsed time is the sum of the sampled time intervals.
• Convection. Convection (with convection speed v) is implemented by shifting the population by one deme away or towards the antibiotic gradient, for negative or positive convection respectively. The convection strength is controlled by performing the shift at a rate 1/|v|, i.e., any time the time elapsed since the last shift is greater than 1/|v|.
For each simulation, we first allow the wild type to reach the steady-state profile c(x). We then introduce one mutant element at position x and run the simulation until either all mutants go extinct, or mutants reach the last deme in the simulation box. No further mutations are allowed in the course of the simulations. The probability of fixation u(x) is then computed as the proportion of the simulations in which a mutant introduced at x reached fixation.

Numerics
Numerical solutions to eqs. 5 and 7 were obtained by evaluating the differential equations using Mathematica's built-in NDSolve routine with the backwards differentiation (BDF) method, with a maximum step size of 1 and a domain size of 1000. Initial conditions were chosen according to the analytical approximations to the steady-state profiles given in the text. This is done to speed up computation and increase numerical stability, but is otherwise inconsequential; starting with different initial conditions leads the same final solution. To obtain steady-state profiles, we solve the full time-dependent problems until the solution no longer changes for longer evaluation times.
To compute the establishment probability u(x) in realistic wild-type population profiles c(x), we first computed the steady-state population density and then used the final profile to compute the (constant in time) local death rate for the mutants b(x) = 1 − c(x). The resulting numerical solutions were integrated numerically using Mathematica's built-in NIntegrate routine to obtain the total death rate B and the rate of adaptation, R, .

S1. Antibiotic concentration profile
In the main text, we suppose an antibiotic gradient varying over a characteristic length scale λ. Here, we present an example for how such an antibiotic gradient may arise. We use a simple reaction-diffusion model for the antibiotic concentration C in a onedimensional system of length L to show that our model for the antibiotic gradient, eq. 1 in the main text, can approximately capture a variety of realistic gradient profiles.
Consider an antibiotic source at position −L/2 and an antibiotic sink at position L/2. The antibiotic concentration at the source is c ∞ and the antibiotic concentration at  Figure S1. Example of realistic growth rate profiles for different convection speeds v C (v C > 0, cyan; v C = 0, black; v C < 0, magenta). Solid lines are solutions to eq. S1, dashed lines are fits with eq. 1 in the main text.
the sink is 0. The drug diffuses with diffusion constant D C and is subject to convection with speed v C . Additionally, the antibiotic may be degraded, e.g., by cells metabolizing the drug or through chemical degradation, at a rate γ. Under these assumption, the steady-state antibiotic concentration profile is described by the reaction diffusion equation The resulting net growth rate s(x) = 1 − 2C(x)/C ∞ is shown for three values of v C and γ = 0 in Fig. S1 (solid lines). The dashed lines are approximations of the form given in eq. 1 in the main text. For our purposes, eq. 1 gives a good enough fit to (simple) realistic antibiotic profiles.

S2. Derivation of survival probability
Let u x (t) denote the probability that a mutation born at lattice site x survives for time t. Denote by a(x) and b(x) the local birth and death rate, respectively, and let v + and v − be the rates to migrate a distance δx (i.e., one lattice site) to the right and to the left, respectively. Then, the u(x, t) after a short time interval satisfies the equation The first term on the right-hand side accounts for the fact that when the initial mutant divides then there are two mutants, and the probability of survival of at least one lineage is 1 minus the square of both lineages disappearing. The second term describes the probability of migrating one lattice site to the left or right, respectively. The first term describes the case of nothing happening in the time interval .
Letting → 0 and performing the Taylor expansion in δx, we obtain eq. 5 in the main text where Note that the sign of the convection term v∂ x u is reversed from that in the equation governing the population density, eq. 5 in the main text.
From the profiles for c(ξ), the total death rate B can directly be computed as B = −3 + 9 + √ 6 ≈ 0.38, as quoted in the main text. Since c(ξ) transition sharply around ξ = 0, the establishment score R = ∞ −∞ c(ξ)u(ξ) can be approximated as 3 ≈ 1.86, which is very close to the numerical result R ≈ 1.91.

S4.
Step-like antibiotic concentration profile with flow In the presence of flow, the equations for the population density becomes For ξ 0, the population density is small such that we can extract the asymptotic behavior from the linear equation. With the boundary condition c(∞) = 0, we get Similarly, for ξ 0, c(ξ) is close to one; introducing Ψ(ξ) = 1 − c(ξ) and linearizing the result equation, we get that Thus, even in the presence of flow, the population density transitions sharply around the antibiotic concentration step. Therefore, we can employ the same approximation as in the no-flow case and approximate c(ξ) ≈ Θ(−ξ). The establishment probability u(ξ) thus approximately satisfies the equation For ξ > 0, the establishment probability approaches 1 exponentially quickly. To see this, we can again expand around 1 and solve the resulting linearized equation. The result is  The results in summarized in Fig. S2. For ξ < 0, the linear term in eq. S15 vanishes, and we have to distinguish two cases, ν > 0 and ν < 0. For ν > 0, we expect a broad profile because co-flow facilitates the establishment of mutations born deep inside the population bulk. Hence, the non-linear term in eq. S15 cannot be neglected. On the other hand, if the co-flow is strong, we can neglect the diffusive term and find asymptotically u(ξ < 0) ∼ 1 in good agreement with the numerical solution (cyan line in Fig. 4b in the main text). For counter-flow, i.e., for ν < 0, it is intuitively clear that the establishment probability must vanish faster than in the no-flow case. Thus, we expect for that u(ξ 0) such that we can neglect the non-linearity. As a result, we get u(ξ < 0) ∼ e −νξ , ν < 0, which also agrees well with the numerical solution (red line in Fig. 4b in the main text).

S5. Cut-off position and finite carrying capacity
Here, we provide a more rigorous derivation of the cut-off position ξ * c . For ν < 0 and m 1, the population density c(ξ, τ ) satisfies the equation Performing the transformation ϕ = ce −νξ/2 we obtain Assume the population is started with an initial population at ξ = −∞. To arrive at the steady-state population density discussed in the main text, ∂ τ ϕ must be positive initially and go to zero as τ → ∞. For the right-hand side of eq. S21 to be positive, we must have s W T (ξ) − ν 2 4 > 0, whence we obtain eq. 15 in the main text. Eq. 17 is obtained by the same arguments. Note that this derivation also gives the general condition ν 2 < 4, i.e., the external flow speed must be smaller than the Fisher wave speed of the population.
We model the effects of a finite carrying capacity K by introducing a cut-off in the net growth rate if the local wild-type population density becomes lower than 1/K. To this end, we replace s W T (ξ) in eq. 5 in the main text with s W T (ξ, K) = s W T (ξ)Θ(c(ξ) − 1/K). (S22) The results are shown in Fig. S3: for fixed ν, smaller K reduce the maximal sustainable local flow speed, such that the cut-off is shifted to the left for smaller K (Fig. S3a). When K is relatively small, the cut-offs can be shifted strongly, especially at strong flow speeds (Fig. S3b).
For an analytical approximation, we use the result that a finite carrying capacity reduced the Fisher wave speed by a factor 1 − π 2 /2(ln K) 2 . Plugging this correction into eq. 16 we get a corrected cut-off position ξ * c , which shows reasonable agreement with our simulations (Fig. S3c, black line).

S6. Small mutation rate estimate
For clonal interference to be negligible, we estimate that the time to fixation of a mutant clone must be much faster than the time it takes for a mutation to arise and establish. In a well-mixed population, this condition can be phrased as µN 1, or, more precisely, τ fix ≈ ln[N s]/s τ establish ≈ 1/(µN s).
In a spatial scenario, the establishment time τ establish is given through a generalization of the well-mixed result, i.e., where we have made the dependence on the carrying capacity K explicit, and the mutation rate µ is given per generation. We have found in the main text that R ∼ in steep gradients and R ∼ λ in shallow gradients. We can estimate the time to fixation by considering a mutation born at a typical distance L from the front of the population; it then has to travel this distance L to take over the whole population. In our model, the distance L can be estimated from the position of the peak of the product u(ξ)c(ξ). This is approximately given by The time it takes to travel this distance depends on the steepness of the gradient: in a steep gradient, mutants have to migrate into region of high antibiotic concentration through random dispersal. We can model their motion as a random walk, and hence, the distance = √ 2Dτ fix such that τ fix = 2 /2D. This gives the criterion for negligible clonal interference as µK 1.
In shallow gradients, the time to fix for a mutant clone is given by the length L divided by the speed v with which it travels, which is roughly v ∼ √ Da 0 . This gives the condition µK 2 λ 2 . (S29)