A Model for Emergence of Multiple Anti-Microbial Resistance in a Petri Torus

TL;DR

This paper presents a lattice model using Ising perceptrons to simulate bacterial interactions with antimicrobial drugs, capturing resistance emergence and highlighting the benefits of drug mixing strategies.

Contribution

It introduces a novel statistical physics lattice model that links theoretical predictions with measurable in vitro experimental parameters.

Findings

Mixing drugs reduces resistance emergence compared to other protocols.

The model reproduces qualitative features of resistance development.

Parameters can be experimentally measured in vitro.

Abstract

This work introduces a new statistical physics lattice model of bacteria interacting with anti-microbial drugs that can reproduce qualitative features of resistance emergence and whose model parameters and outputs can be measured with controlled \textit{in vitro} experiments. The lattice is inhabited by agents modeled by Ising perceptrons. The results show the advantage of mixing drugs among the population compared to other treatment protocols.

Equations

$$q_{ij}^{\mu }=\Theta \left({\Delta }_{ij}^{\mu }\right)p_{ij}^{\mu }\left(1-e^{-{\lambda }_{ij}^{\mu }{\Delta }_{ij}^{\mu }}\right).$$

$${\bar{c}}_{ij}^{\mu }$$

$${\lambda }_{ij}^{\mu }$$

$$p_{ij}^{\mu }$$

$$\bar{c}=\lambda =\frac{x}{1-x},\text{and}p=1-x,$$

$${\sigma }_{+}^2$$

$${\sigma }_{-}^2$$

$$\rho (t)=\frac{1}{K{N}_a}\sum _{k,l}\text{mbox}sgn\left(\Delta {\bar{q}}_{kl}^1(t)\Delta {\bar{q}}_{kl}^2(t)\right),$$

Full text

I Introduction

Anti-microbial resistance (AMR), the resistance of infectious micro-organisms to therapeutic drugs, has been dramatically increasing during the last decades, becoming a threat as dangerous as climate change Woodford and Livermore (2009); Andersson and Hughes (2010) and potentially more urgent. With decreasing rates of anti-microbial drug (AMD) development Charles and Grayson (2004), we risk returning to a defenseless era against infections. The issue is so important that the World Health Organization has proposed a global action plan to address it WHO (2015).

The focus of most studies lies on bacteria, which can develop AMR in several ways. In addition to mutations in their main DNA and the possibility of horizontal gene transfer (HGT) through the intermediate action of phages (bacteria-infecting viruses), they also possess circular DNA fragments, called plasmids, which frequently encode resistance genes Ng et al. (2010). Plasmids can also mutate like the main DNA and they largely contribute to HGT San Millan et al. (2015); Baltrus (2013); Andersson and Hughes (2010) by being exchanged between bacteria of different species or freed into the environment for instance, upon cell death, once the cell’s membrane breaks down. This versatility makes plasmids powerful tools for AMR acquisition.

There are several attempted methods to deal with resistance, although none of them has been enough to alleviate the problem satisfactorily. One common practice is simply to avoid using of some AMD for which resistance has emerged for a (potentially long) amount of time. Although resistance seems to decrease in general with this protocol, its efficiency has been disputed Barbosa and Levy (2000). Evidence shows that the reversal rate is slow Austin et al. (1999); Tan et al. (2011) and, even after some reversal is observed, AMR never goes back to original levels.

Other common treatment protocols make simultaneous use of more than one AMD Bonhoeffer et al. (1997), opening the doors for the emergence of multiple resistance. The most popular are the two-drug protocols called combining, mixing and cycling. Combining AMDs amounts to administering the two drugs simultaneously to all patients. Mixing means that each drug is given separately to a certain fraction of the patients, who keep taking the same drugh during the whole treatment period. Finally, cycling entails the periodic alternation between drugs, each one given to the whole group at the same time.

There is a large literature on modeling AMR using differential equations Bonhoeffer et al. (1997); Alexander et al. (2007); D’Agata et al. (2008); Obolski and Hadany (2012); Ternent et al. (2015). From the perspective of statistical physics, these are effective models, i.e., averages over large populations. They have been used to predict and analyze different aspects of AMR. However, as it is common with effective models, they comprise a large number of adjustable parameters as many of them have to be included ad hoc to allow for a better adjusting of observed features.

In microscopic models, on the other hand, the number of adjustable parameters is usually reduced as they should emerge from microscopic quantities of the system. While effective models can be advantageous from the practical point, once their solution has usually a lower computational cost, microscopic models help in understanding basic principles behind the phenomenon, leading possibly to a better control and refinements or corrections to the effective models.

AMR and evolution of cell cultures have been addressed in different branches of theoretical physics before, but most models were again effective ones Patra and Klumpp (2014); Lambert and Kussell (2015); Bittihn et al. (2017). The use of microscopic models has been far less common Hermsen and Hwa (2010); Greulich et al. (2012). Analytical results from related statistical physics evolutionary models also exist, but AMR is not specifically addressed Tikhonov and Monasson (2017).

Given the importance of the subject, it is surprising that the powerful methods of statistical physics have not been more widely applied to the AMR problem. This paper contributes to fill this gap by introducing a new microscopic model of agents in a lattice. Using fewer parameters than models relying on differential equations, it is capable of reproducing qualitatively main features of AMR. It is clear that simplified models like the present one cannot reproduce perfectly all the subtleties of the biological processes involved in AMR. In this first approximation though, the aim of this work is to reproduce mainly qualitative aspects and understand the fundamental mechanisms at work. Quantitative agreement with experiments requires more thorough studies and model refinement.

In section II we introduce the lattice model representing an artificial bacteria culture in a Petri torus. The model has only three parameters representing a genotype-phenotype map from the cell’s plasmid to their response to AMDs. The introduced model is used to analyze the effect of the single drug protocol in section III, where it is shown that it can reproduce the memory effect we mentioned above. Section IV shows the results for double-drug protocol in which we single out the mixing protocol as the best option to slow down the emergence of resistance. Finally, further discussions and conclusions are presented in section V.

II Artificial Bacteria

The model presented here corresponds to a scenario in which a bacterial culture in a Petri dish is subjected to different treatments with different AMDs. The evolution of the culture in time is then followed in order to provide information about how each treatment affects the dynamical aspects of AMR emergence in the culture.

For convenience, the Petri dish is represented in a simplified way by a square $N\times N$ lattice with periodic boundary conditions (PBC) in both directions, which is going to be called a Petri torus due to the resulting topology. To each site $(i,j)$, $i,j=1,...,N$, a spin variable $\sigma_{ij}$ which is +1 if the site is occupied by a cell and -1 if it is empty. PBC imply $\sigma_{i+N,j}=\sigma_{i,j+N}=\sigma_{ij}$.

The main bacterial life cycle, excluding the action of the AMDs, is modeled by two parameters. At each time step $t$, cells reproduce by choosing one of its four neighboring sites with the same probability. If the site is empty, then there is a probability $r$ for the cell to divide with the newborn bacteria occupying the new site. A natural death probability $d$ for each living cell at time $t$ is used to model all non-drug related processes like other adverse environmental conditions, the natural life-cycle of the cell and a patient’s immune system.

Bacteria are modeled as agents living in the Petri torus subjected to external stimuli provided by local concentrations of the different AMDs. The drugs are administered according to each specific treatment protocol to be analyzed. Pharmaceutical companies usually measure the efficiency of particular AMDs by their Minimum Inhibitory Concentration (MIC), the drug concentration above which bacterial growth stops on average Davey et al. (2015). The MIC is convenient in clinical trials as it avoids the difficulties in isolating the effects of patients’ immune system, but it contains no information about pharmacodynamic properties of the drug (how bacterial growth changes with variations in drug concentrations). It has been proposed that a better proxy is given by the Minimum Bactericidal Concentration (MBC), the concentration that kills at least 99.9% of the bacteria on average French (2006). Here, a variation of the MBC is proposed and used: the concentration below which no cell dies. This definition is much simpler to implement in a probabilistic microscopic model and, because the studied scenario comprises in vitro cultures, it can be actually measured with controlled experiments.

Finally, it is assumed for simplicity that each cell has one single plasmid encoding its AMD response, which is given by the total AMD death probability $q^{\mu}_{ij}$, the probability that the cell occupying the site $(i,j)$ dies if exposed to the local concentration $c_{ij}^{\mu}$ of the $\mu$-th AMD, $\mu\in{\left\{1,...,M\right\}}$, where $M$ is total number of different available AMDs. This probability is modeled by

[TABLE]

where $\Delta^{\mu}_{ij}\equiv c_{ij}^{\mu}-\bar{c}_{ij}^{\mu}$.

There are only three parameters in the above expression, what is already a huge simplification if compared to the models based on differential equations. The MBC $\bar{c}_{ij}^{\mu}\in\left[0,\infty\right)$, enters in the model through the Heaviside step function $\Theta(x)$, which is zero if $x<0$ and 1 otherwise. Its objective is to account for physical and chemical mechanisms having threshold behaviors as, for instance, chemical pumps which can become saturated or membranes whose thickness up to a certain point can prevent the AMDs from entering the cell interior. The maximum death probability $p_{ij}^{\mu}\in{\left[0,1\right]}$ is a scaling that sets the maximum damage a certain AMD can inflict to the cell. Finally, the sensitivity $\lambda_{ij}^{\mu}\in\left[0,\infty\right)$ regulates the increase in cell death with AMD concentration and is related to the actual toxicity of it. All three parameters can, in principle, be obtained from actual designed experiments by measuring changes in bacterial populations.

The genotype of each plasmid will be encoded by a spin chain ${\boldsymbol{\pi}}_{ij}\in{\left\{\pm 1\right\}}^{D}$, where $D$ is an integer. The crucial point in epigenetics is to find an appropriate genotype-phenotype map into the model parameters Stadler and Stadler (2006); Cortini et al. (2016). The map should allow for different genotypes to generate the similar phenotypes and also for learning and forgetting. Learning is required for adaptation in the form of AMR acquisition and forgetting allows acquired adaptations to fade away, as would be the case for reversal of AMR. One of the most studied statistical physics models of biological phenomena with these characteristics is the Ising perceptron Rosenblatt (1958); Engel and van den Broeck (2001), designed to model neuronal responses, with the neuron’s synapses and stimuli both encoded by spin chains.

The perceptron is one of the simplest known learning machines and it is characterized by a function taking a multidimensional vector into a number, called the activation function for biological reasons. This is usually a general function of the scalar product of its synaptic vector, which is a multidimensional parameter encoding the information learned by the perceptron, and the input vector, a vector with the same dimensions as the synaptic vector and which encodes the stimuli provided by the environment to which the perceptron reacts.

In this work, the input vectors correspond to the binary encodings of each AMD into a spin chain of dimension 3D, ${\boldsymbol{A}}^{\mu}={\left({\boldsymbol{\alpha}}^{\mu},{\boldsymbol{\beta}}^{\mu},{\boldsymbol{\gamma}}^{\mu}\right)}\in{\left\{\pm 1\right\}}^{3D}$. The three model parameters are then given by

[TABLE]

which are chosen for being the simplest mappings into the relevant intervals.

Although simple perceptrons cannot approximate general functions, it has been shown that adding one extra layer, corresponding to another set of perceptrons realising an intermediate processing of information between the stimuli and the final response, turns them into universal approximators Cybenko (1989). Variations with several layers, known as deep neural networks, have been successfully used in machine learning applications and recently provided a solution for the long sought problem of creating a computer algorithm capable of playing Go on a level comparable to human masters Silver et al. (2016).

The dimension of the AMD vector was chosen to be three times larger than the plasmid’s one to allow mutations in the latter to simultaneously affect all three model parameters. In nature, each protein usually participates in more than one metabolic process simultaneously. As a consequence, each single mutation might affect more than one of them.

One should notice that, although $p$, $\bar{c}$ and $\lambda$ are being called the parameters of the model, they are actually not numeric adjustable parameters. Once the genotype and the AMD are chosen, they are fully determined. More precisely, the only freedom comes from the choice of the activation functions leading to them. One might call the latter the functional parameters of our model, as it remains open the possibility of choosing more complicated genotype-phenotype maps which might lead to better agreement between observation and theory if needed.

III Single Resistance

In order to obtain the dynamical behaviour of the AMD death probability $q$ for one single drug, the model is simulated in discrete time $t$. Each initial Petri torus occupation is set by putting a cell in each site with probability 1/2. All initial cells carry the same all-ones plasmid ${\boldsymbol{\pi}}={\boldsymbol{1}}\equiv(1,1,...,1)$. To isolate the effect of the AMD, we will use $r=1$ and $d=0$, i.e., cells always reproduce if there is space and do not die unless killed by the AMD. With these two parameters fixed to these values, the dynamics then follows two steps at each $t$:

(1) Reproduction: all living cells are drawn once and only once with the same probability. Then, one of the four neighbors of that cell is chosen with probability 1/4 and, if the corresponding site is empty, the cell generates a child cell on it. The child cell has a probability $m$ of mutation, where one single coordinate of its plasmid, chosen at random with equal probability, is flipped.

(2) AMD Death: in a uniformly random order, one checks whether each cell dies according to the probability $q_{ij}^{\mu}$ for each AMD present locally.

Each run of the simulation consists of $T$ time steps and the information recorded is a double average of the death probability $q^{\mu}_{ij}$ - the average over all living cells and the quenched average over initial configurations.

According to our model, bacteria acquire resistance by (i) increasing the MBC, (ii) decreasing the maximum death probability, (iii) decreasing their sensitivity to the drug or any combination of these. Although the same mechanisms can contribute to different parameters, the larger the plasmid dimension $D$ is, the more uncorrelated we expect them to be.

It is clear that for large dimensions $D$, which is the case in practice, the fitness landscape will be highly AMD dependent. It is reasonable to expect it to b also very complex, with several local minima and maxima which cannot be easily studied only with analytical methods. In order to get some insight into these landscapes, we start our analysis by illustrating how one can use the model to extract useful information in simple cases. One such case is the choice $\alpha_{i}=\beta_{i}=-1$ and $\gamma_{i}=1$, $i=1,...,D$, in which the parameters can be expressed as

[TABLE]

where $x=n/D$ and $n$ is the fraction of -1 coordinates in the plasmid. When $D\rightarrow\infty$, $x$ becomes a continuous real variable in the interval $[0,1]$. In this case it becomes very simple to plot the fitness landscape as a contour plot of $q$ as a function of $x$ and the AMD concentration $c$ as given in fig. 1.

In this particular case, $\bar{c}$ and $\lambda$ are equal by design, which will usually not be the case in practice. Still, qualitatively this is an important situation where they both contribute in different ways for the emergence of AMR. While a higher MBC improves resistance, a higher sensitivity decreases it. This is reflected in the fitness landscape by the maximum of $q$ surrounded by descending profiles both to the left and right. Because this landscape has a only one global maximum, we can more clearly see that resistance will eventually emerge as we move away from it in the $x$ direction, which is always observed in simulations.

For the rest of this work, we use AMDs with $D=70$ where ${\boldsymbol{\alpha}}$ and ${\boldsymbol{\beta}}$ are generated randomly and uniformly. Without loss of generality, we set $\gamma_{i}=1$ to tune the initial value of $p$ to 1. Fig. 2 shows the average death probability $q$ for mutation rates $m=0.4,0.1,0.04$ in a lattice of size $N=100$. Results are averaged over $K=100$ realizations of the same AMD, differing by the initial site occupation, and over $N_{a}=10$ different randomly generated AMDs. Averaging over AMDs drastically increases the deviations from the mean curves as the response to each might be very different for the same plasmid. The local concentrations of AMD are kept constant during simulations, fine-tuned to give an initial $q=0.5$. As expected, larger mutation rates allow for faster adaptation (smaller $q$). The shadowed areas around the mean curves represent deviations from mean plots and, as discussed before, it is clear the dramatic increase in their amplitude as the mutation rate is decreased.

Because the death probability is restricted to the interval $\in[0,1]$, deviations above and below the mean are calculated separately. Let $\bar{q}_{kl}$ be the average death probability for the $k$-th run of the $l$-th AMD. Deviations above and below the mean are respectively

[TABLE]

where $\Delta q_{kl}=q-\bar{q}_{kl},$ $n_{+}$ and $n_{-}$ are the number of cases in which the variation is respectively non-negative or negative.

Deviations become larger for smaller mutation rates $m$ because different AMDs lead to very different rates of emergence. With higher $m$, cells probe larger areas of the fitness landscape before dying, resulting in faster growth of AMR. Higher initial concentrations require greater $m$’s for adaptation to win over extinction. The exact relationship between this threshold and the other parameters of the model requires a better characterization of the phase diagram, which is under study. Realistically, higher concentrations become toxic also to the patients and, although efficient to kill bacteria, might result in potentially fatal side effects.

Fig. 3 reproduces the observed memory effect that stopping the treatment does not restore resistance to its levels before its beginning. The same parameters as before are used, but only for $m=0.4$, averaging over 10 AMDs. Stopping times are written next to the corresponding curves, the bottom one representing a continuing treatment.

IV Double Resistance

While it seems logical to use several AMDs simultaneously, as the probability of being resistant to all of them is surely smaller than to a single one, this creates a multiple selection pressure that might lead to the emergence of resistance to all of them. To study this, we simulate the currently most common two-drug protocols: (i) combining, (ii) mixing and (iii) cycling.

Fig. 4 compares the rate of emergence of resistance when only one AMD (Single) is used against that of the three two-drug protocols. The curves are averages over 100 random initial occupations of the lattice and 10 randomly generated pairs of AMDs. Each AMD is generated independently, its coordinates being $\pm 1$ with equal probabilities. The plot shows only the first AMD of each pair. The results for the second vary slightly quantitatively, but the qualitative behavior is similar.

All protocols were applied using the same initial death probability, $q=0.5$ for the upper plot and $q=0.3$ for the bottom one. $N=100$ and $m=0.4$ for both. Again, shaded areas around curves represent deviations above and below mean values.

For the combining protocol, both AMDs are present at each site simultaneously. Concentrations are site independent and do not vary with time. For the mixing protocol, only one AMD is present at each site. For every site, we choose one AMD of the pair with probability 1/2 and keep its concentration constant in time. All sites with the same AMD have the same drug concentration. For cycling, the simulation starts with the first AMD in each pair present in all sites at the same concentration. After each cycling period $C$, the AMD is changed by the other one in the same pair in all sites. Results are presented for three different periods, $C=10,50$ and 100.

For the chosen parameters, combining shows very little emergence when $q=0.5$. However, this is mainly a result of the double concentration to which cells are exposed. Because the initial $q$ is relatively high for both AMDs, cells die faster than they can adapt. Mixing seems to be the next best, with an apparent asymptotic stop of resistance emergence. Although the cycling protocols do better than the single drug case, they seem to be unable to avoid the emergence of total resistance, which in this case occurs for both AMDs. In particular, this is independent of the cycling period, whose only effect is to create larger oscillations about an average line roughly coinciding with the smaller cycle. For cycling, we see again the memory effect discussed in the introduction: when one AMD is not being used, resistance decreases but not enough to bring a complete recovery to initial levels. It is not clear if there is any influence of the second AMD in the recovery, an interesting question for further studies.

For the plot in fig. 4b, simulations start with $q=0.3$ by appropriately adjusting the AMD concentrations. Differently from the previous case, combining becomes the worst of the multi-drug protocols, worse even than the single-drug one after some time. This change is explained by noticing that with $q=0.3$, cells live long enough to allow AMR to appear through mutations. The data indicates again that, while other protocols result in the emergence of an almost total resistance, mixing saturates at a lower level for the given parameters.

Fig. 5 shows the same graph for the two drugs of the pair for a cycling protocol with period $C=50$. The first AMD in the pair is represented by the line which drops faster in the initial period.

We measured the joint pharmacodynamics of a certain AMD pair by the average relative sign of their variations of $q$

[TABLE]

where ${\mbox{sgn}}\,x=x/|x|$ if $x\neq 0$ and 0 otherwise, $\Delta\bar{q}^{\mu}_{kl}(t)=\bar{q}^{\mu}_{kl}(t)-\bar{q}^{\mu}_{kl}(t-1)$ and $\bar{q}^{\mu}_{kl}(t)$ is the death probability averaged only over lattice sites for the $\mu$-th AMD in the $l$-th pair during the $k$-th simulation run at time step $t$.

This quantity, plotted in fig. 6 for each two-drug protocol (using the same parameters as in fig. 4b), shows how often on average the cells have the same adaptive trend towards resistance to the two drugs in the pair. The value +1 means that both $q$ change in the same direction, i.e., bacteria are either becoming resistant or more susceptible to both AMDs simultaneously. Conversely, -1 means that their adaptive trend is opposite, they are gaining resistant to one of them and losing to the other.

The results indicate that when both AMDs are simultaneously present, resistance changes in the same direction for longer than in cycling, but this is not enough for the latter to outperform mixing. Mixing remains superior because, when reproducing, cells will half of the time spread to a site with a different AMD, to which they are less resistant. Although this seems similar to cycling, the latter has the whole population subjected to the same AMD at every time step.

We now analyse a new protocol, mixed cycling, which in principle should combine the advantages of mixing and cycling: we distribute two AMDs uniformly randomly to the population and then cycle the resulting configuration. Fig. 7 shows a comparison between this protocol and mixing.

The plot, in which the same parameters as in fig. 4b were used and where the cycling period was chosen to be $C=50$, shows the results for the first AMD in the pair. One can see that mixed cycling ends up performing worse than mixing. Unexpectedly, instead of observing a reversal in resistance at each cycling period, the change in AMD has exactly the opposite effect: it increases the AMR emergence rate. The mechanism at work here is subtle, but at the same time very simple. With pure cycling, each period sees adaptation to only one AMD, while the other loses it. Resistance is reversed because no cell is under pressure to adapt to the absent AMD.

On the other hand, in every cycling period of mixed cycling, half of the cells are subjected to one AMD and half to the other. This implies that adaptation to both, although at a smaller rate, is always taking place. Still, that does not explain the whole story as, in the mixing protocol, both AMDs are also present all the time during the runs. The difference is that each cell which has already adapted in mixing, will not suffer any different pressure as the AMD at that site never changes. On the other hand, with the cycling, it is force to re-adapt. That, instead of decreasing AMR, works in its favor. The details of this mechanism need to be studied more carefully and is the subject of our current research.

V Conclusions

The model presented here has many advantages complementary to usual ones based on differential equations. Being a microscopic statistical model, there is larger control over mechanisms for AMR emergence. It also benefits from clarity and versatility while, at the same time, it can reproduce qualitatively actual features of both single and double resistance to arbitrary AMDs using very few parameters.

Additional mechanisms, like exchange or acquisition of plasmids, are easy to include in the model. One exciting possibility is to use machine learning algorithms to encode the structure of actual AMDs and study them. The plasmid genotype can be directly translated to binary code and the macroscopic parameters of the model can be obtained from experiments. Although perceptrons are too simple to approximate general genotype-phenotype maps, it was proven Cybenko (1989) that more complex networks, as deep networks Silver et al. (2016), are universal approximators and can become powerful tools in the search of real new AMDs and evaluation of resistance scenarios.

By applying the model to the most commonly used treatment protocols, we showed that the model (i) is capable of producing single and double resistance and (ii) can reproduce the memory effects observed in practice which indicates that reversal rates in the absence of AMD are inefficient to recover initial levels of resistance.

The simulations that have been carried out indicate that the mixing protocol is the most efficient of the protocols. The reason seems to be deceptively simple: while in all other protocols the whole population is subjected to the same AMD for some period of time, resistance grows faster during these times as all bacteria evolve in the same direction. With mixing, half of the population adapts mainly to one AMD and the other half to another. In addition, whenever bacteria spread, they have a chance of meeting a different AMD than the one in its original site, which changes the adaptation trend. Evidence in favour of this explanation comes from the proposed mixed cycling protocol. One could expect that this would benefit from the advantages of both mixing and cycling. In particular, the reversal in resistance observed at each cycling period. However, because of the mixing in the population, when the cycling period is reached, half of the bacteria are still subjected to the same AMD and, although slower, adaptation continues leading to resistance. On the other hand, contrary to mixing, all bacteria are exposed to both AMD, although in different periods. Because of this, cells at a definite site have contact to both AMDs instead of only one as in mixing, which leads to higher average adaptation to both.

Other statistical physics issues need further analysis as, for instance, finite size effects and the existence or not of true phase transitions to resistant phases. HGT might be an important issue as it is known to change the critical transition point for speciation Park and Deem (2007) and should also affect significantly AMR emergence.

Many qualitative aspects of the model are fairly general and actual experiments are fundamental in testing its predictions as well as identifying aspects in which a more sophisticated modeling is needed in order to allow for a closer approximation of practical situations.

Acknowledgements

I would like to thank Dr J. Neirotti, Dr M.Stich, Dr M.Chli and Dr A. Cheong for their suggestions and discussions.