May 20, 2018, Sunday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Stochastic Models for Bacterial Development of Multidrug Resistance

From MathBio

Jump to: navigation, search

Executive Summary

Here, we present stochastic simulations for the problem of bacterial development of multidrug resistance. We build off of the system presented by Obolski and Hadany [1], which uses an SIR format to examine rates of bacterial mutation when various antibiotic treatment strategies are used on hospitalized patients. We study the system with the therapeutic strategy of mixing, using a MATLAB script that encodes Gillespie's stochastic simulation method [2], in which random events (or reactions) are selected to proceed at random time intervals. We also take advantage of MATLAB's Simbiology Toolbox to perform deterministic, stochastic, and sensitivity analysis. We find that the stochastic solutions in Simbiology closely approximate the deterministic solution, and that emergence of multidrug resistant bacteria increases with the probability of stress-induced mutations.

Biological Phenomena Under Consideration

The model presented describes bacterial antibiotic resistance and the spread of a bacterial disease in a hospital environment. When patients are admitted to a hospital, they often contract bacterial infections unrelated to the reason they entered the hospital. Bacteria develop antibiotic resistance from both random selection and selection due to stress (exposure to antibiotics), and are capable of of developing multidrug resistance using these same mechanisms. Drug resistant strains are becoming increasingly prevalent for many infectious bacteria. Antibiotic resistance can impact patients by leading to increased mortality rates, prolonged hospitalization, and higher costs of treatment[1]. The rate at which bacteria mutate and develop antibiotic resistance far exceeds the rate that new antibiotic drugs can be developed; therefore, antibiotic resistance poses a serious heath problem.

Antibiotic resistance has always been present in bacterial colonies due to random mutation, but since the introduction of antibiotics, the resistance strains have been selected for. Everyone encounters antibiotics on a daily basis whether they realize it or not. Antibiotics are present in livestock feed and most commercially used soap in our homes and offices. Antibiotics are often times over prescribed and misused when they are given to patients. In some counties, antibiotics are available over the counter which leads to even greater misuse and extreme overuse. The high levels of antibiotic use in our society are causing mutated, resistant strains to proliferate and spread more easily.

Figure 1: Number of new beta lactamase enzymes after the introduction of beta lactamase antibiotics. [3]

Antibiotic resistance has been identified as one of the worlds most pressing health issues by the World Health Organization. Research is being performed to better understand antibiotic resistance and design drugs that can overcome resistance. An in vivo study on vancomycin showed 35 mutations in a bacterial genome in only 3 months [4]. Figure 1 shows data for the number of unique \beta -lactamase enzymes since the introduction of \beta -lactamase antibiotics. This study demonstrates the effect of stress induced mutations due to antibiotics. As the antibiotics were prescribed for years and used by patients, different mutated \beta -lactamase strains were more likely to survive and spread their mutation leading to a greater number of unique enzymes that may or may not be susceptible to \beta -lactamase antibiotics.

Heath care professionals are trying to combat drug resistance by using different therapeutic dosing strategies. Cycling is when the drug being used is periodically switched. By cycling the administered drug, bacteria that has become resistant to the first drug is likely to be susceptible to the second drug in the cycle. Mixing is when the patient receives a random antibiotic. Bacteria in this heterogeneous stress environment are then more likely to encounter a drug for which they have not yet developed resistance. Combining is a another tactic in which many drugs are given to one patient at a time, with the hope that all bacterium will be susceptible to at least one of the antibiotics. The model presented below re-examined the mixing strategy using a stochastic model and sensitivity analysis.

Mathematical Model

Deterministic Model

Obolski and Hadany [1] use four ordinary differential equations to describe patient conditions and population groups:

(1) {\frac  {dX}{dt}}=(1-\lambda _{s}-\lambda _{{R1}}-\lambda _{{R2}}-X)m+(\tau +\gamma )S+(\tau \chi _{2}+\gamma )R_{1}+(\tau \chi _{1}+\gamma )R_{2}-\beta X(S+R_{1}+R_{2})

(2) {\frac  {dS}{dt}}=(\lambda _{s}-S)m-(\tau +\gamma )S+\beta XS

(3) {\frac  {dR_{1}}{dt}}=(\lambda _{{R1}}-R_{1})m-(\tau \chi _{2}+\gamma )R_{1}+\beta XR_{1}

(4) {\frac  {dR_{2}}{dt}}=(\lambda _{{R2}}-R_{2})m-(\tau \chi _{1}+\gamma )R_{2}+\beta XR_{2}

Figure 2: Model system for emergence of multidrug resistance [1].

Where the parameters and the values we used in the stochastic models are as follows:

Parameter Value Dimension Comment
X 300 @t=0 Number Uninfected individuals
R_{1} 100 @t=0 Number Individuals infected with strain resistant to antibiotic 1
R_{2} 100 @t=0 Number Individuals infected with strain resistant to antibiotic 2
S 200 @t=0 Number Individuals infected with nonresistant strain
\tau 0.5 1/day Rate at which an individual becomes uninfected due to antibiotic use
\gamma 0.03 1/day Rate at which an individual becomes uninfected due to immune response
\beta 0.9 1/day Rate of bacterial transmissions resulting in infection
m 0.1 1/day Rate of patient turnover
\lambda _{{R1}} 0.1 -- Proportion of infected patients entering the hospital who carry bacteria resistant to antibiotic 1
\lambda _{{R2}} 0.1 -- Proportion of infected patients entering the hospital who carry bacteria resistant to antibiotic 2
\lambda _{{s}} 0.1 -- Proportion of infected patients entering the hospital who carry nonresistant bacteria
\chi _{{1}} 0.5 -- Proportion of patients receiving antibiotic 1
\chi _{{2}} 0.5 -- Proportion of patients receiving antibiotic 2

Additional parameters are introduced to account for mutations that result in antibiotic resistance:

  • \mu _{s}: Rate at which resistance-conferring mutations occur when bacteria are stressed
  • \mu _{r}: Rate at which resistance-conferring mutations occur when bacteria are free from stress
  • \sigma : Relative persistence of antibiotic resistance when there is no direct antibiotic usage, where 0<\sigma <1

To describe the emergence of a bacterial strain that is resistant to both antibiotic 1 and 2, an expression for the frequency of individuals carrying single-drug resistant strains due to stress induced mutation can be written:

E_{{SIM}}(U)=\int _{{t_{0}}}^{{t_{1}}}{\frac  {R_{1}(U)+R_{2}(U)}{R_{1}(U)+R_{2}(U)+S(U)+X(U)}}\,{\mathrm  {d}}t

Where U is the antibiotic treatment strategy (mixing simulated in the stochastic models). The emergence of the multidrug resistance strain can then be written in terms of E_{{SIM}}, the mutation rates (\mu ), and the resistance persistence (\sigma ):

\xi _{{SIM}}(mixing)={\frac  {1}{2}}(\mu _{s}+\mu _{r}\sigma )E_{{SIM}}(mixing)

Stochastic Model

The deterministic ODE model can be written as a set of chemical reactions for stochastic modeling. Reactions that import or export from outside of the compartment were modeled as zero order reactions. The \oslash represents the "null" participation of people outside the hospital in zero order reactions. Zero order reactions were appropriate for hospital intake and discharge because they were assumed to occur at a constant rate (i.e. independent of concentration or population).

Reaction Number Reaction Equation Phenomena Described
1 S{\overset  {m}{\rightarrow }}\oslash Rate of patients with susceptible bacteria leaving hospital
2 \oslash {\overset  {\lambda _{{s}}m}{\rightarrow }}S Rate of patients entering hospital with susceptible bacteria
3 R_{{2}}{\overset  {m}{\rightarrow }}\oslash Rate of patients with bacteria resistant to antibiotic 2 leaving the hospital
4 \oslash {\overset  {\lambda _{{R_{{2}}}}m}{\rightarrow }}R_{{2}} Rate of patients entering the hospital with bacteria resistant to antibiotic 2
5 R_{{1}}{\overset  {m}{\rightarrow }}\oslash Rate of patients with bacteria resistant to antibiotic 1 leaving the hospital
6 \oslash {\overset  {\lambda _{{R_{{1}}}}m}{\rightarrow }}R_{{1}} Rate of patients entering the hospital with bacteria resistant to antibiotic 1
7 \oslash {\overset  {(1-\lambda _{{s}}-\lambda _{{R_{{1}}}}-\lambda _{{{R_{{2}}}}})m}{\rightarrow }}X Rate of patients entering the hospital who are uninfected
8 X{\overset  {m}{\rightarrow }}\oslash Rate of patients who are uninfected leaving the hospital
9 S+X{\overset  {\beta }{\rightarrow }}2S Rate of infection of uninfected patients with susceptible bacteria
10 S{\overset  {\tau }{\rightarrow }}X Rate of patients with susceptible bacteria treated with antibiotic 1 or 2
11 S{\overset  {\gamma }{\rightarrow }}X Rate of immune clearance of patients infected with susceptible bacteria
12 R_{{2}}+X{\overset  {\beta }{\rightarrow }}2R_{{2}} Rate of infection of uninfected patients with bacteria resistant to antibiotic 2
13 R_{{2}}{\overset  {\gamma }{\rightarrow }}X Rate of immune clearance of patients infected with bacteria resistant to antibiotic 2
14 R_{{2}}{\overset  {\tau \chi _{{1}}}{\rightarrow }}X Rate of patients with bacteria resistant to antibiotic 2 being treated with antibiotic 1
15 R_{{1}}+X{\overset  {\beta }{\rightarrow }}2R_{{1}} Rate of infection of uninfected patients with bacteria resistant to antibiotic 1
16 R_{{1}}{\overset  {\gamma }{\rightarrow }}X Rate of immune clearance of patients infected with bacteria resistant to antibiotic 1
17 R_{{1}}{\overset  {\tau \chi _{{2}}}{\rightarrow }}X Rate of patients with bacteria resistant to antibiotic 1 being treated with antibiotic 2

Stochiometric Matrix

The stoichiometric matrix for the system of reactions is:


Reaction Number {\begin{matrix}1\\2\\3\\4\\5\\6\\7\\8\\9\\10\\11\\12\\13\\14\\15\\16\\17\end{matrix}} {\begin{bmatrix}-1&0&0&0\\+1&0&0&0\\0&0&-1&0\\0&0&+1&0\\0&-1&0&0\\0&+1&0&0\\0&0&0&+1\\0&0&0&-1\\+1&0&0&-1\\-1&0&0&+1\\-1&0&0&+1\\0&0&+1&-1\\0&0&-1&+1\\0&0&-1&+1\\0&+1&0&-1\\0&-1&0&+1\\0&-1&0&+1\end{bmatrix}}

Analysis & Results

Stochastic Algorithm

The algorithm for the code was adapted from Gillespie [2] as presented by Rawlings and Ekerdt [5]

1. Initialize. The integer counter must be set to zero and the intial conditions must be set. The stoichiometric matrix has to be entered for all reactions and the reaction probability laws or rate expressions must be entered for all reactions.

2. Compute the reaction probabilities and then compute the total reaction probability (r_{{tot}}).

3. Select two random numbers (p_{{1}}, p_{{2}}) from a uniform distribution on the interval (0,1). Set the time interval until the next reaction as follows:

                                                          {\check  {t}}={\frac  {-ln(p_{{1}})}{{r_{{tot}}}}}
Figure 3: Simple example of how algorithm weights reactions to choose the reaction m that occurs during the time interval. [5]

4. Determine the reaction to take place during this time interval using the reaction probabilities, the total reaction probability, and p_{{2}}. Figure 3 give a simple example of how the reaction probabilities are used to weight the chance that a faster or more favorable reaction will happen during a time interval.

5. Update the simulation time as t(n+1)=t(n)+{\check  {t}} and update the species numbers for a single occurrence of the selected reaction m the occurred during that time interval using the stoichiometric matrix as follows:

                                                 x_{{j}}(n+1)=x_{{j}}(n)+\nu _{{mj}}

Stochastic MATLAB Code

The MATLAB code generated for this modeling can be found in the following PDF: Media:Matlabcode.pdf

Figure 4 shows 20 stochastic simulations in MATLAB.

Figure 4: Twenty stochastic simulations plotted in MATLAB, number of time steps=2000.
Figure 5: 20 stochastic simulations in MATLAB displayed for the short initial time scale to show initial behavior.

Multidrug Resistance Emergence

Figure 6: 20 stochastic simulations in MATLAB for multidrug resistance emergence with mixing antibiotic administration.
Figure 7: Emergence of multiresistant strains, from Obolski and Hadany.

SimBiology Model

Deterministic and Stochastic Solutions

Figure 8: Deterministic solution, with the following initial conditions: S(0) = 200, R2(0) = 100, R1(0) = 100, and X(0) = 300.
Figure 9: Stochastic solution from Simbiology ensemble run. Ten runs with initial conditions as in deterministic solution.

This model was also created and studied using MATLAB's Simbiology Toolbox.

An ODE solver (ode15s) was applied to obtain a deterministic solution. An epidemic is observed, since the number of patients infected with resistant strains R1 or R2 increases before it decreases.

The results from Simbiology's Stochastic Simulation Algorithm (SSA) are shown in Figure 9. Note that the overall behavior observed on these ten runs is very similar to the behavior seen in the deterministic model. We again see epidemic increases in the number of R1 and R2 patients, and the number of uninfected (X) patients quickly drops.

Sensitivity Analysis

Simbiology was used to perform sensitivity analysis on the system. While Obolski and Hadany assumed that the hospitalization time (1/m) and clearance rates (\tau , \gamma ) were identical for all patient classifications, we allowed these values to be independent of one another, such that patient infected with R1 could have a different hospitalization time than an uninfected patient. By using independent rate constants for each reaction, we could evaluate the impact of changing all 17 rate constants in the model to trends in R1, R2, X, and S.

Figure 10 shows the magnitude of the time integral of each patient classification as each of the seventeen rate parameters are changed. The analysis reveals that for both R2 and R1 (not shown), manipulating the discharge parameters in Equation 3 and Equation 5 will have the largest impact on the value of R1 and R2. When the analysis was performed on S, the treating and healing rates in Equations 10 and 11 were found to have the greatest effect. The discharge rate in Equation 8 made the largest difference in the value of X.

Figure 10: Plots from a sensitivity analysis in Simbiology.

The scales on the horizontal axis for each of the three sensitivity plots provides insight into the diversity of the stochastic solutions. Note that the largest value in the R2 plot is nearly 3E4 for the R2 plot, while the X plot caps out at 4E2. In the stochastic ensemble runs, the R1 and R2 solutions do show larger variability, while the X solutions are tighter.

Discussion & Conclusions

The Simbiology analysis shows that the deterministic model used by Obolski and Hadany is a good approximation of the system. Ensemble runs of stochastic simulations do not show solutions that vary widely from the deterministic solution. When we look at the curves for emergence of doubly-resistant bacteria using the stochastic code, we see that, though the shape and magnitude of the curve varies, overall, it is in the same range as the curves presented in the paper [1] For the mixing case, where each patient is randomly given either Antibiotic A or Antibiotic B, there is a low dependence on the persistence of antibiotic resistance. Emergence of multidrug resistant strains depends primarily on the likelihood that a mutation will confer resistance.

Figure 11: Tau or time step as a function of the total reaction probability. As the reaction proceeds the total reaction probability approaches 0 and then time step changes greatly.

Stochastic simulations in MATLAB with the Gillespie code show similar behavior to stochastic simulations from Simbiology's SSA, but there are major differences. Both solution strategies show an increase in infected individuals (S, R1, and R2) close to t=0, then an overall decrease in infected individuals as time progresses. The Simbiology solutions show a faster decrease in the number of individuals infected with susceptible bacteria (as compared to individuals infected with R1 or R2 resistant bacteria), while the Gillespie code shows a slightly slower decay in S. This may be related to the time scales, which are different. The time scale in the MATLAB Gillespie code varies greatly as the time progresses due to the changing Tau. Figure 11 shows how the time step changes as the reaction proceeds. Initially there are many small time steps and as the reaction proceeds the time step gets larger therefore the reactions are sampled less later in time.

It is unclear and would need extensive further analysis to determine whether the MATLAB stochastic algorithm or Simbiology stochastic model more accurately represents the system.

The double resistance emergence stochastic graph, while on a different magnitude, show similar behavior. There seems to be a two order of magnitude difference in emergence between different stochastic simulations. This leads one to believe that it can be very difficult to study and understand antibiotic resistance considering each hospital as a different stochastic environment.

The final presentation of the results and model can be accessed: Media:finalpresentation.pdf


  1. 1.0 1.1 1.2 1.3 1.4 Obolski, Uri, and Lilach Hadany. "Implications of stress-induced genetic variation for minimizing multidrug resistance in bacteria." BMC medicine 10.1 (2012): 89.Online.
  2. 2.0 2.1 Gillespie, Daniel T. "Exact stochastic simulation of coupled chemical reactions." The journal of physical chemistry 81.25 (1977): 2340-2361.
  3. Center for Disease Control.
  4. Julian Davies and Dorothy Davies Microbiol. Mol. Biol. Rev. 2010, 74(3):417. DOI: 10.1128/MMBR.00016-10.
  5. 5.0 5.1 Rawlings, J.B., Ekerdt, J.G. Chemical Reactor Analysis and Design Fundamentals. Nob Hill Publishing, Dec. 2004 Madison Wisconsin.