Estimating time-dependent contact: a multi-strain epidemiological model of SARS-CoV-2 on the island of Ireland

Mathematical modelling plays a key role in understanding and predicting the epidemiological dynamics of infectious diseases. We construct a flexible discrete-time model that incorporates multiple viral strains with different transmissibilities to estimate the changing patterns of human contact that generates new infections. Using a Bayesian approach, we fit the model to longitudinal data on hospitalisation with COVID-19 from the Republic of Ireland and Northern Ireland during the first year of the pandemic. We describe the estimated change in human contact in the context of government-mandated non-pharmaceutical interventions in the two jurisdictions on the island of Ireland. We take advantage of the fitted model to conduct counterfactual analyses exploring the impact of lockdown timing and introducing a novel, more transmissible variant. We found substantial differences in human contact between the two jurisdictions during periods of varied restriction easing and December holidays. Our counterfactual analyses reveal that implementing lockdowns earlier would have decreased subsequent hospitalisation substantially in most, but not all cases, and that an introduction of a more transmissible variant - without necessarily being more severe - can cause a large impact on the health care burden.


Introduction
During an epidemic, behavioural changes are encouraged, and sometimes mandated, to curtail infectious disease transmission. These changes aim to reduce the number of contacts between people broadly (e.g., closure of schools, workplaces, commercial establishments, roads, and public transit; restriction of movement; cancellation of public events; maintenance of physical distances in public) and reduce the chance of infection upon contact (e.g., use of personal protective equipment). Furthermore, tracing and isolating known infectious cases can limit the contact between infectious and susceptible individuals. These actions are collectively referred to as non-pharmaceutical interventions (NPIs) and are often mandated by governments. Slowing the surge of infection (or "flattening the curve") affords an opportunity to reduce infection-induced mortality and morbidity, alleviate health care burden and wait out an epidemic until pharmaceutical solutions (i.e., treatment and vaccines) become available. Implementation of mandated NPIs in historical outbreaks, including during the 1918 influenza pandemic, was crucial for preventing excess death in the United States [1]. NPIs have also been mandated globally during the COVID-19 pandemic.
Mathematical modelling and quantitative analyses of empirical data play a pivotal role in understanding and predicting epidemiological dynamics. Mechanistic epidemiological models have been widely applied to study the dynamics of SARS-CoV-2, and to make predictions of clinical outcomes under alternative scenarios (e.g., an assumed decrease in physical contact [2]). Despite their public health benefits, social distancing measures have been shown to incur high costs in several domains, including in economy [3], mental health [4], and civil liberty [5]. Thus, it is crucial to quantify infection contact, or its derivative quantities like the effective reproductive number, R -to monitor changes in infection burden, achieve desired public health outcomes and improve policy transparency and public engagement. While it is not possible to measure human contact (and its effect on disease transmission) directly, fitting a mathematical model to longitudinal data on observed processes such as reported cases and hospital admissions allows estimation of human contact and its derivatives [6,7].
Many epidemiological models follows a rich tradition of ordinary differential equation (ODE) models [8], which track the spread of infection and often immunity in a population. Specifically, ODE models assume that waiting time processes (such as infectious period and time to hospitalisation) are memoryless, that is to say, that the waiting time until an event (such as recovery and hospitalisation) does not depend on the elapsed time. Seen at the population level, this assumption deduces that times spent by individuals in each compartment are distributed exponentially, implying large individual variability. While mathematically convenient, the lack of memory is unsupported for certain epidemiological processes [9], and empirical evidence indicates other probability distributions with smaller individual variability and nonmonotonic densities (e.g., gamma-, Weibull and log-normal distributions) are better equipped to describe those processes. Previous studies have also demonstrated that quantitative predictions of epidemiological outcomes depend on the assumed probability distribution in a variety of systems [10][11][12][13], including SARS-CoV-2 [7]. As such, it is pertinent to incorporate realistic waiting time distributions, particularly when one aims to obtain quantitative and short-term, rather than qualitative and long-term insights from epidemiological models.
Here we develop a compartmental epidemiological model that accurately predicts inter-individual contact over the first year of the epidemic in the Republic of Ireland (ROI) and Northern Ireland (NI). These neighbouring jurisdictions on the island of Ireland present a compelling contrast due to independent policymaking over a small geographical area. We use a flexible discrete-time approach that incorporates waiting time distributions that reflect accurate assumptions about times spent in each compartment [7].
The COVID-19 pandemic has been characterised by subsequent waves of novel variants with varying disease transmissibility and severity. To separate the effect of human behaviour from the difference in transmissibility of multiple strains, we explicitly model multiple SARS-CoV-2 strains, with differing transmissibilities, seeded in the population at differing times. Our multiple strain model allows consistent estimates of relative human contact between periods even when the dominant variant has changed.
Previous studies of SARS-CoV-2 have estimated the time-dependent contact ratio and its derived quantities using continuous (e.g., basis splines [14]), or piece-wise, discrete functions of time (e.g., consisting of the specified period corresponding to NPI mandates [7]). However, both approaches can be fraught with challenges. On the one hand, it is not obvious to choose the appropriate extent of smoothing of a continuous function, for example, by deciding the number of knots in a basis spline function. On the other hand, abrupt changes imposed by piece-wise functions are at odds with empirical data on human movement during the COVID-19 pandemic [15]. Furthermore, it is difficult to establish a precise definition of the level of an intervention over time as definitions changed over time [16,17]. To address these concerns, we develop an intermediate approach, in which we introduce a prior that allows smoothness in human contact between neighbouring weeks in the absence of information otherwise from empirical data.
In the Irish context, compartmental models have been used by several other research groups to understand the dynamics of the virus, make forecasts of outcomes under various scenarios, and assess economic impacts of policy restrictions [18][19][20][21][22][23]. Our study complements these studies by providing a high-resolution description of the change in human contact over time, comparing the two jurisdictions on the island of Ireland. Leveraging the epidemiological model and estimated parameters, we also perform counterfactual analyses to explore the effects of alternative interventions on cumulative hospitalisation and assess the impact of a novel variant.

Epidemiological model and data
Our multi-strain discrete-time model consists of three types of host compartments (Fig. 1): a susceptible compartment (S) and two infectious compartments for viral strain s, (J s, i and Y s, i where i indicates the infection age, i.e., day since exposure). The compartments J and Y differ in their future clinical outcome: individuals in the components Y eventually get hospitalised while those in J remain out of hospitals. As our primary focus is the inference of NPI in the community, we did not consider within-hospital transmission, recurring hospital admissions of the same patients, or demographic turnover (including death). We ignored the dynamics of recovered hosts who were assumed to have minimal influence on the transmission during the period investigated.
Our model is parameterised by θ, the probability of hospitalisation, Δ s , the daily probability of infection with strain s, and a discrete random variable H that characterises a set of probabilities governing daily transitions to hospital. Δ s is informed by a discrete random variable, Z, that characterises a set of probabilities governing daily transitions into infected compartments, and is defined in Section 2.1.1. Z denotes the time in days from exposure of the infector to exposure of the infectee for a randomly chosen infectee-infector pair (i.e., generation interval), and can be viewed as the average relative contribution of each day to the individual reproduction number. The probability that infection occurs at infector age i, Transmission of infection does not affect the individual's stay in the compartment; however, for transition out of the Y compartments, the event going to hospital at infection age i is conditional on still being in the compartment at infection age i − 1. Given H, the random variable representing time from infection to hospital admission in hospitalised patients, we denote by η i the probability of hospitalisation at infection age i given the individual was still not hospitalised at infection age i − 1, . This is the discrete hazard of hospital admission at infection age i. We use published estimates for Z and H, as described in section 2.1.3 below. Our model assumes that the proportion of infected people hospitalised and the processes governing hospitalisation and recovery over time are constant across strains.

Infection dynamics
We extend a discrete epidemiological modelling framework by Sofonea et al. [7] to accommodate multiple viral strains spreading simultaneously. First, we express the effective density of infectious host Once infected, individuals progress to the next square each day (J and Y), capturing the memory effect of the infection age. After spending n j days, infectious hosts in J are no longer infectious. Alternatively, a fraction, θ of infectious hosts (in Y), is admitted to the hospital with a delay specified by the probabilities η 1 , …, η n y , where η i is the probability that the individual is admitted to hospital on the day i, conditional on their being infectious for i − 1 days. The grey arrows indicate the daily transition of individuals that occurs with probability 1.
population on a given day d that contributes to transmission of the strain s as: is the number of individuals with strain s on day d with infection age i in the community. Multiplying by ζ i and summing over infection ages, can be regarded as the total 'potential for infection' in the community on day d. The effective infectious density, I s (d), is this sum scaled by the contact ratio, c(d), on day d. The contact ratio is the ratio of contact rate on day d to contact rate on day 0. The infectious density is thus a measure of the total amount of transmission in a completely susceptible population.
We allow susceptible hosts encounter the viral strain s with a probability Λ s (d), which as in [7] is assumed to follow the Michaelis-Menten function that saturates with the effective infectious density, I s (d) and the contact rate; under assumptions about initial conditions as in [7] we derive the following: where τ s is the relative transmission advantage of strain s, S 0 the population size, and R 0 1 the basic reproductive number of the original strain.
We define the τ s as the ratio of the basic reproductive numbers of strain s to the original strain R 0 s = τ s R 0 1 .
When a host encounters multiple strains, we model the interaction between strains assuming superinfection with priority determined by order of exposure: i.e., only the first strain that encounters a host establishes infection when the same host subsequently encounters multiple strains. Thus, in the case of two strains, the probability of getting infected with the strain s, Δ s (d), is: where s ′ denotes the non-focal strain. Similarly, expressions can be derived for more than two strains. It follows then that the number of susceptibles on the next day is expressed as: Of those exposed to either viral strain, the proportion θ will develop severe symptoms and eventually be admitted to the hospital (Fig. 1).
For less severe cases that do not result in hospitalisation, J, the infection progresses towards recovery until they are no longer infectious on the day n j : Those that develop severe symptoms, Y, are admitted to hospital with the probability η i on the i-th day following exposure.
It follows then that the number of hospital admissions on day d + 1 equals

Observed longitudinal data
Epidemiological models are often fitted to data on infected caseshowever, case data depends on levels of testing, which varied over time during the COVID-19 pandemic. It can also be problematic to rely on data on deaths -for instance, many deaths in the ROI occurred following outbreaks in care homes, and thus data on deaths may not reflect disease spread in the general community. Biases and uncertainty in estimating the reproductive number arising from such issues are discussed elsewhere [6]. Thus, hospital admission data is likely a better reflection of the community spread of SARS-CoV-2. We used daily COVID-19 hospital admissions in the ROI and NI, reported respectively by the Central Statistics Office COVID data hub for the ROI [24], and the NI Department of Health [25]. One potential caveat of publicly available hospital admission data is their ambiguity on whether infections were acquired in the community or health-care settings. Modelling studies of hospitalacquired SARS-CoV-2 in Germany [26] and England [27] estimate roughly 10 and 20% of hospitalised cases may have originated from transmissions within hospitals: such estimates are not available for Ireland to our knowledge. As our study focuses on community transmission alone, we conduct a sensitivity analysis fitting our model to 80% of reported hospital admission numbers.
Successive invasions of new variants have characterised the COVID-19 pandemic. Our study tracks two strains that circulated in the island of Ireland in the first 12 months of the pandemic: i.e., the original strain (initially detected in Wuhan, China) and the Alpha strain (also known as B.1.1.7., initially detected in Kent, UK). We used publicly accessible data on the frequency of the Alpha strain in the ROI [28], and NI [29], respectively.

Incorporating empirical estimates of waiting time distributions
Linking transitions within and between model compartments are two random variables, each describing a waiting time process. These are the infectious period (generation interval), Z, and the delay between infection exposure and hospitalisation, H. The probability distributions representing these random variables have been estimated elsewhere empirically for SARS-CoV-2 in a global and European context as described below.

Generation interval
The generation interval refers to the time between infection events in a pair of infector and infectee, reflecting the incubation duration and recovery timing. Here, we used the distribution of this interval to model the relationship between the age of infection (i.e., time since exposure) and the infectiousness of the infector. We employed an estimate by Ferretti et al. [30] who found the variation in SARS-CoV-2 generation interval was best described by the Weibull distribution with the mean interval of 5.5 days (shape=3.29 and scale=6.12). We truncated the Weibull distribution at the upper-integer-rounded 99%-quantilewithout this truncation, the discrete model would require infinite timetracking sub-compartments due to a right-unbounded support [0, ∞]. We then discretised the distributions because the dynamics unfold in discrete-time intervals of one day in our model: the upper limit of the discretised distribution corresponds to n j (eq. 6).

Exposure to hospital admission
The waiting time between exposure and hospital admission was estimated as the sum of the incubation period and the delay between symptom onset and hospitalisation. We assumed that the two waiting times were independent due to the absence of evidence otherwise. A meta-analysis of global, but predominately, Chinese data found that the SARS-CoV-2 incubation period was log-normally distributed with parameters μ = 1.63 and σ = 0.50 [31], corresponding to a mean incubation time of 5.78 days (standard deviation of 3.97 days). The distribution of waiting time between symptom onset and hospitalisation was estimated assuming a gamma distribution by Public Health England with a mean of 5.14 days (standard deviation of 4.2 days) [32]. We fitted a gamma distribution to the simulated sum of the two distributions to represent the timing between exposure to infection and hospital admission (shape =4.76 and rate =0.435). As for the generation interval, we discretised the distributions and the upper limit of the discretised distribution corresponds to n y (eq. 8).

Weekly contact ratio
Here, we defined the contact ratio c as the human contact rate relative to the pre-pandemic, pre-intervention baseline (eq. 1 & 2) and estimated this quantity using a piece-wise function consisting of weekly intervals. Specifically, we estimated the ratio in each area a (NI and ROI), per week w (i.e., c a, w ) as a function of ϕ a, w , the log proportional change in the contact ratio from the previous week. We index w from the date of the first public health intervention in either jurisdiction, which took place in ROI on 2020-03-12 (Supporting Information S1: Table S1 & S2); hence the preceding, pre-intervention contact ratios are defined as 1.0.
With this formulation, hierarchical Bayesian inference with priors on the ϕ a, w allows us to estimate the time-varying weekly contact ratios with minimal prior information specific to the modelled system. Specifically, we used a prior ∼ N (0, ε), where ε is a hyperparameter specifying the standard deviation of ϕ, such that c a, w would equal c a, (w− 1) in the absence of signals from epidemiological data ( Table 1). A priori, this formulation avoids over-fitting random weekly variation at the potential risk of smoothing over valid signals of an abrupt change in the weekly contact ratio, for example, following an introduction of lockdown measures. To check for such bias, we examined the extent to which our smoothing approach affects the estimation of sudden changes in the contact ratio, c. We showed that our formulation is unlikely to introduce substantial bias (Supporting Information S2).

Initial conditions
The first case of SARS-CoV-2 on the island of Ireland was identified in NI on 2020-02-27 from an individual travelling back from Northern Italy via Dublin Airport located in ROI (Table S2). Two days later, the first official case in the ROI was also confirmed from a traveller from Northern Italy (Table S1). Initially, most known cases are travel-related, and contact tracing may successfully contain infections. As our model solely tracks community transmission, we started our simulations on the first day that community transmission was detected on the island of Ireland: 2020-03-05 (Table S1). Coincidentally, the exponential growth of confirmed cases appears to have begun around 2020-03-05 in both ROI and NI [24,33]. We account for the uncertainty of the beginning of community transmission by estimating the initial infectious density independently in the two jurisdictions (Table 1). We assume implicitly that the contribution of the travel-related cases is negligible once the infection starts growing exponentially in the community.
The first cases of the Alpha strain were reported in November and December 2020, respectively, in ROI and NI (Tables S1 & S2). Due to high connectivity with the island of Britain, the Alpha strain likely entered the island of Ireland soon after it emerged in England, where the strain was detected in mid-September [34]. By February 2021, the Alpha strain comprised the majority of infections in both ROI and NI. To estimate the date of introduction, we fitted a three-parameter logistic function to the longitudinal data of the Alpha frequency and identified the date on which Alpha cases (frequency of Alpha × known new cases) intersects 1: the date of introduction was estimated as 2020-09-22. Again, we account for the sensitivity of the timing of introduction by estimating the founding infectious density of the Alpha strain, independently in the two areas (Table 1). In our model, viral strains differ only in their transmissibility, τ s .

Fitting
We used a Bayesian approach to fit the above model to two types of longitudinal data from the ROI and NI: daily counts of hospital admissions and the Alpha strain frequency. Model parameters are detailed in Table 1. Hospital admissions per day were modelled as log-normally distributed with standard deviation parameters σ h . We set the probability of hospital admission given infection, θ, to the observed figure in ROI published by the Health Service Executive [42] ( Table 1). The frequency of the Alpha strain was fitted assuming the beta proportion distribution with a standard deviation parameter, σ f . We fitted our model to the data from the first year of the pandemic from the first confirmed case of community transmission on the island, which was detected on 2020-03-05, in ROI, to the end of February 2021. Our modelling period precedes the widespread administration of the full course of vaccination in either jurisdiction: the proportion of fully (twice) vaccinated individuals in ROI and NI was less than 3% and 2% at the end of February 2021, respectively [24,33].
Our model was written in Stan 2.21.2 and fitted through the RStan interface [35]. We fitted the model in parallel in four independent chains, each with 5000 sampled iterations and 1000 warmup iterations. For diagnostics, we confirmed over 400 effective samples and ensured convergence of independent chains (R < 1.1) for all parameters [36]. We assessed the goodness of fit to data using standardised residuals (Supporting Information S3). We also quantified the posterior z-score and posterior contraction to examine the accuracy and precision of posterior distributions and the relative strength of data to prior information [37] (Supporting Information S4).

Counterfactual analyses
Estimating human contacts with a multi-strain model separates the effect of human behaviour from the difference in transmissibility of multiple strains. This separation allows us to leverage the epidemiological model and estimated parameters to simulate an epidemic based Table 1 Description of model parameters and their fixed values, or prior distributions used in Bayesian statistical inference. We assigned an informed prior for R 0 , τ 2 and a generic, weakly informative prior for I s,a (0), ε and measurement error parameters. half-N (0, 1) on data-generating processes consistent with the observed data. In turn, we can modify one part of the fitted modelwhile everything else is constantto conduct counterfactual analyses, which allows us to explore the impact of different factors that affect disease transmission.
Here, we explored two counterfactual scenarios: to examine the effect of lockdown timing; and to isolate the impact of the more transmissible Alpha strain on the hospitalisation outcome.

Effect of lockdown timing
We explored the impact of the timing of lockdown introduction by simulating an epidemic with parameters estimated from the model fitted to the observed data on hospitalisations and strain proportions, but the contact ratios counterfactually shifted earlier by seven and 14 days relative to the actual start of the three lockdowns imposed in the ROI and NI. We then compared counterfactual scenarios and reality by computing the percentage difference of the cumulative hospital admission numbers for the subsequent days under the counterfactual versus the observed scenarios.
Suppose the intervention is to shift the first lockdown date earlier by seven days. Denoting the actual lockdown date d l and the counterfactual contact ratio on day d as c*(d): The counterfactual infectious density on day d follows from eq. 1, and we denote this where the J*, Y* denote the counterfactual numbers in these compartments on day d. From this follows the counterfactuals on day d Λ s *(d), Δ s *(d), and H*(d + 1) from eqs. 2 to 9. For the second and third lockdowns, we assume that the epidemic had proceeded as observed up to the second and third lockdown, respectively.

Impact of a more transmissible variant
We investigated the extent to which the introduction of the more transmissible Alpha strain contributed to public health burden by simulating an alternative epidemic with parameters estimated from the model fitted to the observed data on hospitalisations and strain proportions, but without introducing the Alpha strain in September 2020. We then used the percentage difference to compare the cumulative hospital admission numbers between the counterfactual and real scenarios over time until the end of the modelled period at the end of February 2021.

Estimated contact ratios
We estimated a rapid decline in contact ratios during the first month of the pandemic before a strict lockdown was implemented ( Fig. 3; Tables S1 & S2). The first lockdown started on 2020-03-28 in both jurisdictions (Fig. 3 NI-a & ROI-a). By this date, the estimated contact ratio was already down to about 60% of the pre-pandemic baseline in both jurisdictions. This finding is consistent with pre-lockdown movement alternations reported elsewhere, for example in China, Italy and New York [43]. The changes were likely driven by lighter restrictions that preceded strict lockdowns and spontaneous behavioural changes in response to increasing perceived infection risk (e.g., increased incidence within the social sphere and media coverage) [44]. During the first lockdown, human contact fluctuated only slightly.
In both jurisdictions, the easing of the first lockdown began from 2020-05-18 (Fig. 3 NI-b & ROI-b), and a long period of slow restriction easing took place during the summer months. In the ROI, the estimated contact ratio increased from June and fluctuated between approximately 70-80% of the pre-pandemic baseline in July, August and September. In NI, the contact ratio rose to a peak around the end of July. We detected higher human contact in NI than the ROI in mid-June (indicated by 95% predictive intervals of the difference excluding zero; Fig. 3; bottom panel). Of potential relevance, we note that all nonessential retail outlets were allowed to reopen earlier in NI than the ROI during this period from 2020-06-12 and 2020-06-29, respectively (Fig. 3 NI-c & ROI-c; Tables S1 and S2). Human contact in NI decreased through August but elevated again to about 90% of baseline by the end of September with no parallel increase in the ROI (Fig. 3 NI-d & ROI-d). This period corresponds to the first time primary and secondary teaching resumed in person in both jurisdictions. The increased contact in NI mirrors a trend in detected England where the September schooling reopening led to increased cases, most notably among the teaching staff [45].
Ahead of the second lockdown (Fig. 3 NI-e & ROI-e), the estimated contact ratio declined to about 60% of baseline in both jurisdictions during October. Unlike during the other two lockdowns, the contact ratios tended to increase during the lockdown period in both jurisdictions throughout November (Fig. 3; top and middle panels).
At the beginning of December in the ROI, several mitigation measures were lifted, allowing non-essential businesses, restaurants, cafes and gastro-pubs to open as well as relaxing household gathering restrictions ( Fig. 3 ROI-f; Table S1). This period coincides with an increasing trend in the estimated contact ratio, which reached about 90% of the pre-pandemic baseline the week before Christmas. In NI, on

Fig. 3. Estimated weekly contact ratios in
Northern Ireland (top) and the Republic of Ireland (middle) and differences in the contact ratio between the two jurisdictions (bottom). The three lockdown periods, corresponding to the most strict restrictions in each jurisdiction, are marked in yellow, blue and green, respectively. The letters a-g in black dots correspond to the timing of events described in the main text. The black line, and grey bands correspond to the median, the 50% (dark) and 95% (light) credible intervals. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the other hand, the lockdown remained in place almost two weeks longer (Fig. 3 NI-f; Table S2), and the estimated contact ratio reached a maximum of about 75% of the baseline value before Christmas. Our estimates indicate that human contact was substantially higher in the ROI than NI for two weeks over the Christmas period (indicated by 95% predictive intervals of the difference excluding zero; Fig. 3 NI-g & ROIg): the ROI experienced the highest per capita infection rate in the world during this period. [46]. In both NI and the ROI, the third lockdown introduced in late-December 2020 coincided with the lowest contact ratio ( Fig. 3; green) followed by the first lockdown in late March ( Fig. 3; yellow) and second lockdown in November ( Fig. 3; blue).
The above finding are based on the assumption that 100% of reported hospital admission cases originate in the community. However, recent estimates indicate up to 20% of hospital admission numbers may be attributable to transmissions that take place in health-care settings [27], though the estimates vary across counties and viral strains [26]. Our sensitivity analysisassuming only 80% of reported admission numbers were community-acquireddemonstrate little impact of this assumption on the estimated contact ratio (Supporting Information S5: Fig. S5). The difference in fitted proportion of hospital admission numbers was absorbed by a 25-30% change in the estimated initial effective infectious densities (Supporting Information S5: Fig. S6), which also affects incidence, but the priors for which were less informative than either the contact ratios or R 0 (Supporting Information S4: Fig. S4).

Counterfactual scenarios Effect of lockdown timing
Lockdown measures have been shown effective in reducing the infection burden of SARS-CoV-2, and the timing of introduction is the most significant factor in determining their effectiveness [47,48]. We found that bringing forward the lockdown dates by either seven or 14 days would have substantially reduced the cumulative hospitalisation over the subsequent 50 days from the date of lockdown in most scenarios (indicated by the 95% predictive interval excluding zero; Fig. 4). Of note, we found that a counterfactual simulation to bring forward the second lockdown date by seven days showed a non-conclusive impact on the cumulative hospitalisation in the subsequent 50-day period in either jurisdiction (judged by the 95% predictive interval containing zero; Fig. 4). The second lockdown was preceded by a declining trend in contact ratios while the contact during the lockdown remained relatively higher than the first or third lockdown (Fig. 3).

Impact of a more transmissible variant
Our model estimated that the Alpha strain was approximately 19% more transmissible than the original strain (τ 2 , 95% prediction intervals [16.0,21.8]). It is worthwhile noting that our model does not consider continuous inputs of infection into the island of Ireland, despite the high connectivity among the British Isles. Thus, our estimate of the Alpha transmissibility may be confounded by repeated introductions, for example, from England, where the Alpha strain was first detected. Nonetheless, our estimate is consistent with those from England [34].
To assess the impact of the Alpha strain, which arrived later and is more transmissible than the original strain, we compared the fitted model ( Fig. 5; orange) to a counterfactual simulation without the Alpha strain, in which we assumed the same estimated contact ratio ( Fig. 5; blue). We detected a statistically distinguishable impact of the Alpha strain on the cumulative hospital admissions by earlier January in both jurisdictions -approximately 3.5 months after the initial introduction (indicated by the 95% predictive interval excluding zero; Fig. 5). By the end of February 2021, we show that the Alpha strain was responsible for a 38 and 55% increase in cumulative hospitalisation, in NI and the ROI, respectively (Fig. 5). Our findings demonstrate that an introduction of a more transmissible variant -without necessarily being more severe -can cause a large impact on the health care burden. Fig. 4. Counterfactual analysis demonstrates the effect of lockdown timing on epidemiological outcomes. We examined counterfactual introductions of three lockdowns in Northern Ireland and the Republic of Ireland, assuming that they would have started seven days and 14 days earlier. The percentage difference in cumulative hospital admissions between the counter-factual and factual scenarios is shown. The black line and grey band indicates the median and 95% predictive interval, respectively. The predictive interval signifies 95% of simulated outcomes generated based on samples from the posterior distribution of model parameters and measurement errors.

Conclusion
We developed a multi-strain model of SARS-CoV-2 and estimated timedependent human contact over the first 12 months of the pandemic on the island of Ireland. Unlike many earlier COVID-19 modelling studies that estimate the effective reproductive number of a single strain, our model explicitly incorporates multiple viral strains and focus on estimating contact ratios. An important difference between the contact ratio and the effective reproductive number is that the former is unaffected by changes in virus transmissibility, which is modelled independently. As such, our approach separates the effect of human behaviour from that of the difference in transmissibilities between multiple, co-circulating strains.
Examining the longitudinal patterns and geographical differences in the estimated contact ratios allowed us to identify corresponding policies and events. In addition, we leveraged estimated parameters to conduct counterfactual analyses, in which we examined the role of lockdown timing and a novel variant on cumulative hospitalisation. In a companion paper, we extended the application of the estimated contact ratios to causal inference [49]. Specifically, we used mobility and maskwearing data to independently predict the contact ratios estimated from our epidemiological model described in the current paper and subsequently compared observed hospitalisations with predicted hospitalisations under a counterfactual mask-wearing scenario.
We presented a generic, epidemic model parameterised for SARS-CoV-2 to fit longitudinal hospitalisation data, one of the most reliable and available data types [6]. Of most relevance to COVID-19 at the time of publication, our model lacks human age structure and vaccination: these omissions give rise to certain limitations. For example, hospitalisation risks increase with age while older individuals adjust their behaviour differently from young counterparts [50]. Thus ignoring the age structure may bias our estimate of human contact estimated from hospitalisation data. In addition, the lack of vaccination and associated immunity in our model restricted our scope to the first 12 months of the COVID-19 pandemic. Technically, our model can be extended modularly to relax these assumptions about age structure and vaccination. However, these extensions were outside the scope of this study due to challenges in parameterising these processes reliably. For instance, the output of age-structured models are highly sensitive to assumptions of age-specific contact patterns [51], which likely changed during the epidemic, yet empirical data for time-dependent contact matrices are scarcely available. Behavioural adjustment in response to the pandemic is further complicated by the interaction between age-and sex-specific effects [52]. Furthermore, it is difficult to track and parameterise the state of immunity generated by natural infections from multiple viral strains and multiple vaccine doses using compartmental models.
Finally, our work contributes to the growing COVID-19 modelling literature by providing a transparent Bayesian workflow for fitting a multi-strain epidemic model to longitudinal epidemiological data, which may be readily adapted to modelling SARS-CoV-2 in other jurisdictions and other infectious diseases.

Fig. 5.
Counterfactual analysis shows the extent to which the Alpha strain elevated the burden of hospitalisation. To compute the impact of the Alpha strain, the counterfactual simulation (blue) assumes that the Alpha strain never invaded either jurisdiction. The crosses indicate data, and the coloured bands correspond to 95% predictive intervals of the fitted model (orange) and counterfactual scenario (blue), respectively (left panels). The percentage difference in cumulative hospital admissions between the two scenarios is shown. The black line and grey band indicates the median and 95% predictive interval, respectively (right panels). The predictive interval signifies 95% of simulated outcomes generated based on samples from the posterior distribution of model parameters and measurement errors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Declaration of Competing Interest
None.