A stochastic model for MRSA transmission within a hospital ward incorporating environmental contamination

Summary

Methicillin-resistant Staphylococcus aureus (MRSA) transmission in hospital wards is associated with adverse outcomes for patients and increased costs for hospitals. The transmission process is inherently stochastic and the randomness emphasized by the small population sizes involved. As such, a stochastic model was proposed to describe the MRSA transmission process, taking into account the related contribution and modelling of the associated microbiological environmental contamination. The model was used to evaluate the performance of five common interventions and their combinations on six potential outcome measures of interest under two hypothetical disease burden settings. The model showed that the optimal intervention combination varied depending on the outcome measure and burden setting. In particular, it was found that certain outcomes only required a small subset of targeted interventions to control the outcome measure, while other outcomes still reported reduction in the outcome distribution with up to all five interventions included. This study describes a new stochastic model for MRSA transmission within a ward and highlights the use of the generalized Mann–Whitney statistic to compare the distribution of the outcome measures under different intervention combinations to assist in planning future interventions in hospital wards under different potential outcome measures and disease burden.

INTRODUCTION

Healthcare associated infections (HAIs) are adverse events that can arise during hospitalization. Multidrug-resistant organisms (MDROs), for example methicillin-resistant Staphylococcus aureus (MRSA), are common causes of these HAIs with patients typically becoming colonized with the organism prior to developing an infection. Treatment options for MDROs are becoming increasingly limited due to the relative scarcity in development of new treatments compared to the rate of resistance acquisition [1]. As such, the role of routine infection control and prevention (ICP) practices are of great importance in reducing the occurrence of HAIs.

Intervention studies which typically investigate the effects of one or a combination of interventions in reducing HAIs provide a good first line of evidence for particular interventions to be incorporated into routine ICP practices. These studies also assist in building mathematical model representations of the healthcare setting. Such models then allow for further probing of the effects of the interventions which may not have been feasible or potentially ethical to investigate in a clinical setting but could prove useful in assisting decision-making, particularly when hospital resources are severely limited. The model findings could also provide recommendations for future intervention studies.

Susceptible patients are typically modelled to be colonized (a state which precedes an infection) through a forcing term (referred to as the force of infection) which is a function of the number of colonized patients currently present in the ward as well as the colonized hospital staff in the ward at the time and also contact frequency. As hospital staff are not routinely screened for pathogen colonization [2], obtaining high-quality data on hospital staff has proven difficult.

That said, the most mathematical models consider vector-based cross-transmission between patients and transiently contaminated healthcare workers (HCWs) to be the dominant transmission mechanism for MDROs such as MRSA [3]. Only a small number of papers have considered alternative transmission routes typically by incorporating a constant source (such as in [4]). Even fewer have explicitly modelled environmental contamination as an alternative transmission route [5–10]. However, such models only calibrated the parameter estimates related to the environmental contamination to match observed patient incidence rather than using environmental contamination data.

This paper presents a stochastic model for ward MDRO transmission based on patient dynamics, as patient data are typically more readily available compared to hospital staff, coupled with a time-series model of environmental contamination which was parametrized by environmental contamination data. Due to the low reported prevalence of HCW carriage [11], the small proportion of nosocomial outbreaks attributable to HCWs [2] and the few adverse outcomes reported for HCWs [11], we assumed that transmission is implicitly facilitated by HCWs, who are temporarily contaminated with MRSA through contact with an MRSA-positive patient or environmental contamination, due to the limited mobility of patients, as is also common practice in similar modelling studies [4, 10, 12, 13]. Inclusion of HCWs typically involves substantial simplification of realistic HCW dynamics [8, 14] or substantial additional data collection to account for the heterogeneity between HCWs [15–19] beyond the scope of this study.

The model was run under two settings; the first is based on MRSA dynamics in a developed country (UK [12] and Switzerland [36] study estimates were used here) where MRSA data and parameters are more easily readily sourced, and the second is for a hypothetical scenario where the pathogen is more readily transmitted and not as easy to detect. The second setting could be representative of a novel pathogen in the healthcare setting, a new strain of MRSA that is more virulent than existing strains or perhaps reflective of a resource-poor setting such as in low-income countries [13] where such modelling studies could be of great benefit. The impact of five common healthcare interventions [3] and their various combinations were investigated for six potential outcome measures under both settings separately. Limitations and future directions in model development are provided in the Discussion.

METHODS

Model formulation

The model proposed is for a single ward setting and comprises of: (i) a ward-level patient arrival process; (ii) an individual-based model for patient transitions in the ward; and (iii) a time-series model for the level of environmental contamination.

At any time t, patients in the ward are categorized based on their MRSA status where they can be in the susceptible group [S(t)], the undetected MRSA colonized group [C _xd (t)], the detected with MRSA colonization and undergoing appropriate treatment group [C _d (t)], the undetected MRSA infected group [I _xd (t)], or the detected with MRSA infection and undergoing appropriate treatment group [I _d (t)]. A schematic illustration of the model is provided in Figure 1 with E(t) representing the ward environmental contamination levels.

Fig. 1. Compartmental diagram for the MRSA transmission model incorporating environmental contamination. The solid black lines represent patient transitions between the different states as well as admissions and discharges [only for the S(t) and C _xd (t) compartments]. The red dashed lines denote the contribution from the various compartments to the colonization process while the black dashed lines show the compartments contributing to the evolution of the E(t) compartment.

The model is an example of Discrete Event Simulation (DES), a technique that is widely used in healthcare research [20–22]. While perhaps more commonly used in scheduling problems, DES has also been applied to investigate pathogen transmission [21]. DES provides a flexible modelling approach to represent individual patient transitions during their hospitalization episode, allowing for the inclusion of stochastic variability (important for small population studies such as in a hospital ward) and effects of individual patient information.

Patient admissions into the ward are modelled as a right-censored (at ward capacity M) Poisson process [ $\!A\left(t \right)\sim Po\left( \lambda \right)]$ with a binomial variable to separate arrivals to either susceptibles [AS(t)] or colonized (but not detected, i.e.C _xd ) [AC(t)]. It is assumed that patients cannot be infected on admission (as infected patients are typically isolated or cohorted to reduce transmission risk to other patients). Excess arrivals, beyond the ward capacity M, are assumed to be allocated to a separate ward thus creating the right-censoring in the arrival process.

The likelihood for the admissions at time t can therefore be written as:

$$\eqalign{ P{\rm (}A\left( t \right) & = i,\; AS\left( t \right) = j,\; AC\left( t \right) = i - j\; {\rm \vert} Y\left( {t - 1} \right)) = \; \cr&{\hskip-32pt}\left\{ {\matrix{ {\displaystyle{{\lambda ^i} \over {i!}}{\rm exp}\left\{ { - \lambda} \right\}\left( {\matrix{ i \cr j \cr}} \right)\vartheta ^j \left( {1 - \vartheta} \right)^{i - j}} & {0 \le i \lt Y\left( {t - 1} \right)} \cr {\mathop \sum \limits_{l = Y\left( {t - 1} \right)}^\infty \displaystyle{{\lambda ^l} \over {l!}}{\rm exp}\left\{ { - \lambda} \right\}\left( {\matrix{ l \cr j \cr}} \right)\vartheta ^j \left( {1 - \vartheta} \right)^{l - j}} & {i = Y\left( {t - 1} \right),} \cr}} \right.}$$

where Y(t) is the number of empty beds in the ward at time t and $\vartheta $ is the proportion of admissions that arrive susceptible.

The admissions at time t will then be assigned to the empty beds in the ward but will not undergo the individual patient transitions until the next time point.

The individual-based model, which is for patient transitions in the ward, processes each patient present in the ward at each time point based on the patient's current MRSA status. The following assumptions were used to formulate the individual-based model patient transitions:

1. (1) Each patient can only undergo one transition (discharge, colonization, infection, recovery, detection) per time period.

2. (2) Susceptible patients have to be colonized before developing an infection.

3. (3) Patient colonization will always be undetected when first colonized.

4. (4) Colonized patients will not return to the susceptible state.

5. (5) Undetected colonized patients cannot transition directly to the detected infected state as it counts as two transitions (detection and infection).

6. (6) Detected colonized and infected patients cannot return to the undetected state.

7. (7) Detected colonized patients are placed under the decolonisation treatment and cannot develop an infection.

8. (8) Infected patients only recover to the colonized state, and not to the susceptible state.

9. (9) Detected infected patients are placed under an appropriate treatment which increases their probability of recovery over their infection duration.

10. (10) Undetected infected patients cannot recover as they have not received appropriate treatment yet.

At each time point t, each susceptible patient S can either leave the ward as susceptible with probability p _L , become colonized (but not detected) with probability p _C , or remain susceptible with probability p _s such that p _L  + p _C  + p _s  = 1.

The probability of being colonized is modelled as p _C  = f _E (1 − p _L ), where f _E is an increasing function of E(t), C _xd (t − 1), C _d (t − 1), I _xd (t − 1) and I _d (t − 1). Specifically, the following form for f _E was used

$$f_E \left( t \right) = 1 - \exp \left\{ { - \nu \left( t \right){\rm \Delta} t} \right\},\; $$

where ν(t) = β ₀ + β ₁ C _xd (t − 1) + β ₂ C _d (t − 1) + β ₃ I _xd (t − 1) + β ₄ I _d (t − 1) + β ₅ E(t) is the instantaneous hazard of being colonized, or also known as the force of infection for this model, and $0 \le f_E \left( t \right) \lt 1\; \forall \; t$ . Last, p _S  = (1 − f _E )(1 − p _L ).

Each undetected colonized patient C _xd is detected with probability ρ (assumed to be the screening test sensitivity). Otherwise, the undetected colonized patient can either leave the ward with probability q _L , develop an infection with probability q _I , or remain colonized in the ward with probability q _C such that q _L  + q _I  + q _C  = 1. No additional structure is imposed on these probabilities values as it is assumed that each colonized patients will have the same probability values.

Each detected colonized patient C _d can either leave the ward with probability q _L or remain colonized and detected with probability 1 − q _L . Due to a lack of information to differentiate the probability of leaving for undetected and detected colonized patients, these were assumed to be same. One of the interventions considered (DECOL) increases the probability of leaving for just the detected colonized patients.

Each undetected infected patient I _xd can either be detected with probability ρ or remain undetected with probability 1 − ρ.

Each detected infected patient I _d will have a probability r _C of recovering (transitioning to C _d ) where

$$r_C \left( {t{\rm \vert} \psi, \; ti_k} \right) = \; 1 - \exp \left\{ { - \psi \left( {t - ti_k} \right)} \right\},$$

is an increasing function of the difference of the current time (t) and the time the individual k became infected (ti _k ). In other words, it is assumed that the longer a patient is infected, the more likely the patient will recover at the next time point. An infected patient remains infected with probability 1 − r _C .

By definition, only the (approximate) date that a patient is detected to be colonized or infected is available from hospital surveillance databases. The transition times from susceptible to undetected with MRSA colonization (tc _k ), and subsequently undetected infection (ti _k ) are typically imputed from a range of plausible values between the patient's admission date (a _k ) and first positive screening test date (d _k ) where the full conditional for (tc _k , ti _k ) can be written as

$$\eqalign{& \left( {1 - \rho} \right)^{N_F \left( {ti_k} \right)} \cr & \quad {\rm exp}\left\{ {\mathop \sum \limits_b \log \nu \left( {t_b} \right) - \; \mathop \sum \limits_d S\left( {t_d} \right)\nu \left( {t_d} \right)\left( {t_{d + 1} - t_d} \right)} \right\} \cr & \quad\times q_I {\rm exp}\left\{ { - \; q_I \left( {ti_k - tc_k} \right)} \right\},} $$

where tc _k  < ti _k , N _F (ti _k ) is the number of false-negative screening test results for patient k given ti _k , the b subscript indexes time points where a susceptible patient becomes colonized between tc _k to patient k’s discharge and the d subscript indexes the time points where v(t) changes between a _k and tc _k . The expression can be evaluated for all potential (tc _k , ti _k ) values to obtain a discrete distribution to be used in a Metropolis–Hastings step within a Markov chain Monte Carlo algorithm to impute these unobserved quantities and estimate the remaining model parameters [4, 14, 23].

An autoregressive-moving average time series model with exogenous covariates (ARMAX) [24] is used to describe the environmental contamination levels E(t). The exogenous covariates assumed to be contributing to the levels of environmental contamination at time t are the C _xd and I _xd patients in the ward at time t − 1. It is assumed that detected (colonized and infected) MRSA patients undergo the decolonization treatment which halts shedding from the patient to the environment. The orders of the ARMAX model are determined using the auto.arima() function in the R package forecast [25].

Parameter values

The model parameter values used for the normal burden setting simulations are summarized in Table 1. Additional details of the parametrization are provided in the Supplementary material. The normal burden setting is reflective of MRSA burden in a typical hospital ward in a developed country. These parameters values are also used in the high burden setting simulations with the following modifications:

1. (1) there is an additional factor of two multiplying ν(t);

2. (2) the probability of a colonized patient developing an infection q _I is doubled and q _C is reduced accordingly to ensure q _L  + q _I  + q _c  = 1;

3. (3) there is decreased sensitivity in the screening test, ρ = 0·6;

that is, we assumed that in this setting, the hypothetical pathogen is more likely to colonize susceptible patients, colonized patients more readily develop an infection and it is harder to detect the presence of the pathogen. The high burden setting attempts to mimic either the MRSA dynamics in a developing country [26] or a novel strain of pathogen that is more virulent and less readily detected by routine surveillance.

Table 1. Parameter values for the stochastic model describing multidrug-resistant organisms' transmission in a hospital ward

AR, Autoregressive; MA, moving average.

* Unpublished observations are estimates obtained from fitting a non-homogeneous Poisson process to the data. More details provided in the supplementary material.

There was no available source to estimate the parameter ω which represents the difference between colonized and infected patients on the force of infection. The ω value in the Results section was 1 as a reflection of the lack of information on the parameter. Alternative values of 0·1 and 1·9 were also investigated in the parameter sensitivity analysis (provided in the Supplementary material). We found that the AR, C _xd and C _d outcomes (defined below) were particularly sensitive to a low value of ω (giving a stronger influence to colonized patients) in both normal and high burden settings. Distributions of AR outcome for the different values of ω are provided in Figure 2. Similar plots for the other outcomes and parameters are provided in the Supplementary material.

Fig. 2. Attack rate outcome for normal burden (left plot) and high burden (right plot) settings. The x-axis denotes the baseline, low ω value and high ω value (moving from left to right).

Interventions

Five common intervention strategies were considered in the model investigation below:

1. (1) Not colonized on admission (COA) ( $\vartheta = 1$ ), where all patients who are colonized on admission are assumed to be detected on admission and isolated elsewhere, i.e. universal screening [27].

2. (2) Improved environmental cleaning (ENV), which halved the intercept term in the environmental time-series model (α ₁) [28].

3. (3) Improved contact precaution practices (CP), which decreases ν by a factor of ξ where ξ was set to 0·75 [29].

4. (4) Perfect screening test sensitivity (SENS), where test sensitivity ρ was set to 1 [14].

5. (5) Improved decolonization treatment for colonized patients (DECOL) where the probability for a C _d patient leaving the ward is now q _L  + Δ (with the probability of staying adjusted accordingly) [14].

We considered six outcome measures for the investigations. They are the attack rate (AR) defined as the average of the force of infection ν(t) [14] as well as the cumulative numbers of

• • patients who were colonized on admission (AC),

• • patients who were colonized but not detected (C _xd ),

• • detected, colonized patients (C _d ),

• • patients who were infected but not detected (I _xd ),

• • detected, infected patients (I _d ).

Note that there is a slight abuse of notation where C _xd , C _d , I _xd and I _d refer to the cumulative number of patients in each group for the outcome measures, but the time-varying prevalence of the groups in the model.

Due to the stochastic model formulation, each intervention setting was simulated 1000 times and we compared the distributional differences of the outcomes rather than just point estimates of the outcomes.

Pairs of distributions (denoted generally by X and Y here) were assessed using the generalized Mann–Whitney statistic which estimates the parameter

$$\theta = P\left( {Y \gt X} \right) + \; \textstyle{1 \over 2}\; P\left( {Y = X} \right)\,{\rm using}\,\hat \theta = \; \displaystyle{U \over {mn}},$$

where $U = \; \mathop \sum \nolimits_{i = 1}^m \mathop \sum \nolimits_{j = 1}^n {\open 1}\left( {Y_j \gt X_i} \right) + \; \displaystyle{1 \over 2}{\rm \;} {\open 1}\left( {Y_j = X_i} \right)$ with {Y _j ; j = 1,…, n} and {X _i ; i = 1,…, m} being samples from the Y and X distributions, respectively. Confidence intervals (CIs) for $\hat \theta $ were computed based on method 5 of Newcombe [30].

Following the definition above, values of θ larger than 0·5 indicate that the Y is stochastically larger than X and, conversely, values of θ less than 0·5 indicate X is stochastically larger than Y. For the results below, θ values between 0 and 0·2 (and similarly between 0·8 and 1) are considered strong evidence that the two distributions are substantially different. Intermediate θ values between 0·2 and 0·4 (or 0·6 and 0·8) are assumed to provide weak evidence of a difference between the distributions. Values of θ close to 0·5 (between 0·4 and 0·6) indicate that there is no evidence that the two distributions being compared are dissimilar.

RESULTS

The results for the normal burden setting (Table 2 and Table 3) and high burden setting (Table 4 and Table 5) are summarized below. More detailed comparisons of the interventions combinations for all outcome measures using the generalized Mann–Whitney statistic are provided in the Supplementary material.

The results for the AC, I _xd and I _d outcomes were similar for both the normal and high burden settings, and discussed together here. Results for the AR, C _xd and C _d outcomes are discussed separately for the normal burden setting and high burden setting.

The most important intervention for the AC outcome was the COA intervention which eliminates the possibility of colonized patients being admitted. As such, the COA intervention (and any other intervention combinations which include COA) greatly outperforms interventions of any size which do not include the COA intervention in both settings. Any intervention combination which includes the COA intervention achieved 0 AC, whereas intervention combinations without the COA intervention produced AC distributions with 95% intervals that do not include 0.

The performance of the interventions on the I _d outcome was very similar to that for I _xd since the only transition to I _d is through I _xd , i.e. eliminating I _xd would also eliminate the I _d population. As such, only the results for the I _xd results are discussed for brevity as identical inferences apply to the I _xd outcome. The SENS intervention was the most important intervention for the I _xd outcome as having perfect sensitivity would allow detection of all colonized patients prior to infection developing. As such, the best performing intervention of any size will include the SENS intervention.

However, it should also be noted that the I _xd outcome is generally small for the normal burden setting with even the baseline I _xd having a 95% CI of 0–2 (Table 2).

Table 2. Numerical summaries of output measures for normal burden setting

For explanation of abbreviations in column 1 see ‘Interventions’ section in main text. Values in parentheses are 95% confidence intervals.

Table 3. Summary of intervention combination comparisons for the normal burden setting

CI, Confidence interval.

For explanation of abbreviations in Comparison column see ‘Interventions’ section in main text.

In contrast with the normal burden setting, the SENS intervention (or any combination which includes the SENS intervention) was substantially more favourable in the high burden setting (Table 5). The SENS intervention substantially outperformed all intervention combinations which excluded the SENS intervention here.

Normal burden setting

Table 2 provides the numerical summary of the six outcome measures under the baseline and the various combinations of the five interventions investigated and Table 3 lists the θ comparisons for optimal interventions of different sizes. The baseline scenario refers to the case without any interventions.

There were great improvements in reducing the AR outcome when increasing the number of interventions by up to three with the optimal triplet being {COA, ENV, CP} [2·66 (2·20−3·31) × 10⁻³] (values in parentheses are 95% confidence intervals). This triplet outperformed the best single intervention [CP with AR of 4·32 (3·69−5·05) × 10⁻³] and intervention pair [{COA, CP} with AR of 3·35 (2·88−4·01) × 10⁻³]. The addition of one extra intervention (either DECOL or SENS) did not seem to have a marked effect on the AR distribution [2·50 (2·13−3·02) × 10⁻³ and 2·53 (2·19−2·92) × 10⁻³, respectively]. However, there is a benefit in implementing all five interventions [AR = 2·39 (2·11−2·71) × 10⁻³] compared to just the best three interventions.

For the C _xd outcome, the two best performing pairs [{ENV, CP} and {COA, CP} with C _xd of 17·59 (10–27) and 17·60 (9–28), respectively] performed slightly better compared with the best single intervention [CP with C _xd of 20·78 (12–31)]. A similar performance gain was noted when comparing the best intervention triplet [{COA, ENV, CP} with C _xd  = 14·29 (6–24)] to both the best performing pairs. There does not appear to be substantial changes in the C _xd difference when comparing across the best performing triplet, quartets [{COA, ENV, CP, SENS} and {COA, ENV, CP, DECOL} with C _xd of 13·65 (6–23) and 13·94 (6–23), respectively] and the combination of all interventions [13·44 (6–22)], indicating that there is little gain from considering anything beyond the best performing triplet in reducing the distributional outcome of C _xd for this scenario.

Comparing across different intervention sizes for the C _d outcome, there are notable reductions in support for considering additional numbers of interventions up to the best performing intervention triplet [{COA, ENV, CP} with C _d of 13·96 (6–24)]. The best performing single intervention for the C _d outcome was COA [24·22 (14–36)] and the best performing intervention pair was {COA, CP} (17·21 (9·5–27). There are no discernible differences in the C _d outcome distributions in implementing all five interventions [C _d  = 13·43 (6–22)] or either of the two best performing quartets identified [{COA, ENV, CP, DECOL} and {COA, ENV, CP, SENS} with C _d of 13·32 (6–22) and 13·95 (6–23), respectively] compared to having just the best performing intervention triplet (with θ estimates ranging from 0·46 to 0·50).

High burden setting

The mean and 95% CIs for the six outcome measures across the different intervention combinations considered are listed in Table 4. Compared to the baseline scenario in the normal burden setting (Table 2), we see notable increases in the average AR, C _xd , C _d , I _xd and I _d outcomes but a slight reduction in the AC outcome likely due to the decreased number of admissions overall as colonized and infected patients stay in the ward longer. The comparisons across optimal interventions of different sizes are provided in Table 5 for the high burden setting.

Table 4. Numerical summaries of output measures for high burden setting

For explanation of abbreviations in column 1 see ‘Interventions’ section in main text.Values in parentheses are 95% confidence intervals.

Table 5. Summary of intervention combination comparisons for the normal burden setting

CI, Confidence interval.

For explanation of abbreviations in Comparison column see ‘Interventions’ section in main text.

For the AR outcome in the high burden setting, there is evidence to consider implementing the maximum number of interventions possible (subject to resource constraint) beginning with the CP intervention [12·44 (10·14 − 14·83) × 10⁻³], followed by the SENS intervention [{CP, SENS} with AR of 9·50 (8·35 − 10·79) × 10⁻³], either the COA or ENV intervention [{COA, CP, SENS} with AR of 7·88 (6·77 − 9·14) × 10⁻³ or {ENV, CP, SENS} with AR 7·97 (6·71 − 9·24) × 10⁻³] or both [{COA, ENV, CP, SENS} with AR 6·25 (5·10, 7·53) × 10⁻³], up to all five interventions [5·55 (4·73, 6·46) × 10⁻³]. The reduction in the AR distribution when moving from the best performing quartet to all intervention was not as marked as the other increases in intervention sizes.

Only small gains were obtained from increasing the size of the intervention combinations sequentially for the C _xd outcome. More notable reductions were obtained by moving from the best performing single intervention [CP with C _xd of 45·46 (30–61)] to at least one of the best performing triplets [{ENV, CP, SENS}, {COA, ENV, CP} or {COA, CP, SENS} with C _xd s of 36·57 (23–50), 37·24 (22–53) and 39·21 (26–55), respectively], and similarly from one of the best performing intervention pairs [{ENV, CP}, {CP, SENS} or {COA, CP} with C _xd s of 40·95 (28–55·5), 42·70 (29·5–58) and 43·56 (28–60), respectively] to either the {COA, ENV, CP, SENS} quartet [32·02 (19–46)] or all five interventions [29·95 (17–45)].

For the C _d outcome measure, the results obtained suggest it would be beneficial to consider up to the best performing triplet of interventions [{COA, ENV, CP} with C _d 33·85 (20–49)] subject to resource constraints. The best performing single interventions were COA [53·96 (39–72·5)] and CP [55·58 (39–74)], and the best performing intervention pair was {COA, CP} [39·72 (26–55)]. There was only a slight gain in moving from the best performing triplet to the combination of all interventions [29·95 (17–45)]. The two best performing intervention quartets [{COA, ENV, CP, SENS} and {COA, ENV, CP, DECOL} with C _d ’s of 32·02 (19–46) and 32·80 (19–49), respectively, did not yield C _d distributions substantially different from the best performing triplet.

DISCUSSION

The results obtained from the proposed stochastic model showed that there are differences in the optimal set of interventions depending on the outcome measure of interest as well as the burden setting of the pathogen (as summarized in Table 6).

Table 6. Overall order of importance for the five interventions considered under the normal and high burden setting

For explanation of abbreviations in last two columns see ‘Interventions’ section in main text.

↔ Denotes exchangeability in the order of the interventions and || denotes the optimal sized interventions, i.e. addition of interventions to the right of the || symbol would not affect the associated outcome measure.

For the AC outcome, I _xd and I _d outcome measures where one of the interventions considered eradicated the respective outcome measure (COA for the AC outcome and SENS for both I _xd and I _d ), only that particular intervention was required. This finding, particular for the I _xd and I _d outcome measures, may not be very realistic given that there is always some amount of delay between sample collection and the corresponding action based on the screening results. However, the θ performance measure still showed that in the normal burden setting, eradication of I _xd and I _d was only a slight improvement compared with the other intervention combinations and the baseline on the account of the already low baseline I _xd and I _d prevalence. This is not the case in the high burden setting where eradication of the I _xd and I _d outcomes with the SENS intervention was markedly different from the other intervention combinations which exclude SENS and the baseline scenario. The addition of the aforementioned small delay would have affected all scenarios considered equally and would unlikely have changed the finding in the normal burden setting. It is also unlikely to change the findings in the high burden setting unless the delay was substantive (of the order of days).

The model presented used parameter estimates combined from multiple sources. While it would be ideal if the model parameters were all obtained from one source, this is frequently not the case in such modelling studies where the hypothetical investigations considered typically require some form of data collation from multiple sources in order to fully parametrize the model [5–10]. It could also be argued that this provides such modelling studies with a level of flexibility that could not be obtained from clinical intervention studies. The lack of additional individual patient data for this study also precluded demonstration of the full utility of the individual-based patient transition component in the model. For this application, only the patient transition from I _d to C _d was based on their individual infection times (see expression for r _C ). However, the model can readily include individual-specific covariates into other transition probabilities in the model as well.

There are a number of extensions to the stochastic model proposed here that were not considered. Most of these extensions also involve additional data structures that are not readily available.

One such extension is to generalize the force of infection term such that the colonization threshold is no longer constant [23]. Under the current model formulation, the probability of a patient being colonized is only a function of the current force of infection. However, the generalization proposed in Streftaris & Gibson [23] allows for this transition to also depend on the accumulation of the force of infection terms from a patient's admission date to their colonization date. This quantity is known as the colonization threshold and requires prior knowledge or imputation of the colonization date in order to compute it. This extension is another approach to incorporate patient heterogeneity into the model, specifically related to patient susceptibility.

Another potential extension is to extend the one ward model to a multi-ward model using one of the meta-population models [31, 32] such as the multi-patch models (where each patch represents a ward) or more generally, temporal network models taking into account the fact that the edges between nodes change quite frequently with staff shift changes, and patient admissions and discharges, making the temporal element of the network more important [33, 34]. The high-frequency contact data required for such models have only recently started to be collected [35] and could prove to be a promising research avenue in providing a realistic, detailed representation of hospital pathogen transmission in a ward.

The inclusion of explicit representations of HCWs' roles in the pathogen transmission could be considered in extensions of the model presented here. While having explicit representation of HCWs allows for more realistic investigation of HCW-related interventions, this extension requires either incorporation of additional model assumptions on the HCWs' behaviours, or substantial additional data collection as HCWs are known to be highly heterogeneous population with different HCW categories (e.g. nurses, physicians, technicians) having differing patient contact rates, compliance levels to infection control and prevention practices, and work schedules [15–19]. Moreover, due to the low carriage rates in HCW reported [11], frequent screening of HCWs would be required in order to accurately quantify the temporary contamination status of HCWs, which is associated with high cost and staff time. It is also likely that this extension would require the aforementioned multi-ward extension to realistically capture the impact of HCWs in MRSA transmission as HCWs tend to work across multiple wards.