Transmission thresholds for the spread of infections in healthcare facilities

TL;DR

This paper studies how infections spread in healthcare facilities and introduces a new formula to calculate the basic reproduction number R0, which helps determine if an outbreak can be sustained and how interventions can reduce it.

Contribution

The paper provides a novel formula for R0 in healthcare facilities and shows how statistical properties like patient stay duration affect outbreak potential.

Findings

R0 was estimated at 1.24 for carbapenemase-producing Enterobacteriaceae in long-term acute-care hospitals before intervention.

Weekly surveillance with transmission reduction could lower R0 to 0.85, potentially preventing outbreaks.

Reducing the mean length of stay alone may not reduce R0 if the variance-to-mean ratio remains unchanged.

Abstract

Some infections may be sustained in the human population by persistent transmission among patients in healthcare facilities, including patients colonized with multi-drug-resistant organisms posing a major health threat. A nuanced understanding of facility characteristics that contribute to crossing a threshold for self-sustaining outbreak potential may be crucial to designing efficient interventions for lowering regional disease burden and preventing high-risk infections. Using a mathematical model, we define the facility basic reproduction number R0, where a single facility can sustain an outbreak without ongoing importation under the threshold condition R0 > 1. We define R0 for a general model with heterogeneous patient susceptibility and transmissibility and with generic length-of-stay assumptions, and we provide a software package for numerical calculation of user-defined examples. We estimate R0 using published data for carbapenemase–producing Enterobacteriaceae (CPE) in long-term acute-care hospitals (LTACHs) and the effects of interventions on R0, including surveillance, pathogen reduction treatments, and length-of-stay reduction. In a simple model, R0 is directly proportional to the sum of the mean and variance-to-mean ratio of the length-of-stay distribution. In intervention models, R0 depends on the moment-generating function of the length-of-stay distribution. From the CPE data, we estimated R0 = 1.24 (95% CI: 1.04, 1.45) prior to intervention. Weekly surveillance with 50% transmission reduction of detected patients alone could have reduced R0 to 0.85 (0.72, 0.98), with additional reduction if detected patients could be decolonized. Reducing the mean length of stay does not necessarily reduce R0 if the variance-to-mean ratio is not also reduced. We conclude that R0 > 1 conditions plausibly exist in LTACHs, where CPE outbreaks could be sustained by patients who acquire colonization and subsequently transmit to other patients during the same hospital stay. Our findings illuminate epidemiological mechanisms producing those conditions and their relationship to interventions that could efficiently improve population health.

The spread of infectious diseases between patients admitted to the same healthcare facility can cause large outbreaks of dangerous organisms such as bacteria resistant to available antibiotic treatments. One way to quantify facility transmission is with the basic reproduction number, R0. When a facility patient acquires the organism, R0 is the average number of onward transmissions from that patient, when all other patients are susceptible. There is an important threshold effect at R0=1: when R0<1, a single introduction is unlikely to set off a long-lasting outbreak, but if R0>1, a continual chain of transmissions stemming from one introduction is possible. We provide a novel formula for R0 specific to healthcare facilities and exhibiting this threshold effect. We show that R0>1 likely occurred in facilities with high rates of drug-resistant bacterial infections, and we provide nuanced insights into intervention efforts to push R0<1 in those facilities. Our formulas reveal that R0 depends on statistical properties of healthcare facilities that are underappreciated for their epidemiological importance, such as the variance of patient stay duration. Our publicly available software allows other researchers to analyze models with customizable facility and disease properties and design intervention strategies for efficient reduction of facility outbreak risk.

Figures

Relationships between facility 𝐑0 and equilibrium clinical infection incidence with constant importation.The curve was generated using formulas for Model 4, at values for the surveillance test frequency ranging from 0 to once per week and no pathogen reduction effect of the intervention. Other model parameters were set at the assumed and estimated values described in Section 3.3, fit to data from CPE in LTACHs. The points (circles) represent average surveillance intervals of, from left to right, one week, two weeks, one month, and no surveillance. Facility 𝐑0 was calculated using the formula in Section 3.2.4 with no decolonization, and we calculated infections per admission as the product of the mean length of stay, the per capita clinical detection rate, and the equilibrium prevalence of colonized, non-clinically detected patients, assuming a constant importation rates of 0.01% (solid), 0.1% (dashed), 1% (dotted), and 5% (dash-dot).

Effect of length-of-stay interventions on facility basic reproduction number.Calculations of the facility reproduction number, 𝐑0, for Model 3 with the following fixed values calibrated to pre-intervention CRE data in long-term acute care hospitals (Table 3): transmission rate β= 0.051 per day, clinical detection rate δc= 0.0084 per day, colonization clearance rate γ=1/ 387 per day, and contact precaution effectiveness ϵ=0.5. Length of stay was governed by a 4-parameter mixed exponential–gamma distribution, and specific parameters were varied from the baseline values (Table 3) in a direction that decreased the overall mean length of stay from left-to-right in each of the sub-figures, while holding all other parameters constant. Upper-left: both the exponential rate parameter 𝐫x and the gamma rate parameter 𝐫g were simultaneously increased by the same factor from their baseline values. Upper-right: the exponential rate parameter 𝐫x was increased from its baseline value. Lower-left: the gamma rate parameter 𝐫g was increased from its baseline value. Lower-right: the fraction of patients following the exponential distribution 𝐩x was decreased from its baseline value.

Scatterplot of facility basic reproduction number vs. length of stay (LOS) distribution statistics.Calculations of the facility reproduction number, 𝐑0, for Model 3 with the following fixed values calibrated to pre-intervention CRE data in long-term acute care hospitals (Table 3): transmission rate β=0.051 per day, clinical detection rate δc=0.0084 per day, colonization clearance rate γ=1/ 387 per day, and contact precaution effectiveness ϵ=0.5. The dots represent 𝐑0 results for 1000 sets of values for each of the four parameters governing the mixed exponential–gamma length of stay distribution, drawn from four independent uniform distributions ranging from 60% to 140% of the calibrated values in Table 3. Each subplot displays the same 1000 𝐑0 results, plotted against different statistics of the length of stay distribution for each set of length-of-stay parameters.

Equations

Full text

1 Introduction

Infections caused by multi-drug resistant organisms (MDROs) have limited available treatments and pose a significant health threat [1]. Some MDROs have shown the ability to disseminate rapidly among patients within and between certain healthcare facilities [2], suggesting that efforts to reduce patient-to-patient transmission via inter-facility coordination [3] and/or interventions targeted at high-risk settings [4,5] could provide substantial population benefit. Results from simulations of regional outbreaks and interventions suggest that, in some scenarios, efforts to prevent transmission such as active surveillance for MDRO carriage among high-risk patients could drastically reduce regional vulnerability to explosive outbreaks that would lead to endemic presence of the MDRO [4]. These results point to the possibility that rates of MDRO transmission between healthcare patients could be close to the threshold required to sustain long-lasting chains of person-to-person transmission that are characteristic of successful invasion of a novel transmissible organism in a population.

The quantitative theory of transmission thresholds has generated a rich body of research. Levels of transmission for an organism in a population relative to a threshold are often characterized by the basic reproduction number, $[eqn]$ , the expected total number of transmissions from a typical infected individual in an unexposed population [6]. If $[eqn]$ is larger than the threshold of 1, the organism can invade and produce a sustained outbreak over many generations of transmission. When above the threshold, the value of $[eqn]$ can determine the theoretical level of transmission-reducing intervention effort, such as vaccination, required to control outbreaks [7].

Some studies have applied transmission threshold and reproduction number theory specifically to organisms transmitted primarily in healthcare settings, usually defining $[eqn]$ as the expected number of transmissions from an admitted carrier of an organism during a single stay in a healthcare facility under specific assumptions for a particular pathogen [8–11]. Others have defined this quantity as the admission reproduction number, $[eqn]$ [12–17], or ward reproduction ratio, $[eqn]$ [18], noting that the quantity is less than $[eqn]$ when a single facility or ward stay may represent only a portion of the patient’s potential to transmit during the infectious period.

While the reproduction number in a healthcare facility is partly determined by properties of the organism and its typical course of colonization and infection in humans, it can also be affected by actions the facility takes to prevent patient-to-patient transmission. Of particular interest are scenarios where a facility’s $[eqn]$ for an MDRO under the facility’s standard transmission control efforts, but additional actions, such as active surveillance and contact precautions for detected carriers of an MDRO [19] or decolonization via pathogen reduction approaches [5,20], could potentially move $[eqn]$ below threshold.

Our aims for this work were first to define $[eqn]$ for a broad class of mathematical models of a single healthcare facility and demonstrate its utility as a threshold for facility outbreak potential, and second to demonstrate our model’s utility in assessing intervention components that might push a vulnerable facility below threshold. For this second aim, we use data from a bundled intervention that successfully reduced infections with carbapenemase-producing Enterbacteriaceae (CPE) in long-term acute care hospitals (LTACHs) [21]. Due to their composition of acutely ill patients with long length of stay, LTACHs are high-risk settings for infections with CPE and other high-priority organisms [22]. Efforts to interrupt transmission could provide substantial benefit to LTACHs and other facilities linked by patient exchange [4]. Therefore, understanding the role of individual intervention components for reducing transmission could be crucial to achieving maximum efficiency in regional outbreak prevention.

The overall outline of our study is as follows. First, we formulate differential equations describing patient states in a healthcare facility pertaining to colonization with and susceptibility to an infectious organism as well as admission and discharge. Next, we derive formulas for $[eqn]$ , the expected transmissions from a patient who acquires colonization during a facility stay, in an otherwise fully susceptible facility population and with no ongoing importation. We also derive a formula for the equilibrium facility patient state distribution when the admission state distribution (including the importation rate of colonized patients) is constant. We use this equilibrium equation to estimate model parameters from CPE data from LTACHs in which constant importation and steady state facility prevalence were observed, and we estimate $[eqn]$ using those parameter estimates and the formulas we had derived. Finally, we investigate nuances on the effects of interventions on reducing $[eqn]$ , to provide insights on how similar LTACHs without ongoing CPE importation could substantially reduce invasion risk by keeping $[eqn]$ .

2 Methods

2.1 General model description

First, we defined a general class of mathematical models describing states of patients in a healthcare facility. We split the patients into different states enumerated by integer subscripts $[eqn]$ , and let $[eqn]$ be the probability that a patient is alive, not discharged, and in state $[eqn]$ at time $[eqn]$ after admission. Defining the column vector of probabilities $[eqn]$ containing each $[eqn]$ , we define the system of differential equations

$$[eqn]$$

Here, $[eqn]$ contains state-to-state transition rates and may also contain constant, state-specific facility removal rates by death or live discharge. The hazard function $[eqn]$ is a facility removal rate at time $[eqn]$ of the stay, where the same function applies to every state. The state-specific removal rates and/or the hazard function can be used to calibrate a model to an overall observed patient length of stay distribution. The vector $[eqn]$ is the distribution of states at admission. For this system form, we derive the formula for the equilibrium cross-sectional state distribution of facility patients, under the assumption that patients are continually admitted to the facility at a constant rate with an unchanging admission state probability distribution.

Next, for modeling patient states pertaining to a particular organism that can colonize or infect patients in the facility, we assume the patient states $[eqn]$ are comprised of $[eqn]$ susceptible states $[eqn]$ and $[eqn]$ colonized states $[eqn]$ . The colonized states represent patients carrying the modeled organism. The $[eqn]$ different colonized states could represent any epidemiologically relevant differences between different groups of colonized patients, for example, states of infection at different body sites, different levels of transmissibility to other patients, different hazards of death or discharge, etc. The susceptible states represent patients not colonized by the organism. The $[eqn]$ different susceptible states could be distinguished by different risks of acquiring colonization per level of exposure and/or by other relevant differences.

We use the following notation:

$$[eqn]$$ $$[eqn]$$

In the initial state vector $[eqn]$ describing the admission state distribution, $[eqn]$ is the overall probability of being in a colonized state at admission, $[eqn]$ is the admission state distribution among the susceptible state admissions, with $[eqn]$ , and $[eqn]$ is the admission state distribution among the colonized state admissions, with $[eqn]$ .

In the matrix $[eqn]$ , the terms $[eqn]$ in the diagonal elements are the rates of removal from the facility for patients in each state, and the other terms describe transitions from one state to another. The upper-left $[eqn]$ portion of $[eqn]$ contains elements describing the transitions out of and between the susceptible states. The $[eqn]$ terms are the rates of transitioning from one susceptible state to another, with $[eqn]$ if $[eqn]$ and $[eqn]$ if $[eqn]$ . The $[eqn]$ terms are the rates of transitioning from each susceptible state to a colonized state (rates of acquisition). The $[eqn]$ coefficients describe relative susceptibility to acquisition of patients in each susceptible state, and $[eqn]$ is a baseline acquisition rate that is assumed to depend on the prevalence of patients in the colonized states, who can transmit to the susceptible patients, as follows:

$$[eqn]$$

where the $[eqn]$ coefficients define the relative transmissibility of patients in the different colonized states.

The lower-left $[eqn]$ portion of $[eqn]$ contains rates of transitioning into each colonized state from each susceptible state, with $[eqn]$ . The lower-right $[eqn]$ portion of $[eqn]$ contains elements describing the transitions out of and between the colonized states. For the $[eqn]$ elements in the diagonal we have $[eqn]$ , where the $[eqn]$ terms comprising the upper right $[eqn]$ portion of $[eqn]$ are the rates of transitioning from colonized to susceptible, i.e., the rates of clearing colonization.

For this system, we derive a formula for the facility basic reproduction number $[eqn]$ for given values of the importation state distribution, state transition rates, and relative transmissibility parameters defined above. We derive our equation to match the “basic reproduction ratio” formula in Diekmann et al [6]. Their definition depends on quantifying an individual’s state just prior to acquiring colonization and the expected infectivity of an individual who was in a given state just prior to acquisition at a given prior time, toward a susceptible individual currently in a given state. For our model, the state of a facility patient is defined as the susceptible compartment the patient is in and the time since admission. Hence the $[eqn]$ formula involves both summation over discrete patient states and continuous-time integration of functions involving the time since admission. Details are provided in S1 Text.

We also derive formulas for the equilibrium cross-sectional patient state distribution when the admission rate of patients in each state is constant, as described further in subsequent sections. The model is not designed to represent transient dynamics of a facility patient population starting from a non-equilibrium state distribution. In full generality, such a model requires representing patient states as a multivariate function $[eqn]$ that tracks the number of facility patients in each state who were admitted $[eqn]$ time units ago, at time $[eqn]$ after the start of the dynamic simulation. While studying transient dynamics is outside the scope of this study, we provide an example of a model formulation of this type in S1 Text.

2.2 Example model descriptions

We use our general formula for $[eqn]$ to derive $[eqn]$ formulae for a series of specific model examples.

2.2.1 Model 1: Simple susceptible–colonized model.

We specified a simple version of our general disease model with only two states and one direction of transition between states. The state probability vector $[eqn]$ consists of one susceptible state $[eqn]$ that we rename $[eqn]$ and one colonized state $[eqn]$ that we rename $[eqn]$ . The admission state probability vector in the above notation contains $[eqn]$ and $[eqn]$ :

$$[eqn]$$

We assume that susceptible inpatients acquire colonization at rate $[eqn]$ , and there are no other state transitions or removals other than removal modeled by the hazard function $[eqn]$ . Thus, in the general model notation, we have $[eqn]$ . The state transition matrix $[eqn]$ is then

$$[eqn]$$

For transmission, there is only one colonized state from which transmissions can occur, so there is a single transmission parameter ( $[eqn]$ that we rename $[eqn]$ , and the acquisition rate $[eqn]$ . The time-of-stay-dependent removal hazard $[eqn]$ is left general. The equivalent system in non-matrix form is as follows:

$$[eqn]$$ $$[eqn]$$

2.2.2 Model 2: Clearance of colonization.

To extend Model 1, we add a non-zero clearance rate $[eqn]$ of moving from colonized back to susceptible. In the general model notation, we have $[eqn]$ . The state transmission matrix $[eqn]$ is then

$$[eqn]$$

In non-matrix equation form, we have:

$$[eqn]$$ $$[eqn]$$

2.2.3 Model 3: Clinical detection with contact precautions.

For this model, we extend Model 2 to include two colonized states $[eqn]$ that we rename $[eqn]$ and $[eqn]$ . These two colonized states represent, respectively, patients with no prior detection of their colonization and patients with a prior clinical detection. This type of model is generally used to represent the transition from asymptomatic colonization to the onset of symptomatic, invasive infection that triggers a clinical test for the causative organism. For the admission state probability vector in the above notation, we assume $[eqn]$ , $[eqn]$ , and $[eqn]$ , i.e., all colonized patients are undetected at admission:

$$[eqn]$$

For state transitions among inpatients, we assume that susceptible patients acquiring colonization (at rate $[eqn]$ ) move into the undetected colonized state. Undetected colonized patients move to the susceptible state at clearance rate $[eqn]$ and move to the clinically detected state at clinical detection rate $[eqn]$ , which represent a progression rate from asymptomatic colonization to confirmed clinical infection. Clinically detected patients are assumed to remain in that state until the end of the stay. Thus, in the general model notation, we have:

$$[eqn]$$

The state transmission matrix $[eqn]$ is then

$$[eqn]$$

For the relative transmission rates from patients in the two colonized states, we assume that clinically detected patients transmit at a rate reduced according to the effectiveness $[eqn]$ of contact precautions: $[eqn]$ . Then the acquisition rate $[eqn]$

In non-matrix equation form, Model 3 is written as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

2.2.4 Model 4: Active surveillance and decolonization.

We extend Model 3 to include two simultaneous interventions: active surveillance to detect patients with colonization who have not been detected by clinical test, and decolonization of those detected patients via administration of a decolonizing drug or pathogen reduction agent. Patients detected by active surveillance are assumed to be placed under contact precautions in addition to receiving the pathogen reduction treatment, which increases the rate of clearance above the natural clearance rate.

The state probability vector $[eqn]$ consists of two susceptible states $[eqn]$ that we rename $[eqn]$ and $[eqn]$ and three colonized states $[eqn]$ that we rename $[eqn]$ , and $[eqn]$ . The new states with subscript “sd” represent patients with a positive surveillance test earlier in the facility stay. At admission, we assume that susceptible patients never test positive, and colonized patients test positive with probability $[eqn]$ , which incorporates the probability of receiving an admission surveillance test and the probability that the test avoids a false negative result. For the admission state probability vector in the above notation, we assume $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ :

$$[eqn]$$

For state transitions among inpatients, we assume that undetected colonized patients become surveillance detected at rate $[eqn]$ , representing mid-stay surveillance performed on patients who did not test positive at admission. Surveillance-detected, colonized patients can move to the clinically detected state at the same rate $[eqn]$ as undetected colonized patients. The new susceptible state $[eqn]$ can be reached by surveillance-detected, colonized patients who clear colonization during the stay. Surveillance-positive patients are assumed to remain under contact precautions for their entire remaining stay. For susceptible patients under contact precautions, we assume that those precautions reduced their rate of re-acquisition by the same factor, $[eqn]$ , that they reduce transmissibility. In the general model notation we have the following nonzero elements: $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ . The state transmission matrix $[eqn]$ is then

$$[eqn]$$

For relative transmission rates, we assume that surveillance detected, colonized patients transmit at a rate commensurate with clinically detected patients: $[eqn]$ . Then the acquisition rate $[eqn]$ .

In non-matrix equation form, Model 4 is written as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

2.3 Calibration of Model 3 to pre-intervention CPE in LTACH

As in our prior work [5,23], we calibrated model parameter values to data for CPE in Chicago-area LTACHs [21]. Prior to intervention in those LTACHs, cross-sectional prevalence surveys revealed nearly 50% prevalence of CPE carriers, and there were significant rates of CPE clinical detection (Table 1). The facilities did not routinely perform CPE surveillance before the intervention, thus we used our Model 3 to represent the pre-intervention scenario with no active surveillance. For the facility removal (discharge or death) hazard, we used a function giving rise to a length-of-stay distribution governed by a mixture of an exponential and a gamma distribution, which in prior work we found to provide a parsimonious fit to detailed LTACH length of stay data [4]. The distribution has a total of four parameters: the rate parameter of an exponential distribution, $[eqn]$ , the rate and shape parameters of a gamma distribution, $[eqn]$ and $[eqn]$ , and the portion of patients following the exponential distribution, $[eqn]$ . We calculated the unique values of these four parameters that produce a length-of-stay distribution matching the following statistics reported in the LTACH study [21]: the mean, median, and interquartile range of hospital stay (S1 Text).

Table 1: Data from pre-intervention LTACHs [21].

For the CPE parameters in Model 3, we assumed the CPE clearance rate $[eqn]$ [24], contact precaution effectiveness $[eqn]$ [25], and CPE test sensitivity = 0.85 [26,27]. Because the LTACH importation rate and cross-sectional prevalence of CPE were observed and stable over time, we assumed the Model 3 system was at equilibrium, which means that the acquisition rate, $[eqn]$ , was constant, making Model 3 a linear system. We solved for the value of $[eqn]$ and of the progression rate to clinical detection, $[eqn]$ , that produced a match to the equilibrium cross-sectional CPE carriage prevalence (scaled by the assumed test sensitivity) and CPE clinical detection incidence reported in the source LTACH study (Table 1). Then we calculated the baseline transmission rate, $[eqn]$ , that produced the correct acquisition rate under the calibrated model (Table 2 and S1 Text). Finally, we used the assumed values of $[eqn]$ and $[eqn]$ and the data-calibrated values of $[eqn]$ and $[eqn]$ to calculate $[eqn]$ using the formula we derived for Model 3 (see Results).

Table 2: Model 3 parameters and assumptions.

The value of $[eqn]$ in the importation rate formula is the fraction of patients positive by surveillance test at admission (Table 1). The values of the length of stay parameters were calibrated to match the observed length of stay statistics (Table 1), and the values of $[eqn]$ and $[eqn]$ were calibrated to produce a model equilibrium matching the cross-sectional facility CPE positivity fraction and clinical detection rate data (Table 1).

We accounted for uncertainty arising from certain pre-intervention data that were used for calibration (admission CPE positivity and surveillance adherence fractions), or calibration targets (cross-sectional CPE positivity fraction and clinical detection incidence) in the model. Here, we relied on the reported uncertainty ranges from the source paper [21] (Table 1). Each value with reported uncertainty was derived from a large number of observations [21], justifying the use of a normal distribution approximation to the underlying distribution. We incorporated these assumed parameter distributions into a Latin hypercube sampling scheme [28], which produces random sets of values representative of the multi-variate distribution. We used 1000 combinations of the paper-derived values from the Latin hypercube sampling scheme to create 95% intervals for our results, including $[eqn]$ .

2.4 Application of Model 4 to intervention analysis for CPE in LTACH

We used Model 4 to investigate the potential effects of active surveillance and decolonization for CPE carriage on reducing the facility reproduction number from the pre-intervention estimates. First, we assessed the biweekly surveillance component of the intervention done in the source study [21], which in our prior work we found could largely explain the post-intervention CPE reduction they observed when 50% contact precaution effectiveness was assumed [23]. The surveillance detection rate, $[eqn]$ , was assumed to be the product of the targeted testing rate (once per 14 days), the reported compliance rate (0.954), and the assumed test sensitivity (0.85): $[eqn]$ per day. We tested the effect of increasing the frequency of surveillance from biweekly to weekly and calculated the minimum frequency of surveillance testing required to make $[eqn]$ .

Next, we modeled the effect of combining surveillance with pathogen reduction treatment of surveillance-detected patients, represented by increasing the non-treated rate of clearing colonization among surveillance-detected patients by a factor $[eqn]$ , such that $[eqn]$ . For values of f = 2, 5, 10, 50, 100, and ∞ (immediate decolonization at detection), we calculated $[eqn]$ for weekly and biweekly surveillance, and the surveillance rate required to achieve $[eqn]$ .

2.5 Length of stay intervention analysis

We tested the effects of interventions that would alter one or more of the parameters governing the length of stay distribution on $[eqn]$ , using Model 3 calibrated to LTACH CPE data (Tables 1 and 2) as the starting point. As the overall mean length of stay decreases as either of the rate parameters $[eqn]$ or $[eqn]$ increases, we tested the $[eqn]$ -altering effects of increasing either $[eqn]$ , $[eqn]$ , or both while holding the other parameters constant, to model the effects of efforts to speed up the discharge of LTACH patients. We also tested the effects of altering the $[eqn]$ parameter governing the portion of patients following exponential- or gamma-distributed length of stay in favor of the group with lower mean length of stay.

Finally, we chose random combinations of all 4 parameters governing the length-of-stay distribution and produced scatterplots of the resulting $[eqn]$ (again using Model 3 calibrated to LTACH CPE data) vs. different length-of-stay distribution statistics. We plotted $[eqn]$ against the length of stay mean, standard deviation, variance-to-mean ratio (VMR), and the sum of the mean and VMR, and we calculated Pearson’s correlation coefficient for the correlation between each of those statistics and $[eqn]$ .

3 Results

3.1 General facility basic reproduction number (R0) formula

For the general model described in the Methods section, we derived an expression for the basic reproduction number among facility patients. We express the result using the following components:

$$[eqn]$$

The elements in the matrices $[eqn]$ , $[eqn]$ , and $[eqn]$ were defined within the matrix W in the Methods section. The matrix $[eqn]$ describes transition rates between the susceptible states in the absence of colonized patients, and state-specific removal rates from the facility along the diagonal ( $[eqn]$ terms). The matrix $[eqn]$ contains transition rates from and between the colonized states, including state specific removal rates from the facility along the diagonal ( $[eqn]$ terms). The matrix $[eqn]$ contains the relative transition rates $[eqn]$ from susceptible state $[eqn]$ to colonized state $[eqn]$ upon acquisition. The vector $[eqn]$ contains the transmission rates from each colonized state. The full transition rate from susceptible state susceptible state $[eqn]$ to colonized state $[eqn]$ is $[eqn]$ . The vector $[eqn]$ contains the fractions of admitted, susceptible patients who are in each of the susceptible states at admission. The function $[eqn]$ is the cumulative distribution function of the facility length of stay distribution of patients who are removed via the hazard function $[eqn]$ .

With these elements in hand, we derived the following result for the facility basic reproduction number $[eqn]$ (S1 Text):

$$[eqn]$$

This formula is a dot product of the vector $[eqn]$ , containing the transmissibility per time of each colonized compartment, with a vector of the same length that contains the expected amount of time spent in each colonized compartment between the time of acquiring colonization in the facility and discharge or death.

We provide example models below in which $[eqn]$ can be expressed symbolically, which generally occurs when the eigenvalues of matrices $[eqn]$ and $[eqn]$ can be expressed as formulas in terms of the model parameters. Otherwise, $[eqn]$ can be calculated numerically as follows. In the numerator of the above expression for $[eqn]$ , the vector $[eqn]$ is equivalent to final $[eqn]$ elements of the vector $[eqn]$ , in which

$$[eqn]$$

and $[eqn]$ contains $[eqn]$ rows and $[eqn]$ columns of zeros. We numerically calculate the eigenvalues $[eqn]$ of $[eqn]$ and associated eigenvectors to determine the linear combinations of $[eqn]$ functions that comprise the elements of $[eqn]$ . Then, integrating the product of this vector and $[eqn]$ produces a vector with linear combinations of $[eqn]$ , where $[eqn]$ and $[eqn]$ is the moment generating function of the distribution for which $[eqn]$ is the cumulative distribution function.

For a model with the hazard function $[eqn]$ , the patients’ rate of removal from the facility (by discharge or death) is entirely state-dependent and quantified by the $[eqn]$ values within the matrices. Then the function $[eqn]$ and the basic reproduction number formula (S1 Text) reduces to:

$$[eqn]$$

This result is valid if all patients are guaranteed to eventually be removed from the facility, i.e., there is always a nonzero probability that an inpatient will eventually reach a state $[eqn]$ from which the removal rate $[eqn]$ .

We created an R package, facilityepimath [29], containing a function that numerically calculates $[eqn]$ according to the above formula, for a given set of numerical inputs for $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ (if $[eqn]$ ).

3.2 Example model R0 formulas

3.2.1 Model 1: Simple susceptible–colonized model.

For this model we have:

$$[eqn]$$

and the hazard function $[eqn]$ is left general. Substituting these values in the general formula gives

$$[eqn]$$

In S1 Text we show that $[eqn]$ , and $[eqn]$ , where $[eqn]$ and $[eqn]$ are the mean and variance, respectively, of the distribution for which $[eqn]$ is the cumulative distribution function. Therefore, our result is:

$$[eqn]$$

3.2.2 Model 2: Clearance of colonization.

For this model we have:

$$[eqn]$$

After substitution and some integration (see more detailed steps in S1 Text), the result is:

$$[eqn]$$

In S1 Text demonstrate that, for a real number $[eqn]$ , $[eqn]$ , where $[eqn]$ is the moment-generating function of the distribution for which $[eqn]$ is the cumulative distribution function. For ease of notation, we define the function $[eqn]$ . Our result for this model can then be expressed

$$[eqn]$$

Some commonly used parametric distributions for modeling length of stay provide convenient expressions for the moment-generating function $[eqn]$ within the formula for $[eqn]$ . For example, for length of stay exponentially distributed with rate $[eqn]$ (constant removal hazard), we have $[eqn]$ and $[eqn]$ , which leads to:

$$[eqn]$$

With time to removal gamma distributed with rate parameter $[eqn]$ and shape parameter $[eqn]$ , we have $[eqn]$ and $[eqn]$ , which leads to

$$[eqn]$$

3.2.3 Model 3: Clinical detection with contact precautions.

For this model we have:

$$[eqn]$$

The formula for $[eqn]$ contains $[eqn]$ , which, because $[eqn]$ is a matrix, is the matrix exponential:

$$[eqn]$$

Substituting this and the other expressions in the $[eqn]$ formula, we arrive at the following results after integrating and performing matrix multiplications:

$$[eqn]$$

3.2.4 Model 4: Active surveillance and decolonization.

For this model we have:

$$[eqn]$$

The matrix exponentials $[eqn]$ and $[eqn]$ can be expressed:

$$[eqn]$$ $$[eqn]$$

The values $[eqn]$ and $[eqn]$ are the first and second diagonal elements of the matrix $[eqn]$ . Substituting these and the other expressions in the $[eqn]$ formula, we arrive at the following result after integrating and performing matrix multiplications:

$$[eqn]$$

In the case that $[eqn]$ , the above expression is indeterminate because of division by zero; this corresponds to the matrix $[eqn]$ having a repeated eigenvalue ( $[eqn]$ ) that is degenerate, and the matrix exponential $[eqn]$ includes eigenfunctions $[eqn]$ and $[eqn]$ . In this case, the resulting expression for $[eqn]$ is as follows:

$$[eqn]$$

We also use simplified expressions for other special cases of the modeled intervention. First, with surveillance and contact precautions as the only intervention component, i.e., no decolonization, we have $[eqn]$ and $[eqn]$ simplifies to:

$$[eqn]$$

Second, when decolonization is applied to surveillance-detected patients, and the clearance effect is assumed to be instantaneous ( $[eqn]$ ), we can use the result from Model 2 with $[eqn]$ replaced by the sum of $[eqn]$ and $[eqn]$ . That is,

$$[eqn]$$

3.3 Calibration of Model 3 to pre-intervention CPE in LTACH

To calibrate a model to length of stay data, we previously derived the following general formula for the cumulative distribution function, $[eqn]$ , for the time between admission and facility removal at equilibrium [5,23]:

$$[eqn]$$

The function $[eqn]$ is the cumulative distribution function associated with the time-dependent discharge hazard $[eqn]$ . In models where $[eqn]$ , and facility removal is modeled by state-dependent constant removal terms, the quantity $[eqn]$ is the fraction of patients who are still in the facility at time $[eqn]$ after admission, with admission state distribution $[eqn]$ . Conversely, if there are no facility removal terms within W, $[eqn]$ and the above formula reduces to $[eqn]$ . See our prior study for an example of a model with both nonzero $[eqn]$ (representing live discharge) and nonzero removal elements within W (representing death) [23].

For Model 3, we have $[eqn]$ , and assume a mixed exponential-gamma distribution, with a portion $[eqn]$ of patients following an exponential distribution with rate $[eqn]$ and the rest following a gamma distribution with rate $[eqn]$ and shape parameter $[eqn]$ . The following length of stay statistics reported in the LTACH study were used to estimate the 4 parameters: the mean, median, 25^th^ and 75^th^ percentile of the length of stay distribution. With 4 unknowns fit to those 4 data points (S1 Text), we calculated the unique solution (Table 3).

Table 3: Model calibration and R0 results: CPE LTACH pre-intervention.

For the CPE parameters, under our assumed values for three independently based parameters ( $[eqn]$ , $[eqn]$ per day, $[eqn]$ ), we estimated the pre-intervention values of $[eqn]$ , $[eqn]$ , and $[eqn]$ from CPE data in LTACHs, using Model 3 (Table 3). Our estimate for the baseline transmission rate, $[eqn]$ , was 0.051 per day (95% range 0.043 to 0.059). The estimate of rate of progression to clinical detection, $[eqn]$ , was 0.0084 per day (0.0073 to 0.0097). Using those results and our formula for Model 3 above, we estimated $[eqn]$ = 1.24 (1.04 to 1.45).

3.4 Application of Model 4 to vertical intervention analysis for CPE in LTACH

First, we estimated the effect of the biweekly CPE surveillance component of the intervention reported in the source study. Using the formula for $[eqn]$ from Model 4, with $[eqn]$ (no decolonization effect) and surveillance detection rate $[eqn]$ per day (once per 14 days multiplied by scaling factors for imperfect compliance and imperfect test sensitivity). With these values, and with other assumptions unchanged from the pre-intervention scenario modeled in the previous section, we estimated $[eqn]$ = 0.94 (95% range 0.79 to 1.08). Increasing the frequency of surveillance from biweekly to weekly decreases the estimate to $[eqn]$ = 0.85 (0.72 to 0.98). The minimum frequency of surveillance testing required to make $[eqn]$ was once per 3.2 weeks (1.1 to 23 weeks). We plotted the relationship between $[eqn]$ at different frequencies of surveillance and the number of clinical infections per 1000 admissions at equilibrium, assuming the admission rate of CPE carriers remains constant at 0.01%, 0.1%, 1%, and 5% (Fig 1). When there is continual importation from external sources, reducing $[eqn]$ does not eliminate infections and does not exhibit a noticeable threshold effect when the importation rate is higher. However, when the importation rate is low, the plotted relationship demonstrates an inflection point at the $[eqn]$ threshold, where the infection incidence rises much more steeply with increasing $[eqn]$ when it is above the threshold.

Relationships between facility 𝐑0 and equilibrium clinical infection incidence with constant importation.The curve was generated using formulas for Model 4, at values for the surveillance test frequency ranging from 0 to once per week and no pathogen reduction effect of the intervention. Other model parameters were set at the assumed and estimated values described in Section 3.3, fit to data from CPE in LTACHs. The points (circles) represent average surveillance intervals of, from left to right, one week, two weeks, one month, and no surveillance. Facility 𝐑0 was calculated using the formula in Section 3.2.4 with no decolonization, and we calculated infections per admission as the product of the mean length of stay, the per capita clinical detection rate, and the equilibrium prevalence of colonized, non-clinically detected patients, assuming a constant importation rates of 0.01% (solid), 0.1% (dashed), 1% (dotted), and 5% (dash-dot).

Next, we modeled the effect of combining surveillance with decolonization via pathogen reduction treatment of surveillance-detected patients, represented by increasing the non-treated rate of clearing colonization among of surveillance-detected patients by a factor $[eqn]$ , such that $[eqn]$ . For different values of $[eqn]$ , we calculated $[eqn]$ for weekly and biweekly surveillance, and the surveillance rate required to achieve $[eqn]$ (Table 4). For example, with $[eqn]$ (five-fold increase in clearance rate due to the pathogen reduction treatment), the $[eqn]$ estimate with biweekly surveillance was 0.87 (0.73 to 1.00), an 8% decrease from the estimate of 0.94 with no decolonization effect. The upper bound of the decolonization factor ( $[eqn]$ , or immediate clearance upon detection) achieves $[eqn]$ (0.47 to 0.64) with biweekly surveillance. With immediate clearance, a surveillance interval of once per 10 weeks (6–56 weeks) could potentially achieve $[eqn]$ .

Table 4: R0 and threshold estimates for combined surveillance and decolonization intervention.

3.5 Length of stay intervention analysis

The parameters of Model 3 governing removal (death or discharge) rates of patients were calibrated to LTACH length of stay statistics (Table 1), giving rise to an overall length of stay distribution with mean $[eqn]$ days and standard deviation $[eqn]$ days. With the mixture model we used to govern removal, a fraction $[eqn]$ of admitted patients follow an exponentially distributed length of stay (constant discharge-or-death hazard), with mean and standard deviation $[eqn]$ days. The remaining patients follow gamma-distributed length of stay with shape parameter $[eqn]$ (meaning the removal hazard increases with time of stay), which gives a length of stay distribution with mean $[eqn]$ days and standard deviation $[eqn]$ days. Thus, the patients following the exponential distribution have a slightly higher mean and much higher variance in length of stay than the other group. We refer to the patients following the exponentially distributed length of stay as the high-variance patients and the patients following the gamma-distributed length of stay as the low-variance patients.

We tested the effects of interventions that alter different length of stay parameter values on $[eqn]$ for Model 3 calibrated to pre-intervention LTACH CPE data (Fig 2). We chose four interventions that decrease the mean length of stay from the baseline pre-intervention 33.8 days. First, when increasing the discharge rate of all patients (i.e., increasing the rate parameters $[eqn]$ and $[eqn]$ by the same factor), $[eqn]$ decreases and crosses the $[eqn]$ threshold as the mean length of stay drops below about 26 days. Second, when increasing the discharge rate $[eqn]$ of the high-variance patients only, $[eqn]$ decreases below 1 as the mean length of stay drops to about 28 days. Third, when increasing the discharge rate $[eqn]$ of the low-variance patients, $[eqn]$ does not decrease below 1.2 and in fact begins increasing as the mean length of stay drops below about 29 days and lower. Fourth, when decreasing the fraction $[eqn]$ of patients following the high-variance exponential distribution, $[eqn]$ drops rapidly and falls below 1 as the mean length of stay drops below about 32.6 days. When $[eqn]$ (all patients follow the low-variance gamma distribution with mean 32.0 days), $[eqn]$ 0.88.

Effect of length-of-stay interventions on facility basic reproduction number.Calculations of the facility reproduction number, 𝐑0, for Model 3 with the following fixed values calibrated to pre-intervention CRE data in long-term acute care hospitals (Table 3): transmission rate β= 0.051 per day, clinical detection rate δc= 0.0084 per day, colonization clearance rate γ=1/ 387 per day, and contact precaution effectiveness ϵ=0.5. Length of stay was governed by a 4-parameter mixed exponential–gamma distribution, and specific parameters were varied from the baseline values (Table 3) in a direction that decreased the overall mean length of stay from left-to-right in each of the sub-figures, while holding all other parameters constant. Upper-left: both the exponential rate parameter 𝐫x and the gamma rate parameter 𝐫g were simultaneously increased by the same factor from their baseline values. Upper-right: the exponential rate parameter 𝐫x was increased from its baseline value. Lower-left: the gamma rate parameter 𝐫g was increased from its baseline value. Lower-right: the fraction of patients following the exponential distribution 𝐩x was decreased from its baseline value.

As the results in Fig 2 demonstrate that the mean length of stay did not necessarily correlate with $[eqn]$ , we tested the relationship between $[eqn]$ and other statistics for the length of stay distribution (Fig 3). The dots in the scatterplot represent the results from randomly chosen combinations of the four length-of-stay parameters $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ , each picked from independent uniform distributions ranging from plus-or-minus 40% of the values from the model calibrated to the LTACH data (Table 3). For each of the random set of parameter values, we calculated $[eqn]$ and the mean, standard deviation, variance-to-mean ratio (VMR) and the sum of the mean and VMR. The latter statistic was chosen because of its appearance in our result for the $[eqn]$ formula for Model 1.

Scatterplot of facility basic reproduction number vs. length of stay (LOS) distribution statistics.Calculations of the facility reproduction number, 𝐑0, for Model 3 with the following fixed values calibrated to pre-intervention CRE data in long-term acute care hospitals (Table 3): transmission rate β=0.051 per day, clinical detection rate δc=0.0084 per day, colonization clearance rate γ=1/ 387 per day, and contact precaution effectiveness ϵ=0.5. The dots represent 𝐑0 results for 1000 sets of values for each of the four parameters governing the mixed exponential–gamma length of stay distribution, drawn from four independent uniform distributions ranging from 60% to 140% of the calibrated values in Table 3. Each subplot displays the same 1000 𝐑0 results, plotted against different statistics of the length of stay distribution for each set of length-of-stay parameters.

The length of stay mean produced a weaker correlation with $[eqn]$ (Pearson’s correlation coefficient PCC = 0.857) compared to the correlation between the standard deviation and $[eqn]$ (PCC = 0.986). The VMR produced PCC = 0.860, and the sum of the mean and VMR produced the strongest correlation with $[eqn]$ , with PCC = 0.997 and the relationship appearing nearly linear (Fig 3).

Files with code that produces results in Figs 1–3 and Tables 3 and 4 are available on the facilityepimath GitHub repository at https://github.com/EpiForeSITE/facilityepimath/tree/v0.2.0/inst.

4 Discussion

This study provides a flexible, quantitative framework for calculating the basic reproduction number $[eqn]$ for a wide range of mathematical models of transmissions among patients admitted to healthcare facilities. We showed that such calculations can provide actionable insights for understanding the factors leading to certain facilities reaching $[eqn]$ for high-priority pathogens, creating a setting where a single patient importation can lead to a sustained outbreak. We demonstrated that one set of high-risk facilities plausibly produced $[eqn]$ for CPE, a multi-drug-resistant organism of high concern, and analyzed models that provided nuanced insights into alternate interventions that could reduce $[eqn]$ in such facilities to below threshold.

A facility with $[eqn]$ is at risk not only of explosive outbreaks among its own patients, but also for acting as a transmission amplifier within a regional healthcare network [30]. Patients needing care in a high-risk (and high $[eqn]$ ) facility like an LTACH may be likely to subsequently require care in other hospitals or long-term care facilities, and thus they pose risk of continued spread of pathogens they had acquired in the LTACH. Prior simulation studies of CPE transmission among patients in a regional network of healthcare facilities, including short-stay hospitals, nursing homes and LTACHs, demonstrated that interventions including efforts focused on detecting carriage among LTACH patients could dramatically reduce CPE infections over the whole region [4].

Those prior findings of a dramatic decrease in simulated CPE infections due to LTACH-based intervention are consistent with a theoretical result we produced here: realistic rates of surveillance for CPE carriage and contact precautions for detected CPE carriers in an LTACH can plausibly move $[eqn]$ from above to below the critical threshold of 1. Crossing that threshold can explain the disproportionately high benefit achieved by a moderate increase in mitigation effort observed in simulation studies. A patient population health benefit that increases non-linearly with increasing intervention cost can also manifest in a highly efficient, cost-saving resource allocation, as predicted by a health economic analysis of a model of a CPE decolonization intervention in LTACHs [5]. However, our results in Fig 1 demonstrate one caveat to this insight: a facility with continually high rates of importation may not experience the same threshold-crossing benefit as a facility where colonized patients were more largely self-generated. For example, in a healthcare system with multiple $[eqn]$ facilities linked by patient exchange, a single healthcare facility reducing $[eqn]$ would not experience the full benefits of threshold effect if it still exists in an above-threshold system.

Our result for the effect of CPE surveillance on reducing facility $[eqn]$ was determined by the effectiveness of contact precautions to reduce transmission from detected carriers and rate of mid-stay surveillance only. The rate of admission surveillance does not appear in the $[eqn]$ formula, because we defined $[eqn]$ as the expected number of transmissions from a patient who acquires CPE after admission, between the time of acquisition and discharge or death. While admission surveillance can reduce the chances that otherwise undetected CPE importers will set off a chain of transmissions, nothing short of perfect admission detection and 100% effective contact precautions could entirely prevent a within-facility transmission from occurring. If any transmissions occur in a facility with $[eqn]$ by our definition, admission surveillance would not affect continuing circulation of CPE among in-patients who were each non-carriers at admission.

This insight exhibits why the admission reproduction number (defined $[eqn]$ in many studies), the expected number of transmissions from an admitted colonized patient, is not necessarily equal to $[eqn]$ , the expected number of transmissions from a colonized patient who acquired colonization after admission. If $[eqn]$ , then $[eqn]$ is not the correct threshold condition for characterizing outbreak risk. The implications of this finding are potentially important for facilities considering implementing a surveillance intervention: testing patients during the middle of their stay is the critical component that reduces the chances of a facility being a transmission amplifier. A caveat to this insight is that a broader definition of $[eqn]$ would cover a patient’s entire period of colonization, which might include multiple different stays at healthcare facilities. Admission surveillance, particularly for patients having recently stayed at the same or another high-risk facility, would help reduce $[eqn]$ under that extended model.

For the single-stay facility model, we showed that the biweekly rate of surveillance for CPE reported in a Chicago area LTACH intervention [21] may have been sufficient to reduce $[eqn]$ from above to below threshold, assuming 50% effectiveness of contact precautions at reducing transmission from detected patients, 95% adherence to surveillance testing, and 85% sensitivity of the test. Our $[eqn]$ formula provides a quick method for testing sensitivity of $[eqn]$ estimates to those assumptions as well as the rate of surveillance as an intervention decision variable.

In addition to contact precautions for identified colonized patients, the pathogen reduction intervention we tested represents the administration of treatment to detected carriers that would decolonize the patient or act to increase the rate of clearance by some mechanism, a class of treatments that has recently garnered attention [20]. Our prior modeling study demonstrated the potential for high cost-effectiveness that a drug with rapid decolonization effects could have when administered in a high-transmission facility to even a small portion of colonized patients [5]. Our results here bolster the theoretical underpinnings of that prior result and can help identify scenarios and settings where new drugs could help achieve transmission threshold effects that produces a disproportionately high health benefit. Our model combining surveillance and decolonization demonstrates a synergistic effect of the two strategies when the drug can achieve clearance rapidly.

The feasibility of implementing such a combined surveillance-and-decolonization intervention would depend on the method of pathogen reduction being used to achieve decolonization. For organisms like CPE and other pathogenic bacteria that colonize the gut, treatments with potential pathogen reduction agents are in early stages of conception, development, testing, or approval, such as selective decontamination of the digestive tract [31], bacteriophages [32], or probiotic treatments that could fortify the human microbiome’s natural ability to clear harmful colonization [33]. As the mechanisms, costs, and timings of such treatments become clearer, our model could be refined to produce cost-effectiveness analyses tailored to a particular approach.

As with any epidemiological modeling study attempting to project intervention effects in a real-world setting, our results depend on assumptions that may not be accurate. While we matched our assumptions to published data wherever possible, some key parameter value assumptions did not have strong data-based justification. For example, the relative transmissibility of detected CPE carriers compared to undetected carriers in hospitals is unknown, as the transmitter to an acquisition is rarely identified definitively during CPE outbreaks. In addition, given that the decolonization drug explored in our model is hypothetical, modeling its effects relative to the also-uncertain natural clearance rate is largely based on guesswork. Nevertheless, our equation-based results allow for quick exploration of alternate assumptions, and such sensitivity analyses in this and prior work [5,23] suggest that our findings for these intervention effects are largely robust.

Another major finding of this study is a characteristic property of the length of stay distribution of facility patients that largely governs its relationship to $[eqn]$ . Namely, $[eqn]$ was directly proportional to the ratio of the second raw moment of the length of stay distribution to the mean of the distribution. The second raw moment is the average of the squares of each individual length of stay and is also equal to the sum of the variance and the square of the mean. Thus, the characteristic ratio can also be expressed as the sum of the variance and the variance-to-mean ratio of the distribution. This is an exact result for Model 1, where colonized patients remain colonized from the time of acquisition to discharge or death with no change in their transmissibility. For more complicated models, the proportionality is not exact, as the $[eqn]$ formulas depend on higher moments of the length of stay distribution via the moment-generating function. But we showed in Fig 3 that the characteristic ratio was nevertheless an excellent predictor of changes to $[eqn]$ due to changes in a complicated length of stay distribution model.

Intuition as to why the square of lengths of stay is a characteristic property governing facility $[eqn]$ lies in understanding the length of stay of a typical patient who might acquire colonization while they are an in-patient at a healthcare facility. If a random patient residing in the facility on a given day becomes colonized, their time remaining in the facility before discharge or death is the characteristic time duration governing their subsequent potential to transmit to others. Calculating the average length of stay of a cross-section of patients in the facility at a given time is different than calculating the average length of stay of a set of admitted patients, because longer-stay patients are overrepresented in the cross-section due to their taking up facility bed spaces for longer. Specifically, in a sample of patients in the facility at a given time, patients with a longer length of stay are more likely to be sampled, in direct proportion to their length of stay. For example, a patient with a length of stay twice as long as another was twice as likely to be in the facility on the day of sampling. Therefore, the mean length of stay of the patients sampled is a weighted mean of the actual lengths of stays of all admitted patients, where the weights are the length of stay. The weighted mean of lengths of stay weighted by length of stay reproduces the characteristic quantity: the mean of the square of length of stay divided by the mean.

This finding is an example of well-known result in theoretical ecology of infectious disease [34]. In applications, the sum of the mean and the variance-to-mean ratio of a transmission-relevant quantity being proportional to $[eqn]$ has been shown in a variety contexts, including the number of sexual partners per person infected with sexually transmitted disease [35], the biting rates of mosquitos transmitting disease to humans [36], and household sizes for transmission through a community via housemates [37]. Our finding that this result also applies to facility length of stay for healthcare-associated infections is, to our knowledge, a novel contribution to the field of healthcare transmission modeling.

The practical implications of this result are that both the mean and the variance of the length of stay distribution are important for characterizing facility outbreak risk, and both should be more regularly calculated and considered by facilities undergoing transmission studies. Reported facility length of stay statistics often include only the mean and/or quantiles such as the median or inter-quartile range, which fail to uniquely identify the variance. Individual facilities and facility types vary widely in their length of stay patterns. Aside from LTACHs, other types of long-term care facilities such as nursing homes or skilled nursing facilities may produce different quite different average and variation in length of stay compared to our example, and those patterns may rapidly and unevenly fluctuate due to shifting policies and other factors [38]. Short-stay hospitals generally have a much lower mean length of stay that reduces the risk of sustained transmission after a single importation, but some hospitals could have an unusually high length of stay variance due to a subset of long-stay patients, which would put them at greater risk than the mean length of stay alone would indicate. Our findings produce a unifying framework for analyzing and comparing facilities with widely different and complicated length of stay patterns: the sum of the mean and the variance-to-mean ratio in length of stay correlates with $[eqn]$ and thus with outbreak risk.

Furthermore, we showed that the effect of efforts to reduce lengths of stay as a transmission risk-reduction intervention are not always well characterized by the intervention’s effect on the mean. If the mean length of stay is reduced in a way that simultaneously increases the variance-to-mean ratio, it is possible that facility outbreak risk could increase, as measured by $[eqn]$ . For LTACHs and similar facilities that specialize in care for high acuity patients, it might be possible to characterize the condition and care needs of incoming patients to estimate not only their expected length of stay, but also the uncertainty (i.e., variance) in days of care they will require. If a facility can control the composition of admitted patients and/or target length of stay reduction efforts on certain groups, efforts to reduce the fraction of the high-variance group of patients or to increase their discharge rates would be favored for reducing facility outbreak risk.

Our general model has capabilities for investigating additional nuances of length of stay effects, including use of state-dependent discharge or death rates via the $[eqn]$ removal terms. To reduce unnecessary complexity in our examples, we set $[eqn]$ and focused on using the state-independent removal hazard h(t) to model and calibrate the length of stay distribution. Our assumption ignores potential correlation between patient states and discharge and death hazard. For example, patients experiencing symptomatic infection could have a reduced or eliminated discharge hazard and a higher death rate, and patients with more serious underlying conditions may have both a lower discharge hazard and a higher per-exposure susceptibility than other patients. Calibrating a model with these details to length of stay statistics would require more care but is possible under our general framework, as shown in our prior work [5,23].

We have released an R package, facilityepimath [29], that includes a function to calculate $[eqn]$ numerically, using the general equation we derived, for a wide variety of models. The package contains example code reproducing the results for the models we investigated in this manuscript, as well as example demonstrations of how other users can use the function to calculate $[eqn]$ for any model that fits our broad general framework, including models with nuanced assumptions governing facility length of stay. The package also provides other functions, e.g., to calculate the equilibrium of a model for a given, constant rate of importation of infectious individuals to the facility, which is useful for calibrating models to observed levels of admission and cross-sectional prevalence, as done in this work, and for assessing the impacts of transmission-reducing interventions under ongoing importation [5,23].

While our formula applies to a wide range of possible models, there are limitations to what can be assumed. For contact assumptions affecting transmission, the model form requires that all susceptible compartments have the same relative exposure to individuals in different colonized compartments. Therefore, our results do not apply to segregated contact models in which, for example, the compartments represent patients in different wards of a healthcare facility with different within-ward vs. cross-ward contact rates. Our model also does not include the effects of transmission during an episode of infectiousness extending beyond a single facility stay. Extensions to multiple facility stays via readmission to the same facility and transmission in multiple different facilities or in the community outside a healthcare facility are potential subjects of future work. Our model parameter estimates from CPE data were possible because the facilities observed steady state levels of importation and cross-sectional prevalence, which made our equilibrium assumption reasonable. Fitting parameters to facilities observing more transient patient state dynamics would require additional model components and solution techniques, as would studying the transient dynamics toward a new equilibrium after introducing an intervention. With a deterministic framework, our model also does not quantify stochastic variation in outbreaks stemming from a single introduction. In a stochastic model simulation, a sustained chain of transmission after a single introduction would have a much higher probability of occurring when $[eqn]$ , but would not be guaranteed. With high variability in factors affecting transmission, a facility with $[eqn]$ might experience many introductions leading to only small outbreaks before a large outbreak eventually takes hold [39].

In conclusion, we have derived a novel, general mathematical result for theoretical quantities that have important implications for study and control of high-priority infectious disease health threats. Our findings for healthcare facility transmission models could be applicable to other at-risk settings with similarly rapid population turnover. Locations with above-threshold transmission conditions can pose a significant risk to population health, while at the same time presenting opportunities for targeted interventions that are highly efficient. Seeking a nuanced understanding of those conditions and their relationship to intervention mechanisms should be a component of strategic planning for public health efforts.

Supporting information

S1 TextSupplementary methods and results.(PDF)