A Stochastic Model of Nosocomial Epidemics in Hospital Intensive Care Units

A stochastic differential equation model is developed to portray an epidemic of antibiotic resistant infections in a hospital intensive care unit. The dynamical behaviors of the solutions of this model are analyzed and numerical simulations are given to illustrate the dynamics of the solutions.


Introduction
Nosocomial (hospital-acquired) infections caused by antibiotic resistant bacteria are a major global public health problem.Numerous factors contribute to the emergence and spread of these bacterial infections in hospital settings.To fully understand the impact of different factors, various mathematical models have been developed [3,4,6,8,13,14].Many of these models are formulated as ordinary differential equations (ODEs).These models divide patients and health care workers (HCW) into different compartments, such as infected or uninfected patients and contaminated or uncontaminated HCW.The rate of change of each compartment is described by an ODE under the assumption that the interactions among different groups are deterministic.These models have been applied to the population-level analysis of the spread of nosocomial epidemics [4,5,14].
One drawback of ODE models is that they cannot directly reflect randomness in epidemic events.For instance, the parameters in ODE models should be viewed as averages.Consequently, ODE models can only describe average behavior.Hence there is a need to formulate randomness more precisely in nosocomial models.This is especially true for nosocomial models in hospital subunits, such as an intensive care unit (ICU), which are usually very small, and where randomness may have large influence.Thus, continuous-time Markov chain models (CTMC) and individual based models (IBM) have been used to model nosocomial Corresponding author.Email: wangmin@rowan.edu

Derivation of the stochastic model
In [13] an antibiotic resistant infection epidemic in an ICU is treated with both IBM and ODE approaches.The epidemic population was divided into the following seven compartments: • uninfected patients (P U ); • patients infected with a nonresistant bacterial strain not on antibiotics specific to this strain (P N off ); • patients infected with a nonresistant bacterial strain on antibiotics specific to this strain (P N on ); • patients infected with the resistant strain of this bacteria (P R ); • uncontaminated HCW (H U ); • HCW contaminated with the nonresistant bacterial strain (H N ); • and HCW contaminated with the resistant strain (H R ).
Then a system involving seven ODEs was derived; see [13, Eq. ( 1)-( 7)].That system was further simplified based on certain assumptions and reduced to a model involving three ODEs; the reader is referred to [13, Section 3] for the details.
In this section, we first use the idea similar to [13] to build an ODE model.To simplify the problem, we combine P N off and P N on defined above as a single compartment: patients infected with a nonresistant bacterial strain (P N ).All other compartments remain the same.The transmission among the compartments is described by  Based on Figure 2.1, a system involving six ODEs is derived as follows: where , are all positive constant parameters.The model assumptions and meanings of parameters may be summarized as follows (the reader is referred to [13] for the details).
1.There are N P patients and N H HCWs in the ICU.The units of time are days.
2. All patients who exit the ICU are immediately replaced by uninfected patients.The exit rate of an uninfected patient is not specified as another uninfected patient replaces such a patient immediately.
3. T N is the average time of the length of stay (LOS) of a patient infected with the nonresistant strain.It is assumed this value is additional to any time already spent in the ICU as an uninfected patient.1/T N is interpreted as the exit rate of a patient infected with the nonresistant strain.
4. T R is the average time of the LOS of a patient infected with the resistant strain.It is assumed this value is additional to any time already spent in the ICU as a patient uninfected plus time spent infected with the nonresistant strain.1/T R is interpreted as the exit rate of a patient infected with the resistant strain.
5. T V is the average time (in days) of a patient-HCW visit.N H /(N P T V ) is the average number of visits received by a patient per day.
6.During each visit of a patient by an uncontaminated HCW there is a probability ω N of HCW contamination by a patient infected with a nonresistant strain and a probability ω R of HCW contamination by a patient infected with a resistant strain.A contaminated HCW remains contaminated only for the subsequent next visit to a patient.

7.
During each visit of an uninfected patient by a HCW contaminated with the nonresistant strain, there is a probability π N that the patient is infected with the nonresistant strain.
During each visit of a patient infected with the nonresistant strain by a HCW contaminated with the resistant strain, there is a probability π R that the patient is infected with the resistant strain.
Implicitly assumed in the model is that patients infected with the non-resistant strain, P N (t), are prescribed antibiotics and the resistant strain can only infect patients on antibiotics.The justification for this assumption is that antibiotics provide a favorable environment for the resistant strain to infect a patient and, in the absence of antibiotics, the resistant strain cannot establish an infection.Additionally, within-host mutation from the non-resistant strain to the resistant strain is assumed to be sufficiently rare, so that infected patients (on antibiotics) only become resistant through exposure to the circulating resistant strain via contaminated HCW.Since the populations of patients and HCW remain constant in time, equations (2.1) and (2.4) can be eliminated, and the following system is obtained: ) Further since the time-scale of patient-HCW visits is much faster than patient turnover, it is assumed in [13] that (H1) the HCW compartments are in a quasi-steady state, i.e., Then, by (2.5), (2.6), and (H1), Hence, System (2.7)-(2.10) is simplified to two equations: ) It is clear that (H1) plays an important role to derive System (2.12), (2.13).However, (H1) is a simplification of large scale uncertainty in random events.Thus, it is more realistic to allow H N and H R to vary from their quasi-steady states.We therefore weaken assumption (H1) by introducing random perturbations into (2.11): where Ḃ is a white noise, and σ 1 ≥ 0 and σ 2 ≥ 0 are the intensities of the noise.Thus, by (2.7), (2.8), and (2.14), we obtain the following SDE model: ) where .17) .18) and B(t) is a standard Brownian motion.
Remark 2.1.The solutions P N and P R of System (2.15), (2.16) are two continuous stochastic processes defined on a probability space (Ω, A, P ).Let T be an interval in time.Then both The normal convention is that the variable ω is often suppressed; see for example [1] for the details.In the sequel, we will use P N (t) and P R (t) to denote the solutions of System (2.15), (2.16).ω is only included as needed.
Other parameters in the model, such as transmission probabilities, contact or removal rates, may also vary randomly.However, we leave this for future research and only consider random perturbations of the quasi-steady states for contaminated HCW.The inclusion of random perturbations in these terms reflects the inherent uncertainty in the quasi-steady state assumption.The fast dynamics associated with HCW contamination and decontamination which allow for the quasi-steady state assumption, also can lead to high-frequency noise when stochasticity is included.Indeed, in multiple types of stochastic models of nosocomial infections, HCW contamination is highly variable on small intervals of time [11,12,14].Thus, it seems that equations in (2.14) are an effective way to introduce random perturbations, while also allowing for analysis of a reduced model afforded by the quasi-steady state assumption.
Proof.Since f 1 , f 2 , g 1 , and g 2 are locally Lipschitz continuous, for any initial value (P ++ , there exists a unique local solution (P N (t), P R (t)) on [0, τ e ), where τ e is the explosion time.Let k 0 > 0 be large enough so that P For any k ≥ k 0 , define the stopping time by and τ k is increasing.Let τ ∞ = lim k→∞ τ k .Then we have τ ∞ ≤ τ e .If we can prove that τ ∞ = ∞ a.s., so is τ e .Hence System (2.15), (2.16) has a unique global solution which remains in R 2  ++ a.s.Now we will prove τ ∞ = ∞ a.s.Assume the contrary.There exist T > 0 and ε ∈ (0, 1) such that P {τ ∞ ≤ T} > ε.
Hence there exists k 1 ≥ k 0 such that Let Ω k = {τ k ≤ T}.For any ω ∈ Ω k , at least one of the follows must hold: It is easy to see that V 1 (x, y) ≥ 0 on R 2 ++ .Furthermore, (3.1) implies that where F and G are defined by with f 1 , f 2 , g 1 , and g 2 defined by (2.17)-(2.20).
After some computation, we have where By using the fact that x ≤ 2(x + 1 − ln(x)), we have where C 2 = max{C 0 , 2C 1 }.Then for any t 1 ≤ T, we have where a ∧ b := min{a, b}.Hence, Taking the expectation on both sides, we have Then by the Gronwall inequality, we have . By (3.2) and (3.3), for k ≥ k 1 , we have Next, we consider the stability of the trivial solution (0, 0) of (2.15), (2.16).In the sequel, let P N (t) = P N (t; R ) be the solutions of (2.15) and (2.16) starting from the initial value (P for any (P Proof.For any initial value (P ++ , by Theorem 3.1, System (2.15), (2.16) has a unique solution (P N , P R ) that remains in R ++ with probability 1.By (2.15), (2.16), we have where Then by (3.4)-(3.7)and It ô's formula, • T N is the average time of LOS of a patient infected with the nonresistant strain; • N H /(N P T V ) is the average number of visits received by a patient per day; • ω N is the probability of an uncontaminated HCW contamination by a patient infected with a nonresistant strain during each visit; • π N is the probability that an uninfected patient is infected with the nonresistant strain during each visit by a HCW contaminated with the nonresistant strain.
Therefore by (3.4), R 0 relates to these 4 factors.Indeed, R 0 can be interpreted as the number of contacts between HCW and a patient infected with nonresistant strain during the infectious period, T N N H /(N P T V ), multiplied by the probability HCW contamination and subsequent transmission during visit with next patient in wholly susceptible population, ω N π N .Furthermore, Theorem 3.2 enlightens us that one way to control the spread of the epidemic in an ICU is to reduce the values of these four factors, especially the probabilities ω N and π N .
Remark 3.4.By Theorems 3.1 and 3.2 and their proofs, we know that the noise intensities σ 1 and σ 2 do not affect the positivity of solutions and the almost surely exponential stability of the trivial solution of (2.15), (2.16) when R 0 < 1.

Numerical simulations
In this section, we use numerical simulations to verify the results obtained in Section 3. The values of the parameters of (2.15), (2.16) and the initial values in the following examples are chosen for illustration purpose and are not from actual ICU data.
Example 4.1.Example 4.2.We also investigate the effect of varying the intensities σ 1 and σ 2 when R 0 < 1 in Figure 4.3.By Theorems 3.1 and 3.2, the positivity and convergence to the disease-free equilibrium are not affected, but observe in the figures that the amplitude of random oscillations in the transient dynamics increases with σ 1 and σ 2 .We have also carried out numerical simulations when R 0 > 1. Results show that the solutions will approach certain stationary solutions, which depend on another reproduction number R 01 .This work will appear in another paper.

Conclusions
In this paper we derive a SDE model of an antibiotic resistant infection epidemic in a hospital ICU.The positivity of solutions and almost surely exponential stability of the trivial solution are proved, when the reproduction number R 0 < 1.These results are then illustrated by numerical simulations.These behaviors are consistent with the results obtained in [13].This concurrence suggests that we may control the spread of nosocomial epidemics by adjusting parameters such as π N , π R (the probabilities of patient infection due to patient-HCW visits) or ω N , ω R (the probabilities of HCW contamination due to patient-HCW visits).Our work demonstrates that SDE models are useful for investigating nosocomial epidemics, since SDE models can better reflect the uncertainty and randomness that occur in actual ICU settings.

Figure 4 . 1 :Figure 4 . 2 :
Figure 4.1: (a) Sample path of P N (t) (blue) and P R (t) (green) with the initial condition (19, 11) when R 0 < 1.(b) Plot of the random perturbations about the HCW quasi-steady state for this example, i.e.H N (t) (blue) and H R (t) (green) in equation 2.14.The parameter values are stated in the text for Example 4.1.
5, and σ 2 = 0.2.It is easy to see that R 0 = 0.96 < 1.Then by Theorems 3.1 and 3.2, all the solutions of System (2.15), (2.16) remain in R 2 ++ and the trivial solution is almost surely exponentially stable in probability.Simulation is performed by Matlab.A sample path with the initial condition (19, 11) is given in Figure 4.1.Phase portraits of 12 solutions with different initial conditions are given in Figure 4.2.