Bacterial motility can govern the dynamics of antibiotic resistance evolution

Spatial heterogeneity in antibiotic concentrations is thought to accelerate the evolution of antibiotic resistance, but current theory and experiments have overlooked the effect of cell motility on bacterial adaptation. Here, we study bacterial evolution in antibiotic landscapes with a quantitative model where bacteria evolve under the stochastic processes of proliferation, death, mutation and migration. Numerical and analytical results show that cell motility can both accelerate and decelerate bacterial adaptation by affecting the degree of genotypic mixing and ecological competition. Moreover, we find that for sufficiently high rates, cell motility can limit bacterial survival, and we derive conditions for all these regimes. Similar patterns are observed in more complex scenarios, namely where bacteria can bias their motion in chemical gradients (chemotaxis) or switch between motility phenotypes either stochastically or in a density-dependent manner. Overall, our work reveals limits to bacterial adaptation in antibiotic landscapes that are set by cell motility.

Line 39: Perhaps a little subjective, but this line sounds a little overly dramatic to me.I therefore initially took the manuscript a little less serious.(feel free to keep it in, but I am providing a reader perspective here) 'Yet, despite these important realizations, there is surprisingly little evidence about how motility contributes to bacterial adaptation, namely to the evolution of antibiotic resistance, which is ranked among other major threats such as climate change and terrorism.'Line 162: 'At high motility (ν > δ, Supplementary Movie 1), the wild-type is spread across space (SI Theorem 3), competing with mutants for space even in regions where the wild-type cannot grow.' I do not fully understand this, yet this is important for much of the manuscript.Do the authors mean that at these regions the wild-type cannot grow, but can be maintained and therefore contributes to K, and therefore is able to compete for space with mutants?Can the authors explain/ discuss? Figure 1b: Blue and red are only explained at panel b, but are already described at panel a.Therefore suggest to move up to panel a.Line 247: Replace in 'this' case for in 'the latter' case for clarification.
Line 323: 'This result stems from the fact that the slower phenotype is naturally selected for in the absence of phenotypic switching, given that the slower phenotype is exposed to antibiotics less often and maximises the local fitness of the wild-type, and stochastic switching levels out the numbers of both phenotypes' Clarification needed: Why would the slower phenotype be less exposed to antibiotics?Line 368: 'As a result, the population has decreased antibiotic exposure and adaptation rate.'Yet, if motility would be sufficiently fast, then the antibiotic exposure may be averaged at short timescale.Correct?(I do understand that average-step sizes are not in the staircase model) but would that change dynamics by a different (likely larger) mutational supply, and (likely lower) cost of resistance ?The authors refer to these aspects in the last part of the discussion, and I wonder whether time averaged concentrations would therefore change the effect of motility as described in line 368.
Reviewer #2: Remarks to the Author: Piskovsky and Oliveira used a stochastic mathematical model to study the interaction of cell motility and bacterial adaptation in a spatially heterogeneous environment.In particular, they use an adapted version of the staircase model of Hermsen et al (a lattice of genotypes and spaces with increasing antibiotic concentrations) and show that increasing the migration rate drastically changes the dynamics of the original model, since not only evolutionary time (defined as the earliest time a more resistance mutant is produced) but also ecological time (defined as the earliest time when a mutant population becomes larger than the wild-type in the overlap region) becomes relevant for adaptation.They further depict how other characteristics of active bacterial motility, i.e. phenotype switching between high and low motility and chemotaxis, can change adaptive rates in more complex models.They describe 3 different motility regimes and show that high motility can lead to population extinction in a sink-like environment.They conclude that motility can also limit adaptation, an effect which has been overseen in the original staircase model.
The modeling/analysis in this work is of high quality, and its predictions could be a good foundation to assess the effect of motility on resistance evolution experimentally.However the biological relevance of the model remains unclear.In addition, some of the outcomes of the analysis are somehow trivial.

Major points
1.It is unclear what the relationship between ν and δ represents in the physical world and when, in practice, should we expect transitions between high and low motility scenarios.What are the relevant units (in the real world) and how can they be compared, especially in a continuous space?How to know (in a real experiment) which regime applies?
2. The simulation results need error bars, e.g. in Figure 2. It is unclear when the decrease/increase in adaptation rate is actually significant, especially when analytics predict a steep reduction.
3. When evaluating the effect of cell-density dependent switching on adaptation rates, the authors consider all cells in a single compartment to switch their phenotype.This leads to dynamics where swarming populations (in the fast to slow scenario) might not switch phenotype ever again and just spread into all the compartments at low density, which is a rather trivial outcome.Would the author's conclusions still hold true in a model where the random switching rate s is expressed as a function of cell density N, and where not all cells but only a subpopulation in a compartment undergoes this switch (similarly to Fig. 3).This would also much better reflect phenotypic heterogeneity in biological systems, where both phenotypes usually coexist in the same space.4. Growth in each compartment is dependent on the presence of all cells in this compartment, whether they are dividing or not.In nature, not growing cells might not compete fully for resources with growing cells.The authors have taken this assumption from the original staircase model, which is adequate for showing differences with the original model, but to make the current model biologically more meaningful it would be important to explore how differential competition between growing and non growing cells might interfere with the dynamics they have described.
5. The title should indicate that this work is exclusively numerical/theoretical.Minor points 6.Several important parameters, such as the mutation rate (which is set to be different by 3 magnitudes between increase and decrease of resistance without justification) are identical to the original work by Hermsen et al..While this is adequate to show differences with the original staircase model, it would be good if they could at least comment on this choice, and in the best case also look at its interaction with motility.
7. The authors should comment on the ratio of different parameters that can lead to the different motility regimes they have described, and how these ratios relate to data from biological/experimental systems, and if they are biologically meaningful.8. Mutant fitness (calculated as r(1-Nx/K)-δ) is actually expected growth dependent on cell density, and does in my opinion not really correspond to the usage of this term in experimental work.The fitness gradient therefore rather represents the steepness of the population distribution.Another term or a comment in the text would be needed to make this less misleading.9.The graphical representations of the model, especially the depiction of the heat maps (usually in Fig. Xb on the left) would benefit from more clarity.x and y scales should be uniformous within a figure.The comments on Figure 5 ("increase time" and "increase motility") in combination with arrows are misleading for the understanding of the Figure .Reviewer #3: Remarks to the Author: # Review: Piskovsky & Oliveira This manuscript reports results on (versions of) a mathematical / computational model.It aims to study the effects of several aspects of motility on the rate of antibiotic resistance evolution in an environment that contains a concentration gradient of an antibiotic.
The model used is based on the so-called Staircase Model published in 2012 in PNAS (Hermsen *et al*), which assumes a one-dimensional environment that is subdivided into compartments with increasing drug concentrations.In the article, the analysis of the model was restricted to a particular parameter regime: it assumed that the carrying capacity was large, the migration rate between neighboring compartments low, and mutation rates lower still.The current manuscript expands this analysis by studying the regime of fast migration/motility.Next, the authors expand the model by including stochastic switching between a fast and a slow motility state, or a density-dependent motility rate.
As I see it, the main findings are presented in the first four figures.Figure 1 demonstrates that that the behavior of the high-motility regime is qualitatively different from that of the low-motility regime: it results in smoother density profiles, which reduces the growth rate of mutants with a higher level of resistant emerging in the system and hence increases the "ecological waiting time".Figure 2 shows that, as a consequence, the adaptation rate depends non-monotonously on the motility rate, with a positive dependence in the low motility regime but a negative dependence in the high-motility regime.Figure 3 presents results of the model in which individual cells switch stochastically between two motility rates.The main point is that this extended model behaves very similarly to the original model with some effective migration rate.The same holds for the model with density-dependent motility, as presented in Figure 4.
Apart from these main points presented in full in the main text, quite a lot of additional work has been shifted to supplementary figures, where the effect of a fitness cost, chemotaxis, HGT, different types of antibiotics, etc. are also examined, and to the supplementary text, which contains detailed mathematical derivations and analysis.
Generally, I think the analysis is sound and the conclusions seem warranted.The text is clearly and attractively written and the figures are carefully designed.I only have minor comments.
## Comments on content 1. Around the same time when the Staircase Model was published, Greulich et al (PRL, 2012) published a similar model.An important difference between Hermsen et al and Greulich et al was that Hermsen et al focused on a series of relatively isolated compartments, whereas Greulich et al considered continuous space.It seems to me that this difference disappears in the high-motility regime.While the authors cite Greulich et al, they do not compare their results to that paper in much detail.I think they should.2. Section D and Figure 5 of the Results section discusses the effect of migration on the density profile from a dynamic systems perspective.Although the mathematics presented here and in the SI are fun (for some readers), the ultimate biological messages are very intuitive and can easily be summarized in a few sentences, as the authors in fact do in the last paragraph of this section.I would advice to move this part to the supplement.3. How did the authors calculate the adaptation rates of, say, Fig. 2b in their simulations?If I understand well, they used the equation $a_R = 1/(\mathbb{E}T_{evo}) + \mathbb{E}T_{eco})$ and measured $T_{evo}$ and $T_{eco}$ in their simulations.I would argue it makes more sense to measure the actual time it took in the simulations to take a certain number of steps.The authors argue that this should equal the equation above, but this is only approximately so.In the current setup, this assumption remains untested.4. In introduction and discussion you claim that "cell motility has largely been overlooked".I think this is claim is a bit dubious given that all mathematical models that you cite do include motility, albeit at a fixed and relatively low rate.Perhaps you can make this claim more specific.## Comments on figures 1.In several figures, genotype is indicated with a blue or yellow color scale.Even though genotype is restricted to integers, the legend belonging to the color scale shows integers at the *border* between the colors.This is confusing.2. In the three panels of Fig. 3b, the lines indicating the borders between the regimes seem to be inconsistent.In particular, in the heat maps the horizontal and vertical black lines are not at $\nu_i = \delta = 10^{-1}$.3. I would prefer that the parameters used in the figures are mentioned in the figure caption.

Comments from reviewers REVIEWER 1
Review Bacterial motility as a driver of antibiotic resistance evolution by Vit Piskovsky and Nuno M. Oliveira.This manuscript describes the e↵ect of motility of bacterial cells in the adaptation to spatial gradients of antibiotics.
It builds on empirical and theoretical literature on the spatial evolution of antibiotic resistance, and extends that with an investigation of (heterogeneous) bacterial motility.The findings are relevant and significant, because motility/ migration is relevant in general for a wide range of (infectious) (micro)organisms, and in this context specifically for the rate of adaptation to antibiotics.
The work extensively investigates e↵ect of several relevant motility scenarios, and even extends those to match other potentially relevant biological aspects, such as HGT.The addition of these scenarios is in line with the (main) conclusion of this manuscript.
As far as I can judge the methods are su ciently well described and the methodology is sound.I therefore only have very few minor comments on this manuscript.

A1. Comments
Line 39: Perhaps a little subjective, but this line sounds a little overly dramatic to me.I therefore initially took the manuscript a little less serious.(feel free to keep it in, but I am providing a reader perspective here) 'Yet, despite these important realizations, there is surprisingly little evidence about how motility contributes to bacterial adaptation, namely to the evolution of antibiotic resistance, which is ranked among other major threats such as climate change and terrorism.'This paragraph is now rephrased in the revised manuscript to "Despite these important realizations, we know little about how cell motility contributes to bacterial adaptation, namely to the evolution of antibiotic resistance, which is a major public health concern [8]."Line 162: 'At high motility (⌫ > , Supplementary Movie 1), the wild-type is spread across space (SI Theorem 3), competing with mutants for space even in regions where the wild-type cannot grow.'I do not fully understand this, yet this is important for much of the manuscript.Do the authors mean that at these regions the wild-type cannot grow, but can be maintained and therefore contributes to K, and therefore is able to compete for space with mutants?Can the authors explain/discuss?
In our model, we assume logistic growth in each compartment.Importantly, even in regions where wild-type cells cannot grow, they can be continually imported from neighbouring regions and, consequently, occupy space and contribute towards the carrying capacity K.As a result, the presence of wild-type cells decreases the growth of mutant cells.We have modified line 162 so that this idea is now clearer.As another reviewer also asked about this issue, and in particular about the plausibility of assuming equal competitive ability between growing and non-growing cells, we have now considered the e↵ect of non-growing cells that compete at a reduced rate with growing cells (new Supplementary Text S4, new Supplementary Fig. 4, shown here as Fig. A1).Fig. A1 shows that di↵erential competition between growing and non-growing cells does not a↵ect the key conclusions of our work.We now mention this e↵ect in Results section B, and in the Discussion.We had a typo in the caption of the figure, namely (black dots) instead of (blue and red dots).To clarify, dots in panel a represent cells, and we do not explain the meaning of specific colours.The meaning of the colours (i.e., why we distinguish between blue and red dots) only becomes important when we explain the adaptation process.For this reason, we postpone the explanation of the colour scheme to panel b of the figure.| Di↵erential competition between growing and non-growing cells.In our original model, growing and non-growing cells are assumed to contribute equally to the carrying capacity.Here, we consider the e↵ect of non-growing cells contributing to the carrying capacity at a reduced rate ↵ 2 [0, 1] when compared to growing cells (↵ = 0 corresponds to no contribution, ↵ = 1 corresponds to the same contribution as in the original model).Therefore, the division rate of growing cells (g x) is max(0, r(1 (N g x + ↵N n x )/K), where N g x is the number of growing cells above the staircase (g x) and N n x is the number of non-growing cells below the staircase (g < x) in a given spatial compartment x.The heatmap of adaptation rate on the (⌫, ↵) plane explains how the reduction in competitive strength of non-growing cells ↵ changes the relationship between motility rate ⌫ and adaptation rate aR.The data shows that there is no significant change in the relationship between adaptation rate and motility.A detailed description of this model can be found in Supplementary Text S4.
Line 247: Replace in 'this' case for in 'the latter' case for clarification. Replaced.
Line 323: 'This result stems from the fact that the slower phenotype is naturally selected for in the absence of phenotypic switching, given that the slower phenotype is exposed to antibiotics less often and maximises the local fitness of the wild-type, and stochastic switching levels out the numbers of both phenotypes' Clarification needed: Why would the slower phenotype be less exposed to antibiotics?
The slower phenotype is less exposed to antibiotics because slower cells are less likely to migrate into regions with large antibiotic concentrations, and their population will be concentrated in regions where they can divide.This argument is consistent with the theory developed by Alan Hastings in [A1], where the author shows that spatial heterogeneity alone (such as the one introduced by a stable antibiotic gradient) selects against the evolution of dispersal/migration/motility.In particular, the author argues that slower phenotypes are better at staying in regions with high fitness, such as low antibiotic concentration in our case.Moreover, we note that our equation for the wildtype profile (S16) is of the same form of equation ( 1) in Hastings's paper if the switching rate s is 0. Also, our Fig.1d shows that, at lower motility, the population is concentrated in compartments where it can divide.In this sense, the slower phenotype is less exposed to antibiotics.We modified the original line 323 to make this idea clearer: "This result stems from the fact that the slower phenotype is naturally selected for in the absence of phenotypic switching, given that the slower phenotype spends more time in regions of low antibiotic concentration where it can proliferate, ...".
Line 368: 'As a result, the population has decreased antibiotic exposure and adaptation rate.'Yet, if motility would be su ciently fast, then the antibiotic exposure may be averaged at short time-scale.Correct?(I do understand that average-step sizes are not in the staircase model) but would that change dynamics by a di↵erent (likely larger) mutational supply, and (likely lower) cost of resistance?The authors refer to these aspects in the last part of the discussion, and I wonder whether time averaged concentrations would therefore change the e↵ect of motility as described in line 368.
To clarify, line 368 was used in the context of density-dependent motility, with cells moving fast at low density (rate ⌫ L ) and slowly at high density (rate ⌫ H ). In this case, the vast majority of the population reaches high density, moves slowly and resides in low antibiotic concentrations (Fig. 4d, lower right panel).Therefore, the time-averaged antibiotic concentrations of the entire population are approximated by the time-averaged antibiotic concentrations of the high-density phenotype, which are low and independent of the time-averaged antibiotic concentrations of the few fast cells.This pattern is reflected in the idea of e↵ective motility, where the system behaves as if only the high-density phenotype was present (with motility ⌫ H ). It should be noted also that line 368 is used in a very specific context, namely to highlight exceptions to the rule of e↵ective motility: "First, when cells switch from fast to slow above the threshold S, the population becomes very stable as the fast cells at the front quickly return to the slow population bulk (Fig. 4c).As a result, the population has decreased antibiotic exposure and adaptation rate (Fig. 4b)."Therefore, line 368 simply meant to highlight that the population adapts slower than a population with single motility ⌫ H in the original staircase model.Regarding the e↵ects of resistance costs and mutational supply, if we understood correctly, the reviewer asks if an increased mutation rate µ f and decreased cost of resistance c can change this comparison between the original model and the model with density-dependent motility.To clarify this issue, we would like to note that in the density-dependent motility model, we assumed c = 0 and thus resistance costs cannot be further reduced.Moreover, the mutation rate does not a↵ect the reduction in cell number at the population front of the wild-type profile, which can be concluded from the absence of µ f parameter in the wild-type profile equations (S6) and (S29).Therefore, we believe that reducing c and or increasing µ f do not change the idea expressed in our original line 368.
We would like to thank the reviewer for this question because it made us realize one aspect that we did not consider originally.While we studied the e↵ects of resistance costs (Supplementary Fig. 2), this question made us wonder about the e↵ects of mutation rates µ f .We have now varied the mutation rate µ f (new Supplementary Fig. 3, Fig. A2 here) and showed that our key conclusions are not a↵ected by it for realistic values of mutation rate, and that the low, high and deadly motility regimes capture well the evolutionary dynamics.Piskovsky and Oliveira used a stochastic mathematical model to study the interaction of cell motility and bacterial adaptation in a spatially heterogeneous environment.In particular, they use an adapted version of the staircase model of Hermsen et al (a lattice of genotypes and spaces with increasing antibiotic concentrations) and show that increasing the migration rate drastically changes the dynamics of the original model, since not only evolutionary time (defined as the earliest time a more resistance mutant is produced) but also ecological time (defined as the earliest time when a mutant population becomes larger than the wild-type in the overlap region) becomes relevant for adaptation.
They further depict how other characteristics of active bacterial motility, i.e. phenotype switching between high and low motility and chemotaxis, can change adaptive rates in more complex models.They describe 3 di↵erent motility regimes and show that high motility can lead to population extinction in a sink-like environment.They conclude that motility can also limit adaptation, an e↵ect which has been overseen in the original staircase model.
The modeling/analysis in this work is of high quality, and its predictions could be a good foundation to assess the e↵ect of motility on resistance evolution experimentally.However the biological relevance of the model remains unclear.In addition, some of the outcomes of the analysis are somehow trivial.and how can they be compared, especially in a continuous space?How to know (in a real experiment) which regime applies?
Our main purpose was to develop theory on how cell motility a↵ects bacterial adaptation in spatially-heterogeneous environments and we use the evolution of antibiotic resistance as a case study.While our goal was to provide a general theory that was not meant to recapitulate any experimental setup in particular, we understand the concern of the reviewer and the importance of relating our model results to experimental studies.We now created a new supplementary section in the paper that addresses this issue (new Supplementary Text S7).In particular, we now considered what we call a visiting number, a non-dimensional number that captures the antibiotic variability experienced by an average cell during its lifetime.Specifically, the visiting number V is defined as the number of regions that di↵er in antibiotic concentrations (on MIC scale, MIC = minimal inhibitory concentration) and that are visited by an average cell during its lifetime.This number can be computed for both our model and di↵erent experimental setups, and we now examine the visiting number in our model and in two works, [A5] and [A6], that studied experimentally the evolution of antibiotic resistance in heterogeneous environments.In our model, the visiting number is V = ⌫/ as ⌫ is the number of compartments visited per unit of time and 1/ is the timescale of bacterial lifespan.In experiments, , where v is the characteristic cell speed, t is the doubling time and l is the length-scale over which drug concentrations vary on MIC scales.This expression follows as v/l estimates the number of di↵erent MIC scales visited per unit of time and t characterises the timescale of bacterial lifespan.The visiting number can be used to identify the adaptation regime since Fig. 2b shows that the low motility regime corresponds to ⌫ < (i.e., V < 1) and the high motility regime to ⌫ > (i.e., V > 1).To relate our model to experiments, we computed the visiting numbers in the experiments [A5, A6] and noted the predicted motility regime in Table I (new Supplementary Table I).To further confirm such prediction, we compared the qualitative features of the experimental evolutionary dynamics with the results of our model in Fig. 1d.In [A5] (V = 120), the resident mutants quickly invaded the entire environment, which is characteristic of the high motility regime in Fig. 1d.In [A6] (V = 0.1846), the adaptation dynamics was located at the population front and multiple strains with di↵erent resistance levels coexist, which is characteristic of the low motility regime in Fig. 1d.This comparison suggests that the key results of our model can be compared with experimental studies.Notably, our new analysis predicts that bacteria in [A5] and [A6] were evolving in very di↵erent motility regimes, which may explain why the authors observed di↵erent adaptation rates in their studies.Find more details in Discussion and in Supplementary Text S7.
2. The simulation results need error bars, e.g. in Figure 2. It is unclear when the decrease/increase in adaptation rate is actually significant, especially when analytics predict a steep reduction.
We now added error bars to Fig. 2b and showed that they are smaller than the symbols in Fig. 2b, see Fig. A3, consistent with what was found in [A11].We now clarified this observation in the caption of Fig. 2. To clarify, the central limit theorem guarantees a small size of the errors for a su ciently large number of simulations n (see a detailed description of computation below).In . The equation for estimating the adaptation rate ( 5) is Defining a R in this way, we can notice that the central limit theorem implies that Using the delta method for g(x) = x 1 , we can see that Therefore, we can use the sample moments Using this formula, we plotted the standard deviation as error bars centred at the symbols that denote the estimated mean value of the adaptation rate a R (Fig. A3).We can see that the error bars are so small that they are only visible after su cient zooming into the figure, which is ensured by taking su ciently large n = 50.The standard deviations associated with estimating the adaptation rate aR by finite number n = 50 of simulations can be computed using the central limit theorem.We plotted the standard deviation as error bars centred at the symbols that denote the estimated mean value of the adaptation rate aR and that are used in Fig. 2b.
3. When evaluating the e↵ect of cell-density dependent switching on adaptation rates, the authors consider all cells in a single compartment to switch their phenotype.This leads to dynamics where swarming populations (in the fast to slow scenario) might not switch phenotype ever again and just spread into all the compartments at low density, which is a rather trivial outcome.Would the author's conclusions still hold true in a model where the random switching rate s is expressed as a function of cell density N, and where not all cells but only a subpopulation in a compartment undergoes this switch (similarly to Fig. 3).This would also much better reflect phenotypic heterogeneity in biological systems, where both phenotypes usually coexist in the same space.
For simplicity, we modelled phenotypes implicitly (i.e., all cells in a given compartment have the same phenotype).
One of the advantages of this modelling is that it reduces computational time.However, to address the comment of the reviewer, we now show that our model can be derived as a limit of a model suggested by the reviewer, where both phenotypes are explicitly tracked in each compartment (Fig. A4a).In particular, we show that the evolutionary dynamics is governed by the same e↵ective motility in both models and our key conclusions hold as originally presented (Fig. A4b-f).The alternative modelling technique is now described in Results section C, new Supplementary Fig. 8 and new Supplementary Text S6.For completeness, we give below a detailed account of the correspondence between the two models of density-dependent switching (Supplementary Text S6).
Details on the correspondence between the models.First, we construct the explicit model suggested by the reviewer.
In the explicit model, we consider the same setup as in the model for stochastic switching (Fig. 3a) where two phenotypes of di↵erent motility ⌫ 1,2 are considered and switch at rate s.To introduce density-dependent motility, we modify the switching rates s between the motility phenotypes P (Fig. A4a).In particular, at a given spatial position x with N x cells, phenotype 1 switches into phenotype 2 at a rate and phenotype 2 switches into phenotype 1 at a rate The total switching rate s 1!2 + s 2!1 = 2s is kept constant and the phenotypes switch -times more likely in the direction 1 ! 2 than 2 ! 1 at low density (N x < S), while they switch -times more likely in the direction 2 ! 1 than 1 ! 2 at high density (N x > S).To identify phenotype 1 with the low-density phenotype (⌫ 1 = ⌫ L ) and phenotype 2 with the high-density phenotype (⌫ 2 = ⌫ H ), we restrict to  1.Moreover, for switching to be important, we assume that s .We can notice that the implicit model of density-dependent switching corresponds to s ! 1 and !0. Therefore, the implicit model is expected to exhibit the same dynamics as the explicit model.We show this correspondence in Fig. A4b-f.The adaptation regime is again governed by the same e↵ective motility ⌫ ± = ⌫ H/L (Fig. A4b-e) and the wild-type profiles are very similar and di↵er only in predicted phenotypic proportions at individual spatial positions between the explicit and implicit models (Fig. A4f).Both models give rise to the same mutant fitness (red bars), which further explains the similarity of the evolutionary dynamics in the two alternative models.(c,e) Adaptation in source-like environment.Adaptation rate heatmap on the (⌫L, ⌫H ) plane for di↵erent switching thresholds S in the source-like environment, where ⌫L,H are the low-density/high-density motility rates.The (⌫L, ⌫H ) plane can be partitioned into di↵erent combinations of adaptation regimes and its diagonal corresponds to a population of a single motility.The diagonal separates slow-to-fast and fast-to-slow switching combinations, which di↵er in relative motility at low-to-high density.Bacteria generically adapt as if all cells had the same e↵ective motility, which corresponds to the intersection of the level sets (white dashed lines) with the diagonal.This e↵ective motility generically matches the low-density (resp.high-density) motility at high (resp.low) threshold S. At low threshold S, this generic rule has one exception in the implicit model.The adaptation rate is reduced in the mixed motility combination of the fast-to-slow switching case, as all cells in the overlap region switch to fast phenotype and are unlikely to start a growing mutant colony there.(b,d) Adaptation in sink-like environment.Same as panels (c,e), but for a sink-like environment.Importantly, the same e↵ective motility governs the adaptation dynamics as in panels (c,e).At low threshold S, this generic rule of e↵ective motility has three exceptions.First, the adaptation rate is reduced in the mixed motility combination of fast-to-slow switching.Second, the adaptation rate is reduced for slow-to-fast switching when the environment is sink-like and the high-density motility phenotype moves above the critical motility.Third, a deadly motility regime (grey) can only occur if the initial population is of low density and the low-density motility phenotype moves above the critical motility.In the explicit model, the critical motility is slightly increased as the explicit presence of a few slow cells with a high-density phenotype can prevent extinction.(f ) Wild-type profiles in the implicit and explicit model.The wild-type profiles are similar and di↵er only in predicted phenotypic proportions.Note that the mutant fitness (red bars) is similar in both models, which further explains why both models exhibit the same adaptation dynamics.
4. Growth in each compartment is dependent on the presence of all cells in this compartment, whether they are dividing or not.In nature, not growing cells might not compete fully for resources with growing cells.The authors have taken this assumption from the original staircase model, which is adequate for showing di↵erences with the original model, but to make the current model biologically more meaningful it would be important to explore how di↵erential competition between growing and non-growing cells might interfere with the dynamics they have described.
We already answered a similar question from reviewer 1 and we now consider growing and non-growing cells with di↵erent competitive ability.In short, our model assumes logistic growth and thus, even in regions where wild-type cells cannot grow, they can be continually imported from neighbouring regions, and consequently occupy space, which contributes towards the carrying capacity K. To understand the e↵ect of di↵erential competition between growing and non-growing cells suggested by the reviewer, we now considered a model where only a proportion ↵ of the nongrowing cells contributes to the carrying capacity, see Fig. A1 (Supplementary Text S4, Supplementary Fig. 4).Fig. A1 shows that di↵erential competition between growing and non-growing cells does not a↵ect our key conclusions.
5. The title should indicate that this work is exclusively numerical/theoretical.
Our intention was to o↵er a title that can be understood as a new hypothesis, and not to highlight the specific approach we used to test the hypothesis.We followed closely the titles of important works in the field, for example "On the rapidity of antibiotic resistance evolution facilitated by a concentration gradient," "Mutational Pathway Determines Whether Drug Gradients Accelerate Evolution of Drug-Resistant Cells," "Acceleration of emergence of bacterial antibiotic resistance in connected microenvironments" or "Spatiotemporal microbial evolution on antibiotic landscapes."Some of these works are mathematical and others are experimental, and they do not distinguish the specific methodology used.Therefore, we are keen to keep the title as it is.
A3. Minor points 6.Several important parameters, such as the mutation rate (which is set to be di↵erent by 3 magnitudes between increase and decrease of resistance without justification) are identical to the original work by Hermsen et al.While this is adequate to show di↵erences with the original staircase model, it would be good if they could at least comment on this choice, and in the best case also look at its interaction with motility.
We have assumed the mutation rate used by Hermsen et al., indeed to make our results on the role of bacterial motility comparable to theirs, and because the choice of mutation rate was consistent with experimental evidence [A2-A4].However, we do appreciate the concern of the reviewer and we now added a new figure where we show the e↵ect of varying the mutation rate on the evolutionary dynamics of our model.See Supplementary Fig. 3.The reviewer can also see this figure in our answer to a question posed by reviewer 1, Fig. A2.Essentially, we show that for realistic mutation rates, our conclusions hold.
7. The authors should comment on the ratio of di↵erent parameters that can lead to the di↵erent motility regimes they have described, and how these ratios relate to data from biological/experimental systems, and if they are biologically meaningful.
We show in Fig. 2 that the low motility regime appears for ⌫ < and the high motility regimes appears for ⌫ > .
Therefore, the ratio ⌫/ describes the motility regime in our model.Moreover, in our answer to Major Points question 1 of the reviewer, we show that this ratio can be interpreted as a visiting number V and we relate this number to experimental studies [A5, A6].This new metric then suggests that the motility regimes described in our work are biologically meaningful.8. Mutant fitness (calculated as r(1 N x /K) ) is actually expected growth dependent on cell density, and does in my opinion not really correspond to the usage of this term in experimental work.The fitness gradient therefore rather represents the steepness of the population distribution.Another term or a comment in the text would be needed to make this less misleading.
In our work, mutant fitness follows the classic Malthusian fitness, which is frequently used to experimentally determine bacterial fitness [A12].To clarify, Malthusian fitness is defined as the net growth rate during the exponential phase of a population.In population genetics, for example, the alternative definition of absolute fitness (the expected number of adult o↵spring per generation) is often used [A13].As Malthusian fitness is often considered a better metric for models with continuous time [A14], and our model is formulated in continuous time, we use the definition of the Malthusian fitness.We added a comment to clarify this definition to Results section A, where we define mutant fitness.
We should note that there is no density-dependence in our definition of fitness in the usual sense of dependence on the density of the phenotype whose fitness is measured.This is because the term N x in the mutant fitness r(1 N x /K) corresponds to the total number of wild-type cells rather than mutant cells.The wild-type cells contribute to the carrying capacity K and modify the environment experienced by mutant cells, and therefore, modify the net growth of mutant cells during their exponential growth phase.9.The graphical representations of the model, especially the depiction of the heat maps (usually in Fig. Xb on the left) would benefit from more clarity.x and y scales should be uniformous within a figure.The comments on Figure 5 ("increase time" and "increase motility") in combination with arrows are misleading for the understanding of the Figure.
We now changed the colour maps to avoid confusion between phenotypes (blue/yellow) and the colours of the colour map.We also made x and y scales uniformous within each Figure .Finally, we removed the arrows from the the comments in Figure 5, to make the plots clearer.

REVIEWER 3
This manuscript reports results on (versions of) a mathematical / computational model.It aims to study the e↵ects of several aspects of motility on the rate of antibiotic resistance evolution in an environment that contains a concentration gradient of an antibiotic.
The model used is based on the so-called Staircase Model published in 2012 in PNAS (Hermsen et al ), which assumes a one-dimensional environment that is subdivided into compartments with increasing drug concentrations.
In the article, the analysis of the model was restricted to a particular parameter regime: it assumed that the carrying capacity was large, the migration rate between neighboring compartments low, and mutation rates lower still.The current manuscript expands this analysis by studying the regime of fast migration/motility.Next, the authors expand the model by including stochastic switching between a fast and a slow motility state, or a density-dependent motility rate.
As I see it, the main findings are presented in the first four figures.Figure 1 demonstrates that that the behavior of the high-motility regime is qualitatively di↵erent from that of the low-motility regime: it results in smoother density profiles, which reduces the growth rate of mutants with a higher level of resistant emerging in the system and hence increases the "ecological waiting time".Figure 2 shows that, as a consequence, the adaptation rate depends nonmonotonously on the motility rate, with a positive dependence in the low motility regime but a negative dependence in the high-motility regime.Figure 3 presents results of the model in which individual cells switch stochastically between two motility rates.The main point is that this extended model behaves very similarly to the original model with some e↵ective migration rate.The same holds for the model with density-dependent motility, as presented in Figure 4.
Apart from these main points presented in full in the main text, quite a lot of additional work has been shifted to supplementary figures, where the e↵ect of a fitness cost, chemotaxis, HGT, di↵erent types of antibiotics, etc. are also examined, and to the supplementary text, which contains detailed mathematical derivations and analysis.
Generally, I think the analysis is sound and the conclusions seem warranted.The text is clearly and attractively written and the figures are carefully designed.I only have minor comments.Greulich and colleagues focus on the shape of the antibiotic gradient and various mutational pathways rather than the e↵ect of cell motility on the evolution of antibiotic resistance.However, some of their results can be directly compared with our work.For example, the authors report that "for the nonuniform drug distribution [Fig.2(b)], t (the average time to resistance) varies non-monotonically as a function of the (gradient) steepness, with a minimum at 0.01" and this result is comparable to the results in our Fig.2b if the gradient steepness is replaced by motility.
To clarify, if motility increases, bacteria explore the steepness of the gradient faster, which explains the link between our results and theirs.As we discuss the relationship between motility and gradient shape already by comparing our work with [A15, A16], we now added the work of Greulich et al to this list of studies.
2. Section D and Figure 5 of the Results section discusses the e↵ect of migration on the density profile from a dynamic systems perspective.Although the mathematics presented here and in the SI are fun (for some readers), the ultimate biological messages are very intuitive and can easily be summarized in a few sentences, as the authors in fact do in the last paragraph of this section.I would advice to move this part to the supplement.
We agree that the biological messages of Section D and Figure 5 of the Results section are relatively intuitive, but this section has an important purpose, to bring together key ideas from all variants of the staircase model.The concept of density profiles is a crucial idea that runs through the entire paper.In section A, the density profiles distinguish the low and high motility regimes and give rise to mutant fitness profiles.In section B, it is discussed how the combination of density profiles and mutant fitness profiles can be used to understand the evolutionary dynamics, namely the adaptation rate.In section C, the density profiles are used to understand the impact of phenotypic switching on the evolutionary dynamics.Therefore, section D primarily serves as a section that connects ideas from previous sections.The loss of density profiles at a critical motility is therefore a particularly important feature that arises in all variants of the model (Fig. 2, 3, 4).We have now also added a new movie where we compare the dynamics  We now changed the lines in Fig. 3b to be consistent.However, notice that one of these lines describes the surface of ⌫ 1 + ⌫ 2 = 2 , as indicated in the figure.In log-log plot, this surface has the shape of a line with a smooth corner rather than a smooth line perpendicular to ⌫ 1 = ⌫ 2 (as is usual in plots with linear scale), which explains why these lines end at ⌫ i = 2 rather than ⌫ i = as suggested by the reviewer.
3. I would prefer that the parameters used in the figures are mentioned in the figure caption.
We understand the concern of the reviewer, and indeed we considered the suggested option initially.However, we realized that this option would significantly increase the length of the caption of each figure.Our figures typically have several panels and parameters would need to be specified for each panel separately.As it can be appreciated in our section dedicated to the description of parameters, Supplementary Text S10, such descriptions are lengthy.
Therefore, we favour the format originally submitted.

A6. Comments on text and notation
1. What is the rationale behind the term "overlap region"?Overlap between what?
We use the term overlap region to denote the only spatial compartment where the regions where mutant cells can grow and regions where wild-type cells cannot grow overlap.We modified the introduction of the term in Results section A to make this idea clearer.
We note that previous work used the term sink instead of overlap region [A11, A20].The term sink comes from the ecological literature where source/sink denotes a compartment with a positive/negative net growth rate, respectively.
In particular, the overlap region is the only sink occupied by the wild-type population when motility is low (Fig. 1d, Theorem 3).However, when motility is high, the wild-type occupies all sink compartments x > R and the terminology used in previous work that focused on low motility only [A11, A20] cannot be used for the overlap region x = R + 1.Using the source/sink notation, the overlap region is the only compartment where the wild-type sink and mutant source overlap.
2. In several places (both in the main text and in the supplement), assumptions are made that are not explicitly mentioned.As mentioned above, the equation a R = 1/(ET evo + ET eco ) is only approximately true.Also, I don't think it is mentioned anywhere that the analytical results rely on the assumption that K is large (at least, I think they do) and that the mutation rates are low.(Without these assumptions, I do not expect the mean-field equations and quasi-steady-state assumptions to work.)Please check the text for implicit assumption.
We now explicitly mentioned all assumptions in the main text.We note that as a reply to one of the questions of reviewer 1, we now considered the e↵ect of di↵erent mutation rates, and showed that higher mutations rates do not a↵ect the key conclusions of our work (Fig A2 here, new Supplementary Fig. 3).However, this comment made us realize that we had implicitly assumed that the first mutant in the overlap region always survives.Therefore, we now consider the mutant survival probability in our analytics (new Supplementary Text S8) and explicitly address this assumption.While the incorporation of mutant survival probability did not a↵ect our conclusions, it improved our analytical results when the environment is sink-like and motility is close to the critical motility (Fig. 2b).We thank the reviewer for this helpful comment that allowed us to make this improvement.
3. The notation of the expected value in the equation a R = 1/(ET evo + ET eco ) in the main text is not introduced.
We now explicitly mentioned the notation for expected values.Furthermore, we added upper indices on T evo and T eco to clarify the resistance state R, as explained in Comments on content, question 3. 4. On line 195, the authors write "low motility (⌫ < ) accelerates the adaptation rate".Similar phrases appear later in the manuscript, including on line 216 ("high motility (. . . ) decelerates the adaptation rate").These phrases are confusing because in fact the adaptation rate is much lower in the low-motility regime than in the high-motility regime.I would rephrase this.
The idea that heterogeneous environments can accelerate the evolution of antibiotic resistance has been introduced and discussed in the literature, for example in [A5, A16, A21, A22], and we use the same terminology when discussing the net e↵ect of bacterial motility on bacterial adaptive evolution.We show that bacterial motility can increase the adaptation rate in heterogenous environments at lower rates (i.e., accelerates bacterial evolution), consistent with the ideas published in the literature.However, we wanted to highlight that for high rates of motility, motility can actually decrease bacterial adaptation rate (i.e., decelerates bacterial evolution).It should be noted that we compare local trends, within identical heterogeneous environments that di↵er infinitesimally in cell motility.We now highlight that our comparison is local to avoid confusion whenever it is not clear if we refer to global or local trends.
6. line 328: What is meant by "stochastic switching levels out the numbers of both phenotypes"?

Figure 1d :
Figure 1d: Nature of the red bars is unclear.Therefore suggest to add description to caption.

Figure 1b :
Figure 1b: Blue and red are only explained at panel b, but are already described at panel a.Therefore suggest to move up to panel a.

Figure 1d :
Figure 1d: Nature of the red bars is unclear.Therefore suggest to add description to caption.Added.

Figure A2 |
FigureA2| E↵ect of mutation rate.The forward mutation rate µ f (towards increasing resistance) is varied in the original model.The heatmap of adaptation rate on the (⌫, µ f ) plane explains how the mutation rate µ f changes the relationship between motility rate ⌫ and adaptation rate aR.Mutation rate below µ f < 10 4 does not a↵ect the relationship between adaptation rate and motility: adaptation rate increases at low motility, decreases at high motility, and no adaptation occurs at deadly motility (grey).This mutation rate threshold corresponds to the probability of µ f / = 10 3 resistance mutations per cell division, which is very high compared to experimental estimates of 10 6 and 10 9 mutations per cell division[A2-A4].

A2. Major points 1 .
It is unclear what the relationship between ⌫ and represents in the physical world and when, in practice, should we expect transitions between high and low motility scenarios.What are the relevant units (in the real world) Figures include the original staircase model as a special case, errors are expected to be negligible in all other Figures.We now explicitly list all errors in the new Source Data file and clarify this e↵ect in Methods.For completeness, we add below the description of the computation.Detailed computation of the error bars.Let us denote the observed times between two adaptation steps in a simulation i by T R i = T R compute the standard deviation associated with estimating the adaptation rate a R by finite number n of

Figure A3 |
Figure A3| Standard deviations of the estimated adaptation rate aR.The standard deviations associated with estimating the adaptation rate aR by finite number n = 50 of simulations can be computed using the central limit theorem.We plotted the standard deviation as error bars centred at the symbols that denote the estimated mean value of the adaptation rate aR and that are used in Fig.2b.

Figure A4 |
Figure A4| Further e↵ects of density-dependent motility.(a) Explicit model.In the main text, we model motility phenotypes implicitly as all cells at a given spatial position have the same motility (implicit model).Alternatively, motility phenotypes can be modeled explicitly at each spatial position x (explicit model).In this case, motility phenotypes are represented by a new dimension in the staircase model and switch at a large rate (s ) with a bias towards low-density (resp.high-density) phenotype at low (resp.high) density (characterised by ⌧ 1), where low/high density is determined by comparing the density of cells Nx at a spatial position x to the switching threshold S. See Supplementary Text S6 for details.(c,e) Adaptation in source-like environment.Adaptation rate heatmap on the (⌫L, ⌫H ) plane for di↵erent switching thresholds S in the source-like environment, where ⌫L,H are the low-density/high-density motility rates.The (⌫L, ⌫H ) plane can be partitioned into di↵erent combinations of adaptation regimes and its diagonal corresponds to a population of a single motility.The diagonal separates slow-to-fast and fast-to-slow switching combinations, which di↵er in relative motility at low-to-high density.Bacteria generically adapt as if all cells had the same e↵ective motility, which corresponds to the intersection of the level sets (white dashed lines) with the diagonal.This e↵ective motility generically matches the low-density (resp.high-density) motility at high (resp.low) threshold S. At low threshold S, this generic rule has one exception in the implicit model.The adaptation rate is reduced in the mixed motility combination of the fast-to-slow switching case, as all cells in the overlap region switch to fast phenotype and are unlikely to start a growing mutant colony there.(b,d) Adaptation in sink-like environment.Same as panels (c,e), but for a sink-like environment.Importantly, the same e↵ective motility governs the adaptation dynamics as in panels (c,e).At low threshold S, this generic rule of e↵ective motility has three exceptions.First, the adaptation rate is reduced in the mixed motility combination of fast-to-slow switching.Second, the adaptation rate is reduced for slow-to-fast switching when the environment is sink-like and the high-density motility phenotype moves above the critical motility.Third, a deadly motility regime (grey) can only occur if the initial population is of low density and the low-density motility phenotype moves above the critical motility.In the explicit model, the critical motility is slightly increased as the explicit presence of a few slow cells with a high-density phenotype can prevent extinction.(f ) Wild-type profiles in the implicit and explicit model.The wild-type profiles are similar and di↵er only in predicted phenotypic proportions.Note that the mutant fitness (red bars) is similar in both models, which further explains why both models exhibit the same adaptation dynamics.

A4. Comments on content 1 .
Around the same time when the Staircase Model was published, Greulich et al (PRL, 2012) published a similar model.An important di↵erence between Hermsen et al and Greulich et al was that Hermsen et al focused on a series of relatively isolated compartments, whereas Greulich et al considered continuous space.It seems to me that this di↵erence disappears in the high-motility regime.While the authors cite Greulich et al, they do not compare their results to that paper in much detail.I think they should.

TABLE I |
Experimental examples of di↵erent adaptation regimes explored in our work.
Figure A5| Alternative definition of adaptation rate.(a)Adaptationrate between consecutive founder states.The adaptation process in the staircase model corresponds to jumps between resistance states R !R + 1 that happen at rate aR.The resistance states R are characterised by founder states (a first mutant appears in the overlap region) and winner states (mutants outgrow wild-type in the overlap region), which are separated in time by evolutionary T R evo and ecological T R eco waiting times.In the main text, we define the adaptation rate aR as the rate at which consecutive founder states appear, which is a definition used in previous work[A11].(b)Adaptationrate between consecutive winner states.As an alternative to adaptation rate between consecutive founder states, one can consider adaptation rate as the rate at which consecutive winner states appear.(c)Adaptationrate as a function of motility rate.The same description as in Fig.2b, but the definition of adaptation rate used is according to winner states.Both definitions of adaptation rate lead to the same key results: adaptation rate increases at low motility, decreases at high motility, and no adaptation occurs at deadly motility.Moreover, simulations (dots) of the adaptation rate match appropriate analytical computation (lines).