A real-time regional model for COVID-19: Probabilistic situational awareness and forecasting

The COVID-19 pandemic is challenging nations with devastating health and economic consequences. The spread of the disease has revealed major geographical heterogeneity because of regionally varying individual behaviour and mobility patterns, unequal meteorological conditions, diverse viral variants, and locally implemented non-pharmaceutical interventions and vaccination roll-out. To support national and regional authorities in surveilling and controlling the pandemic in real-time as it unfolds, we here develop a new regional mathematical and statistical model. The model, which has been in use in Norway during the first two years of the pandemic, is informed by real-time mobility estimates from mobile phone data and laboratory-confirmed case and hospitalisation incidence. To estimate regional and time-varying transmissibility, case detection probabilities, and missed imported cases, we developed a novel sequential Approximate Bayesian Computation method allowing inference in useful time, despite the high parametric dimension. We test our approach on Norway and find that three-week-ahead predictions are precise and well-calibrated, enabling policy-relevant situational awareness at a local scale. By comparing the reproduction numbers before and after lockdowns, we identify spatially heterogeneous patterns in their effect on the transmissibility, with a stronger effect in the most populated regions compared to the national reduction estimated to be 85% (95% CI 78%-89%). Our approach is the first regional changepoint stochastic metapopulation model capable of real time spatially refined surveillance and forecasting during emergencies.

Conserving the number of people in a location requires a balance in people not moving, people travelling out from the location and people travelling to it. In terms of the mobility matrices, this means the following: For a given location i and time t, all incoming transitions are summed up (column-wise), and this must equal the sum of outgoing transitions at time t + 1 (row-wise). Assume now that we want to simulate ahead in time, using only one mobility matrix for all times t. When we use the same matrix all the time, then Eq (S1) can be rewritten as Eq (S2) (where we drop t): For all sites i k:k̸ =i which means that the matrix row-sums are equal to the column-sums. A symmetric matrix trivially fulfils this requirement. However, this is not the only possibility, and skewed mobility matrices are important due to commuting patterns.

S1
We find the matrix Z that preserves the population of each region by solving the following Quadratic Programming problem: subject to the linear constraints ∀i k:k̸ =i For the predictions in this paper, we used regularised mobility data from the dates in Table A

S1.2 Hospitalisation risk and length of stay
To calculate the hospitalisation probability conditional on infection per age group, we started with the estimated hospitalisation probabilities in [1]. Note that the compartmental model did not have age compartments, as they would make the model computationally intractable. On the other hand, the probability of requiring hospitalisation was highly dependent on age. We corrected the hospital admission probabilities to account for age-dependencies in transmission. The disease burden did not just follow the age demography, and we needed to quantify the proportion of infected individuals per age group. To this end, we used the age distribution of the positive cases as a proxy for the age distribution of COVID-19 infections. Test data a commonly available data source. For Norway, we obtained test data from the MSIS registry. To account for a changing age profile of infected throughout S2 the pandemic, we used monthly estimates from May 1, 2020 in the national model, reported in Table B. In the regional model we used daily estimates. Before May 1, 2020, we do not use the test data. Hence, because the demographic age profile in each county was available, we computed the county-wise probabilities of hospital admission as a mixture of the probabilities per age group in [1] prior to May 1, 2020, each multiplied with the proportion of the age group in the counties.
Compared to the population under study in [1], Norway has had a practice of not admitting elderly with COVID-19 living in nursing homes to hospitals. Therefore, we reduced the hospitalisation risk for the 80+ age group by 15% which corresponds to the percentage of non-nursing home residents . The total hospitalisation, h i,t for time, t and region i is given by where p ijt is the fraction of positive tests in age group j in region i at time t and S j are the baseline probabilities from [1] corrected by the percentage of non-nursing home residents. The time-scale of t used was monthly in the national model and daily in the regional model (see section below). When we predicted future hospitalisations in this paper, we did not use the future hospitalisation risks as it would require the use of impending positive test data that would not be available in a real-time setting. Instead, we used the latest available risk estimates for the relevant dates.

S1.2.1 Regional hospitalisation risks
The age distribution of the positive cases varied between counties and changed over time. We estimated the daily hospitalisation proportion of positive tests in each age group, county-wise and nationally, over a moving window over the last 28 days. We used the observed proportion of hospitalised cases in each age group. We added a small number in both the number of hospitalisations and the number of positive tests in each age group to avoid dividing by zero. Specifically, a pseudo count of 0.001 was added to the number of individuals hospitalised in each age group before dividing by the total number of individuals testing positive plus the number of age groups times 0.001. The country-specific proportions were adjusted using national data in those regions when few tested positive in a given period. This was done by calculating a weighted fraction: the county-specific S3  [2,3]. The two first detected cases of the alpha variant imported to Norway occurred on December 27, 2020 [4]. Since then, Norway experienced an increasing proportion of cases attributed to the alpha variant, which eventually became predominant. We adjusted for a gradually increasing hospitalisation risk by multiplying by a logistic increase in hospitalisation risk between January 1, 2021 and March 31, 2021, assuming a total increase of 60%, resulting in the following logistic function B 117 (t) = 1000 exp(−0.073t) + 1.6 · 10 exp(0.073t) 1000 exp(−0.073t) + 10 exp(0.073t) (S6) whose slope was estimated from the test data, and t denotes number of days from January 1st 2021.
The equation S7 shows how the regional hospitalisation probabilities are calculated from May 1st 2020 taking into account the take over curve of the B 117 variant showed in S6. In the national model, we assumed an increase of 9% in January, 50% in February and 60% from March 2021 corresponding to the monthly averages of the logistic increase.

S1.2.3 Length of stay in hospital and time to hospital admission
To compute hospital incidence, we needed to know the time from symptom onset to hospital admission. Moreover, the predictions of hospital prevalence and the number of occupied intensive care treatment beds have been of interest, requiring estimates of the length of stay (LOS) in hospital for COVID-19 patients and the proportion needing intensive care treatment. The LOS in hospital for COVID-19 patients may change over time. We therefore allowed changepoints in the hospitalisation parameters. In addition, there may be regional differences, and we thus also included the possibility for separate parameters for certain counties.

S5
These periods will here be called time-space periods. The parameters for the timespace periods were estimated separately.
The aim was to estimate distributions for the time from symptom onset to hospitalisation, the LOS for patients with and without ventilator treatment, and for the proportion of hospitalised cases needing ventilator treatment. The LOS for patients needing ventilator treatment was divided into three periods: before, during and after ventilator treatment to obtain estimates of the number of required ventilators per day. For patients with more than one COVID-19 hospitalisation episode, the LOS was set equal to the sum of all admission periods.
We used a proportional hazard Cox model to estimate the parameters for each time-space period, assuming right censoring in time. The survival curves were estimated with time-space periods as an explanatory variable, using the survfit function in the R-package survival. Specifically, one minus the cumulative distribution function of the negative binomial distributions was fitted to each survival curve by minimising the sum of squared errors, using the function optim in R.
To estimate the hospitalisation parameters we used individual-level data from the Norwegian Surveillance System for Communicable Diseases (MSIS) and the Norwegian Intensive Care and Pandemic Registry (NIPaR) [5] as explained in the main text. By April 12, 2021, NIPaR included 4.248 hospitalised individuals, of whom 735 needed mechanical respirator treatment, while MSIS contained 103.461 positive test results. We include three changepoints, dividing the data into four periods -before June 2020, June 2020-December 2020, January 2021-February 2021, and March 2020-April 12, 2021.
Except for the county of Oslo, there were no significant differences in the regional COVID-19 hospitalisation parameters. The hospitals in Oslo increased the use of mechanical ventilators from January 2021. Therefore, we obtained separate estimates for Oslo in the last two time periods. When estimating the time spent under ventilator treatment, we merged the periods January 2021-February 2021 and March 2021-April 12, 2021, to reduce the level of censoring. In a separate project, we are trying to find the causes of these changes, see also [6]. Figs A, B, C, D and E show for each time-space period a histogram of time to event together with an estimated negative binomial distribution fitting the data in red. The parameters of the negative binomial distribution (mu and size) together with the number of patients (N) in the data are shown in the upper right corner of each panel.

S1.2.4 Time from symptom onset to hospitalisation
When estimating time from symptom onset to the first hospitalisation, patients with symptom onset date equal to the COVID-19 test date were disregarded. In Norway, this registration practice is common in situations where symptom onset is unknown. We assumed that these cases were missing at random and hence did not introduce bias. All patients included in this analysis were hospitalised, so no right censoring was needed for this calculation. Fig A shows the real data and the fitted distribution in red.

S1.2.5 Time from hospitalisation to discharge
The time from hospitalisation to discharge was estimated separately for two groups, those who required ventilator treatment and those who did not. Right censoring occurred in the last time-space period from March 2021; 15 and 65 patients from Oslo and the rest of Norway, respectively, were still in hospital at the end of that period, out of a total of 476 patients from Oslo and 960 from the rest of Norway. The number of patients included in each time-space period is shown in Fig  B together with the real data and the fitted distribution in red.
For those needing mechanical ventilator treatment, we considered three time intervals: time from admission to start of ventilator treatment (Fig C), time on ventilator ( Fig D) and time to discharge after ventilator treatment (Fig E). Both the estimation of the distribution of the time on ventilator and time after ventilator treatment included censoring events in the last period from January 2021. 13 patients from Oslo and 23 patients from the rest of Norway were censored when estimating time on ventilator, out of a total of 132 patients from Oslo and 207 patients from the rest of Norway. For the time after ventilator treatment, the censored patients were 17 and 28 for the same periods, out of a total of 120 for Oslo and 170 for the rest of Norway.
The proportion of patients needing ventilator treatment has changed over time.

S1.3 Parameter assumptions
A lot is uncertain about the natural history of SARS-CoV-2. Hence, there are many unknown parameters in the model that our method estimates (or calibrates) from data. We calibrate the amplification factor and the transmission parameters β i t , which are directly related to the reproductive numbers. We also calibrate the parameters π 0 and π 1 in the model for the detection of positive tests, and the delay d between entering the infectious symptomatic or asymptomatic class and a positive test. As described above, we estimate the hospitalisation parameters from individual-level data on hospitalised cases.
The remaining parameters related to the natural history of COVID-19 are assigned fixed values or distributions based on values found in the literature. All the parameter choices and values are provided in Table C.

S1.4 Regional baseline model
To compare the performance of the regional model with a regional baseline or null model, we implemented a simple baseline model. The model uses the last 21 S8 days of the data (hospitalisation or confirmed cases) and fits a negative binomial regression model with only an intercept term. The regression model is then used to make predictions for the next 21 days as the in the other models. To ensure some variation in the data if there were no hospitalisations during the training period, we added one hospitalisation on the first day of the training data.

S2.1 Estimated national reproduction numbers
In addition to Table 2 in the main text, we also provide a visualisation of the estimated national reproduction numbers in Fig F.

S2.2 Validation results
The validation scores are provided in Table D. Table D: Validation results comparing the regional, national and regional baseline model. National Regional Regional baseline Week

S2.3 Fit to test and hospitalisation data
When estimating the transmissibility, we choose to use both the hospitalisation data from the beginning of the epidemic and the test data from May 2020, as both contain information.

S2.3.1 Fit to data up to April 2021
The fit to the regional hospitalisation data until April 2021, along with three-weeks ahead predictions, is provided in Figs G-H, for the model allowing region-specific transmissibilities. The fit to the test data is provided in Figs I-J. Note that the test data are more volatile than the hospitalisation incidence due to flare-ups of local outbreaks that become controlled. For this reason, test data are not well captured by the constant step-function in the transmissibility unless the dates for the changepoints are tailored to the dates of the outbreaks (see for example Agder). However, the test data do contain information on the transmissibility, in particular for the more recent days. We note that for the first days of prediction, our hospital incidence predictions seem to cover the real data that are later observed, but that we start overpredicting over time. This is due to local restrictions implemented in Oslo on March 17, 2021, and national restrictions implemented on March 25, 2021, which have reduction effects that we have not captured with our most recently estimated effective reproduction number, which is active from March 2, 2021. This is particularly clear in the fit to the largest region in Norway, Viken, and in the national aggregate.
Note that we do not capture all the uncertainty in the test data. This is likely because we calibrate to a 7-days-backwards moving average of the daily test data. Moreover, the binomial assumption could be too simple, and a beta-binomial distribution would allow more variance, at the cost of an additional parameter.

S2.3.2 Fit to data up to March 2021
We include also the fit to the test data and hospitalisation data for the calibration up to March 2021. Here, no changepoints later than February 8, 2021 were used. The fit to the regional hospitalisation data until March 1, along with three-weeks ahead predictions, is provided in Figs K-L, for the model allowing region-specific transmissibilities. The fit to the test data is provided in Figs M-N. The predictions of future hospitalisations are better here than for April 1, 2021, as the intervention measures implemented later in March take more time before the effect shows up in the data. We note however that we are underestimating the hospitalisation a bit, in particular in Viken.
We repeat the calibration up to March 1, 2021 using only hospitalisation data. The fit to the hospital incidence is provided in Figs O-P. If we consider Viken, we note that the hospitalisation predictions are better in this setting when only calibrating to the hospitalisation data, than when we calibrate to both the hospitalisa-S10 tion data and the test data, when we underpredicted the hospitalisation incidence (Figs K-L). However, the variance is also larger.

S2.3.3 Fit to data up to November 1, 2020
We include also the fit to the test data and hospitalisation data for the calibration up to November 1, 2020. Here, no changepoints later than October 1, 2020 were used. The fit to the regional hospitalisation data until November 1, 2020, along with three-weeks ahead predictions, is provided in Figs Q-R, for the model allowing region-specific transmissibilities. The fit to the test data is provided in Figs S-T. We note that for the future hospitalisations, we are underestimating the cases in Oslo.

S2.3.4 Fit to data up to October 1, 2020
We include also the fit to the test data and hospitalisation data for the calibration up to October 1, 2020. Here, no changepoints later than September 1, 2020 were used. The fit to the regional hospitalisation data until October 1, 2020, along with three-weeks ahead predictions, is provided in Figs U-V, for the model allowing region-specific transmissibilities. The fit to the test data is provided in Figs W-X. We note the low level of hospitalisation incidence in Norway in October 2020.

S2.3.5 Fit to data up to September 1, 2020
We include also the fit to the test data and hospitalisation data for the calibration up to September 1, 2020. Here, no changepoints later than August 1, 2020 were used. The fit to the regional hospitalisation data until September 1, 2020, along with three-weeks ahead predictions, is provided in Figs Y-Z, for the model allowing region-specific transmissibilities. The fit to the test data is provided in Figs AA-AB. In September 2020, there were few COVID-19 cases in Norway.

S2.3.6 Fit to data up to August 1, 2020
We include also the fit to the test data and hospitalisation data for the regional calibration up to August 1, 2020. This is the calibration of the past which is considered fixed in the later regional calibrations. Here, no changepoints later than July 1, 2020 were used. The fit to the regional hospitalisation data until August 1, 2020 is provided in Figs AC-AD, for the model allowing region-specific transmissibilities. The fit to the test data is provided in Figs AE-AF. We note S11

S2.5 Computational time and convergence
In particular for the regional calibration, the parameter space is very high-dimensional and correlated in both time and space, which makes the calibration exercise difficult. The calibration took one week running on 1400 cores. The run-time for a single sample is in seconds, what mainly takes time is to find parameters which improve the error, not the disease spread simulation itself. Fig AH shows a typical example of the evolution of the estimated distributions of the parameters throughout the calibration rounds.

S2.6 Policy compliance and reduced mobility
The mobility data are indicators for compliance with mobility reducing policies, like travel restrictions and working from home. Fig AI shows the daily number of trips out from each of Norway's counties. There was a clear reduction in mobility during the March 2020 lockdown of approximately 50%. Thereafter, mobility increased steadily until the summer. In July, geographic differences in mobility corresponded to summer vacations. Later in 2020, Norway experienced high, regionally varying mobility, with peaks around vacations. Late 2020 and early 2021, mobility decreased again during the second and third waves. The systematic mobility drop on weekends indicated commuting to work even though teleworking was recommended. We compared mobility to March 1, 2020, as we do not have mobility data before the pandemic, which also means that we cannot compare general seasonalities in mobility volume.
For some counties (e.g., Troms og Finnmark), the mobility recovered to the levels prior to the first lockdown, while this was not the case for others (Oslo and Viken).

S2.7 Sensitivity to mobility implementation
We performed a sensitivity analysis to ensure robustness to the choice of mobility rule for moving people between regions in the model. To allow for comparison of the spread of the pandemic between different regions, we utilised a national model for this analysis, assuming common fixed reproductive values for each region.
We ran the model from February 17, 2020 to October 18, 2021 with two different R 0 values of 1.1 and 1.5 in order to understand the importance of the mobility rule with different infectivity levels in the pandemic. We use the actual imported cases in the simulations. To understand the importance of these mobility rules in a S14 more interpretable scenario of limited foreign seeding, we in addition ran simulations with each rule starting with 100 cases in Oslo on day 1 (February 17, 2020) and no seeding thereafter.
We ran 100 simulations of the model assuming different mobility rules, when considering the problem of how to move individuals from location A to location B. We denote by "visitors" the individuals who are present in B and have A as their home location. We denote by "host" the individuals who are present in location B and have B as their home location. We denote by "others" the individuals who are present in location B and have other locations than A or B as their home location. The rule applied in the paper will therefore hereafter be referred to as V isitors − Host−Others (V HO), as we first send the individuals here defined as "visitors", then the individuals here defined as "host", and lastly the individuals here defined as "others". If there are not enough individuals present in location B, we add susceptibles to the host population of B. As a control experiment, we also ran another set of 100 simulations for the rule applied in the paper, hereafter referred to as V isitors−Host−Others2 (V HO2). Hence, VHO2 is a repetition of VHO to assess the variability between two runs with the same rule. We additionally explored four other mobility rules defined as follows: 1. Random (R): We randomly sample from all people in a location, without distinguishing between visitors and host population. Then, we add extra people if necessary.
2. V isitors − Random (V R): We send back visitors from location B first, and then randomly sample from the remaining individuals who are present in A and add extra people if necessary.
3. V isitors − Others − Host (V OH): We send back visitors from location B first, and then select from the other visitors to the location. Finally, we select from the host population and add extra people if necessary.
4. V HORegularised (Reg): We used the VHO mobility rule applied in the paper to move people between two locations, but we used regularised matrices instead of the original origin-destination matrices. This ensures that we do not have to add extra people into the model at any point.
To better understand the effect of regularisation of the mobility matrices, we additionally quantified the absolute daily difference between inter-county mobility in the regularised matrices and the original mobility matrices from February 17, 2020 to September 27, 2021 (Fig AM). We observe that the differences are generally higher during the summer and autumn seasons than the winter and spring seasons.
We determined the date and size of the peak, i.e. the maximum number of S15 infected cases on a national and county level for each of the simulations. We define infected cases as the number of infectious, infectious asymptomatic, and pre-symptomatic infectious cases who are residents in the county of interest and are thus considered a part of the host population of that county. We then plotted the mean of the peak size for the 100 simulations nationally, (Fig AJ), and countywise ( Fig AK).We also computed the mean number of days until 1% of each county's resident population was infected (Fig AL). We only conducted regional analyses for a scenario with limited foreign seeding, as it is easier to interpret the effect of the mobility rules without additional seeding throughout the pandemic.
Using the metrics of the national and county-specific peak sizes and time to 1% infected, we observe very few statistically significant differences between the rules. However, the random rule (R) has a higher national peak that is almost statistically significant, which is expected, as infections tend to spread more with random movement of people [8]. Though, we do not find this behaviour in the county-specific analysis. In Oslo county, the peak is lower when moving visitors and hosts randomly (VR) compared to those obtained with the other mobility rules. This finding may be explained by the strong asymmetry in the matrices from this particular county. With the regularised matrices (Reg), the peaks tend to be lower than with the other rules, though not significantly. One reason for this observation is that no extra people are added to the model, which may give rise to fewer infected people at the peaks. Given that the differences between the results from the different mobility rules were almost always statistically insignificant, the sensitivity analysis indicates that our simulations are robust to the choice of mobility rule.

S2.8 Effect of mobility on exponential growth rates
We illustrate the effect of mobility on the exponential growth by simulating a set of scenarios where we increase the reproduction number in Oslo county while setting the reproduction number in every other county equal to one. We simulate these scenarios with and without mobility included in the model. In Fig AN we see that when the difference in reproduction number is large between counties and with substantial mobility, there is a significant decrease in the exponential growth rate for a given reproduction number when we include mobility. When there are large differences between reproduction numbers, one could easily underestimate the reproduction number in the county with higher transmission by 0.5. The effect is driven by a large number of infected persons travelling to other counties, while S16 few infected people travel back to Oslo since the prevalence of infected are much lower in other counties.

S2.9 Comparison of regional reproduction numbers with alternative approaches
A much simpler approach to estimating the (regional) reproduction numbers based on case data is the commonly used method implemented in the R-package EpiNow2 [9]. We run EpiNow2 with the same changepoints as in the regional model. We use an incubation period based on [10] and a generation time with mean 7 days standard deviation of 3.5 days matching the generation time in the regional model. Using a method like EpiNow2, we lose a lot of details, like for example changes in imported cases over time, the connection between the different regions over time due to mobility, and changes in testing criteria and capacity. It is also not straightforward to include multiple case data sources at the same time, and hence we show the results from EpiNow2 when calibrated to the number of positive tests only. The resulting reproduction numbers for Oslo and Viken are provided in Fig AO and AP, respectively. In general, the overall patterns are the same. However, we see that the discrepancy is large for the early part of the pandemic. As discussed earlier, this is likely because the test data were not reliable in this period due to frequent changes in testing criteria and capacity. For Oslo, the estimates from our model framework are in general higher, as we have an outflow of infectious individuals from Oslo due to mobility (see also the previous section). The opposite effect is found for Viken. For both Viken and Oslo, we note that the estimates for early autumn 2020 from EpiNow2 are higher than our estimates, which is likely due to a large number of imported cases coinciding with border reopening, which is not taken into account in EpiNow2.

S2.10 Sensitivity analysis of the probability of asymptomatic infection
In order to evaluate whether the estimated reproduction numbers were robust to the probability of developing asymptomatic infection (p a ), which was set to 0.4 in the main analyses, we re-calibrated the model with a national reproduction number, with data until April 1, 2021 (like in section "The National picture" in the main text) with p a = 0.3 and p a = 0.5. Fig AQ presents a comparison of the S17

S2.11 Verification test for Split-SMC-ABC
In order to study the performance of the Split-SMC-ABC calibration method suggested in this paper, we have performed one simple test with synthetic data, generated with a certain set of known, true underlying parameters. We generated synthetic data (Fig AR) from the beginning of the pandemic up until October 1, 2020 assuming that the reproductive number was the same in all regions except Oslo. We used four changepoints, resulting in a total of ten reproductive numbers: five for Oslo and five for the other counties (see Table H for more details). The remaining parameters of the model (π 0 ,π 1 , delay between infectious and testing) and the amplification factor, were assumed fixed.
In the Split-SMC-ABC, we applied a split on July 25, 2020, running the first calibration using data 15 days into the future, as in the main analysis. The assumed prior distributions were the same as those used in the main analysis.
As shown in Fig AS and Fig AT, both the presplit and postsplit stages converged asymptotically to the ground truth values for most parameters. One exception is R 1 in Oslo, where the method seems to converge, but to a number which is somewhat higher than the ground truth. Note, however, that we are calibrating to one single realisation of the stochastic disease spread, and hence it might be that the underlying truth is not the best fitting parameter. S18 S3 Implementation     Observed (orange dots) and simulated hospitalisation incidence for each region, along with 95%, 90%, 75%, and 50% credibility intervals. The regions are ordered according to population, with the capital Oslo as second. The incidence for the next three weeks not used in the calibration is provided in red. Observed (orange dots) and simulated hospitalisation incidence for each region, along with 95%, 90%, 75%, and 50% credibility intervals. The regions are ordered according to population. In the lower right panel, we have the national fit for the model assuming the same transmissibility in all counties. The incidence for the next three weeks not used in the calibration is provided in red.   Observed (orange dots) and simulated hospitalisation incidence for each region, along with 95%, 90%, 75%, and 50% credibility intervals. The regions are ordered according to population, with the capital Oslo as second. The incidence for the next three weeks not used in the calibration is provided in red. Observed (orange dots) and simulated hospitalisation incidence for each region, along with 95%, 90%, 75%, and 50% credibility intervals. The regions are ordered according to population. In the lower right panel, we have the national fit for the model assuming the same transmissibility in all counties. The incidence for the next three weeks not used in the calibration is provided in red. Observed (orange dots) and simulated hospitalisation incidence for each region, along with 95%, 90%, 75%, and 50% credibility intervals. The regions are ordered according to population, with the capital Oslo as second. The incidence for the next three weeks not used in the calibration is provided in red.