Forecasting local hospital bed demand for COVID-19 using on-request simulations

Accurate forecasting of hospital bed demand is crucial during infectious disease epidemics to avoid overwhelming healthcare facilities. To address this, we developed an intuitive online tool for individual hospitals to forecast COVID-19 bed demand. The tool utilizes local data, including incidence, vaccination, and bed occupancy data, at customizable geographical resolutions. Users can specify their hospital’s catchment area and adjust the initial number of COVID-19 occupied beds. We assessed the model’s performance by forecasting ICU bed occupancy for several university hospitals and regions in Germany. The model achieves optimal results when the selected catchment area aligns with the hospital’s local catchment. While expanding the catchment area reduces accuracy, it improves precision. However, forecasting performance diminishes during epidemic turning points. Incorporating variants of concern slightly decreases precision around turning points but does not significantly impact overall bed occupancy results. Our study highlights the significance of using local data for epidemic forecasts. Forecasts based on the hospital’s specific catchment area outperform those relying on national or state-level data, striking a better balance between accuracy and precision. These hospital-specific bed demand forecasts offer valuable insights for hospital planning, such as adjusting elective surgeries to create additional bed capacity promptly.


Introduction
The speed of community transmission of virulent pathogens during epidemics or pandemics of infectious diseases has the potential to overwhelm healthcare systems [1].Hospital beds are usually in limited supply, because most hospitals run at nearfull capacity under normal circumstances.A sudden increase in patients needing hospital care might thus not be met with the required available beds if the growth in demand is unexpected.Reliable forecasting of bed demand during epidemics is therefore essential for hospitals to free up beds in time.
However, epidemic trajectories can differ considerably regionally [2], and because patients are usually admitted to hospitals in their own region [3], the bed demand can widely vary between hospitals as well.Consequently, nationwide forecasts of bed demand may not reflect the individual hospital's situation, and may fail to predict overload problems.We, therefore, developed a bed demand forecasting model [4] based on local epidemic trajectories in the hospital's individual catchment area and its current bed occupancy.
To be widely applicable, this model needs to be able to produce bespoke forecasts for each hospital in a country for which basic levels of data are available.This capability requires a flexible forecasting system because of the vast number of potential combinations of parameter choices for all individual hospitals, in particular concerning the possible choices for catchment areas per hospital, the number of currently occupied beds, and patients' lengths of stay on various wards.Such a variety of parameters enables the extreme forecast flexibility required to produce ad-hoc reports and forecasts for all situations of specific hospitals.
We made the use of this model accessible for non-technical users by creating an online interactive dashboard.This platform allows individual users (i.e., hospital managers) to enter the data and parameters related to their specific context and requirements and produce ad-hoc bed occupancy estimates and predictions.End user-specific forecasting of bed demand is a unique concept, as on-request forecasting is often avoided because of the computational challenges of multiple concurrent user modelling platforms.
Here we show how some of these computational challenges can be solved.Careful streamlining of the model and its technical implementation can help produce epidemic forecasts in a reasonable time, allowing users to explore the potential future epidemic trajectory without being restricted to the developer's viewpoint.
We will first discuss the forecast model itself, its data and parameter requirements, as well as the different modules it contains.Each module is designed to guide the user through the steps and assumptions needed to produce the forecast.Then we will discuss the technical implementation of the online dashboard and the technologies that ensure usability for multiple concurrent users.Finally, we will present performance statistics of the dashboard collected between December 2021 and March 2022.

General model structure
The complete forecast model consists of 5 main modules: 1) Data loading and nowcasting, 2) Reproduction number estimates, 3) Vaccination coverage forecasting, 4) Incidence model, and 5) Care path model.The calculations done in the first two modules rely only on the loaded data and user input, while modules 3-5 also build further on the results of each previous module (see figure 1).This hierarchical structure implies that any change to parameters in a module updates the outputs in all dependent modules.However, by guiding the user through the modules in sequential order, the number of required computations can be drastically reduced.

Data structure and definitions
The model needs data on infection incidence, vaccine doses, bed occupancy in general wards and intensive care units, and a simulation window for which a forecast is required.Although we developed this model for hospitals in Germany, based on data provided by the Robert Koch Institute [5] and the Deutsche Interdisziplinäre Vereinigung für Intensiv-und Notfallmedizin (DIVI) [6], the model can be applied to any country that has the following data available (preferably if in remote-accessible, machine-readable format): 1. incidence, I g (t) 2. administered vaccinations, V d,g (t) 3. bed occupancy on general wards B GW,g (t) and ICU B ICU,g (t) where t denotes the time in days since a reference date before the start of the pandemic, g denotes a geographical subdivision (in the case of Germany, this is a district, i.e., Landkreis or Stadtkreis), and d denotes the dose of vaccine administered (1st, 2nd, or 3rd/booster).Incidence and vaccination data are needed for each day of the period or reference, while merely the value at the start date of the forecast is needed for bed occupancy.
The user defines the catchment area of interest, selecting geographical areas from the available list (K all ), into the catchment set (K u ).After this, the input data is summed over the selected areas: The start of the simulation is denoted as T S , for which the incidence, number of vaccinations, and bed occupancy are taken from the (last data-point of the) data, the first forecasted day is thus T S + 1.The length of the forecast is defined as L; the last forecasted date is therefore given as T E = T S + L. The size of any class of individuals in the model is given in absolute numbers of individuals; thus S(0) = N , with N being the population size.
Because many of the parameters in the dashboard's models are continuously distributed (as Exponential, Gamma, or Weibull), and the model uses discrete time steps, we need to discretise the parameter distributions.For each continuous distribution, with probability density function (PDF) f c (t), the discretised PDF is defined as User input for all distributions is given in terms of the mean µ and standard deviation σ, and distribution parameters are then determined by moment matching.For a gamma distribution, this means that f (x, α, β) = f (x, α = µ 2 /σ 2 , β = µ/σ 2 ), for a Weibull distribution, f (x, k, λ) = f (x, k = ( σ µ ) −1.086 , λ = µ/Γ(1 + 1/k)), and for an exponential distribution f (x, λ) = f (x, 1/µ).

Vaccination model
The vaccination model consists of two parts.In the first, we forecast the number of people vaccinated over the length of the forecast for each of the doses (1st, 2nd, and 3rd/booster dose), while in the second part we convert these administered vaccinations into population-level protection against transmission.
The first doses (V 1 (t)) are assumed to be administered at the same rate as during the last observed week: The second dose is assumed to be administered at a fixed time delay after the first dose.This time delay ∆T V is extracted from the data as the maximum time for which holds.The second doses are then a direct reflection of the number of first doses ∆T V days ago, The daily administered third/booster doses are assumed to be a continuation of the mean daily third doses over the previous week: This is in line with the forecasting of the first doses, with the exception that the cumulative number of third doses is not allowed to exceed the cumulative number of second doses a certain delay (∆T B ) ago, such that has to hold.

Population protection
Each dose is assumed to have a time-dependent additive effect on the populationwide protection against infection.At the individual level, G d (τ ) denotes the proportion of individuals no longer susceptible to infection τ days after the administration of dose d.This function thus serves to simulate the delay between vaccine administration and maximum vaccine protection at the individual level.We assume no waning of immunity, such that G d (∞) = E d , with E d the vaccine effectiveness against transmission elicited by doses d.
The additive effect on vaccine effectiveness of each additional dose is defined as consequently, the additive increase in an individual's protection due to an extra dose can be written as: G d (τ ) is input as the cumulative distribution function (CDF) of a normal distribution with mean µ(G 1 ) = 15, µ(G 2 ) = 15, µ(G 3 ) = 7, and standard deviation σ(G 1 ) = 3.8, σ(G 2 ) = 6.5, σ(G 3 ) = 3.8 as default values.
The population protection at time t is then given by, which is used in both the R forecasting step and the incidence model.

R estimation and forecasting
In order to inform the transmission process part of the incidence model, we need to forecast the reproduction number of the pathogen over the model timeframe.This forecast is based on the observed time-varying effective reproduction number (R e (t)), and implemented either as the static continuation of the most recent observation, or forecasted using an ETS (Error, Trend, Seasonal) model [9] based on the last 100 days of observations, depending on the user's preference.Furthermore, we provide the option to incorporate the effect of a Variant of Concern (VoC) on the development of the R 0 (t).In all cases, we use estimated values of R e (t) from the available data up to the simulation start date T S , R e (t ≤ T S ), and forecasted values after the simulation start date (R e (t > T S )).
We estimate the time-varying effective reproduction number R e (t) using the Cori et al.R estimation method [7] as implemented in the EpiEstim R package (version 2.2-4), with a given non-parametric Serial Interval (SI).This SI distribution denotes the number of days between equal disease events (e.g., onset of symptoms) of directly connected cases.The default SI is assumed to be Gamma distributed with a mean of 5 days and a standard deviation of 4.9.These values were chosen to reflect estimates of mean SI from multiple studies [8] with a relatively high standard deviation.
By keeping track of the cumulative number of infected individuals, as well as the previously calculated population protection (G P (t)), we can calculate R 0 (t) from R e (t), given that R e is related to R 0 through the susceptible population:

ETS model
To produce a dynamic forecast of R 0 (t), we used an ETS model which underlies an exponential smoothing method [9].In general, exponential smoothing methods generate point forecasts of time series using weighted averages.
For an observed time series y 1 , y 2 , ..., y n the point forecast for h periods ahead is y t+h .When the forecast is based on all data up to time t, it is denoted as y * t+h|t .For a method with an additive trend, the point forecast results from the level of the series l and the growth b at time t: y * t+h|t = l t + hb t which gives the forecast equation.While l t and b t are based on the level and trend (slope) at t − 1 and are given by the smoothing equations: This approach is also known as Holt's linear method.Values for the initial states l 0 and b 0 as well as the smoothing parameters α and β * are estimated from the observed data, with 0 < α, β * < 1. Increasing α gives more weight to the more recent observations and less weight to observations that lie longer ago.In our case, α and β * were set to 0.25 and 0.15, respectively.
While the exponential smoothing methods generate point forecasts, the underlying statistical models, in addition, generate prediction intervals.To generate such an ETS forecast model, we need to specify the probability distribution for the residual at time t, e t .For additive errors, we assume the residuals are normally and independently distributed with mean 0 and variance σ 2 , notated as e t = ε t ∼ NID(0, σ 2 ).Expressing the error as with β = αβ * .These expressions are referred to as the state-space models as they include an equation for the observation and equations for the unobserved states (level, trend).Each observation y t in the ETS model therefore includes a component for the level (or smoothed value), an error component (E), a trend component (T), and seasonality (S).In our case, we use the ETS(A,A,N) model, which has additive errors, additive trend and no seasonality; We use a log(R 0 (t)) time series of a 100 days, i.e. log(R 0 (T S − 99...T S )), as input for the ETS model.The logarithm serves to forecast the relative changes in R 0 and avoids negative R 0 forecasts.Per model iteration, we pick a random quantile for the prediction interval and use this as the trajectory for R 0 (t > T S ).

Variant of Concern implementation
To model the impact of new pathogen variants emerging in the population [10], we need to consider how a variant of concern (VoC) has an advantage over the background viral population.We considered two possible mechanisms: 1) through higher transmissibility, increasing the base R 0 (t), and 2) by immune evasion, effectively increasing the size of the susceptible class (S(t)).
In case of increased transmissibility, we increase the basic reproduction number of the background viral population R 0,B (t) by the relative added fitness advantage Note that the previously calculated R 0 (t) is the average of the background and VoC reproduction numbers weighted by the proportion of cases caused by either (respectively ρ B (t) and ρ V (t)), We can calculate R 0,B (t), used as the general reference reproduction number in the model, from R 0 (t) as or we can directly calculate this from the observed R e (t) value .
The proportion of cases caused by the VoC (ρ V (t)) over time can be described with a sigmoidal function.A similar approach to estimating fitness advantages for new variants has been used by others [11].This is done by assuming both the background population and the VoC have a static growth rate unaltered by added immunity through new cases (i.e., the relation between their growth rates is fixed over time).We can calculate these growth rates (r) as a function of R 0 , using the mean serial interval (σ) [12] This means that the number of infections from both background variants and the VoC can be described as respectively The proportion VoC is .
By defining the total number of infected individuals as I T (0) = I V (0) + I B (0), and using .
This shows that we only need the relative added fitness advantage (A) in combination with a starting proportion of the variant (ρ V (t = t ρ )) at a certain reference date (t ρ ) to describe the proportion VoC over time.In the dashboard code, the reference date for the proportion of VoC does not overlap with the starting date of the incidence, vaccination, and bed occupancy data.In other words, with I(t I = 0) and ρ V (t ρ = 0), t I = t ρ , because ρ V (t I = 0) = 0, making a precise estimate of ρ V (t) impossible, or at least impractical.We, therefore, set t ρ to a later time point at which the ρ V is identifiable.
The added relative fitness advantage of the VoC can then be estimated from data about the proportion VoC among all isolates per week in the given catchment area, available through for instance genomic surveillance efforts.The user can then set the reference date, starting proportion, and added relative fitness advantage based on such analyses to be used in the forecast.

Immune evasion
In case the VoC has an advantage over the background variant through (partial) immune evasion, the number of individuals susceptible to the VoC (S V (t)) would be greater than for the background variant (S B (t)), giving it a higher R e (t) even if the R 0 (t) would remain the same.
We define immune evasion ( ) as the proportion of individuals immune against infection with the background variant that are not immune against infection with the VoC.Thus, the number of individuals immune to infection with the VoC scales linearly with the number of individuals immune to infection with the background variant, if very few or no infections with the VoC have occured yet.
To correctly assess the added fitness advantage at the VoC reference point (t ρ ), we need to calculate the increase in the number of individuals susceptible to infection with the VoC relative to those susceptible to the background variant.For simplicity, we assume that the cumulative number of infections with the VoC is negligible relative to the cumulative number of infections with the background variants at t ρ , that vaccination is aimed at the background variants, and immunity against the VoC is only determined by cross-immunity (1 − ) from immunity against the background variants.
Given that R e,B (t) = R 0 (t)S B (t)/N and R e,V (t) = R 0 (t)(1+A)S V (t)/N , the relative added advantage by immune evasion, given the pre-existing immunity against the background variant, is and the total added advantage is As A * is also the observed fitness advantage, the user will need to make an assumption about what part of this advantage is caused by immune evasion or by increased transmissibility on the reference date.Although the reference date (t ρ ) and the start of the simulation (T S ) are not necessarily the same, since the user can define t ρ independently of T S , from here on, we only consider the case where T S = t ρ , for clarity and brevity.

Incidence model
The incidence forecast model is based on a stochastic implementation of a Susceptible-Infected-Recovered (SIR) model.The model describes the number of new cases as a function of the population still susceptible (S(t)) to infection and of infection pressure (P (t)) exerted by all those currently infectious, and the total population size (N ).The size of the components is given as absolute numbers of individuals; thus S(0) = N .The infection process is governed by the forecasted basic reproduction number R 0 (t ≥ T S ) computed in the preceding modules of the framework.
If no VoC is defined, the model describes the dynamics of a single variant (the background variant).Conversely, the model converts into a two-strain model if a VoC is defined, keeping track of the infected individuals for either variant as well as the total number of individuals susceptible to either variant.

Single variant model
The size of the susceptible class is determined by both the cumulative number of infected individuals and the population level protection through vaccination, as defined before, The infection pressure exerted on S(t) is then given as the weighted number of infected individuals in the preceding serial interval, These are thus all past cases that are able to cause new cases on the current day t, weighted by the probability that they do so i days after they were reported themselves.Given the current R 0 (t) and P (t), the expected mean total number of newly infected individuals on day t+1 is I(t + 1) = R 0 (t)P (t).The transmission rate per infected individual (The probability of infecting each of the other individuals in the population in a unit of time), is given as β(t) = R 0 (t)/δN .Since we use infection pressure P (t) calculated over the entire SI, we can set δ to 1, such that β(t) = R 0 (t)/N .
The infection pressure exerted on each individual within the population is then calculated as which for low numbers of infected individuals (and therefore low P (t)) can be simplified to P I (t) = βP (t).The newly infected individuals (I(t+1)) are then randomly picked from a binomial distribution I(t + 1) = Binom(p = P I (t), N = S(t)).

VoC implementation
If a VoC is defined, we need to track the size of both the susceptible and the infected class of individuals for both the background variant and the VoC.This is first done retrospectively, for t ≤ T S , to determine the size of the classes at T S .
With an estimate of the VoC proportion over time (ρ V (t ≤ T S )), the past infections caused by the VoC can now be calculated as where ... denotes rounding to the nearest integer.The number of infections with the background variants thus follows as I B (t) = I(t) − I V (t).This also gives us the exact size of the susceptible class for both the background variant and the VoC Note that in the absence of new cases, the susceptible classes still need to be calculated forward in time before the incidence model is actually run, as forward vaccinations reduce their size.
The forward, prospective, model of the infection process, for t > T S , is then performed for both the background variants and VoC in the same way as with a single variant, splitting I(t), P (t), P I (t), and S(t) in a background variant ( * B ) and VoC ( * V ) class.The final epicurve then combines both variants again,

Care path (Within-hospital) model
To forecast bed occupancy in the hospital, the main endpoint of this work, we first simulate the path of each single patient admitted to the hospital through the general ward, ICU, and step-down units (care path model), as described in Donker et al ([4]).Then we join this model with the predicted incidence to forecast the hospital admission rate and the bed occupancy.
We use the following parameters to model the care path of individual patients through the hospital.
• Length of stay (LoS) distribution on general wards (GW) • LoS distribution on intensive care units (ICU) • LoS distribution on step-down units (SDU) • proportion of patients being transferred from GW to ICU • proportion of patients being transferred from ICU to SDU It is possible to add the proportion of patients directly admitted to ICU (H I ) to this list, but this is usually estimated from the data (see below).
Once all admissions and discharges are determined, we track the total number of patients on each ward type on each day of the forecast, presenting the final result of the forecast.Note that, despite modelling the general ward and step-down unit separately, we combine both wards into a single number of beds occupied on the general ward in any other step, because the occupied beds in both cases are usually reported as part of the same type (i.e., occupied non-ICU beds).

Admission rate
To estimate the proportion of reported cases being admitted to hospital, we need to compare the incidence of reported cases to the total number of patients in hospital at a given moment (prevalence).This can be achieved by taking the length of stay of admitted patients into account.We create a complete Length of Stay distribution L(t a ), defined as the proportion of patients still in hospital t a days after admission by simulating the care paths of a large number of patients (N P = 10000) with the same admission date using these parameters.For each day after admission, we then track how many patients are still in hospital (N H (t a )), resulting in the Length of Stay distribution, The admission rate can then be directly estimated using the current number of occupied beds in general wards B GW and in ICU B ICU which we assume to remain stable after the start of the simulation: The same method delivers the distribution of admitted patients on the GW L GW (t a ) and on the ICU L ICU (t a ).Note that L(t a ) = L GW (t a ) + L ICU (t a ), and while L(0) = 1, L GW (0) ≤ 1 and L ICU (0) ≤ 1.To be exact, L ICU (0) = H I , reflecting the directly to ICU admitted patients, and The admission rate can also be calculated for the separate wards, .
We can calculate the expected number of ICU beds under the assumption that no patients enter the ICU directly. .
The surplus of ICU patients (B ICU (t) − B * ICU (t)) is then caused by the direct admission of patients to ICU, and calculated as the ICU surplus over the total number of occupied beds.
Similar to α(t), we assume that α GW (t), α ICU (t), and H I (t) remain stable after T S .

Admission rate and VoC
In general, the assumption that the admission rate is stable over time reasonably holds unless the dynamics of the disease suddenly change.For now, we consider one such case: the occurrence of a variant with significant immune evasion, causing infections in individuals with immunity against the background variants through infection or vaccination (here called "reinfections", but it includes breakthrough infections to a large extent).In such a case, immunity in these individuals may protect against severe disease (to some degree).A sudden increase in the proportion of infections in these exposed individuals may lower the admission rate drastically.
We first need to calculate the number of reinfections.To do this, we track the proportion of individuals that have any immunity against any or either variant, assuming the second dose of the vaccine gives enough protection to reduce hospitalisation rates significantly, The number of reinfections is then determined by comparing the susceptible class (S) with the previous immunity class (U ).The sizes of these classes differ because of non-complete protection from vaccination and cross-immunity between strains.
The proportion of reinfections among all infections is thus The admission rate can then be adjusted to consider the effect of reinfections.This is done by adjusting future admission rates relative to the observed admission rate at the start of the simulation.First, we calculate the risk of admission for immunologically naive individuals (the "Primary hospitalisation risk"), π, at the start of the simulation given the relative admission risk for secondary infections γ.The admission rate over time then develops with X * (t) This admission rate then determines the number of individuals admitted to the hospital, entering the care path model, pulled from the binomial distribution Binom(p = α * (t), N = I(t)).
The care path model then depends on determining lengths of stays on each ward and movements between ward types for each individual admitted patient sampled from their respective distributions.
However, the care path model needs to determine the lengths of stays for the patients already in hospital at T S , that is B GW (T S ) and B ICU (T S ).To simulate their future discharge and transfer events, we create an admission record using a uniform number of patients per day for the 100 days preceding T S and simulate their care paths.
We then select all patients present in the GW or ICU at T S , creating the "Current" patient population in the hospital.From the current population, we randomly select N GW = B GW (T S ) and N ICU = B ICU (T S ) patients to recreate the discharges for the current population.

Dashboard server implementation
The main backbone of the model is written in R (version 4.1.2[13]) as an R-Shiny [14] (version 1.6) app running using shiny-server on an 8-core, 16GB ram, Ubuntu (version 20.04.3) VM server.An IP-hashed load balancer divides traffic over eight separate instances of shiny-server, reducing the number of concurrent users per shiny-server.We measure the performance of the model based on this infrastructure.
The most computationally demanding part of the care path model was written in the Julia programming language [15] to improve performance.Julia is rapidly gaining momentum as a tool for scientific computing due to higher performance compared to languages like R or Python without compromising ease of use.The care path model written in Julia takes approximately 1 hundredth of the processing time of the R equivalent.We use JuliaConnectoR (v.1.0.0.9009) [16] to allow communication between the Julia and R code bases.One of the drawbacks of Julia is that the functions require just-in-time compilation the first time they are used; the overhead related to compilation time exceeds the running time of the entire care path model.To address this issue, we generated a fully compiled version of the model code [17], virtually eliminating loading time.Once the main server is started, a concurrent Julia session is created and is shared between all shiny-server sessions.

User interface
The user interface of the dashboard is divided into sections (tabs) following the main modules of the model, in combination with a side panel showing the basic controls: the catchment area selection, the forecast start date, the number of simulation runs, and the simulation length in days.On each tab, parameter choices can be made that affect the current and further tabs, but not any previous tabs.In this way, we prevent unnecessary re-running of the simulations.
Tabs are ordered as follows: 1) Incidence, showing the daily reported number of cases with and without the nowcast.2) Vaccination, showing the cumulative number of vaccine doses administered, as well as the forecast of vaccinations.It also includes the option to change the assumptions on the future vaccinations (number of first doses and booster doses per day, and the minimum delay between 2nd dose and booster).3) Effective R, showing the time-varying R estimation based on the EpiEstim package.This tab also includes the option to use the ETS model and the option to include a variant of concern (VoC), together with the needed VoC parameters.4) Incidence forecast, showing the results of the incidence model, with the option to include the reporting pattern related to the weekdays.5) Bed forecast, reporting the results of the within-hospital care path model.The model is run using data on the current number of occupied beds in the hospital of interest, inputted either manually or by uploading a specifically structured file.6) Parameter choices.This tab includes all parameters included as default values in the other tabs, which should only be changed by users with more advanced knowledge of the underlying model.

Performance
Between 23-12-2021 and 07-04-2022, the online dashboard was visited 2352 times (counted as the number of unique sessions) by 575 users (counted as the number of unique machines by cookie UUID).On average, sessions lasted 22 minutes, but the distribution of session durations is heavily right-skewed, ranging between 0 seconds (i.e., immediately closed by the user) and 12 hours, and a median of 1.6 minutes (IQR 0.3 -5.9 minutes).There were active sessions during 471 hours of the 2523 hours.
For 219 hours, the dashboard served multiple concurrent users, specifically more than three users for 44.5 hours, to a maximum of six users on 23-12-2021 and 19-01-2022.Users primarily selected hospital catchments consisting of a single district (n=2507), followed by 44 districts (N=1426), coinciding with the number of districts in the state of Baden-Württemberg, and 412 districts (N=807), indicating all districts in Germany.The incidence model was called 6875 times and took 2.7s on average (median 1.4s IQR 0.8-2.5s),with a strong dependency on the total number of run-days (number of runs × length of the simulation, see figure 3).The care path model was called 5735 times and took 4.2s on average (median 1.9s IQR 0.9-4.1s).
As with the incidence model, the run time of the care path model depended strongly on the number of run-days.Additionally, it depends on the total number of beds occupied at the start of the simulation.On 157 occasions, users uploaded a file specifying the number of beds occupied in a specific hospital.This uploading was done by 24 unique users, of which two were responsible for the vast majority of uploading instances (74 and 33 times).The other forecasts were thus based on the number of beds occupied in each of the districts as automatically loaded from the central database.

Conclusion
This work shows that hospital-specific on-request production of bed demand forecasts during pandemics is possible.The developed on-request forecasting tool offers users maximum flexibility by allowing the definition of hospital-specific assumptions and/or starting conditions, while at the same time basing all forecasts on the same model.The combination of user-defined inputs within an existing model framework allows public health experts to adapt advanced forecasting methods to their local setting without having to build such models from scratch.Such an approach has the potential to make mathematical modelling of infectious diseases available to a much larger group of stakeholders, each with their specific set of questions to be answered.
Although our dashboard had a relatively modest number of users during the study duration, we are confident that it can cope with many more users, given the low proportion of time concurrent users were present.Nonetheless, there is certainly a limit to the traffic capacity of the dashboard, and particularly peak demand may be challenging.For this reason, we avoided publicly promoting our dashboard, relying instead on our network of hospitals and other healthcare organisations to circulate it.However, the dashboard was openly and freely available to any user over the entire period.Furthermore.the code is freely available and can be installed and run locally if hospitals need their own instance of the dashboard.
In conclusion, on-request forecasts are a fundamentally different method of providing infectious disease forecasts compared to the more common static reports.We showed how our dashboard implementation solves the specific set of challenges related to open-ended dynamic forecasting platforms, in particular computational requirements.The dashboard was successfully used by local healthcare providers, hospitals, and healthcare policymakers to evaluate incidence and hospital bed occupancy in Germany during the 2020-2022 COVID-19 pandemic.We argue that these on-request forecasts are much more helpful in informing stakeholders at a local level where health management decisions, such as cancelling elective surgeries, directly affect the bed capacity.This way, the pandemic or epidemic response can be driven in near real-time on the level where it matters most.
The views and opinions expressed herein are the authors' own and do not necessarily state or reflect those of ECDC.ECDC is not responsible for the data and information collation and analysis and cannot be held liable for conclusions or opinions drawn.

Figure 1 :
Figure 1: Model structure showing each of the modules, the data carried over between modules, and the possible user input in each module.

Figure 2 :
Figure 2: The user interface of the on-request COVID-19 bed demand forecasting model.A) Side panel with basic controls, B) Reported incidence, C) Effective R forecast, D) Vaccination forecast, E) Incidence forecast, and F) Bed occupancy forecast.

Figure 3 :
Figure 3: The run-time of B) the incidence model and B) the care path model, as a function of the total number of run-days and bed-run-days, respectively.
The other authors have no conflicts of interest to declare.