Abstract
Multidrugresistant Klebsiella pneumoniae is an increasing cause of infant mortality in developing countries. We aimed to develop a quantitative understanding of the drivers of this epidemic by estimating the effects of antibiotics on nosocomial transmission risk, comparing competing hypotheses about mechanisms of spread, and quantifying the impact of potential interventions. Using a sequence of dynamic models, we analysed data from a oneyear prospective carriage study in a Cambodian neonatal intensive care unit with hyperendemic thirdgeneration cephalosporinresistant K. pneumoniae. All widelyused antibiotics except imipenem were associated with an increased daily acquisition risk, with an odds ratio for the most common combination (ampicillin + gentamicin) of 1.96 (95% CrI 1.18, 3.36). Models incorporating genomic data found that colonisation pressure was associated with a higher transmission risk, indicated sequence type heterogeneity in transmissibility, and showed that withinward transmission was insufficient to maintain endemicity. Simulations indicated that increasing the nursepatient ratio could be an effective intervention.
Introduction
Infections with multidrugresistant Enterobacteriaceae constitute a major threat to public health in all regions of the world (Schwaber and Carmeli, 2008; Theuretzbacher, 2017) and extended spectrum βlactamase (ESBL) producing and carbapenemresistant Enterobacteriaceae have been prioritised by the World Health Organization as pathogenic bacteria in need of novel therapeutics (World Health Organization, 2014; World Health Organization, 2017). These organisms pose the highest risk to subgroups of patients such as those undergoing surgery, requiring invasive devices, and neonates (Peleg and Hooper, 2010; Goldmann, 1981). Antimicrobial resistance is of particular concern in developing (lower and middleincome) countries where the estimated per capita mortality from drugresistant bacteraemia is far greater than in highincome countries and where lastline antibiotics may be unavailable or unaffordable (Lim et al., 2016; Zaidi et al., 2005). The high risk of difficulttotreat nosocomial infections threatens to undermine patient confidence in developing world hospitals and health systems (Dondorp et al., 2018).
In this paper, we focus on Klebsiella pneumoniae, which is amongst the most clinically important multidrugresistant Gramnegative bacteria in developing country settings (FoxLewis et al., 2018; Zellweger et al., 2017; Musicha et al., 2017). Genomic studies characterising the population structure of K. pneumoniae have revealed a complex consisting of three separate species (K. pneumoniae, K. quasipneumoniae, K. variicola) that are indistinguishable by culture or standard biochemical assays (Holt et al., 2015). Where these isolates remain undifferentiated by molecular assay, we refer to them as K. pneumoniae sensu lato (in the broad sense).
While there is currently widespread concern about Gramnegative bacteria as an emerging threat due to high levels of plasmidborne resistance (Nordmann et al., 2011), pathogens such as K. pneumoniae have been considered a major problem in nosocomial settings for over half a century (Yow, 1955). Research during the 1960s identified a number of drivers of Klebsiella colonisation and infection in hospital settings including invasive devices (Mertz et al., 1967), environmental contamination (Kresky, 1964), introduction from the community (Kessner and Lepper, 1967), persontoperson transmission (Weil et al., 1966), endogenous selection from antibiotic pressure (Selden et al., 1971) and transient carriage on the hands of healthcare workers (Adler et al., 1970). In a comprehensive review of this body of work, Montgomerie concluded that ‘The likely means of transmission of Klebsiella is via the hands of hospital staff members’ (Montgomerie, 1979).
More recently, carriage studies in highincome hospital settings in temperate regions have used wholegenome sequencing to show the critical importance of asymptomatic carriage for understanding the epidemiology of K. pneumoniae, establishing a firm link between gastrointestinal carriage and clinical infection (Martin et al., 2016; Gorrie et al., 2017; Gorrie et al., 2018). In contrast to studies in highincome countries where multidrugresistant K. pneumoniae is typically rare, prospective carriage studies in hospitalised paediatric populations in developing countries in Africa and Asia have reported hyperendemic levels of ESBLproducing Enterobacteriaceae including K. pneumoniae (Andriatahina et al., 2010; Roberts et al., 2019; Founou et al., 2019; Turner et al., 2016) consistent with the greater burden of disease due to these organisms in lower income settings (Lim et al., 2016; Musicha et al., 2017).
Despite their clinical importance, there are major gaps in our knowledge of the epidemiology of multidrugresistant K. pneumoniae. First, while a number of recent investigations of localised hospital outbreaks in highincome settings have provided evidence that longterm environmental contamination of sinks and other sites may play a role (Clarivet et al., 2016; Decraene et al., 2018; Mathers et al., 2018), the relative importance of persistently contaminated point sources versus patienttopatient transmission in endemic settings remains unclear. Second, though it is widely assumed that antibiotic exposures play an important role in selecting for multidrugresistant K. pneumoniae (Baker et al., 2018), such effects have not previously been quantified at the patient level in a way that disentangles antibiotic effects from general exposure to the hospital environment. Third, attempts to quantify transmissibility and how this varies by sequence type (ST) are lacking. Fourth, the impacts of other factors that might plausibly affect transmission including staffing levels, infant breast feeding and use of probiotics have not been explored.
To address these knowledge gaps, we used data collected from a yearlong prospective observational carriage study in a Cambodian neonatal intensive care unit and analysis methods which build upon a previously described datadriven stochastic model (Forrester and Pettitt, 2005). We fit four models with logit link functions to estimate the impact of covariates on the daily risk of acquisition of 3GCR K. pneumoniae s.l.. As these models are unable to identify the force of infection, and genomic data show the ward 3GCR K. pneumoniae s.l. to be a highly heterogeneous bacterial community, we then fit five linear transmission models to estimate the force of infection for acquisition of each ST. Finally, we use the estimated parameters to provide modelbased assessments of the potential impact of hypothetical control measures. We define "acquisition" as the first detection of the organism within an infant following an initial negative swab on admission, which may indicate transmission from other colonised infants or hospital staff, or endogenous selection as a result of antibiotic pressure.
Results
Descriptive epidemiological data
Over the yearlong observation period, there were consistently high rates of patient carriage of third generation cephalosporinresistant (3GCR) K. pneumoniae sensu lato. Of 333 infants admitted to the neonatal unit, 121 of 289 (42%) were found to be colonised on the first swab taken within 48 hours of ward admission. A further 21 out of 44 (48%) were positive on the first swab that was taken more than 48 hours after admission. Overall, 109/191 (57%) infants who initially screened negative for 3GCR K. pneumoniae s.l. became positive during their stay in the neonatal unit. Almost all 3GCR K. pneumoniae s.l. isolates were ESBL producers (1412/1423; 99%), and only 5/1423 (0.35%) were resistant to imipenem. Cocolonisation with 3GCR E. coli was observed in 52 infants on their first swab, and a further 102 infants became cocolonised with both resistant organisms during their stay on the neonatal unit. Full details on the study population, including bloodstream infections and mortality have been reported previously (Turner et al., 2016), and a summary is provided in Table 1.
The daily counts of infants known to be colonised with 3GCR K. pneumoniae s.l. (Figure 1A), shows no clear trend but large stochastic fluctuations, which are expected given the small size of the ward (eight beds) and frequent discharges of patients and introduction of colonised infants (imported carriage). A representation of swabbing interval outcomes, as used in the models, is shown in Figure 1B. While the median length of stay was four days, the distribution is highly skewed with a tail of longstaying patients (Figure 1C). The most frequently used antibiotic combination was ampicillin with gentamicin which was used empirically to treat suspected sepsis in infants admitted from the community and was taken on one fifth (19%) of patient days on the ward. This was followed by imipenem, which was taken on 12% of patient days. Imipenem was typically used when culture results showed nonsusceptibility to firstline antibiotic choices and empirically in infants with suspected hospitalacquired infection. All other antibiotic combinations were used at a much lower frequency (Figure 1D).
Factors associated with carriage acquisition
Infants were prospectively screened for the organism whilst on the ward by culture of rectal swabs on selective media. For infants with negative cultures for 3GCR K. pneumoniae s.l. on admission to the ward, the outcome (acquisition of the organism) was therefore interval censored between subsequent rectal swabs that were taken a median of 2 days apart (IQR 1, 3 days). In total there were 400 swab outcomes (either negatives or a first positive swab) over 864 patient days from 191 infants with a negative culture at entry.
Four models were fitted to the interval censored swab data to determine factors associated with daily risk of carriage acquisition. The best performing single intercept model (model A; Table 2) with the lowest WAIC (see Methods) considered exposure to antibiotics in the previous 96 hours and did not include a term for colonisation pressure (i.e. daily perpatient acquisition risk did not depend on the number of other patients who were colonised on a given day). Models with a 48 hour antibiotic exposure period (model B) and those that included a colonisation pressure term (model C) showed slightly worse fits.
The covariates associated with reduced daily risk of acquisition were breast milk feeding (odds ratio [OR] 0.69 [95% CrI 0.35, 1.41]) and increasing the number of nurses, for instance three nurses in the ward was associated with an OR of 0.55 (95% CrI 0.15, 1.77) relative to zero nurses (baseline). Male sex was also associated with a reduced risk of acquisition (OR 0.68 [95% CrI 0.43, 1.04]). Other covariates such as taking a probiotic (OR 0.88 [95% CrI 0.55, 1.41]), a severe condition (OR 1.10 [95% CrI 0.55, 2.13]), prior colonisation with 3GCR E. coli (OR 1.07, [95% CrI 0.65, 1.75]), and age at admission (10 days compared with zero days, OR 0.95 [95% CI 0.73, 1.24]) had ORs distributions centred closer to unity. As anticipated (though, to our knowledge, not previously shown), antibiotics taken (intravenously, with the exception of cloxacillin) within the past 96 hours were mostly associated with an increased risk of colonisation with 3GCR K. pneumoniae s.l.. Ampicillin (OR 1.77 [95% CrI 0.88, 3.26]), ampicillin with gentamicin (OR 1.96 [95% CrI 1.18, 3.36]), ceftriaxone (OR 1.85 [95% CrI 0.68, 4.54]), oral cloxacillin (OR 1.49 95% CrI 0.47, 4.02]), and cloxacillin with gentamicin (OR 1.94 [95% CrI 0.55, 5.66]) were all associated with an increased risk of acquisition. Only intravenous imipenem (OR 1.01 [0.40, 2.30]) had a posterior distribution centred near unity, consistent with the carbapenem sensitivity of 3GCR K. pneumoniae found in this setting. See Figure 2A for odds ratio posterior distributions.
Using the covariate posterior distributions, we also estimated the probability of colonisation for each of the 864 patient days where patients were at risk for acquiring 3GCR K. pneumoniae s.l.. The median daily probability of first acquisition of 3GCR K. pneumoniae s.l. for an infant was estimated from the best fitting model as 0.15. There is considerable variability in the risk of acquisition between patient days and the medians of the posterior probability distribution ranges nearly eightfold from 0.047 to 0.35 (Figure 2B).
The risk of becoming colonised with 3GCR K. pneumoniae s.l. is cumulative over an infant’s length of stay and varies in response to interventions, such as consumption of antibiotics. We show the cumulative risk of first acquisition under two scenarios: 1) where an initially four day old, breastfed, female infant remains in the ward for eight days without taking antibiotics or probiotics; and 2) where an infant with the same characteristics is prescribed ampicillin with gentamicin from day three onwards. The median cumulative risk of acquiring 3GCR K. pneumoniae s.l. after eight days for the first scenario is 0.57 (80% CrI 0.42, 0.71) and for the second scenario is 0.75 (80% CrI 0.64, 0.85). Although the median cumulative risk between the two scenarios diverges the longer the infants are in the neonatal unit, the uncertainty also increases with time (Figure 2C).
Swab sensitivity
We estimated the sensitivity of rectal swabs for detecting 3GCR K. pneumoniae s.l. by examining the swabs that followed a positive from the same patient. There were 936 such swabs which were taken from a patient after at least one swab positive for 3GCR K. pneumoniae s.l., and 90 (9.6%) of these were negatives. Under the assumption that all negatives following a positive culture are false negatives, the false negative rate posterior median was 0.096 (95% CrI 0.078, 0.12) and the posterior median swab sensitivity was 0.90 (95% CrI 0.88, 0.92). Under the second assumption that three or more consecutive negative swabs following a positive culture represent a true decolonisation event, there were 72 false negatives, giving a false negative rate of 0.073 (95% CrI 0.058, 0.091) and a swab sensitivity of 0.93 (95% CrI 0.91, 0.94).
Model assessment and comparison
The measures of Markov chain convergence showed high effective sample sizes (>400) and $\hat{R}$<1.01, indicating that the chains had run for long enough and had mixed well (see Methods and Figure 2—figure supplement 1). Model assessment was performed with a posterior predictive check; we estimated the probability of acquisition for each of the 400 swabbing intervals and binned these probabilities into groups defined by the quintiles. We then calculated the expected number of colonisation events in each of the five groups and compared these with the observed number of acquisitions. Within each of these groups, the posterior median of the predicted number of acquisitions was close to the observed number of events, and the observed values were always within the 80% CrI of the model estimates (Figure 2D). The results from fitting risk factor model A with alternative prior distributions are shown in Figure 2—figure supplement 2; substantially reducing the variance of the priors for the model intercept and covariates had a negligible effect on the posterior parameter estimates.
When the intercept was permitted to vary by study month in a hierarchical model (risk factor model D; Table 2), little variation was observed between months; the median posterior baseline probability ranged from 0.20 to 0.23 with wide credible intervals. As these models did not include a colonisation pressure term, the intercept incorporated timevarying changes in the underlying intensity of transmission. The low variance in the monthly intercepts therefore suggests a relatively constant force of infection over the 12month study period. In models where we included a colonisation pressure term (risk factor model C; Table 2) this was found to have a slightly negative slope for acquisition of 3GCR K. pneumoniae s.l. (OR 0.96 [95% CrI 0.86, 1.09]). This is surprising, as if patienttopatient transmission was occurring, we would usually expect the force of infection to increase with the colonisation pressure (Bonten, 2012). The finding therefore suggested that one of the following three possibilities was true: i) patienttopatient transmission was not occurring at a high frequency in this ward; ii) patienttopatient transmission was occurring but, because of the continually high wardlevel prevalence, variations in the force of infection could not be identified; iii) patienttopatient transmission was occurring but exposure to the presence of two or more colonised patients presented a similar risk for acquisition as exposure to one. We therefore used K. pneumoniae s.l. wholegenome sequence data to help determine the most plausible scenario.
Klebsiella wholegenome assemblies
We examined wholegenome assemblies of 317 3GCR K. pneumoniae s.l. isolates cultured from rectal or environmental swabs in the neonatal unit over a four month period (see Methods). A phylogeny based on kmer distances between assemblies is shown in Figure 3A. Of note is the highly diverse and structured nature of the pathogen population, in contrast to one dominated by a clonal expansion of a single lineage. Overall 62 distinct sequence types were identified in our collection of isolates. The species identified from culture as K. pneumoniae s.l. consists of three distinct subpopulations that meet the criteria for separate species. We isolated all three species from infants in the cohort (K. pneumoniae n = 219, K. quasipneumoniae subspecies similipneumoniae n = 95, K. variicola n = 3), and found diversity similar to that observed in a global collection of K. pneumoniae isolates (Holt et al., 2015), suggesting that the diversity accumulated within a Cambodian neonatal unit over four months is comparable to the diversity of K. pneumoniae globally. Many STs were characterised by only a single carriage isolate, suggestive of importations that were not subsequently transmitted to other patients (Figure 3B). The STs with the largest number of carriage isolates were ST334 (K. quasipneumoniae; n = 29), ST101 (K. pneumoniae; n = 21), ST1074 (K. pneumoniae; n = 21) and ST45 (K. pneumoniae; n = 17). The most frequent bla_{ESBL} genes in the wholegenome assemblies were CTXM15 (201/317; 63%), followed by CTXM14 (43/317; 13%) and CTXM63 (38/317; 12%).
Transmission models for sequence types
We fitted mechanistic models representing different transmission processes to the ST swab data. Within the four month period where sequence data were available, there were 171 events for first acquisition of a 3GCR K. pneumoniae s.l. ST among 150 infants. Among transmission models 1–3 that were initially tested, the model with the best fit to data by WAIC was transmission model 2 (see Table 3 and Methods), which has an intercept (α), representing a constant risk of acquisition, and a slope (β) which scales the risk of acquisition for each infant colonised with a given ST in the ward (i.e. it accounts for colonisation pressure). The model estimated the values for α as 0.0019 (95% CrI 0.0015, 0.0023) and β as 0.0097 (95% CrI 0.0075, 0.012). We then fitted transmission model 2 with a random effect term, where β was permitted to vary by ST and the underlying distribution of β was assumed to follow a beta distribution, with shape hyperparameters α and β (transmission model 4). We also fitted transmission model 2 with a random effect term where α was permitted to vary by ST (transmission model 5).
Transmission model 1 has a constant colonisation pressure by ST which was not linked to the daily number of colonised individuals with that ST; this model showed a substantially worse fit to the data by WAIC (Table 3). Transmission model 3, which included terms for colonisation pressure and contamination in the hospital environment also failed to improve the model fit, and the high estimate for λ (4.7 95% CrI [2.0, 7.0]) suggests that contamination left by previously colonised infants decays rapidly to background levels, with an estimated environmental half life of 3.6 hours (95% CrI 2.4, 7.6 hours; Table 4). Varying the α parameter by ST (transmission model 5) resulted in a substantially worse fit to data, suggesting there was not enough information in the model to differentiate STspecific background rates of colonisation.
The central estimates of the force of infection by ST from transmission model 4 are shown in Figure 3C, with estimates from the four most frequent STs highlighted in colour. The uncertainty around these parameter estimates for the four STs is shown in Figure 3D. The daily intervalcensored colonisation incidence for the most frequent STs are clustered in time, generally emerging and reaching extinction in the ward within a matter of weeks, suggestive of importation and subsequent patienttopatient transmission. The incidence and estimated force of infection for the four most frequent STs over a period in the study where all 3GCR K. pneumoniae s.l. isolates were sequenced are shown in Figure 3E and F. Parameter estimates for all transmission models are shown in Table 3. Model fitting diagnostics showed that Markov chains had converged satisfactorily (see Methods and Figure 3—figure supplement 1).
Sequence Type SNP diversity
We mapped reads from isolates in two STs with the highest estimated force of infection (ST45 and ST101; see Figure 3) to ST consensus reference genomes. The 21 ST101 carriage isolates had a mean read depth ranging from 25x to 125x (median 57x), and the 17 ST45 carriage isolates had a mean read depth ranging from 20x to 67x (median 56x). All the carriage isolates in both STs had >90% of the genome covered by >5x coverage. We then called and filtered SNPs (see Methods) to determine if the relatedness of carriage isolates within STs was consistent with recent persontoperson transmission. The pairwise number of variants in ST101 isolates between infants ranged from 15 to 68 SNPs (median 30), which was comparable to the variation seen within infants in ST101 (from 18 to 38 SNPs; median 28). Similarly in ST45 the pairwise SNP differences between infants varied from 13 to 223 (median 125), which was comparable to the withinhost ST diversity (from 55 to 212 SNPs; median 124). Therefore the SNP diversity observed within and betweenhosts was very similar for both ST45 and ST101.
Simulations with agentbased models
We used the posterior parameter estimates obtained from model fitting for forward simulations using a dynamic agentbased model in order to evaluate the potential impact of interventions (see Methods). We first estimated the wardlevel reproduction number (R_{A}) for 3GCR K. pneumoniae s.l. by simulating the introduction of a single colonised patient into a ward of eight susceptible patients. For all patients, length of stay was sampled from the empirical length of stay distribution (Figure 1C) and colonised patients had a transmission potential sampled from the posterior hyperparameter distribution from transmission model 4 (Table 3). The median of the R_{A} distribution was 0.65% and 95% of values fell between 0.36 and 1.09. The distribution of R_{A} values is shown in Figure 4A.
We then simulated the impact of interventions to reduce the rate of 3GCR K. pneumoniae s.l. acquisition by combining parameter estimates for colonisation pressure from transmission model 4 with the marginal effect of modifiable covariates from risk factor model A (Table 2). In the first intervention scenario, we varied the proportion of infants given an oral probiotic, in addition to varying the proportion of infants that were colonised on entry (imported cases; 5% or 40%). We used as an outcome the proportion of infants susceptible to 3GCR K. pneumoniae s.l. on admission that remained uncolonised on discharge. When the importation rate was high (40% colonised on entry; similar to our study population, see Table 1), setting the proportion of infants taking the probiotic to be 0%, 50% or 100% resulted in the median proportion remaining uncolonised as 0.54 (95% interval 0.34, 0.72), 0.56 (95% interval 0.38, 0.72) and 0.59 (95% interval 0.39, 0.75) respectively. In the lower importation setting (5% of infants colonised on entry), setting the proportion of infants taking the probiotic to be 0%, 50% or 100% resulted in the median proportion remaining uncolonised as 0.80 (95% interval 0.54, 0.93), 0.82 (95% interval 0.60, 0.93) and 0.84 (95% 0.61, 0.94), respectively.
In the second intervention scenario, we varied the proportion of breast milk fed infants between 25%, 50% and 90%, in addition to varying the proportion of infants that were categorised as imported cases. When the importation rate was high (40%), altering the proportion of infants breast fed between 25%, 50% and 90% resulted in the median proportion remaining uncolonised as 0.48 (95% interval 0.25, 0.69), 0.51, (95% interval 0.31, 0.70), and 0.56 (95% interval 0.38, 0.72), respectively. When the importation rate was low (5%) altering the proportion of infants breast fed between 25%, 50% and 90% resulted in the median proportion remaining uncolonised as 0.73 (95% interval 0.40, 0.92), 0.77 (95% interval 0.48, 0.92) and 0.84 (95% 0.61, 0.94) respectively.
In the third scenario, we simulated either three, four or eight nurses in the ward each day, corresponding to infant:nurse ratios of roughly 3:1, 2:1 and 1:1 respectively. Again, we examined this effect in settings with different importation rates. Varying the number of nurses in the ward between three, four and eight in the high importation setting resulted in median proportions of initially uncolonised infants who remained uncolonised throughout their neonatal unit stay of 0.54 (95% interval 0.36, 0.71), 0.61 (95% interval 0.39, 0.78) and 0.81 (95% interval 0.21, 0.97) respectively. In the lower importation setting, varying the infant:nurse ratio resulted in the median proportions remaining uncolonised of 0.81 (95% interval 0.56, 0.93), 0.86 (95% interval 0.62, 0.95) and 0.96 (95% interval 0.32, 0.99), respectively. Of all the simulated interventions therefore, increasing the number of nurses on the ward had the largest impact on reducing colonisation rates. The distributions of the outcome variables from all simulations are shown in Figure 4B, C and D.
Discussion
In a hyperendemic, developing country hospital setting, we analysed the transmission dynamics of 3GCR K. pneumoniae s.l. prospectively over one year. We compared the support for competing hypotheses about modes of spread, quantified effects of antibiotic exposures as drivers of the epidemic, evaluated risk factors for transmission, and forward simulated to evaluate the potential benefits of interventions.
We found that carriage strains of 3GCR K. pneumoniae s.l. among neonates constituted a highly diverse population, with considerable intrahost variation within STs. There were frequent acquisitions or detection of resistant K. pneumoniae s.l. that were closely related to strains carried by other infants on the ward at the time of acquisition (Figure 3E). Withinhost diversity of these lineages was similar to betweenhost diversity in potentially linked cases. Moreover, once genomic information was considered, models incorporating colonisation pressure as a risk factor for acquisition showed substantially better fits to data than models without colonisation pressure (Table 3). Taken together with the lack of persistent environmental contamination (Smit et al., 2018) and the lack of improvement in model fit when longterm environmental contamination was considered, these findings add support to the view that patienttopatient transmission (much of which is likely to be mediated by contacts with healthcare workers) is the main driver of resistant K. pneumoniae acquisition, and that colonised patients represent the primary reservoir. Notably, this result was only apparent when we incorporated genomic data into our models to estimate transmission parameters by ST.
We are aware of one other prospective study of ESBLproducing Enterobacteriaceae colonisation dynamics from a developing country (Bonneault et al., 2019). This study considered ESBLproducing Enterobacteriaceae colonisation in a neonatal unit in Madagascar over a period of six months, and found similar rates of acquisition with ESBLproducing Enterobacteriaceae to those seen here. Though sequencing data were not used, fitting transmission models to data provided evidence of patienttopatient and healthcare workerpatient transmission, particularly for ESBLproducing K. pneumoniae and the estimated daily transmission parameter (0.05; 0.008, 0.14) is similar to our estimates for the most common STs (Figure 3D).
Two other prospective carriage studies of K. pneumoniae have used whole genome sequencing to identify possible transmission events: while a one year study in an adult intensive care unit in Australia found five epidemiologically plausible intrahospital transmission chains (Gorrie et al., 2017), a one year carriage study in two geriatric wards (also in Australia) found no evidence of patienttopatient transmission chains (Gorrie et al., 2018). Other studies using genomic epidemiology to study resistant K. pneumoniae transmission have been performed in adult wards in highincome settings (Snitkin et al., 2012; Haller et al., 2015; Snitkin et al., 2017). While such studies have also provided support for patienttopatient transmission playing an important role, such retrospective investigations made in response to reported outbreaks of multidrugresistant K. pneumoniae might not reflect typical patterns of transmission. A study using proximity sensors to investigate transmission of ESBLproducing Enterobacteriaceae over four months in a long term care facility in France found stronger support for persontoperson transmission as the main route of acquisition for K. pneumoniae, though the evidence for persontoperson transmission of ESBLproducing E. coli was weaker (Duval et al., 2019). This suggests that separate mechanisms may drive the transmission of different ESBLproducing organisms and that this should be considered when analysing carriage data from multiple species of Enterobacteriaceae.
Under the assumption that patienttopatient transmission was driving the epidemic, we calculated that the rates of transmission within the ward were insufficient to maintain endemic transmission (i.e. the wardlevel reproduction number, R_{A}<1) (Cooper et al., 2004). These findings are comparable to a study which determined the relative force of infection between ESBLproducing E. coli and other ESBLproducing Enterobacteriaceae in 13 European intensive care units, finding that the latter (mainly K. pneumoniae) had a transmission rate almost three times greater than the former, and that the singleadmission reproduction number for both classes of organisms were well below one (Gurieva et al., 2018). Our central estimate of R_{A} is, however, substantially higher (0.65 compared with 0.17); this difference in estimated transmission potential may reflect differences in stafftopatient ratios, different standards of hygiene and infection control, different patterns of antibiotic use, and differences in the patient population (Dondorp et al., 2018). The findings from both studies indicate that repeated importation into the unit is needed to maintain endemicity of resistant Klebsiella. Imported cases may be acquired from other wards within the same hospital, other hospitals within the referral network, or community transmission. When rates of importation are high, as observed in this setting where up to 43% of infants may have been colonised on initial ward admission (Table 1), even effective interventions are limited in how many acquisition events they can prevent due to a high underlying colonisation pressure (Figure 4). The public health implications are therefore that control measures should be coordinated regionally and targeted to the wider hospital patient referral network (Donker et al., 2012; Ciccolini et al., 2013).
Amongst the most important findings was the consistent association between a patient’s antibiotic exposure and an increased risk of acquiring 3GCR K. pneumoniae s.l. or detection of the organism due to withinhost selection. To our knowledge, the role of antibiotics as drivers of carriage dynamics has not previously been explored with appropriate methods to account for timedependent antibiotic exposures in Gramnegative bacteria, but there are reasons for believing it is likely to be a key mechanism through which antibiotics select for ESBL K. pneumoniae (Tedijanto et al., 2018). With the exception of imipenem, for which there was no association with the risk of acquisition for these predominantly carbapenemsusceptible bacteria (1418/1423 K. pneumoniae s.l. isolates sensitive to imipemen; 99.6%), effect sizes were similar for different antibiotic combinations (median posterior OR ~2; Figure 2A). The narrower credible interval for the ampicillin + gentamicin combination (OR 1.96, 95% CrI 1.18, 3.36) reflects the much higher usage of this antibiotic combination compared to others (Figure 1D). These effects are consistent with hypotheses about antibiotic therapy leading to reduced microbiome diversity and subsequently increasing the risk of colonisation with drugresistant bacteria, which face less competition from fitter, sensitive strains (Lipsitch et al., 2000). Microbiome analysis has shown that the greatest dysbiosis following antibiotic therapy in healthy adults is four days after treatment starts (Palleja et al., 2018), lending confidence to our model comparison selecting a 96 hour exposure period over a 48 hour period (Table 2). While we cannot differentiate the effects of antibiotics on i) increasing the susceptibility of an infant to acquisition from another infant, versus ii) endogenous selection for resistant bacteria within that infant (Lipsitch and Samore, 2002) in the risk factor models, the small estimate for α relative to β in transmission model two suggests that background selection plays a relatively small role in acqusition compared with persontoperson transmission (Table 3).
Breast milk feeding was associated with a reduced risk of colonisation with 3GCR K.pneumoniae s.l. though uncertainty was large (Figure 2A). This accords with our understanding of the development of a healthy gut microbiome in the early stages of life, which can be adversely affected by the replacement of breast milk by formula (Bäckhed et al., 2015), and the protective effect of a diverse microbiome that competes against potentially pathogenic bacteria (Langdon et al., 2016). In this population, breastfeeding rates were high (90%) though our simulations showed that dropping the proportion to 50% or 25%, as seen in other developing world populations (Lauer et al., 2004), could increase the proportion becoming colonised during admission by around 5% and 8% respectively.
The finding that the oral probiotic Lactobacillus acidophilus was not strongly protective against acquisition was disappointing in light of earlier results that suggested a possible effect in slowing rates of acquisition of ESBLproducing Enterobacteriaceae (Turner et al., 2016). This negative result was shown by the odds ratio, which was located close to zero in the risk factor models for acquisition (Figure 2A) and supported by the forward simulations which showed only a 4% median decrease in infant colonisation rates when 100% of infants were prescribed probiotics (Figure 4B). Evidence is still limited about the value of probiotics for neonates. One large randomised trial in rural India reported a beneficial effect from a symbiotic preparation (combining a probiotic, Lactobacillus plantarum, with fructooligosaccharide) in preventing sepsis in infants (Panigrahi et al., 2017). Another randomised trial in preterm infants in England found no benefit from the probiotic Bifidobacterium breve BBG001 in preventing necrotising enterocolitis, blood culture positive sepsis or death before hospital discharge (Costeloe et al., 2016). There is also evidence that some probiotics given after antibiotic consumption can impair and delay the recovery of normal gut flora in humans (Suez et al., 2018).
An increased number of nurses on the ward was negatively associated with the risk of acquiring 3GCR K. pneumoniae s.l.. Our simulation studies using a dynamic agentbased model showed that increasing the nurse:infant ratio from 1:3 to 1:1 could reduce the number of infants becoming colonised in both high and low importation settings by about a quarter. These findings are consistent with results from large observational studies. For example, (Rogowski et al., 2013), in a retrospective cohort study in 67 neonatal units in the USA found a strong association between neonatal unit understaffing and an increased nosocomial infection rate (where understaffing was defined as a nursepatient ratio below US guidelines for the patient acuity level). Two distinct mechanisms might account for such an association. First, an imbalance between workload and staffing levels may lead to reduced attention to basic infection control measures such as hand hygiene, as has been reported in several studies (Pittet et al., 2006). Second, as nurse:patient ratios decrease a lower proportion of patient contacts will be cohorted as each nurse will need to contact more patients in a shift, substantially increasing the potential for crosstransmission (Archibald et al., 1997; Austin et al., 1999; Hugonnet et al., 2004).
From our model estimates, we observed considerable variation in the risk of acquisition per patient day, with the median posterior probability varying nearly eightfold from 0.047 to 0.35, showing that even neonates within a single hospital ward constitute a highly heterogeneous population. This challenges the assumptions of compartmental models for nosocomial transmission that treat patients as broadly homogeneous in their risk of acquiring drugresistant bacteria (Grundmann and Hellriegel, 2006; van Kleef et al., 2013; Domenech de Cellès et al., 2013).
The effect of colonisation pressure was not identifiable when we considered the total number of 3GCR K. pneumoniae s.l. isolates on a given day, but became identifiable when we considered transmission by ST. From our analysis of the Klebsiella carriage population, we observed three species (K. pneumoniae, K. quasipneumoniae and K. variicola) along with a considerable number (62) of STs. This indicates that Klebsiella in carriage did not reflect a single population, but rather repeated introductions of diverse isolates which were then either spread around the ward over a number of weeks or were not transmitted and became locally extinct.
Comparison of different transmission models strongly supported the inclusion of a colonisation pressure term to account for patienttopatient transmission. There was also support for a hierarchical model where transmissibility varied by ST, though the STs we estimated to have the highest transmissibility have not been highlighted as dominant in other settings (Wyres and Holt, 2016) suggesting that the ST composition within a region may reflect adaptation to local pressures (Stoesser et al., 2015b). The pairwise SNP diversity in two of the major STs (ST45 and ST101) was greater than expected, but the withinhost and betweenhost diversity was similar in both cases, suggesting that transmission within these STs is biologically plausible. Withinhost diversity has the potential to hinder the reconstruction of transmission networks (Worby et al., 2014; Didelot et al., 2016), and our results here highlight the importance of capturing withinhost diversity in sequencing studies. While withinhost diversity was particularly high in ST45, all acquisitions of this ST occurred within a 19 day window (over a possible four month period when carriage isolates were sequenced) with overlapping colonised patient stays, strongly suggesting that cases were epidemiologically related (Figure 3E). A previous analysis of a subset of the genomic data identified closely related clusters suggestive of transmission, though with a smaller number of pairwise SNPs than we observed (Smit et al., 2018). This discrepancy is likely due to differences in methodology as we mapped isolates to STspecific reference genomes, which results in a greater proportion of the genome being callable compared with mapping very diverse isolates to a single reference. We cannot entirely rule out less parsimonious explanations for the temporal clustering of STs, such as a transient increase of certain STs in the water supply, though the persontoperson transmission route, mediated by healthcare workers, is most strongly supported by our models.
An important strength of the study is that we considered asymptomatic carriage, rather than clinical isolates, and collected detailed patientlevel data from infants who became colonised as well as from infants who remained uncolonised. This allowed us to quantify the factors driving the epidemic, which would not have been possible if we had considered only clinical isolates. Also, by using a prospective design rather than a reactive exploration of an outbreak or unusual cluster of cases our findings should be more representative of typical patterns of transmission. The use of an inferrential approach that accounted for the intervalcensored nature of the data and the use of whole genome sequencing were also key factors in developing a quantitative understanding of the transmission dynamics.
Our study has limitations. We sequenced isolates of 3GCR K. pneumoniae s.l. from a four month period (rather than the full 12 months of the carriage study) due to resource constraints. Such constraints also prevented us from sampling mothers and other family members of infants, and healthcare workers in close contact with infants. Greater density of withinhost sampling (Stoesser et al., 2015a; Wymant et al., 2018; Lees et al., 2019), and inclusion of longread sequencing to investigate plasmid transmission (Conlan et al., 2014) would have also provided a more complete epidemiological picture. Although our analysis accounts for patient heterogeneity and interval censoring, it assumes that a positive culture taken within 48 hours of admission indicates that the patient was colonised on admission, and that cultures had 100% sensitivity for detecting carriage. Such restrictions could all be relaxed using a data augmentation framework (Cooper et al., 2008), but validated implementations of such an approach allowing inclusion of an arbitrary number of covariates are not currently available. Antibiotics were not prescribed randomly, but according to the clinical judgement of the clinicians and in response to the pathology of the patient. It is therefore possible that the effects of antibiotics in the model estimates may be confounded by the clinical severity of the infant. We did attempt to mitigate this by including a covariate for clinical severity, which links to the use of invasive devices and we note that, with the exception of imipenem, all antibiotics show similar effects to each other on the risk of first acquisition.
In summary, this study provides strong evidence for persontoperson within ward transmission of 3GCR K. pneumoniae s.l. mediated by healthcare workers, estimates key epidemiological parameters, quantifies the role of different antibiotics in driving the epidemic, and highlights interventions with the potential to contribute to control efforts.
Materials and methods
Epidemiological and microbiological data
Request a detailed protocolWe used data collected prospectively from a neonatal intensive care unit in a children’s hospital in Siem Reap, Cambodia between September 2013 (when the neonatal unit newly opened) and September 2014. The study protocol required rectal swabs to be taken within 48 hours of admission and subsequently every two to three days. We assumed that any patient testing positive on their first swab taken within 48 hours of admission was positive on arrival. Patients who had their first swab taken >48 hours from admission but tested negative were included in the analysis, while those that tested positive were omitted. Infants that were readmitted to the ward and had been colonised prior to first discharge were assumed to still be colonised and were thus excluded. Details of the microbiological treatment of rectal swabs, including resistance assays, have been published previously (Turner et al., 2016).
Probability of acquisition
View detailed protocolWe estimated the daily probability of colonisation with 3GCR K. pneumoniae s.l. using a discrete time model with time steps of one day. Each day in the ward a previously uncolonised patient can become colonised. As rectal swabs are not taken on every day of a patient’s stay the outcome is interval censored: we know that a negative swab followed by a positive swab indicates that a patient became colonised on some day between the two swabs, but not on which day. If the probability of becoming colonised on day i for patient j is p_{ij}, given the patient is uncolonised at the start of the day, then the probability of remaining uncolonised is (1p${}_{i\mathit{}j}$). In interval k for patient j consisting of N_{kj} days, then the probability of remaining uncolonised is:
Therefore the probability of becoming colonised (v_{kj}) is the complement:
The outcome for patient j in interval k, Y_{kj }∈ {0,1}, as the patient either becomes colonised (1) or remains uncolonised (0) with 3GCR Klebsiella. Therefore the likelihood is given by:
Risk factor models for carriage acquisition
Request a detailed protocolThe daily probability of becoming colonised (p_{ij}) is related by the logit link function to a linear function of covariates:
Where x_{1}, x_{2}, x_{3} … is a vector of predictors (data) and β_{1}, β_{2}, β_{3} … is a vector of slopes (parameters) that are to be estimated. The intercept α can be a single parameter, or permitted to vary over m periods of time in a random effects model. When using such a random effects model α was assumed to be normally distributed, with a mean µ and standard deviation σ, which are themselves parameters with their own prior distributions. The prior distributions used in the analysis are shown in Table 2. Results from models with alternative prior distributions are shown in Figure 2—figure supplement 2.
Fourteen standard covariates included in all risk factor models were: exposure to the six most common combinations of antibiotics, taken intravenously unless otherwise stated (ampicillin, ampicillin + gentamicin, cloxacillin (oral), ceftriaxone, cloxacillin + gentamicin, and imipenem) within the past 48 or 96 hours; whether the infant was breast milk fed; if the infant recieved an oral probiotic on entry (Lactobacillus acidophilus) for prevention of necrotising enterocolitis; sex; born prematurely (before the 37th week of pregnancy); severity (defined as either i. requiring ventilation, ii. requiring continuous positive airway pressure or iii. requiring inotopes); and if already colonised with 3GCR E. coli. These explanatory variables were treated as binary (0/1). We also included the age in days on first admission to the NU, and the daily number of nurses on the ward. An additional covariate included in one of the risk factor models (see below) was a term for colonisation pressure, which is an integer value representing the number of individuals known to be colonised with 3GCR K. pneumoniae s.l. on that day. Covariates were recorded for every day the infant was present in the neonatal unit and data were treated as complete. We considered the following models:
Single intercept with standard covariates (14). Exposure to antibiotics considered if taken within the past 96 hours.
Single intercept with standard covariates (14). Exposure to antibiotics considered if taken within the past 48 hours.
Single intercept with standard covariates (14) plus an additional covariate for colonisation pressure. Exposure to antibiotics considered if taken within the past 96 hours.
Variable intercept by study month that uses partial pooling (hierarchical model). Standard covariates (14) and exposure to antibiotics considered if taken within the past 96 hours.
Rectal swab/culture sensitivity
Request a detailed protocolSwab sensitivity was estimated from the number of negative swabs following a positive swab, i) under the assumption that all negative results following a positive swab were false negatives and ii) under the assumption that three or more consecutive negative swabs following a positive represented a true decolonisation event. Posterior distributions of the false negative rate and swab sensitivity were estimated using a conjugate beta prior; beta(α=1, β=7) (Bolker, 2008).
Pathogen sequencing and bioinformatics
Request a detailed protocolWe wholegenome sequenced 317 cultured isolates identified morphologically as 3GCR K. pneumoniae s.l. from i) rectal swabs from all colonised patients within a four month period of the study and ii) twice weekly swabs from seven environmental surfaces around the ward (6 sinks and one computer keyboard) within the same time frame. Sequencing was performed with the Illumina HiSeq 2500 platform, producing 150 basepair pairedend reads. The reads were trimmed for adapter sequence using TrimGalore (v0.4.4) before assembly with Unicycler (v0.4.5) (Wick et al., 2017), contigs <1 kilobase were discarded. Distances between assemblies were calculated using mash (v1.1) (Ondov et al., 2016) and a phylogeny constructed with mashtree (v0.33). Sequence types (STs) were identified using Kleborate 0.2.0 (Wyres et al., 2016). Variant calling was performed within STs by mapping reads to ST consensus reference genomes (5.32 Mbp) with SMALT (v.0.7.6) https://www.sanger.ac.uk/science/tools/smalt0 (parameters x y 0.85 r 1 j 100 i 800), before sorting and removing unpaired mate reads and technical duplicates from binary alignment files with samtools (v.1.8). Single nucleotide variants (SNPs) were called by piping output from samtools mpileup into bcftools (v.1.8) and were conservatively filtered to remove SNPs within 50 bp of indels, with a read depth <10x and >200x, a mapping quality score <30 or read quality score <100. Repetitive regions were identified with nucmer (v.3.1), phage regions with PHASTER (Arndt et al., 2016), and recombination with ClonalFrameML (v.1.11) (Didelot and Wilson, 2015), and these regions subsequently masked from SNP calling.
Transmission models incorporating sequence type data
Request a detailed protocolWe assessed withinward transmission of 3GCR K. pneumoniae s.l. STs under the assumption that individuals remain colonised with a given ST for the duration of their stay until discharge (Birgand et al., 2013). We fitted data to five linear transmission models.
Where the daily risk of acquiring any ST is constant (intercept only).
A constant term plus a covariate for colonisation pressure ($\beta $), where the explanatory variable (n${}_{i\mathit{}c}$) is the number of other patients colonised with ST c in the neonatal unit on day i.
As (2) with an additional term for the contribution of environmental contamination to transmission ($\gamma $), whereby the number of cases with ST c on an earlier day i’ (n${}_{{i}^{\mathrm{\prime}}\mathit{}c}$) are assumed to leave some trace in the environment that decays exponentially over time at a rate given by $\lambda $.
A hierarchical version of (2) which permits the transmission parameter $\beta $ to vary by ST.
A hierarchical version of (2) which permits the intercept $\alpha $ to vary by ST.
The probability p of colonisation for individual j on day i with ST c for the respective models are:
We opted to fit transmission models on the linear, rather than logistic, scale as we consider a linear increase in the force of infection with the number of colonised infants to be a more realistic assumption than the, initially, exponential increase that results from a logistic transformation.
Statistical model fitting
Request a detailed protocolWe fitted the statistical models using Hamiltonian Markov chain Monte Carlo in Stan (version 2.17.3) within the R environment (v. 3.4.3). Prior distributions were selected to be weakly informative normal distributions for the riskfactor models (McElreath, 2018), see Table 2. For the transmission models, beta distributions were used as priors. in the case of the hierarchical transmission models (4 and 5) the scale and location parameters are themselves hyperparameters with their own priors. Studies have shown the half life of carbapenemresistant K. pneumoniae inoculated onto hospital surfaces to be around 12 hours (Weber et al., 2015); this corresponds to λ=1.39 in our model. We therefore gave λ in transmission model 3 a more informative prior, normal(µ=1, σ=2), to give weight to biologically plausible estimates. Prior distributions for all transmission models are shown in Table 3. Results from a less informative prior distribution for λ are shown in Figure 3—figure supplement 2. Chains were run for a varying number of iterations depending on the number of parameters to estimate, though with a minimum of 12,000 iterations over four chains, including burnin. The GelmanRubin statistic ($\hat{R}$) was used as a diagnostic, where values <1.01 indicate chains have converged and additionally posterior chains were visually inspected for convergence. Effective sample sizes (ESS) for both the centre of the posterior distribution (bulk) and the ends of the distribution (tail) was ensured to be >400 (Vehtari et al., 2019). Model comparison was performed with widely applicable information criterion (WAIC) (Vehtari et al., 2017). We used 95% credible intervals (CrIs) as a measure of uncertainty around posterior parameter distributions and posterior medians as the central estimate. Uncertainty in parameter estimates is represented by the posterior distributions, and the estimated probability that a particular parameter is within a certain range is given by the area under the curve of that parameter’s marginal posterior distribution within that range.
Agentbased forward simulations
Request a detailed protocolWe forward simulated the impact of interventions using an agentbased model implemented in Python (v. 2.7.15) and hosted at https://github.com/tc13/wardinfectionABM (copy archived at https://github.com/elifesciencespublications/wardinfectionABM). In brief, we model in discrete time steps a ward containing a fixed number of beds and where patients sample a length of stay and colonisation status on entry. Simulations to estimate the ward reproduction number R_{A} introduced a single colonised individual into a full ward of susceptible patients. The probability of any uncolonised patient acquiring resistant Klebsiella from the index patient on day one is p_{ijc}, where p_{ijc} is sampled from a beta distribution with shape hyperparameters α and β from transmission model 4 (Table 3). This probability changes on subsequent days based on the number of other infants that become colonised and start transmission, and the lengths of stay. The simulation ends when the index patient is discharged.
For simulations that explored the impact of interventions, the effect of each of the covariates of interest (probiotic consumption, breast milk feeding and number of nurses on the ward) were obtained from risk factor model A (Table 2 and Figure 2A). The marginal effects of each risk factor were transformed into a logodds ratio and used to modify the daily risk of acquisition sampled from a beta distribution with shape hyperparameters α_{s} and β_{s}. For n_{c} colonised patients on the ward on day j, the probability of an uncolonised patient with covariate k acquiring 3GCR K. pneumoniae s.l. is 1(1 p_{k})^{ncj}. The simulation ends after 365 days. To reduce variability between model runs for estimation of R_{A} and the interventions, simulations were run 100 times with each posterior parameter sample and the mean outcome value obtained. As we used 2000 posterior samples from each parameter estimated by model fitting, this resulted in a total of 200,000 model runs for each simulated scenario.
Data availability
Code for reproducing the statistical model fitting and anonymised patient data are available at https://github.com/tc13/transmissiondynamicsklebsiella (copy archived at https://github.com/elifesciencespublications/transmissiondynamicsklebsiella). The code for the agent based model and parameter values for forward simulations are available at https://github.com/tc13/wardinfectionABM (copy archived athttps://github.com/elifesciencespublications/wardinfectionABM). Short read sequence data is available from NCBI under accession numbers PRJNA395864 and PR600JEB24970. Genomic data can be visualised in the Microreact project https://microreact.org/project/BV05TIXkU.

NCBI BioProjectID PRJEB24970. Whole genome sequencing of Klebsiella pneumoniae collected from patients and environment in a hospital in Cambodia.

NCBI BioProjectID PRJNA395864. Klebsiella pneumoniae on a Cambodian neonatal unit.
References

Patient density, nursetopatient ratio and nosocomial infection risk in a pediatric cardiac intensive care unitThe Pediatric Infectious Disease Journal 16:1045–1048.https://doi.org/10.1097/0000645419971100000008

PHASTER: a better, faster version of the PHAST phage search toolNucleic Acids Research 44:W16–W21.https://doi.org/10.1093/nar/gkw387

Duration of colonization by extendedspectrum βlactamaseproducing Enterobacteriaceae after hospital dischargeAmerican Journal of Infection Control 41:443–447.https://doi.org/10.1016/j.ajic.2012.05.015

Transmission routes of ExtendedSpectrum BetaLactamaseProducing Enterobacteriaceae in a neonatology ward in MadagascarThe American Journal of Tropical Medicine and Hygiene 100:1355–1362.https://doi.org/10.4269/ajtmh.180410

Infection prevention in a connected world: the case for a regional approachInternational Journal of Medical Microbiology 303:380–387.https://doi.org/10.1016/j.ijmm.2013.02.003

Singlemolecule sequencing to track plasmid diversity of hospitalassociated carbapenemaseproducing EnterobacteriaceaeScience Translational Medicine 6:254ra126.https://doi.org/10.1126/scitranslmed.3009845

An augmented data method for the analysis of nosocomial infection dataAmerican Journal of Epidemiology 168:548–557.https://doi.org/10.1093/aje/kwn176

A Large, Refractory Nosocomial Outbreak of Klebsiella pneumoniae CarbapenemaseProducing Escherichia coli Demonstrates Carbapenemase Gene Outbreaks Involving Sink Sites Require Novel Approaches to Infection ControlAntimicrobial Agents and Chemotherapy 62:e01689.https://doi.org/10.1128/AAC.0168918

Withinhost evolution of bacterial pathogensNature Reviews Microbiology 14:150–162.https://doi.org/10.1038/nrmicro.2015.13

ClonalFrameML: efficient inference of recombination in whole bacterial genomesPLOS Computational Biology 11:e1004041.https://doi.org/10.1371/journal.pcbi.1004041

Use of stochastic epidemic modeling to quantify transmission rates of colonization with methicillinresistant Staphylococcus aureus in an intensive care unitInfection Control & Hospital Epidemiology 26:598–606.https://doi.org/10.1086/502588

Antimicrobial resistance in invasive bacterial infections in hospitalized children, Cambodia, 20072016Emerging Infectious Diseases 24:841–851.https://doi.org/10.3201/eid2405.171830

Bacterial colonization and infection in the neonateThe American Journal of Medicine 70:417–422.https://doi.org/10.1016/00029343(81)907828

Gastrointestinal carriage is a major reservoir of Klebsiella pneumoniae infection in intensive care patientsClinical Infectious Diseases 65:208–215.https://doi.org/10.1093/cid/cix270

Mathematical modelling: a tool for hospital infection controlThe Lancet Infectious Diseases 6:39–45.https://doi.org/10.1016/S14733099(05)70325X

The transmissibility of AntibioticResistant Enterobacteriaceae in intensive care unitsClinical Infectious Diseases 66:489–493.https://doi.org/10.1093/cid/cix825

Nursing resources: a major determinant of nosocomial infection?Current Opinion in Infectious Diseases 17:329–333.https://doi.org/10.1097/01.qco.0000136931.83167.d2

Epidemiologic studies on gramnegative bacilli in the hospital and communityAmerican Journal of Epidemiology 85:45–60.https://doi.org/10.1093/oxfordjournals.aje.a120674

Control of gramnegative bacilli in a hospital nurseryAmerican Journal of Diseases of Children 107:363–369.https://doi.org/10.1001/archpedi.1964.02080060365007

Fast and flexible bacterial genomic epidemiology with PopPUNKGenome Research 29:304–316.https://doi.org/10.1101/gr.241455.118

Antimicrobial use and antimicrobial resistance: a population perspectiveEmerging Infectious Diseases 8:347–354.https://doi.org/10.3201/eid0804.010312

BookStatistical Rethinking: A Bayesian Course with Examples in R and StanCRC Press.

A hospital outbreak of Klebsiella pneumonia from inhalation therapy with contaminated aerosol solutionsThe American Review of Respiratory Disease 95:454–460.https://doi.org/10.1164/arrd.1967.95.3.454

Epidemiology of Klebsiella and hospitalassociated infectionsClinical Infectious Diseases 1:736–753.https://doi.org/10.1093/clinids/1.5.736

Global spread of Carbapenemaseproducing EnterobacteriaceaeEmerging Infectious Diseases 17:1791–1798.https://doi.org/10.3201/eid1710.110655

Recovery of gut Microbiota of healthy adults following antibiotic exposureNature Microbiology 3:1255–1265.https://doi.org/10.1038/s4156401802579

Hospitalacquired infections due to gramnegative BacteriaNew England Journal of Medicine 362:1804–1813.https://doi.org/10.1056/NEJMra0904124

Evidencebased model for hand transmission during patient care and the role of improved practicesThe Lancet Infectious Diseases 6:641–652.https://doi.org/10.1016/S14733099(06)706004

Antimicrobialresistant Gramnegative colonization in infants from a neonatal intensive care unit in ThailandJournal of Hospital Infection 103:151–155.https://doi.org/10.1016/j.jhin.2019.04.004

Nurse staffing and NICU infection ratesJAMA Pediatrics 167:444–450.https://doi.org/10.1001/jamapediatrics.2013.18

Nosocomial Klebsiella infections: intestinal colonization as a reservoirAnnals of Internal Medicine 74:657.https://doi.org/10.7326/00034819745657

Tracking a hospital outbreak of carbapenemresistant Klebsiella pneumoniae with wholegenome sequencingScience Translational Medicine 4:148ra116.https://doi.org/10.1126/scitranslmed.3004129

Colonization with Enterobacteriaceae producing ESBLs in children attending preschool childcare facilities in the lao people's Democratic RepublicJournal of Antimicrobial Chemotherapy 70:1893–1897.https://doi.org/10.1093/jac/dkv021

Global antimicrobial resistance in Gramnegative pathogens and clinical needCurrent Opinion in Microbiology 39:106–112.https://doi.org/10.1016/j.mib.2017.10.028

High prevalence of Antimicrobialresistant Gramnegative colonization in hospitalized cambodian infantsThe Pediatric Infectious Disease Journal 35:856–861.https://doi.org/10.1097/INF.0000000000001187

Practical bayesian model evaluation using leaveoneout crossvalidation and WAICStatistics and Computing 27:1413–1432.https://doi.org/10.1007/s1122201696964

Carbapenemresistant Enterobacteriaceae: frequency of hospital room contamination and survival on various inoculated surfacesInfection Control & Hospital Epidemiology 36:590–593.https://doi.org/10.1017/ice.2015.17

Nosocomial infection with Klebsiella type 25New England Journal of Medicine 275:17–22.https://doi.org/10.1056/NEJM196607072750104

Unicycler: resolving bacterial genome assemblies from short and long sequencing readsPLOS Computational Biology 13:e1005595.https://doi.org/10.1371/journal.pcbi.1005595

PHYLOSCANNER: inferring transmission from within and BetweenHost pathogen genetic diversityMolecular Biology and Evolution 35:719–733.https://doi.org/10.1093/molbev/msx304

Identification of Klebsiella capsule synthesis loci from whole genome dataMicrobial Genomics 2:e000102.https://doi.org/10.1099/mgen.0.000102

Klebsiella pneumoniae population genomics and AntimicrobialResistant clonesTrends in Microbiology 24:944–956.https://doi.org/10.1016/j.tim.2016.09.007

Clinical significance of rising incidence of infections due to GramNegative bacilliPostgraduate Medicine 17:413–419.https://doi.org/10.1080/00325481.1955.11708216

A current perspective on antimicrobial resistance in southeast asiaJournal of Antimicrobial Chemotherapy 72:2963–2972.https://doi.org/10.1093/jac/dkx260
Decision letter

Miles P DavenportReviewing Editor; University of New South Wales, Australia

Neil M FergusonSenior Editor; Imperial College London, United Kingdom

Miles P DavenportReviewer; University of New South Wales, Australia

Jon ZelnerReviewer; UMICH

Edward J FeilReviewer; Univesity of Bath, United Kingdom
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This manuscript addresses the important question of transmission of drug resistance in a hospital setting. It analyses a dataset of children regularly swabbed for Klebsiella pneumoniae after admission to an intensive care unit in Cambodia. Such data is inherently complex and difficult to interpret. The authors use mathematical modelling approach to investigate the role of a number of factors such as colonization pressure, antibiotic treatment, and breastfeeding. The results provide insights into the spread of drug resistant bacteria in an important setting for the development and transmission of drug resistant bacteria.
Decision letter after peer review:
Thank you for submitting your article "Transmission dynamics and control of multidrugresistant Klebsiella pneumoniae in neonates in a developing country" for consideration by eLife. Your article has been reviewed by four peer reviewers, including Miles P Davenport as the Reviewing Editor, and the evaluation has been overseen by a Neil Ferguson as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Jon Zelner; Edward J Feil.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
This manuscript uses a modelling approach to investigate acquisition of drug resistant bacteria in a hospital setting. The authors study various risk factors for 'acquisition', and find that antibiotic use was the only significant effect. The strength of this study is that it looks at a very clinically important organism in the exact setting in which it is problematic. The value of the study is tempered by the fact that much of the screening data was already published, along with more coarsegrained analysis of the role of antibiotic on acquisition. The modelling and analysis is therefore a central factor in the study, but also falls short in several ways that need to be addressed. The reviewers expressed that there are several possibilities to improve the analysis that need to be explored.
Essential revisions:
1) the Bayesian approach to model fitting of the risk factor and the transmission models, uses arbitrary priors, that do not truly reflect prior probability. Despite the authors claim that they were selected to be weakly informative, most prior distributions cover a range that make it possible that they are significantly influencing the posterior. A straightforward way of conveying this to the reader is to plot the prior and the posteriors together. This is a key shortcoming because if the model estimates are driven by these arbitrary choices, much of the remaining paper comes into question. This is certainly a criticism that can be addressed by changing the fitting approach to remove the influence of the priors, or by switching to another approach, maximum likelihood or proportional hazard models.
2) the results of the model are not presented carefully, many of the risk factors have credibility intervals that would not generally be considered significant, yet are presented as if they are (e.g.: "covariates associated with reduced daily risk of acquisition," but nothing in this paragraph is significant; "suggests that contamination left by previously colonized infants decays rapidly to background levels", when there is no support for the model with this parameter and the λ estimate is heavily influenced by the prior with mean 1 and variance 2). This is also true for the results presented in Figure 4B, C, D and the lengthy discussion in Discussion section, which fails to accurately convey the possibility that these factors play no role in colonization is also consistent with the data. The large uncertainty in most parameter estimates likely reflects the challenge of fitting 109 observations with a 14parameter model. This needs to be stated and the limitations emphasized.
3) the analysis of the transmission models appears to have a bias toward finding model 4 to be the best because there is a single value of α for all the sequence types, but there are different average colonization probabilities across the sequence types (because some had many more observed colonization events than others). This variability it likely to be picked up, at least to some extent, in the ST specific betas. A fairer test is whether ST specific β fits better than a ST specific alphas. The temporal clustering of ST in Figure 3F suggest this may be true, but the approach seem to have been biased to find an effect of β. Are there explanations for the temporal clustering of the ST other than transmission in the NU, such a transient increased in the community or in the water supply etc.
4) A concern with this analysis relate to the way the different infection models are described. Specifically, the differentiation between a 'risk factor model' and a 'transmission model'. It is not clear from the main text that the transmission models considered are essentially hierarchical regression models in which the impact of an increasing number of colonized patients is considered to have an additive rather than a multiplicative effect on individual risk. This is not a critique of the modeling, which is generally appropriate and welldone. But it is difficult to read the manuscript and understand the distinction between these two sets of analyses. In addition, the authors should clarify earlier that the nonsignificant, negative OR for colonization pressure may be due to this fact as a way of framing the need for an additive model that more able to directly account for the impact of transmission on daily infection risk.
5) It is not clear to me why the transmission models could not incorporate some of the covariates used in the first set of risk factor models. For example, age at admission or the number of nurses on the unit could be used to modulate infection risk, i.e. be used to calculate an individual risk ratio that is then multiplied by the hazard of infection, e.g. exp(\ζ X_{i})*(\α +\β n_{i}), where X_{i} is a set of individual or dayspecific risk factors, and \ζ is a vector of logrisk ratios. This would seem to make the most of the data and incorporate both the additive effect of colonization pressure with the riskmodulation of individual level factors.
https://doi.org/10.7554/eLife.50468.sa1Author response
Summary:
This manuscript uses a modelling approach to investigate acquisition of drug resistant bacteria in a hospital setting. The authors study various risk factors for 'acquisition', and find that antibiotic use was the only significant effect. The strength of this study is that it looks at a very clinically important organism in the exact setting in which it is problematic. The value of the study is tempered by the fact that much of the screening data was already published, along with more coarsegrained analysis of the role of antibiotic on acquisition. The modelling and analysis is therefore a central factor in the study, but also falls short in several ways that need to be addressed. The reviewers expressed that there are several possibilities to improve the analysis that need to be explored.
We thank the four reviewers, including the reviewing and senior editors for their comments and criticisms of our article. We have taken these points into consideration and conducted a number of additional analyses to confirm the robustness of our results to different modelling decisions (e.g. prior distributions). For other criticisms we consider it more appropriate to include these as study limitations in the Discussion.
While the individual points are dealt with below, we would like to respectfully challenge the reviewers’ reference to statistical significance, e.g. in the Summary and Essential revisions 2 and 4. The limitations and misuse of significance testing has been widely discussed by the statistical community, please see “American Statistical Association [ASA] statement on Statistical Significance and pvalues”^{1}, and discussions by Gelman^{2,3}, McElreath^{4,5} and Spiegelhalter^{6}. In line with current statistical thinking which recommends that statistical “bright lines” (such as p<0.05) are avoided, we make no reference to “significant” or “nonsignificant” results in our article. Rather we consider the full posterior distribution to be informative about parameter estimates, in addition to interpreting our results in the context of the existing literature. In this analysis we are not setting up null models to be refuted, but rather estimating biomedically relevant parameters based on our data and the model structure. We do not consider any of our models to be correct, but instead provide a number of possible models and use an information criterion (and our background biological knowledge) as a guide for the most appropriate model. This is in line with best practices in ecological modelling that have been established for some time^{7,8}. We hope therefore that the reviewers will consider our analyses and the response to their comments in this light.
Essential revisions:
1) the Bayesian approach to model fitting of the risk factor and the transmission models, uses arbitrary priors, that do not truly reflect prior probability. Despite the authors claim that they were selected to be weakly informative, most prior distributions cover a range that make it possible that they are significantly influencing the posterior. A straightforward way of conveying this to the reader is to plot the prior and the posteriors together. This is a key shortcoming because if the model estimates are driven by these arbitrary choices, much of the remaining paper comes into question. This is certainly a criticism that can be addressed by changing the fitting approach to remove the influence of the priors, or by switching to another approach, maximum likelihood or proportional hazard models.
We agree that plotting the prior and posterior distributions together is a useful method to visualise the effect of the prior distribution on the result, and we have done this as suggested (see Figure 2—figure supplement 2 and Figure 3—figure supplement 2). Our priors were chosen on the basis of best practice guidelines provided by leading statisticians, including the creators of the STAN fitting software^{9}. The most important validation of our model is the posterior predictive distribution (Figure 2D), which shows that the Risk Factor Model closely predicts the observed number of colonisation events.
While the choice of prior distributions remains an area of active debate, it is important to note that the prior will always influence the posterior. Prior distributions are an important tool for i) incorporating prior beliefs or information into the model, ii) preventing the model from exploring regions of parameter space that are biologically implausible, and iii) avoiding overfitting. Indeed, the use of weakly informative priors is directly analogous to the use of standard regularisation methods in frequentist analysis. Methods from frequentist analysis mentioned by the reviewers (maximum likelihood, proportional hazards models) are not somehow “more objective” and contain their own, although often less explicitly stated, assumptions. Using diffuse or uniform (“flat”) priors with a large range in Bayesian inference will typically give similar results to maximum likelihood. However, the use of “flat” or uniform priors is not generally recommended^{10}.
We do, however, fully agree that there is value in exploring sensitivity of results to the choice of prior distributions. Following this criticism of the reviewers, we have therefore explored different prior distributions for Risk Factor Model A to test the sensitivity of the result to the prior. Varying the priors for Risk Factor Model A to normal distributions with much lower variance, as suggested, results in only small changes to the posterior distributions (in Figure 2—figure supplement 2); WAIC is also unaffected.
2) the results of the model are not presented carefully, many of the risk factors have credibility intervals that would not generally be considered significant, yet are presented as if they are (e.g.: "covariates associated with reduced daily risk of acquisition," but nothing in this paragraph is significant; "suggests that contamination left by previously colonized infants decays rapidly to background levels", when there is no support for the model with this parameter and the λ estimate is heavily influenced by the prior with mean 1 and variance 2). This is also true for the results presented in Figure 4B, C, D and the lengthy discussion in Discussion section, which fails to accurately convey the possibility that these factors play no role in colonization is also consistent with the data. The large uncertainty in most parameter estimates likely reflects the challenge of fitting 109 observations with a 14parameter model. This needs to be stated and the limitations emphasized.
In line with current statistical thinking we think it is inappropriate to dichotomise results into significant/ nonsignificant (see comment above and statement from American Statistical Association^{1}), and so consider our presentation and discussion of the posterior odds ratio distributions appropriate. We disagree that there is “no support” for transmission model 3; there is less support according to WAIC, but due to the interest in colonisation with resistant K. pneumoniae resulting from contact with contaminated surfaces in hospitals^{11}, we consider it is appropriate to discuss the estimates from this model. The results from our analysis should not be interpreted from a frequentist perspective; posterior distributions inherently account for uncertainty in the parameter estimates. For instance if 6% of the posterior odds distribution for a covariate is <1, then there is a 94% probability that the covariate increases the risk of acquisition / detection. We have made this clearer in the text; in the Materials and methods subsection “Statistical Model Fitting” we add “Uncertainty in parameter estimates is represented by the posterior distributions, and the estimated probability that a particular parameter is within a certain range is given by the area under the curve of that parameter’s marginal posterior distribution within that range.”
To address the reviewers’ concerns that the estimate of λ in transmission model 3 is heavily influenced by the prior distribution, we have refit the model using a less informative prior (Figure 3—figure supplement 2. The model estimate of the K. pneumoniae halflife on surfaces changes from 3.6 hours with our initial prior distribution (normal(1,2)) to 2.5 hours with a less informative prior (normal(0,5)). Our interpretation of the results from these models is qualitatively the same (“K. pneumoniae surface contamination decays quickly to background levels”), and altering the prior distribution had no effect on WAIC. The only published value from experimental studies for resistant K. pneumoniae halflife on surfaces gave a value of ~12 hours^{12}, which is equivalent to a parameter estimate of 1.39 in our model. Our initial prior distribution (normal(1,2)) therefore gives more weight to biologically plausible areas of parameter space, and we continue to use the result from this model in the main text. We have expanded the text to better explain our prior choice for this model– see Materials and methods subsection “Transmission Models Incorporating Sequence Type Data”.
3) the analysis of the transmission models appears to have a bias toward finding model 4 to be the best because there is a single value of α for all the sequence types, but there are different average colonization probabilities across the sequence types (because some had many more observed colonization events than others). This variability it likely to be picked up, at least to some extent, in the ST specific betas. A fairer test is whether ST specific β fits better than a ST specific alphas. The temporal clustering of ST in Figure 3F suggest this may be true, but the approach seem to have been biased to find an effect of β. Are there explanations for the temporal clustering of the ST other than transmission in the NU, such a transient increased in the community or in the water supply etc.
We agree that it would be interesting to explore variance in the α parameter and if this results in a better fit to data by WAIC than variance in β. We have therefore added a model which allows α to vary by sequence type (transmission model 5, described in Materials and methods). The WAIC from this model is much higher than from transmission models 24 (1793, compared with 1739, 1740 and 1733), indicating a poor fit to data (shown in Table 3).
We agree that other explanations for the temporal clustering of STs are possible, though less parsimonious. Overall the estimate of α (background contamination in the absence of colonised infants in the ward) is very low in all models, suggesting that effects other than persontoperson transmission are minimal in spreading resistant K. pneumoniae around the ward. Transmission from the community may be driving imported cases into the ward, and the variability of the K. pneumoniae population that we observe in the sequence data (Figure 3A) may reflect a large, diverse population of community K. pneumoniae that has not yet been investigated. We have added a sentence to the Discussion that “We cannot entirely rule out less parsimonious explanations for the temporal clustering of STs, such as a transient increase of certain STs in the water supply, though the persontoperson transmission route is most strongly supported by our models.”
4) A concern with this analysis relate to the way the different infection models are described. Specifically, the differentiation between a 'risk factor model' and a 'transmission model'. It is not clear from the main text that the transmission models considered are essentially hierarchical regression models in which the impact of an increasing number of colonized patients is considered to have an additive rather than a multiplicative effect on individual risk. This is not a critique of the modeling, which is generally appropriate and welldone. But it is difficult to read the manuscript and understand the distinction between these two sets of analyses. In addition, the authors should clarify earlier that the nonsignificant, negative OR for colonization pressure may be due to this fact as a way of framing the need for an additive model that more able to directly account for the impact of transmission on daily infection risk.
We agree that there is a need to clarify the analysis earlier in the manuscript and frame the reasoning for our additive transmission models. We have amended the text at the end of the Introduction to better introduce the different types of statistical models that follow in the analysis; “We fit four models with logit link functions to estimate the impact of covariates on the daily risk of acquisition or detection of 3GCR K. pneumoniae s.l.. As these models are unable to identify the force of infection, and genomic data show the ward 3GCR K. pneumoniae s.l. to be a highly heterogeneous bacterial community, we then fit five linear transmission models to estimate the force of infection for acquisition or detection of each ST.”
5) It is not clear to me why the transmission models could not incorporate some of the covariates used in the first set of risk factor models. For example, age at admission or the number of nurses on the unit could be used to modulate infection risk, i.e. be used to calculate an individual risk ratio that is then multiplied by the hazard of infection, e.g. exp(ζ X_{i})*(α +β n_{i}), where X_{i} is a set of individual or dayspecific risk factors, and ζ is a vector of logrisk ratios. This would seem to make the most of the data and incorporate both the additive effect of colonization pressure with the riskmodulation of individual level factors.
We did not include additional covariates when fitting the transmission models, as the K. pneumoniae sequence data was only taken within a four month window of the yearlong study. In this period there were 171 first acquisition/detection events for different STs, however given the number of ST atrisk patient days was >48,000, the events were much sparser than the 109 acquisition events over 864 patient days in the risk factor model and are likely to have wider uncertainty intervals. As the underlying data are different in the risk factor and transmission models and as the models use different linkfunctions, the estimates of covariates will not be exactly the same and we considered that presenting two sets of odds ratios for the risk factors might confuse the reader.
We performed analysis similar to that described by the reviewer above to obtain the parameters for the individualbased model; we varied the force of infection inferred from the transmission models with the effect of covariates inferred from the risk factor models, by converting all parameters onto the logodds scale to create a combined linear function and then using the inverselogit function to obtain a probability of acquisition. This is described in the Materials and methods as “For simulations that explored the impact of interventions, the effect of each of the covariates of interest (probiotic consumption, breast milk feeding and number of nurses on the ward) were obtained from risk factor model A (Table 2 and Figure 2A). The marginal effects of each risk factor were transformed into a logodds ratio and used to modify the daily risk of acquisition sampled from a β distribution with shape hyperparameters α_{s} and β_{s}.”
References
1) Ronald L. Wasserstein and Nicole A. Lazar (2016) The ASA Statement on pValues: Context, Process, and Purpose, The American Statistician, 70:2, 129133, DOI: 10.1080/00031305.2016.1154108
2) Blakeley B. McShane, David Gal, Andrew Gelman, Christian Robert and Jennifer
L. Tackett (2019) Abandon Statistical Significance, The American Statistician, 73:sup1, 235245, DOI: 10.1080/00031305.2018.1527253
3) Gelman, A. and Stern, H., 2006. The difference between “significant” and “not significant” is not itself statistically significant. The American Statistician, 60(4), pp.328331.
4) McElreath, R., 2018. Statistical rethinking: A Bayesian course with examples in R and Stan. Chapman and Hall/CRC.
5) Smaldino, P.E. and McElreath, R., 2016. The natural selection of bad science. Royal Society open science, 3(9), p.160384.
6) Matthews, R., Wasserstein, R. and Spiegelhalter, D., 2017. The ASA’s pvalue statement, one year on. Significance, 14(2), pp.3841
7) Hilborn, R. and Mangel, M., 1997. The ecological detective: confronting models with data. Princeton University Press
8) Bolker, B.M., 2008. Ecological models and data in R. Princeton University Press
9) https://github.com/standev/stan/wiki/PriorChoiceRecommendations
10) Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. and Rubin, D.B., 2013. Bayesian data analysis. Chapman and Hall/CRC.
11) Freeman, J.T., Nimmo, J., Gregory, E., Tiong, A., De Almeida, M., McAuliffe, G.N. and Roberts, S.A., 2014. Predictors of hospital surface contamination with Extendedspectrum βlactamaseproducing Escherichia coli and Klebsiella pneumoniae: patient and organism factors. Antimicrobial resistance and infection control, 3(1), p.5
12) Weber, D.J., Rutala, W.A., Kanamori, H., Gergen, M.F. and SickbertBennett, E.E., 2015. Carbapenemresistant Enterobacteriaceae: frequency of hospital room contamination and survival on various inoculated surfaces. infection control and hospital epidemiology, 36(5), pp.590593
https://doi.org/10.7554/eLife.50468.sa2Article and author information
Author details
Funding
Wellcome (106698/Z/14/Z)
 Nicholas PJ Day
Medical Research Council (MR/K006924/1)
 Ben S Cooper
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank the staff and patients’ families at the Angkor Hospital for Children for their participation in the study. We also thank Anastasia Hernandez Koutoucheva and Patrick Musicha for their comments on the manuscript.
Ethics
Human subjects: Written consent was obtained from mothers before study enrolment. The study was reviewed and approved by the Angkor Hospital for Children Institutional Review Board (1055/13 AHC) and the University of Oxford Tropical Ethics Committee (104713).
Senior Editor
 Neil M Ferguson, Imperial College London, United Kingdom
Reviewing Editor
 Miles P Davenport, University of New South Wales, Australia
Reviewers
 Miles P Davenport, University of New South Wales, Australia
 Jon Zelner, UMICH
 Edward J Feil, Univesity of Bath, United Kingdom
Publication history
 Received: July 25, 2019
 Accepted: November 26, 2019
 Accepted Manuscript published: December 3, 2019 (version 1)
 Version of Record published: January 23, 2020 (version 2)
Copyright
© 2019, Crellen et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,615
 Page views

 209
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.