Modeling the effect of antibiotic exposure on the transmission of methicillin-resistant Staphylococcus aureus in hospitals with environmental contamination

1.   Introduction

Nosocomial infections caused by antibiotic-resistant bacteria are a major threat to global public health today. According to the Centers for Disease Control and Prevention (CDC) ^[1]:" Each year in the United States, at least 2 million people become infected with bacteria that are resistant to antibiotics, and at least 23,000 people die each year as a direct result of these infections." Furthermore, CDC identifies methicillin-resistant Staphylococcus aureus (MRSA), a Gram-positive bacterium, as one of the most common causes of hospital-acquired infections, especially in newborns and intensive-care unit patients. MRSA infections are usually treated by antibiotics, however, due to overprescribing and misprescribing, MRSA has been resistant to multiple commonly used antibiotics, which makes MRSA infections harder to be treated and even causes life-threatening cases in intensive care units.

An observation that patients with MRSA are about 64$\%$ more possible to die than patients with a non-resistant form of the infection in hospitals was revealed by a World Health Organization report in April 2014 ^[2]. In fact, some studies have observed a clear association between antibiotic exposure and MRSA isolation (Dancer ^[3], Tacconelli ^[4], Tacconelli et al. ^[5]). It has been shown that patients with prior antibiotic exposure are vulnerable to skin infections and are more likely to be colonized by MRSA. Furthermore, due to antibiotic exposure, patients may have a lower probability of successful treatment, a lengthier stay in hospitals and extra costs of treatment. Hence, it is necessary to take antibiotic exposure into account in modeling MRSA infections in hospitals. It is also found that under certain circumstances MRSA is capable of surviving for days, weeks or even months on environmental surfaces such as door handles, healthcare facilities, health-care workers' gloves (Boyce et al. ^[6], Dancer ^[7]). So environmental contamination is also a necessarily essential factor when we study the transmission of MRSA in hospitals.

In order to understand the diverse factors contributing to infections of antibiotic-resistant bacteria such as MRSA in hospitals, various models have been proposed and studied, see for example Austin and Anderson ^[8], Austin et al. ^[9], Bergstrom et al. ^[10], Bonhoeffer et al. ^[11], Bootsma et al. ^[12], Browne et al. ^[13], Browne and Webb ^[14], Chamchod and Ruan ^[15], Cooper et al. ^[16], D'Agata et al. ^[17,18], Hall et al. ^[19], Huang et al. ^[20,21], Lipsitch et al. ^[22], Smith et al. ^[23], Wang et al. ^[24], Wang and Ruan ^[25], Wang et al. ^[26,27], Webb ^[28], Webb et al. ^[29]. We refer to survey papers of Bonten et al. ^[30], Grundmann and Hellriegel ^[31], Temime et al. ^[32], van Kleef et al. ^[33] and the references cited therein on modeling antimicrobial resistance. These studies showed quantitatively how infection control measures such as hand washing, cohorting, and antibiotic restriction affect nosocomial cross-transmission. It is observed that the direct transmission via the hands of health-care workers (HCWs) is a crucial factor in the transmission of MRSA. In particular, D'Agata et al. ^[34] developed models to investigate the impact of persistent gastrointestinal colonization and antibiotic exposure on transmission dynamics of vancomycin resistant enterococci (VRE). Chamchod and Ruan ^[15] proposed models to investigate the effect of antibiotic exposure on the transmission of MRSA in hospitals. Mathematical models have also been developed to study the effect of environmental contamination on the spread of antibiotic-resistant bacteria in hospitals (Browne and Webb ^[14], McBryde and McElwain ^[35], Plipat et al. ^[36], Wang and Ruan ^[25], Wang et al. ^[26,27], Wolkewitz et al. ^[37]). Especially, Wang et al. ^[26,27] used both deterministic and stochastic models to focus on exploring the interaction between volunteers and their environment in the hospital system with a special care pattern in China. However, to the best of our knowledge, the combined effects of antibiotic exposure and environmental contamination have not been studied. This is the motivation for the current study.

In our previous studies (Wang et al. ^[24] and Wang and Ruan ^[25]) on nosocomial infections of MRSA in the emergency ward and respiratory intensive care unit in Beijing Tongren Hospital, Beijing, China, data on HCW, volunteers, patients, and environmental contamination were obtained. Wang et al. ^[24] constructed a mathematical model to determine the role of volunteers in the prevalence and persistence of MRSA in Beijing Tongren Hospital. Wang and Ruan ^[25] studied the effect of environmental contamination on the transmission MRSA in Beijing Tongren Hospital. Based on the data in ^[24,25], in this article, we first develop a deterministic ordinary differential equations model (ODE) to investigate the combined effects of antibiotic exposure and environmental contamination on the transmission dynamics of MRSA in hospitals. When there is no admission of colonized patients, we study the steady-states and estimate the basic reproduction number. Numerical simulations and sensitivity analysis are also provided. However, in hospital subunits, where the population is usually small, randomness may matter. Then we formulate a stochastic differential equations model (SDE) to study the transmission dynamics of MRSA that are not well described by the deterministic ODE model. Numerical simulations show that the average of multiple stochastic outputs aligns with the ODE output.

2.   Deterministic and stochastic models

2.1.   The deterministic model

The patients, healthcare workers (HCWs), and free-living bacteria in the environment in hospitals are divided into the following seven compartments (see Figure 1):

$P_u(t)$ – number of uncolonized patients without antibiotic exposure at time $t;$

$P_{uA}(t)$ – number of uncolonized patients with antibiotic exposure at time $t;$

$P_c(t)$ – number of colonized patients without antibiotic exposure at time $t;$

$P_{cA}(t)$ – number of colonized patients with antibiotic exposure at time $t;$

$H_u(t)$ – number of uncontaminated healthcare workers at time $t;$

$H_c(t)$ – number of contaminated healthcare workers at time $t;$

$B_e(t)$ – density of the free-living bacteria in the environment at time $t.$

We define antibiotic exposure as having received antibiotics within one month on admission or receiving any antibiotics treatment currently in hospital ^[38]. We also assume that there is no cross-infection between patients, the hospital is always fully occupied, and the bacteria in the environment are uniformly distributed. Patients are in four compartments according to their status: uncolonized without or with antibiotic exposure ($P_u$ and $P_{uA}$, respectively), colonized without or with antibiotic exposure ($P_c$ and $P_{cA}$, respectively). Patients are recruited at a total rate $\Omega(t)$ from any of these four compartments with the corresponding fractions $\theta_u$, $\theta_{uA}$, $\theta_c$, and $\theta_{cA}$, respectively, and can be discharged from any of these four compartments at corresponding rates $\gamma_u$, $\gamma_{uA}$, $\gamma_c$, and $\gamma_{cA}$, respectively. We calculate the discharge rate as the reciprocal of the length of stay specific for each compartment.

Based on the assumption that the hospital is always fully occupied, it is reasonable to see that the inflow of patients is equivalent to the outflow of patients, that is $\Omega(t) = \gamma_uP_u+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA}$, which results in a constant population size of patients $N_p = P_u+P_{uA}+P_c+P_{cA}$. Note that the total number of HCWs is also a constant $N_h = H_u+H_c$. Besides, patients without antibiotic exposure would move to patients with antibiotic exposure at an antibiotic prescribing rate of $\epsilon$ per day ^[39]. Since we assume that there is no cross-infection between patients, uncolonized patients without antibiotic exposure $P_u$ can be colonized either by direct contacting contaminated HCWs, $\alpha_p\beta_p(1-\eta)P_uH_c$, or indirect transmission via free-living bacteria in the environment, $\kappa_pP_uB_e$. A similar process occurs as uncolonized patients with antibiotic exposure move to colonized patients with antibiotic exposure, $\alpha_p\beta_{pA}(1-\eta)P_{uA}H_c+\kappa_{pA}P_{uA}B_e$. $\alpha_p$ is the contact rate per day per person, $\beta_p\ (\beta_{pA})$ is the probability of colonization after a contact with contaminated healthcare worker for $P_u\ (P_{uA})$, $\eta$ is the compliance rate with the hand hygiene, and $\kappa_p\ (\kappa_{pA})$ is the colonization rate from the environment for $P_u (P_{uA})$. Meanwhile, HCWs move from uncontaminated state to contaminated state either by contacting colonized patients (without or with antibiotic exposure), $\alpha_p\beta_h(1-\eta)P_cH_u+\alpha_p\beta_{hA}(1-\eta)P_{cA}H_u$, or by the indirect transmission via free-living bacteria in the contaminated environment, $\kappa_hH_uB_e$, where $\kappa_h$ is the contamination rate from the environment for HCWs. Because of hygiene standard of HCWs in hospitals, contaminated HCWs can move to uncontaminated HCWs, $\mu_c H_c$, in which $\mu_c$ is the decontamination rate for the HCWs. Even though free-living bacteria in the environment can survive for months, they cannot reproduce themselves in hospitals due to lack of enough reproduction conditions. Hence shedding bacteria into the environment from colonized patients, $\upsilon_pP_c+\upsilon_{pA}P_{cA}$, serves as an important source of transmission. As well contaminated HCWs contaminate the environment wherever they touch at the rate of $\upsilon_hH_c$. Here $\upsilon_p$, $\upsilon_{pA}$, and $\upsilon_h$ are the corresponding contaminated rates to the environment from $P_c$, $P_{cA}$, and $H_c$. $\gamma_b$ is the cleaning/disinfection rate of the environment.

Due to antibiotic exposure, patients may be more likely to have thrush, skin rashes, and gastrointestinal symptoms, and they have a higher probability of colonization, a lower probability of successful treatment, a lengthier stay, and a larger contamination rate to HCWs and the environment. We therefore assume that $\beta_{pA} \geq \beta_p$, it was estimated that uncolonized patients with antibiotic exposure $P_{uA}$ is 1.67 times more likely to be colonized than uncolonized patients without antibiotic exposure $P_u$ ^[39,15]. We also assume that $\beta_{hA} \ge \beta_h$ (a higher contamination rate to HCWs), $\gamma_{cA} \le \gamma_c \le \gamma_{uA}\le \gamma_u$ (a lengthier stay or a lower discharge rate), $\upsilon_{pA} \ge \upsilon_p$ (a higher contamination (shedding) rate to the environment). Besides, of new admission, the fraction of patients having antibiotic exposure ($\theta_{uA}+\theta_{cA}$) is assumed to be 0.38 ^[40,15]. As this article is a continuation of our previous studies on modeling the effect of antibiotic exposure ^[15] and impact of environmental contamination in Beijing Tongren Hospitals ^[24,25], we adapt most parameter values from these papers.

Detailed parameter values are listed in Table 1. From the flowchart shown in Figure 1, we formulate an ordinary differential equations model describing the transition between compartments as follows:

with initial conditions $P_u(0) = P_u^0$, $P_{uA}(0) = P_{uA}^0$, $P_c(0) = P_c^0$, $P_{cA}(0) = P_{cA}^0$, $H_u(0) = H_u^0$, $H_c(0) = H_c^0$, $B_e(0) = B_e^0$ specified at time 0.

2.2.   The stochastic model

We know that one disadvantage of deterministic models is that they cannot directly reflect randomness in epidemic events. For nosocomial models in hospital subunits where randomness may matter, there is a need to formulate randomness more precisely. The natural extensions of a deterministic ordinary differential equations model are usually two types of stochastic setting, continuous-time Markov chains (CTMCs) and stochastic differential equations (SDEs), where the time variable is continuous, but the state variables are discrete or continuous, respectively. In the formulation of CTMCs, forward Kolmogorov differential equations for the transition probability density functions can be derived and they, in turn, lead directly to the SDEs. Even though it is difficult to find analytical solutions for multivariate processes, SDEs are useful to numerically simulate stochastic realizations (sample paths) of the process. It is believed that the SDEs are easier to solve numerically than the Kolmogorov differential equations and faster than simulating sample paths of the CTMCs model for multivariate processes ^[42]. Thus, we develop a CTMCs model, an SDE model and its simulations in the following (^[43,25]).

By the assumption $P_u+P_{uA}+P_c+P_{cA} = N_p, H_u+H_c = N_h, \forall t\geq0$, we have the process $(P_c, P_{uA}, P_{cA}, H_c, B_e)$ in $\mathbb{R}^5$ with $P_u(t) = N_p-P_{uA}-P_c-P_{cA}$ and $H_u(t) = N_h-H_c$. These five variables have a joint probability denoted by

with $s \ge 0, j \ge 0, k \ge 0, 0 \le s+j+k \le N_p, 0 \le m \le N_h$ and $n \ge 0$. Assume that $\triangle t > 0$ is sufficiently small, the transition probabilities associated with the stochastic process are defined for a small period of time $\triangle t > 0$ as follows:

where $i_1, i_2, i_3, i_4, i_5\in \{ -1, 0, 1 \}$, Hence the transition probability is as follow:

We must choose the time step $\triangle t$ sufficiently small. In our case it is too complicated to express the transition matrix. Instead, we still are able to express the probabilities $p_{(s, j, k, m, n)}(t+\triangle t)$ by using the Markov property:

Naturally, a system of forward Kolmogorov differential equations can be derived:

We now try to develop a stochastic differential equations model (SDE) from the deterministic epidemic model (2.1). The system has five variables with a joint probability defined by:

with $s, j, k = 0, ..., N_p, m = 0...N_h$, and $n\geq 0$, with transition probabilities given in (2.2). Let $X(t) = (P_c(t), P_{uA}(t), P_{cA}(t), H_c(t), B_e(t))^T$ with infinitesimal

We express the infinitesimal mean matrix $f(X(t), t)$ as follows:

where

$e_c = \theta_c(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})+\alpha_p\beta_p(1-\eta)(N_h-P_c-P_{uA}-P_{cA})H_c+\kappa_pP_uB_e-\gamma_cP_c-\epsilon P_c$,

$e_{uA} = \theta_{uA}(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})-\alpha_p\beta_{pA}(1-\eta)P_{uA}H_c-\kappa_{pA}P_{uA}B_e-\gamma_{uA}P_{uA}\\ +\epsilon (N_h-P_c-P_{uA}-P_{cA})$,

$e_{cA} = \theta_{cA}(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})+\alpha_p\beta_{pA}(1-\eta)P_{uA}H_c+\kappa_{pA}P_{uA}Be-\gamma_{cA}P_{cA}+\epsilon P_c$,

$e_h = \alpha_p\beta_h(1-\eta)P_c(N_h-H_c)+\alpha_p\beta_{hA}(1-\eta)P_{cA}(N_h-H_c)+\kappa_h(N_h-H_c)B_e-\mu_cH_c$,

$e_b = \upsilon_pP_c+\upsilon_{pA}P_{cA}+\upsilon_hH_c-\gamma_bB_e$.

and also the infinitesimal variance matrix $\Sigma (X(t)t)$:

where

$\delta_c = \theta_c(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})+\alpha_p\beta_p(1-\eta)(N_h-P_c-P_{uA}-P_{cA})H_c+\kappa_pP_uB_e+\gamma_cP_c+\epsilon P_c,$

$\delta_{uA} = \theta_{uA}(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})+\alpha_p\beta_{pA}(1-\eta)P_{uA}H_c+\kappa_{pA}P_{uA}B_e+\gamma_{uA}P_{uA}\\ +\epsilon (N_h-P_c-P_{uA}-P_{cA}),$

$\delta_{cA} = \theta_{cA}(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})+\alpha_p\beta_{pA}(1-\eta)P_{uA}H_c+\kappa_{pA}P_{uA}B_e+\gamma_{cA}P_{cA}+\epsilon P_c,$

$\delta_h = \alpha_p\beta_h(1-\eta)P_c(N_h-H_c)+\alpha_p\beta_{hA}(1-\eta)P_{cA}(N_h-H_c)+\kappa_h(N_h-H_c)B_e+\mu_cH_c, \\ \delta_b = \upsilon_pP_c+\upsilon_{pA}P_{cA}+\upsilon_hH_c+\gamma_bB_e$.

It is easy to check that $\delta_c, \delta_{uA}, \delta_{cA}, \delta_h, \delta_b$ are all nonnegative. Hence we have a matrix $G$ satisfying $GG^T = \Sigma$, where $G$ is a 5$\times$ matrix to order $\triangle t$,

where

$a_1 = \theta_c(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})+\alpha_p\beta_p(1-\eta)(N_h-P_c-P_{uA}-P_{cA})H_c+\kappa_pP_uB_e$,

$a_2 = \gamma_cP_c$,

$a_3 = \epsilon P_c$,

$a_4 = \theta_{uA}(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA})+\epsilon (N_h-P_c-P_{uA}-P_{cA})$,

$a_5 = \kappa_{pA}P_{uA}Be+\gamma_{uA}P_{uA}$,

$a_6 = \alpha_p\beta_{pA}(1-\eta)P_{uA}H_c$,

$a_7 = \theta_{cA}(\gamma_u(N_h-P_c-P_{uA}-P_{cA})+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})+\kappa_{pA}P_{uA}B_e$,

$a_8 = \gamma_{cA}P_{cA}$,

$a_9 = \alpha_p\beta_h(1-\eta)P_c(N_h-H_c)+\alpha_p\beta_{hA}(1-\eta)P_{cA}(N_h-H_c)+\kappa_h(N_h-H_c)B_e$,

$a_{10} = \mu_cH_c$,

$a_{11} = \upsilon_pP_c+\upsilon_{pA}P_{cA}+\upsilon_hH_c$,

$a_{12} = \gamma_bB_e$.

Then the stochastic differential equations have the following form:

More precisely,

where $W_1, \cdots, W_{12}$ are twelve independent Wiener processes. Next we are able to run stochastic simulations.

3.   Mathematical analysis of the deterministic model

In this section we provide detailed analysis of the deterministic ODE model (2.1).

3.1.   Positivity and invariance of solutions

Based on the biological background of model (2.1), we only consider solutions of model (2.1) starting at $t = 0$ with nonnegative initial values:

Lemma 1. If $P_u^0, P_{uA}^0, P_c^0, P_{cA}^0, H_u^0, H_c^0, B_e^0 \ge 0$, then $(P_u(t), P_{uA}(t), P_c(t), P_{cA}(t), H_u(t), H_c(t), B_e(t))$ the solutions of model (2.1) are nonnegative for all $t \ge 0$ and ultimately bounded. In particular, if $P_u^0, P_{uA}^0, P_c^0, P_{cA}^0, H_u^0, H_c^0, B_e^0 > 0$, then the solutions $(P_u(t), P_{uA}(t), P_c(t), P_{cA}(t), H_u(t), H_c(t), B_e(t))$ are also positive for all $t \ge 0$.

Proof. Firstly, by the continuous dependence of solutions with respect to initial values, we only need to prove that when $P_u^0, P_{uA}^0, P_c^0, P_{cA}^0, H_u^0, H_c^0, B_e^0 > 0$, $(P_u(t), P_{uA}(t), P_c(t), P_{cA}(t), H_u(t), H_c(t), B_e(t))$ the solutions are also positive for all $t \ge 0$. That is, the solutions remain in the positive cone if the initial conditions are in the positive cone of $\mathbb{R}^7$. Set

Clearly, $m(0) > 0$. Assuming that there exists a $t_1 > 0$ such that $m(t_1) = 0$ and $m(t) > 0$ for all $t \in [\ 0, t_1)$.

If $m(t_1) = P_u(t_1)$, from the first equation of model (2.1) it follows that $\frac{dP_u}{dt} \ge -(\alpha_p\beta_p(1-\eta)H_c(t)+\kappa_pB_e(t)+\gamma_u+\epsilon)P_u$ for all $t \in [\ 0, t_1)$. Since $H_c(t) > 0, B_e(t) > 0$ for all $t \in [\ 0, t_1)$, we have

which leads to a contradiction. Similar contradictions can be deduced in the cases of $m(t_1) = P_{uA}(t_1), m(t_1) = P_c(t_1), m(t_1) = P_{cA}(t_1), m(t_1) = H_u(t_1), m(t_1) = H_c(t_1), m(t_1) = B_e(t_1)$. Hence, the solutions remain in the positive cone if the initial conditions are in the positive cone $\mathbb{R}^7$.

Secondly, let $T(t) = P_u(t)+P_{uA}(t)+P_c(t)+P_{cA}(t)+H_u(t)+H_c(t)+B_e(t)$. Then

where $N_p = P_u(t)+P_{uA}(t)+P_c(t)+P_{cA}(t)$ and $N_h = H_u(t)+H_c(t)$, which implies that

So $B_e(t)$ is bounded by a fixed number

Let $N = N_p+N_h+M$, we have

Thus, the solutions remain bounded in a positive cone of $\mathbb{R}^7$, and the system induces a global semiflow in the positively invariant set of $\mathbb{R}^7$. This completes the proof.

Remark 2. Denote set G as follows

Then Lemma 1 implies that G is a positively invariant set with respect to model (2.1).

3.2.   Basic reproduction number

When $\theta_c$ = 0, $\theta_{cA}$ = 0, that is, there are no colonized patients admitted into hospitals, model (2.1) has a unique infection-free equilibrium (IFE) which is defined by

We derive the basic reproduction number $R_0$ for model (2.1) by using the techniques in Diekmann et al. ^[44] and van den Driessche and Watmough ^[45], which involves linearizing the original nonlinear ordinary differential equations at the infection-free equilibrium. Re-order the components of $E_0$ as

and set

where $\Omega = (\gamma_uP_u+\gamma_cP_c+\gamma_{uA}P_{uA}+\gamma_{cA}P_{cA})$,

Since $\theta_c$ = 0, $\theta_{cA}$ = 0, then we can derive that

Therefore, let $\omega_1 = \frac{\upsilon_p\gamma_{cA}+\upsilon_{pA}\epsilon}{\gamma_b\gamma_{cA}(\gamma_c+\epsilon)}$, $\omega_2 = \frac{\upsilon_{pA}}{\gamma_b\gamma_{cA}}$, $\omega_3 = \frac{\upsilon_h}{\gamma_b\mu_c}$, we have

and the basic reproductive number is defined as the spectral radius of $FV^{-1}$:

where

By Theorem 2 in van den Driessche and Watmough ^[45], we have the following theorem:

Theorem 3. If $R_0 < 1$, then the infection-free equilibrium $E_0$ is locally asymptotically stable; If $R_0 > 1$, then $E_0$ is unstable.

Moreover, from the proof of Theorem 2 in van den Driesshce and Watmough ^[45] or the proof of Lemma 2.1 in Wang and Zhao ^[46], we have the following observation: Denote

Let $s(J_1)$ be the maximum real part of the eigenvalues of $J_1$. Since $J_1$ is irreducible and has non-negative off-diagonal elements, $s(J_1)$ is a simple eigenvalue of $J_1$ with a positive eigenvector. Then we have the following corollary:

Corollary 4. There hold two equivalences:

3.3.   Vanishing of infection

The existence and stability of the infection-free equilibrium $E_0$ indicates that the MRSA infection is vanishing.

Theorem 5. If $R_0 < 1$, then the infection-free equilibrium $E_0$ is globally asymptotically stable.

Proof. From Theorem 3 we know that $E_0$ is locally asymptotically stable. Now we prove the global attractivity of the infection-free equilibrium $E_0$.

By the first equation of model (2.1), non-negativity of the solutions and previous assumptions, we get

Since $\gamma_{uA} = \max\{\gamma_{uA}, \gamma_c, \gamma_{cA}\}$, it implies that

So $\forall \delta > 0$, there exists $t_1 > 0$, such that $P_u \le N^*+\delta,$ for all $t\ge t_1$.

Similarly, by the third equation of model (2.1), non-negativity of the solutions and previous assumptions, we get

that is,

Then $\forall \delta > 0$, there exists $t_2 > 0$, such that $P_{uA} \le N_p-N^*+\delta,$ for all $t\ge t_2$.

Let $T = \max \{ t_1, t_2\}.$ If $t > T$, since $\theta_c = \theta_{cA} = 0$, then

Considering the following auxiliary system:

Define

It follows from Corollary 4 that if $R_0 < 1$, then $s(J_1(0)) < 0$. Since $s(J_1(\delta))$ is continuous for small $\delta$, so there exists $\delta$ small enough such that $s(J_1(\delta)) < 0$. Thus there is a negative eigenvalue of $s(J_1(\delta))$ with a positive eigenvector. Obviously if $t \to \infty$, then $\tilde{ P_c}, \tilde{P_{cA}}, \tilde{H_c}, \tilde{B_e } \to 0$. Then by the comparison principle we get

Therefore, $E_0$ is globally attractive when $R_0 < 1$. This completes the proof.

3.4.   Uniform persistence

Uniform persistence of system (2.1) demonstrates that all components of the dynamical model have positive lower bounds which in turn indicates that MRSA infection persists in the hospital.

Theorem 6. If $R_0 > 1$, then model (2.1) is uniformly persistent.

Proof. We first define

It can be seen that both $X$ and $X_0$ are positively invariant with respect to model (2.1). Clearly, $\partial X_0$ is relatively closed in $X$. Lemma 1 implies that model (2.1) is point dissipative, which implies that the solutions of model (2.1) admit a global attractor. Then we define

Now we prove that

For any point $\varphi_0 = (P_u(0), P_c(0), P_{uA}(0), P_{cA}(0), H_u(0), H_c(0), B_e(0))$ in $M_{\partial}$, we suppose that one of $P_c(0), P_{cA}(0), H_c(0), B_e(0)$ is not zero, that is to say, $\varphi_0 \notin \{ (P_u, 0, P_{uA}, 0, H_u, 0, 0):P_u \ge 0, P_{uA}\ge0, H_u = N_h \}$. Without loss of generality, we suppose that $P_c(0) = 0, P_{cA}(0) = 0, H_c(0) = 0, B_e(0) > 0$. By the second, fourth, and sixth equations, we have

Thus, there exists $\delta_0 > 0$, if $0 < t < \delta_0$ then $P_c(t) > 0, P_{cA}(t) > 0, H_c(t) > 0, B_e(t) > 0$, which imply that $\varphi_0 \notin \partial X_0$. we will get the similar result for other cases ($P_c(0) > 0$, or $P_{cA}(0) > 0$, or $H_c(0) > 0$). Thus $\varphi_0 \notin M_{\partial}$. This gives us a contradiction. Hence $\varphi_0 \in \{ (P_u, 0, P_{uA}, 0, H_u, 0, 0):P_u \ge 0, P_{uA}\ge0, H_u = N_h \}$. So $M_{\partial} \subseteq \{ (P_u, 0, P_{uA}, 0, H_u, 0, 0):P_u \ge 0, P_{uA}\ge0, H_u = N_h \}$. Obviously we have $\{ (P_u, 0, P_{uA}, 0, H_u, 0, 0):P_u \ge 0, P_{uA}\ge0, H_u = N_h \} \subseteq M_{\partial}$, therefore, $M_{\partial} = \{ (P_u, 0, P_{uA}, 0, H_u, 0, 0):P_u \ge 0, P_{uA}\ge0, H_u = N_h \}$. Let $\varphi_0$ be an initial value. Clearly there is only one equilibrium $E_0 = (N^*, 0, N_p - N^*, 0, N_h, 0, 0)$ in $M_{\partial}$, so $\cup_{\varphi_0 \in M_{\partial}}\omega(\varphi_0) = E_0$. Therefore, $\{E_0\}$ is a compact and isolated invariant set in $\partial X_0$.

Next we claim that there exists a positive constant $\ell$ such that for any solution of model (2.1), $\Psi_t(\varphi_0), \varphi_0 \in X_0$, we have

where $d$ is a distant function in $X_0$. We construct by contradiction so that we suppose the claim is not true. Then $\limsup_{t \to\infty} d(\Psi_t(\varphi_0), E_0) \le \ell$ for any $\ell > 0$, namely, there exists a positive constant $T$, such that $N^*-\ell \le P_u(t) \le N^*+\ell, \ P_c(t) \le \ell, \ N_p-N^*-\ell \le P_{uA}(t) \le N_p-N^*+\ell, \ P_{cA}(t) \le \ell, \ N_h-\ell \le H_u(t) \le N_h+\ell, \ H_c(t) \le \ell, \ B_e(t) \le \ell,$ for any $t > T$. While $t > T$, we have,

Consider the following auxiliary system:

we define

For $R_0 > 1$, by Corollary 4, we have $s(J_1(0)) > 0$. Since $s(J_1(\ell))$ is continuous for small $\ell$, so there exists a positive constant $\ell$ small enough such that $s(J_1(\ell)) > 0$. Thus, there is a positive eigenvalue of $s(J_1(\delta))$ with a positive eigenvector. It is easy to see if $t \to \infty$, then $\tilde{ P_c}, \tilde{P_{cA}}, \tilde{H_c}, \tilde{B_e } \to \infty$. Then by the comparison principle we get

This contradicts our assumption and completes the proof of the claim.

The claim implies that $\{E_0\}$ is an isolated invariant set in $X$ and $W^s(E_0) \cap X_0 = \emptyset$. Therefore, system (2.1) is uniformly persistent if $R_0 > 1$ by Theorem 1.3.1 in ^[47]. This completes the proof.

4.   Numerical simulations of the deterministic model

In this section we present numerical simulations on the deterministic model and sensitivity analysis of the basic reproduction number in terms of model parameters.

Our deterministic model is simulated for 365 days. Data ^[24] collected in Beijing Tongren Hospital, where a total of 23 beds were in the emergency ward and were always fully occupied, are used to estimate the initial values of patients and healthcare workers. We assume an initial bacteria density being $1000\ ACC/cm^2$ as comparable to the measurement scale obtained by Bogusz et al. ^[48]. With the initial values $(P_u^0, P_{uA}^0, P_c^0, P_{cA}^0, H_u^0, H_c^0, B_e^0) = (4, 6, 7, 6, 17, 6, 1000)$ and parameter values shown in Table 1, we simulate the following outcomes: numerical solutions of the deterministic model (2.1), prevalences of colonized patients without or with antibiotic exposure, and the basic reproduction number $R_0$. Simulations are also performed to evaluate the effect of various interventions on changing the prevalence of colonized patients and $R_0$.

4.1.   Behavior of the model

Using the baseline parameters in Table 1, Figures 2 and 3 give the behaviors of solutions to the deterministic model (2.1), which imply that 36$\%$ of patients are colonized with MRSA with antibiotic exposure, and 4$\%$ are colonized without antibiotic exposure; Figure 4 shows that, while with no admission of MRSA-positive patients ($\theta_c = \theta_{cA} = 0$), 21$\%$ of patients are colonized with MRSA with antibiotic exposure and 3$\%$ are colonized without antibiotic exposure; while with no admission of patients with history of antibiotic exposure ($\theta_{uA} = \theta_{cA} = 0$), 27$\%$ of patients are colonized with MRSA with antibiotic exposure, and 7.5$\%$ are colonized without antibiotic exposure; while with no admission of patients with history of antibiotic exposure and MRSA-positive ($\theta_{uA} = \theta_c = \theta_{cA} = 0$), 14$\%$ of patients are colonized with MRSA with antibiotic exposure, and 3.5$\%$ are colonized without antibiotic exposure. Hence, to control hospital infections, we may need to reduce the proportion of colonized patients ($\theta_c$ and $\theta_{cA}$) at admission by increasing the detection and isolation of the admitted MRSA patients and also reduce the proportion of uncolonized patients with antibiotic exposure ($\theta_{uA}$) by strengthening the public education about how to use antibiotics properly in the community.

4.2.   Basic reproduction number

In the case where colonized patients are admitted into hospitals, the infections will always persist. When $\theta_c = 0, \theta_{cA} = 0$, that is no colonized patients are admitted into hospital, the infection-free equilibrium (IFE) is defined to be $E_0 = (P_u, P_c, P_{uA}, P_{cA}, H_u, H_c, B_e) = (N^*, 0, N_p - N^*, 0, N_h, 0, 0)$ where $N^*$ = $\frac{\theta_u\gamma_{uA}N_p}{\theta_{uA}\gamma_u+\theta_u\gamma_{uA}+\epsilon}$. By parameters listed in Table 1, the basic reproduction number is estimated to be 1.2860, which means that the infections are persistent. We want to reduce $R_0$ to below unity by some interventions. Here we perform some simulations to evaluate the effect of the following interventions in reducing the prevalence of colonized patients with or without antibiotic exposure, and $R_0$: (1) Prescription rate of antibiotics $\epsilon$; (2) Hand hygiene compliance of HCWS $\eta$; (3) The discharge rate for colonized patients with or without antibiotic exposure $\gamma_c$ and $\gamma_{cA}$, respectively; (4) Environmental cleaning rate $\gamma_b$; and (5) Decontamination rate of HCWs $\mu_c$.

The predicted effects of individual interventions on reducing the prevalence of MRSA and the reproduction number $R_0$ are shown in Figure 5. Figure 5A shows that increasing the compliance rate of hand hygiene for HCWs, $\eta$, from 0.4 (baseline) to 1, just reduces $R_0$ from 1.2860 to 1.2197, and reduces the prevalence of colonized patients with or without antibiotic exposure by 4.51$\%$ (from 20.56$\%$ to 16.04$\%$) and 0.54$\%$ (from 2.45$\%$ to 1.91$\%$), respectively. When antibiotic prescribing rate is reduced from 0.12 (baseline) to 0 (no antibiotic use), we get a result in around 19$\%$ reduction in the prevalence of colonized patients with antibiotic exposure, while a little increase and then decrease in the prevalence of colonized patients without antibiotic exposure, and a change from 1.2860 to 0.9251 in $R_0$ in Figure 5B. We investigate the discharge rate (i.e., the reciprocal of the length of stay) of colonized patients without antibiotic exposure $\gamma_c$, and with antibiotic exposure $\gamma_{cA}$, respectively, in Figures 5C-5D. When the discharge rate of $P_c$ is increased from the baseline value 0.06 to 0.2 (i.e., the length of stay of $P_c$ is decreased from 16.6 days to 5 days), $R_0$ reduces to 1.1308, and the prevalence of $P_{cA}$ and $P_c$ reduces by 9.58$\%$ (from 21.08$\%$ to 11.50$\%$) and 1.72$\%$ (from 2.57$\%$ to 0.85$\%$), respectively. Especially, we notice that if we decrease the discharge rate of $P_{cA}$ a little bit from the baseline value 0.055, there are dramatic increases in both $R_0$ and the prevalence of $P_{cA}$. However, many studies show that colonized patients with antibiotic exposure $P_{cA}$ usually lead to a lengthier stay ^[15], which in turns makes the situation worse. Hence, the rapid and efficient treatment of colonized patients, especially those with antibiotic exposure, is key in controlling MRSA infections. Furthermore, we find that improving environmental cleaning rate $\gamma_b$ is the most effective intervention from Figure 5E. When the environmental cleaning rate is increased from 0.7 (baseline) to 1, we are able to decrease the prevalence of $P_{cA}$ and $P_c$ from 20.56$\%$ to 1.99$\%$ and from 2.45$\%$ to 0.21$\%$, respectively, and successfully reduce $R_0$ to below unity. Figure 5F shows that decontamination rate of HCWs has little effect.

Furthermore, Figure 6 presents a direct simulation of the stability transition at $R_0 = 1$ to support our above conclusion that when $R_0 < 1$ the infection-free equilibrium is globally asymptotically stable and when $R_0 > 1$ the infection is uniformly persistent.

Observing that individual invention is hard to reduce $R_0$ to below unity, we examine the effects of combined interventions (Figure 7). When we decrease antibiotic use and in the meanwhile increase the discharge rate of $P_{cA}$, we reduce $R_0$ to below unity efficiently (Figure 7(b)). Similar result occurs when combining the increased environmental cleaning rate and decreased discharge rate of $P_{cA}$ (Figure 7(f)).

4.3.   Sensitivity analysis

Latin hypercube sampling (LHS) method is used to engage a sensitivity analysis ^[49,50]. Partial rank correlation coefficients (PRCCs) are calculated for the following nine parameters against the prevalence of colonized patients and $R_0$ over time: discharge rate for colonized patients with antibiotic exposure $\gamma_{cA}$, environmental cleaning rate $\gamma_b$, probability of colonization for $P_{uA}$ after a contact with a contaminated HCW $\beta_{pA}$, probability of contamination for HCW after a contact with a colonized patient with antibiotic exposure $\beta_{hA}$, hand hygiene compliance rate $\eta$, decontamination rate of HCWs $\mu_c$, shedding rate to the environment by colonized patients with antibiotic exposure $\upsilon_{pA}$, antibiotic prescribing rate $\epsilon$, contamination rate from environment for uncolonized patients with antibiotic exposure $\kappa_{pA}$. We also test for significant PRCCs for the above parameters to evaluate which parameters are essential to our model. Since we find that the PRCC values vary little after about 100 days, it is reasonable and efficient for us to just study the PRCC values on this specific day 100 (Figure 8). Figure 8(d) implies that the first four parameters have the most impact on the outcome of $R_0$, which are the environmental cleaning rate $\gamma_b$, shedding rate to the environment by colonized patients with antibiotic exposure $\upsilon_{pA}$, contamination rate from the environment for uncolonized patients with antibiotic exposure $\kappa_{pA}$, and antibiotic prescribing rate $\epsilon$. From Figures 8(a)-8(c), we illustrate the PRCC values of the nine examined parameters and corresponding p-values for different outcome parameters for $P_c$, $P_{cA}$, and $B_e$. All simulations are done by MATLAB and input parameters are assumed to be normally distributed, due to the lack of present data concerning distribution functions, as shown in Table 2.

5.   Stochastic simulations

Finally we run some numerical simulations of the stochastic model. Using the baseline parameters in Table 1 and initial values $(P_u^0, P_{uA}^0, P_c^0, P_{cA}^0, H_u^0, H_c^0, B_e^0) = (4, 6, 7, 6, 17, 6, 1000)$, we run 100 stochastic simulations of the SDE model (2.3). In Figure 9, the blue curves, red curves, and shaded regions represent the averages of 100 runs, outputs of the deterministic model (2.1), and $90\%$ bound of 100 runs, respectively. It is shown that the average of stochastic runs is consistent with the outcome of the deterministic model; however, randomness does make a difference in a single stochastic case. Also, we study the effect of environmental cleaning rate $\gamma_b$ and antibiotic prescribing rate $\epsilon$ on the number of colonized patients in Figures 10 and 11, respectively. Compared with Figure 9, we increase environmental cleaning rate, $\gamma_b = 1$, in Figure 10, which shows that increasing environmental cleaning rate can reduce the average number of colonized patients in the SDE model. Similarly, compared with Figure 9, we reduce the antibiotic use to an extreme case, $\epsilon = 0$, in Figure 11, which reduces the average number of colonized patients greatly in the SDE model.

Next, we consider the case with no admission of MRSA-positive patients, i.e., $\theta_c = \theta_{cA} = 0$, so $R_0$ can be calculated. We have proved that when $R_0 < 1$, MRSA infections will go to extinction in the deterministic model (2.1). By choosing different parameter values to make $R_0 < 1$ in both Figures 12 and 13, it is shown that MRSA infections do go to extinction in the deterministic model. However, the average number of colonized patients in the SDEs model persists, even though it is small, in Figures 12 and 13. In the above simulations, the number of colonized patients was treated as a real number rather than an integer. In order to intuitively see a difference between deterministic and stochastic models for a small population, in Figure 14, we round the real number to the nearest integer. This shows that with no admission of colonized patients, there are still reinfections after the number of colonized patients drops to zero, which indicates that the free-living bacteria may persist in the environment and later be transmitted back to the patients. That is to say, the free-living bacteria in the environment may be able to cause a later outbreak even though the infections have been died out from patients. Hence, hospitals should pay attention to environmental cleaning strategies to prevent MRSA infections.

Furthermore, in a relatively large population, where we consider initial values

Figure 15 shows that the average of multiple stochastic runs is consistent with the deterministic outcome.

6.   Discussion and conclusions

We developed a comprehensive study of MRSA infections in hospitals, which includes crucial factors such as antibiotic exposure and environmental contamination. Both deterministic and stochastic mathematical models were developed to study the transmission dynamics of MRSA infections in hospitals, including uncolonized patients without and with antibiotic exposure, colonized patients without and with antibiotic exposure, uncontaminated and contaminated health-care workers, and free-living MRSA. Under the assumption that there is no admission of the colonized patients, the basic reproduction number $R_0$ was calculated. It was shown that when $R_0 < 1$ the infection-free equilibrium is globally asymptotically stable, and when $R_0 > 1$ the infection is uniformly persistent.

For the deterministic model, numerical simulations were performed to demonstrate the behavior of the solutions, the effect of several interventions on reducing the prevalence of colonized patients and the basic reproduction number, and the dependence and sensitivity of the basic reproduction number of various parameters. Until recently, control strategies focus on the direct transmission between HCWs and patients, however, our results strongly supported that the environmental cleaning is the most effective intervention. When we increased the environmental cleaning rate, we could decrease the prevalence of colonized patients greatly and successfully reduce $R_0$ to below unity in Figure 5. Sensitively analysis in Figure 8 also showed that the environmental cleaning rate $\gamma_b$, shedding rate to the environment by colonized patients with antibiotic exposure $\upsilon_{pA}$, contamination rate from the environment for uncolonized patients with antibiotic exposure $\kappa_{pA}$, and antibiotic prescribing rate $\epsilon$ had remarkable impacts on the number of colonized patients and $R_0$. Even though it is difficult to quantify the environmental cleaning, we suggest that hospitals should try to use more effective cleaning products, improve monitoring strategies such as providing feedback to cleaning teams, and even use new technology (cleaning robots) to supplement the manual cleanings. Besides, it was shown that a higher discharge rate is associated with a lower prevalence of MRSA. The rapid and efficient treatment of colonized patients, especially those with antibiotic exposure, is key in controlling MRSA infections. However, the discharge rate depends on the time required for treatment and cannot be arbitrarily modified at will, which makes the control of MRSA infections challenge. In the cases of outbreaks, hospitals should try proper isolation of those colonized patients. Also, screening and isolating colonized patients at admission are important control strategies. Our study also emphasized the importance of effective antimicrobial stewardship programs in reducing antibiotic usage both in hospitals and communities.

For the stochastic model, numerical simulations were also carried out to study the behavior of the stochastic model and the effect of antibiotic prescribing rate $\epsilon$, and environmental cleaning rate $\gamma_b$, on the number of colonized patients, respectively. Moreover, we chose different parameter values to make $R_0 < 1$ and found that MRSA infections go to extinction in the deterministic model; however, the average number of colonized patients in multiple stochastic runs persisted in Figures 12 and 13. In order to intuitively see a difference between deterministic and stochastic models for a small population, in Figure 14, we rounded the real number to the nearest integer. This shows that the free-living bacteria in the environment may be able to cause a later outbreak even though the infections have been died out from patients. Hence, hospitals should pay attention to environmental cleaning strategies to prevent MRSA infections.

In the proposed model, the heterogeneity in infection risk among different types of wards was omitted and free-living bacteria were assumed to be uniformly distributed for the sake of simplicity; however, heterogeneity should be taken into account in the future work for a more realistic consideration ^[12,14]. Also, model parameters were regarded as a constant, however, it is not necessarily true in the real world. In applications, parameters in the model need to be inferred from noisy data. Recent works have developed several methods to infer parameters in models with complex and flexible structures ^[51,52,53], which should be potentially considered in future works. In addition, in current model the complex contact pattern was simplified to a full mixing, however, it may be useful to take into account the contact network structure in future works to represent heterogeneity in population behavior, locations and contact patterns ^[54,55]. In current model setting, patients were assumed to be hospitalized and under antibiotic treatments because of other diseases, then for those patients with antibiotic treatment, they were just more likely to be colonized by MRSA, and then had a lower probability of successful treatment, a lengthier stay in hospitals and extra costs of treatment. That is to say, we ignored the difference between colonized patients and infected patients since treatment should be needed for infected patients. A more comprehensive work may be done in the future ^[15]. Furthermore, an extension should be specifically compared with this proposed model and its conclusions and applied to actual MRSA infection data in hospitals (^[25,57]).

Acknowledgments

This research was partially supported by the University of Miami Provost's Research Award (UM PRA 2019-409). We would like to thank Dr. Xi Huo and three anonymous reviewers for their helpful comments and suggestions which helped us to improve the presentation of the paper.

Conflict of interest

All authors declare no conflicts of interest in this paper.