A Multicompartment SIS Stochastic Model with Zonal Ventilation for the Spread of Nosocomial Infections: Detection, Outbreak Management, and Infection Control.

In this work, we study the environmental and operational factors that influence airborne transmission of nosocomial infections. We link a deterministic zonal ventilation model for the airborne distribution of infectious material in a hospital ward, with a Markovian multicompartment SIS model for the infection of individuals within this ward, in order to conduct a parametric study on ventilation rates and their effect on the epidemic dynamics. Our stochastic model includes arrival and discharge of patients, as well as the detection of the outbreak by screening events or due to symptoms being shown by infective patients. For each ventilation setting, we measure the infectious potential of a nosocomial outbreak in the hospital ward by means of a summary statistic: the number of infections occurred within the hospital ward until end or declaration of the outbreak. We analytically compute the distribution of this summary statistic, and carry out local and global sensitivity analysis in order to identify the particular characteristics of each ventilation regime with the largest impact on the epidemic spread. Our results show that ward ventilation can have a significant impact on the infection spread, especially under slow detection scenarios or in overoccupied wards, and that decreasing the infection risk for the whole hospital ward might increase the risk in specific areas of the health-care facility. Moreover, the location of the initial infective individual and the protocol in place for outbreak declaration both form an interplay with ventilation of the ward.

INTRODUCTION

The risk of acquiring nosocomial infections is a recognized problem in health‐care facilities worldwide (Harbarth, Sax, & Gastmeier, 2003). While the transmission routes for some diseases are well documented, the precise mode of transmission is uncertain for many infections, particularly for those pathogens that cause health‐care–acquired infections (HCAIs). Although it is probable that the majority of transmission occurs via contact routes (Sax et al., 2009), there is increasing recognition that the air plays an important role in disease spread (Li et al., 2007). Understanding the role that ventilation airflow plays in the dispersion of infectious microorganisms is tantamount to assessing exposure to pathogens and hence infection risk. This study aims to provide an analytical link between airborne hospital infection spread, ventilation design, and outbreak management.

Ventilation has been found to have a significant impact on the distribution of infectious material in hospital settings. Examples include Influenza A (e.g. H5N1 and H7N9) (Reed & Kemmerly, 2009), Mycobacterium tuberculosis (Escombe et al., 2010), measles (Atkinson, Chartier, Pessoa‐Silva, Jensen, & Li, 2009), and norovirus (Teunis et al., 2008). One of the most infamous examples occurred in 2003 during the severe acute respiratory syndrome (SARS) outbreak in Hong Kong. Analysis of airflow patterns and outbreak data demonstrated that ventilation routes were critical in the short‐ and long‐range spread of aerosolized coronavirus (Li, Huang, Yu, Wong, & Qian, 2005). Ventilation is recognized as an important infection control approach in health‐care design, with strategies such as mechanical ventilation and pressure zoning set out in international (Atkinson et al., 2009) and national guidance (Department of Health, 2007).

Evaluating the influence of ventilation on infection risk typically applies models such as the Wells–Riley equation (Riley, Murphy, & Riley, 1978) or a dose–response approach (Sze & Chao, 2010) to estimate the influence of ventilation on the number of new cases of an infection. Liao, Chang, and Liang (2005) presented a probabilistic transmission dynamic model to assess indoor airborne infection risks and Ko, Burge, Nardell, and Thompson (2001) and Ko, Thompson, and Nardell (2004) developed models for tuberculosis spread incorporating a zonal ventilation model. A number of authors have also looked at control strategies, including Wein and Atkinson (2009), who modeled infection control measures for pandemic influenza; Brienen, Timen, Wallinga, Van Steenbergen, and Teunis (2010), who analyzed the effect of mask use on the spread of influenza; and King, Noakes, and Sleigh (2015), who developed a stochastic model to link airborne and contact transmission. It is also worth mentioning the recent work by Carruthers et al. (2018), where a zonal ventilation model similar to the one considered in this article is linked to a dose–response approach to estimate the risk of infection after an accidental release of bacteria Francisella tularensis in a microbiology laboratory.

While these studies enable some understanding of the influence of the environment on transmission, they do not consider relationships between ventilation parameters and the progression and control of an infection outbreak. In an earlier study, it was demonstrated that the Wells–Riley model could be coupled to an epidemic model to relate ventilation rate and transmission in a fully mixed environment (Noakes, Beggs, Sleigh, & Kerr, 2006). In later work, a zonal air distribution and a stochastic formulation (Noakes & Sleigh, 2009) was considered, and cost benefits of ventilation from an energy and infection risk perspective were explored (Noakes, Sleigh, & Khan, 2012).

The model presented in this article is constructed on a scenario defined in previous work (Noakes & Sleigh, 2009), where the role played by the airflow during a nosocomial outbreak is assessed by linking a deterministic zonal ventilation model with an stochastic epidemic model using a computational approach. While the previous approach enabled exploration of the basic interaction between the ventilation and the outbreak, there are a number of limitations: We refer the reader to Keeling and Ross (2007) where the limitations of analyzing this type of epidemic processes by simulation are discussed in more detail, and where the benefits of following exact analytical approaches instead are highlighted.

The epidemic dynamics are represented through a simple epidemic model, not accounting for relevant factors such as the discharge and admission of patients, or the detection and declaration of the outbreak.

Results reported by Noakes and Sleigh (2009) have high variability, which is related to the fact that they were obtained by means of stochastic simulations of the epidemic process.

The large number of parameters associated with each ventilation scenario makes it difficult to identify, from stochastic simulations, the specific factors of the ventilation air distribution that facilitate or mitigate epidemic spread.

Our aim here is to show how this zonal ventilation model can be linked to more complex stochastic epidemic models for the spread of nosocomial pathogens, while accounting for patients, admission and discharge, and different outbreak detection and declaration hypotheses. We show how to implement exact analytical procedures for computing summary statistics of the outbreak (statistics measuring outbreak infectiousness), which by means of a perturbation analysis enables identification of specific characteristics of the ventilation setting that are crucial for the spread or control of the infection.

Finally, we carry out a comprehensive numerical study of six ventilation strategies for a hypothetical hospital ward in order to identify particular ventilation characteristics that may promote or inhibit spread of airborne nosocomial infections. Our results explore the interplay between ward ventilation, location of patients, ward overoccupancy, and outbreak detection management.

THE MODEL

Noakes and Sleigh (2009) propose a stochastic model that links a zonal ventilation model with epidemic dynamics by means of an SI model for the infection spread among individuals in a hospital setting. The combined model can be thought of as an adaptation of the Wells–Riley model (Riley et al., 1978; Sze & Chao, 2010), where each susceptible individual is infected with a per capita rate, that is proportional to the number I of infective individuals. Here, Q is the room ventilation rate (m3/s), p is the pulmonary ventilation (breathing) rate (m3/s), and q is the unit of infection (quantum) as described by Riley et al. (1978). In a standard stochastic epidemic model, where individuals are either susceptible or infective, individuals do not recover. This can be expressed in terms of a continuous‐time Markov chain (CTMC) , where the random variable represents the number of susceptible individuals at time t, and the aim is to compute transient probabilities for any , where N is the number of individuals in the population. These probabilities satisfy the Kolmogorov differential equations:

Interevent times are assumed to be exponentially distributed and the stochastic process is simulated using the Gillespie algorithm (Gillespie, 1976). In particular, this algorithm generates samples of the interevent times by using the fact that

and then updates the state of the system according to the probabilities of each possible event occurring.

A Zonal Ventilation Model for Linking Airflow Dynamics and Infection Rates

In Noakes and Sleigh (2009), this Wells–Riley process is adapted to investigate ventilation scenarios in a hypothetical hospital ward split in M ventilation zones. The air is assumed to be uniformly mixed within each zone; however, there is incomplete mixing between the zones and unbalanced zone boundaries allow for the effect of directional flow to be examined. In particular, the per capita infection rate for susceptible individuals at zone k is defined as , where is the concentration of infectious material at zone k and is the pulmonary rate of these individuals. We note that this concentration could depend on the number of infected individuals, , in every zone , due to airflow. In Noakes and Sleigh (2009), the spatial distribution of infectious material is represented through the differential equation: where is the volume of zone k, is the generation rate of infectious quanta, represents the extract ventilation rate in zone k, and and amount to the volume flow of air between zones k and j. Moreover, each interzonal flow rate represents the sum of two contributions: where β0 is a global mixing rate and is an additional contribution representing net flow across the zonal boundary, from zone k to zone j.

Equation (2) leads to a ventilation matrix that characterizes the ventilation air distribution under study, representing ventilation in a hospital ward divided into M ventilation zones.

By assuming steady‐state conditions for airflow, and taking into account Equation (2), Noakes and Sleigh (2009) propose to link infection rates , with the ventilation matrix V as follows: This means that per capita infection rates for susceptible individuals at zones depend on how many infective individuals there are in any zone at any given time, computed in a specialized manner (Equation (3)) that takes into account the ventilation distribution through matrix V. Once this procedure for computing infection rates is proposed, stochastic simulations for the epidemic dynamics are carried out in Noakes and Sleigh (2009) by following steps (i)–(iv) (Noakes & Sleigh, 2009), which assume exponentially distributed interevent times and make use of the property depicted in Equation (1).

We note here that Equation (3) means considering per capita infection rates at each zone as functions of the number of infectives within each ventilation zone at the hospital ward. In Section 2.2, we now go on to exploit this, to link the zonal ventilation model with a multicompartment model with detection, to evaluate the infection spread dynamics within the hospital ward until detection of the outbreak. Instead of carrying out stochastic simulations, we present an exact approach for analyzing a summary statistic of the outbreak: the total number of infections occurring until the outbreak ends or is detected and declared. This exact approach does not only allow us to compute this quantity of interest, but it also allows one to carry out a sensitivity analysis on the model parameters, so that the impact that different characteristics of the ventilation setting has on this summary statistic can be evaluated.

A Multicompartment Stochastic Model for the Infection Spread Dynamics

At the epidemic level, we assume that patients in each zone i are discharged at rate , so that amounts to the average length of stay (LOS) of patients in zone i. Discharges are immediately replaced by new admitted patients, a reasonable approximation for hospital wards under high demand (Pelupessy, Bonten, & Diekmann, 2002; Wolkewitz, Dettenkofer, Bertz, Schumacher, & Huebner, 2008). Moreover, we consider that the nosocomial outbreak will go undetected by health‐care workers for some time, and incorporate this fact into our model by considering that each infected individual in zone i can be discovered/detected at some rate . The reciprocal represents the average time until some symptoms arise that alert health‐care workers to a patient's infection, or the average time until the infected individual is detected through screening policies put in place at this hospital ward. Fig. 1 represents the epidemic dynamics for an individual in zone i.

Figure 1

Individual epidemic dynamics for a patient in zone i. Event represents the discharge of the infective patient in zone i, immediately replaced by a susceptible one.

This leads to a multicompartment SIS epidemic model that can be described as a CTMC , where represents the number of infective individuals in zone j at time , defined over the space of states . is the total number of patients in zone j, leading to patients in the hospital ward. State Δ represents that the nosocomial outbreak has been detected and declared by health‐care workers by the first detection of an infected patient in the hospital ward. We note that absorbing state (0, …, 0) represents the end of the outbreak (lack of infective individuals), due to patients' discharge (i.e., if all the patients infected by the pathogen are discharged before the outbreak is actually detected). We consider Δ also as an absorbing state in this process, since we are only interested in the dynamics of the process until the end or declaration of the outbreak, and the transitions (obtained from diagram in Fig. 1) described in Table I. We note that, according to our comments in Section 2.1, is a function of the state , representing the per capita infection rate of susceptible individuals in zone j when we have infective individuals within the ward, computed from Equation (3) for each .

Table I

Transitions or Events Occurring in the Stochastic Process , and Corresponding Rates

Event	Stochastic Transition	Rate
Infection of a patient in zone j	(i_1,⋯,i_j,⋯,i_M)→(i_1,⋯,i_j+1,⋯,i_M)	λ_j(i_1,⋯,i_M)(N_j−i_j)
Discharge of an infective patient in zone j	(i_1,⋯,i_j,⋯,i_M)→(i_1,⋯,i_j−1,⋯,i_M)	γ_ji_j
Detection of the outbreak	(i_1,⋯,i_M)→Δ	∑j=1Mδ_ji_j

METHODOLOGY

Our interest is in analyzing the infectious potential of an outbreak in a given hospital ward, for different ventilation configurations and outbreak detection hypotheses. We measure the infectious potential in terms of the following discrete random variable: which can take values in . R is used here as a measure of how well ventilation can act as a preventive infection control strategy that is in place even before the actual detection of the outbreak by health‐care workers occurs and reactive strategies can be implemented.

Analyzing R

For an initial state , our aim is to compute probabilities for ; that is, the probability distribution of R for some initial state . We can compute these probabilities from a system of linear equations, which is obtained by a first‐step argument. In particular, by proposing notation we get for and any , and with boundary conditions , for all ; a detailed explanation on how Equations (4) and (5) are obtained is in the Appendix. This means that probabilities for (Equation (4)) can be computed by solving a system of linear equations.1 Once these are in hand, probabilities for can be computed by solving the system of linear equations given by Equation (5), which also consists of  equations.

Algorithm 1 in the Appendix computes probabilities for any , . It works sequentially, computing probabilities for for increasing values of and .

Local Sensitivity Analysis

Our analysis allows one to identify the most important characteristics of the ventilation scenario, regarding the infectious potential of the outbreak until detection, by means of computing partial derivatives of the form with respect to ventilation parameters . We note that, for an initial state , , so that Partial derivatives can be computed from direct differentiation of Equations (4) and (5). In Equations (4) and (5), the only quantities that depend on parameter are infection rates and probabilities . Thus, we get for , and any . Partial derivatives can then be computed from the equations above by following arguments similar to those in Algorithm 1 in the Appendix. In order to solve these equations, one needs to have in hand values of (previously computed from Algorithm 1), as well as derivatives . These derivatives can be straightforwardly obtained from Equation (3) as:

where represents the element‐by‐element partial derivative of matrix V with respect to parameter θ (Gómez‐Corral & López‐García, 2018).

Spread until the Dth Individual Detection

As outlined above, declaration of the outbreak is identified with the first detection of an infective patient, where each patient is detected in zone j at rate . If detection of an infective patient occurs because this patient shows symptoms, outbreak declaration might require several () patients showing some common symptoms, since for some nosocomial pathogens, associated symptoms are quite common and pass unnoticed (Ekkert, 2015). For example, norovirus causes gastrointestinal symptoms such as nausea, vomiting, or diarrhea that are common to many diseases and conditions. The National Guidelines on the Management of Outbreaks of Norovirus Infection in Health‐Care Settings (National Disease Surveillance Centre, 2003), issued by the National Disease Surveillance Centre in Ireland, requires for patients to show these symptoms in a hospital ward for a potential norovirus outbreak declaration. Once the outbreak has been declared, control strategies such as immediate cleaning and decontamination, frequent handwashing, or cohorting of affected patients are recommended.

Thus, our interest in this subsection is to analyze the summary statistic R when the detection of the outbreak requires D patients to show symptoms, for some value , and results in the subsections above can be seen as the particular case . We define the augmented process , where the increasing variable amounts to the number of detected patients up to time . We consider that the outbreak is declared once , and the space of states of this CTMC is given by: Thus, , with state Δ representing outbreak declaration (i.e., the detection of the Dth infected patient). Events occurring in this process, at different rates, are described in Table II.

Table II

Transitions or Events Occurring in the Stochastic Process , and Corresponding Rates

Event	Stochastic Transition	Rate
Infection of a patient in zone j	(i_1,⋯,i_j,⋯,i_M,d)→(i_1,⋯,i_j+1,⋯,i_M,d)	λ_j(i_1,⋯,i_M)(N_j−i_j)
Discharge of an infective patient in zone j	(i_1,⋯,i_j,⋯,i_M,d)→(i_1,⋯,i_j−1,⋯,i_M,d)	γ_ji_j
Detection of an infective patient	(i_1,⋯,i_M,d)→(i_1,⋯,i_M,d+1), 0≤d≤D−2	∑j=1Mδ_ji_j
Outbreak declaration	(i_1,⋯,i_M,D−1)→Δ	∑j=1Mδ_ji_j

Our arguments in the subsections above can be adapted for process . For example, Equations (4) and (5) become for any , and with boundary conditions , and for all and . 1 represents a function that takes value 1 if A is satisfied, and 0 otherwise. An adapted version of Algorithm 1, not reported here for the sake of brevity, allows for an efficient solution of this system.

RESULTS

We consider here the hypothetical hospital ward in Fig. 2 in Noakes and Sleigh (2009), and ventilation settings according to Table 1 in Noakes and Sleigh (2009). This hypothetical hospital ward consists of three six‐bedded bays connected through a corridor. Each bay is split into two ventilation zones, and each ventilation zone contains three patients. The corridor is split into three ventilation zones and in the first instance has no patients, so that the ward has patients. Ventilation rate over the whole ward was , which equated to an air change rate of 3 AC/h. Diagrams of the ventilation flows and specific rates within and between each zone for each setting are given in Fig. 2. Moreover, we set and for all patients at all zones and for every ventilation scenario (Noakes & Sleigh, 2009), and assume an average LOS for each patient  days.

Figure 2

Diagrams of the ventilation settings considered, and corresponding ventilation matrices V (for ). Gray arrows: ventilation between zones (, m3/min); black arrows: ventilation supply and extract to the ward; circled values: extract ventilation rates (, m3/min)

If we order the ventilation zones as the ventilation matrix is given by

and ventilation settings in Table 1 in Noakes and Sleigh (2009) lead to the ventilation matrices reported in Fig. 2.

In Sections 4.1–4.3 and 4.6, we consider that outbreak declaration occurs after one patient shows symptoms, with each patient showing symptoms after an average time (i.e., for all ). Alternative outbreak detection and declaration hypotheses are discussed in Sections 4.4 and 4.5, while the impact of parameter q in our numerical results is explored in Section 4.7.

Impact of Ventilation Setting on Spread Dynamics

In Fig. 3, we plot the probability mass function of R versus different values of the global mixing rate β0, the average time at which each infective patient shows symptoms, and for ventilation scenarios A–F. For these results, it is assumed that an infective patient in zone 1a starts the outbreak, and we report in Table III the mean values E[R] computed for these distributions.

Figure 3

Probability mass function of the number R of infections until the end or detection of the outbreak, for ventilation settings A–F. Declaration of the outbreak occurs once one patient shows symptoms, and each infective patient shows symptoms after an average time . Global mixing rate . Initial infective located in zone 1a. Red stars identify the best ventilation settings in terms of E[R], while black circles identify the worst ones; see Table III.

Table III

Mean Number E[R] of Infections until End or Declaration of the Outbreak, for Scenarios in Fig. 3

β0	δ^−1	A	B	C	D	E	F
9	12 h	1.75	1.54	1.24	1.75	1.74	1.50
	48 h	5.01	4.34	3.87	5.20	4.86	4.50
27	12 h	1.88	1.76	1.61	1.86	1.87	1.77
	48 h	5.59	5.25	5.03	5.59	5.53	5.37

Ventilation setting C can be identified in Fig. 3 and Table III as the best one, while settings A, D, and E are identified as the worst ones depending on the detection parameter δ and the global mixing rate β0. We note that ventilation setting C has significant extract ventilation at the initially infected zone 1a, so that the airflow is directed from 1a outward to the hospital ward. On the other hand, ventilation setting D represents a well‐mixed ward ( for all i and k) with no extract ventilation at zone 1b (), which might favor the spread of pathogens from 1a toward other zones within the ward, leading to more infections occurring until outbreak detection.

In general terms, worse scenarios can be identified for values h and , where the long‐tailed distribution of R for these scenarios in Fig. 3 indicates that large outbreaks occur with significant probability. Larger differences among ventilation settings are also found for value h. Thus, our results suggest that ventilation of the ward should be of special concern for pathogens that have longer infectious asymptomatic periods, or in hospital wards with more limited surveillance policies. It is also clear that the average individual detection time has a higher impact on the infection spread than the specific ventilation setting in the ward, so that outbreak detection seems to dominate ventilation regarding infection spread.

Dependence on Location of Initial Infective

In Fig. 4, we plot analogous results to those in Fig. 3 when the infective patient starting the outbreak is located in zones , for h and . Corresponding mean values E[R] are reported in Table IV. We note that zones 3a and 3b are equivalent to zones 1a and 1b, for all ventilation settings in Fig. 2, and thus we do not test them. For zones near the corridor (i.e., 1b and 2b), ventilation setting B is identified as the best one, while D is identified as the worst one. We note that ventilation setting B has no extract ventilation in zones 1b and 2b (), but it directs the airflow instead toward corridor areas. In this setting B, corridor areas have no patients and significant ventilation, with airflow unbalance from bays to corridor areas acting in practice as an infection control measure. Thus, our results suggest that the spread control ability of a given ventilation setting depends on the location of the patient starting the outbreak as well as the airflow direction. However, from results in Fig. 4 and Table IV, ventilation setting D seems to perform poorly regardless of the initial infective location, suggesting that some ventilation settings might be inadvisable regardless of this location (i.e., if this location is unknown).

Figure 4

Probability mass function of the number R of infections until the end or detection of the outbreak, for ventilation settings A–F. Declaration of the outbreak occurs once one patient shows symptoms, and each infective patient shows symptoms after an average time h. Initial infective located in zones . Global mixing rate

Table IV

Mean Values E[R] for Distributions in Fig. 4, That Is, for Different Locations of the Initial Infective

Location	A	B	C	D	E	F
1a	1.75	1.54	1.24	1.75	1.75	1.50
1b	1.74	1.28	1.57	1.93	1.50	1.73
2a	1.79	1.57	1.27	1.79	1.78	1.54
2b	1.78	1.32	1.62	1.98	1.54	1.78

When focusing on a particular location for the initial infective, comments above are supported by the sensitivity analysis on the ventilation parameters. For example, in Tables V and VI we report, for ventilation parameters , partial derivatives , and elasticities for ventilation settings B and D, h, and an infective patient starting the outbreak in zone 1b. We note that while dimensionless elasticities are useful for comparison purposes, they equal zero if parameter θ is zero.

Table V

Partial Derivatives and Elasticities of E[R] with Respect to Ventilation Parameters, for Setting B

Partial Derivatives						Elasticities
β0	0.0230	β_Q2b,c2	−0.0141	Q_o,1b	−0.0890	β0	0.2899	β_Q2b,c2	−0.1042	Q_o,1b	0
β_Q1a,1b	−0.0201	β_Qc2,2b	0.0243	Q_o,c1	−0.0398	β_Q1a,1b	−0.1483	β_Qc2,2b	0	Q_o,c1	−0.2876
β_Q1b,1a	0.0207	β_Qc2,c3	0.0019	Q_o,2a	−0.0209	β_Q1b,1a	0	β_Qc2,c3	0	Q_o,2a	0
β_Q1b,c1	−0.0382	β_Qc3,c2	−0.0008	Q_o,2b	−0.0274	β_Q1b,c1	−0.2735	β_Qc3,c2	0	Q_o,2b	0
β_Qc1,1b	0.0301	β_Q3a,3b	−0.0032	Q_o,c2	−0.0221	β_Qc1,1b	0	β_Q3a,3b	−0.0240	Q_o,c2	−0.1615
β_Qc1,c2	0.0097	β_Q3b,3a	0.0060	Q_o,3a	−0.0120	β_Qc1,c2	0	β_Q3b,3a	0	Q_o,3a	0
β_Qc2,c1	−0.0041	β_Q3b,c3	−0.0081	Q_o,3b	−0.0159	β_Qc2,c1	0	β_Q3b,c3	−0.0601	Q_o,3b	0
β_Q2a,2b	−0.0056	β_Qc3,3b	0.0141	Q_o,c3	−0.0132	β_Q2a,2b	−0.0420	β_Qc3,3b	0	Q_o,c3	−0.0971
β_Q2b,2a	0.0106	Q_o,1a	−0.1089			β_Q2b,2a	0	Q_o,1a	0

Note: Average individual detection time h, global mixing rate , and an initially infective patient starting the outbreak in zone 1b.

Table VI

Partial Derivatives and Elasticities of E[R] with Respect to Ventilation Parameters, for Setting D

Partial Derivatives						Elasticities
β0	−0.0032	β_Q2b,c2	−0.0198	Q_o,1b	−0.1746	β0	−0.0151	β_Q2b,c2	0	Q_o,1b	0
β_Q1a,1b	0.0184	β_Qc2,2b	0.0241	Q_o,c1	−0.0872	β_Q1a,1b	0	β_Qc2,2b	0	Q_o,c1	−0.1356
β_Q1b,1a	−0.0248	β_Qc2,c3	0.0046	Q_o,2a	−0.0411	β_Q1b,1a	0	β_Qc2,c3	0	Q_o,2a	−0.1278
β_Q1b,c1	−0.0212	β_Qc3,c2	−0.0026	Q_o,2b	−0.0718	β_Q1b,c1	0	β_Qc3,c2	0	Q_o,2b	0
β_Qc1,1b	0.0124	β_Q3a,3b	0.0053	Q_o,c2	−0.0590	β_Qc1,1b	0	β_Q3a,3b	0	Q_o,c2	−0.0918
β_Qc1,c2	0.0167	β_Q3b,3a	−0.0074	Q_o,3a	−0.0293	β_Qc1,c2	0	β_Q3b,3a	0	Q_o,3a	−0.0912
β_Qc2,c1	−0.0090	β_Q3b,c3	−0.0140	Q_o,3b	−0.0516	β_Qc2,c1	0	β_Q3b,c3	0	Q_o,3b	0
β_Q2a,2b	0.0074	β_Qc3,3b	0.0171	Q_o,c3	0.0434	β_Q2a,2b	0	β_Qc3,3b	0	Q_o,c3	−0.0675
β_Q2b,2a	−0.0102	Q_o,1a	−0.1017			β_Q2b,2a	0	Q_o,1a	−0.3165

Note: Average individual detection time h, global mixing rate , and an initially infective patient starting the outbreak in zone 1b.

Regime B requires airflow toward the corridor in order to expel pathogens from zone 1b, since . Thus, rates , , , and correspond to significantly large negative elasticities reported in Table V (i.e., increasing the values of these rates would lead to decreasing values of E[R]). Global mixing rate β0 has a significant impact (large positive elasticity) favoring disease spread, since increasing the value of β0 represents increasing the rate at which pathogens flow among all zones, instead of flowing specifically toward the extract ventilation areas (corridors in this setting).

According to results in Table VI, ventilation setting D could be significantly improved by increasing extract ventilation (especially in zones 1a, 1b, and c1), as well as increasing airflow from 1b to 1a and to c1. This is directly related to the fact that, since there is no extract ventilation in zone 1b, infectious material in this zone can only be expelled by directing it toward adjacent zones 1a and c1.

Decreasing Hospital Ward Infection Spread Risk Might Increase Risk at Specific Bays

It is clear that the number, R, of infections occurring until the end or detection of the outbreak can be split according to where these infections actually occur as where is the number of infections occurring at bay j. Although probabilities can be analytically computed by adapting arguments in Section 3, details are omitted here for the sake of brevity, and results in Table VII are obtained from 106 stochastic simulations of the process.

Table VII

Mean Values of E[R], , , and for , h and Ventilation Settings A, D, and E; Initially Infective Individual in Zone 1a

Regime	E[R]	E[R(1)]	E[R(2)]	E[R(3)]
A	5.01	2.65	1.34	1.02
D	5.20	2.56	1.47	1.17
E	4.86	2.69	1.25	0.92

In Table VII, we report values of for , h, ventilation settings A, D, and E, and an infective patient starting the outbreak in zone 1a. Results suggest that epidemic spread can be limited by switching ward ventilation from setting D to A, and further containment is obtained by switching to ventilation setting E. However, infection risk in bay 1 (in terms of ) behaves contrarily; although the global hospital ward infection risk (in terms of E[R]) is lower for setting E, this is at the expense of expelling pathogens from the infected zone 1a toward zones 1b and c1, and thus posing a greater risk to patients in bay 1.

Interplay with Detection Management

As explained in Section 3.3, detection and declaration of an outbreak in the hospital ward may require several patients showing symptoms, and not only one. In Fig. 5, we plot analogous results to those in Fig. 3 when declaration of the outbreak occurs after patients show symptoms (each after average time ). The corresponding mean values of E[R] are reported in Table VIII.

Figure 5

Probability mass function of the number R of infections until the end or detection of the outbreak, for ventilation settings A–F. Declaration of the outbreak occurs once two patients show symptoms, and each infective patient shows symptoms after an average time . Global mixing rate . Initial infective located in zone 1a.

Table VIII

Mean Values of E[R] for Scenarios in Fig. 5

β0	δ^−1	A	B	C	D	E	F
9	12 h	3.20	2.78	2.33	3.26	3.14	2.79
	48 h	8.42	7.31	6.78	8.79	8.10	7.72
27	12 h	3.52	3.28	3.07	3.50	3.49	3.34
	48 h	9.37	8.87	8.63	9.39	9.26	9.10

We note that values in Table VIII are significantly larger than those in Table III, since outbreak declaration takes longer to occur, allowing for more infections to take place. This increase is significantly larger than the differences that can be observed, in Table III, between different ventilation settings, suggesting again that detection policy is likely to dominate ventilation as an infection control strategy. Under slow detection scenarios ( h), we observe in Fig. 5 a clear bimodality for the distribution of R. Thus, our model predicts that under slow detection, a two‐output situation can be expected: either the initially infective patient is discharged before infecting any other patient (so that ), or this patient infects a second patient, leading to a large outbreak (represented by the second mode in Fig. 5).

Screening at Admission

In Sections 4.1–4.4, we analyze infection spread under the assumption that each individual in zone j is detected (by showing symptoms) at rate , with for all j, and where the outbreak is detected and declared after one (or several) infective patients are detected. This leads to the contribution in Equations (4) and (5). However, if the detection of the outbreak is due instead to the screening of the newly admitted patient who starts the outbreak, and results of this screening arrive after an average time , then the outbreak is detected at rate δ, and one needs to replace by δ in Equations (4) and (5).

Under this hypothesis, we plot in Fig. 6 the probability mass function of the number R of infections until the end or detection of the outbreak, when the results of this screening (and thus, the declaration of the outbreak) arrive after an average time . Corresponding mean values of E[R] are reported in Table IX. If results arrive after h, ventilation has a less significant impact on the nosocomial spread, and low values of E[R] are reported in Table IX. The number of infections until the end or detection of the outbreak proportionally increases with the delay in obtaining the screening results. In particular, for h significant differences in E[R] can be noticed among the different ventilation settings, and a marked bimodality can be observed in Fig. 6.

Figure 6

Probability mass function of the number R of infections until the end or detection of the outbreak, for ventilation settings A–F. Declaration of the outbreak occurs after an average time , which is independent of the number of infectives and represents a screening event. Global mixing rate . Initial infective located in zone 1a.

Table IX

Mean Values E[R] for Scenarios in Fig. 6

δ^−1	A	B	C	D	E	F
4 h	1.12	0.92	0.66	1.11	1.11	0.89
8 h	2.66	2.13	1.63	2.72	2.59	2.16
12 h	4.24	3.43	2.74	4.36	4.11	3.54
24 h	8.15	6.91	5.94	8.36	7.91	7.17

Ventilation and Overoccupancy

In this subsection, our aim is to shed some light on the interplay between ventilation, nosocomial spread, and overoccupancy of the hospital ward. We represent hospital ward overoccupancy by locating three additional patients at the corridor areas; in particular, we set and locate one additional patient in each of the zones. This practice is common in U.K. hospitals during times of high demand. We assume that the outbreak is detected and declared after the first patient shows symptoms, each patient showing symptoms after an average time h. For interzonal mixing , we report in Table X the mean number R of infections until the end or detection of the outbreak, for ventilation settings A–F, and for the initial infective patient being located in zones .

Table X

Mean Values E[R] for a Hospital Ward with Overoccupancy

Location of Initial Infective	A	B	C	D	E	F
1a	2.01	1.82	1.37	2.04	2.02	1.72
c1	1.95	1.14	1.95	2.19	1.71	1.92
c2	1.99	1.18	2.00	2.23	1.75	1.96
c3	1.95	1.14	1.95	2.19	1.71	1.92

Note: Patients show symptoms after an average time h, and we set .

We first note that results in Table X are exactly the same for an initially infective individual being located in zones c1 and c3. This is explained by noting that bays 1 and 3 are completely symmetric for all ventilation settings as noted in Section 4.2; see diagrams and matrices in Fig. 2. On the other hand, when the initial infective patient is located in zone 1a, results in Table X can be compared to those in Table III for . For zone 1a, values of E[R] are larger in Table X than in Table III; if the nosocomial outbreak is initiated by an infective patient in zone 1a, more infections during this outbreak should be expected under ward overoccupancy. This might not be only related to having more patients in the ward under overoccupancy (21 instead of 18), but also to the potential of patients in the corridor to act as infection links between bays. For example, under overoccupancy, an infective individual in zone 1a might infect individuals in bay 2 by, as a first step, infecting individuals in the corridor areas. These people then might more easily infect individuals in bay 2, before being discharged, due to being in closer proximity and depending on the particular ventilation setting in place in the ward.

Infection dynamics related to the scenario above highly depend on the particular ventilation setting under study, which can be noticed by inspecting rows corresponding to zones in Table X. While an individual in zone 1a has a larger infectious potential (in terms of E[R]) than individuals located in the corridor when ventilation settings B and E are in place, this is not the case for ventilation settings C, D, and F, and these infectious potentials are comparable under ventilation setting A, which represents a well‐mixed scenario. Our results then indicate that overoccupancy leads in general to higher airborne spread risks, and that this increase can be especially significant depending on the specific ventilation in place.

The Unit of Infection

We note that parameter q is highly pathogen dependent, ranging from quanta/min for rhinovirus to quanta/min for measles (Noakes & Sleigh, 2009). We perform a parametric analysis by varying q to assess the sensitivity of our conclusions, and report expected infections E[R] until the end of the outbreak in Table XI. In particular, we are interested in the mean number E[R] of infections if the outbreak is detected on the first patient showing symptoms, each patient showing symptoms after an average time h, and where we consider m3/min. We note that, as expected, increasing values of q lead to increasing mean number E[R] of infections. However, this does not seem to affect the relative infectiousness of ventilation setting C, which is identified as the best scenario regardless of the value of q. On the other hand, less advantageous ventilation schemes are dependent on the value of q. For example, ventilation setting A and E can be identified as the worst for quanta/min, while setting D can be seen as the worst one for quanta/min.

Table XI

Mean Number E[R] of Infections until End or Declaration of the Outbreak, for m3/min and When Outbreak Declaration Occurs After the First Patient Shows Symptoms

q	A	B	C	D	E	F
0.1	0.41	0.36	0.27	0.38	0.42	0.34
0.5	1.75	1.54	1.24	1.75	1.74	1.50
1.0	3.06	2.66	2.26	3.14	3.00	2.69
10.0	11.43	10.44	10.16	11.75	11.13	10.92
50.0	15.54	15.12	15.05	15.66	15.41	15.35

Note: Each patient shows symptoms after an average time h. Different values of q (quanta/min) considered.

DISCUSSION AND CONCLUSIONS

In this work, we link a zonal ventilation model for the generation and airborne spread of infectious material within a hospital ward with a multicompartment SIS Markovian model for the infection of patients within this ward. Our model incorporates the possibility of considering a wide range of ventilation settings, the discharge and arrival of patients within the ward, as well as different hypotheses regarding how outbreak detection and declaration occurs. Moreover, it allows us to explore the interplay between ventilation, outbreak management, ward overoccupancy, and the location of the infective patient starting the outbreak.

Our results suggest that detection time dominates ventilation when the variable of interest is the number of infections occurring before the declaration or end of the outbreak, with longer detection times leading to significantly more infections happening. Longer detection times can arise when analyzing pathogens with long infectious asymptomatic periods, when declaration of an outbreak requires for several patients to show symptoms, or when this declaration depends on screening events for which results take longer to arrive. The interplay between ventilation of the hospital ward and location of the initially infective patient starting the outbreak implies that recommendations on where to locate potentially infected (e.g., newly admitted) patients in a given hospital ward could be issued depending on the ventilation in place in the ward. Our model also predicts that decreasing the infection spread risk in the hospital ward can sometimes come at the expense of increasing the risk in particular areas of the ward.

Similar models have already been considered in the literature for linking zonal ventilation scenarios with epidemic spread models (Carruthers et al., 2018; Ko et al., 2004; Noakes & Sleigh, 2009), where epidemic dynamics are usually analyzed by means of stochastic simulations. To the best of our knowledge, this is the first time that this link is carried out by defining in detail the CTMC for the infection spread, where infection rates at each ventilation zone are in fact functions of the number of infectives in each zone at any given time, and where represents in fact a state of the CTMC under study. This detailed mathematical construction allows for the analytical computation of summary statistics (such as R in this work), and for carrying out a local sensitivity analysis that allows one to identify the particular factors of each ventilation setting having the most significant impact on the infection spread.

It should be noted that the primary objective in this study is to demonstrate this detailed mathematical analysis and how it can be applied to evaluate the relative influence of different parameters. The model is applied to a hypothetical hospital ward, which, while it is representative of multibed ward environments in many hospitals, is a very simplified model of reality. The results demonstrate that the ventilation flow settings may influence the dispersion of airborne pathogens and hence the risk of transmission; however, these should be interpreted with caution. We assume a steady‐state ventilation scenario with the flow pattern replicated exactly between neighboring bays 1, 2, and 3. In reality, the flows will not be exactly identical for every bay, and other factors such as heat sources and movement of people will alter the mixing with and between zones. In particular, corridor ventilation often has a directional flow due to wider spacing of ventilation supply/extract grilles, which may hinder or improve the control of infection. However, the analysis we have carried out gives some clear insight into why particular directional flows influence risk, and the relative importance of detection strategies, ventilation control, and occupancy.

It is clear that some of our conclusions could be highly dependent on the hospital ward structure, and therefore the flexibility of our methodology comes into play. It can be applied to any hospital ward of interest by appropriately adapting the corresponding ventilation matrix V. Although carrying out a detailed mathematical analysis of a number of potentially different hospital ward structures is out of the scope of this article, we include a short numerical study of an alternative hospital ward in the Supplementary Material. The aim of this is twofold: (i) to show how our methodology can be easily implemented for a different hospital ward to that in Fig. 2 by just adapting the ventilation matrix V, illustrating how this matrix varies with hospital ward structure; and (ii) to show that while some of our conclusions might be hospital ward structure dependent, others seem to be valid for a wide range of hospital ward structures (e.g., detection dominates ventilation as well in this alternative hospital ward).

In this article, we go beyond the SI epidemic model in Noakes and Sleigh (2009), proposing an SIS‐type model, which allows us to incorporate patients' arrival/discharge and outbreak detection and declaration. This is similar to the model recently proposed by López‐García and Kypraios (2018), as a unified framework for modeling the spread of nosocomial infections. We note that this epidemic model structure would be especially relevant for pathogens with no or short (i.e., negligible compared to the average patient's LOS) noninfectious or latent periods, and where the infectious period is long enough so that recovery of patients does not occur before discharge (or detection). Depending on the hospital ward under analysis and the average LOS of patients in this ward, this could be the case for influenza or norovirus. Pathogens with nonnegligible incubation periods (e.g., 7–21 days for measles) might require more complex stochastic epidemic model structures such as the SEIR (susceptible–exposed–infective–recovered). On the other hand, when analyzing hospital wards with longer average patient LOS, so that individuals may become infected and recover during their stay, SIRS‐type epidemic models would be required to represent the recovery of patients () before discharge () occurs. We note here that in principle, the methodology outlined in Section 3 can be extended to any of these compartmental‐based epidemic models for the disease spread dynamics, where the link between the deterministic zonal ventilation model for the airflow dynamics and the stochastic epidemic model for disease spread dynamics would still be as in Section 2. In a similar way, more complex epidemic model structures could allow one to study the infection spread dynamics after outbreak detection and declaration occurs. In this article, we have focused instead on the impact of ventilation on disease spread until the end or declaration of the outbreak. Considering these alternative compartmental‐based epidemic model structures could be the aim of future work.

Finally, we note that when carrying out our analysis, the main computational effort lies in solving systems of linear equations, where the number of equations is determined by the number of states of the corresponding CTMC, given by Equation (6) in our model. Limitations of our approach are then of a computational nature, since highly complex epidemic models (here, a multicompartment SIS stochastic model with detection) linked to large hospital wards split in many different ventilation zones (here, zones with three empty zones and six zones containing three patients each) would lead to an intractable number of equations, and stochastic simulation approaches would prevail.

DATA, SOFTWARE, AND REPRODUCIBILITY

Computer codes (in Python) in order to reproduce our numerical results are available López‐García, King, and Noakes (2019).

Fig. S1. Diagrams of the ventilation settings considered, and corresponding ventilation matrices V.

Fig. S2. Probability mass function of the number R of infections until the end or detection of the outbreak, for ventilation settings SA‐SD.

Table S1. Mean Values of E[R] for Scenarios in Fig. S2

Click here for additional data file.