Virus-host interactions shape viral dispersal giving rise to distinct classes of travelling waves in spatial expansions
Michael Hunter, Nikhil Krishnan, Tongfei Liu, Wolfram Möbius, Diana Fusco
aa r X i v : . [ phy s i c s . b i o - ph ] O c t Virus-host interactions shape viral dispersal giving rise to distinct classes of travellingwaves in spatial expansions
Michael Hunter, Tongfei Liu,
1, 2
Wolfram M¨obius,
3, 4 and Diana Fusco ∗ Cavendish Laboratory, University of Cambridge, Cambridge, CB3 0HE, United Kingdom Department of Physics, University of Oxford, Oxford, OX1 3PJ, United Kingdom Living Systems Institute, University of Exeter, Exeter, EX4 4QD, UK Physics and Astronomy, College of Engineering,Mathematics and Physical Sciences, University of Exeter, Exeter, EX4 4QL, UK (Dated: October 26, 2020)Reaction-diffusion waves have long been used to describe the growth and spread of populationsundergoing a spatial range expansion. Such waves are generally classed as either pulled, where thedynamics are driven by the very tip of the front and stochastic fluctuations are high, or pushed, wherecooperation in growth or dispersal results in a bulk-driven wave in which fluctuations are suppressed.These concepts have been well studied experimentally in populations where the cooperation leadsto a density-dependent growth rate. By contrast, relatively little is known about experimentalpopulations that exhibit a density-dependent dispersal rate.Using bacteriophage T7 as a test organism, we present novel experimental measurements thatdemonstrate that the diffusion of phage T7, in a lawn of host
E. coli , is hindered by steric interactionswith host bacteria cells. The coupling between host density, phage dispersal and cell lysis caused byviral infection results in an effective density-dependent diffusion rate akin to cooperative behavior.Using a system of reaction-diffusion equations, we show that this effect can result in a transition froma pulled to pushed expansion. Moreover, we find that a second, independent density-dependent effecton phage dispersal spontaneously emerges as a result of the viral incubation period, during whichphage is trapped inside the host unable to disperse. Our results indicate both that bacteriophagecan be used as a controllable laboratory population to investigate the impact of density-dependentdispersal on evolution, and that the genetic diversity and adaptability of expanding viral populationscould be much greater than is currently assumed.
I. INTRODUCTION
Spatial range expansions are ubiquitous in nature,from the expansion of invasive plant species, through themigration of ancient human populations, to the rangeshifts of many organisms to higher altitudes and lati-tudes due to climate change [1–8]. One of the hallmarksof spatial expansions is the rapid loss of genetic diver-sity due to the enhanced fluctuations at the front [9, 10].This effect can, however, be significantly mitigated inthe presence of density-dependent growth [11, 12], suchas an Allee effect [13], or density-dependant dispersal,where individuals in highly dense patches tend to dis-perse more quickly [14]. In particular, it has recentlybeen shown theoretically that the ratio between the de-terministic velocity of the front and that of its linearizedapproximation is sufficient to classify the expansions inthree distinct types of travelling waves, nominally pulled,semi-pushed and fully-pushed, which respectively exhibitqualitatively distinct behaviors in the decay of heterozy-gosity, the stochastic wandering of the front position, andthe probability distribution of the most recent commonancestor [14, 15].Because density-dependent growth can play such a cru-cial role in the evolutionary dynamic of a population, it ∗ E-mail: [email protected] has been extensively investigated in both naturally occur-ring range expansions in animals, such as the invasions ofboth Eurasian gypsy moths and house finches in NorthAmerica [16, 17], and in laboratory microbial model sys-tems, where the expansion dynamics in populations ofthe budding yeast
Saccharomyces cerevisae transitionfrom pulled to pushed as growth becomes more coop-erative [18], with a corresponding preservation of geneticdiversity [19]. In comparison, relatively little is knownabout the population dynamic of experimental systemsthat exhibit density-dependant dispersal, even if it hasbeen documented in several natural populations [20] andthe transition to pushed waves has been theoretically pre-dicted [14].One laboratory system that has been hypothesized toundergo density-dependent dispersal is bacteriophage ex-panding in a bacterial lawn. The crowded bacterial en-vironment is thought to hinder phage diffusion becauseof steric interactions, resulting in a density-dependentdiffusion rate due to the coupling between the host andthe viral population densities [21]. Direct experimentalquantification of this density-dependent diffusion is, how-ever, very limited [22], and its consequence on the frontpopulation dynamic mostly unknown.Here, we address the open questions of (i) whetherand how the rate of phage diffusion depends on the den-sity of surrounding bacteria, (ii) under what conditionstransitions to semi-pushed and fully-pushed expansionscan occur, and (iii) what role density-dependent diffusionplays in this. We first design an experimental protocolto measure the dependence of the phage diffusion rate onthe density of surrounding bacteria. We then constructa system of reaction-diffusion equations to determine thephage front velocity, demonstrating that transitions toboth semi-pushed and fully-pushed waves can occur. Wefind that the presence and location of these transitionsare controlled by two independent effects that alter thedensity-dependent diffusion of the virus: the first is asso-ciated with the steric interactions with the surroundingbacteria, while the second spontaneously emerges fromthe viral infection dynamic, which prevents a viral parti-cle from diffusing during infection of the host.Taken together, our results identify bacteriophages asa controllable laboratory model system to investigatethe role of density-dependent dispersal in evolution andprovide a quantitative explanation of the physical pro-cesses that control the phage population dynamic duringa range expansion. Going beyond phages, our findingssuggest that a broad range of viruses may expand viapushed travelling waves and, consequently, may be muchmore adaptable then previously thought.
II. EXPERIMENTAL MEASUREMENTS OFDENSITY-DEPENDENT DISPERSAL INCOLIPHAGE T7
Starting with Yin and McCaskill [21], it is usually rec-ognized that bacteria can act as a barrier to phage dif-fusion, resulting in a diffusion rate that depends on thebacterial density. This dependence is indeed necessaryin phage expansion models to correctly reproduce thenon-monotonicity of front velocity observed as a functionof bacterial density [21, 23–25]. Despite its importance,the dependence of phage diffusion rate on bacterial den-sity has never been measured experimentally, possiblybecause direct observation of phage diffusion in a densebacterial environment using optical microscopy is exceed-ingly challenging.To address this problem, we designed an experimentalsetup that leverages the different interactions betweenphage and two bacterial strains, one resistant and onesusceptible to the virus. Phage droplets were inoculatedon a uniform background of resistant
E. coli of tunabledensity placed on top of an agar plate. Since phage ad-sorption is assumed to be prevented in this strain [26, 27],the resistant bacteria serve as passive barrier to phagedispersal. To probe the phage location, several dropletsof susceptible
E. coli were placed at different distancesfrom the phage droplets. The time ∆ t required by thephage to travel the distance r between a viral droplet anda close-by susceptible bacteria droplet was monitored invivo by tracking the appearance of clearings in the sus-ceptible droplets (Fig. 1a and b).By gathering statistics over many droplet-dropletpairs, we were able to first confirm that the relationshipbetween distance travelled and mean first passage time is consistent with diffusive behavior for the whole range ofbackground densities tested (Fig. 1b), and then calculatethe rate of phage diffusion D as a function of backgroundbacterial density (Fig. 1c).Building on previous efforts to account for density-dependence in plaque models [21], we fit our data usingFricke’s Law [24, 28], which describes the diffusion of asolute through a suspension of spheroids [29]: D = 1 − b bη D ; b = BB max , (1)where b indicates the fraction of bacteria B relative toa maximum value B max and η accounts for the shapeof the cells: spherical cells correspond to η =2, while E. coli cells have previously been determined to corre-spond to η =1.67 [24]. Our experimental data allow forthe first time to estimate the two fitting parameters re-quired by Frickes’s law in this context: the free diffusionrate D (i.e. the diffusion rate in the absence of sur-rounding bacteria), and the bacterial density B max atwhich diffusion is expected to be completely halted. Weestimate D = 4 . ± . µ m /s, which is in good agree-ment with the rate of 4 µ m /s previously determined byOuchterlony double immunodiffusion in 10 g/l agar ofphage P22 (similar size and shape of T7) [30, 31]; and B max = 2 . ± . µ m − , which is consistent with thetypical dimensions of an E. coli cell (assuming
E. coli cells are approximately 0.5 × µ m, we would expect a1 µ m cross section to contain between 1 and 4 closelypacked cells, depending on their orientation and defor-mation). III. MODELLING PLAQUE GROWTH:DENSITY-DEPENDENT DIFFUSION ANDADSORPTION TO INFECTED CELLS
To investigate whether the phage expansion on a bac-terial lawn occurs as a pulled or a pushed wave, andto uncover the role of host density-dependence, we com-pare the actual front velocity with the velocity c F of thecorresponding linearized system, as their ratio has beenshown to be sufficient to determine the wave class in sin-gle species range expansions [15]. To this end, we developa mathematical model that accommodates the density-dependent diffusion we have experimentally measured.We model the spatial dynamics of bacteriophageplaque growth by considering the interactions betweenthree populations: viruses (phage) V , uninfected hostbacteria B and infected host bacteria I , similar to[21, 23–25, 32–35]. The process may be summarised as V + B rate −−→ α I delay −−−→ τ βV, (2)where β is the burst size, α is the rate of adsorption, and τ is the lysis time. Δ t Background Bacteria (b)
Host BacteriaPhage r (c)(a) Background Bacterial Density [ (cid:1) m -2 ] s - ] Experiments
Experiments r =4 D Δ t +const Δ t [s] r [ (cid:0) m ] D= B/B max B/ (cid:2) B max D B FIG. 1. Experiments show how rate of phage diffusion is reduced by surrounding bacteria. (a): The basis of the experimentalset-up, consisting of a droplet of phage and host bacteria, separated by a distance r . After time ∆ t , a plaque begins to form inthe host bacterial droplet (starred region). (b): The full experimental set-up, consisting of many phage-host droplet pairs ontop of a lawn of phage resistant bacteria of variable density. The presence of chloramphenicol in the plate media ensures thatthe background bacterial density is constant over the course of the experiment (see Methods). An example plot of r against∆ t data, with linear fit, for a resistant bacteria density of 0.36 µ m − is shown. (c): The diffusion rate obtained as a functionof resistant bacteria density (i.e. for several instances shown in (b)), fit with Fricke’s Law (Eq. 1). As the model is deterministic, without loss of general-ity, we describe these populations with a set of reaction-diffusion equations in 1D, similar to those examined byJones et al. [32]: ∂B∂t = − αV B, (3a) ∂I∂t = αV B − αV t − τ B t − τ , (3b) ∂V∂t = ∂∂x (cid:18) D ∂V∂x (cid:19) − αV B − α ∗ V I + βαV t − τ B t − τ , (3c)where V , B and I indicate the concentration of the pop-ulation as a function of space and time. The subscript isused to indicate that those terms are evaluated at time t − τ . D is the density-dependent diffusion coefficient ofthe phage, which we experimentally fit to Fricke’s law inthe previous section (Eq. 1). α ∗ = α or α ∗ = 0, depend-ing on whether adsorption to previously infected hosts isallowed or prevented, respectively. We assume that thehost bacteria are motionless and that adsorption to un-infected hosts always leads to successful infection whileneglecting desorption.Our model introduces two ingredients that are biolog-ically and physically relevant, and that are expected toaffect the front dynamic. First, in contrast to previouswork [21, 23–25, 33–35], where the diffusion coefficient D only depends on the initial bacterial concentration B ( b = B B max in Eq. 1), we allow D to vary in time andspace according to the local bacterial density ( b = B + IB max in Eq. 1), resulting in faster diffusion inside the phageclearing (Fig. 2a). Secondly, we allow for the possibil-ity that phage can adsorb to previously infected cells( − α ∗ V I term in Eq. 3c), as is the case for phage T7.The presence or absence of these two effects generatesfour model variants that are summarized in Fig. 2a: Uni-form vs. Variable Diffusion model (UDM vs. VDM), andadsorption vs. non-adsorption to infected cells (+ vs. -).In line with previous studies, we cast the equations us-ing dimensionless variables. We measure concentrationsin terms of the initial bacterial density B , time in unitsof τ , and length in units of L = p D ( B ) τ (the distancediffused within a lysis time at the front). This results inthe following set of dimensionless variables: B ≡ B/B , I ≡ I/B , V ≡ V / ( β − B , t ≡ t/τ , x ≡ x/L and K ≡ ατ B . Consequently, c = c p τ /D , where c and c are the dimensionless and dimensional velocity of theexpansion front, respectively (Fig. 2a).In these units, the UDMs are characterized by a con-stant dimensionless diffusion coefficient D = Dτ /L = 1by definition, while the VDMs exhibit a dimensionlessdensity-dependent diffusion coefficient of the form: D = 1 − f ( B + I )1 + f ( B + I ) /η . f /η − f , (4)where f = B /B max is the initial fraction of bacteria.It is important to note that D from Eq. 4 is alwaysgreater than or equal to 1, and can therefore be inter- (a) (b) UDM-UDM+ VDM-VDM+ P o p u l a t i o n D e n s i t y Position 0 0.2 0.4 0.6 0.8 1 c f=B /B max c [ mm / h r ] UDM+UDM-VDM+VDM-PulledSemi-PushedFully-Pushed K max =1.8 V Adsorption to Infected CellsAllowed Prevented
Position DD V a r i a b l e D i ff u s i o n U n i f o r m D i ff u s i o n c FIG. 2. (a): A sketch of the population concentrations B , I and V as a function of location at the expansion front (the preciselocation is not important here, as the qualitative shape of the fronts remain constant during the expansion). The front ispropagating with dimensionless velocity c to the right. The dimensionless width ∆ x I , characterising the width of the infectedregion is given by the difference in position of the uninfected ( B ) and infected fronts ( B + I ). The differing diffusion andadsorption behaviours explored lead to four different model variants in this work. Variants either have a Uniform or Variablediffusion rate (UDM or VDM respectively, black line in (a)), and adsorption to previously infected cells either does (+), ordoes not (-) occur, leading to the four model variants (UDM+, UDM-, VDM+ and VDM-). (b): Dimensionless front velocity c as a function of bacteria fraction f , with shaded regions indicating the different expansion types. Error bars on the velocitiesare smaller than the symbols. Inset also shows the dimensional velocity c . Model parameters are chosen to represent typicalT7 expansions with β = 50, τ =18 mins, and αB max =0.1 min − [21, 32, 35], corresponding to K max = 1 . preted as a “boost” in diffusion that the VDMs exhibitin the bulk of the plaque in comparison to the correspond-ing UDMs. This boost mathematically describes the de-crease in steric interactions between phage and bacteriadue to the lysis of the host as the viral infection proceeds(black line in Fig. 2a).In terms of these variables, our model (Eqs. 3) be-comes: ∂B∂t = − K ( β − V B, (5a) ∂I∂t = K ( β − V B − K ( β − V t − B t − , (5b) ∂V∂t = ∂∂x (cid:18) D ∂V∂x (cid:19) − KV B − K ∗ V I + βKV t − B t − , (5c)where K = ατ B and K ∗ = α ∗ τ B .As our goal is to determine whether the travellingwaves are either pulled or pushed, we will require the so-lution to the linearised approximation of the model. Toachieve this we expand the model (Eqs. 5) about the tipof the front where ( V , B , I ) ≈ (0, 1, 0), keeping only lin-ear terms. This results in the following set of equations: ∂B∂t = − K ( β − V , (6a) ∂I∂t = K ( β − V − K ( β − V t − , (6b) ∂V∂t = ∂∂x (cid:18) D ∂V∂x (cid:19) − KV + βKV t − . (6c) From Eqs. 5 three natural parameters emerge: the di-mensionless adsorption coefficient K = ατ B , the burstsize β and the dimensionless diffusion coefficient D . Inthe UDMs, D = 1, leaving K and β as the only twoparameters of the model. By contrast, in the VDMs, D is a function of B (Eq. 4), which entangles the effectof initial bacterial density on K and D . To decoupleadsorption and diffusion, we define a set of three newindependent parameters that we will use in the follow-ing to analyse the model variants: the initial fraction ofbacteria f = B /B max , the maximum dimensionless ad-sorption coefficient K max = ατ B max ( K = f K max ), andthe burst size β . In the linearised approximation (Eqs. 6) c F is the same for all four model variants and, therefore,depends only on the dimensionless adsorption coefficient K = f K max and the burst size β - see Methods. IV. FROM PULLED, TO SEMI-PUSHED TOFULLY PUSHED
By numerically solving the PDE system in Eqs. 5, weobtain the front velocity c and compare it with the ve-locity c F of the linearised approximation of the model(Eqs. 6, see Methods for details). In addition to front ve-locity, we also determine the characteristic width of theinfection region ∆ x I (Fig. 2a), which we will discuss lateron.The transitions between different travelling waveregimes are then determined from the ratio cc F accord-ing to Ref. [15]: (i) pulled waves for cc F = 1, (ii) semi-pushed wave for 1 < cc F < √ , (iii) fully pushed wavesfor cc F ≥ √ (see Methods). We point out here that thetransition between semi-pushed and fully pushed waveshas been uncovered and investigated only for single vari-able range expansions, so far. The viral model presentedhere is more complex because of the coupling betweenthe dynamics of different populations (bacterial and vi-ral) and because of the presence of a time delay. Asa result, the demographic noise in our model may differfrom that in Ref. [15]. While this does not affect the tran-sition between pulled and pushed waves, the distinctionbetween semi-pushed and fully pushed waves might, inprinciple, be different. It would be interesting to investi-gate the nature of stochastic fluctuations in viral modelsin future work.An example of these transitions in the different modelvariants for a set of infection parameters typical of T7is shown in Fig. 2b. Under these conditions, we observethat the UDM+ exhibits a pulled wave for the full rangeof initial bacterial fraction, while the UDM-, the VDM+and the VDM- waves become increasingly more pushedas f increases. In terms of dimensional velocity, the dif-ference between the model variants is minimal (inset inFig. 2b), justifying why these effects have gone unnoticedin past theoretical work that aimed at predicting exper-imental phage front speeds. A. Wave Transitions are Very Sensitive toVirus-Host Interactions
To generalise our findings and fully characterise theorigin and nature of the transitions in front dynamic forthe different model variants, we extend our investigationto a broader range of parameter values, by varying K max (and β , see Appendix B) about the typical parametersused in Fig. 2b.Fig. 3 shows the type of expansion that occurs in eachof the models as a function of f and K max . The resultsclearly indicate that the presence or absence of density-dependent diffusion and adsorption to infected cells candramatically alter the type of travelling wave undergoneby phage, with the UDM+ being the only model result-ing in a pulled population wave for the whole range ofparameters explored. In the following, we will provide aphysical interpretation for these observations, by identi-fying two independent effects that alter phage dispersalin a density-dependent fashion. B. Decreased Steric Effects due to Host LysisPromote the Transition to Pushed Waves at HighBacterial Densities
The first effect can be appreciated by comparing thephase diagram of the UDM+ to that of the VDM+, andis a direct consequence of the variable diffusion coefficientthat our model explicitly introduces in Eq. 4.In the VDM+, transitions to semi-pushed and fullypushed waves occur at high values of f with very weak de-pendence on K max . This results from the boost in phage VDM- f = B /B max K m a x = α τ B m a x VDM+ f = B /B max K m a x = α τ B m a x UDM+ f = B /B max K m a x = α τ B m a x UDM- f = B /B max K m a x = α τ B m a x Pulled Semi-Pushed Fully-Pushed
FIG. 3. Phase diagrams showing the expansion types for thefour model variants as a function of bacterial density f andmaximum dimensionless adsorption coefficient K max - burstsize β =50 throughout. Lines in the UDMs, and data pointsin the VDMs indicate the parameter combinations for whichnumerical integration was performed, and velocities calcu-lated. These values are interpolated to estimate the transitionboundaries between different classes of travelling waves (yel-low lines). In the UDM+, we do not observe pushed transi-tions, while in the VDM+ transitions occur at approximatelyconstant bacteria fractions. In the UDM-, as K and β arethe only free parameters, transitions occur at specific valuesof K . In the VDM-, the transitions are heuristically approx-imated as linear relationships with gradient m and intercept a ( f = mK max + a ). diffusion that occurs in the bulk of the plaque as hostcells lyse and steric effects decrease. Because the boostincreases with increasing difference in bacterial densitybetween the front ( B = B/B = 1) and the back ( B = 0),higher initial bacterial density f will generate a strongerboost. Beyond a given threshold, controlled exclusivelyby f , the phage behind the propagating front will dis-perse sufficiently fast to be able to catch up with thefront and generate a semi-pushed or even a fully-pushedwave. For what follows, it is useful to name this explicitboost to diffusion D exp , which is mathematically identi-cal to the dimensionless diffusion coefficient in Eq. 4, andreaches its maximum in the bulk of the plaque where nobacteria are left (Fig. 2a and dashed red line in Fig. 4b). C. A second “Implicit” Density-DependentDiffusion Emerges from the Viral InfectionDynamics
Since the UDMs lack the explicit density-dependentdiffusion, the appearance of transitions to pushed regimes UDM+UDM-VDM+VDM- ((cid:3)(cid:4) (cid:5)(cid:6)(cid:7)i(cid:8) (cid:9)(cid:10)(cid:11) (cid:12)(cid:13)(cid:14))(cid:15)(cid:16)(cid:17) E(cid:18)(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:24) (cid:25)(cid:26)(cid:27)(cid:28)(cid:29)(cid:30) on Diffus (cid:31) !"
Steric InteractionsImpli c t on Diffus ,-./
Infection Dynamic timetime ψ
012 (cid:215) (1+ ψ ) (1+ ψ ) D exp D exp UDM+ UDM- VDM+ VDM- ψ FIG. 4. (a): An illustration of both explicit and implicit effects to phage diffusion. Due to the explicit effect, phage diffusionis hindered by steric interactions with bacterial hosts, while the implicit effect hinders phage diffusion by trapping the virusfor a period τ during which it cannot disperse. (b): Proxies for the diffusive behaviour in each of the model variants plottedas a function of position across the expansion front. The base diffusion rate D imp + (Eq. 8) in the UDM+ (black solid line)is modified either by the term 1 + ψ (blue arrow and blue dashed line) in the UDM- (blue solid line), which accounts for thenow unhindered diffusion in the region of infected cells, or by an additional term D exp (red arrow and red dashed line) in theVDM+ (red solid line) which accounts for the hindrance due to steric effects. Both modifications occur in the VDM- (magentaline). Faint grey lines indicate the different front profiles from the model variants used to calculate the diffusion rate profiles(see Methods). (c): Diffusion profiles for three representative cases are shown, all highlighted in comparison to the semi-pushedtransition boundaries for the models shown in Fig. 3: (i) high K max and low f , where the UDM- and VDM- are pushed; (ii)intermediate f and K max , where only the VDM- is pushed; (iii) low K max and high f , where the VDM+ and VDM- are pushed.Red and blue arrows highlight the shift from the UDM+ to the VDM+ and UDM- respectively at the position where the viralprofile is approximately 3/4 times its steady-state. in the UDM- may seem surprising (Fig. 3). To under-stand the origin of these transitions, it is helpful to con-sider the effects of the parameter K = ατ B that con-trols the transition. Adsorption and incubation (quan-tified by the parameter K ) are not only key for the ef-fective growth rate of the phage population, but also forthe effective dispersal rate of the phage, as they con-trol the time and the probability that phage particlesare “trapped” in a host cell, unable to disperse. As K increases, either more phage adsorb to host cells atthe front per unit of time, or they are trapped there for longer , resulting in a hampered dispersal of the phage atthe front (Fig. 4a).To quantify this reduced dispersal, we consider a sys-tem of point-like phage particles diffusing across a field ofcompletely permeable “sticky” obstacles, mimicking hostbacteria that trap phage for a time τ . A simple mean-field analytical argument (see Methods), supported bytwo-dimensional Monte Carlo simulations (Appendix A),demonstrates that the particles in this system exhibit ahindered diffusion D compared to their free diffusion D ,such that DD = D imp = 11 + bK max , (7)where b is the local bacterial density (similarly to Eq. 1). Because of the dependence on the local bacterial density,the phage behind the front will diffuse faster relative tothe phage at the front, and might be able to catch upand contribute to the advancement of the expansion.When adsorption to infected cells occurs (UDM+),infected cells trap phage as much as uninfected cells( b = ( B + I ) /B max in Eq. 7) resulting in a dimensionlessdiffusion coefficient of the form D imp + = 1 + K B + I ) K . (8)By contrast, when adsorption to infected cells is pre-vented (UDM-), phage can no longer become trappedin the infected region behind the front ( b = B/B max inEq. 7), so that D imp − = 1 + K BK . (9)Comparison between the two expressions (black vs.blue lines in Fig. 4b, Methods) shows that preventingadsorption to infected cells is equivalent to a boost inimplicit diffusion in the infected region just behind thefront, which can be approximated to D imp − D imp + = 1 + ψ ≈ IK BK , (10)(blue dashed line in Fig. 4b). This boost is sufficient toshift the fast diffusing phage closer to the expanding frontand, if K is sufficiently large, to generate a transition topushed waves.It is important to point out that this implicit density-dependent diffusion emerges spontaneously from the vi-ral infection dynamics (common to most viruses), whereinfecting viruses trapped in the host cannot contributeto the advancement of the front until they are releasedfrom the host. As a consequence, unlike the explicitdensity-dependent diffusion, this effect cannot be eas-ily accommodated into the diffusion coefficient of ourmodel, as it does not act independently of the infectionand growth processes. Indeed, an alternative interpreta-tion for this mechanism can be provided by a density-dependent death rate. In our model, adsorption to pre-viously infected hosts is equivalent to phage ‘death’, as itresults in the permanent loss of these phage. Going fromthe case where adsorption to infected cells occurs to thecase where it does not (from + to - models) will then leadto an increase in net growth rate in the region of infectedcells, which lies at intermediate viral density (Fig. 2a).The result is a higher net growth rate at intermediatepopulation densities similar to what an Allee effect wouldgenerate in a mono-species expansion [12, 13]. D. Implicit and Explicit Density-DependentDiffusions Act Independently with MultiplicativeEffects
Because the implicit and explicit boosts to diffusiondiscussed above have different physical origins and arecontrolled by different parameters ( K and f , respec-tively), they play significant roles in different regions ofparameter space. The implicit boost that results from alack of adsorption to infected cells, encoded in (1 + ψ ), isstronger at large K max , where more phages are trappedby hosts for a longer period of time. Instead, the explicitboost caused by steric interactions, encoded in D exp , isdominant at low K max . The ratio of the two effects overparameter space is shown in Appendix C, Fig. 11.Extending the analytical argument with which we de-fined the implicit boost to diffusion, we can show that, toa first approximation, explicit and implicit effects act in-dependently over a basal diffusion coefficient (see Meth-ods). As a consequence, preventing adsorption to in-fected cells corresponds to multiplying the diffusion co-efficient by 1 + ψ (from + to - models, blue arrows inFig. 4). Similarly, including steric effects corresponds tomultiplying the diffusion coefficient by D exp (from UDto VD models, red arrows in Fig. 4). As a result, we canwrite the dimensionless diffusion coefficient of the VDM-,which exhibits both effects, as D V DM − = (1 + ψ ) D exp D imp + , (11)where all the terms are calculated with respect to B and I from the VDM- simulations. D ( V = V m a x ) D VDM- D VDM+ D UDM- D UDM+ (1+ ψ ) (1+ ψ ) D exp : D exp FIG. 5. Dimensionless diffusion rates in each of the modelsdetermined at the front position where the phage populationis 3/4 times the steady state population V max , plotted as afunction of both f and K max . Contour lines indicate levelsof constant D . The behaviour of the contours qualitativelymatches that of the transition boundaries in Fig. 3. In contrast to the other models, this function dependsnon-trivially on K max and f , making it challenging tofind a simple parameter combination that controls thetransitions to pushed waves. Nonetheless, we see that thediffusion coefficient determined at 3/4 times the steadystate phage population is able to qualitatively capture thebehavior of the transition lines in all models (Fig. 5) andit explains why the VDM- approaches the UDM- and theVDM+ at high and low K max , respectively, where eithereffects dominate (Fig. 3 and Fig. 4c). While the phagediffusion at a specific population density is not sufficientto recapitulate pushed behavior, which by definition de-pends on the whole wave dynamic, Fig. 5 illustrates thatregions in parameter space with similar effective diffusionwithin a model correspond to similar types of expansions.This supports the idea that the density-dependent diffu-sion, whether implicit or explicit, is the key ingredientthat leads to transitions to pushed waves. E. Reduced Burst Size Results in a More PushedExpansion
Fig. 3 illustrates the transition to pushed expansionsfor a fixed burst size β . Because burst size affects thegrowth rate of the phage and, consequently, its expan-sion velocity, it is natural to wonder whether it also hasany effect on the transitions between expansion types.Interestingly, we find that the general shape of the tran-sitions for the model variants does not change when de-creasing β by a factor of 2 (Fig. 8), in line with the ideathat density-dependent dispersal is the main driver toalterations in wave dynamics. The exact location of thetransitions, however, are slightly affected by the burstsize (Fig. 9). The similarity between the transition po-sitions as a function of burst size and the correspondingvelocity c F (Fig. 10) suggest that decreasing burst sizeslows down the very front of the wave facilitating thetransition to pushed regimes. A more detailed discussionof the influence of burst size can be found in Appendix B. V. DISCUSSION
In this work, we have demonstrated experimentallythat the diffusion of phage in a bacterial lawn is hinderedby steric interactions with the host bacterial cells, re-sulting in a density-dependent diffusion coefficient whichwe have quantified. Going beyond current descriptionsof plaque growth, which have considered host density-dependence only for setting a constant diffusion rate pa-rameter, we construct a reaction-diffusion model of thephage-bacteria system that explicitly incorporates a dif-fusion rate that depends on local host density, and there-fore varies in time and space. We show that this ef-fect can lead to the transition to pushed waves at highhost densities. We also show that a second, independentdensity-dependence in diffusion emerges implicitly fromthe underlying viral dynamics, whereby phage are unableto diffuse during replication within the host. We find thatwhen adsorption to infected host cells is prevented, thiseffect can also lead to the transition to pushed waves.Together, this indicates that bacteriophage offer an ex-cellent experimental system to study the effect of density-dependent diffusion on expansion dynamic.
A. Eco-Evolutionary Consequences ofDensity-Dependent Dispersal
The transition from a pulled wave to a pushed wavehas traditionally been associated with increased co-operativity between individuals, quantified by density-dependent growth, or more recently, density-dependentdispersal [13–15]. By analogy, the density-dependence inphage diffusion can be interpreted as an emergent co-operativity, which stems from the fact that as phagework together towards cell lysis, they remove bacterialobstacles, indirectly favouring the dispersal of neighbor-ing phage. The fact that the diffusion is dynamicallychanged as phage replicate could lead to interesting eco-logical feedback. Ecological feedback on diffusion hasbeen shown in other contexts to theoretically lead to pat-tern formation, and in some cases help maintain geneticdiversity and mitigate the risk of extinctions [36]. In-deed, density-dependent dispersal has been shown to bea key ingredient in a generic route to pattern formationin bacterial populations [37].It has also been shown that pushed dynamics in range expansions have significant consequences for the evolu-tion of the population. In pulled expansions, the highsusceptibility to stochastic fluctuations results in inef-ficient selection which can lead to the accumulation ofdeleterious mutations, known as expansion load [38, 39].This has been demonstrated previously in experimen-tal microbial populations [40–42]. By contrast, in fullypushed waves, the genetic drift is reduced [12], poten-tially leading to more efficient selection and a reducedrisk of accumulating deleterious mutations.Even in the case of semi-pushed waves, where fluctu-ations are high, the growth and ancestry processes arelocalised behind the tip of the front, thereby suppress-ing the strong founder effect found in pulled expansions[14, 15]. Our results therefore suggest that expandingpopulations of viruses might be subject to more efficientselection compared to other species, e.g. microbes, andbe better equipped to adapt to new environments.
B. Relevance to Real Viral Species
We find that the transition to a pushed wave canoccur due to two separate effects: an explicit density-dependent diffusion rate, caused by steric interactionsbetween the phage and the host bacteria, which is domi-nant in crowded host environments, and an implicit hin-drance to the diffusion of the phage population at thefront caused by the viral infection dynamics. The evo-lutionary benefit of reducing genetic drift is expected tobe strongest in populations that experience both effectsand where adsorption to infected host is limited (VDM-). Some bacteriophage have mechanisms that preventadsorption to already infected cells, usually by blockingreceptor sites post-infection [43]. Bacteriophage T5 pro-duces a lipoprotein (Llp) that is expressed at the be-ginning of infection, preventing superinfection by block-ing its own receptor site (FhuA protein), and protectingnewly produced phage from inactivation by binding tofree receptors released by lysed cells [44, 45]. Similarmechanisms are also well documented in several temper-ate phage. Phage ΦV10 possesses an O-acetyltransferasethat modifies the specific ΦV10 receptor site (the O-antigen of
E. coli
O157:H7) to block adsorption [46].Similarly,
Pseudomonas aeruginosa prophage D3 mod-ifies the O-antigen of LPS on the host surface to preventadsorption of the many phage that bind to the O-antigen[47]. This is similar to other
Pseudomonas prophagewhich encode for twitching-inhibitory protein (TiP) thatmodifies the type IV pilus on the
P. aeruginosa , prevent-ing further adsorption [48, 49].Indeed, mechanisms to prevent superinfection by pre-venting adsorption to infected cells have been observed inviruses beyond bacteriophage [50]. For instance, cells re-cently infected with Vaccinia virus VacV (the live vaccineused to eradicate smallpox) express two proteins that re-pel super-infecting virions, resulting in plaques that growfour-fold faster than would be expected by replication ki-netics alone [51]. Our results show that even in absenceof explicit steric effects, pushed expansions can occur ifadsorption to infected hosts is prevented (UDM-), as isthe case for VacV, suggesting that pushed waves mightbe more widespread than previously thought among dif-ferent viral systems.As much as it may be evolutionarily advantageous forthe viral population to undergo a pushed expansion, itwould be just as advantageous for the host cells to pre-vent this from happening. From the host’s point of view,it would likely be beneficial if infected cells were able toadsorb as many virions as possible. According to ourmodel (VDM+), this would reduce the effect of the im-plicit density-dependent hindrance to diffusion, leavingonly the explicit effect to act at very high host concen-trations.
C. Life History Parameters in Spatial Settings
Our experiments have shown that the rate of phagediffusion strongly depends on the host environment and,in particular, can dramatically differ from liquid culturemeasurements, where even at high overnight densities of ∼ cells/ml, the volume fraction occupied by the cellsis ∼ . x I during one lysis time interval. As a result, the di-mensionless infection region width ∆ x I is equivalent tothe dimensionless velocity c , as predicted in all the mod-els and confirmed by the numerics (Fig. 6). By utilisingphage engineered to result in fluorescent infected cells[54], imaging a growing plaque of such phage over timein both fluorescent and bright-field channels should yieldinformation about the distributions of infected and un-infected bacteria (see Fig. 2a), and enable determination UDM+UDM-VDM-VDM+1.6 ;<= >?@
ABC
DFG
FIG. 6. Dimensionless width ∆ x I shown to be equal to c forall models across the range of f and K max investigated inFig. 3, corresponding to an expansion that travels the widthof the infected region every lysis time. We attribute the smalldiscrepancy observed at lower values to the limited conver-gence of the front profile to its steady-state because of trade-offs between precision and computational cost. of the dimensional equivalent of ∆ x I . By simultaneouslymeasuring the velocity of the expansion, the lysis timein solid media and its variation during the course of anexpansion could be estimated. D. Inferring the Fisher Velocity in ExperimentalSystems
While this work provides theoretical predictions andphysical insights regarding the transition from pulled topushed waves in phage expansions, it also shows thatphage plaques represent a well-controlled model systemto investigate these different population dynamics exper-imentally. To test our results experimentally, however,we need to determine the velocity ratio c/c F . In contrastto simulations, the velocity of the linearized system c F cannot be directly measured in experiments, but it couldbe inferred by measuring the parameters that contributeto it.One possible route could consist in collapsing thegrowth and spread of the phage into a single reaction-diffusion equation, namely: ∂n∂t = ∂∂x (cid:18) D ( n ) ∂n∂x (cid:19) + nr ( n ) , (12)where n represents the density of phage, while r ( n ) and D ( n ) are respectively the density dependent growth anddiffusion rate of phage, that would implicitly describe theprocesses of diffusion, adsorption, replication and lysis.The Fisher velocity of this model is c F = 2 p r (0) D (0),where r (0) and D (0) are the growth and diffusion rate atlow density, respectively [55]. In principle, both of these0quantities could be measured experimentally, without theneed to determine the full dependence of r ( n ) and D ( n )on the complex underlying processes.Importantly, front velocity is but one of the dynamicalproperties affected by a pushed transition. It has beenshown, for instance, that the loss of genetic diversity dur-ing the expansion behaves qualitatively different in pulledversus pushed waves [15]. By controlling the density ofhost bacteria in a lawn one could tune the strength ofthe emergent co-operativity in an experimental expand-ing phage population. Either by utilising fluorescent la-belled phage [56, 57], or by periodically sequencing theexpansion, it should be possible to track the loss of ge-netic diversity over time. This would allow for the classi-fication of expanding phage populations as either pulledor pushed, without the need to map experiments to pre-cise models of front velocity. VI. MATERIALS AND METHODSBacterial Strains
Five strains of
E. coli were involved in this work. Thefirst strain,
E. coli
BW25113 (CGSC
E. coli eWM43 [26]. This strain wasfurther transformed with plasmid pAK501 (Addgene
E. coli eMTH43. The fourthstrain used,
E. coli ∆ waaC , is resistant to phage infec-tion through deletion of the waaC gene, the product ofwhich is involved in the production of lipopolysaccharide,the recognition of which is essential for the adsorption ofphage [27]. This strain was transformed previously witha plasmid expressing mCherry to yield the final strain E.coli eWM44 [26]. The strains eMTH43 and eWM44 werethe two strains used respectively as the susceptible andresistant host in our experiments.
Bacteriophage T7
The phage used in the study is the obligately lyticbacteriophage T7. The phage was originally obtained asan aliquot from the wild-type stock of the Richardsonlab (Harvard Medical School, Boston, MA). To preparestocks of this phage, phage were added to an exponen-tially growing liquid culture of BW25113, and incubatedat 37 ◦ C until clear. The lysate was then mixed withNaCl to a final concentration of 1.4 M. This was thenspun down to remove cell debris, and the resulting su-pernatant was stored at 4 ◦ C. Sample Preparation
To measure the diffusion coefficient of phage, 96 wellplate sized omni-plates, containing 35 ml of 20 g/l agar(VWR Chemicals), with 25 g/l LB (Invitrogen) - NaClconcentration 10 g/l - and 15 µ g/ml chloramphenicol(Sigma-Aldrich) were prepared and kept at room temper-ature for 2 days. The presence of chloramphenicol in theplate prevents growth of the background strain eWM44,so to maintain its density constant over the course of theexperiment. Plates were then refrigerated if they were tobe used at a later date. Overnight liquid cultures of E.coli were grown from single colonies at 37 ◦ C in 25 g/lLB with either 100 µ g/ml ampicillin (Sigma-Aldrich) or15 µ g/ml chloramphenicol (Sigma-Aldrich) for eWM44and eMTH43 respectively.To create the background lawn of bacteria, the opticaldensity at 600 nm (OD ) of the eWM44 culture wasmeasured, and diluted with 25 g/l LB to obtain the de-sired density (calculated on the basis that OD = 0 . cells/ml). A 500 µ l droplet of this culturewas then spread with glass beads (radius 4 mm) acrossthe surface of the agar until dry. This process (a 500 µ ldroplet spread with fresh beads) was repeated a furthertwo times to achieve as uniform a distribution as pos-sible. The plate was then left for a further 10 minutesbefore proceeding to the next step.10 ml of eMTH43 overnight culture was spun down andre-suspended in fresh 25 g/l LB, so as to give an OD reading of 0.50 if diluted hundredfold. 2 µ l droplets of theconcentrated culture were then pipetted onto the lawn ofeWM44 in a grid like pattern, spaced approximately 1 cmapart, and left to dry (approximately 15 mins after thelast droplet was pipetted). Each plate contained approx-imately 60 host droplets. The plate was then incubatedat 37 ◦ C for 1 hour.10 µ l of stock T7 phage (10 pfu/ml) was diluted in 100 µ l of black food dye. 1 µ l droplets of this dilution werepipetted onto the surface of the agar, in the gaps betweenthe previously pipetted droplets of eMTH43, and left todry. A schematic of the resulting set-up can be seen inFig. 1a. Data Acquisition
The plates were imaged using a Zeiss Axio Zoom.V16stereo microscope equipped with a Zeiss PlanApo Z0.5x/0.125 FWD 114 mm objective. Images of the sam-ple were taken every 20 minutes for a period of 24 hours.During the imaging period, the sample was kept with itslid on at 37 ◦ C using an ibidi Multi-Well Plate HeatingSystem.1
Data Analysis
All of the images were analysed using Fiji (v1.52h), anopen source distribution of ImageJ focused on scientificimage analysis [59, 60].The time ∆ t necessary for a clearing to appear in thedroplets of eMTH43, and the separation r between thepoint at which a clearing forms and the nearest phagedroplet was recorded (Fig. 1). By measuring r and ∆ t over many droplet-droplet pairs, one can measure themean first-passage time t (the mean time taken for thefirst phage to diffuse a fixed distance r ), and calculate therate of phage diffusion by fitting the relationship [61]: r = 4 Dt + constant , (13)where D is the diffusion coefficient, and the constantarises due to the delay between the arrival of the phageand the formation of a visible plaque. Numerical Solutions of Reaction-Diffusion Model
We consider the scaled equations (Eq. 5) on a finite in-terval of length L D with homogeneous Neumann bound-ary conditions. Throughout we used L D = 120 and amaximum possible time of t max = 50. Initially, we set B = 1 over the whole interval, V = 1 for x ≤
2, and V = 0 elsewhere. There are initially no infected bacteria( I = 0). Solutions are determined on a mesh of uniformspace and time, with divisions of dt = 0 .
1, and dx = 0 . dx = 0 . c of the front is determined by tracking the midpoint ofthe bacterial wave (i.e. B = 0 .
5) over time. For pulledfronts, the spreading velocity is known to demonstrate apower law convergence to an asymptotic value [62]. In thecase where a steady spreading velocity was not reached,the spreading velocity was given by the asymptotic valuewhich produced the best power law fit to the data.The transitions between pulled, semi-pushed and fullypushed have been found to occur at specific wave veloc-ities with respect to the linearised (Fisher) velocity c F -the velocity determined solely by the linear dynamics atthe tip of the front [15]. Pulled expansions spread witha velocity equal to the velocity of the linearised model c = c F , while pushed expansions spread faster [11]. Thetransition between semi-pushed and fully pushed occursat a velocity of cc F = √ , with waves below this ve-locity being semi-pushed, and above this velocity beingfully pushed [15]. These thresholds have been shownto be robust to the details of the population dynam-ics [14], though their appropriateness for multi-speciesexpansions requires further investigations. We here usethe same values for illustration purposes.Even though strictly speaking pulled expansions occurwhen the spreading velocity is equal to the velocity of the linearised model ( c = c F ), due to errors in determiningthe velocity over a finite time (i.e. errors due to thepower law fit when determining the asymptotic velocity),we conservatively consider velocities within 1% of thelinearised velocity as corresponding to pulled expansions.A length scale characterising the width of the infectedregion ∆ x I was also computed by tracking the separationbetween the midpoint on the wave of uninfected bacteria( B = 0 . B + I =0.5) over time. An average was takenover the final 20 time points for the reported value of∆ x I . Linearised Solution of Reaction-Diffusion Model
To determine the transition between pulled, semi-pushed and fully pushed regimes, the solution to the lin-earised model is required. To achieve this, we first lookfor travelling wave solutions to Eq. 5 in the co-movingcoordinate z ≡ x − ct where c is the dimensionless frontvelocity.As the components approach their limiting concen-trations at the leading edge of the front, the linearisedform of the model becomes valid, and so, following pre-vious work [21, 24], we assume the concentrations takethe form V = a exp( − λz ), B = 1 − a exp( − λz ) and I = a exp( − λz ) where λ is a dimensionless width pa-rameter, and a , a and a are positive constants.Substituting into the linearised version of the model(Eqs. 6) and writing in matrix notation yields: K ( β − − λc K ( β − − e − λc ) 0 − λcλ − λc + K ( β e − λc −
1) 0 0 a a a = 0 . (14)To find non-trivial solutions the determinant of the ma-trix is set to zero, leading to the characteristic equation: λ − λc + K ( β e − λc −
1) = 0 . (15)As we are assuming here that the front is pulled, thefront propagates with the minimum possible velocity [62]: c = min λ> [ c ( λ )] . (16)By implicitly differentiating Eq. 15 with respect to λ , andsetting dc/dλ = 0 according to Eq. 16, this leads to thesecond characteristic equation:2 λ − c − Kβc e − λc = 0 . (17)The dimensionless spreading velocity c is given as theunique solution to both Eq. 15 and Eq. 17 which wesolved numerically for each set of parameters.2 Analytical Model of “Implicit” Density Dependence
To develop a simple mean-field analytical model to de-scribe the effect of the underlying viral dynamics on thephage diffusion, we imagine phage diffusing across a lawnof “sticky” penetrable disks. These disks are used to rep-resent host bacteria cells that are able to adsorb phagefor a period equivalent to the lysis time, after which thephage desorb and continue to diffuse. The disks do notpose a hindrance to the phage through steric interactions.In this set-up, phage diffuse through a series of discretesteps, where phage move a certain distance with eachstep. In any given step, the probability that a phagewill become adsorbed to one of the host bacteria is p α φ ,where p α represents the probability of adsorbing when aphage encounters a host, and φ represents the fraction ofall space occupied by the host bacteria. This representsa Poisson point process, where events (adsorption) occurcontinuously and independently. Therefore, the numberof steps that a phage takes before becoming adsorbed toa host is exponentially distributed with mean t ads = p α φ .Consequently, over the period of time T = t ads + τ s , where τ s is the lysis time (in steps), the phage will only haveon average actually moved for t ads of that time.For long times (over many adsorption/desorptionevents), this process can be thought of as a hindereddiffusion process with relative diffusion coefficient equalto DD = D imp = t ads t ads + τ s = 11 + p α τ s φ , (18)where D is the free diffusion coefficient, and D imp is therelative density-dependent diffusion coefficient resultingfrom the hindrance posed by the underlying viral dy-namics, which we have termed the “implicit” density-dependence.This can be re-written in terms of the parameters usedin the main text as DD = D imp = 11 + p α τ s φ = 11 + AbK max , (19)where A is a scaling parameter given by: A = p α τ s φατ B max b . (20)To compare the parameters used in the two descriptions,we can consider that in our mean-field model φ = b . Wecan then use the fact that the term αB max determines therate at which phage are adsorbed when in contact withbacteria (as all space is filled with bacteria at B max ),and so, like p α τ s , the term αB max τ measures the thetotal probability that phage will be adsorbed over a lysistime, assuming that the phage are always in contact withbacteria (i.e. p α τ s = αB max τ ). This leads to a value for A of: A = p α τ s φατ B max b = p α τ s ατ B max = 1 . (21) Consequently, the implicit density dependence can bewritten equivalently in terms of either parameters as: DD = D imp = 11 + p α τ s φ = 11 + bK max . (22) Analytical Model Predicts Multiplicative Effects ofSteric Interactions and Infection Dynamic
The model introduced above can be modified to ac-count for the presence of steric effects. In the absence ofadsorption i.e., if excluded-volume interactions were theonly hindrance to diffusion, the average fraction of stepssuccessfully “jumped” by phage compared to the totalattempted would be 1 − φ , and we can define a relativediffusion coefficient D exp = 1 − φ . Although this is an ap-proximation as it does not take into account the fact thatjumps may be correlated, which is why it deviates fromthe more precise Fricke’s equation, it helps extending theanalytical model in the previous section.If we now introduce adsorption, as explained in theprevious section, the average number of steps taken byphage before adsorbing to an obstacle will be t ads = p α φ .In the presence of steric interactions, only a fraction 1 − φ of these steps will be successful, so that t succ = − φp α φ .Thus, on average, the success rate of jumping if bothadsorption and steric interactions are taken into accountwill be: DD = D exp + imp = t succ t ads + τ s = 1 − φ p α τ s φ , (23)which is equivalent to the product D exp D imp , indicatingthat when both implicit and explicit effects are present,the total behaviour can be expressed as the productof both effects individually. Because in our model thefunctional form is imposed on top of the PDE system,the same argument holds if we replace the simplified D exp = 1 − φ with the more precise Fricke’s law param-eterised by our experiments, resulting in Eq. 11. Implicit ‘Boost’ to Diffusion
This section will derive how the implicit diffusion ratein the UDM- can be thought of as a boost to the implicitdiffusion coefficient in the UDM+. Here, let the implicitslow-down to diffusion in the UDM+ and UDM- be de-noted by D imp + and D imp − , respectively. Following ourprevious derivations, these rates are given by: D imp + − = 11 + b + − K max , (24)where b + = B + IB max and b − = BB max , owing to the fact thatinfected cells are either adsorbing or non-adsorbing in thetwo models.3So as to compare to the dimensionless set of parametersused in the model (where D = 1 at the expansion front),we re-scale the implicit coefficients as: D imp + − = 1 + f K max ρ + − f K max , (25)where ρ + = B + I and ρ − = B .If we then compare the ratio of the two coefficients, itcan be seen that D imp − D imp + ≈ f K max Bf K max B + I ) f K max f K max (26)= 1 + ( B + I ) f K max Bf K max (27)= 1 +
If K max
Bf K max . (28)The approximation arises from the assumption that thebacterial curves B and I are the same in both expansions.This is supported by the observation that the density pro-files for B and I are very similar when compared acrossmodels (Fig. 4b). Therefore we can see that we can write D imp − in terms of D imp + as D imp − ≈ (1 + ψ ) D imp + ; ψ = If K max
Bf K max . (29) Diffusion Profiles
Fig. 4b shows the proxy diffusion rates of each of themodel variants as a function of position across the ex-pansion front. To generate this, the population profiles V , B and I of each of the model variants were taken atthe final time step of the numerical solution (so as to beas close to the steady state as possible), and then alignedso that the half max points of the density profiles of un-infected bacteria ( B = 0 .
5) coincide. These populationcurves were then used to determine the proxy dimension-less diffusion coefficients: D UDM + = D imp + , (30a) D UDM − = (1 + ψ ) D imp + , (30b) D V DM + = D exp D imp + , (30c) D V DM − = (1 + ψ ) D exp D imp + . (30d) ACKNOWLEDGMENTS
We acknowledge help from R. Majed in performingbacterial transformations. DF thanks Kirill Korolev and Oskar Hallatschek for helpful discussions. MH acknowl-edges studentship funding from EPSRC under grantnumber EP/R513180/1. TL acknowledges college grantfrom St Edmund Hall. This work was performed usingresources provided by the Cambridge Service for DataDriven Discovery (CSD3) operated by the Universityof Cambridge Research Computing Service, provided byDell EMC and Intel using Tier-2 funding from the Engi-neering and Physical Sciences Research Council (capitalgrant EP/P020259/1), and DiRAC funding from the Sci-ence and Technology Facilities Council.
Appendix A: Monte Carlo Simulations of Diffusion
Monte Carlo simulations of point-like phage diffusingacross a lawn of penetrable “sticky” disks, representingadsorbing host bacteria in the absence of steric effectswere performed. The simulations are carried out in a 2Dbox with periodic boundary conditions, populated withcircular obstacles of fixed diameter d obs , representing alawn of host bacteria cells. Obstacles are generated withrandom center co-ordinates within the boundary of thelawn until a fraction φ of the area of the lawn is occupied(obstacles may overlap). We use sticky penetrable obsta-cles to mimic the implicit hindrance to diffusion due toinfection dynamics. A hundred of point-like tracers (rep-resenting phage) are placed on the lawn and diffuse indiscrete steps. At each step, every free tracer proposes anew random co-ordinate within a circular region of radius r step of its current position. If there is no obstacle at thenew co-ordinate, the tracer jumps to the new position. If,however, there is an obstacle at the new co-ordinate, thetracer jumps to the new position and, with probability p α , adsorbs to that obstacle. Upon successful adsorption,the tracer remains bound to the obstacle for a number ofsteps τ s (representing the viral incubation period).The size of the lawn is 1000 d obs , and r step = 5 d obs .We use 100 phages for each simulation, and run 100 sim-ulations. The simulations are run for 3 million steps, toobserve equilibration of the diffusive behavior. After abrief initial transient period, we plot the mean squaredisplacement of the tracer particles in each simulation asa function of time and fit a linear function. The slope ofthis function is equivalent to 4 D where D is the diffusioncoefficient. The final diffusion coefficient is estimatedby calculating the average and the standard deviation ofthe diffusion coefficients over the 100 independent sim-ulations. The relative diffusion coefficient D imp is thendetermined by dividing this results with the theoreticalfree diffusion coefficient given by the simulation parame-ters.These results show that the hindrance posed by theunderlying viral dynamic is equivalent to an “implicit”density-dependent reduction in diffusion D imp (Fig. 7).We find quantitative agreement with the mean-field ana-lytical model discussed in the main text (Eq. 7 and Meth-ods), both in terms of the dependence on φ , and that the4 Analytical Model
HIJ =2 KLMNO =200 0.2 0.4 0.6 0.800.20.40.60.81
FIG. 7. Relative effective diffusion coefficient obtained fromMonte Carlo simulations of point-like phage diffusing througha field of penetrable “sticky” obstacles. These obstacles resultin a reduction in the effective diffusion coefficient of the phagein agreement with the mean-field analytical model of implicitdensity-dependence Eq. 22 (see Methods). Simulation dataoffset slightly for clarity. behaviour does not depend on p α and τ s independently,but rather depends only on the product p α τ s . This is inagreement with the fact that the behaviour in the UDMsdepends only on K = ατ B . Appendix B: The Impact of Burst Size
In Fig. 8 we present phase diagrams such as in Fig. 3,but for a burst size β =20 instead of β =50. It can be seenby comparing both Figures that the qualitative behaviourof the transitions remains unchanged: in the UDM-, tran-sitions are characterised by constant values of K s and K p ; in the VDM-, transitions are approximately straightlines characterised by gradients m s and m p , and inter-cepts a s and a p ( f = m s,p K max + a s,p ); in the VDM+,transitions are largely independent of K max , and occurat critical values f s and f p .As it seems the behaviour can be characterised in thesame manner when burst size is changed, we simplify ourexamination by focusing only on how the burst size altersthese specific characteristic parameters. Rather than at-tempting to produce the whole phase diagram for eachof the models at various β , as this is very computation-ally intensive when β is either small or large, we insteadchoose a specific K max value, and for various values of β ,we vary f at this K max value. From this, the parameters K s,p and f s,p as a function of β can be easily obtained(i.e. the parameters describing the UDM- and VDM+transitions respectively).To simplify our investigation of the behaviour in theVDM-, and limit the number of computationally inten-sive calculations required, we assume that m s,p are con-stant with burst size, and using data from one specific f = B /B max K m a x = α τ B m a x f = B /B max K m a x = α τ B m a x f = B /B max K m a x = α τ B m a x f = B /B max K m a x = α τ B m a x VDM-VDM+UDM+ UDM-Pulled Semi-Pushed Fully-Pushed
FIG. 8. Phase diagrams showing the expansion types forthe four model variants as a function of f and K max , β =20throughout. As can be seen by comparison to Fig. 3, thequalitative behaviour of the transitions in each of the modelsremains the same, and can be characterised using the sameparameters: K s,p for the UDM-; f s,p for the VDM+; m s,p and a s,p for the VDM-. As before, lines in the UDMs, anddata points in the VDMs indicate the parameter combinationsfor which numerical integration was performed, and velocitiescalculated. Transition boundaries (yellow lines) are inferredfrom the data points calculated. K max value we calculate how a s,p vary as a result. Inthe two cases where the full phase diagrams were com-puted, the transition gradients were calculated as m s = − . , − . m p = − . , − . β = 50 ,
20, indicating agreement to within 2 σ and1 σ respectively.We find that while the general shape of the transitionsfor the model variants does not depend on β (Fig. 8),the exact location of the transitions are affected (Fig. 9).The dependence of all of the transition parameters on β is similar across models. Above β ≈
40, the parametersexhibit only a weak dependence on burst size, whereaswhen β decreases below this value, the transition pa-rameters also decrease, increasing the parameter rangeof K max and f in which we observe a pushed wave.This behaviour qualitatively matches the dependenceof the spreading velocity of the linearised model c F onburst size (Fig. 10). While f and K max determine theability of phage in the bulk to catch up to the front andcontribute to the dynamics, either due to explicit or im-plicit hindrance to diffusion, β only contributes to thephage growth rate and, as a result, the velocity of thefront. At lower values of β , the spreading velocity isgreatly reduced as the limited number of phage releasedat the tip after each lysis event struggle to clear the hostcells around them, allowing the phage in the back to catchup more easily, regardless of the mechanism, and con-tribute to the expansion. As burst size is increased how-5 (a) (b) (c)20 40 60 80 100 20 40 60 80 10020 40 60 80 100 0.50.6 FIG. 9. Behavior of the critical parameters describing the location of the transitions (Fig. 3 and Fig. 8) as a function of burstsize β . (a): Critical values K s and K p in the UDM-. (b): Critical values a s and a p in the VDM-, calculated from transitionlocations at K max =2.2, assuming the gradients of the transitions m s and m p are approximately constant when varying β , tomaintain computational feasibility. (c): Critical values f s and f p in the VDM+, similarly calculated from transition locationsat K max =2.2. For β >
50 we did not observe transitions to fully pushed waves in the parameter regime of f ≤ .
95 investigated.Extending the parameter regime was computationally unfeasible. In each case, to determine the error bars, we assume a 1%error in the model velocities as before, and shift the velocities accordingly, to determine the resultant shift in the transitionparameters.
FIG. 10. Spreading velocity of the linearised model for K =1 . β (in the linearised model, K and β are the only independent parameters). ever, the opposite is true, although the velocity gains thatcome with increased β become increasingly marginal, asthe uninfected host within the vicinity of recently lysedcells become saturated with newly released phage. Appendix C: Ratio of Explicit to Implicit Effects
Here we present a plot of the ratio of the explicitboost to diffusion D exp , when transitioning from UDMto VDM, to the implicit boost to diffusion (1 + ψ ) whentransitioning from UDM+ to UDM-, or from VDM+ toVDM-, as a function of f and K max . This is obtained bydetermining the strength of each effect at the front posi-tion where the phage population is 3/4 times the steadystate population V max . It can be seen in Fig. 11 that the implicit boost is dominant at large K max , while theexplicit boost dominates at low K max . PQRSTU /(1+ ψ ) - 1 D exp FIG. 11. The ratio of the explicit boost to diffusion D exp tothe implicit boost to diffusion (1 + ψ ), as a function of f and K max . [1] R. I. Colautti and S. C. H. Barrett, Rapid adaptation toclimate facilitates range expansion of an invasive plant.,Science (New York, N.Y.) , 364 (2013).[2] L. L. Cavalli-Sforza, P. Menozzi, and A. Pi-azza, Demic expansions and human evolution,Science , 639 (1993).[3] N. A. Rosenberg, J. K. Pritchard, J. L. Weber, H. M.Cann, K. K. Kidd, L. A. Zhivotovsky, and M. W.Feldman, Genetic Structure of Human Populations,Science , 2381 LP (2002).[4] S. Ramachandran, O. Deshpande, C. C. Roseman,N. A. Rosenberg, M. W. Feldman, and L. L.Cavalli-Sforza, Support from the relationship ofgenetic and geographic distance in human popula-tions for a serial founder effect originating in Africa.,Proceedings of the National Academy of Sciences of the United States of America , 15942 (2005).[5] A. R. Templeton, Out of Africa again and again,Nature , 45 (2002).[6] I.-C. Chen, J. K. Hill, R. Ohlem¨uller, D. B. Roy,and C. D. Thomas, Rapid range shifts of speciesassociated with high levels of climate warming.,Science (New York, N.Y.) , 1024 (2011).[7] G. M. Hewitt, Some genetic consequences of iceages, and their role in divergence and speciation,Biological Journal of the Linnean Society , 247 (1996).[8] G. Hewitt, The genetic legacy of the quaternary ice ages,Nature , 907 (2000).[9] O. Hallatschek and D. R. Nelson,Gene surfing in expanding populations,Theoretical Population Biology , 158 (2008).[10] S. Klopfstein, M. Currat, and L. Excoffier, The Fate ofMutations Surfing on the Wave of a Range Expansion,Molecular Biology and Evolution , 482 (2006).[11] W. van Saarloos, Front propagation into unstable states,Physics Reports , 29 (2003).[12] L. Roques, J. Garnier, F. Hamel, andE. K. Klein, Allee effect promotes diver-sity in traveling waves of colonization.,Proceedings of the National Academy of Sciences of the United States of America , 8828 (2012).[13] W. C. Allee, Co-Operation Among Animals,American Journal of Sociology , 386 (1931).[14] G. Birzu, S. Matin, O. Hallatschek, and K. S. Ko-rolev, Genetic drift in range expansions is very sen-sitive to density dependence in dispersal and growth,Ecology Letters , 1817 (2019).[15] G. Birzu, O. Hallatschek, and K. S. Korolev, Fluc-tuations uncover a distinct class of traveling waves,Proceedings of the National Academy of Sciences of the United States of America , E3645 (2018).[16] D. M. Johnson, A. M. Liebhold, P. C. Tobin, and O. N.Bjørnstad, Allee effects and pulsed invasion by the gypsymoth, Nature , 361 (2006).[17] R. R. Veit and M. A. Lewis, Dispersal, popula-tion growth, and the allee effect: Dynamics ofthe house finch invasion of eastern North America,American Naturalist , 255 (1996).[18] S. R. Gandhi, E. A. Yurtsev, K. S. Korolev, andJ. Gore, Range expansions transition from pulledto pushed waves as growth becomes more coop-erative in an experimental microbial population,Proceedings of the National Academy of Sciences , 6922 (2016).[19] S. R. Gandhi, K. S. Korolev, and J. Gore,Cooperation mitigates diversity loss in aspatially expanding microbial population,Proceedings of the National Academy of Sciences of the United States of America , 23582 (2019).[20] E. Matthysen, Density-dependent dispersal in birds andmammals, Ecography , 403 (2005).[21] J. Yin and J. S. McCaskill, Replication of Virusesin a Growing Plaque - A Reaction-Diffusion Model,Biophysical Journal , 1540 (1992).[22] L. J. Alvarez, P. Thomen, T. Makushok,and D. Chatenay, Propagation of flu-orescent viruses in growing plaques,Biotechnology and Bioengineering , 615 (2007).[23] L. You and J. Yin, Amplification andSpread of Viruses in a Growing Plaque,Journal of Theoretical Biology , 365 (1999).[24] J. Fort and V. M´endez, Time-delayed spread ofviruses in growing plaques, Physical Review Letters ,10.1103/PhysRevLett.89.178101 (2002).[25] V. Ortega-Cejas, J. Fort, V. M´endez, and D. Campos,Approximate solution to the speed of spreading viruses,Physical Review E , 031909 (2004).[26] W. M¨obius, A. W. Murray, and D. R.Nelson, How Obstacles Perturb PopulationFronts and Alter Their Genetic Structure,PLOS Computational Biology , e1004615 (2015).[27] U. Qimron, B. Marintcheva, S. Tabor, and C. C.Richardson, Genomewide screens for Escherichiacoli genes affecting growth of T7 bacteriophage,Proceedings of the National Academy of Sciences , 19039 LP (2006).[28] H. Fricke, A mathematical treatment of the electric con-ductivity and capacity of disperse systems I. The electricconductivity of a suspension of homogeneous spheroids,Physical Review , 575 (1924).[29] J. Crank and E. P. J. Crank, The Mathematics of Diffusion , Oxford science pub-lications (Clarendon Press, 1979).[30] D. Stollar and L. Levine, Two-dimensional immunodiffu-sion, Methods in Enzymology , 848 (1963).[31] H. W. Ackermann, Classification of tailed enterobacteriaphages, Pathologie-biologie , 359 (1976).[32] D. A. Jones, H. L. Smith, H. R. Thieme, and G. R¨ost,On spread of phage infection of bacteria in a petri dish,SIAM Journal on Applied Mathematics , 670 (2012).[33] D. R. Amor and J. Fort, Virus infec-tion speeds: Theory versus experiment,Physical Review E , 061905 (2010).[34] D. R. Amor and J. Fort, Cohabitation reac-tion–diffusion model for virus focal infections,Physica A: Statistical Mechanics and its Applications , 611 (2014).[35] V. de Rioja, J. Fort, and N. Isern, Frontpropagation speeds of T7 virus mutants,Journal of Theoretical Biology , 112 (2015).[36] H. J. Park and C. S. Gokhale, Eco-logical feedback on diffusion dynamics,Royal Society Open Science , 181273 (2019).[37] M. E. Cates, D. Marenduzzo, I. Pagonabarraga, andJ. Tailleur, Arrested phase separation in reproducingbacteria creates a generic route to pattern formation,Proceedings of the National Academy of Sciences of the United States of America , 11715 (2010). [38] D. Panja, Effects of fluctuations on propagating fronts,Physics Reports , 87 (2004).[39] S. Peischl, M. Kirkpatrick, and L. Excoffier, Expansionload and the evolutionary dynamics of a species range.,The American naturalist , 81 (2015).[40] L. Bosshard, I. Dupanloup, O. Tenaillon, R. Bruggmann,M. Ackermann, S. Peischl, and L. Excoffier, Accumula-tion of Deleterious Mutations During Bacterial RangeExpansions, Genetics , 669 (2017).[41] T. Matsuyama, K. Komatsu, A. Nakahara,and M. Matsushita, Experimental Investiga-tion on the Validity of Population Dynam-ics Approach to Bacterial Colony Formation,Journal of the Physical Society of Japan , 1205 (1994).[42] A. Giometto, A. Rinaldo, F. Carrara, andF. Altermatt, Emerging predictable fea-tures of replicated biological invasion fronts,Proceedings of the National Academy of Sciences of the United States of America , 297 (2014).[43] S. J. Labrie, J. E. Samson, and S. Moineau,Bacteriophage resistance mechanisms,Nature Reviews Microbiology , 317 (2010).[44] V. Braun, H. Killmann, and C. Herrmann, Inactivationof FhuA at the cell surface of Escherichia coli K-12 by aphage T5 lipoprotein at the periplasmic face of the outermembrane, Journal of Bacteriology , 4710 (1994).[45] I. Pedruzzi, J. P. Rosenbusch, and K. P. Locher, In-activation in vitro of the Escherichia coli outer mem-brane protein FhuA by a phage T5-encoded lipoprotein,FEMS Microbiology Letters , 119 (1998).[46] L. L. Perry, P. SanMiguel, U. Minocha, A. I. Terekhov,M. L. Shroyer, L. A. Farris, N. Bright, B. L. Reuhs,and B. M. Applegate, Sequence analysis of Escherichiacoli O157:H7 bacteriophage Φv10 and identification of aphage-encoded immunity protein that modifies the O157antigen, FEMS Microbiology Letters , 182 (2009).[47] G. J. Newton, C. Daniels, L. L. Burrows,A. M. Kropinski, A. J. Clarke, and J. S. Lam,Three-component-mediated serotype conversionin Pseudomonas aeruginosa by bacteriophage D3,Molecular Microbiology , 1237 (2004).[48] J. Bondy-Denomy, J. Qian, E. R. Westra, A. Buck-ling, D. S. Guttman, A. R. Davidson, and K. L.Maxwell, Prophages mediate defense againstphage infection through diverse mechanisms,ISME Journal , 2854 (2016).[49] S. van Houte, A. Buckling, and E. R. Westra, Evolu-tionary Ecology of Prokaryotic Immune Mechanisms., Microbiology and molecular biology reviews : MMBR , 745 (2016).[50] D. Walsh and M. H. Naghavi, Exploitation of Cy-toskeletal Networks during Early Viral Infection,Trends in Microbiology , 39 (2019).[51] V. Doceul, M. Hollinshead, L. Van Der Linden, and G. L.Smith, Repulsion of superinfecting virions: A mechanismfor rapid virus spread, Science , 873 (2010).[52] J. Yin, Evolution of bacteriophage T7 in a growingplaque., Journal of Bacteriology , 1272 LP (1993).[53] M. Choua and J. A. Bonachela, Ecological andevolutionary consequences of viral plasticity,American Naturalist , 346 (2019).[54] L. Vidakovic, P. K. Singh, R. Hartmann, C. D. Nadell,and K. Drescher, Dynamic biofilm architecture confersindividual and collective mechanisms of viral protection,Nature Microbiology , 26 (2018).[55] R. A. Fisher, The Wave of Advance of AdvantageousGenes, Annals of Eugenics , 355 (1937).[56] E. J. Slootweg, H. J. H. G. Keller, M. A. Hink,J. W. Borst, J. Bakker, and A. Schots, FluorescentT7 display phages obtained by translational frameshift,Nucleic Acids Research , e137 (2006).[57] S. Huang, H. Wang, C. A. Carroll, S. J. Hayes, S. T.Weintraub, and P. Serwer, Analysis of proteins stainedby Alexa dyes, Electrophoresis , 779 (2004).[58] A. Kaczmarczyk, J. A. Vorholt, and A. Francez-Charlot, Cumate-inducible gene expression systemfor sphingomonads and other Alphaproteobacteria,Applied and Environmental Microbiology , 6795 (2013).[59] J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig,M. Longair, T. Pietzsch, S. Preibisch, C. Rueden,S. Saalfeld, B. Schmid, J.-Y. Tinevez, D. J. White,V. Hartenstein, K. Eliceiri, P. Tomancak, and A. Car-dona, Fiji: an open-source platform for biological-imageanalysis, Nature Methods , 676 (2012).[60] C. T. Rueden, J. Schindelin, M. C. Hiner, B. E. DeZonia,A. E. Walter, E. T. Arena, and K. W. Eliceiri, ImageJ2:ImageJ for the next generation of scientific image data,BMC Bioinformatics , 529 (2017).[61] G. Klein, Mean first-passage times ofBrownian motion and related problems,Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences , 431 (1952).[62] U. Ebert and W. Van Saarloos, Front propaga-tion into unstable states: Universal algebraic con-vergence towards uniformly translating pulled fronts,Physica D: Nonlinear Phenomena146