Disentangling the force of infection of SARS-CoV-2 in Dutch long-term care facilities

TL;DR

This study analyzed how SARS-CoV-2 spread in Dutch long-term care facilities, tracking infection rates and sources over time.

Contribution

The study introduces a Bayesian model to estimate infection sources and transmission rates in LTCFs during the pandemic.

Findings

The force of infection was highest in December 2020 and lowest in June 2021.

Vaccination reduced transmission rates, but the general population remained a key source of infections.

Resident susceptibility increased in late 2021, partly due to the Delta variant.

Abstract

During the COVID-19 pandemic, residents of long-term care facilities (LTCFs) were disproportionately affected. To inform decision-making around interventions, we quantified the SARS-CoV-2 infection risk for residents and the relative contribution of different infection sources. We estimated the force of infection (FOI) experienced by Dutch LTCF residents over time and quantified the contributions of residents, LTCF healthcare workers (HCWs), and the general population.

Case data were obtained by Municipal Health Services as part of the Dutch national surveillance program. During the study period (1 October 2020 to 10 November 2021), testing policies included symptom-based testing, exposure-based testing, and facility-wide serial testing. We used a data augmentation approach to include uncertainty in the timing of infection, while taking account of different testing policies. We constructed a Bayesian generalized linear model to estimate group-specific transmission rate parameters and contributions to the FOI experienced by residents.

During the study period 36,877 cases were registered among residents and 19,676 among HCWs. The total daily FOI towards residents was highest in December 2020 (1.7*10−3, 95% CI: 1.5*10−3 – 1.9*10−3) and lowest in June 2021 (1.1*10−5 95%CI: 7.6*10−6 – 1.7*10−5). Transmission rate parameters and FOI declined as COVID-19 vaccination rollout started in residents, HCWs, and the older general population (February-May 2021). Most resident infections in spring and summer 2021 (April to August) were attributable to infections in the general population. The relative contribution of the general population to the FOI decreased in July 2021 when vaccination was available population-wide. In October-November 2021, transmission rate parameters and FOI increased again. We observed an increase in residents’ susceptibility to infection in this period, which was only partially explained by the emergence of the Delta variant.

Relative contributions to the FOI in LTCF residents varied substantially over time. COVID-19 vaccination appears to have been effective in reducing SARS-CoV-2 transmission towards residents, although other factors such as seasonality or non-pharmaceutical interventions may also have contributed to this. Such estimates of temporal trends in contributions to the FOI in LTCF residents can help inform future intervention efforts.

The online version contains supplementary material available at 10.1186/s12889-025-25015-x.

Figures

Descriptive results. A Epidemic curve based on reporting date in residents and HCW. Vertical lines and letters indicate study periods. B Correlation between daily number of reported resident and HCW cases. The fitted linear regression model is shown with its corresponding R^2^ value. C Epidemic curve based on reporting date in the general population. Vertical lines and letters indicate study periods. D Correlation between daily number of reported resident and general population cases. The fitted linear regression model is shown with its corresponding R^2^ value. E Cumulative density plot of ou

Model fits and outcomes for the best-fitting model. A Comparison between the daily number of observed infections (in the data) and predicted infections (from the model). Values of the 7-day average are shown. Dark lines indicate mean values. Shaded areas show 95% prediction intervals obtained by generating 1000 predictions by randomly sampling from the posterior distributions and augmented datasets. B Relative susceptibility to infection in residents. The horizontal line at y=1 indicates no difference compared to period A. C Daily relative force of infection from each LTCF group towards reside

Full text

Introduction

Long-term care facilities (LTCFs), which house amongst the most vulnerable people in our societies, were disproportionately affected during the Coronavirus disease (COVID-19) pandemic. LTCF residents faced an increased risk of high COVID-19 impact, due to a combination of risk factors for respiratory infections, such as frailty, comorbidities, crowded living conditions, and high-frequency contact with health care workers (HCWs) [1, 2]. The lack of appropriate protective equipment for HCWs in the early stages of the COVID-19 pandemic further contributed to this. The impact is reflected in high incidence and disease burden, peaking in the high age groups [3]. Seroprevalence estimates among healthcare workers (HCWs) in Dutch LTCFs also indicated high Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2) infection rates relative to the total HCW population (16.9% in LTCF HCWs vs. 7.5% in total HCW population, after the first wave) [4]. In the first half of 2020, deaths among LTCF residents accounted for 37–66% of all COVID-19 related deaths across several EU/EEA countries [5]. Estimations of the infection fatality rate for LTCF residents in this period are high, with estimates around 22% [6].

In response to the many outbreaks in LTCFs, the WHO published a policy brief with recommendations for managing SARS-CoV-2 transmission in LTCFs [7]. Recommendations included the implementation of and adherence to infection prevention and control standards in LTCFs as well as prioritization of testing, contact tracing, and monitoring of spread in both residents and health care workers [7]. In line with these recommendations, governments introduced non-pharmaceutical interventions (NPIs) such as visitor bans, individual movement restrictions, quarantine policies, masking policies, and a plethora of testing strategies [8]. However, NPIs might also affect physical and mental well-being of residents and their families and therefore require a careful balance and judicious implementation [9].

NPIs can be categorized into two groups: those focused on reducing the introduction of infections into the LTCF and those focused on reducing onwards transmission within the LTCF after introduction. The potential impact of a specific NPI depends on the contribution of introductions relative to within-LTCF transmission to the overall case load. Here, we consider within-LTCF transmission as the subsequent infections occurring after an introduction, between and among HCWs and residents. LTCFs are generally open for people to enter and leave, albeit sometimes with restrictions such as visiting hours or locked units where residents cannot exit freely for their own safety. SARS-CoV-2 outbreaks in LTCFs were found to be correlated to transmission patterns of SARS-CoV-2 in the general population [10–13]. Also, the number of HCWs at the facility, after controlling for LTCF size, was identified as a risk factor for infections among residents [14], and staff and residents were shown to be frequently infected with genetically similar SARS-CoV-2 viruses [15], indicating transmission between HCWs and residents. Similarly, the positive association between multibedded rooms as well as high occupancy rates with infection risk indicates a role for resident-to-resident transmission [16]. Such positive association has also been observed for other respiratory infections [17]. Despite evidence of transmission originating from these three groups (i.e., general population, HCWs, and residents), a quantitative understanding of their contribution to transmission and the role of interventions therein is lacking.

In this study, we aimed to disentangle the contributions of residents, HCWs, and the general population to the force of infection (FOI) experienced by LTCF residents (i.e., daily rate at which susceptible LTCF residents get infected). First, we described the characteristics of the COVID-19 epidemic in LTCFs in the Netherlands. We then used a Bayesian modelling framework including data augmentation to estimate the transmission rate parameters of each group towards residents in the period between October 2020 and November 2021. We also explored the changes in these rates over time and discuss how these relate to changes in the epidemiological context.

Data and methods

Setting

A total of 3.2 million SARS-CoV-2 positive tests were recorded in the Netherlands in the years 2020 and 2021 [18], with an estimated infection seroprevalence of 26% in November 2021 [19]. The first detection of SARS-CoV-2 infection in the Netherlands was on 27 February 2020 [20]. COVID-19 vaccine rollout started in HCW and LTCF residents ahead of the general population. We distinguished five periods based on progress in vaccination campaigns. The chosen dates do not represent exact starting dates of vaccination, but approximate the moment when immunity started to develop in the relevant population because of vaccination. Period A (1 October 2020–25 January 2021): pre-vaccination period. Period B (26 January 2021–28 February 2021) and period C (1 March 2021–30 May 2021): the periods after the first (B) and second (C) vaccination round for the primary series in LTCF residents and HCW. Period D (31 May 2021–30 September 2021): general population was getting vaccinated. Period E (1 October 2021–10 November 2021): period after vaccination of the general population. Data from before October 2020 was excluded due to limited accessibility to testing.

Dutch LTCFs for older people with high-level care needs can range from small-scaled homelike nursing homes to large residential care centers that have both nursing home and care home/service flat facilities. Typically, LTCFs have common areas for meals and socializing, and the rooms are predominantly one-bedded. LTCFs offer specialized wards, mostly either somatic wards for physically disabled residents, psychogeriatric wards for residents with dementia, or geriatric rehabilitation wards. In addition to society-wide control measures, such as physical distancing, hygiene measures, and self-quarantine, specific measures were implemented in LTCFs [21]. These included isolation and quarantine measures and visitor restrictions. Testing policies during the study period (1 October 2020-10 November 2021) included symptom-based testing, exposure-based testing, and facility-wide serial testing, although adherence to these measures varied [22].

Data collection

Data were obtained as part of the national surveillance program, through contact tracing by the Municipal Health Service (GGD: Gemeentelijke Gezondheidsdienst). Data from 24 out of 25 GGD regions were included. We used two types of data: (1) cases in the general population were obtained from aggregated data with total number of cases (i.e., positive RT-PCR test) per GGD region per day and (2) cases in HCW and residents were obtained from an individual-level dataset with data about cases associated with LTCFs for elderly care. A case was classified as associated with an LTCF when this person appeared in an LTCF-associated cluster (i.e., most likely setting of infection based on contact-tracing information). The individual-level dataset included: sex, age, GGD region (using an anonymized ID), symptom onset date, date of positive test result, date of notification, all-cause mortality (yes/no and date), whether they lived or worked in an LTCF, and a cluster identifier (i.e., pseudonymized LTCF identifier if the LTCF was the most likely setting of infection). A large LTCF organization with multiple facility locations was viewed as consisting of multiple LTCFs, such that an LTCF refers to a specific LTCF location in this study.

Additional data sources included the age distribution of the general population (Statistics Netherlands (CBS) [23]) and the bed capacity of LTCFs obtained through Patiëntenfederatie Nederland [24]. This information was available for 26% of LTCF locations and was linked to LTCFs by matching postal codes before further pseudonymization of the study dataset.

Resident and HCW cases associated with LTCFs (N = 56,553) were included in descriptive analyses. Individuals who were recorded both as resident and HCW (n = 137) were classified as resident. In these analyses, we used a Kruskal-Wallis rank sum test to compare outbreak sizes between study periods and a Pearson’s chi-squared test to compare the number of LTCFs between regions.

Statistical model

We quantified the force of infection experienced by LTCF residents. The FOI ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\lambda\:$$\end{document} ) denotes the rate at which susceptible individuals become infected per day. Assuming infections follow a Poisson process, the probability of a susceptible individual to acquire infection during a period \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\varDelta\:t$$\end{document} (1 day) is related to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\lambda\:$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:P=1-{e}^{-\lambda\:\varDelta\:t}$$\end{document}

Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\lambda\:\:$$\end{document} depends on the fraction of individuals who are infectious ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\frac{I}{N}$$\end{document} , where I is the number of infectious individuals and N is the population size) and on the transmission rate parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\beta\:$$\end{document} : the average number of daily new infections per infectious individual, such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\lambda\:=\beta\:\frac{I}{N}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\beta\:$$\end{document} can be considered the product of the contact rate and the probability of a contact between an infectious and susceptible individual to result in infection. The latter is related to infectivity and susceptibility of the infectee and recipient, respectively.

We distinguished the contribution to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\lambda\:$$\end{document} by three groups of infectious individuals in LTCF: (i) residents, (ii) HCWs, and (iii) general population. HCW were classified as such when they self-identified as HCW in an LTCF during contact tracing. The general population coming into contact with residents consisted of visitors, non-healthcare LTCF staff, and possibly HCWs in case of missing information regarding profession. We assumed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\beta\:$$\end{document} to differ between these groups because of differences in contact rates and infection probability upon a contact between an infectious and susceptible individual. Adding this into Eq. 1, the daily probability for a resident to become infected by either of three groups is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:P=1-{e}^{-{\sum\:}_{i=1}^{3}\left({\beta\:}_{i}\:\frac{{I}_{i}}{{N}_{i}}\:\text{*}\:{\text{F}\text{r}}_{\text{i}}\right)\varDelta\:t}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)$$\end{document} is the transmission rate parameter for group i, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{I}_{i}$$\end{document} is the number of infectious people for group i, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{N}_{i}$$\end{document} is the population size for group i. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\text{F}\text{r}}_{i}\:$$\end{document}  represents the fraction of overall resident contacts that are with individuals of group i. Inclusion of this fraction ensures that group-specific contact rates do not exceed the total contact rate of residents. We set these fractions ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:F{r}_{i}$$\end{document} ) to be 1/3 such that the sum across the three groups equals 1. The transmission rate parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)$$\end{document} and contact fraction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:F{r}_{i}$$\end{document} are not identifiable individually. We therefore report their product. We explored the impact of using unequal contact fractions (section ‘Sensitivity Analysis’). We allowed group-specific transmission rate parameters to differ between periods (see Supplementary Material Methods).

Applying Eq. 2 to all LTCFs (f) and for all time points (t), gives the likelihood function with θ representing the parameters, z representing the data (i.e. number of resident cases (Cf, t) and number of residents that can get infected (Sf, t) per LTCF per time point) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{P}_{f,t}$$\end{document} the probability for a resident to become infected:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\:P\left(z|\theta\:\right)=\prod\:_{f=1}^{F}\prod\:_{t=1}^{T}{{P}_{f,t}\left(\theta\:\right)}^{{C}_{f,t}}\\&*\:{(1-{P}_{f,t}(\theta\:\left)\right)}^{{S}_{f,t}-{C}_{f,t}}\end{aligned}$$\end{document}

Inference

To estimate the transmission rate parameters by group and over time, we fitted a generalized linear model (GLM) to the data in a Bayesian framework. To this end, we linearized Eq. 2 using a complementary log-log transformation. The dependent variable of this model represents the number of cases among residents (C) divided by the total resident population that can get infected (S) (i.e., the total population minus ongoing infections), which is the maximum likelihood estimate for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:P$$\end{document} from Eqs. 2 and 3. Models with and without interaction terms between group and period were compared using the Widely Applicable Information Criterion (WAIC) [25]. The following equation represents the full model, including all interaction terms:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}\:\text{log}\left(-\text{log}\left(\frac{C}{S}\right)\right)={c}_{0}&+{c}_{{HCW}}*HCW\:\\&+\:{c}_{gen\:pop}*GEN\:POP\\&+{c}_{p}*PERIOD\\&+{c}_{gen\:pop\_period}\text{*}PERIOD\text{*}\:GEN\:POP\\&+{c}_{HCW\_period}\text{*}PERIOD\text{*}\:HCW\\&+\:\text{O}\text{F}\text{F}\text{S}\text{E}\text{T}\end{aligned}$$\end{document}

with c denoting the regression coefficients. Here, c_0_ denotes the intercept, with residents (res) and period A being the reference categories. Other coefficients relate to HCW (hcw), the general population (gen pop), period (p), and the interactions of period with the general population (gen pop_period) and with HCW (HCW_period). The mean prevalence, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\text{l}\text{o}\text{g}\left(\sum\:_{i=1}^{3}(\frac{{I}_{i}}{{N}_{i}}\text{*}{\text{F}\text{r}}_{\text{i}}\right)\:)$$\end{document} , denotes the offset variable [26]. The capitalized explanatory variables for each group consist of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:\frac{{I}_{i}}{{N}_{i}}\text{*}F{r}_{i}/\sum\:_{i=1}^{3}(\frac{{I}_{i}}{{N}_{i}}\text{*}{\text{F}\text{r}}_{\text{i}}\left)\right),$$\end{document} with i denoting the group and residents in period A acting as the reference group.

The daily prevalence estimates for residents and HCW were calculated at LTCF-level, whereas the daily prevalence in the general population was calculated at regional level and applied to all LTCFs within that region. For the general population, we assumed an underreporting of cases in the general population of about 50% and therefore multiplied reported case numbers by two (see Supplementary Material methods for more information and the section ‘Sensitivity Analysis’ for an exploration of the impact of this assumption). The data described the timing of symptom onset or diagnostic test result, not when infections occurred and how long they persist, which is needed to calculate prevalence. We used data augmentation techniques [27, 28] to account for uncertainty in the dates of infection within the inference framework thereby creating 100 different datasets. For the aggregated number of cases detected in the general population, no data augmentation was performed, but rather we imputed the estimated day of infection. Lastly, for those cases with missing information on LTCF identifier we imputed this based on the LTCF associated with the cluster they were part of. For those clusters that were associated with multiple LTCFs, the LTCF was randomly selected amongst these for each of the datasets to account for this uncertainty. Further details on the applied augmentation techniques are described in Supplementary Material Methods.

For each of 100 augmented datasets, we fitted Eq. 4 in a Bayesian GLM framework with a Hamiltonian Monte Carlo sampling algorithm of 2000 iterations, discarding the first 1000 as burn-in period. Uniform priors across any real number were used for the regression coefficients. A gamma distributed prior (alpha = 0.01, beta = 0.01) was used for the overdispersion parameter and a beta distributed prior (alpha = 1, beta = 1) for the zero-inflation parameter. We calculated WAIC values [25] for each model run (i.e., each dataset) and averaged these to assess model fit. WAIC values were used to compare model variants assuming binomial, beta-binomial, zero-inflated binomial, and beta-binomial zero-inflated distributions and provide insight into the extent to which superspreading events played a role in LTCFs (i.e., resulting in overdispersion and/or zero-inflation in the data). Further analyses were done on pooled results from all model runs. All analyses were conducted in R version 4.2.1., model fitting was done using the brms package [29].

Period-specific transmission rate parameters from each group towards residents ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)$$\end{document} ) can be derived from the regression coefficients. For example, the transmission rate parameter from the general population towards residents in period C was calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{e}^{{(c}_{0}+{c}_{gen\:pop}+{c}_{period\_c}+{c}_{gen\:pop\_period\_c})}$$\end{document} . These regression coefficients can also be used to estimate the susceptibility of residents during a specific period relative to period A and the infectivity of HCW and the general population relative to residents. Relative resident susceptibility for period i was estimated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{e}^{{c}_{p\_i}}$$\end{document} . Relative infectivity of the general population was estimated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{e}^{({c}_{gen\:pop}+{c}_{gen\:pop\_period})}$$\end{document} and similarly for HCW: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{e}^{({c}_{HCW}\:+\:{c}_{HCW\_period})}$$\end{document} . In this analysis, we corrected for an increased transmission success of the Alpha (1.7 times higher [30]) and Delta (2 times higher [31]) variants, relative to the circulating variant (Supplementary Material Methods). By doing so, we aimed to disentangle temporal changes in susceptibility due to new variants from external factors such as natural immunity and vaccination. Uncorrected results are presented in supplementary Fig. 11.

Sensitivity analyses

We conducted three sensitivity analyses. Firstly, estimating the number of new cases requires information on the number of residents that can get infected for each day at each LTCF (represented by S in Eq. 4). This number was defined as the LTCF capacity minus the number of infectious residents. We thus do not explicitly correct for the build-up of natural immunity. Rather, the build-up of immunity is reflected in the estimates of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)$$\end{document} : a lower probability for an infectious contact to result in infection (p) reflects a less susceptible resident population, e.g. because of natural or vaccine-induced immunity. In a sensitivity analysis, we explored the impact of the assumption regarding the size of the resident population that can get infected. For this, we modelled immunity after infection as perfect and lasting, thereby reducing the estimated population that can get infected over time. In this analysis, the number of residents that can get infected was defined as the LTCF capacity minus the cumulative number of infected residents in the LTCF. Vaccine-induced immunity, and changes over time therein, remains reflected in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)$$\end{document} estimates. Additionally, we conducted a second sensitivity analysis to assess the impact of underreporting of infections in the general population by increasing the underreporting rate to 75%. Finally, in a third sensitivity analysis, we explored the impact of assuming an unequal distribution of resident contacts across the three groups ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\text{F}\text{r}}_{i}$$\end{document} ) where 30% of resident contacts are with other residents, 60% with HCW, and 10% with the general population.

Results

Descriptive analyses: LTCF-associated cases

Over the course of the study period (1 October 2020 to 10 November 2021) a total of 36,877 SARS-CoV-2 cases were registered among LTCF residents and 19,676 among HCWs (Fig. 1A). These cases were distributed over 1787 LTCF locations in the Netherlands. The median age of residents was 86 years (interquartile range (IQR): 80–91 years) and the median age of HCWs was 46 years (IQR range: 30–56 years). Both the resident and HCW populations consisted mainly of women, with 91% of HCW and 67% of residents. Of the 5791 recorded deaths in LTCF-associated cases, 6 (0.1%) were among HCWs. Survival of infected residents differed between periods, with the highest survival in periods C and D (Fig. 1G). The number of cases among residents and HCWs was strongly correlated (R^2^ = 0.80) (Fig. 1B). This correlation was lowest in period E (R^2^ = 0.17). With on average 2,295 cases per week, case numbers were highest in period A. Cases among LTCF residents were weakly correlated to the number of detected infections in the general population (R^2^ = 0.27) (Fig. 1C & D). This correlation was strongest in period E (R^2^ = 0.83) and weakest in period D (R^2^= −0.01). When comparing a lead and a lag of one, two and three weeks, we found that a one week lead of LTCF cases compared to the general population showed the strongest correlation (R^2^ = 0.32).

Descriptive analyses: outbreaks among residents in LTCFs

Resident cases for whom the LTCF identifier (68% of residents, n = 25,086) was recorded could be grouped into outbreaks. Outbreaks were defined as three or more cases in the same LTCF notified within 14 days of each other. Of detected resident cases that could be matched to a specific LTCF, 88% were part of an outbreak. Overall, a total of 1984 outbreaks were detected in 1297 (73%) LTCFs, with 55% (n = 980) of LTCFs reporting more than one outbreak during the study period (> 2: 23%, > 3: 6.7%). The maximum number of outbreaks reported in a single LTCF was six, which was reported by four LTCFs. Numbers of outbreaks per week across LTCFs varied over time, with most outbreaks per week occurring in period A (77 per week) and fewest in period C (7 per week). The median outbreak size was 7 (IQR 4–13) and lasted for a median of 14 days (IQR 7–23).

Across all outbreaks, 63.0% were considered large outbreaks (i.e., > 5 cases), and 8.8% of outbreaks resulted in more than 30% of residents becoming infected (Fig. 1E). There was no correlation between outbreak size and LTCF size (R^2^ = 0.02). Outbreak size varied significantly by period, with the largest median outbreak size observed in period A (median size = 9) and the lowest in period C (median size = 4) (Kruskal-Wallis rank sum test, p < 0.001) (Fig. 1F). Fig. 1. Descriptive results. A Epidemic curve based on reporting date in residents and HCW. Vertical lines and letters indicate study periods. B Correlation between daily number of reported resident and HCW cases. The fitted linear regression model is shown with its corresponding R^2^ value. C Epidemic curve based on reporting date in the general population. Vertical lines and letters indicate study periods. D Correlation between daily number of reported resident and general population cases. The fitted linear regression model is shown with its corresponding R^2^ value. E Cumulative density plot of outbreak sizes. F Outbreak sizes per period. Dots represent the median values. G Mortality hazard rate relative to period A

Model-based analyses: estimation of group-specific transmission rates

A subset of the individual-level data that contained bed capacity data (comprising 466 of the 1787 LTCFs; 26.1%) was used to estimate the contributions of different groups to the FOI in LTCF residents over time. For this subset, the median bed capacity of LTCFs was 71.5 (IQR 42–114). The total number of residents across these LTCFs was 39,150. LTCFs with available bed capacity were not equally distributed across the 24 GGD regions, with the number of LTCFs per region ranging from 4 to 38 (Pearson’s chi-squared test: p < 0.001).

We used a Bayesian approach to estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)\:$$\end{document}  for each group towards residents and quantify the contribution of each group to the overall FOI in residents. We compared model variants assuming different underlying distributions in the data and found a beta-binomial distribution to provide the best representation of the data (Supplementary Fig. 3 for a comparison of WAIC estimates). A model including a period term with only the significant interaction terms, which were the general population with period C and with period D, resulted in the lowest WAIC and was retained for further analysis (Supplementary Fig. 4).

The model presented a good fit to the observed data (Fig. 2A, Supplementary Fig. 5) although the sizes of the largest outbreaks were underestimated (Supplementary Fig. 6). After correcting for the increased transmission potential of Alpha and Delta variants, we found that resident susceptibility was highest in period A and lowest in period C (62% reduction from period A), after which it increased to 57% of the initial susceptibility in period E (Fig. 2B). Infectivity of HCWs was 32% lower compared to residents. Infectivity of the general population varied between periods, and was lowest in period C (59% lower than residents) (Fig. 2D). We used the uncorrected estimates to calculate the transmission rate parameter multiplied by the contract fraction ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)\text{*}{\text{F}\text{r}}_{i})\:$$\end{document} from each group towards residents (Supplementary Table 2). This value was 0.09 (95%CI 0.08–0.09) in period A, which declined as COVID-19 vaccines were rolled out in period C (0.06 95%CI: 0.04–0.08), after which it increased again up to 0.10 (95%CI: 0.09–0.11) in period E. A similar initial decrease was seen in HCWs with a 19% reduction between period A and B, and a 22% reduction between period B and C. The largest observed difference was a 71% reduction between period B and C in the general population. The values were robust to assumptions regarding the fractions of resident contact with each of the three group (Supplementary Table 5). The effective reproduction number of resident-to-resident transmission (calculated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)*{Fr}_{i}*\:\:$$\end{document} * duration of infectious period (Supplementary Material Methods)*) was below 1 during the full study period, suggesting that resident-to-resident transmission was self-limiting.

The overall daily FOI towards residents was highest in December 2020 (1.710^ −3^, 95%CI: 1.510 ^−3^–1.910 ^−3^) and lowest in June 2021 (1.110^ −5^ 95%CI: 7.610 ^−6^–1.710 ^−5^) (Supplementary Fig. 10). We inferred that, in the pre-vaccination period (period A), 42.1% (95%CI 34.6% − 48.8%) of the FOI experienced by residents could be attributed to fellow LTCF residents (Fig. 2C). The contribution of the general population during this period was similar (40.8% 95%CI 35.1% − 46.9%). The contribution of residents to the FOI rapidly declined with the rollout of vaccination, with only 1 in 4 (23.7% 95%CI 11.6%−37.4%) infections in residents likely to result from another resident in period C.

HCW contributed least to the FOI across all periods, with a highest contribution of 19.7% (95%CI 13.2%−26.5%) in period B to the lowest in period C (5.8% 95%CI 1.7% − 11.2%). Like residents, the contribution of HCWs to the FOI declined with the rollout of vaccination, but for HCWs this did not return to pre-vaccination levels. The general population was the main source of infections in residents during period C (70.6%) and D (60.2%), when the absolute FOI was lowest and the uncertainty in estimates larger. These are the periods where residents and HCW had been vaccinated, but vaccination of the general population was still being rolled out.

Sensitivity analyses

When comparing the assumption of permanent, complete immunity in residents to assuming imperfect immunity in residents (default), we found that the values of the transmission rate parameter multiplied by the contact fraction ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)*{Fr}_{i}$$\end{document} ) were on average 24% higher (supplementary Table 3). Similarly, the mean FOI was on average 20% higher (supplementary Fig. 13). However, the relative FOI, relative resident susceptibility and relative infectivity of the other groups were robust to assumptions around immunity (supplementary Fig. 12). Assumptions regarding the extent of infection underreporting in the general population affected the transmission rate parameters and relative infectivity related to this group, but none of the other outcome measures (supplementary Table 4, supplementary Figs. 15 & 16). Variation in the distribution of resident contacts across the three groups affected estimates of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)$$\end{document} , but estimates of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\:{\beta\:}_{i}\left(t\right)*{Fr}_{i}$$\end{document} remained the same (supplementary Table 5), as did the relative and absolute FOI (supplementary Fig. 17 C, supplementary Fig. 18). The estimates for infectivity of the general population and HCW relative to residents were sensitive to assumptions around the contact fractions (supplementary Fig. 17D).

Fig. 2. Model fits and outcomes for the best-fitting model. A Comparison between the daily number of observed infections (in the data) and predicted infections (from the model). Values of the 7-day average are shown. Dark lines indicate mean values. Shaded areas show 95% prediction intervals obtained by generating 1000 predictions by randomly sampling from the posterior distributions and augmented datasets. B Relative susceptibility to infection in residents. The horizontal line at y=1 indicates no difference compared to period A. C Daily relative force of infection from each LTCF group towards residents. Dark lines indicate mean values. Shaded areas show 95% prediction intervals obtained by generating 1000 predictions by randomly sampling from the posterior distributions and datasets. D Relative infectivity compared to residents. The horizontal line at y=1 indicates an infectivity equal to that of residents. Estimates for the infectivity of the general population differed in periods C and D compared to all other periods as these were significant interaction terms in the statistical model

Discussion

We used a Bayesian GLM fitted to LTCF resident case data from October 2020 to November 2021 to quantify the contribution of residents, HCWs, and the general population to the force of infection in LTCF residents. We described characteristics of the COVID-19 epidemic in Dutch LTCF settings, with a focus on outbreaks among LTCF-associated cases. We observed the highest incidence rates as well as largest outbreaks in the pre-vaccination period. Vaccination rollout was followed by a stark reduction in the average estimated infection susceptibility of residents. This reduction was robust to strong assumptions on the buildup of natural immunity, indicating that this reduction could not be explained by natural immunity alone. This protection decreased over time, which could partly be explained by increased transmission success of newly emerging Alpha and Delta variants. HCWs contributed least to the FOI towards residents over the study period, based on the information available regarding who is a LTCF HCW. Contributions from residents and the general population were similar in the pre-vaccination period, however, the general population formed the main source of infections in LTCF residents in the period after vaccination was rolled out in LTCFs, while vaccination was not yet widely available for the general population. This suggests a substantial reduction in within-LTCF transmission after vaccination. Resident-to-resident effective reproduction numbers were below 1 during the whole study period, meaning that resident-to-resident transmission was, on average, self-limiting.

The observed temporal trends in FOI contributions could be explained by several processes. During early spring 2021 (period B & C) Alpha became the dominant strain, followed by Delta in summer 2021 (period D & E) [32]. We detected no increase in transmission in LTCFs during the emergence of Alpha, which coincided with the rollout of vaccination in residents and HCWs. Estimates of temporal trends in residents’ susceptibility suggest that vaccine-induced immunity in the resident population was sufficient to counter the increased transmission success of the Alpha variant. Resident-directed transmission rate parameters increased roughly six months after the first round of vaccination. This increase could be due to a combination of waning immunity [33] and the Delta variant becoming the dominant strain, which was associated with increased transmission success [34] such as through immune escape [35]. After correcting for the estimated transmission advantage of the Delta variant, estimates of residents’ susceptibility still increased towards the end of the study period, but remained low compared to those of the pre-vaccination period. These findings suggest that both waning immunity and the new variant played a role. Additionally, trends in transmission rate parameters could also be a consequence of seasonality [36], variation in case detection rates, and variation in contact between residents and all three groups due to additional control measures. Policies and the levels of adherence to control measures could have differed between LTCFs, which makes it difficult to study the effect of these interventions when looking at overall results. If detailed information on interventions and adherence to control measures would be available at LTCF or regional level, stratified analyses could shed light on their impact. Similar to observed trends in transmission rate parameters, mortality rates also started increasing towards the end of the study period, indicating that protection against mortality had also waned, possibly due to an increased mortality risk associated with the Delta variant, although its impact on mortality remains uncertain [37–39]. Mortality was recorded as all-cause mortality in the present study and likely suffered from underreporting due to the absence of a notification requirement. By estimating mortality rates relative to period A, we aimed to reduce the impact of this, although temporal changes in underreporting can bias our findings. Another data limitation that might have affected our results is the incompleteness of LTCF identifier for residents. While the model showed a good fit to population-level case trends, fits at LTCF-level showed an underestimation of case numbers for LTCFs with large outbreaks. This could be due to infected residents with missing LTCF identifier being incorrectly assigned to a different LTCF during data augmentation (see Supplementary Information Methods) or due to suboptimal linking of LTCF identifier with its capacity information.

Our results indicate that a large proportion of infections in residents originated from the general population. Several other studies discussed the contribution of different groups to transmission in LTCFs, with significant variation in their findings. Genetic analyses of a large outbreak in a Dutch LTCF in the early phase of the pandemic indicated multiple introductions into the resident population, with a limited role for within-LTCF transmission [40]. Similarly, an outbreak in a Belgian nursing home originating from a single visitor showed a significantly increased infection risk in residents who attended the visitor’s event and limited onward transmission to those not attending [41]. Two model-based analyses from the United Kingdom [42, 43] also suggested the importance of external introductions, specifically by HCWs. Our results suggest a relatively small contribution of HCWs to the FOI towards residents. However, underreporting of infections in HCW is likely, e.g., due to misclassification of HCW status. An unknown part of the transmission from the general population should therefore probably be ascribed to HCWs or to other staff. Additionally, due to missing information about in which LTCF each HCW worked, we imputed this information when possible. This misclassification contributed to the uncertainty in HCW’s contribution to the FOI. Moreover, it has been shown that adherence to the reactive testing policy was higher for residents than for HCW, which indicates that asymptomatic infections in HCW were more likely to be missed [22]. It is also important to note that interventions such as masking were in place during the study period which could have contributed to the limited role of HCWs to transmission. Relative contributions to resident infections may have been different during the first wave, which was not in our study period, when fewer interventions were implemented.

Interventions are generally targeted towards reducing introductions into an LTCF, reducing transmission within an LTCF, or reducing the impact of infection. We showed that the estimated susceptibility of residents declined when vaccination in residents and HCWs was rolled out. This indicates that vaccination contributed to reducing introduction risk, limiting within-LTCF transmission, and reducing the impact of infections. This observed reduction in susceptibility, rather than in infectivity, aligns with several studies showing that virus shedding was similar in vaccinated and unvaccinated individuals [44–46]. We only found statistical support for a reduction in infectivity in the general population relative to residents in periods C and D (March-October 2021). While the general population tended to be more infectious than residents in periods A, B, and E, they were found to be less infectious than residents in periods C and D. The reason for this is unclear, but could be due to behavioral changes related to vaccination status or higher vaccine efficacy as the general population was generally younger [47, 48]. While the temporal trends of relative infectivity were robust to assumptions around the distribution of resident contacts, its absolute values were sensitive to this distribution. After the vaccination rollout in LTCFs, the general population was the dominant source of infections in residents and resident-directed Rs dropped below one. The effectiveness of control strategies depends on how important the targeted infection source is for the overall infection load in residents. The variation in dominant sources of infection, therefore, implies that the most effective control strategy for reducing infections in LTCFs varied over the course of the epidemic.

While we applied our model framework to LTCFs in the Netherlands, this approach could be easily adapted to other countries or similar settings where housing and healthcare are combined. Using a Bayesian approach combined with data augmentation allowed for a rigorous quantification of uncertainty and we showed a good model fit to the data. During a large-scale outbreak, dates of infection and case characteristics, such as an LTCF identifier, are often not observed or recorded. Data augmentation allows to explicitly account for these uncertainties in the analysis, making results robust to the impact of missing information. The availability of contact survey information for LTCF residents made it possible to adjust the prevalence in the general population to the prevalence in resident contacts (supplementary Fig. 2). While many other studies are limited by a lack of information on asymptomatic infections, the reactive testing policy in LTCFs resulted in less underreporting and allowed for detection of asymptomatic infections, although the adherence to this policy varied [22]. Asymptomatic individuals have been shown to contribute to transmission [49] but are challenging to detect, particularly in LTCFs, making this an important strength of this study. This was particularly important in later periods, when a larger share of infections lead to mild symptoms due to immune protection and are therefore more likely to be missed using symptom-based testing. Case affirmation rates among the general population are uncertain and may have varied due to variation in testing behavior over time and by age. To account for this, we assumed an underreporting proportion of 50% and restricted the study to the period where testing was available to everyone at sufficient capacity. Sensitivity analyses showed that deviations in this proportion affect transmission rate parameter estimates for this population, but does not affect relative and absolute contributions to the resident-directed FOI from this group. A potential limitation is the incompleteness of data on LTCF capacity. However, there is no indication that this was related to regional or within-LTCF SARS-CoV-2 dynamics as this information predates the epidemic.

Conclusions

When designing interventions to reduce transmission it is important to characterize the epidemiological situation to inform the most effective strategy and avoid unnecessary harmful side-effects. We observed substantial variation over time in relative contributions of residents, HCWs, and the general population towards resident infections in LTCFs. This suggests that the effectiveness of potential interventions, which were focused on preventing introductions or within-LTCF transmission, varied over the course of the outbreak. This is valuable information that can be used for decision-making, such as the scheduling of vaccination rounds or the implementation of visitor restrictions to lower introduction risk. The ability to generate such insights relies heavily on robust surveillance systems, highlighting the importance of continued investments in these efforts. The application of advanced statistical techniques to routinely collected data opens opportunities to tailor interventions to the current epidemiological context and improve data-driven decision-making.

Supplementary Information

Supplementary Material 1.