Loading [MathJax]/jax/output/SVG/jax.js
Research article

Identification of a novel functional miR-143-5p recognition element in the Cystic Fibrosis Transmembrane Conductance Regulator 3’UTR

  • Received: 01 September 2017 Accepted: 08 February 2018 Published: 23 February 2018
  • MicroRNAs (miRNAs) are small non-coding RNAs involved in regulation of gene expression. They bind in a sequence-specific manner to miRNA recognition elements (MREs) located in the 3’ untranslated region (UTR) of target mRNAs and prevent mRNA translation. MiRNA expression is dysregulated in cystic fibrosis (CF), affecting several biological processes including ion conductance in the epithelial cells of the lung. We previously reported that miR-143 is up-regulated in CF bronchial brushings compared to non-CF. Here we identified two predicted binding sites for miR-143-5p (starting at residues 558 and 644) on the CFTR mRNA, and aimed to assess whether CFTR is a true molecular target of miR-143-5p. Expression of miR-143-5p was found to be up-regulated in a panel of CF vs non-CF cell lines (1.7-fold, P = 0.0165), and its levels were increased in vitro after 20 hours treatment with bronchoalveolar lavage fluid from CF patients compared to vehicle-treated cells (3.3-fold, P = 0.0319). Luciferase assays were performed to elucidate direct miRNA::target interactions and showed that miR-143-5p significantly decreased the reporter activity when carrying the wild-type full length sequence of CFTR 3’UTR (minus 15%, P = 0.005). This repression was rescued by the disruption of the first, but not the second, predicted MRE, suggesting that the residue starting at 558 was the actual active binding site. In conclusion, we here showed that miR-143-5p modestly but significantly inhibits CFTR, improving the knowledge on functional MREs within the CFTR 3’UTR. This could lead to the development of novel therapeutic strategies where miRNA-mediated CFTR repression is blocked thereby possibly increasing the efficacy of the currently available CFTR modulators.

    Citation: Chiara De Santi, Sucharitha Gadi, Agnieszka Swiatecka-Urban, Catherine M. Greene. Identification of a novel functional miR-143-5p recognition element in the Cystic Fibrosis Transmembrane Conductance Regulator 3’UTR[J]. AIMS Genetics, 2018, 5(1): 53-62. doi: 10.3934/genet.2018.1.53

    Related Papers:

    [1] Avner Friedman, Najat Ziyadi, Khalid Boushaba . A model of drug resistance with infection by health care workers. Mathematical Biosciences and Engineering, 2010, 7(4): 779-792. doi: 10.3934/mbe.2010.7.779
    [2] Hermann Mena, Lena-Maria Pfurtscheller, Jhoana P. Romero-Leiton . Random perturbations in a mathematical model of bacterial resistance: Analysis and optimal control. Mathematical Biosciences and Engineering, 2020, 17(5): 4477-4499. doi: 10.3934/mbe.2020247
    [3] Michele L. Joyner, Cammey C. Manning, Brandi N. Canter . Modeling the effects of introducing a new antibiotic in a hospital setting: A case study. Mathematical Biosciences and Engineering, 2012, 9(3): 601-625. doi: 10.3934/mbe.2012.9.601
    [4] Zongmin Yue, Yuanhua Mu, Kekui Yu . Dynamic analysis of sheep Brucellosis model with environmental infection pathways. Mathematical Biosciences and Engineering, 2023, 20(7): 11688-11712. doi: 10.3934/mbe.2023520
    [5] Tianfang Hou, Guijie Lan, Sanling Yuan, Tonghua Zhang . Threshold dynamics of a stochastic SIHR epidemic model of COVID-19 with general population-size dependent contact rate. Mathematical Biosciences and Engineering, 2022, 19(4): 4217-4236. doi: 10.3934/mbe.2022195
    [6] Robert E. Beardmore, Rafael Peña-Miller . Rotating antibiotics selects optimally against antibiotic resistance, in theory. Mathematical Biosciences and Engineering, 2010, 7(3): 527-552. doi: 10.3934/mbe.2010.7.527
    [7] Damilola Olabode, Libin Rong, Xueying Wang . Stochastic investigation of HIV infection and the emergence of drug resistance. Mathematical Biosciences and Engineering, 2022, 19(2): 1174-1194. doi: 10.3934/mbe.2022054
    [8] Xiaxia Kang, Jie Yan, Fan Huang, Ling Yang . On the mechanism of antibiotic resistance and fecal microbiota transplantation. Mathematical Biosciences and Engineering, 2019, 16(6): 7057-7084. doi: 10.3934/mbe.2019354
    [9] Tsanou Berge, Samuel Bowong, Jean Lubuma, Martin Luther Mann Manyombe . Modeling Ebola Virus Disease transmissions with reservoir in a complex virus life ecology. Mathematical Biosciences and Engineering, 2018, 15(1): 21-56. doi: 10.3934/mbe.2018002
    [10] Jing Jia, Yanfeng Zhao, Zhong Zhao, Bing Liu, Xinyu Song, Yuanxian Hui . Dynamics of a within-host drug resistance model with impulsive state feedback control. Mathematical Biosciences and Engineering, 2023, 20(2): 2219-2231. doi: 10.3934/mbe.2023103
  • MicroRNAs (miRNAs) are small non-coding RNAs involved in regulation of gene expression. They bind in a sequence-specific manner to miRNA recognition elements (MREs) located in the 3’ untranslated region (UTR) of target mRNAs and prevent mRNA translation. MiRNA expression is dysregulated in cystic fibrosis (CF), affecting several biological processes including ion conductance in the epithelial cells of the lung. We previously reported that miR-143 is up-regulated in CF bronchial brushings compared to non-CF. Here we identified two predicted binding sites for miR-143-5p (starting at residues 558 and 644) on the CFTR mRNA, and aimed to assess whether CFTR is a true molecular target of miR-143-5p. Expression of miR-143-5p was found to be up-regulated in a panel of CF vs non-CF cell lines (1.7-fold, P = 0.0165), and its levels were increased in vitro after 20 hours treatment with bronchoalveolar lavage fluid from CF patients compared to vehicle-treated cells (3.3-fold, P = 0.0319). Luciferase assays were performed to elucidate direct miRNA::target interactions and showed that miR-143-5p significantly decreased the reporter activity when carrying the wild-type full length sequence of CFTR 3’UTR (minus 15%, P = 0.005). This repression was rescued by the disruption of the first, but not the second, predicted MRE, suggesting that the residue starting at 558 was the actual active binding site. In conclusion, we here showed that miR-143-5p modestly but significantly inhibits CFTR, improving the knowledge on functional MREs within the CFTR 3’UTR. This could lead to the development of novel therapeutic strategies where miRNA-mediated CFTR repression is blocked thereby possibly increasing the efficacy of the currently available CFTR modulators.


    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.

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

    Figure 1.  Flowchart of the model consisted of uncolonized patients without antibiotic exposure $ P_u(t) $, uncolonized patients with antibiotic exposure $ P_{uA}(t) $, colonized patients without antibiotic exposure $ P_c(t) $, colonized patients with antibiotic exposure $ P_{cA}(t) $, uncontaminated healthcare workers $ H_u(t) $, contaminated healthcare workers ($ H_c(t) $), and free-living bacteria in the environment $ B_e(t) $.

    $ 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:

    Table 1.  Parameters and descriptions.
    SymbolDescriptionValueReference
    $ \epsilon_0 $Antibiotic prescription rate (day$ ^{-1} $)0.12[39][41]
    $ \epsilon_1 $Magnitude of change of antibiotic prescription rate (no dimension)0.25[56]
    $ \theta_u $Proportion of $ P_u $ on admission (day$ ^{-1} $)0.617 [15][39]
    $ \theta_{uA} $Proportion of $ P_{uA} $ on admission (day$ ^{-1} $)0.349 [15][40]
    $ \theta_c $Proportion of $ P_c $ on admission (day$ ^{-1} $)0.003 [15][39]
    $ \theta_{cA} $Proportion of $ P_{cA} $ on admission (day$ ^{-1} $)0.031 [15][40]
    $ \gamma_u $Discharge rate of $ P_u $ (day$ ^{-1} $)0.2 [15]
    $ \gamma_{uA} $Discharge rate of $ P_{uA} $ (day$ ^{-1} $)0.2[15]
    $ \gamma_c $Discharge rate of $ P_c $ (day$ ^{-1} $)0.06[39]
    $ \gamma_{cA} $Discharge rate of $ P_{cA} $ (day$ ^{-1} $)0.055[15][39]
    $ \gamma_b $Disinfection (cleaning) rate of environment (day$ ^{-1} $)0.7[25]
    $ \alpha_p $Contact rate (day$ ^{-1} $ person$ ^{-1} $)0.0435 [25]
    $ \beta_p $Probability of colonization for $ P_u $ after a contact with $ H_c $ (no dimension)0.42 [25]
    $ \beta_{pA} $Probability of colonization for $ P_{uA} $ after a contact with $ H_c $ (no dimension)0.42*1.67[39][15]
    $ \beta_h $Probability of contamination for HCW after a contact with $ P_c $ (no dimension)0.2 [25][15]
    $ \beta_{hA} $Probability of contamination for HCW after a contact with $ P_{cA} $ (no dimension)0.25[15]
    $ \eta $Hand hygiene compliance with HCWs (no dimension)0.4[25]
    $ \mu_c $Decontamination rate of HCWs (day$ ^{-1} $)24 [25]
    $ \upsilon_p $Shedding rate to environment from $ P_c $ (day$ ^{-1} $ person$ ^{-1} $ ACC/$ cm^2 $)235[25]
    $ \upsilon_{pA} $Shedding rate to environment from $ P_{cA} $ (day$ ^{-1} $ person$ ^{-1} $ ACC/$ cm^2 $)470[26][39]
    $ \upsilon_h $Contamination rate to environment by $ H_c $ (day$ ^{-1} $ person$ ^{-1} $ ACC/$ cm^2 $)235[25]
    $ \kappa_p $Colonization rate from environment for $ P_u $ (day$ ^{-1} $ (ACC/$ cm^2 $)$ ^{-1} $)0.000004[25]
    $ \kappa_{pA} $Colonization rate from environment for $ P_{uA} $ (day$ ^{-1} $ (ACC/$ cm^2 $)$ ^{-1} $)0.000005 [15][25]
    $ \kappa_h $Colonization rate from environment for $ H_u $ (day$ ^{-1} $ (ACC/$ cm^2 $)$ ^{-1} $)0.00001[25]
    $ N_p $Total number of patients23[25]
    $ N_h $Total number of HCWs23[25]

     | Show Table
    DownLoad: CSV
    $ dPudt=θu(γuPu+γcPc+γuAPuA+γcAPcA)αpβp(1η)PuHcκpPuBeγuPuϵPu,dPcdt=θc(γuPu+γcPc+γuAPuA+γcAPcA)+αpβp(1η)PuHc+κpPuBeγcPcϵPc,dPuAdt=θuA(γuPu+γcPc+γuAPuA+γcAPcA)αpβpA(1η)PuAHcκpAPuABeγuAPuA+ϵPu,dPcAdt=θcA(γuPu+γcPc+γuAPuA+γcAPcA)+αpβpA(1η)PuAHc+κpAPuABeγcAPcA+ϵPc,dHudt=αpβh(1η)PcHuαpβhA(1η)PcAHuκhHuBe+μcHc,dHcdt=αpβh(1η)PcHu+αpβhA(1η)PcAHu+κhHuBeμcHc,dBedt=υpPc+υpAPcA+υhHcγbBe $ (2.1)

    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.

    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

    $ p_{(s, j, k, m, n)}(t) = \Pr(P_c(t) = s, P_{uA}(t) = j, P_{cA}(t) = k, H_c(t) = m, B_e(t) = n) $

    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:

    $ p(s+i1,j+i2,k+i3,m+i4,n+i5);(s,j,k,m,n)(t)=Pr[(Pc(t+t),PuA(t+t),PcA(t+t),Hc(t+t),Be(t+t))=(s+i1,j+i2,k+i3,m+i4,n+i5)|(Pc(t),PuA(t),PcA(t),Hc(t),Be(t))=(s,j,k,m,n)], $

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

    ${\mathit{\boldsymbol{p_{(s+i_1, j+i_2, k+i_3, m+i_4, n+i_5);(s, j, k, m, n)}(\triangle t)}}}\\ ={{θc(γu(Npsjk)+γcs+γuAj+γcAk)+αpβp(1η)(Npsjk)m+κp(Npsjk)n}t(i1,i2,i3,i4,i5)=(1,0,0,0,0)γcst(i1,i2,i3,i4,i5)=(1,0,0,0,0)ϵst(i1,i2,i3,i4,i5)=(1,0,1,0,0){θuA(γu(Npsjk)+γcs+γuAj+γcAk)+ϵ(Npsjk)}t(i1,i2,i3,i4,i5)=(0,1,0,0,0)(κpAjn+γuAj)t(i1,i2,i3,i4,i5)=(0,1,0,0,0)αpβpA(1η)jmt(i1,i2,i3,i4,i5)=(0,1,1,0,0){θcA(γu(Npsjk)+γcs+γuAj+γcAk)+κpjn}t(i1,i2,i3,i4,i5)=(0,0,1,0,0)γcAkt(i1,i2,i3,i4,i5)=(0,0,1,0,0){αpβh(1η)s(Nhm)+αpβhA(1η)k(Nhm)+κh(Nhm)n}t(i1,i2,i3,i4,i5)=(0,0,0,1,0)μcmt(i1,i2,i3,i4,i5)=(0,0,0,1,0)(υps+υpAk+υhm)t(i1,i2,i3,i4,i5)=(0,0,0,0,1)γbnt(i1,i2,i3,i4,i5)=(0,0,0,0,1)0otherwise. $ (2.2)

    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:

    $ p(s,j,k,m,n)(t+t)=p(s1,j,k,m,n)(t)[θc(γu(Nps+1jk)+γcs+γuAj+γcAk)+αpβp(1η)(Nps+1jk)m+κp(Nps+1jk)n]t+p(s+1,j,k,m,n)(t)γc(s+1)t+p(s+1,j,k1,m,n)(t)ϵ(s+1)t+p(s,j1,k,m,n)(t)[θuA(γu(Npsj+1k)+γcs+γuAj+γcAk)+ϵ(Npsj+1k))]t+p(s,j+1,k,m,n)(t)[κpA(j+1)n+γuA(j+1)]t+p(s,j+1,k1,m,n)(t)αpβpA(1η)(j+1)m+p(s,j,k1,m,n)(t)[θc(γu(Npsjk+1)+γcs+γuAj+γcA(k+1))+κpjn]t+p(s,j,k+1,m,n)(t)γcA(k+1)t+p(s,j,k,m1,n)(t)[αpβh(1η)s(Nhm+1)+αpβhA(1η)k(Nhm+1)+κh(Nhm+1)n]t+p(s,j,k,m+1,n)(t)μc(m+1)+p(s,j,k,m,n1)(t)(υps+υpAk+υhm)t+p(s,j,k,m,n+1)(t)γb(n+1)t+(t). $

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

    $ dps,j,k,m,ndt=p(s1,j,k,m,n)[θc(γu(Nps+1jk)+γcs+γuAj+γcAk)+αpβp(1η)(Nps+1jk)m+κp(Nps+1jk)n]+p(s+1,j,k,m,n)(t)γc(s+1)+p(s+1,j,k1,m,n)(t)ϵ(s+1)+p(s,j1,k,m,n)(t)[θuA(γu(Npsj+1k)+γcs+γuAj+γcAk)+ϵ(Npsj+1k))]+p(s,j+1,k,m,n)(t)[κpA(j+1)n+γuA(j+1)]t+p(s,j+1,k1,m,n)(t)αpβpA(1η)(j+1)m+p(s,j,k1,m,n)(t)[θc(γu(Npsjk+1)+γcs+γuAj+γcA(k+1))+κpjn]+p(s,j,k+1,m,n)(t)γcA(k+1)+p(s,j,k,m1,n)(t)[αpβh(1η)s(Nhm+1)+αpβhA(1η)k(Nhm+1)+κh(Nhm+1)n]+p(s,j,k,m+1,n)(t)μc(m+1)+p(s,j,k,m,n1)(t)(υps+υpAk+υhm)+p(s,j,k,m,n+1)(t)γb(n+1). $

    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:

    $ p_{(s, j, k, m, n)}(t) = \Pr\{ P_c(t) = s, P_{uA}(t) = j, P_{cA}(t) = k, H_c(t) = m, B_e(t) = n \} $

    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

    $ \triangle X(t) = (\triangle P_c(t), \triangle P_{uA}(t), \triangle P_{cA}(t), \triangle H_c(t), \triangle B_e(t))^T. $

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

    $ E(\triangle X(t)|X(t) = (eceuAecAeheb) \triangle t = f(X(t), t)\triangle t, $

    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) $:

    $ E(X(t)(X(t))T|X(t))=(δc0ϵPc000δuAαpβpA(1η)PuAHc00ϵPcαpβpA(1η)PuAHcδcA00000δh00000δb)t=Σ(X(t),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 $,

    $ G = (a1a2a3000000000000a4a5a600000000a300a6a7a8000000000000a9a10000000000000a11a12), $

    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:

    $ dX(t) = f(X(t), t)dt+G(X(t), t)dW(t). $

    More precisely,

    $ {dPc(t)=ecdt+a1dW1a2dW2a3dW3,dPuA(t)=euAdt+a4dW4a5dW5a6dW6,dPcA(t)=ecAdt+a3dW3+a6dW6+a7dW7a8dW8,dHct=ehdt+a9dW9a10dW10,dBe(t)=ebdt+a11dW11a12dW12. $ (2.3)

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

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

    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:

    $ P_u^0 \ge 0, P_{uA}^0 \ge 0, P_c^0 \ge 0, P_{cA}^0 \ge 0, H_u^0 \ge 0, H_c^0 \ge 0, B_e^0 \ge 0. $

    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

    $ m(t) = \min\{P_u(t), P_{uA}(t), P_c(t), P_{cA}(t), H_u(t), H_c(t), B_e(t)\}, \forall t \gt 0. $

    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

    $ 0 = P_u(t_1) \ge P_u^0 \exp(-\int_{0}^{t_1}(\alpha_p\beta_p(1-\eta)H_c(s)+\kappa_pB_e(s)+\gamma_u+\epsilon)ds) \gt 0, $

    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

    $ dT(t)dt=dBe(t)dt=υpPc+υpAPcA+υhHcγbBeυpNp+υpANp+υhNhγbBe(t), $

    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

    $ B_e(t) \le \frac{(\upsilon_pN_p+\upsilon_{pA}N_p+\upsilon_hN_h)}{\gamma_b}(1-e^{-\gamma_b t})+B_e^0 e^{-\gamma_b t}. $

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

    $ M = \frac{(\upsilon_pN_p+\upsilon_{pA}N_p+\upsilon_hN_h)}{\gamma_b}+B_e^0. $

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

    $ T(t) = P_u(t)+P_{uA}(t)+P_c(t)+P_{cA}(t)+H_u(t)+H_c(t)+B_e(t) \le N. $

    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

    $ G: = \{(P_u, P_{uA}, P_c, P_{cA}, H_u, H_c, B_e) \in \mathbb{R}^7_{+} : P_u+P_{uA}+P_c+P_{cA}+H_u+H_c+B_e \le N )\}. $

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

    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

    $ 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), \;\; N^* = \frac{\theta_u\gamma_{uA}N_p}{\theta_{uA}\gamma_u+\theta_u\gamma_{uA}+\epsilon}. $

    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

    $ E_0 = (P_c, P_{cA}, H_c, B_e, P_u, P_{uA}, H_u) = (0, 0, 0, 0, N^*, N_p-N^*, N_h) $

    and set

    $ \mathcal{F} = (αpβp(1η)PuHc+κpPuBeαpβpA(1η)PuAHc+κpAPuABeαpβh(1η)PcHu+αpβhA(1η)PcAHu+κhHuBe0000), $
    $ \mathcal{V} = (γcPc+ϵPcθcΩγcAPcA[ϵPc+θcAΩ]μcHcγbBe(υpPc+υpAPcA+υhHc)αpβp(1η)PuHc+κpPuBe+γuPu+ϵPuθuΩαpβpA(1η)PuAHc+κpAPuABe+γuAPuA[ϵPu+θuAΩ]αpβh(1η)PcHu+αpβhA(1η)PcAHu+κhHuBeμcHc), $

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

    $ \mathcal{V}^- = (γcPc+ϵPcγcAPcAμcHcγbBeαpβp(1η)PuHc+κpPuBe+γuPu+ϵPuαpβpA(1η)PuAHc+κpAPuABe+γuAPuAαpβh(1η)PcHu+αpβhA(1η)PcAHu+κhHuBe), \;\; \mathcal{V}^+ = (θcΩϵPc+θcAΩ0υpPc+υpAPcA+υhHcθuΩϵPu+θuAΩμcHc). $

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

    $ F = (00αpβp(1η)NκpN00αpβpA(1η)(NpN)κpA(NpN)αpβh(1η)NhαpβhA(1η)Nh0κhNh0000), $
    $ V = (γc+ϵ000ϵγcA0000μc0υpυpAυhγb). $

    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

    $ FV^{-1} = (ω1κpNω2κpNω3κpN+αpβp(1η)NμcκpNγbω1κpA(NpN)ω2κpA(NpN)ω3κpA(NpN)+αpβpA(1η)(NpN)μcκpA(NpN)γbω1κhNh+βhγcA+βhAϵγcA(γc+ϵ)αp(1η)Nhω2κhNh+αpβhA(1η)NhγcAω3κhNhκhNhγb0000) $

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

    $ R0=sp(FV1)=α3α1+α1+α2 $ (3.1)

    where

    $ α1=((α3627μ3cγ3cA+α5+α4)2α33+α3627μ3cγ3cA+α5+α4)13,α2=α63μcγcA,  α3=α73μcγcA+α269μ2cγ2cA,  α4=α6α76μ2cγ2cA,α5=(βpκpAβpAκp)N(NpN)Nhα2p(1η)2[ ω1βhA(γc+ϵ)ω2(βhγcA+βhAϵ)] 2μcγcA(γc+ϵ),α6=μcγcA(ω1κpN+ω2κpA(NpN)+ω3κhNh),α7=Nhα2p(1η)2[βhAβpA(NpN)+βhγcA+βhAϵγc+ϵβpN]+(NpN)Nhαp(1η)[ω2κhβpAγcA+ω3κpAβhAμc]+NNhαp(1η)[ω1βpκhγcA+ω3κpμcβhγcA+βhAϵγc+ϵ]. $

    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

    $ J_1 = F-V = (γcϵ0αpβp(1η)NκpNϵγcAαpβpA(1η)(NpN)κpA(NpN)αpβh(1η)NhαpβhA(1η)NhμcκhNhυpυpAυhγb). $

    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:

    $ R_0 \lt 1 \iff s(J_1) \lt 0; \;\;\; R_0 \gt 1 \iff s(J_1) \gt 0. $

    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

    $ \frac{dP_u}{dt} \le \theta_u[\gamma_uP_u+\gamma_{uA}(N_p-P_u)]-\gamma_uP_u-\epsilon P_u. $

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

    $ \frac{dP_u}{dt} \le \theta_u\gamma_{uA}N_p-(-\theta_u\gamma_u+\theta_u\gamma_{uA}+\gamma_u+\epsilon)P_u = \theta_u\gamma_{uA}N_p-(\theta_{uA}\gamma_u+\theta_u\gamma_{uA}+\epsilon)P_u. $

    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

    $ \frac{dP_{uA}}{dt} \le (1-\theta_u)[\gamma_u(N_p-P_{uA})+\gamma_{uA}P_{uA}]-\gamma_{uA}P_{uA}+\epsilon(N_p-P_{uA}), $

    that is,

    $ \frac{dP_{uA}}{dt} \le (\theta_{uA}\gamma_u+\epsilon)N_p-(\theta_{uA}\gamma_u+\theta_u\gamma_{uA}+\epsilon)P_{uA}. $

    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

    $ {Pc(t)αpβp(1η)(N+δ)Hc+κp(N+δ)BeγcPcϵPc,PcA(t)αpβpA(NpN+δ)Hc+κpA(NpN+δ)BeγcAPcA+ϵPc,Hc(t)αpβh(1η)NhPc+αpβhA(1η)NhPcA+κhNhBeμcHc,Be(t)υpPc+υpAPcA+υhHcγbBe. $ (3.2)

    Considering the following auxiliary system:

    $ {~Pc(t)=αpβp(1η)(N+δ)˜Hc+κp(N+δ)˜Beγc˜Pcϵ˜Pc,~PcA(t)=αpβpA(NpN+δ)˜Hc+κpA(NpN+δ)˜BeγcA˜PcA+ϵ˜Pc,~Hc(t)=αpβh(1η)Nh˜Pc+αpβhA(1η)Nh˜PcA+κhNh˜Beμc˜Hc,~Be(t)=υp˜Pc+υpA˜PcA+υh˜Hcγb˜Be. $ (3.3)

    Define

    $ J_1(\delta) = (γcϵ0αpβp(1η)(N+δ)κp(N+δ)ϵγcAαpβpA(1η)(NpN+δ)κpA(NpN+δ)αpβh(1η)NhαpβhA(1η)NhμcκhNhυpυpAυhγb). $

    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

    $ \lim\limits_{t \to\infty} P_c = 0, \lim\limits_{t \to\infty} P_{cA} = 0, \lim\limits_{t \to\infty} H_c = 0, \lim\limits_{t \to\infty} B_e = 0. $

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

    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

    $ X = \{(P_u, P_c, P_{uA}, P_{cA}, H_u, H_c, B_e): P_u \ge 0, P_c \ge 0, P_{uA} \ge 0, P_{cA} \ge 0, H_u \ge0, H_c\ge0, B_e\ge0 \}, $
    $ X_0 = \{(P_u, P_c, P_{uA}, P_{cA}, H_u, H_c, B_e) \in X: P_c \gt 0, P_{cA} \gt 0, H_c \gt 0, B_e \gt 0 \}, \;\; \partial X_0 = X\backslash X_0. $

    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

    $ M={(Pu(0),Pc(0),PuA(0),PcA(0),Hu(0),Hc(0),Be(0)):(Pu(t),Pc(t),PuA(t),PcA(t),Hu(t),Hc(t),Be(t))X0,t0}. $

    Now we prove that

    $ M_{\partial} = \{ (P_u, 0, P_{uA}, 0, H_u, 0, 0):P_u \ge 0 , P_{uA}\ge0, H_u = N_h \}. $

    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

    $ \frac{dP_c(0)}{dt} \ge \kappa_pP_u(0)B_e(0) \gt 0; \frac{dP_{cA}(0)}{dt} \ge \kappa_{pA}P_{uA}(0)B_e(0) \gt 0; \frac{dH_c(0)}{dt} \ge \kappa_hH_u(0)B_e(0) \gt 0. $

    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

    $ \limsup\limits_{t \to\infty} d(\Psi_t(\varphi_0), E_0) \ge \ell, $

    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,

    $ {Pc(t)αpβp(1η)(N)Hc+κp(N)BeγcPcϵPc,PcA(t)αpβpA(NpN)Hc+κpA(NpN)BeγcAPcA+ϵPc,Hc(t)αpβh(1η)(Nh)Pc+αpβhA(1η)(Nh)PcA+κh(Nh)BeμcHc,Be(t)υpPc+υpAPcA+υhHcγbBe. $ (3.4)

    Consider the following auxiliary system:

    $ {~Pc(t)=αpβp(1η)(N)~Hc+κp(N)~Beγc~Pcϵ~Pc,~PcA(t)=αpβpA(NpN)~Hc+κpA(NpN)~BeγcA~PcA+ϵ~Pc,~Hc(t)=αpβh(1η)(Nh)~Pc+αpβhA(1η)(Nh)~PcA+κh(Nh)~Beμc~Hc,~Be(t)=υp~Pc+υpA~PcA+υh~Hcγb~Be. $ (3.5)

    we define

    $ J_1(\ell) = (γcϵ0αpβp(1η)(N)κp(N)ϵγcAαpβpA(1η)(NpN)κpA(NpN)αpβh(1η)(Nh)αpβhA(1η)(Nh)μcκh(Nh)υpυpAυhγb). $

    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

    $ \lim\limits_{t \to\infty} P_c = \infty, \; \lim\limits_{t \to\infty} P_{cA} = \infty, \; \lim\limits_{t \to\infty} H_c = \infty, \; \lim\limits_{t \to\infty} B_e = \infty. $

    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.

    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 $.

    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.

    Figure 2.  Solutions of uncolonized patients without or with antibiotic exposure ($ P_u(t), P_{uA}(t) $) and colonized patients without or with antibiotic exposure ($ P_c(t), P_{cA}(t) $) of deterministic model (2.1) with 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 $ \theta_u = 0.617, \theta_{uA} = 0.349, \theta_c = 0.003, \theta_{cA} = 0.031 $ on admission. All parameter values are given in Table 1.
    Figure 3.  (a) Prevalence of colonized patients with or without antibiotic exposure; (b) density of bacteria ($ ACC/cm^2 $) in the environment of deterministic model (2.1) with 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 $ \theta_u = 0.617, \theta_{uA} = 0.349, \theta_c = 0.003, \theta_{cA} = 0.031 $ on admission.
    Figure 4.  Model behaviors with various colonization ratios upon admission. (a) Prevalence of colonized patients without antibiotic exposure; (b) prevalence of colonized patients with antibiotic exposure (c) density of bacteria ($ ACC/cm^2 $) in the environment of deterministic model (2.1) with 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)$.

    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.

    Figure 5.  Effects of individual interventions on the prevalence of colonized patient without antibiotic exposure (dashed lines), colonized patients with antibiotic exposure (dashed-dot lines) and the basic reproduction number $ R_0 $ (solid lines). The following interventions are investigated: (A) compliance with hand hygiene; (B) antibiotic prescribing rate; (C) discharge rate of colonized patients without antibiotic exposure $ P_c $; (D) discharge rate of colonized patients with antibiotic exposure $ P_{cA} $; (E) environmental cleaning rate; and (F) decontamination rate of HCWs.

    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.

    Figure 6.  Model behaviors with different values of $ R_0 $. (a) Prevalence of colonized patients with or without antibiotic for an arbitrary choice of $ R_0 < 1 $; (b) prevalence of colonized patients with or without antibiotic for an arbitrary choice of $ R_0 > 1 $.

    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)).

    Figure 7.  Effects of two interventions on the basic reproduction number $ R_0$.

    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.

    Figure 8.  (a)–(c) PRCCs of the nine parameters for $ P_c, P_{cA}, B_e $ when t = 100 day; (d) PRCCs for $ R_0 $ when $ \theta_c = \theta_{cA} = 0 $. All the parameters come from Latin Hypercube sampling.
    Table 2.  Variables evaluated in the sensitivity analysis.
    Symbol Distribution reference
    $ \gamma_{cA} $ $ N(0.055, 0.005) $ estimated by [26]
    $ \gamma_b $ $ N(0.7, 0.2) $ estimated by [26]
    $ \beta_{pA} $ $ N(0.43, 0.1) $ estimated by [26]
    $ \beta_{hA} $ $ N(0.2, 0.05) $ estimated by [26]
    $ \eta $ $ N(0.4, 0.1) $ estimated by [26]
    $ \mu_c $ $ N(24, 5) $ estimated by [26]
    $ \upsilon_{pA} $ $ N(470,150) $ estimated by [26]
    $ \epsilon $ $ N(0.12, 0.02) $ estimated by [26]
    $ \kappa_{pA} $ $ N(0.000005, 0.0000006) $ estimated by [26]

     | Show Table
    DownLoad: CSV

    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.

    Figure 9.  One hundred runs of the SDE model (2.3) with parameter values shown 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) $. The blue curves represent the averages of 100 runs in each compartment, the red curves are the outputs of deterministic model (2.1), and the shaded regions represent $ 90\% $ bound of 100 SDE model simulations.
    Figure 10.  One hundred runs of the SDE model (2.3) with $ \gamma_b $ = 1, other parameter values shown 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) $. The blue curves represent the averages of 100 runs in each compartment, the red curves are the outputs of deterministic model (2.1), and the shaded regions represent $ 90\% $ bound of 100 SDE model simulations.
    Figure 11.  One hundred runs of the SDE model (2.3) with $ \epsilon $ = 0, other parameter values shown 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) $. The blue curves represent the averages of 100 runs in each compartment, the red curves are the outputs of deterministic model (2.1), and the shaded regions represent $ 90\% $ bound of 100 SDE model simulations.

    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.

    Figure 12.  One hundred runs of the SDE model (2.3) with $ R_0 < 1 $ ($ \gamma_b $ = 1, $ \theta_c = \theta_{cA} = 0 $), other parameter values shown 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) $. The blue curves represent the averages of 100 runs in each compartment, the red curves are the outputs of deterministic model (2.1), and the shaded regions represent $ 90\% $ bound of 100 SDE model simulations.
    Figure 13.  One hundred runs of the SDE model (2.3) with $ R_0 < 1 $ ($ \epsilon $ = 0, $ \theta_c = \theta_{cA} = 0 $, $ \gamma_b = 2 $), other parameter values shown 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) $. The blue curves represent the averages of 100 runs in each compartment, the red curves are the outputs of deterministic model (2.1), and the shaded regions represent $ 90\% $ bound of 100 SDE model simulations.
    Figure 14.  Behavior of the SDE model (2.3) (simulation at one time) with $ R_0 < 1 $ ($ \theta_c = \theta_{cA} = 0 $, $ \gamma_b = 1.5 $), other parameter values shown 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) $. The black curves are the outcomes of the deterministic model (2.1).

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

    $ (P_u^0, P_{uA}^0, P_c^0, P_{cA}^0, H_u^0, H_c^0, B_e^0) = (34, 16, 17, 16, 17, 6, 1000), $

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

    Figure 15.  One hundred runs of the SDE model (2.3) with initial values $ (P_u^0, P_{uA}^0, P_c^0, P_{cA}^0, H_u^0, H_c^0, B_e^0) = (34, 16, 17, 16, 17, 6, 1000) $ and parameter values shown in Table 1. The blue curves represent the averages of 100 runs in each compartment, the red curves are the outputs of deterministic model (2.1), and the shaded regions represent $ 90\% $ bound of 100 SDE model simulations.

    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]).

    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.

    All authors declare no conflicts of interest in this paper.

    [1] Lee RC, Ambros V (1993) An extensive class of small RNAs in Caenorhabditis elegans. Science 294: 862–864.
    [2] Chandra S, Vimal D, Sharma D, et al. (2017) Role of miRNAs in development and disease: Lessons learnt from small organisms. Life Sci 185: 8–14. doi: 10.1016/j.lfs.2017.07.017
    [3] Rupaimoole R, Calin GA, Lopez-Berestein G, et al. (2016) miRNA Deregulation in Cancer Cells and the Tumor Microenvironment. Cancer Discov 6: 235–246. doi: 10.1158/2159-8290.CD-15-0893
    [4] Szymczak I, Wieczfinska J, Pawliczak R (2016) Molecular Background of miRNA Role in Asthma and COPD: An Updated Insight. Biomed Res Int 2016: 7802521.
    [5] Chen JQ, Papp G, Szodoray P (2016) The role of microRNAs in the pathogenesis of autoimmune diseases. Autoimmun Rev 15: 1171–1180. doi: 10.1016/j.autrev.2016.09.003
    [6] Kiriakidou M, Nelson PT, Kouranov A, et al. (2004) A combined computational-experimental approach predicts human microRNA targets. Genes Dev 18: 1165–1178. doi: 10.1101/gad.1184704
    [7] Lai EC, Tam B, Rubin GM (2005) Pervasive regulation of Drosophila Notch target genes by GY-box-, Brd-box-, and K-box-class microRNAs. Genes Dev 19: 1067–1080. doi: 10.1101/gad.1291905
    [8] Vodicka P, Pardini B, Vymetalkova V, et al. (2016) Polymorphisms in Non-coding RNA Genes and Their Targets Sites as Risk Factors of Sporadic Colorectal Cancer. Adv Exp Med Biol 937: 123–149. doi: 10.1007/978-3-319-42059-2_7
    [9] Nossent AY, Hansen JL, Doggen C, et al. (2011) SNPs in microRNA binding sites in 3'-UTRs of RAAS genes influence arterial blood pressure and risk of myocardial infarction. Am J Hypertens 24: 999–1006. doi: 10.1038/ajh.2011.92
    [10] Wang X, Jiang H, Wu W, et al. (2017) An Integrating Approach for Genome-Wide Screening of MicroRNA Polymorphisms Mediated Drug Response Alterations. Int J Genomics 2017: 1674827.
    [11] Thomson DW, Bracken CP, Goodall GJ (2011) Experimental strategies for microRNA target identification. Nucleic Acids Res 39: 6845–6853. doi: 10.1093/nar/gkr330
    [12] Gillen AE, Gosalia N, Leir SH, et al. (2011) MicroRNA regulation of expression of the cystic fibrosis transmembrane conductance regulator gene. Biochem J 438: 25–32. doi: 10.1042/BJ20110672
    [13] Megiorni F, Cialfi S, Dominici C, et al. (2011) Synergistic post-transcriptional regulation of the Cystic Fibrosis Transmembrane conductance Regulator (CFTR) by miR-101 and miR-494 specific binding. PLoS One 6: e26601. doi: 10.1371/journal.pone.0026601
    [14] Hassan F, Nuovo GJ, Crawford M, et al. (2012) MiR-101 and miR-144 regulate the expression of the CFTR chloride channel in the lung. PLoS One 7: e50837. doi: 10.1371/journal.pone.0050837
    [15] Oglesby IK, Chotirmall SH, McElvaney NG, et al. (2013) Regulation of cystic fibrosis transmembrane conductance regulator by microRNA-145, -223, and -494 is altered in ΔF508 cystic fibrosis airway epithelium. J Immunol 190: 3354–3362. doi: 10.4049/jimmunol.1202960
    [16] Amato F, Seia M, Giordano S, et al. (2013) Gene mutation in microRNA target sites of CFTR gene: a novel pathogenetic mechanism in cystic fibrosis? PLoS One 8: e60448. doi: 10.1371/journal.pone.0060448
    [17] Ramachandran S, Karp PH, Osterhaus SR, et al. (2013) Post-transcriptional regulation of cystic fibrosis transmembrane conductance regulator expression and function by microRNAs. Am J Respir Cell Mol Biol 49: 544–551. doi: 10.1165/rcmb.2012-0430OC
    [18] Viart V, Bergougnoux A, Bonini J, et al. (2015) Transcription factors and miRNAs that regulate fetal to adult CFTR expression change are new targets for cystic fibrosis. Eur Respir J 45: 116–128. doi: 10.1183/09031936.00113214
    [19] World Health Organization. The molecular genetic epidemiology of cystic fibrosis: report of a joint meeting of WHO/ECFTN/ICF(M)A/ECFS; June 19, 2002; Genoa, Italy.
    [20] The Clinical and Functional TRanslation of CFTR (CFTR2), 2017. Available from: http://cftr2.org.
    [21] McKiernan PJ, Greene CM (2015) MicroRNA Dysregulation in Cystic Fibrosis. Mediators Inflamm 2015: 529642.
    [22] Oglesby IK, Bray IM, Chotirmall SH, et al. (2010) miR-126 is downregulated in cystic fibrosis airway epithelial cells and regulates TOM1 expression. J Immunol 184: 1702–1709. doi: 10.4049/jimmunol.0902669
    [23] Oglesby IK, Vencken SF, Agrawal R, et al. (2015) miR-17 overexpression in cystic fibrosis airway epithelial cells decreases interleukin-8 production. Eur Respir J 46: 1350–1360. doi: 10.1183/09031936.00163414
    [24] Mirković B, Murray MA, Lavelle GM, et al. (2015) The Role of Short-Chain Fatty Acids, Produced by Anaerobic Bacteria, in the Cystic Fibrosis Airway. Am J Respir Crit Care Med 192: 1314–1324. doi: 10.1164/rccm.201505-0943OC
    [25] Livak KJ, Schmittgen TD (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25: 402–408. doi: 10.1006/meth.2001.1262
    [26] Agarwal V, Bell GW, Nam JW, et al. (2015) Predicting effective microRNA target sites in mammalian mRNAs. Elife 4.
    [27] Betel D, Koppal A, Agius P, et al. (2010) Comprehensive modeling of microRNA targets predicts functional non-conserved and non-canonical sites. Genome Biol 11: R90. doi: 10.1186/gb-2010-11-8-r90
    [28] Miranda KC, Huynh T, Tay Y, et al. (2006) A pattern-based method for the identification of MicroRNA binding sites and their corresponding heteroduplexes. Cell 126: 1203–1217. doi: 10.1016/j.cell.2006.07.031
    [29] Clément T, Salone V, Rederstorff M (2015) Dual luciferase gene reporter assays to study miRNA function. Methods Mol Biol 1296: 187–198. doi: 10.1007/978-1-4939-2547-6_17
    [30] Karimi L, Mansoori B, Shanebandi D, et al. (2017) Function of microRNA-143 in different signal pathways in cancer: New insights into cancer therapy. Biomed Pharmacother 91: 121–131. doi: 10.1016/j.biopha.2017.04.060
    [31] He M, Zhan M, Chen W, et al. (2017) MiR-143-5p Deficiency Triggers EMT and Metastasis by Targeting HIF-1α in Gallbladder Cancer. Cell Physiol Biochem 42: 2078–2092. doi: 10.1159/000479903
    [32] Wu XL, Cheng B, Li PY, et al. (2013) MicroRNA-143 suppresses gastric cancer cell growth and induces apoptosis by targeting COX-2. World J Gastroenterol 19: 7758–7765. doi: 10.3748/wjg.v19.i43.7758
    [33] FDA expands approved use of Kalydeco to treat additional mutations of cystic fibrosis, FDA News Release, 2017. Available from: https://www.fda.gov/NewsEvents/Newsroom/PressAnnouncements/ucm559212.htm
    [34] Fajac I, De Boeck K (2017) New horizons for cystic fibrosis treatment. Pharmacol Ther 170: 205–211. doi: 10.1016/j.pharmthera.2016.11.009
    [35] Clancy JP, Rowe SM, Accurso FJ, et al. (2012) Results of a phase IIa study of VX-809, an investigational CFTR corrector compound, in subjects with cystic fibrosis homozygous for the F508del-CFTR mutation. Thorax 67: 12–18. doi: 10.1136/thoraxjnl-2011-200393
    [36] Boyle MP, Bell SC, Konstan MW, et al. (2014) A CFTR corrector (lumacaftor) and a CFTR potentiator (ivacaftor) for treatment of patients with cystic fibrosis who have a phe508del CFTR mutation: a phase 2 randomised controlled trial. Lancet Respir Med 2: 527–538. doi: 10.1016/S2213-2600(14)70132-8
    [37] Wainwright CE, Elborn JS, Ramsey BW, et al. (2015) Lumacaftor-Ivacaftor in Patients with Cystic Fibrosis Homozygous for Phe508del CFTR. N Engl J Med 373: 220–231. doi: 10.1056/NEJMoa1409547
    [38] Zarrilli F, Amato F, Morgillo CM, et al. (2017) Peptide Nucleic Acids as miRNA Target Protectors for the Treatment of Cystic Fibrosis. Molecules 22: 1144. doi: 10.3390/molecules22071144
    [39] Cantin AM (2016) Cystic Fibrosis Transmembrane Conductance Regulator. Implications in Cystic Fibrosis and Chronic Obstructive Pulmonary Disease. Ann Am Thorac Soc 13 Suppl 2: S150–S155.
    [40] Solomon GM, Raju SV, Dransfield MT, et al. (2016) Therapeutic Approaches to Acquired Cystic Fibrosis Transmembrane Conductance Regulator Dysfunction in Chronic Bronchitis. Ann Am Thorac Soc 13 Suppl 2: S169–S176.
  • This article has been cited by:

    1. Selenne Banuelos, Hayriye Gulbudak, Mary Ann Horn, Qimin Huang, Aadrita Nandi, Hwayeon Ryu, Rebecca Segal, 2021, Chapter 6, 978-3-030-57128-3, 111, 10.1007/978-3-030-57129-0_6
    2. Danfeng Pang, Yanni Xiao, Xiao-Qiang Zhao, A cross-infection model with diffusive environmental bacteria, 2022, 505, 0022247X, 125637, 10.1016/j.jmaa.2021.125637
    3. Iqbal Ahmad, Hesham A. Malak, Hussein H. Abulreesh, Environmental antimicrobial resistance and its drivers: a potential threat to public health, 2021, 27, 22137165, 101, 10.1016/j.jgar.2021.08.001
    4. Mohammad Irfan, Alhomidi Almotiri, Zeyad Abdullah AlZeyadi, Antimicrobial Resistance and Its Drivers—A Review, 2022, 11, 2079-6382, 1362, 10.3390/antibiotics11101362
    5. Danfeng Pang, Yanni Xiao, Xiao-Qiang Zhao, A cross-infection model with diffusion and incubation period, 2022, 27, 1531-3492, 6269, 10.3934/dcdsb.2021316
    6. Suttikiat Changruenngam, Charin Modchang, Dominique J. Bicout, Modelling of the transmission dynamics of carbapenem-resistant Klebsiella pneumoniae in hospitals and design of control strategies, 2022, 12, 2045-2322, 10.1038/s41598-022-07728-w
    7. Cara Jill Sulyok, Lindsey Fox, Hannah Ritchie, Cristina Lanzas, Suzanne Lenhart, Judy Day, Mathematically modeling the effect of touch frequency on the environmental transmission of Clostridioides difficile in healthcare settings, 2021, 340, 00255564, 108666, 10.1016/j.mbs.2021.108666
    8. Kamal Raj Acharya, Jhoana P Romero-Leiton, Elizabeth Jane Parmley, Bouchra Nasri, Identification of the elements of models of antimicrobial resistance of bacteria for assessing their usefulness and usability in One Health decision making: a protocol for scoping review, 2023, 13, 2044-6055, e069022, 10.1136/bmjopen-2022-069022
    9. Shyni Unni Kumaran, Lavanya Rajagopal, Manavaalan Gunasekaran, Sensitivity assessment of optimal control strategies and cost-effectiveness analysis of a novel Candida Auris environmental transmission model in intensive care facilities, 2024, 595, 00225193, 111931, 10.1016/j.jtbi.2024.111931
    10. Reuben Iortyer Gweryina, Muhammadu Yahaya Kura, Timothy Terfa Ashezua, Optimal control model for the infectiology of staphylococcus aureus with dual transmission pathways, 2024, 14, 26667207, 100364, 10.1016/j.rico.2023.100364
    11. Bilal Aslam, Rubab Asghar, Saima Muzammil, Muhammad Shafique, Abu Baker Siddique, Mohsin Khurshid, Muhammad Ijaz, Muhammad Hidayat Rasool, Tamoor Hamid Chaudhry, Afreenish Aamir, Zulqarnain Baloch, AMR and Sustainable Development Goals: at a crossroads, 2024, 20, 1744-8603, 10.1186/s12992-024-01046-8
    12. Saeed Shafait, Shazia Nisar, Kinza Nawabi, Hassan Riaz, Ayesha Masood, Mehtab Ahmed, Comparison of Frequency of Pathogenic Micro-Organisms Causing Bloodstream Infections in Patients Admitted at Tertiary Care Hospital Rawalpindi, 2024, 2790-9352, 115, 10.54393/pjhs.v5i07.1435
    13. Danfeng Pang, Yanni Xiao, Xiao-Qiang Zhao, A periodic reaction-diffusion model of hospital infection with crowding effects, 2024, 539, 0022247X, 128487, 10.1016/j.jmaa.2024.128487
    14. Saleh Mohammed Al-maaqar, Abdulaziz Radhi S. Al Johni, Nasser A. Al-Tayyar, Jafar Abdullah Alhamad, Abdullah A. Khan Ghyathuddin, Wael A. Alsubhi, Ammar AL-Farga, Nahid Kamal Eldin, Hala Mohammad Marouf, Mohsen A. Khormi, Bdellovibrio bacteriovorus: a “magic weapon” against bacterial pathogens, 2025, 75, 1869-2044, 10.1186/s13213-025-01794-x
  • Reader Comments
  • © 2018 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5014) PDF downloads(1077) Cited by(10)

Figures and Tables

Figures(2)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog