Assisted percolation of slow-spreading mutants in heterogeneous environments

Environmental heterogeneity can drive genetic heterogeneity in expanding populations; mutant strains may emerge that trade overall growth rate for an improved ability to survive in patches that are hostile to the wild type. This evolutionary dynamic is of practical importance when seeking to prevent the emergence of damaging traits. We show that a sub-critical slow-spreading mutant can attain dominance even when the density of patches is below their percolation threshold and predict this transition using geometrical arguments. This work demonstrates a phenomenon of ''assisted percolation'', where one sub-critical process assists another to achieve super-criticality.

Pesticides are used to control crop pests, antimicrobials to eliminate microbes, and cancer drugs to contain tumours.The emergence of mutants that are resistant to these agents is a major concern in all these scenarios [1][2][3].The double-edged sword of use and control on the one hand and loss of efficacy through the emergence of resistance is widely acknowledged, but we know relatively little about the role spatial structure plays in the dynamics of resistance emergence.Research on the effects of compartmentalisation, for example in the human body [4], gradients mimicking such compartmentalisation [5], and mosaic application, for example on fields of crops [6], have all emphasised the role of a reservoir with low concentrations of the control agent.Therein, mutants can originate that are then selected for in regions with a higher concentration.How these mutants then spread in a complex environment is not understood.
We here address this question from a theoretical physics perspective, for a two-dimensional environment with isolated patches that can be thought of as being protected by the control agent.We generalize the Type C variant [7] of the Eden model [8], a lattice-based model for growth which is suitable for spread in heterogeneous environments, intrinsically incorporates stochasticity, and whose computational efficiency matches the requirement to investigate large systems and many replicates (Figs.1a-b and S1, and expanding on previous work [9][10][11][12][13], see Supplementary Material).
Incorporating environmental heterogeneity, mutation, and selection in a two-species Eden model, we consider a hexagonal lattice of sites.Unoccupied sites can be either 'background' or 'patch site' varieties, whilst occupied sites are either either 'Wild-Type' (WT) or 'Mutant' (M) pathogens (Fig. 1a).During each time step an occupied site on the population frontier is selected to reproduce, at which point offspring of the same type as the parent are placed in a random unoccupied neighbouring site; only M are able to reproduce into a patch site, to represent their resistance to the control measure.During reproduction, there is a finite probability µ of mutation from WT to M (Figs.1a and S1a).Selection of WT or M sites to reproduce occurs with probability proportional to their 'fitness' which is taken to be 1 for WT and F < 1 for M, modelling the cost of resistance [14] (Fig. 1b).
In finite-width spreading fronts, M will eventually come to dominate as this is the only absorbing state.However, the timescale can vary dramatically: Even in the absence of patch sites, one can distinguish a supercritical phase of fast fixation and a sub-critical phase of exponentially slow fixation [15].Below criticality, small clusters of M appear, but typically die out before coalescing with others.Increasing M fitness or mutation rate causes these clusters to grow in size or frequency, respectively, to the point where multiple coalescence events can occur and the M population becomes supercritical (Fig. S1b).For a similar model with a flat expanding front (Fig. S2a), the dynamics fall into the directed percolation universality class [16].The roughness inherent to the Eden Model we use as the basis for our modified model greatly perturbs us from the DP universality class.Kuhr et al [15] performed phenomenological analysis on the Eden Model to determine the phase boundary at µ ≈ p * (1 − F ) 1.4 with p * ≈ 0.407 [15].For a given mutation rate, we can define the critical fitness F c (µ).As a general result, M dominates quickly if F > F c (µ).
We incorporate macroscopic environmental heterogeneity into this model by arranging treated patch sites into circular patches of fixed radii, either randomly placed (Fig. 1b) or in a lattice arrangement.This type of heterogeneity differs from that of previous work by acting asymmetrically on different genotypes [12]: They act as hard boundaries for WT [17], but are transparent to M. For random placements, continuum percolation theory gives a critical threshold of φ • c ≈ 0.68 [18] for the ratio φ of area covered by patches to total area; for φ > φ • c the patches themselves form a percolating cluster allowing only M lineages to survive.Therefore, for a given mutation rate, M dominates either if F > F c or φ > φ • c .The example shown in Fig. 1c demonstrates, however, that these are merely sufficient conditions.Rapid domination of M can occur with both fitness and patch area ratio being significantly lower than the critical thresholds.We aim to examine the full structure of the phase diagram sketched in Fig. 1d.
Close examination of simulations such as that presented in Fig. 1c and Video S1 reveals the mechanisms driving M dominating.When small M clusters intersect with a treated patch they spread through it and emerge from the other side ahead of the faster spreading WT population that is forced to take a longer route around the patch acting as an obstacle.If the 'escape region' beyond a treated patch is large enough, it will intersect with another patch and M population growth will continue.This is an effect of "assisted percolation" as we will explore later.To determine the boundary of the fast fixation phase, it is therefore necessary to compute (i) the expected size of escape regions, and (ii) the effective between-patch percolation process.
At first glance, the statistics of lone M clusters are important to this problem.These statistics are remarkably complex; to our knowledge, only the scaling behaviours have been determined for a square lattice in the literature [15].However, as demonstrated in Fig. 2a, the vast majority of isolated M clusters in the regime considered here are much smaller than the typical size of the patches (area ∼ 10 3 ).We have undertaken further analytical work in determining the dimensions of the lone M clusters for an equivalent model with a flat front (expanding on previous work [15,[19][20][21][22][23][24][25][26][27], see Supplementary Material).
How often these clusters lead to invasion of a patch depends on mutation rate and fitness.Fig. 2b illustrates that the probability to invade a patch increases linearly and then saturates with increasing mutation rate as expected; similarly, a higher fitness results in higher probability of invasion.To abstract from both the sizes of isolated clusters and their abundance, we focus on the consequences of individual patches being invaded.In this way, we capture the long-term behaviour of the 'thermodynamic limit' of a large system of patches with rare mutations.
While clusters are small relative to the patch's radius when outside a patch, they can spread unimpaired through the patch, leading to large domains within  patches, as seen in Fig. 2c.This domain may eventually become trapped within or escape the patch.Examples of the latter are shown in Fig. 3a, where the invading M domain spreads upward through the patch, and is able to escape before it can be headed off by WT.We expect the existence of the escape region and, if applicable, its height to depend on the patch invasion angle α (which acts as the starting point for a race between M and WT strains), M's fitness F (the relative speed of M), as well as stochastic effects.In fact, simulations with invasions seeded at different locations of the patch's boundary show that the median of escape region height increases with fitness F and decreases with invasion angle, i.e., it is largest if the invasion occurs at the bottom of the patch (Fig. 3c,d).
To understand this dependence quantitatively we turn to geometric arguments.Previous work has characterised front shape of a population encountering an obstacle in the absence of mutations and if front speed is the same everywhere outside the obstacle [11].There, the front shape was determined as the set of all points that can be reached within a given time.Here we aim to find the point along the symmetry axis which is reached at the same time by WT expanding around the patch (with relative speed 1) and M expanding through the patch (with relative speed F ), Figs.3d and S3a,b.Measured in units of patch radius, we find that the typical maximum extent of escape regions Λ (F, α) solves the following equation (for |α| < π 2 ): Full details of the derivation are given in the Supplementary Material.If a real, positive solution does not exist, this means that the M were cut off immediately and did not escape, thus Λ = 0.The numerical solution of Eq. 1 describes the simulation data well when varying fitness F or patch invasion angle α (Fig. 3c,d).In the following, we will limit the discussion to patches invaded at the bottom and consequently the tallest escape regions, given that we expect the majority of patch invasions to occur around α ≈ 0.
Having developed an understanding of the escape region from a single patch, we can examine the macrostructure emergent in a system of many randomly distributed patches.Three expansions for mutation rate µ = 10 −3 and with varying fitness F and patch area ratio φ are displayed in Fig. 4a.For low F and φ, patches and/or escape regions rarely overlap, while when either of these values crosses a threshold overlaps appear to lead to a growth of the fraction of M and ultimately domination of the front.To predict the patch area ratio at which this transition takes place, we estimate the percolation threshold for patches including escape regions (Figure S4): We treat the escape region (of rescaled height Λ) as a deformation to the patch shape, elongating in the direction of motion of the population front.A simple but effective heuristic is to consider M populations entering at the base of each patch, and treat the deformation as approximately elliptical.The percolation threshold for this system of ellipses can be found by rescaling the vertical direction to deform ellipses into disks.Conversely, the percolation threshold for a system of randomly distributed disks can be used to approximate the percolation threshold for the system of patches with escape regions, which is a function of fitness F and patch area ratio φ and which we denote by φ * c (F ).We reuse φ • c for the percolation threshold of disks and obtain (see Supplementary   Material for details): where Λ(F, α = 0) is the rescaled escape region height determined by Eq. 1.As expected, for vanishing escape region height, Λ = 0, we obtain φ * c = φ • c , i.e., M will only dominate if patches themselves percolate.In this argument we demonstrate that super-criticality can be achieved via the dynamics of a sub-critical M population being perturbed by the presence of a sub-critical area ratio of patches, hence the term "assisted percolation".
To test how well Eqs. 1 and 2 capture the transition from subcritical to supercritical regime, we simulated the system 50 times for a wide range of fitness values F and patch area ratios φ with patch radius R = 50 and computed the probability P D with which M fixes at the front conditional on invasion of one patch (Fig. 4, see Supplementary Material for details).To ensure that we only study the fate of a single mutation, the mutation rate is set to 0 after a mutation occurs.To ensure that we have studied the case where M has invaded a patch, as the isolated cluster grows we keep track of the number of M on the population frontier: If this value ever exceeds the diameter of a patch, we can be confident that a patch has been invaded.If a cluster collapses before this threshold is met, the simulation is re-run for the same distribution of patches.The transition region is characterised by P D being distinct from 0 and colors.φ * c (F ), indicated as a black line, indeed captures this transition region very well.This means that not only the approximations made, but also the description of macrostructures interacting with each other capture the dynamics of the system very well.
Motivated by wanting to further explore the applicability of these geometric arguments, and to develop a symmetrical patch distribution which can be designed to inhibit M domination, we considered patches organised on a hexagonal lattice.For a given patch radius R, the patch area ratio φ is a function of the lattice constant (the separation between the centres of adjacent patches).A sufficient condition for M dominating is φ to be larger than φ contact ≈ π 2 √ 3 , at which point patches are in contact and thus not leave a path for WT to propagate, with the approximation capturing lattice artifacts.Fig. 5a demonstrates that the phase transition profile permits rapid M domination below each of these thresholds, and we again computed the probability P D of M to dominate the front (Fig. 5b) following invasion of one isolated patch (motivated by [18,28,29]: See Supplementary Material for details).
To further characterise the region of fitness F and patch area ratio φ within which M dominates the front quickly, we consider two different rationales, one valid on short, the other on longer time scales.For short times, we ask whether an escape region can lead to invasion of an adjacent downstream patch.For long times, we expect M to dominate if it can propagate faster vertically through the lattice of patches than WT.This transition can be determined by comparing the path length of WT snaking around patches while M passing straight through, similar to the computation of escape region height above, which yields the analytical result: Both of these approaches are demonstrated in Fig. 5b (and fully described by Fig. S5 as well as and Supplementary Material).As Fig. 5b illustrates, the transition computed numerically for the shortterm argument and the transition based on the analytical long-term argument yielding Eq. 3 adequately match the simulation results.
Comparing Fig. 4b to Fig. 5b demonstrates that the choice of patch distribution strongly affects the phase transition.Tackling the complex optimisation problem of preventing M domination for given patch area ratio would be a natural next step in the translation of this work to an applied setting.One would also need to incorporate a finite mutation rate, which we anticipate will perturb these results.
We mapped the question of how mutants spread in a complex environment of control agents to a modified Eden model with mutations in a heterogeneous environment.In the analysis, we incorporated results from dis-parate analyses of the Eden model (mutation-selection balance in the absence of patches [15] and the perturbation of a single-strain front in the presence of obstacles [11]).Our observation that two sub-critical processes can combine to achieve super-criticality may have wider relevance: For example, one may consider the interaction between vaccine deployment and the emergence of vaccine-escape variants in epidemiology.The presence of a similar assisted percolation dynamic in social contact networks could have wide ramifications for the deployment of disease intervention strategies.
In this paper, we have mapped the question of how mutants spread in a complex environment of control agents to a modified Eden model with mutations in a heterogeneous environment.In doing so, we have demonstrated the existence of the novel dynamic of assisted percolation in a generalised version of a popular surface growth model.Although the model we have presented here is the first example of assisted percolation that we are aware of, we speculate that the phenomenon might be relevant to a range of other systems.In particular we expect that further examples may be found in the field of complex networks where, for example, a weak signal might achieve long-range transmission through a sub-critical set of amplifying nodes; potentially important applications to epidemiology and social dynamics are not hard to imagine. 1 Generalised Eden model To investigate the dynamics of mutants spreading in a heterogeneous environment, we employ a generalisation of the Eden Model [1], type C [2].The classical Eden model is a lattice model for surface growth [1,2].Sites can be either occupied or empty.During the growth process in the Eden model of type C, a site with at least one empty neighboring site is randomly chosen to occupy or propagate one neighboring site.This mimics replication on a population front and suggests a way to implement mutations during this propagation process.The original Eden model has been generalised to include occupied sites of different types, interpreted as genotypes in an expanding population [3,4].Modification of the lattice allows one to investigate the effects of environmental heterogeneity, such as inaccessible sites [5], 'disorder sites' of intermediate growth rate [6], or surface curvature [7].
In this work, we use a hexagonal lattice on which two different genotypes, denoted as wild type (WT, blue) and mutant (M, red), can propagate (Fig. S1a).M can progress onto all unoccupied sites on the lattice -the lattice is homogeneous from their perspective.WT, however, is barred from occupying hostile 'patch sites' denoted in dark grey, and the environment is heterogeneous from WT's perspective.These patch sites are arranged into circles of radius R (in units of lattice spacing), with area ratio φ.During expansion, a site with neighbors that can be occupied (denoted in Fig. S1a by black dotted sites) is randomly chosen.In this process, WT sites are chosen with relative weight 1 while M sites are chosen with weight ('fitness') F < 1.One of the neighboring sites (black arrows in Fig. S1a) is chosen to be converted into the same type as the parent site.However, in the process of converting a neighboring site, a WT site may 'mutate', i.e., convert the site into M instead of WT with probability µ.We do not consider back mutations, i.e., conversion of M sites into WT sites.The system is initialised with a linear front of WT at the bottom of the system and periodic boundary conditions are applied to the boundaries on the sides.
In this model and with this initial condition, a population of WT and M expands upwards on a lattice with regions of M forming (Fig. S1b).Due to the absence of back mutations, a front being composed of only M is the one and only absorbing state.How quickly this state is reached depends on the parameters of the system (Fig. S1b).

Random distribution of patches
Throughout this work, we consider circular patch area ratio φ as control parameter.This parameter relates to the number density n in an infinitely large system by [8,9] φ = 1 − e −nπR 2 .
We chose the number of patches to be placed as the largest integer smaller or equal to number density times domain size.Patches were placed uniformly within the domain while taking care of periodic boundary conditions.

Patches arranged in hexagonal lattice
When arranging patches in a hexagonal lattice, the system width needed to be adjusted slightly to ensure periodic boundary conditions of the system.

Emergence of mutations and conditions upon which simulations ended
Specific observables required us to slightly modify the model and define conditions upon which a given simulation ends: • Cluster Size Distribution, Fig. 2a in main text: A single M is seeded to occur on a specific site when it is occupied by WT, and µ = 0 throughout the simulation.This is to ensure that only a single M cluster is monitored.The simulation ends when the front consists entirely of WT again.
• Probability of patch invasion, Fig. 2b in main text: For a single patch and a µ > 0, an initially WT population frontier engulfs a patch.If M is placed within a patch, the repeat reports an invasion success.This processes is repeated 10000 times.
• Escape Region Height, Figs.2c,d in main text: Here, a single M is seeded to occur at a specific invasion angle, α, just outside the patch.The simulation ends when there are no M on the population frontier and is repeated a 1000 times.
• Probability of M dominating, phase diagrams in Figs. 4 and 5 in main text.For a specific M fitness, F , and patch area ratio, φ, over 50 repeats the mutation rate, µ, is finite until a mutation occurs, at which point µ is set to 0. Each repeat terminates after one of the following conditions is met: (i) When the entire population frontier consists of M (M strain has dominated system).(ii) When the entire population frontier consists of WT (M strain has gone extinct).(iii) When the height of the population frontier exceeds a limit: in the context of a procedurally-generated vertically expanding system with fixed width, this corresponds to stopping a simulation that takes too long to end.During the simulations, the number of M on the population frontier is recorded for each time step.If the simulation terminates for either of the latter two reasons, this record is checked: if the maximum number of M on the population frontier does not exceed 2R, this implies that no mutant ever invaded a patch.In this case, this run is discarded, and the simulation is repeated for the same patch distribution: these 'sub-repeats' continue until the simulation terminates with M having invaded a patch, or some limit is reached (this limit, 50, is to ensure that an unfortunate patch distribution does not result in a system where patches result in the WT being cut off before M can reasonably appear).
3 Cluster size distribution in a flat front model Constructing an analytical description of the rough front cluster size distribution is difficult.The roughness of the front in the Eden model, characterised by standard deviation of front advancement or height, is well understood in the context of the KPZ universality class [10,11].The implementation of mutation and selection perturbs the model away from the KPZ dynamics [11]; to understand the effect of the M clusters, Kuhr et al. phenomenologically extracted the dependence of cluster correlation lengths, corresponding to their vertical and horizontal extensions, on fitness and mutation rate as well as width of the domain or lattice [11].We are not aware of an analytical description of the probability density function for cluster sizes.
In order to better understand the cluster size distribution, we implemented a version of our model with a flat front and without patches (Fig. S2a).The system is updated row-by-row within the same time step, corresponding to a linear Domany-Kinzel model [12].The simplification of a flat front allows us to obtain an analytical result for the flat front cluster height distribution, which we compare to the rough-front equivalent determined from simulations (Fig. 2a in main text).This height distribution of so-called isotropic clusters, clusters in which both edges have equal probabilities of moving inwards and equal probabilities of moving outwards, has been derived previously [12,13] to be Eq.S2 below.It will be re-derived here as starting point to derive the distribution of the maximum widths of clusters.

Height Distribution
We consider mutation rate to be so low that clusters emerge from single mutations and that no other mutations occur in the context of this cluster.A cluster of interest then starts from a flat frontier of WT with a single M site of fitness F .The relative probability of a site on the next row being selected to be M depends on the two neighbors it has on the current or original row, see Fig. S2a.One of these neighbors is to be the parent of the site in question: If both are of the same type, then the site will be the same type as well (cases 1 and 3 in Fig. S2a).If they are different (case 2 in Fig. S2a), the probability of the site being occupied by M is F 1+F , where F is the relative weight that M will be the parent site in analogy to the rough-front model described above.
Let us denote the width of the M cluster on row r, i.e., the number of M sites in that row, as w r .The Markovian transition probabilities p wr,wr+1 (probability of the width of M on subsequent rows changing from w r to w r+1 ) are: p 0,0 = 1 means that once a cluster terminates it stays terminated.The initial condition is w 0 = 1: A cluster starts from a single site originating from a mutation (Fig. S2b).
Note that Eq.S1 is analogous to a random walk starting from w 0 = 1 with absorbing boundary condition at 0. The probability to step right (w

and finally to step left (w
The height H of the cluster is defined as the row r for which w r = 0 for the first time.Note that with w 0 = 1 this means there are H rows that have at least one M site.A cluster of height H involves H random walk steps (with 'stationary steps' included).
We first consider clusters that have an odd height H = 2n + 1 (n ∈ N 0 ) and in which there are no stationary steps present (we shall define the number of stationary steps in this random walk as S).The Hth step must be a left step in order for the cluster to terminate.Given that we start with w 0 = 1 and demand w 2n = 1, the remaining sequence of length 2n consists of left steps and right steps, but there can never be a point at which the total number of left steps outnumbers the number of right steps because that would mean the cluster terminated before reaching height H.This corresponds to the number of Dyck paths of length 2n, known to be the n th Catalan number C n [14]: Thus, the probability of getting to state w H = 0 from state w 0 = 1 (with w r<H > 0) without making any stationary steps is: Let us now consider the possibility of stationary steps in clusters of height H = 2n + 1.The number of left steps needs to be one larger than the number of right steps.One can only substitute pairs of left and right steps by stationary steps without changing cluster height.Thus, the number of stationary steps needs to be even.For example, a cluster of height H = 2n + 1 may consist of n left steps, n − 1 right steps, and 2 stationary steps.The probability to obtain such a cluster is the product of the probability to have two stationary steps p = b 2 , the probability to have a cluster of height H = 2n − 1 without stationary steps (see above), and a combinatorial factor for embedding the two stationary steps within the cluster ( 2n 2 ): Following the same argument we obtain the probability for a cluster of height H = 2n + 1 with 2k stationary steps): The total probability for the cluster to be of height H is then where in the last step we used b 2 /ac = 4 and a computer algebra system [15].
Expressing the results in terms of fitness F and cluster height H results in The mean of the cluster size follows to be A similar sequence of arguments for clusters of even height H = 2n yields again Eq.S2.This analytical result matches the simulations for many individual clusters for fitness F = 0.9 as shown in Fig. 2a of the main text.(Flat front simulations were implemented and performed independently from the Eden model described above.)

Width Distribution
To obtain the distribution of width W for clusters of height H, we follow the same strategy, but now consider subsets of Dyck paths.We again first restrict our derivation to clusters of odd height, i.e., H = 2n + 1 with n ∈ N 0 .
In the absence of stationary steps, a cluster of height H is described by Dyck paths of length 2n of which there are C n (C n again being the Catalan number).Of those, the number of paths that explore a width W (1 ≤ W ≤ n) from the origin is given by T (n, W ) [16,17] with , n ≥ 1 .
If the Dyck path explores a width W , then the corresponding cluster has a maximum width of W + 1.
In the absence of stationary steps (S = 0) the probability to obtain a cluster of width W and height H is thus C n P (H|S = 0) , with P (H|S = 0) derived above.To obtain P (W, H) we need to consider stationary steps and sum over all possible numbers of stationary steps.It is easy to see that Generally, the number of stationary steps needs to be even, S = 2k, as above, where the maximum value of k is n, which corresponds to a cluster of only stationary steps terminated by left step.Generalising the case of S = 2 and summation leads to Using the definitions of a, b, and c, substituting the index of summation to m with k = n − m, and reversing the sum, we obtain: Combining Eqs.S2 and S3, we obtain the distribution for the width of clusters for given height: A similar argument for clusters with even height H leads to Eq. S3 as well.
We compare this analytical expression with simulation results in Fig. S2b for F = 0.9.

Height of escape region as function of patch invasion angle
We compute the height of the escape region by equating the time it takes M to propagate through and beyond a patch with the time it takes WT to propagate around and beyond the same patch as sketched in Fig. S3.

Invasion below the patch's horizontal symmetry axis
Let us first consider invasion of a patch (radius R, centered at origin) below the patch's horizontal symmetry axis (Fig. S3a).The location of the invasion point I can be parameterized by an angle θ.To determine the maximum height of the escape region, we equate the time taken by M to propagate from I to the tip of the escape region, D, with the time WT requires to reach D. A single patch is encompassed by a linear front below the horizontal symmetry axis and the shortest path from the initial front to the patch, AB, is perpendicular to the front [5].The path then follows the perimeter of the patch, BC, before following the tangent from C to D.
Note that D, the point M and WT reach after the same amount of time, is located on the vertical axis of symmetry, because the WT path shown in blue together with a symmetric partner on the other side cut off M at point D. We assume for now that D is located above the patch, i.e., that an escape region exists.The length of the path taken by M can be obtained immediately: The path taken by WT is subdivided into three parts: The straight line AB, the circular arc BC, and the straight line CD.Point C is thereby determined by the fact that CD is tangent to the circle defining the patch.From top to bottom we obtain: • CD = OD 2 − R 2 by virtue of the right-handed triangle OCD, Without loss of generality we can set the front speed of WT to 1; the front speed of M then equals fitness F .The time for M to reach D is the time for WT to reach D is We can solve for the rescaled escape region height Λ defined as OD /R − 1 by equating T M with T W T : In the main text we use the invasion angle α = π 2 − θ instead of θ resulting in

Invasion above the patch's horizontal symmetry axis
A relation analogous to Eq. S3 can be derived for the case of invasion occurring above the patch's horizontal symmetry axis with M propagating from I to D in a straight line and WT propagating from I to D first along the perimeter of the circle to C and then following the tangent to D (Fig. S3b).We obtain: Equating both times, using rescaled escape height Λ = OD /R − 1 and invasion angle α = ρ + π 2 , yields:

Vanishing escape region height
For some values of F , θ, and ρ there may no real positive solution to Eq. S3 or Eq.S4: these correspond to cases where the M does not emerge from the patch before being cut off by WT.This corresponds to a non-existent escape region and thus Λ = 0.

Maximum escape region height
The angle of maximum escape region height can be found by implicit differentiation.We first consider Eq.S3: With ∂Λ ∂α = 0 we obtain which can be solved for Λ numerically.
No such extremum can be found for Eq.S4.Heuristically, we identify the extremum at α = 0 as local maximum of rescaled escape region height Λ (see Fig. 3d in main text).

Patches organised randomly
We expect M to quickly dominate the front, i.e., the system to be supercritical, if escape regions overlap with patches so that large clusters of these macrostructures form spanning the system.The transition from sub-criticality to super-criticality will occur at a critical patch area ratio φ * c , which itself depends on fitness F .Given the nature of the system with a front mainly propagating from bottom to top, we expect most invasions to occur at the bottom of the patch (invasion angle α = 0).In this case, macrostructures consist of the patch and an symmetric escape region whose height depends on fitness F (Fig. S4, left).Without loss of generality, we can choose R = 1, in which case escape region height is given by Λ (Eq.S5).The macrostructure thus has width 2 and height 2 + Λ.We approximate these macrostructures by ellipses of semi-minor axis 1 and semi-major axis (2 + Λ)/2 (Fig. S4, centre).Note that the ellipses are aligned.Percolation is independent of rescaling the system, for example when rescaling the height of the system such that ellipses become circles (Fig. S4, right).Figure S4: Description of the approach and approximation used to determine whether, for given fitness F and patch area ratio φ, the system displays sub-criticality (M extinction) or super-criticality (M domination).(Left) The patches are randomly distributed.With each patch invaded from below, fitness-dependent escape regions form.(Center) The macrostructure consisting of patch and escape region is approximated by an ellipse.(Right) Rescaling the height of the system allows one to recover the case of overlapping disks for which the percolation threshold is known.The system is sub-critical below and super-critical above the percolation threshold.
Let us denote the critical patch area ratio as φ * c .In the thermodynamics limit, this is related to the critical number density of patches by n * c by see, e.g., Refs.[8,9].Because the number density is independent of shape, n * c is also the critical number density of ellipses that approximate the macrostructures.Rescaling system height by 2/(2 + Λ) turns ellipses into circles.In this process, number density increases by a factor (2 + Λ)/2.The critical number density n • c of the discs representing the macrostructures after rescaling is therefore which in turn is related to the corresponding critical area ratio by Taken together, we can express critical area ratio of patches φ * c assisting percolation of the slow-spreading mutants by the critical area ratio of discs in two dimensions, φ • c ≈ 0.68 [8]: where we explicitly denote that φ * c depends on fitness F because the (rescaled) escape region height Λ depends on F .Note that Eq.S6 reduces to φ * c (F ) = φ • c for Λ(F ) = 0 as expected.When F is sufficiently small such that no escape region is expected to arise, the transition line becomes vertical as seen in Fig. 4b of the main text.
Eq. S6 represents the predicted transition in Fig. 4b of the main text which matches the simulated transition very well.This may be surprising given the approximations made: (i) invasions occurring at the bottom of patches, (ii) escape regions growing deterministically, (iii) vertical alignment, (iv) approximation by ellipses.Potentially, some of approximations may cancel each other or are negligible when averaging over many patches.
6 Patches organised in a hexagonal lattice When exploring the spread of M in an environment with regularly distributed patches, we consider the hexagonal lattice not only because it has the maximum patch area ratio φ of all the Bravais lattices [18], but it also has the most degrees of symmetry.In the hexagonal lattice, the separation between the surfaces of adjacent patches, S, is related to the patch area ratio by φ = 2πR 2 /( √ 3 (2R + S) 2 ).
Postulating that a critical area ratio φ c exists, one can easily obtain two bounds.First, for S = 0 neighboring patches are in contact and there is no path for WT to propagate vertically and therefore: Secondly, when WT can propagate vertically unhindered by obstacles, isolated M clusters will become enclosed by WT for any F < 1.The presence of these vertical paths depends upon the orientation of the lattice relative to the population front.For the case where nearest-neighbor patches are horizontally adjacent (the 'horizontal alignment', as seen in Fig. S5a, upper), vertical channels occur for S ≥ 2R.For the case where nearest neighbor patches are vertically adjacent (the 'vertical alignment', as seen in Fig. S5a, lower) vertical channels are present for S ≥ 4 √ 3 − 2 R. Thus: ≤ φ c for vertical alignment .
In the following we only consider the horizontal alignment of the lattice of patches.
We describe two approaches for more accurately estimating the critical patch area ratio φ c above which M dominates quickly: First, a 'short-term' consideration asking if invasion of a single patch leads to the invasion of an immediate downstream patch.Second, a 'long-term' consideration in which we compare the speed at which WT propagates around patches with the speed by which M propagates vertically.Note that φ c is a function of fitness F ; conversely, we can define a critical fitness φ c depending on patch area ratio φ or surface separation S.

Short-term success of M
We here ask whether invasion of a patch by M (from the base of the patch) will lead to a patch invasion in the next row of patches.To that end, we consider whether the escape region envelope overlaps with the patch in the next row as illustrated in Fig. S5b.We thereby approximate the envelope by connecting the points E (where WT and M arrive at the same time at the patch's boundary) and F (the point determining escape region height considered above).Note that this is based on the observation that the escape region has the shape of a triangle, rather than a rigorous derivation of escape region shape.
As starting condition, we consider the WT front to be approximately flat along the line AG between the patches.Previous work [5,19] predicts that the front would be perturbed due to patches below; we neglect this effect as well as stochasticity-induced roughness.
Let us first derive the location of the point E parameterized by angle β = ∠DOE.The paths for WT to reach point E from point A can be decomposed into four segments with (arc) lengths as follows: Equating the time to reach E, taking into account that M is propagating with relative speed F , we yield: To find F, located at the tip of the escape region on the axis of symmetry, we perform a similar procedure as in the case of an isolated patch (Section 4), but this time taking into account the effect of neighboring patches, i.e., the less direct path from A to D: For given fitness F and separation S, solving Eq.S8 numerically for β and subsequently Eq.S9 for OF determines both the location of E and F. We then determine numerically whether EF intersects with the neighboring patch as is the case in Fig. S5b.
In summary, making a number of assumptions outlined, we can predict whether invasion of a patch leads to invasion of a neighboring patch downstream for given fitness F and surface separation S (or equivalently patch area ratio φ), leading to critical area ratio φ c (F ) or critical fitness F c (φ).

Long-term success of M
The analysis above does not take into account the long-term time evolution of the system.If M spreads from one patch to another downstream, M does not necessarily need to invade further patches because in the next row the invasion occurs at a different invasion angle (from the bottom right instead of the bottom in Fig. S5b).This motivates the need to consider a more robust evaluation with focus on the long-term time evolution of the system, sketched in Fig. S5c, which asks whether WT following the winding paths around patches or M propagating through patches with slower relative speed F are overall faster.
The path length of WT from one row to the next around patches is circular segments ÃB and CD as well as the segment BC, which forms part of the inner tangent.The length of these paths has been considered above and we obtain: The length of the direct path from one row to the next is given by WT can outpace the M despite the longer path length if Thus, the critical fitness F LT c , depending on surface separation S using the long-term consideration, is 7 Captions for Supplementary Videos    When WT invades an empty site, it has a mutation probability µ = 0.001 of mutating into M of fitness F = 0.9.Only M are able to invade the circular patches of radius R = 40.Patches are randomly placed such that patch area ratio φ ≈ 0.7.M rapidly dominates the population frontier, indicating that the system is super-critical for this parameter pair (F, φ).Note that this is because patches are in the continuum percolation regime, creating a spanning cluster of patches that head off WT.

FIG. 1 .
FIG. 1.(a) Rules underlying the modified Eden model.The possible reproduction events of the next simulation step are represented by arrows.See also Fig. S1a.(b) Visualisation of the effects of increasing fitness F of M individuals and patch area ratio φ.See also Fig. S1b.(c) M dominates in an expanding front in the presence of treated patches (dark circles) in which only M (dark/red) survive.See Video 1 for a depiction of the full simulation.(d) Sketch of phase diagram, for sufficiently high fitness F > Fc and patch area φ > phic dominance of M is expected.Parameters for simulation in panel (c) lie outside of these regions, indicating a richer structure to the phase diagram in which assisted percolation takes place.

FIG. 2 .
FIG. 2. (a) Cluster height distributions for the Eden model (rough front) and a flat-front model, both obtained from simulations, together with the analytical result for flat fronts for F = 0.9.Inset: Typical clusters in our model.(b) Probability of an isolated patch being invaded as a function of mutation rate for different fitness values F .Black line represents a linear relationship.(c) Two examples for how clusters invade a patch with µ = 5 × 10 −4 , F = 0.9: (case I) Invasion from the bottom and (case II) invasion from the bottom left.In both cases M inside the patch lags the WT front outside.

FIG. 3 .
FIG. 3. (a) Continuing the evolution of the patch invasions in Fig. 2c with emphasis on the escape region.See Videos S2 and S3 for a depiction of the full simulation.(b) Sketch of the deterministic escape region height, found by equating the time taken for M to pass through the patch with the time taken for WT to pass around the patch to the same point.(c) Median of the normalised escape region height as a function of fitness F for patch invasion angle α = 0. Black data points: Simulation results, red line: Geometrical prediction.(d) Like panel (c) but as a function of patch invasion angle α for fitness F = 0.95.
FIG. 4. (a) Three snapshots of simulations for different values of fitness F and patch area ratio φ.See Videos S4-S6 for depictions of the full simulations.(b) Grid heat map of the probability of M dominating the front determined by simulations.The black line indicates the prediction of the boundary by Eq. 2.

FIG. 5 .
FIG. 5. (a) Three snapshots of simulations for different combinations of φ and F .See Videos S7-S9 for depictions of the full simulations.(b) Grid heat map of the probability of M dominating the front, determined by simulation.Dotted black line indicates a numerical prediction generated by consideration for short-time success, solid black line indicates analytical prediction generated by consideration for long-term success.Dashed lines indicate the mutation-fitness and patch coalescence phase transitions.

1
Figure S1: (a) (Inset) Rules underlying the modified Eden model.Each site is either a 'patch site' or part of a 'background' site.Occupied sites can be in of type 'wild type' (WT) or 'mutant' (M).Upon reproduction into an unoccupied neighboring site mutation from WT to M is possible (but not the reverse).Patch sites can only be turned into M sites.(Main panel) Section of the hexagonal lattice including the front of occupied sites.Sites that can convert a neighboring unoccupied site are marked with black dots.For sites marked 1-4, arrows indicate which neighboring unoccupied sites can be turned into occupied sites.(b) Snapshots of simulations for different sets of parameters for mutation rate µ, patch area ratio φ, and fitness F .(Top) µ = 10 −3 , φ = 0.4, F = 0.9 : M and WT coexist, clusters of M typically disappear before merging into larger and larger clusters.(Bottom, left) µ = 2 × 10 −2 , φ = 0.4, F = 0.9 : Increasing the mutation rate results in more frequent clusters coalescing, the absorbing state of all M at the front is reached quickly.(Bottom, center) µ = 10 −3 , φ = 0.4, F = 0.95 : Increasing the mutant fitness results in the clusters getting larger, nearby clusters coalesce, the absorbing state of all M at the front is reached quickly.(Bottom, right) µ = 10 −3 , φ = 0.7, F = 0.90 : Increasing the patch area ratio provides a 'geographical benefit' to M, the absorbing state of all M at the front is reached quickly.

Figure S2 :
Figure S2: (a) Flat front model where each row is updated based on the preceding row.The fate of an unoccupied site depends only upon the two occupied neighbors below it.Three individual cases are highlighted and described in the text.(b) Illustration of cluster width as a random walk with initial condition w 0 = 1 and absorbing boundary condition w = 0 (red circle).Cluster width (yellow circle) is shorthand for maximum width.(c) Probability of a cluster of height H having width W (left: H = 10, center: H = 20, right: H = 30).Black crosses indicate simulations, red lines the predictions by Eq.S3.

Figure S3 :
Figure S3: Paths of M through patches and WT around patches to the tip of the escape region.Patch invasion occurs at point I, the escape region extends up to point D, with escape region height |OD| − R. (a) Invasion occurring below the patch's horizontal symmetry axis, (b) invasion occurring above the patch's horizontal symmetry axis.

Figure S5 :
Figure S5: (a) Sketch of hexagonal lattice of patches aligned horizontally (upper) and vertically (lower).(b) Paths of WT (cyan) and M (red) taken in the short-term estimate of whether M will dominate the front, for details see Section 6.1.(c) Paths of WT (cyan) and M (red) taken in the long-term estimate of whether M will dominate the front, for details see Section 6.2.

2 (
length of circular arc, difference of ∠O O O and ∠BO O), CD = ÃB (because of symmetry), DE = βR (length of circular arc).M reaches point E via a different path, first propagating as WT to point H, mutating, and following the secant to point E. The corresponding paths lengths are: S) − R and HE = R 2 + 2 sin β.

Video 1 (
relating to Fig 1c): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 1000 × 1000 lattice points.When WT invades an empty site, it has a mutation probability µ = 0.001 of mutating into M of fitness F = 0.95.Only M are able to invade the circular patches of radius R = 50.Patches are placed such that about half of the system is covered by patches φ = 0.5.The system is expected to result in rapid M domination of the front if F ≥ 1 or φ >≈ 0.68, but here we see rapid M domination in spite of neither of these conditions being satisfied.Video 2 (relating to Fig2ci and 3ai): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 400 × 400 lattice points.When WT invades an empty site, it has a mutation probability µ = 0.0005 of mutating into M of fitness F = 0.9.Only M are able to invade the circular patch of radius R = 50.Here, a cluster of M invades the patch from the bottom, resulting in an 'escape region' forming vertically above the patch.Video 3 (relating to Fig2cii and 3aii): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 400 × 400 lattice points.When WT invades an empty site, it has a mutation probability µ = 0.0005 of mutating into M of fitness F = 0.9.Only M are able to invade the circular patch of radius R = 50.Here, a cluster of M invades the patch a significant distance away from the bottom of the patch; this results in an 'escape region' forming vertically above the patch, which is significantly smaller than the one in Fig.2ci(Video S2).

Video 4 (
relating to Fig 4ai): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 2000 × 2000 lattice points.When WT invades an empty site, it has a mutation probability µ = 0.001 of mutating into M of fitness F = 0.9.Only M are able to invade the circular patches of radius R = 40.Patches are randomly placed such that patch area ratio φ ≈ 0.4.WT persists on the population frontier for the entirety of the simulation, indicating that the system is sub-critical for this parameter pair (F, φ).

Video 5 (
relating to Fig 4aii): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 2000 × 2000 lattice points.When WT invades an empty site, it has a mutation probability µ = 0.001 of mutating into M of fitness F = 0.98.Only M are able to invade the circular patches of radius R = 40.Patches are randomly placed such that patch area ratio φ ≈ 0.4.M rapidly dominates the population frontier, indicating that the system is super-critical for this parameter pair (F, φ).

Video 6 (
relating to Fig 4aiii): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 2000 × 2000 lattice points.

Video 7 (
relating to Fig 5ai): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 2000 × 2000 lattice points.When WT invades an empty site, it has a mutation probability µ = 0.001 of mutating into M of fitness F = 0.9.Only M are able to invade the circular patches of radius R = 40.Patches are placed on a hexagonal lattice.The surface separation of the nearest-neighbor patches are chosen such that patch area ratio φ = 0.4.WT persists on the population frontier, with M clusters being finite, indicating that the system is sub-critical for this parameter pair (F, φ).

Video 8 (
relating to Fig 5aii): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 2000 × 2000 lattice points.When WT invades an empty site, it has a mutation probability µ = 0.001 of mutating into M of fitness F = 0.98.Only M are able to invade the circular patches of radius R = 40.Patches are placed on a hexagonal lattice.The surface separation of the nearest-neighbor patches are chosen such that patch area ratio φ = 0.4.We see that M rapidly dominates the population frontier, indicating that the system is super-critical for this parameter pair (F, φ) Video 9 (relating to Fig 5aiii): Time evolution of the modified Eden model, with initially flat front of WT in a domain of 2000 × 2000 lattice points.When WT invades an empty site, it has a mutation probability µ = 0.001 of mutating into M of fitness F = 0.9.Only M are able to invade the circular patches of radius R = 40.Patches are placed on a hexagonal lattice.The surface separation of the