Identifying determinants of persistent MRSA bacteremia using mathematical modeling

Abstract

Persistent bacteremia caused by Staphylococcus aureus (SA), especially methicillin-resistant SA (MRSA), is a significant cause of morbidity and mortality. Despite susceptibility phenotypes in vitro, persistent MRSA strains fail to clear with appropriate anti-MRSA therapy during bacteremia in vivo. Thus, identifying the factors that cause such MRSA persistence is of direct and urgent clinical relevance. To address the dynamics of MRSA persistence in the face of host immunity and typical antibiotic regimens, we developed a mathematical model based on the overarching assumption that phenotypic heterogeneity is a hallmark of MRSA persistence. First, we applied an ensemble modeling approach and obtained parameter sets that satisfied the condition of a minimum inoculum dose to establish infection. Second, by simulating with the selected parameter sets under vancomycin therapy which follows clinical practices, we distinguished the models resulting in resolving or persistent bacteremia, based on the total SA exceeding a detection limit after five days of treatment. Third, to find key determinants that discriminate resolving and persistent bacteremia, we applied a machine learning approach and found that the immune clearance rate of persister cells is a key feature. But, fourth, when relapsing bacteremia was considered, the growth rate of persister cells was also found to be a key feature. Finally, we explored pharmacological strategies for persistent and relapsing bacteremia and found that a persister killer, but not a persister formation inhibitor, could provide for an effective cure the persistent bacteremia. Thus, to develop better clinical solutions for MRSA persistence and relapse, our modeling results indicate that we need to better understand the pathogen-host interactions of persister MRSAs in vivo.

Fig 1. Dynamics of in vitro bacterial growth treated with antibiotics, vancomycin.

(A) Schematic diagram of the mathematical model for in vitro growth. The diagram depicts the two state variables as blue ovals, normal and persister SA, whose abundances are controlled by growth rate, death rate by vancomycin (Van), and switch rates. The mathematical model is able to recapitulate not a phenomenon of resistance, but an in vitro persistence against antibiotics. (B) Heatmaps showing the number of normal, persister and total SA at 0, 24, and 120 h post-vancomycin addition as a function of switching rate sw_N2P and sw_P2N. The pre-simulations were run with N(0) = 1 colony-forming unit per arbitrary unit (cfu/a.u.) and P(0) = 0 cfu/a.u. without antibiotics for 8 h, then vancomycin was added (0 h post-vancomycin addition). The parameter values used are shown in Table 1. (C) The number of total, normal and persister SA graphed as a function of time for each combination of sw_N2P and sw_P2N.

Fig 2. Schematic diagram of the mathematical model of in vivo MRSA populations and ensemble modeling.

(A) The diagram depicts the two state variables as blue ovals, whose abundances are controlled by growth, death by antibiotics, clearance by the immune system (Im), as well as switch rates. Antibiotic treatment, vancomycin (Van), affects growth rates of both normal and persister and killing of the normal but not persister bacterium. The two differential equations are shown below the diagram. Van and Im are a binary value of 0 or 1 indicating absence or presence. (B) Algorithm for ensemble modeling to extract the sets of parameters (ensemble) using two criteria. All the models extracted had similar minimum inoculum doses to establish infection between around 100 to 1,000 colony-forming unit per arbitrary unit (cfu/a.u.). The parameter values randomized or fixed during the ensemble modeling are shown in Table 2.

Fig 3. Simulation of the models with vancomycin therapy.

(A) The selected parameter sets were used for simulation with vancomycin treatment and classified into resolving (RB) and persistent bacteremia (PB) as described in the algorithm. Vancomycin was administered when normal SA (N) exceeded 10⁴ cfu/a.u. The models were judged as PB when the number of total SA (T) at 5 days post-treatment exceeded a detection limit under an assumption where positive bacteremia is observed when the number of SA at infection foci exceed a certain number (see text in details). The detection limit of bacteremia is set as 0.1 cfu/a.u. in this study under an assumption where relapse bacteremia is less in RB (see Fig 3B). The relapse bacteremia was determined when vancomycin was re-administered triggered by regrowth of SA after 1^st vancomycin treatment. Four cases of typical simulations, RB no relapse, RB with relapse, PB no relapse and PB with relapse, are shown in the bottom. (B) Bar plots show the number of models with or without relapse for RB (left) and PB (right) over a range of detection limit from 0.001 to 10 cfu/a.u. A dotted line indicates the total number of models.

Fig 4. Comparison of models between resolving and persistent bacteremia.

(A, B) Line plots of the number of total, normal, and persister SA over time (0–60 days) for resolving (RB) and persistent bacteremia (PB). The detection limit of bacteremia was set as 0.1 cfu/a.u. and 2,479 and 2,135 models were classified into RB and PB, respectively. In RB, 52 of models showed relapsing bacteremia and they were also included in the analysis. (C) Two-dimensional scatter plot of switching rates. (D) Violin plots and box plots of the parameter distributions of RB and PB. Dotted lines indicate the range of randomized parameter values.

Fig 5. Key determinants to distinguish resolving and persistent bacteremia.

(A) Pearson correlation coefficients were calculated between each feature and the binary response variable (0: RB, 1: PB). The absolute values of the coefficients are shown. (B) Weights calculated by QPFS methods are shown. (C) The sole ranked feature by QPFS, c_P, was used to build a logistic regression model. For the comparison, all parameter and all parameter except for c_P were applied to the classification model. Best accuracy, ROC curve, and area under ROC curve are graphed for each model. (D) Two-dimensional scatter plot for c_P and sw_N2P, which showed 2^nd high correlation with target response (Fig5 A) with probability density for RB and PB.

Fig 6. Key determinants to distinguish resolving, persistent and relapsing bacteremia.

(A) Line plots of the number of total SA over time (0–180 days) for resolving bacteremia (RB) without relapse, persistent bacteremia (PB) without and with relapse bacteremia. (B) Violin plots and box plots of the parameter distributions. Dotted lines indicate the range of randomized parameter values. (C) Weights calculated by QPFS methods are shown. The ranked features by QPFS, c_P and g_P were used to build a multinomial logistic regression model. For the comparison, all parameter and all parameter except for c_P and g_P were applied to the classification model. The overall accuracy and area under ROC curve for each class are graphed. In the classification, an equal number of datasets in each class were used by reducing the number of data in major classes, RB without relapse and PB with relapse, to the same number of minor class, 940 of PB without relapse. (D) Two-dimensional scatter plot for c_P and g_P with probability density of each value.

Fig 7. Pharmacological strategies for persistent and relapsing bacteremia.

(A) Three possible pharmacological strategies, persister killer, persister reverter and persister formation inhibitor, are shown on the schematic diagram of in vivo MRSA model. (B) Doughnut charts show the number of models for resolving bacteremia (RB) without relapse, RB with relapse, persistent bacteremia (PB) without relapse, and PB with relapse by the treatment of three drug candidates in the presence of vancomycin. The pharmacological strengths of persister killer and persister reverter are expressed as a kinetic rate (h^-1). During the simulation, persister killer and persister reverter were administrated in conjunction with vancomycin. Persister formation inhibitor was treated from the beginning of the simulation with different initial numbers of persister cells, 1 or 0 cfu/a.u.