Transmission thresholds for the spread of infections in healthcare facilities

Transmission thresholds for the spread of infections in healthcare facilities | medRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-P4HH5NV'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Transmission thresholds for the spread of infections in healthcare facilities View ORCID Profile Damon J.A. Toth , View ORCID Profile Karim Khader , View ORCID Profile Christopher Mitchell , Matthew H. Samore doi: https://doi.org/10.1101/2025.02.21.25322698 Damon J.A. Toth 1 Division of Epidemiology, Department of Internal Medicine, University of Utah School of Medicine , Salt Lake City, Utah, United States of America 2 Department of Veterans Affairs Salt Lake City Healthcare System , Salt Lake City, Utah, United States of America 3 Department of Mathematics, University of Utah , Salt Lake City, Utah, United States of America Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Damon J.A. Toth For correspondence: Damon.Toth{at}hsc.utah.edu Karim Khader 1 Division of Epidemiology, Department of Internal Medicine, University of Utah School of Medicine , Salt Lake City, Utah, United States of America 2 Department of Veterans Affairs Salt Lake City Healthcare System , Salt Lake City, Utah, United States of America Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Karim Khader Christopher Mitchell 4 Department of Mathematics, Tarleton State University , Stephenville, Texas, United States of America Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Christopher Mitchell Matthew H. Samore 1 Division of Epidemiology, Department of Internal Medicine, University of Utah School of Medicine , Salt Lake City, Utah, United States of America 2 Department of Veterans Affairs Salt Lake City Healthcare System , Salt Lake City, Utah, United States of America Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Data/Code Preview PDF Abstract Some infections may be sustained in the human population by persistent transmission among patients clustered 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 infections among high-risk patients. Using a mathematical model, we define the facility basic reproduction number R 0 , where a single facility can sustain an outbreak without ongoing importation under the threshold condition R 0 > 1. We define R 0 for a general model with heterogeneous patient states of 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 R 0 using published data for carbapenemase–producing Enterobacteriaceae (CPE) in longterm acute-care hospitals (LTACHs) and the effects of interventions on R 0 , including surveillance, pathogen reduction treatments, and length-of-stay reduction. In a simple model, R 0 is directly proportional to the sum of the mean and variance-to-mean ratio of the length-of-stay distribution. In intervention models, R 0 depends on the moment-generating function of the length-of-stay distribution. From the CPE data, we estimated R 0 = 1.24 (95% CI: 1.04, 1.45) prior to intervention. Weekly surveillance with 50% transmission reduction of detected patients alone could have reduced R 0 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 R 0 if the variance-to-mean ratio is not also reduced. We conclude that R 0 >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. 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 interfacility 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, R 0 , the expected total number of transmissions from a typical infected individual in an unexposed population ( 6 ). If R 0 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 R 0 can determine the theoretical level of transmissionreducing 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 R 0 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, R A ( 12 - 17 ), or ward reproduction ratio, R w ( 18 ), noting that the quantity is less than R 0 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 R 0 >1 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 R 0 below threshold. Our aims for this work were first to define a facility reproduction number for a broad class of mathematical models 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. 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 i , and let x i ( t ) be the probability that a patient is alive, not discharged, and in state i at time t after admission. Defining the column vector of probabilities x ( t ) containing each x i ( t ), we define the system of differential equations Here, W contains state-to-state transition rates and may also contain state-specific facility removal rates by death or live discharge. The hazard function h ( t ) is a facility removal rate at time t of the stay, where the same function applies to every state – this function can be used to calibrate a model to an overall observed patient length of stay distribution. The vector x a 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 the dynamics of a particular organism that can colonize or infect patients in the facility, we assume the patient states x are comprised of n susceptible states ( S 1 ,…, S n ) and m colonized states ( C 1 ,…, C m ). The colonized states represent patients carrying the modeled organism. The m 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 n 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: In the initial state vector x a describing the admission state distribution, p a is the overall probability of being in a colonized state at admission, ( θ 1 ,…, θ n ) is the admission state distribution among the susceptible state admissions, with ∑ θ i =1, and ( k 1 , … k m ) is the admission state distribution among the colonized state admissions, with ∑ k i =1. In the matrix W , the terms ( ω 1 , …, ω n + m ) 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 n × n portion of W contains elements describing the transitions out of and between the susceptible states. The s ij terms are the rates of transitioning from one susceptible state to another, with s jj = ∑ i ≠ j s ij if n >1 and s 11 =0 if n = 1. The s j α terms are the rates of transitioning from each susceptible state to a colonized state (rates of acquisition). The s j coefficients describe relative susceptibility to acquisition of patients in each susceptible state, and α 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: where the ( β 1 , …, β m ) coefficients define the relative transmissibility of patients in the different colonized states. The lower-left m × n portion of W contains rates of transitioning into each colonized state from each susceptible state, with s j = ∑ a ij . The lower-right m × m portion of W contains elements describing the transitions out of and between the colonized states. For the c jj elements in the diagonal we have c jj = ∑ i ≠ j c ij + ∑ r ij , where the r ij terms comprising the upper right n × m portion of W 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 R 0 for given values of the importation state distribution, state transition rates, and relative transmissibility parameters defined above, using the procedure outlined in Diekmann et al. ( 6 ). We also derive R 0 formulae for a series of specific model examples. 2.2 Example model descriptions 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 x consists of one susceptible state ( S 1 ) that we rename S and one colonized state ( C 1 ) that we rename C . The admission state probability vector in the above notation contains θ 1 = 1 and k 1 = 1: We assume that susceptible inpatients acquire colonization at rate α , and there are no other state transitions or removals other than removal modeled by the hazard function h ( t ). Thus, in the general model notation, we have s 1 =1, a 11 =1, r 11 =0, ω 1 = ω 2 =0. The state transition matrix W is then For transmission, there is only one colonized state from which transmissions can occur, so there is a single transmission parameter ( β 1 ) that we rename β , and the acquisition rate α = βC . The time-of-stay-dependent removal hazard h ( t ) is left general. The equivalent system in non-matrix form is as follows: 2.2.2. Model 2: Clearance of colonization To extend Model 1, we add a non-zero clearance rate γ of moving from colonized back to susceptible. In the general model notation, we have s 1 = 1, a 11 = 1, r 11 = γ, ω 1 = ω 2 =0. The state transmission matrix W is then In non-matrix equation form, we have: 2.2.3 Model 3: Clinical detection with contact precautions For this model, we extend Model 2 to include two colonized states ( C 1 , C 2 ) that we rename C and C cd . 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 θ 1 = 1, k 1 = 1, and k 2 = 0, i.e., all colonized patients are undetected at admission: For state transitions among inpatients, we assume that susceptible patients acquiring colonization (at rate α ) move into the undetected colonized state. Undetected colonized patients move to the susceptible state at clearance rate γ and move to the clinically detected state at clinical detection rate δ c , 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: The state transmission matrix W is then 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 ε of contact precautions: β 1 = β, β 2 = β (1 − ε ). Then the acquisition rate α = β ( C +(1 − ε ) C cd ). In non-matrix equation form, Model 3 is written as follows: 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 x consists of two susceptible states ( S 1 , S 2 ) that we rename S and S s d and three colonized states ( C 1 , C 2 , C 3 ) that we rename C,C s d , and C cd . 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 Π a , 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 θ 1 = 1, θ 2 = 0, k 1 = 1 − Π a , k 2 = Π a , k 3 = 0: For state transitions among inpatients, we assume that undetected colonized patients become surveillance detected at rate δ s , 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 δ c as undetected colonized patients. The new susceptible state S s d 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, (1 − ε ), that they reduce transmissibility. In the general model notation we have the following nonzero elements: s 1 = a 11 = 1, s 2 = a 22 = 1 − ε, r 11 = γ, c 21 = δ s , c 31 = δ c , r 32 = γ d , c 32 = δ c , ω 1 = ω 2 = ω 3 = ω 4 = ω 5 = 0. The state transmission matrix W is then For relative transmission rates, we assume that surveillance detected, colonized patients transmit at a rate commensurate with clinically detected patients: β 1 = β, β 2 = β 3 = β (1 − ε ). Then the acquisition rate α = β ( C + (1 − ε )( C s d + C cd )). In non-matrix equation form, Model 4 is written as follows: 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, r x , the rate and shape parameters of a gamma distribution, r g and k , and the portion of patients following the exponential distribution, p x . 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 (Supplementary Material). View this table: View inline View popup Download powerpoint Table 1. Data from pre-intervention LTACHs ( 21 ) For the CPE parameters in Model 3, we assumed the CPE clearance rate γ =1 / (387 days) ( 24 ), contact precaution effectiveness ε =0.5 ( 25 ), and CPE test sensitivity = 0.85 ( 26 , 27 ). Because the LTACH importation rate and crosssectional prevalence of CPE were observed and stable over time, we assumed the Model 3 system was at equilibrium, which means that the acquisition rate, α , was constant, making Model 3 a linear system. We solved for the value of α and of the progression rate to clinical detection, δ c , 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, β , that produced the correct acquisition rate under the calibrated model ( Table 2 and supplementary material). Finally, we used the assumed values of γ and ε and the data-calibrated values of β and δ c to calculate R 0 using the formula we derived for Model 3 (see Results). View this table: View inline View popup Download powerpoint Table 2. Model 3 parameters and assumptions 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 R 0 . 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, δ s , 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): δ s ≈ 0.058 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 R 0 < 1. 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 f > 1, such that γ d = fγ . For values of f = 2, 5, 10, 50, 100, and ∞ (immediate decolonization at detection), we calculated R 0 for weekly and biweekly surveillance, and the surveillance rate required to achieve R 0 < 1. 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 R 0 , using Model 3 calibrated to LTACH CPE data ( Table 1 and Table 2 ) as the starting point. As the overall mean length of stay decreases as either of the rate parameters r x or r g increases, we tested the R 0 -altering effects of increasing either r x , r g , 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 p x 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 R 0 (again using Model 3 calibrated to LTACH CPE data) vs. different length-of-stay distribution statistics. We plotted R 0 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 R 0 . 3. Results 3.1 General facility basic reproduction number (R 0 ) 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: The elements in the matrices S, A , and C were defined within the matrix W in the Methods section. The matrix S describes transition rates between the susceptible states in the absence of colonized patients, and state-specific removal rates from the facility along the diagonal ( ω i terms). The matrix C contains transition rates from and between the colonized states, including state specific removal rates from the facility along the diagonal ( ω i terms). The matrix A contains the relative transition rates a ij from susceptible state j to colonized state i upon acquisition. The vector β contains the transmission rates from each colonized state. The full transition rate from susceptible state susceptible state j to colonized state i is a ij ( β 1 C 1 + … + β m C m ). The vector θ contains the fractions of admitted, susceptible patients who are in each of the susceptible states at admission. The function F ( t ) is the cumulative distribution function of the facility length of stay distribution of patients who are removed via the hazard function h ( t ). With these elements in hand, we derived the following result for the facility basic reproduction number R 0 (Supporting Information): This formula is a dot product of the vector β , 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 R 0 can be expressed symbolically, which generally occurs when the eigenvalues of matrices S and C can be expressed as formulas in terms of the model parameters. Otherwise, R 0 can be calculated numerically as follows. In the numerator of the above expression for R 0 , the vector is equivalent to final m elements of the vector , in which and 0 i×j contains i rows and j columns of zeros. We numerically calculate the eigenvalues λ i of and associated eigenvectors to determine the linear combinations of functions that comprise the elements of .Then, integrating the product of this vector and (1 − F ( t )) produces a vector with linear combinations of K ( λ i ), where K ( x ) = (1 − M ( x )) ⁄ x and M ( x ) is the moment generating function of the distribution for which F is the cumulative distribution function. For a model with the hazard function h ( t )=0, the patients ’ rate of removal from the facility (by discharge or death) is entirely state-dependent and quantified by the ωi values within the matrices. Then the function F ( t )=0and the basic reproduction number formula (supplementary material) reduces to: 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 i from which the removal rate ω i >0. We created an R package, facilityepimath [ https://epiforesite.github.io/facilityepimath/ ] containing a function that numerically calculates R 0 according to the above formula, for a given set of numerical inputs for S, C, A, β , θ , and M (if h ( t ) ≠ 0). 3.2 Example model R 0 formulas 3.2.1 Model 1: Simple susceptible–colonized model For this model we have: and the hazard function h ( t ) is left general. Substituting these values in the general formula gives In the supplementary material we show that , and , where μ and σ 2 are the mean and variance, respectively, of the distribution for which F ( t ) is the cumulative distribution function. Therefore, our result is: 3.2.2. Model 2: Clearance of colonization For this model we have: After substitution and some integration (see more detailed steps in the supplementary material), the result is: In the supplementary material we demonstrate that, for a real number , where M is the moment-generating function of the distribution for which F is the cumulative distribution function. For ease of notation, we define the function K ( λ ) = ( M ( λ ) − 1) ⁄ λ . Our result for this model can then be expressed Some commonly used parametric distributions for modeling length of stay provide convenient expressions for the moment-generating function M within the formula for K . For example, for length of stay exponentially distributed with rate r (constant removal hazard), we have μ = 1 ⁄ r and K ( λ ) = 1 ⁄( r − λ ), which leads to: With time to removal gamma distributed with rate parameter r and shape parameter k , we have μ = k ⁄ r and K ( λ ) = ((1 − λ / r ) − k − 1) ⁄ λ , which leads to 3.2.3 Model 3: Clinical detection with contact precautions For this model we have: The formula for R 0 contains e C t , which, because C is a matrix, is the matrix exponential: Substituting this and the other expressions in the R 0 formula, we arrive at the following results after integrating and performing matrix multiplications: 3.2.4 Model 4: Active surveillance and decolonization For this model we have: The matrix exponentials e S t and e C t can be expressed: The values c 11 and c 22 are the first and second diagonal elements of the matrix C . Substituting these and the other expressions in the R 0 formula, we arrive at the following result after integrating and performing matrix multiplications: In the case that δ s − ( γ d − γ ) = 0, the above expression is indeterminate because of division by zero; this corresponds to the matrix C having a repeated eigenvalue ( c 11 = c 22 = δ s − δ c − γ ) that is degenerate, and the matrix exponential e C t includes eigenfunctions and . In this case, the resulting expression for R 0 is as follows: 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 γd = γ and R 0 simplifies to: Second, when decolonization is applied to surveillance-detected patients, and the clearance effect is assumed to be instantaneous ( γ d → ∞), we can use the result from Model 2 with γ replaced by the sum of γ and δ s . That is, 3.3 Calibration of Model 3 to pre-intervention CPE in LTACH For the length of stay distribution, we assume a mixed exponential-gamma distribution, with a portion p x of patients following an exponential distribution with rate r x and the rest following a gamma distribution with rate r g and shape parameter k . 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 (Supplementary Methods), we calculated the unique solution ( Table 3 ). View this table: View inline View popup Download powerpoint Table 3. Model calibration and R 0 results: CPE LTACH pre-intervention For the CPE parameters, under our assumed values for three independently based parameters ( ε =0.5, γ =1 / 387 per day, σ =0.85), we estimated the pre-intervention values of β, δ c , and R 0 from CPE data in LTACHs, using Model 3 ( Table 3 ). Our estimate for the baseline transmission rate, β , was 0.051 per day (95% range 0.043 to 0.059). The estimate of rate of progression to clinical detection, δ c , was 0.0084 per day (0.0073 to 0.0097). Using those results and our formula for Model 3 above, we estimated R 0 = 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 R 0 from Model 4, with γ d = γ (no decolonization effect) and surveillance detection rate δ s =0.058 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 R 0 = 0.94 (95% range 0.79 to 1.08). Increasing the frequency of surveillance from biweekly to weekly decreases the estimate to R 0 = 0.85 (0.72 to 0.98). The minimum frequency of surveillance testing required to make R 0 < 1 was once per 3.2 weeks (1.1 to 23 weeks) ( Table 4 , top row). View this table: View inline View popup Download powerpoint Table 4: R 0 and threshold estimates for combined surveillance and decolonization intervention 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 f > 1, such that γ d = fγ . For different values of f , we calculated R 0 for weekly and biweekly surveillance, and the surveillance rate required to achieve R 0 < 1 ( Table 4 ). For example, with f =5 (five-fold increase in clearance rate due to the pathogen reduction treatment), the R 0 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 ( f → ∞, or immediate clearance upon detection) achieves R 0 = 0.56 (0.47 to 0.64) with biweekly surveillance. With immediate clearance, a surveillance interval of once per 10 weeks (6 to 56 weeks) could potentially achieve R 0 < 1. The code producing the results in Table 3 and Table 4 is available as an installed file in our facilityepimath R package [ https://github.com/EpiForeSITE/facilityepimath/tree/main/inst ]. 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 μ =33.8 days and standard deviation σ = 28.1 days. With the mixture model we used to govern removal, a fraction p x =0.580 of admitted patients follow an exponentially distributed length of stay (constant discharge-or-death hazard), with mean and standard deviation μ x = σ x = 35.1 days. The remaining patients follow gamma-distributed length of stay with shape parameter k >1 (meaning the removal hazard increases with time of stay), which gives a length of stay distribution with mean μ g = 32.0 days and standard deviation σ g = 13.4 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 R 0 for Model 3 calibrated to pre-intervention LTACH CPE data ( Figure 1 ). 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 r x and r g by the same factor), R 0 decreases and crosses the R 0 < 1 threshold as the mean length of stay drops below about 26 days ( Figure 1A ). Second, when increasing the discharge rate r x of the high-variance patients only, R 0 decreases below 1 as the mean length of stay drops to about 28 days ( Figure 1B ). Third, when increasing the discharge rate r g of the low-variance patients, R 0 does not decrease below 1.2 and in fact begins increasing as the mean length of stay drops below about 29 days and lower ( Figure 1C ). Fourth, when decreasing the fraction p x of patients following the high-variance exponential distribution, R 0 drops rapidly and falls below 1 as the mean length of stay drops below about 32.6 days. When p x = 0 (all patients follow the low-variance gamma distribution with mean 32.0 days), R 0 = 0.88. Download figure Open in new tab Figure 1. Effect of length-of-stay interventions on facility basic reproduction number Calculations of the facility reproduction number, R 0 , for Model 3 with the following fixed values calibrated to pre-intervention CRE data in longterm 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–gammadistribution, 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 r x and the gamma rate parameter r g were simultaneously increased by the same factor from their baseline values. Upper-right: the exponential rate parameter r x was increased from its baseline value. Lower-left: the gamma rate parameter r g was increased from its baseline value. Lower-right: the fraction of patients following the exponential distribution p x was decreased from its baseline value. As the results in Figure 1 demonstrate that the mean length of stay did not necessarily correlate with R 0 , we tested the relationship between R 0 and other statistics for the length of stay distribution ( Figure 2 ). The dots in the scatterplot represent the results from randomly chosen combinations of the four length-of-stay parameters p x , r x , r g , and k , 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 R 0 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 R 0 formula for Model 1. Download figure Open in new tab Figure 2. Scatterplot of facility basic reproduction number vs. length of stay (LOS) distribution statistics Calculations of the facility reproduction number, R 0 , for Model 3 with the following fixed values calibrated to pre-intervention CRE data in longterm 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 R 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 R 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 R 0 (Pearson ’ s correlation coefficient PCC = 0.857) compared to the correlation between the standard deviation and R 0 (PCC = 0.986). The VMR produced PCC = 0.860, and the sum of the mean and VMR produced the strongest correlation with R 0 , with PCC = 0.997 and the relationship appearing nearly linear ( Figure 2 ). The code producing Figures 1 and 2 is available as installed files in our facilityepimath R package [ https://github.com/EpiForeSITE/facilityepimath/tree/main/inst ]. 4. Discussion This study provides a flexible, quantitative framework for calculating the basic reproduction number R 0 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 R 0 > 1 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 R 0 >11 for CPE, a multi-drug-resistant organism of high concern, and analyzed models that provided nuanced insights into alternate interventions that could reduce R 0 in such facilities to below threshold. A facility with R 0 > 1 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 ( 29 ). Patients needing care in a high-risk (and high R 0 ) 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 R 0 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 ). Our result for the effect of CPE surveillance on reducing facility R 0 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 R 0 formula, because we defined R 0 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 R 0 > 1 by our definition, admission surveillance would not affect continuing circulation of CPE among in-patients who were each non-carriers at admission. Thus, our finding provides a potentially important insight 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 R 0 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 R 0 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 R 0 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 R 0 formula provides a quick method for testing sensitivity of R 0 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 ). 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 ( 30 ), bacteriophages ( 31 ), or probiotic treatments that could fortify the human microbiome ’ s natural ability to clear harmful colonization ( 32 ). 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. 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 R 0 . Namely, R 0 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 R 0 formulas depend on higher moments of the length of stay distribution via the moment-generating function. But we showed in Figure 3 that the characteristic ratio was nevertheless an excellent predictor of changes to R 0 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 R 0 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 suggests that measures identifying each of the first two moments (e.g. mean and variance or mean and standard deviation) 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 second moment. 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 R 0 . 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. We have provided an R package, facilityepimath, that includes a function to calculate R 0 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 R 0 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. 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. Data Availability All data produced are available online at https://github.com/EpiForeSITE/facilityepimath/ Supplementary Methods and Results Length of stay model calibration to LTACH data We assume that patient length of stay is a random variable with cumulative distribution function F ( t ). The source LTACH study ( 21 ) provides the mean and 25 th , 50 th , and 75 th percentiles of the overall length of stay distribution ( μ, l 25 , l 50 , l 75 ). For fitting Model 3 to data, we have four equations constraining the length of stay distribution: We assume a mixed exponential-gamma distribution for the length of stay distribution, with a portion p x of patients following an exponential distribution with rate r x and the rest following a gamma distribution with rate r g and shape parameter k . For this distribution: We apply these two functions to the four equations above and numerically solve for the four unknown parameters p x , μ x , μ g , and k . From the LTACH study ( 21 ), we have μ = 33.8 days, l 25 = 16 days, l 50 = 28 days, and l 75 = 43 days, which leads to: Equilibrium calibration to LTACH CPE data As described in the main text, we calibrated the patient acquisition rate α at equilibrium and the per-capita clinical detection rate δ c of a colonized patient to data describing the cross-sectional prevalence of CPE carriage and the overall CPE clinical detection incidence rate observed at the LTACHs ( 21 ). We assumed that Model 3 (main text) was at equilibrium with the observed constant CPE admission positivity rate. Therefore, we required a formula for the equilibrium prevalence of the three patients states of Model 3 at the given initial conditions (admission states). Here we derive a formula for the equilibrium of a general model, then apply the formula to Model 3. Let x ( t ) be the vector of probabilities that a patient is alive, not discharged, and in each state x i at time t after admission. We have the system Here, W contains state-to-state transition rates and state-specific facility removal rates during the stay. The hazard function h ( t ) is the removal rate at time t of the stay, where the same function applies to every state. The vector x a is the distribution of states at admission. The solution to this equation is expressed If we express the hazard function as that of a probability distribution: h ( t ) = f ( t ) ⁄ (1 − F ( t )), where f and F are the probability density function and cumulative distribution function of the same distribution, then we can express x ( t ) as If there are no discharge or death rates within the W matrix, then F is the cumulative distribution function of the facility length of stay distribution. If we assume that patients are continually admitted to the facility at a constant rate of qpatients per time, ad infinitum , then the expected number of patients in the facility between stay-time t 1 and t 2 in each state is . So, the expected total number of inpatients in each state is given by the vector And the expected total number of patients in the facility is N ∗ = 1 T n ∗ . Then, the equilibrium cross-sectional state distribution within the facility is x ∗ = n ∗ ⁄ N ∗ , or To evaluate the integral in this expression, we must consider the integral for the eigenfunctions φ that comprise the elements of the matrix exponential e W t . If the state dependent removal rates ω i are all nonzero, the eigenvalues λ i of W are negative and real. The eigenfunctions associated with non-degenerate eigenvalues then take the form . We use integration by parts: Here, f and M are the probability density function and moment-generating function of the distribution associated with the removal hazard, respectively. For ease of notation, we define the function If there are any repeated, negative eigenvalues that are degenerate, there would be eigenfunctions of the form for integer r ≥1. We note that if then , where K ( r ) is the r th derivative of K , defined above. Finally, for models where any state-dependent removal rates ω i are set to zero, there may be an eigenvalue of W equal to zero, in which case there will be at least one eigenfunction φ i ( t ) that is constant. Then there will be a component of the expression equal to the mean μ of the distribution for which F ( t ) is the cumulative distribution function, as . In summary, the vector is constructed by calculating e W t x a and replacing the terms in each vector element with K ( r ) ( λ i ) and the constant terms with μ . For Model 3, we have the following components: To calculate the matrix exponential e W t we require the eigenvalues of W . The first two eigenvalues are the eigenvalues of the upper-left 2-by-2 submatrix, and the third eigenvalue is 0. The first two eigenvalues λ 1 and λ 2 are: Then, the matrix exponential is: For the function K we require the moment generating function M of the mixed exponential–gamma distribution that was calibrated to the LTACH length of stay distribution as described above: Then, Using the above formulas, we solved for the values of α and δ c that produce an equilibrium result that meets the following two criteria: The cross-sectional fraction of patients who are in one of the colonized states at equilibrium is the cross-sectional positivity f divided by the assumed surveillance test sensitivity σ : The clinical detection rate in the facility at equilibrium is d : The values of f, σ , and d , as well as the other parameters aside from α and δ c within the formulas for C ∗ and above are determined by the data and/or by assumption as described in Tables 1 and 2 in the main text. Therefore, we have two equations to solve simultaneously for the unique solution of α and δ c . Facility basic reproduction number formula derivation For the general system defined in the main text, we will show that the facility basic reproduction number Where To derive our equation we show that it is consistent with a “ basic reproduction ratio ” formula in Diekmann et al ( 6 ). Their definition depends on a measure of an individual ’ s “ h-state ”, η , just prior to acquiring colonization. First, they define the quantity S ( η ): the distribution of h-states among susceptible individuals at equilibrium in the absence of any colonized patients. Then they define the quantity A( τ,ξ,η ), which is the expected infectivity of an individual who was in h-state η just prior to acquisition time τ ago, toward a susceptible individual currently in h-state ξ . In our model, the h-state is characterized by both the individual ’ s susceptibility compartment, s ∈ (1,…, n ), and time since admission, t ∈ (0,∞). So we have S ( η ) = S ( s η , t η ) and A( τ,ξ,η ) = A( τ,s ξ , t ξ , s η , t η ). The time-since-admission component is only relevant for the infectivity of the currently colonized individual represented by η , whose probability of removal (death or discharge) from the facility (thus no longer infective toward other facility patients) between t and t + τ might depend on t . However, the infectivity toward the susceptible individual represented by ξ depends only on that individual ’ s current compartment s , via the susceptibility factor ψ s (the probability of being removed while susceptible will be incorporated in the S function). Therefore, A does not depend on t ξ , and for simplicity of notation we let t = t η and seek expressions for the functions S ( s η , t ) and A( τ,s ξ , s η , t ). The infectivity A will be the susceptibility ψ sξ times the sum of infectivity values β c of each possible colonized compartment c , weighted by the probabilities p ( c,s η , t,τ ) of being in colonized compartment c time τ after acquisition at stay-time t from compartment s η and still in the facility. So we have: This expression for A satisfies a separability condition in Diekmann et al., i.e. it can be written as a product of two functions Where Therefore, the following formula from Diekmann et al. for the basic reproduction ratio applies: For simplicity we let s = s η and specify the R 0 formula as First, from our equilibrium results above we have that S ( s,t ) is the s th element of the column vector: Next, we show that p ( c,τ,s,t ) is the element in the c th row, s th column of the matrix where A 0 is the matrix A above with each column scaled by its sum. The s th column of A 0 represents the distribution of colonized compartments initially entered when acquiring from compartment s , and e C τ projects the initial compartment forward τ time units according to the colonized state transition dynamics embedded in C . The expression (1 − F ( t + τ )) ⁄ (1 − F ( t )) is the probability that removal does not occur before time-of-stay t + τ given that the stay was longer than t (the acquisition time). So is the element in the c th row, s th column of the matrix Noting that the sum of column s of the matrix A is ψ s , we have that is the element in the c th row, s th column of the matrix Substituting into the R 0 formula above, we have We use the substitution z = τ + t and rearrange terms to get Integrate the numerator by parts using the following substitutions: ,so In the numerator of this expression for R 0 , the integral with respect to t can be calculated similarly to integral in the denominator, discussed above. I.e., if all state-dependent removal rates are nonzero, the elements of the vector are linear combinations of functions with negative eigenvalues can be calculated similarly to integral i of the S and C matrices and integer r ≥0. Then, integrating the product of this vector and (1 − F ( t )) produces a vector with linear combinations of K ( r ) ( λ i ), where K ( x ) = (1 − M ( x )) ⁄ x and M ( x ) is the moment generating function of the distribution for which F is the cumulative distribution function. If there are state-dependent removal rates equal to zero, then the vector may contain a term linear in t . In which case we evaluate where μ and σ are the mean and standard deviation of the distribution for which F is the cumulative distribution function. Examples No time-dependent removal hazard For the model with removal hazard h ( t ) = 0, and thus F ( t ) = 0, it is convenient to return to the following expression for R 0 within the derivation above: Plugging in F ( t ) = 0, we have The following are true for invertible matrices C and S : Hence Model 1: Simple susceptible–colonized model Model 2: Clearance of colonization With length of stay exponentially distributed with rate r : With length of stay gamma distributed with rate parameter r and shape parameter k Model 3: Clinical detection with contact precautions When applying this model to LTACH data, we assume a mixed exponential-gamma distribution for the length of stay, with a portion p x of patients following an exponential distribution with mean μ x and the rest following a gamma distribution with mean μ g and shape parameter k . For this distribution: Model 4: Active surveillance and decolonization When δ s − ( γ d − γ ) ≠ 0: Numerical calculation of R 0 As demonstrated in the examples above, when the matrix exponentials e S t and e C t can be expressed symbolically for a particular model, the R 0 formula we derived can be used directly to find a formula for R 0 in terms of the model parameters. We can also calculate R 0 numerically using the following procedure. First, we create the following matrix M : With S and C being n × n and m × m matrices, respectively, the 0 in the upper-right corner of M represents an n × m block of zeros. The 0 in the column vector θ M represents a column vector of m zeros. The system can be decoupled into two linked systems: Solved sequentially we get: Returning to our R 0 formula: We can substitute in x s and x c to arrive at an alternate, equivalent expression: We can numerically calculate the numerator and denominator of that expression using the following procedure: Numerically solve for the eigenvalues and eigenvectors of the matrix M . Place the eigenvalues in the diagonal matrix Λ and the eigenvectors in the corresponding columns of the matrix V , then calculate: Then the integral consists of elements 1 through n of , and consists of elements n through n + m of . The product V − 1 θ M is calculated by solving the system Vx = θ M for x , which does not require calculating V − 1 . If the matrix M has repeated eigenvalues, the above procedure fails when there is an incomplete set of linearly independent eigenvectors with which to construct a nonsingular matrix V . In this case, our code attempts to calculate generalized eigenvectors to replace the duplicated eigenvectors in the matrix V . If v is an eigenvector of M associated with a repeated, degenerate eigenvalue λ , then a generalized eigenvector g solves the equation ( M − λ I ) g = v . To calculate a solution for g numerically, we first compute the Moore-Penrose inverse, also known as the pseudoinverse, ( M − λ I ) + , of the matrix ( M − λ I ). The Moore-Penrose inverse B + of a matrix B satisfies BB + B = B . It can be shown that g = B + v is a solution to Bg = v . If we start with the equation to be solved and multiply both sides on the left by BB + , we get BB + Bg = BB + v . Using the definitional property of B + , the left-hand side is equivalent to Bg , so Bg = B ( B + v ), which shows that g = B + v solves the equation. We use a numerical method for the singular value decomposition of the matrix ( M − λ I ) to calculate ( M − λ I ) + and then calculate a generalized eigenvector g = ( M − λ I ) + v . Then, we construct a Jordan matrix J to replace Λ and a matrix P containing eigenvectors and generalized eigenvectors such that M = PJP −1 . Then, the solution , and the component K ( J ) includes elements above the diagonal involving derivatives of the function K evaluated at the repeated eigenvalue(s). Acknowledgements The authors acknowledge Prabasaj Paul for his role in inspiring this work and for helpful insights during its development. References 1. ↵ World Health Organization . Antibiotic Resistance: Global Report on Surveillance . Geneva : WHO Press ; 2014 . 2. ↵ Won SY , Munoz-Price LS , Lolans K , Hota B , Weinstein RA , Hayden MK , et al. Emergence and rapid regional spread of Klebsiella pneumoniae carbapenemase-producing Enterobacteriaceae . Clin Infect Dis . 2011 ; 53 ( 6 ): 532 – 40 . OpenUrl CrossRef PubMed 3. ↵ Slayton RB , Toth D , Lee BY , Tanner W , Bartsch SM , Khader K , et al. Vital Signs: Estimated Effects of a Coordinated Approach for Action to Reduce Antibiotic-Resistant Infections in Health Care Facilities - United States . MMWR Morb Mortal Wkly Rep . 2015 ; 64 ( 30 ): 826 – 31 . OpenUrl CrossRef PubMed 4. ↵ Toth DJA , Khader K , Slayton RB , Kallen AJ , Gundlapalli AV , O’Hagan JJ , et al. The potential for interventions in a long-term acute care hospital to reduce transmission of carbapenem-resistant Enterobacteriaceae in affiliated healthcare facilities. Clin Infect Dis . 2017 . 5. ↵ Toth DJA , Samore MH , Nelson RE . Economic Evaluations of New Antibiotics: The High Potential Value of Reducing Healthcare Transmission Through Decolonization . Clin Infect Dis . 2021 ; 72 ( Suppl 1 ): S34 – S41 . OpenUrl PubMed 6. ↵ Diekmann O , Heesterbeek JA , Metz JA . On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations . J Math Biol . 1990 ; 28 ( 4 ): 365 – 82 . OpenUrl CrossRef PubMed Web of Science 7. ↵ Hethcote HW . The mathematics of infectious diseases . SIAM Review . 2000 ; 42 ( 4 ): 599 – 653 . OpenUrl CrossRef PubMed 8. ↵ Cooper BS , Medley GF , Scott GM . Preliminary analysis of the transmission dynamics of nosocomial infections: stochastic and management effects . J Hosp Infect . 1999 ; 43 ( 2 ): 131 – 47 . OpenUrl CrossRef PubMed Web of Science 9. Austin DJ , Bonten MJ , Weinstein RA , Slaughter S , Anderson RM . Vancomycin-resistant enterococci in intensive-care hospital settings: transmission dynamics, persistence, and the impact of infection control programs . Proc Natl Acad Sci U S A . 1999 ; 96 ( 12 ): 6908 – 13 . OpenUrl Abstract / FREE Full Text 10. Lipsitch M , Bergstrom CT , Levin BR . The epidemiology of antibiotic resistance in hospitals: paradoxes and prescriptions . Proc Natl Acad Sci U S A . 2000 ; 97 ( 4 ): 1938 – 43 . OpenUrl Abstract / FREE Full Text 11. ↵ Cooper BS , Stone SP , Kibbler CC , Cookson BD , Roberts JA , Medley GF , et al. Systematic review of isolation policies in the hospital management of methicillin-resistant Staphylococcus aureus: a review of the literature with epidemiological and economic modelling . Health Technol Assess . 2003 ; 7 ( 39 ): 1 – 194 . OpenUrl CrossRef PubMed 12. ↵ Cooper BS , Medley GF , Stone SP , Kibbler CC , Cookson BD , Roberts JA , et al. Methicillin-resistant Staphylococcus aureus in hospitals and the community: stealth dynamics and control catastrophes . Proc Natl Acad Sci U S A . 2004 ; 101 ( 27 ): 10223 – 8 . OpenUrl Abstract / FREE Full Text 13. Bootsma MC , Diekmann O , Bonten MJ . Controlling methicillin-resistant Staphylococcus aureus: quantifying the effects of interventions and rapid diagnostic testing . Proc Natl Acad Sci U S A . 2006 ; 103 ( 14 ): 5620 – 5 . OpenUrl Abstract / FREE Full Text 14. Bootsma MC , Wassenberg MW , Trapman P , Bonten MJ . The nosocomial transmission rate of animal-associated ST398 meticillin-resistant Staphylococcus aureus . Journal of the Royal Society, Interface / the Royal Society . 2011 ; 8 ( 57 ): 578 – 84 . OpenUrl 15. Gurieva T , Bootsma MC , Bonten MJ . Successful Veterans Affairs initiative to prevent methicillin-resistant Staphylococcus aureus infections revisited . Clin Infect Dis . 2012 ; 54 ( 11 ): 1618 – 20 . OpenUrl CrossRef PubMed 16. Donker T , Wallinga J , Slack R , Grundmann H. Hospital networks and the dispersal of hospital-acquired pathogens by patient transfer . PLoS One . 2012 ; 7 ( 4 ): e35002 . OpenUrl CrossRef PubMed 17. ↵ Hetem DJ , Westh H , Boye K , Jarlov JO , Bonten MJ , Bootsma MC . Nosocomial transmission of community-associated methicillin-resistant Staphylococcus aureus in Danish Hospitals . J Antimicrob Chemother . 2012 ; 67 ( 7 ): 1775 – 80 . OpenUrl CrossRef PubMed Web of Science 18. ↵ McBryde ES , Pettitt AN , McElwain DL . A stochastic mathematical model of methicillin resistant Staphylococcus aureus transmission in an intensive care unit: predicting the impact of interventions . J Theor Biol . 2007 ; 245 ( 3 ): 470 – 81 . OpenUrl CrossRef PubMed Web of Science 19. ↵ Rubin MA , Samore MH , Harris AD . The Importance of Contact Precautions for Endemic Methicillin-Resistant Staphylococcus aureus and Vancomycin-Resistant Enterococci . JAMA . 2018 ; 319 ( 9 ): 863 – 4 . OpenUrl CrossRef PubMed 20. ↵ Mangalea MR , Halpin AL , Haile M , Elkins CA , McDonald LC . Decolonization and Pathogen Reduction Approaches to Prevent Antimicrobial Resistance and Healthcare-Associated Infections . Emerg Infect Dis . 2024 ; 30 ( 6 ): 1069 – 76 . OpenUrl PubMed 21. ↵ Hayden MK , Lin MY , Lolans K , Weiner S , Blom D , Moore NM , et al. Prevention of colonization and infection by Klebsiella pneumoniae carbapenemase-producing enterobacteriaceae in long-term acute-care hospitals . Clin Infect Dis . 2015 ; 60 ( 8 ): 1153 – 61 . OpenUrl CrossRef PubMed 22. ↵ Munoz-Price LS . Long-term acute care hospitals . Clin Infect Dis . 2009 ; 49 ( 3 ): 438 – 43 . OpenUrl CrossRef PubMed Web of Science 23. ↵ Toth DJA , Khader K , Beams A , Samore MH . Model-based Assessment of the Effect of Contact Precautions Applied to Surveillance-detected Carriers of Carbapenemase-producing Enterobacteriaceae in Long-term Acute Care Hospitals . Clin Infect Dis . 2019 ; 69 ( Suppl 3 ): S206 – S13 . OpenUrl PubMed 24. ↵ Zimmerman FS , Assous MV , Bdolah-Abram T , Lachish T , Yinnon AM , Wiener-Well Y. Duration of carriage of carbapenem-resistant Enterobacteriaceae following hospital discharge . Am J Infect Control . 2013 ; 41 ( 3 ): 190 – 4 . OpenUrl CrossRef PubMed 25. ↵ DalBen MF , Teixeira Mendes E , Moura ML , Abdel Rahman D , Peixoto D , Alves Dos Santos S , et al. A Model-Based Strategy to Control the Spread of Carbapenem-Resistant Enterobacteriaceae: Simulate and Implement . Infect Control Hosp Epidemiol . 2016 ; 37 ( 11 ): 1315 – 22 . OpenUrl PubMed 26. ↵ Adler A , Navon-Venezia S , Moran-Gilad J , Marcos E , Schwartz D , Carmeli Y. Laboratory and clinical evaluation of screening agar plates for detection of carbapenem-resistant Enterobacteriaceae from surveillance rectal swabs . J Clin Microbiol . 2011 ; 49 ( 6 ): 2239 – 42 . OpenUrl Abstract / FREE Full Text 27. ↵ Mathers AJ , Poulter M , Dirks D , Carroll J , Sifri CD , Hazen KC . Clinical microbiology costs for methods of active surveillance for Klebsiella pneumoniae carbapenemase-producing Enterobacteriaceae . Infect Control Hosp Epidemiol . 2014 ; 35 ( 4 ): 350 – 5 . OpenUrl CrossRef PubMed 28. ↵ Marino S , Hogue IB , Ray CJ , Kirschner DE . A methodology for performing global uncertainty and sensitivity analysis in systems biology . J Theor Biol . 2008 ; 254 ( 1 ): 178 – 96 . OpenUrl CrossRef PubMed Web of Science 29. ↵ Lin MY , Lyles-Banks RD , Lolans K , Hines DW , Spear JB , Petrak R , et al. The importance of long-term acute care hospitals in the regional epidemiology of Klebsiella pneumoniae carbapenemase-producing Enterobacteriaceae . Clin Infect Dis . 2013 ; 57 ( 9 ): 1246 – 52 . OpenUrl CrossRef PubMed 30. ↵ Wittekamp BHJ , Oostdijk EAN , Cuthbertson BH , Brun-Buisson C , Bonten MJM . Selective decontamination of the digestive tract (SDD) in critically ill patients: a narrative review . Intensive Care Med . 2020 ; 46 ( 2 ): 343 – 9 . OpenUrl CrossRef PubMed 31. ↵ Gencay YE , Jasinskyte D , Robert C , Semsey S , Martinez V , Petersen AO , et al. Engineered phage with antibacterial CRISPR-Cas selectively reduce E. coli burden in mice . Nat Biotechnol . 2024 ; 42 ( 2 ): 265 – 74 . OpenUrl CrossRef PubMed 32. ↵ Tosh PK , McDonald LC . Infection control in the multidrug-resistant era: tending the human microbiome . Clin Infect Dis . 2012 ; 54 ( 5 ): 707 – 13 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted February 23, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about medRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Transmission thresholds for the spread of infections in healthcare facilities Message Subject (Your Name) has forwarded a page to you from medRxiv Message Body (Your Name) thought you would like to see this page from the medRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Transmission thresholds for the spread of infections in healthcare facilities Damon J.A. Toth , Karim Khader , Christopher Mitchell , Matthew H. Samore medRxiv 2025.02.21.25322698; doi: https://doi.org/10.1101/2025.02.21.25322698 Share This Article: Copy Citation Tools Transmission thresholds for the spread of infections in healthcare facilities Damon J.A. Toth , Karim Khader , Christopher Mitchell , Matthew H. Samore medRxiv 2025.02.21.25322698; doi: https://doi.org/10.1101/2025.02.21.25322698 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Epidemiology Subject Areas All Articles Addiction Medicine (568) Allergy and Immunology (863) Anesthesia (300) Cardiovascular Medicine (4435) Dentistry and Oral Medicine (444) Dermatology (382) Emergency Medicine (608) Endocrinology (including Diabetes Mellitus and Metabolic Disease) (1509) Epidemiology (15227) Forensic Medicine (30) Gastroenterology (1124) Genetic and Genomic Medicine (6597) Geriatric Medicine (668) Health Economics (997) Health Informatics (4534) Health Policy (1368) Health Systems and Quality Improvement (1613) Hematology (540) HIV/AIDS (1264) Infectious Diseases (except HIV/AIDS) (15916) Intensive Care and Critical Care Medicine (1103) Medical Education (623) Medical Ethics (146) Nephrology (667) Neurology (6599) Nursing (346) Nutrition (998) Obstetrics and Gynecology (1144) Occupational and Environmental Health (957) Oncology (3332) Ophthalmology (974) Orthopedics (369) Otolaryngology (420) Pain Medicine (436) Palliative Medicine (130) Pathology (663) Pediatrics (1693) Pharmacology and Therapeutics (691) Primary Care Research (711) Psychiatry and Clinical Psychology (5447) Public and Global Health (9230) Radiology and Imaging (2198) Rehabilitation Medicine and Physical Therapy (1370) Respiratory Medicine (1196) Rheumatology (593) Sexual and Reproductive Health (712) Sports Medicine (530) Surgery (712) Toxicology (99) Transplantation (289) Urology (265) (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'a0033593dcdf3fe2',t:'MTc3OTUzMDM4MQ=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();