Modeling spatial evolution of multi-drug resistance under drug environmental gradients

Abstract

Multi-drug combinations to treat bacterial populations are at the forefront of approaches for infection control and prevention of antibiotic resistance. Although the evolution of antibiotic resistance has been theoretically studied with mathematical population dynamics models, extensions to spatial dynamics remain rare in the literature, including in particular spatial evolution of multi-drug resistance. In this study, we propose a reaction-diffusion system that describes the multi-drug evolution of bacteria, based on a rescaling approach (Gjini and Wood, 2021). We show how the resistance to drugs in space, and the consequent adaptation of growth rate is governed by a Price equation with diffusion. The covariance terms in this equation integrate features of drug interactions and collateral resistances or sensitivities to the drugs. We study spatial versions of the model where the distribution of drugs is homogeneous across space, and where the drugs vary environmentally in a piecewise-constant, linear and nonlinear manner. Applying concepts from perturbation theory and reaction-diffusion equations, we propose an analytical characterization of average mutant fitness in the spatial system based on the principal eigenvalue of our linear problem. This enables an accurate translation from drug spatial gradients and mutant antibiotic susceptibility traits, to the relative advantage of each mutant across the environment. Such a mathematical understanding allows to predict the precise outcomes of selection over space, ultimately from the fundamental balance between growth and movement traits, and their diversity in a population.

Keywords: environmental heterogeneity; infectious diseases; multi-drug antibiotic resistance; selection; spatial evolutionary dynamics.

Figure 1:. Example of outcomes among two strains for constant gi over space.

A. Competitive exclusion. In this example, strains start at uniform distribution over space, with $g_1>g_2$, hence dynamics lead to a traveling wave solution for $f_1(z,t)$ and $f_2(z,t)$ with strain 1 traveling at faster speed and ultimately being the winner everywhere over long time. B. (Neutrally-stable) coexistence at 50:50 because the mutants start at equal total abundances and $g_1=g_2$. C. (Neutrally-stable) coexistence different from 50:50 because mutants start at different total abundances and $g_1=g_2$. D. (Neutrally-stable) coexistence different from 50:50 because mutants start at equal total abundances with $g_1=g_2$, but their initial distribution over space favours one of them that starts at higher abundance in the center of the domain.

Figure 2:. Coexistence example of 2 strains everywhere in space for space-dependent gi(z) which are mutually symmetric about L/2.

A. In this case, one strain is better-adapted in the first half of the domain, the other strain is better-adapted in the second-half of the domain with the selective advantages exactly counterbalanced. A. For low diffusion, the two strains coexist such that each strain dominates in frequency in the part of the domain where it experiences a relatively higher growth rate, maintaining a high-degree of spatial segregation in the system. C. As diffusion increases, the coexistence frequencies become more similar and tend towards 1/2 in both halves of the domain. D. Eventually, for very high-diffusion, the growth variation starts to matter less and less, and the two strains tend to the same frequency everywhere, leading to a uniformly homogeneous spatial distribution of diversity over space.

Figure 3:. Competitive exclusion everywhere in space, but the ultimately winning strain depends on parameters of gi(z) variation.

A. In this piece-wise growth rate example, the $g_1(z)$ and $g_2(z)$ are such that the mean growth rates for both strains are the same $\overline{g_1}={\int }_0^1{g}_1(z)dz={\int }_0^1{g}_2(z)dz=\overline{g_2}$ for $b=1/2$. Yet, even with equal spatially-averaged growth rates, the strain with the central advantage will be the winner. When $b$ changes, the final winner is a result of $b$ as well as $(\mathit{max}(g)-\mathit{min}(g))$ magnitude. B. In this example, the winner can be overturned by modulating the width of the interval where $g_1(z)>g_2(z)$, while keeping the shape of the two functions. We assume the growth rates are non-monotonic functions of space, represented by a concave and a convex parabola with vertices near the middle of the domain: $g_1(z)=m-\sigma (z-L/2{)}^2+h$ and $g_2(z)=m+2\sigma (z-L/2{)}^2$ with $m=0.3,\sigma =0.4,D=0.015$ and $h$ varied. The critical value of $h$ for overturning the final outcome is $h=0.04$. Mutant 1 loses if $h<0.04$ but it wins if $h>0.04$, when its fitness advantage in the center of the domain is sufficiently high to compensate for its disadvantage near the boundaries. This cannot be predicted with the mean growth rate difference $\overline{g_1}-\stackrel{-}{g_2}$ but can be predicted with ${\lambda }_1$ difference for mutants 1 and 2.

Figure 4:. The basis for the atlas of multi-drug resistance evolution patterns over space.

A. The four canonical mutant types for resistance phenotypes to two drugs, distributed in the $(\alpha ,\beta )$ space of rescaling parameters: blue - fully resistant to drug 1 and sensitive to drug 2; red - fully resistant to drug 2 and sensitive to drug 1; purple - intermediate resistance to both drugs; brown - wild-type, sensitive to both drugs. B. The 3 drug fitness landscapes used: synergistic (left) independent (center) and antagonistic (right), as specified in Eqs. 16. These drug landscapes will be used to give rise to $g_i(z)=G\left({\alpha }_{i}x(z),{\beta }_{i}y(z)\right.$) as a function of two drug variation over space $x(z)$ and $y(z)$. The relative fitnesses of the strains are hence dependent both on drug variation over space and on the details of the underlying growth landscape $G$.

Figure 5:. An atlas for 2-drug resistance evolution in space under spatial heterogeneity.

We consider only four available mutants each with different resistance phenotypes to two drugs, distributed in the $(\alpha ,\beta )$ space of rescaling parameters: blue - fully resistant to drug 1 and sensitive to drug 2; red - fully resistant to drug 2 and sensitive to drug 1; purple - intermediate resistance to both drugs; brown - wild-type, sensitive to both drugs. We considered a diffusion coefficient of $D=0.01$; the spatial equilibrium is obtained numerically by considering the system at $t=1000$. For all the simulations, we considered the same initial distributions with 99% wild type mutants and the remaining 1% distributed equally among the three resistant mutants. Without loss of generality, the initial distributions of each mutant were shaped as the function $\sin \left(\frac{\pi z}{L}\right)$, so that the homogeneous Dirichlet boundary conditions were respected. The growth landscapes were as specified in Equation 16. The interaction strength is fixed at $q=0.5$ both in the case of synergistic and antagonistic interaction. For more drug gradient scenarios, under a more complex drug interaction profile and two diffusion rates see Supplementary Figures S1–S3.

Figure 6:. Validating selection predictions based on λ1 ranking among several competing mutants.

We illustrate a model simulation under the linear drug gradients in A, with 10 multi-drug resistant mutants varying in $\left({\alpha }_i,{\beta }_i\right)$ traits (B), growing $\left(g_i(z)\right.$ in C) and spreading over space with diffusion coefficient $D=0.01$. The ${\lambda }_1$ values (Eq. 14) in $\mathrm{D}$, match very well with the spatial selection dynamics observed numerically (E). Initial conditions (99% vs 1%: WT vs. all mutants) were assumed equal for all strains, satisfying the boundary conditions $n_i(z,0)\sim \mathit{sin}(\pi z)$.

Figure 7:. Selection outcomes for periodic drug regimes leading to periodic growth rates over space.

A. The periodic variation of drug 1 and drug 2 over space, keeping the total amount of each drug equal. The periodic function parameters, under conservation of total drug, are: $k_1=A_1=0.5,T_1=2$ and $k_2=A_2=0.72,T_2=0.4$. Further we show mutant growth rates and selection outcome under: B. synergistic drug interactions; C. independent drugs; C. antagonistic drug interactions; E. even more antagonistic drugs. The first column shows resulting growth rates $g_i(z)$ for each mutant following the linear $G$ function combinations in Eqs. 16, with $q=0.5$ (first three rows), and $q=1.5$ in the last row, depicting the case of stronger drug antagonism. The second column shows associated final fitnesses of the 3 mutants over space, computed on the basis of the principal eigenvalue. Assumed diffusion coefficient is $D=0.01$.

Figure 8:. Fitness ranking among 3 mutants for periodic drug regimes, as a function of spatial period of drug 2.

A. Synergistic drug interaction. B. Independent drug action C. Antagonistic drug interaction. The strength of interaction when assumed, was q = 0.5 and G(x, y) were specified as in Eq. 16. The periodic variation of drug 1 x(z) was held fixed, while drug 2 concentration y(z) over space was varied by varying the period T₂. These parameters were fixed: k₁ = A₁ = 0.5, T₁ = 0.8 and k₂ = A₂ = 0.5 before normalization, which then leads to total conservation of drug 1 and drug 2, fixed amount=1 for each spatial period of drug 2 T₂. Final fitnesses of the 3 mutants over space, computed on the basis of the principal eigenvalue. Assumed diffusion coefficient is D = 0.01. In blue: single-resistance to drug 1, in red: single-resistance to drug 2, in purple: double resistant mutant with intermediate resistance to each drug.

Figure 9:. Drug-resistance selection outcomes for periodic 2 drugs as a function of their spatial periods T1 and T2.

A. Synergistic drugs. B. Independent drugs. C. Antagonistic drugs. Shaded blue region: single-resistance to drug 1 has higher fitness, shaded red region: singleresistance to drug 2 has higher fitness, shaded purple region: double resistant mutant with intermediate resistance to each drug has the higher fitness. The periodic variations of drug 1 and drug 2 over one-dimensional space $z\in [\mathrm{0,1}]$ are constructed in such way as to keep the total amount of each drug equal to 1. The periodic function parameters are initially specified as: $k_1=A_1=k_2=A_2$ for any combination of periods $T_1$ and $T_2$, and then immediately scaled by the integral of the periodic function over space, to obtain a total amount of drug equal to unity in each case. Assumed diffusion coefficient is $D=0.01$. The growth functions of each mutant over space are obtained following Eqs. 16 together with the assumption that a mutant with traits $\left({\alpha }_i,{\beta }_i\right)$ experiences the two drugs at concentrations ${\alpha }_{i}x$ and ${\beta }_{i}y$. The interaction strength is fixed at $q=0.5$ both in the case of synergistic and antagonistic interaction. In the case of synergistic/antagonistic interaction, the effect is to decrease/increase the growth rate of bacteria relative to the simple additive effect of the two drugs. See Figure S4 for the analogous figure under a more complex drug interaction function, highlighting the sensitivity to fitness landscape.