
MBW:Stochastic Models for Bacterial Development of Multidrug ResistanceFrom MathBioContentsExecutive SummaryHere, 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 stressinduced mutations. Biological Phenomena Under ConsiderationThe 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 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 lactamase enzymes since the introduction of 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 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 lactamase antibiotics.
Mathematical ModelDeterministic ModelObolski and Hadany ^{[1]} use four ordinary differential equations to describe patient conditions and population groups: (1) (2) (3) (4) Where the parameters and the values we used in the stochastic models are as follows:
Where 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 , the mutation rates (), and the resistance persistence ():
Stochastic ModelThe 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 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).
Stochiometric MatrixThe stoichiometric matrix for the system of reactions is: Reaction Number Analysis & ResultsStochastic AlgorithmThe 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 (). 3. Select two random numbers (, ) from a uniform distribution on the interval (0,1). Set the time interval until the next reaction as follows: 4. Determine the reaction to take place during this time interval using the reaction probabilities, the total reaction probability, and . 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 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: Stochastic MATLAB CodeThe MATLAB code generated for this modeling can be found in the following PDF: Media:Matlabcode.pdf Figure 4 shows 20 stochastic simulations in MATLAB.
Multidrug Resistance EmergenceSimBiology ModelDeterministic and Stochastic SolutionsThis 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 AnalysisSimbiology was used to perform sensitivity analysis on the system. While Obolski and Hadany assumed that the hospitalization time (1/m) and clearance rates (, ) 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. 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 & ConclusionsThe 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 doublyresistant 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.
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 final presentation of the results and model can be accessed: Media:finalpresentation.pdf References
Appendix 