Slow spatial migration can help eradicate cooperative antimicrobial resistance in time-varying environments

TL;DR

It is shown that when the environment and the deme composition vary on the same timescale, the joint effect of slow migration and fluctuations can help eradicate drug resistance by speeding up and enhancing the extinction probability of resistant bacteria.

Abstract

Antimicrobial resistance (AMR) is a global threat and combating its spread is of paramount importance. AMR often results from a cooperative behaviour with shared drug protection. Microbial communities generally evolve in volatile, spatially structured settings. Migration, space, fluctuations, and environmental variability all have a significant impact on the development and proliferation of AMR. While drug resistance is enhanced by migration in static conditions, this changes in time-fluctuating spatially structured environments. Here, we consider a two-dimensional metapopulation consisting of demes in which drug-resistant and sensitive cells evolve in a time-changing environment. This contains a toxin against which protection can be shared (cooperative AMR). Cells migrate between demes and connect them. When the environment and the deme composition vary on the same timescale, strong population bottlenecks cause fluctuation-driven extinction events, countered by migration. We investigate the influence of migration and environmental variability on the AMR eco-evolutionary dynamics by asking at what migration rate fluctuations can help clear resistance and what are the near-optimal environmental conditions ensuring the quasi-certain eradication of resistance in the shortest possible time. By combining analytical and computational tools, we answer these questions by determining when the resistant strain goes extinct across the entire metapopulation. While dispersal generally promotes strain coexistence, here we show that slow-but-nonzero migration can speed up and enhance resistance clearance, and determine the near-optimal conditions for this phenomenon. We discuss the impact of our findings on laboratory-controlled experiments and outline their generalisation to lattices of any spatial dimension.

Figures (14)

• Figure 1: Microbial community model. Panel A: Eco-evolutionary dynamics in an isolated deme ($m=0$) subject to constant antimicrobial input rate and intermediate environmental switching (Model & Methods). Top: Illustrative temporal evolution when the environment switches between mild ($K_+=12$) and harsh ($K_-=6$) environments at rates $\nu_{\pm}\lesssim 1$, with cooperation threshold $N_{\text{th}}=3$ (Model & Methods). Resistant microbes (blue, $R$) produce a resistance enzyme that locally inactivates the drug (green shade) at a metabolic cost. When $N_{R}\geq N_{\text{th}}$, the drug is inactivated in the entire deme and sensitive cells (red, $S$) benefit from the protection at no cost (e.g., bottom-left green shade). The fraction of $S$ thus increases (solid red arrow). When $N_{R}< N_{\text{th}}$, the drug hampers the spread of $S$ (top-right red crosshairs) while $R$'s remain protected and thrive (blue arrow). In the mild environment (left, $K=K_+$), $N_{R}\to N_{\text{th}}$, whereas $N_S\to K_+-N_{\text{th}}$ (solid red arrow). Similarly, in the harsh environment (right, grey background, $K=K_-$), we still have $N_R\to N_{\text{th}}$ while $N_S\to K_-- N_{\text{th}}$ (blue arrow). $K$ is assumed to switch suddenly between $K_+$ and $K_-$ (environmental variability), driving the deme size ($N=N_S+N_R$) that fluctuates in time (Model & Methods, see Supplementary Fig \ref{['fig:DynEnvSwitch']} and Sec. \ref{['Sec:single-deme_PDMP']}). When $\nu_{\pm}\lesssim 1$ (intermediate switching), the deme experiences bottlenecks at every mild ($K=K_+$) to harsh ($K=K_-$) switch. When $K_+/K_-\gtrsim N_{\text{th}}$hernandez2023coupled (Model & Methods), demographic fluctuations may cause the extinction of $R$ cells after each bottleneck (curved red arrow). Bottom: Stochastic realisation of $N_S$ (red) and $N_R$ (blue) in a deme vs. time, with parameters $N_{\rm th}=40$ (dashed), $K_{+}=400$, $K_{-}=80$, $\nu_{+}=0.075$, and $\nu_{-}=0.125$ (Model & Methods). White/grey background indicates mild/harsh environment. Population bottlenecks (white-to-grey) enforce transient $N_R$ dips (blue arrows) promoting fluctuation-driven $R$ eradication (red arrow) hernandez2023coupled (Model & Methods). Panel B: Eco-evolutionary metapopulation dynamics; legend and parameters are as in A. The metapopulation is structured as a (two-dimensional) grid of connected demes, all with carrying capacity $K(t)\in \{K_-,K_+\}$ given by \ref{['eq:K(t)']}. Each $R$ and $S$ cell can migrate onto a neighbouring deme at rate $m$ (curved thin arrows, Model & Methods). Owed to local fluctuations of $N_R$, drug inactivation varies across demes (different shades). Bottlenecks can locally eradicate $R$, e.g. in deme ($*$), but migration from a neighbouring deme ($\dagger$) can rescue resistance (curved thin blue arrow). Resistance is fully eradicated when no $R$ cells survive across the entire grid (curved dashed red arrow).
• Figure 2: The eradication of $R$ cells depends on the bottleneck strength and migration rate. The shared parameters in all panels are $\nu=1$, $\delta=0.75$, $L=20$, $a=0.25$, $s=0.1$, $N_\text{th}=40$, and $K_-=80$ (Model & Methods) with migration according to \ref{['eq:Mig']}. Other parameters are as listed in Table \ref{['tab:sim_params']}. Panel A: Heatmap of the probability $P(N_R(t)= 0)$ of total extinction of $R$ (resistant) cells as a function of bottleneck strength, $K_+/K_-$, and migration rate $m$ at time $t=500$. Each $(m, K_+/K_-)$ value pair represents an ensemble average of ${\cal R}=200$ independent simulations, where we show the fraction of realisations resulting in complete extinction of $R$ (resistant) microbes after 500 microbial generations (standard error of the mean in $P(N_R(t=500)=0)$ below $4\%$; see Supplementary Sec. \ref{['subsubsec:wald_interval']}). The colour bar ranges from light to dark red, where darkest red indicates complete $R$ extinction in all 200 simulations at time $t=500$, $P(N_R(t=500)=0)=1$. The green line is the theoretical prediction of Eq. \ref{['eq:mc']} and the white dashed vertical line indicates an axis break separating $m=0$ and $m=10^{-5}$ (Model & Methods). The black and white annotated letters point to the specific $(m,K_+/K_-)$ values used in the outer panels. Panels B-H: Typical example trajectories of the fraction of demes $\rho_{S}(t)$ without $R$ cells, defined by Eq. \ref{['eq:frac_SR']} and corresponding to the fixation of $S$ in the metapopulation. (The fraction of demes without $S$ cells, $\rho_{R}(t)$, is vanishingly small and unnoticeable.) Here, $\rho_{S}(t)$ is shown as a function of time (microbial generations) for the $(m, K_+/K_-)$ value pairs indicated in Panel A (see Supplementary Sec. \ref{['Sec:Model.Subsec:Comp']}).
• Figure 3: A closer look to individual demes: Migration and intermediate environmental switches shape local eradication of $R$ cells. Example eco-evolutionary dynamics of the metapopulation in a single simulation realisation. Parameters are $K_+=2000$, $\nu=0.1$, $\delta=0.5$, and $m=0.001$, with density-dependent migration according to Eq. \ref{['eq:Mig']}; other parameters are as in Table \ref{['tab:sim_params']}. Panels A-F: Snapshots of the $20\times20$ metapopulation at six microbial generation times $t\in\{70, 95, 105, 125, 300, 499\}$. Red pixels indicate $R$-free demes (containing only $S$ cells) and pink pixels are demes where $R$ and $S$ cells coexist. The two demes, $\vec{u}$ and $\vec{v}$, whose time composition is tracked in Panels G and H are indicated by a black border. Panel A shows the metapopulation a few generations after an environmental bottleneck. From panels A to B no bottleneck occurs, and many $S$-only demes are recolonised by $R$ cells (many red pixels become pink). Between B and C, the metapopulation experiences a bottleneck causing a burst of local $R$ extinctions (with burst of randomly located red pixels, see also the spike of $\rho_S(t=105)$ in Panel I). Panel D: Pink clusters spread across the grid due to the migration of $R$ cells causing many recolonisation events ($\rho_S(t)$ in Panel I decreases for $t\in [105,125]$). Panels E-F: After a sequence of bottlenecks starting at $t\approx220$, the number of $S$-only demes increases overwhelmingly across the grid ($\rho_S(t\lessapprox 220) \rightarrow 1$ in Panel I), and resistance persists only in a few demes where $R$ and $S$ coexist. See Supplementary Sec. \ref{['sec:Movies']} Movie 3 for a video of the full spatial metapopulation dynamics for this example realisation and its detailed description. Panels G-H: Temporal evolution of the fraction of resistant cells $N_{R}(\vec{u},t)/N(\vec{u},t)$ and $N_{R}(\vec{v},t)/N(\vec{v},t)$ in the example demes $\vec{u}$ and $\vec{v}$ indicated as highlighted pixels in Panels A-F. Green bands indicate periods in the harsh environment (where $K_-=80$); harsh periods shorter than 1 microbial generation are not shown (Supplementary Sec. \ref{['subsubsec:harsh_env']}). Each transition from white background to a green band indicates an environmental bottleneck. The deme $\vec{u}$ of panel G first exhibits $R/S$ coexistence, followed by fluctuation-driven $R$ eradication at $t\simeq70$ due to environmental bottlenecks. In Panel H, similar dynamical development is followed by the restoration of resistance through recolonisation of the deme by $R$ cells, as indicated by the blue spikes at long times ($t\approx350$, Discussion). Panel I: Temporal evolution of the fraction $\rho_S(t)$ of demes without $R$ cells (red pixels within Panels A-F, see Eq. (\ref{['eq:frac_SR']})). From left to right, the dashed vertical lines indicate the corresponding snapshot times in Panels A-F. Green background areas as in Panels G-H.
• Figure 4: Near-optimal conditions for resistance clearance: Slow migration can speed up and enhance the eradication of $R$ cells. Temporal evolution of the heatmap showing the probability $P(N_R(t)= 0)$ of $R$ extinction as a function of bottleneck strength, $K_+/K_-$, and migration rate $m$ (implemented according to Eq. \ref{['eq:Mig']}) at $t=200$ (Panel A), $t=300$ (Panel B), $t=400$ (Panel C), and $t=500$ (Panel D) with environmental switching rate $\nu=0.1$ and bias $\delta=0.5$; other parameters are as in Table \ref{['tab:sim_params']}. As in Fig \ref{['fig:KvsD_nu1_delta0.75']}A, each $(m,K_+/K_-)$ value pair is an ensemble average over 200 independent metapopulation simulations and the $P(N_R(t)= 0)$ colour bar ranges from light to dark red indicating the fraction of simulations that have eradicated $R$ cells at each snapshot in time (standard error of the mean in $P(N_R(t)=0)$ is below $4\%$; see Supplementary Sec. \ref{['subsubsec:wald_interval']}). The green and dashed white lines represent the theoretical prediction of Eq. \ref{['eq:mc']} and an eye-guiding axis break, respectively (as in Fig \ref{['fig:KvsD_nu1_delta0.75']}A). The golden lines in Panels D-E show $K_+/K_-=\frac{\nu}{mK_-}$, with $P(N_R(t)= 0)\approx 1$ in the (upper) region between the golden and green lines, according to Eq. \ref{['eq:opt_m_0']}. The grey horizontal lines in Panels A-E indicate the example bottleneck strength used in Panel E. Panel E: Probability of $R$ extinction $P(N_R(t)= 0)$ as a function of migration rate $m$ at bottleneck strength $K_+/K_-=70.7$ for $t=100, 200, 300, 400, 500$ (bottom to top). Solid lines (full symbols at $m=0$) show results averaged over 200 realisations; shaded areas (error bars at $m=0$) indicate binomial confidence interval computed via the Wald interval (see Supplementary Sec. \ref{['subsubsec:wald_interval']}). Panel F: 90th and 95th percentile ($\tau_{90}$ and $\tau_{95}$ respectively) of $R$ eradication times as function of the migration rate with a bottleneck strength $K_+/K_-=400$ (see Supplementary Sec. \ref{['subsubsec:percentile']}). Panel F shows a single minimum at $m^* \approx 10^{-3.5}$ corresponding to $\tau_{90/95}(m^*)=\tau(m^*)=t^*\approx 240-270$.
• Figure S1: Microbial dynamics in an isolated deme subject to a switching environment.(A) Example evolution of the number of $S$ (red curve) and $R$ cells (blue) in an isolated deme (in microbial generation time, see Sec. \ref{['Sec:Model.Subsec:Comp']}) for a slow-switching environment. Parameters are $N_{\rm th}=45$ (blue dashed line), $K_{+}=400$, $K_{-}=80$, $\nu=0.001$, $\delta=0.25$, $s=0.1$, and $a=0.25$. Here, the slow-switching environment stays in the mild state ($K=K_{+}$), where $N_{S}\approx K_{+}-N_{\text{th}}$ while $N_{R}\approx N_{\text{th}}$ (Model & Methods and Sec. \ref{['Sec:single-deme_PDMP']}). The dynamics is thus the same as in a static environment with carrying capacity $K=K(t=0)$. (B) Example realisation as in panel (A) for an intermediate-switching environment with $\nu=0.75\sim s$. The environment switches back and forth between harsh ($K=K_{-}$, grey background shade) and mild states ($K_{+}$, white shade). Frequent environmental bottlenecks (white-to-grey) are accompanied by a sequence of transient $N_R$ dips leading to the fluctuation-driven clearance of resistance (red arrow, Model & Methods and Sec. \ref{['Sec:single-deme_PDMP']}). (C) Same as in panels (A,B) for a fast-switching environment with rate $\nu=10$. Environmental variations is so rapid that the carrying capacity self averages, with $K\approx\mathcal{K}$ and $N\approx \mathcal{K}\gg 1$ (not shown), $N_R\approx N_{{\rm th}}$ and $N_R\approx \mathcal{K}-N_{{\rm th}}$ (Sec. \ref{['Sec:single-deme_PDMP']}). (D) Bimodal quasi-stationary probability density of the total population $N=N_S+N_R$ in an isolated deme, sampled from $10^4$ realisations run for $t=100$ microbial generations, with the same parameters as in (A) (black, Model & Methods and Sec. \ref{['Sec:single-deme_PDMP']}). The solid green line shows the stationary PDMP density $\rho(N)$ given by Eq. \ref{['eqS:NPDMP']} (Sec. \ref{['Sec:single-deme_PDMP']}). The dashed vertical green line indicates the average $\langle K(t)\rangle=(K_++K_-)/2+\delta(K_+-K_-)$, which is close to the average deme size $\langle N\rangle$ (Model & Methods). (E) As in (D), with the same parameters as in panel (B). The dashed green vertical line shows the PDMP approximation of the average deme size, computed using Eq. \ref{['eqS:NPDMP']} according to $\int_{K_-}^{K_+} N\rho(N)~dN$. (F) Same as in (D) and (E), with the same parameters as in panel (C). The dashed vertical green line indicates ${\cal K}=2K_{-}K_{+}/\left[\left(1+\delta\right)K_{-}+\left(1-\delta\right)K_{+}\right]$, which is close to $\langle N\rangle$ and to its PDMP approximation (Sec. \ref{['Sec:single-deme_PDMP']}).