Modeling MRSA decolonization: Interactions between body sites and the impact of site-specific clearance

Nasal carriage of MRSA is the most persistent and most affected by the decolonization protocol

For each body site, we estimated the persistence of colonization by calculating the probability that the site would be colonized in the next time step when it was colonized in the previous time step, assuming other sites were not colonized. We found that in the education group, the estimated persistence was highest in the nares (Figure 2a). In the decolonization group, application of the protocol decreased persistence at each site, with the largest decrease in the nares[21] (Figure 2a). We also estimated the relapse probability, i.e., the probability of a site getting colonized when it was not colonized in the previous time step, assuming other sites were not colonized. However, this did not seem to be significantly affected by the decolonization protocol (Figure 2b). The number of samples from wounds was relatively small in the analysis, which caused the estimated parameters for the wound to have considerably wider CIs than estimates for the other sites. Repeating the analysis only on patients who reported adherence to the protocol led to almost identical results since most (70%) of the data were collected from adherent patients (Supplementary Materials).

• Download figure
• Open in new tab

Figure 2 Estimated (a) persistence and (b) relapse probability of MRSA colonization by body site.

(a) Persistence probabilities are defined as the probability that a site will be colonized in the next time step when it is colonized in the current time step. (b) Relapse probabilities defined as the probability of a site getting colonized when it was not colonized in the previous time step. Persistence and relapse probabilities are calculated assuming other sites are not colonized. The calculation of the posterior distributions is explained in the Methods. Means and 90% credible intervals (CI) are represented in the figure by diamonds and lines, respectively.

The decolonization protocol reduces the transmission of MRSA between body sites

We quantified the strength of dependencies of MRSA colonization between body sites using the increase in the log-odds of colonization in the next time step caused by colonized sites in the previous time step (Figure 3). First, as expected, the most important determinant of colonization at any site was whether the site itself was colonized in the previous time step. Second, the nares was the only site linked to all the other sites in the education group, highlighting it as the likely source from which MRSA spreads to the other sites. Third, the decolonization protocol decreased the impact of the nares on the throat and removed the impact of the nares on wounds (Figure 3). Except for a minor increase in the impact of the nares on the skin, other dependencies between body sites were not affected by the treatment. These findings indicate two mechanisms responsible for the efficiency of the decolonization protocol: 1) the protocol decreased persistence in the nares, which is a central hub for colonization (Figure 2a), and 2) the protocol weakened links between body sites, reducing self-inoculation (Figure 3).

• Download figure
• Open in new tab

Figure 3 Estimated interactions among body sites in the (a) education and (b) decolonization groups.

The red edges in the graphs show the increase in the log-odds of a site being colonized in the following time step if the linked site is colonized. Note there were no decreases in the log-odds. The edge thickness represents the expected change in the log-odds, and edges are excluded from the graph if the possibility for no effect is included in the 90% credible intervals (CI). Calculation of the log-odds is explained in the Methods. Tables show the mean and 90% CIs for all relations. We note that the number of wound samples was relatively small, which caused some wound-related parameters to have considerably wider CIs.

The CHMM accurately predicts the decrease in MRSA carriage over the study period

The key metric of success for decolonization is the extent to which it reduces MRSA carriage. To examine the CHMM’s predictive accuracy, we simulated patient trajectories using parameters from the estimated models and compared the decrease in colonization in our model-based simulations to the observed decrease from the CLEAR trial data. The predicted reduction of the carriage was from 66% to 48% (90% CI is 45%-51%) in the education group and down to 26% (23%-28%) in the decolonization group (Figure 4). These predictions accurately match the reductions observed in the trial. However, a discrepancy between predictions and observations was seen in intermediate visits. Especially in the decolonization group, the observed carriage rate decreases more rapidly than predicted in the beginning and then more slowly. This outcome may reflect heterogeneity in treatment response among the study population: some patients responded very fast to the treatment, others more slowly.

• Download figure
• Open in new tab

Figure 4 The observed decrease in MRSA carriage over time compared to the decrease predicted by the model for actual and hypothetical treatments.

In the trial, study subjects were in the education group or the decolonization group, which consisted of applying mupirocin to the nares, chlorhexidine mouthwash (CHG Oral) to the throat, and chlorhexidine body washes (CHG Skin) to the skin and, if present, wound. The predictions for hypothetical treatments assumed either that only a subset of medications was applied or that the decolonization protocol was used with an additional optimized intervention that resulted in immediate clearance of the corresponding target site. Dotted lines and shaded regions represent the mean and 90% credible intervals (CI) of simulations, and the solid lines represent values observed in the data.

Treatments on nares, skin, and wound account for 97% of the effect of full decolonization

The importance of the nares has previously been noted in data from an intensive care unit[22]. However, the added value of decolonization efforts that focus on body sites other than the nares is not well understood. To quantify the impact of decolonization efforts at each individual site, we predicted the efficiency of hypothetical simplified treatments, consisting of only a subset of the decolonization protocol (Figure 4). Nasal mupirocin decreased the estimated carriage from 66% to 31% (29%-34%), compared to 26% (23%-28%) of the full decolonization (Figure 4). Other single treatments were far less effective, reducing the carriage to 47% (44%-50%) with chlorhexidine mouthwash (CHG Oral) or 44% (41%-47%) with chlorhexidine bodywash (CHG Skin), a minor improvement compared to the 47% (44%-50%) of education alone. The best combination of two treatment components includes mupirocin and CHG Skin, i.e., leaving out CHG Oral, decreased carriage to 27% (25%-30%), corresponding 97% (90%-99%) of the decrease caused by the full decolonization protocol. In contrast, the combination of CHG Skin + CHG Oral and mupirocin + CHG Oral decreased the carriage only to 43%(40%-46%) and 30%(27%-33%), respectively (Suplementary Materials).

Enhancing clearance at the nares is predicted to achieve the largest gain in overall decolonization success

We used our model to predict the impact of future improvements to components of the decolonization protocol. We assumed optimal treatment components, defined as treatments that resulted in immediate clearance of a specific site, could be applied on top of the existing protocol. The results were consistent with the previous finding of the nares’ central role in colonization and its clearance. In particular, the estimated final carriage with an optimized clearance of MRSA nasal carriage was 14% (13%-16%), while optimizing the throat or skin/wound treatments decreased the carriage to 21% (19%-23%) or 20% (18%-22%), respectively (Figure 4). Optimizing treatments both on nares + throat has the potential to reduce carriage further down to 10% (9%-11%) compared to 11% (10%-12%) of the nares + skin/wound, and 16% (14%-18%) of the throat + skin/wound (Suplementary Materials).

Isolate Collection

MRSA isolates were collected as part of the CLEAR (Changing Lives by Eradicating Antibiotic Resistance) Trial. The trial was designed to compare the impact of a repeated decolonization protocol plus education on general hygiene and environmental cleaning to education alone on MRSA infection and hospitalization [3]. Study subjects in the trial were recruited from hospitalized patients based on an MRSA positive culture or surveillance swabs. After recruitment, swabs were obtained from different body parts of subjects (nares, skin, throat, and wound) around the time of hospital discharge (ENRL) and at 1,3,6 and 9 months (V1-V4, respectively) following the initial visit. The decolonization lasted only for 6 months, and consequently, we only modeled visits until V3 in this study. We note that some enrolled study subjects, despite a history of MRSA colonization or infection, did not have swabs positive for MRSA at the first time point (ENRL). More details about sample collection and data preprocessing can be found in the [3].

Hidden Markov Model

Markov models are a well-known tool for time series analysis. Hidden Markov Models (HMM) are a special case where Markov process is observed indirectly through noisy observations. Discrete-time discrete-state Hidden Markov model with a latent sequence π, π_t ∈ {1, …, K} and observations x, x_t ∈ {1, …, L} for t ∈ {1, …, τ}, is defined as where K is number of latent states, τ is the number of time steps, and L the number of possible observed states. Parameters of the model are the initial state probability π₀, transition probability T, and emission probability E, defined as

Note that the rows of T and E are probability distributions by definition. We will jointly denote the model parameters as θ = {π₀, T, E}. If θ is known, then we can estimate the latent states π conditional on parameters θ and observations x using the forward-filtering backward-sampling algorithm.

Forward-Filtering Backward-Sampling

Filtering is the estimation of the current hidden state by using all observations so far, p (π_t| x_1:t) [12]. This is known as forward-filtering and it can be represented as follows:

We define α (π_t) = p (π_t | x_1:t) and substitute that into Equation 6, which yields the following recursive equation known as α-recursion or forward recursion;

The recursion starts with α (π₁) = p (π₁| x₁) = p (x₁| π₁) p (π₁) /p (x₁). This shows that the filtered distribution α(.) is propagated through to the next time-step, where it acts like a ‘prior’. In other words, at each time step the calculated posterior becomes the new prior for the next time-step [26].

To estimate the posterior distribution of latent states, we sample from the joint distribution of the latent sequence, p(π_1:τ | x_1:τ):

According to Equation 8, to be able to compute the posterior distribution of latent states, we need ‘time-reversed’ transitions p(π_t |π_t+1, x_1:t), which are obtained using α-recursions from the forward filtering. The simulation starts from the end of the sequence and proceeds backwards recursively. Starting point is

For each time step t, the unnormalized sampling probability can be calculated using the previously sampled latent state , and the can be sampled using following:

This procedure is known as forward-filtering backward-sampling[26].

Learning Model Parameters

Conditionally on the latent states and observations the parameters θ = {π₀, T, E} of the model can be estimated. The posterior distribution is given by

Factorizing the likelihood term in Equation 14 as in Equation 1, we can get the following result:

Since, π₀, rows of the T, and rows of the E are probability distribution, [.]’s inside the Equation 16 corresponds to Multinomial distribution with parameters ,and where corresponds to the initial state counts, T^counts corresponds to the state to state transition counts, and E^counts corresponds to the state to observation emission counts. Since the Dirichlet is conjugate prior for the Multinomial distribution, we set a Dirichlet prior for θ = {π₀, T, E} with hyperparameters , and , respectively. Because of the conjugacy, the resulting posterior in Equation 14 is also a Dirichlet distribution:

Inference of HMM with MCMC Sampling

There are several ways to estimate the parameters of HMM when both hidden states and model parameters θ are unknown, for example the expectation-maximization algorithm, variational Bayes, or MCMC sampling. In this paper, we will use MCMC. The goal of MCMC is to produce draws from the posterior p(π, θ | x). The Gibbs sampler is an MCMC algorithm suitable for high dimensional problems, such as the HMM, and it samples parameters one-by-one using their distributions conditional on the other parameters. In the HMM, it will alternate between sampling the model parameters θ conditional on the latent sequence π from p(θ | π, x), and the latent sequence π conditional on the model parameters θ from p(π | θ, x). The algorithm iterates the following steps N times, which gives N draws from the posterior:

1. Sample the model parameters from Equation 17, Equation 18, and Equation 19, respectively. (Note that they are independent of each other given latent sequence.)

2. Sample the latent sequence of HMM using the updated model parameters θ using the Forward-Filtering Backward-Sampling algorithm.

Before calculating the posteriors, the first N₀ samples are discarded as a convergence (warm-up) period of the Markov chain. The plate diagram of the HMM is given in the Figure 5a, and the algorithm is described in Algorithm 1.

• Download figure
• Open in new tab

Figure 5 Plate diagram of (a) Hidden Markov Model (HMM) and (b) Coupled Hidden Markov Model (CHMM).

Here, for clarity, we illustrate only one chain. π_t represents the unobserved true states, and x_t represents the noisy observations at time t. The faded observation nodes correspond to missing values. Plain parameters are hyperparameters of the model. Edges from π_t to x_t represent emission probabilities. In the CHMM, the next state in one site depends on the previous states of the other sites. It is illustrated by π^’. Therefore, transition probabilities change at each time step.

Algorithm 1

Hidden Markov Model

• Download figure
• Open in new tab

Coupled Hidden Markov Model

In the ordinary HMM, transition parameters do not change. However, in the CHMM, transitions in one chain are affected by other chains. In principle, it would be possible to define a single joint HMM where the latent state represents the latent states of all individual chains jointly and adapt the solution for the HMM described above. However, this is inefficient where there are many chains because the number of states would grow as O(K^C), where K is the number of hidden states in a chain and C is the number of chains. On the high level, our key idea is that the transition matrix of each chain is modeled conditionally on the states of the other chains (and not as a single large joint transition matrix). A similar formulation is considered by [15]; however, the details are different. Note that a transition matrix T for a specific chain will not be time-independent anymore but changes at each time step depending on the states of the other chains. The graphical representation of the model (showing just two chains) is given in Figure 1, and the plate diagram of the CHMM is given in Figure 5b. There are O(KC) parameters in this formulation, which is much more efficient for the increased number of chains.

Defining the transition matrix is a critical part of the algorithm, and it is essential to carry as much information about the other chains as possible. We model the dependencies between the chains with parameters β. Assume there are C chains and let ℂ denote the set of chains. Further, assume that each chain can be in one of 𝕂 possible states, and let K denote the set of states. Finally, assume that there is a baseline state such that if a chain is in that state, it does not affect other chains (in our application, this state corresponds to the absence of MRSA colonization in the respective body site). The transition matrix for the chain ĉ, at time t, denoted as , is defined by;

The transition matrix is obtained in Equation 21 by applying a row-wise softmax operator σ_row to the unnormalized transition matrix . Parameter corresponds to an intercept matrix and it specifies the transition matrix of the target chain ĉ when all other chains are in the baseline state . Parameter represents the impact of chain c on the target chain ĉ and it is added to whenever chain c was in state in the previous time step. We will denote all the parameters for target chain ĉ as β^[ĉ].

In Equation 20, the unnormalized transition probability U is an additive function of the latent states of the other chains; hence we call it the Additive-CHMM. We also implemented a simple formulation than the Additive-CHMM, which we call Or-CHMM. In the Or-CHMM, the target chain is not affected by the other chains individually; instead, there the target chain is affected whenever any of the other chains is in some state different from the baseline state, and the unnormalized transition matrix is given by:

In our application we have K = 2, corresponding to the presence and absence of colonization in a given body site. In this case, Equation 22 has a simple form: If all other chains are in state , the output is only , and if any of the other chains is in a state other , then the output is the sum of and .

Design and Interpretation of the β parameters

In the CHMM, the transition matrix T for each chain is modeled as a function of β parameters, where each β_k is a matrix of the same size as the transition matrix. Therefore, the β parameters are a list of matrices. Since the rows of a transition matrix T are assumed independent, we also model β parameters such that the rows of those matrices are independent. However, to avoid redundant parameters, the rows on each row of β are assumed to sum to zero. In practice, we sample the first K −1 parameter on a given row and set the Kth element to equal the negative of the sum of all the other parameters. This constraint ensures that the mapping from β to T is one-to-one.

To give an interpretation to the parameters β, we start with Equation 21 and write the softmax for a single element of a transition matrix T (i, j): from which it follows after straightforward algebra:

In our application K = 2 and sums of the rows of β_k are assumed equal to 0, so consequently also the rows in the U_t sum to zero, which gives us: where ¬j refers to the other element on the row that is not j. Therefore, in this 2-dimensional case, all β₀ values correspond to half of the log-odds of the respective transition probability, and similarly, parameters β_k representing the interactions between the chains correspond to half of the change in the log-odds because of the presence of colonization in the other chain.

Adaptive Metropolis-Hastings Within Gibbs Sampling

So far, we have described the design specifications of the CHMM. Since transitions are changing at each time step, it is not feasible to calculate their sufficient statistics as in Equation 16. To resolve this, we sample the β^[ĉ]* using a Metropolis-Hastings(MH) step within the Gibbs sampler because conditional on the latent states of all the chains, the transition parameters are fully determined by the β^[ĉ]. Parameters. To get a draw from , where states of the other chains, we use the following steps;

1. Make a proposal β^[ĉ]* ~ q β^[ĉ]

2. Transform β^[ĉ]* samples to transition probabilities according to latent states of the other chains, π′

3. Accept the proposal β^[ĉ]* with probability

where

The first quantity on the right-hand side of the Equation 28 can be directly calculated with the transition parameters for corresponding chain. We used the Gaussian proposal in which the previously sampled β^[ĉ] is the mean such that q (β^[ĉ]* | β^[ĉ])= N(β^[ĉ]* | β β^[ĉ], c² Σ). As a variance parameter, we used a fixed diagonal covariance matrix during the warm-up period and then used the covariance matrix estimated from the previous samples. We initially set as metropolis jumping scaling factor since it is theoretically the most efficient scaling factor[27], where d is the dimension of the sampling. Then we adaptively scaled the scaling factor c as described in the [28]. Since this proposal distribution is symmetric (i.e., Normal), q(β^[ĉ] |44β^[ĉ]* ) and q(β^[ĉ]* |β^[ĉ])in Equation 27, cancel out each other.

To sample the latent sequences we need to modify the Forward-Filtering Backward-Sampling algorithm[18], since we want samples from instead of . Because in our application the within chain dependencies are stronger than between chain dependencies, we will use the Forward-Filtering Backward-Sampling with a modified model where the outgoing edges from the current chain are neglected as the proposal distribution q, and adjust for the discrepancy between the proposal and the true conditional distribution by using the following MH acceptance probability[18]: where p(.) distribution on the right-hand side of the Equation 29 is the joint posterior probability of the latent sequences, which is straightforward to calculate conditionally on the current CHMM parameters, and q(.) is the probability of samples for particular chain that can be calculated during the backward sampling process (for details, see [18]). Algorithm 2 summarizes the whole algorithm. To validate the implementation, we used the algorithm to estimate parameters in simulated data. The results in Supplementary Figure 1 show that the algorithm estimated the parameters correctly and yielded well-calibrated posterior distributions.

Implementation details

We set the initial covariance matrix for the metropolis proposal as 0.01 ×I, where I is the identity matrix, which corresponds to a step size giving the optimal acceptance rate of ≈23% [29]. We set the prior of β₀ as N (β₀ | 0, 5), which is almost uninformative so that the estimates are not affected strongly by the prior. For the rest of the β parameters, denoted by β_k, we used sparsity encouraging Horseshoe prior with mean and scale parameters are 0 and 1, respectively (see Suplementary Materials for details). We used a uniform prior on the initial state probabilities π₀, and very weak Dirichlet priors for the rows of emission probabilities E such that we set the value to 0.9 if the latent state emits the same observation and set the value to 0.1/(n −1), n is the number of states, if the latent state emits a different observation. In such a formulation, except for the initialization, the prior has a negligible effect because it is summed with observation counts during inference. We drew 15,000 MCMC samples, and we set the warm-up length as 7,500. Posterior probabilities are calculated using the remaining MCMC samples. HMM and CHMM implicitly assume that time intervals between observations are the same, which is not the case in our data. Therefore, during the model training, we assumed that there are missing observations at 2, 4, and 5 months after enrollment.

Algorithm 2

Coupled Hidden Markov Model

• Download figure
• Open in new tab