COVID-19 Delta Wave Caused Early Overburden of Hospital Capacity in the Bulgarian Healthcare System in 2021

We develop and apply our methodology to estimate the overburdening of hospitals in Bulgaria during the upcoming delta surge. We base our estimations on an exponential risk model from the UK. Still, the methodology is generally applicable to all risk models, depending on age. Our hypothesis is that during the delta wave in Bulgaria, the system experienced a burden from late August due to decreased capacity. This will explain most of the excess mortality during the wave. We estimate the number of people from the active cases in need of hospitalization and intensive care.


Introduction
Coronavirus disease 2019 (COVID-19) continues to be a worldwide threat to mankind, requiring specialists, policymakers, and governments to address several issues that extend far beyond the health and well-being effects of this pandemic [1]. While the pandemic's immediate health consequences play out, we need research and actions to be reconfigured to reduce risk, establish continuity, and enhance resilience for potential recovery [1].
We are witnessing conflicts in the ultimate social determinant of health as a result of the pandemic, which is producing a wide variety of challenges ranging from health policy constraints to worries about the distribution and access to healthcare [2]. Reorienting healthcare resources to COVID-19 management has reduced the capacity of already overburdened and underfunded health systems to handle other disease loads. Furthermore, discontinuing regular treatments and interventions and follow-up and immunization programs lead to outbreaks of avoidable transmissible illnesses, increased cancer incidence, and the number of challenging medical conditions in the late stages [1].
Furthermore, COVID-19 has grown difficult even for the most enduring healthcare systems. To stop the spread of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) coronavirus, extraordinary measures were implemented worldwide. Lockdowns appear to be one of the only successful interventions, although at the expense of economic slowdowns and restrictions on human liberties [2].
However, the current (fourth) wave in Bulgaria, typical of respiratory diseases, startled mankind once more and did more serious harm than the first ones. As a result, the optimism has started to fade. Importantly, developing herd immunity to control a pandemic through natural infection will undoubtedly result in millions of fatalities [3]. Furthermore, the delta variant of SARS-CoV-2 also raised people's concerns, including healthcare workers [4].
At the outset of the pandemic, in this period of despair and lack of control over the illness, the discovery and implementation of COVID-19 vaccines provided us a glimpse of hope. It appears that one of the most effective public-health interventions-vaccinationmay be used to tackle the pandemic [5].
Our hypothesis is that many people who need hospitalizations do not receive adequate and in-time hospital care due to changed admission rules. Therefore, the risk of hospitalization grows exponentially with age. We know that the share of people over 60 years is the main predictor for the ratio of hospitalized patients to active cases. The latter ratio should have at least a linear or slightly convex relationship with the first one. It should not tend to a constant value before the first ratio reaches 1 s. This is a basis for our methodology that allows us to qualitatively judge whether there is hospital overburdening and approximately estimate the need for hospital beds. Furthermore, our models for predicting deaths from cases by age groups show that mortality risk is also exponential, which reinforces the validity of the hospitalization risk model (deaths should be proportional to cases) and also shows a slight lag of 0 to 7 days between discovery and death for the 60+ age groups [10]. This small lag makes proper hospitalization in time crucial for controlling the mortality risk. If the distribution of new cases by age is f (t), where t is time in a year (age), we can obtain the distributions of the share of hospitalizations by age f h (t) by multiplying it with the function g(t): (1) gives the risk of hospitalization per age in percentage [11], and we need to normalize (2).

Methodology for
This distribution varies very little with new cases day by day. This is due to the non-stationarity of the process and can be estimated in various stages of the pandemic, for example, for different dominating strains, by taking the cumulative cases for the period. The distribution for the total cumulative cases has a very stable value, because it reflects stable population demographics and contact patterns across ages that lead to this distribution of infected people, as a general picture of the hospitalizations. Since the data are not differentiated by sex, we expect to have bimodal distributions. We chose generalized Gaussian function (3) and estimated it with adjusted R 2 = 0.984, standard error (sum of squared residuals) SSE = 0.0002386, and root mean squared deviation RMSE = 0.008918.
The age groups in the data are 0-19, 20-29, 30-39, 40-49, 50-59, 60-69, 70-79, 80-89, and 90+. We take it from the Open Data Portal [10] and we evaluate by taking the middle of the intervals, with the points 10, 25, 35, 45, 55, 65, 75, 85, and 95 years of age ( Figure 1). To get our distribution, we need to integrate (3):  To treat (3) as probability density distribution, we need to scale with 10.4255, since for the appropriate age interval 10.4255 (4) To get our distribution, we need to integrate (3): This is the share of all individuals hospitalized up to 3 September 2021. So, for 46,0691 cases, this gives us 60,626.1 hospitalizations (Figure 2). With this distribution, the share among 60+ is 66.66%, which is cumulative for 3.09%. The same share in the new cases is 52%. Thus, we can see disproportionally many patients over the age of 60 with the model of [11]. With this distribution, the share among 60+ is 66.66%, which is cumulative for 3.09%. The same share in the new cases is 52%. Thus, we can see disproportionally many patients over the age of 60 with the model of [11]. The risk of hospitalization is exponential-e 0.044t [11], which allows us to use the population above a certain age threshold and its ratio to the whole population (in the new cases). Getting the mean risk from the integral of (1) divided by the length of the interval (0,100) gives us a mean risk of 18.2843% for the general population group. The mean risk for under 60 is 2.96%, which is under 1/6 of the total risk. For that reason, the ratio of people over 60% is enough to determine the hospitalization risk with good precision, to avoid using the whole distribution, which needs to be obtained from data for each day if we follow the approach, we described in Section 2.1.

2.2.2.
Modeling the Dependency between the Ratio of People over 60, the Active Cases, the Hospitalized Patients, and Patients in ICU There will be a fixed share of hospitalizations for a specific fixed distribution of new cases by age. Instead of distribution, we use the ratio of people over 60 years. In this scenario, the model for the number of patients H(a, r 60+ ), depending on active cases a and the ratio r 60+ , is linear: Such a linear model can be estimated with linear regression, given certain essential factors: The number of hospitalized patients is cointegrated with the number of active cases and the ratio of people over 60 years. This model can be estimated for various data ranges and would give different values (m). To fit this model, first, we need to prepare our data.

Preparation of Data
We take active cases from the Open Data Portal [12]. The ratio of people over 60 we calculate for every new day. We differentiate the cumulative cases by age to obtain the daily cases by age, containing multiple periods. We analyzed frequency spectrums of active cases and this ratio to filter that noise and chose 21 days for the moving mean ( Figure 3). We used 25 July2021 as the starting point of modeling when the ratio was lowest. The delta variant became dominant in Bulgaria (officially at the beginning of August, but data are delayed 7-14 days). This period of monotonic increase allows us to apply linear regression on cointegrated time series. Such a linear model can be estimated with linear regression, given certain essential factors: The number of hospitalized patients is cointegrated with the number of active cases and the ratio of people over 60 years. This model can be estimated for various data ranges and would give different values ( ). To fit this model, first, we need to prepare our data.

Preparation of Data
We take active cases from the Open Data Portal [12]. The ratio of people over 60 we calculate for every new day. We differentiate the cumulative cases by age to obtain the daily cases by age, containing multiple periods. We analyzed frequency spectrums of active cases and this ratio to filter that noise and chose 21 days for the moving mean ( Figure 3). We used 25 July2021 as the starting point of modeling when the ratio was lowest. The delta variant became dominant in Bulgaria (officially at the beginning of August, but data are delayed 7-14 days). This period of monotonic increase allows us to apply linear regression on cointegrated time series.

Testing for Cointegration
We used the Engle-Granger cointegration test in Octave for the ratio r 60+ (t) and hospitalization cases H(t), and they are cointegrated with p-value = 0.001 < 0.05. For hospitalized and active cases (moving averages with the period of 21 days), the Johansen cointegration test [13] in Octave gives cointegration with a p-value = 0.001 < 0.05. Furthermore, we know in advance that hospitalizations are a certain percentage of new cases and total hospitalizations from active cases, depending mostly on that ratio from our analysis, which gives a theoretical background to our modeling efforts.

Linear Regression with PCA-Hospitalizations
The correlation between the ratio r 60+ (t) and active cases a(t) is very high-95%, due to the stage of growing hospitalizations and growing ratio of 60+. Nonetheless, we need both variables, since we want to check how the coefficient k h , representing the ratio between hospitalizations and active cases, changes with the increase of r 60+ (t). For that reason, we employ principal component analysis and make the regression with the transformed data. We first make a linear model for the whole interval to assess the fit. Then, we will use it by removing data points backward in time-from 3.09 15 points backward as a compromise between the number of observations needed for regression and the number of estimated coefficients.
The fitted model for hospitalizations with data from 25.07 to 3.09 has almost perfect characteristics (Table 1), as shown in Figure 4, including adjusted R 2 = 1. Here, k 60+ = 2583.2 and k a = 0.13027.   The model here says that on average, 13.027% of active cases become hospitalized above some number of patients, which depend linearly on the ratio . This is a linear approximation of a nonlinear relationship, which is similar to the approximation a tangent line makes to a curve around the point of tangency, or a plane, that is tangent to a surface-which is our case. We are constructing a plane in three dimensions-one is the active cases, the other is the ratio, and the third dimension is , in (7). This is a local The model here says that on average, 13.027% of active cases become hospitalized above some number of patients, which depend linearly on the ratio r 60+ . This is a linear approximation of a nonlinear relationship, which is similar to the approximation a tangent line makes to a curve around the point of tangency, or a plane, that is tangent to a surface-which is our case. We are constructing a plane in three dimensions-one is the active cases, the other is the ratio, and the third dimension is H(a, r 60+ ) in (7). This is a local approximation of a surface. We can reconstruct the surface from a series of linear approximations-making this linear regression for m, m-1, m-2 . . . m-14 points and obtaining the values for k a c , which represent the proportion of hospitalized patients and how it changes with r 60+ .
Obtaining a Model for the Dependency between r 60+ and k a c With the increase of r 60+ , the coefficients k a increase toward a fixed value. The exponential nature of the risk with age suggests the lowest possible growth (as a limit) to be linear or the slowly aperiodic behavior of the process-the increase of the ratio of 60+ leads almost linearly to a value in an interval e 0.044 * 60 , e 0.044 * 100 = [14.01%, 81.45%] with an average of 38.3169% (7a). A simple way to understand this is to imagine that the risks for people over and under 60 (7b) are uniformly distributed with age, and their expected risks fully represent them. Then, the combination of risks depends linearly on the ratio of people over 60. We get an approximation for the expected risk of the distribution of new cases (8). Thus, the choice of the general model is founded on knowledge from the risk model per age and is not just from the data. This is important, since the data contain risk and hospital capacity overburdening information. The theory allows us to distinguish them by justifying the selected method of extrapolation of data.
The general model is (10)-for hospitalizations and for ICU admissions alike, it is aperiodic, but its usage is as an almost linear model, with very low convexity (small exponential constant b) To calculate the results for hospitalization-we repeat the linear regression fit with decreasing time windows-3.09-25.07, 2.09-25.07 . . . .-15 consecutive times. Then, we fit model (8). For each of these points, we obtain good regressions with R 2 = 1 ( Figure 5) and similar characteristics such as p-values and F statistics to those in Table 1. Then, we fit (9) with these data, and we obtain the model (10). It is applicable for r 60+ > 0.05: With that model, we can calculate for different numbers of active cases and different values of the simple model (11), as shown in Figure 6.

ICU Beds
The exact same methodology is applied to ICU beds. We assume that the exponential risks of hospitalizations have a fixed proportion of cases needing intensive care and that this risk also grows exponentially. Thus, the general linear regression for the whole period is as good as it is for hospitalizations.
By using the same approach as for hospitalizations, we can fit (9) (Figure 7) and obtain (12) (Figure 8): With that model, we can calculate for different numbers of active cases and different values of r 60+ the simple model (11), as shown in Figure 6.
Healthcare 2022, 10, 600 8 of 13 With that model, we can calculate for different numbers of active cases and different values of the simple model (11), as shown in Figure 6.

ICU Beds
The exact same methodology is applied to ICU beds. We assume that the exponential risks of hospitalizations have a fixed proportion of cases needing intensive care and that this risk also grows exponentially. Thus, the general linear regression for the whole period is as good as it is for hospitalizations.
By using the same approach as for hospitalizations, we can fit (9) (Figure 7) and obtain (12) (Figure 8)

ICU Beds
The exact same methodology is applied to ICU beds. We assume that the exponential risks of hospitalizations have a fixed proportion of cases needing intensive care and that this risk also grows exponentially. Thus, the general linear regression for the whole period is as good as it is for hospitalizations. By using the same approach as for hospitalizations, we can fit (9) (Figure 7) and obtain (12) (Figure 8):

Projections for Hospital Capacity Overburden, Based on That Methodology
Estimating the total number of people left out of the hospital system requires additional methodology. By comparing these models to the actual data, we can assess the difference between the number of expected hospitalizations and the number of actual hospitalizations. The same can be done for ICU patients. First, the overburdening of the system can be assessed visually from Figures 9 and 10, which are the already fitted models with added data from 4.09 to 21.09. Second, the expected number of hospitalized patients and patients in the ICU can be calculated from the fitted models and the number of active cases (all data are with moving average with a period of 21 days) (Figures 11 and 12).

Projections for Hospital Capacity Overburden, Based on That Methodology
Estimating the total number of people left out of the hospital system requires additional methodology. By comparing these models to the actual data, we can assess the difference between the number of expected hospitalizations and the number of actual hospitalizations. The same can be done for ICU patients. First, the overburdening of the system can be assessed visually from Figures 9 and 10, which are the already fitted models with added data from 4.09 to 21.09. Second, the expected number of hospitalized patients and patients in the ICU can be calculated from the fitted models and the number of active cases (all data are with moving average with a period of 21 days) (Figures 11 and 12).

Projections for Hospital Capacity Overburden, Based on That Methodology
Estimating the total number of people left out of the hospital system requires additional methodology. By comparing these models to the actual data, we can assess the difference between the number of expected hospitalizations and the number of actual hospitalizations. The same can be done for ICU patients. First, the overburdening of the system can be assessed visually from Figures 9 and 10, which are the already fitted models with added data from 4.09 to 21.09. Second, the expected number of hospitalized patients and patients in the ICU can be calculated from the fitted models and the number of active cases (all data are with moving average with a period of 21 days) (Figures 11 and 12).

Discussion
The present methodology has several limitations connected to the data. First and foremost, the exponential risk model we used is obtained from data in the UK, where the numbers of tests per 100,000 are several times higher. The system of contact tracing and isolation is much more thorough. The data, released in the open portal from the Ministry of Health, indicate that new hospitalizations follow the shape of the exponential model but are almost double, which is expected because Bulgaria has a double case fatality ratio of 4%, instead of the average 2%, which indicates that we find half of the infected in comparison with many other countries. However, our modeling is done so that we do not depend on that discrepancy between the model and actual new hospitalizations-we assess the ratio of the active cases and active hospitalizations, in which the errors in the numerator and denominator cancel out. Secondly, we assume that 21 days of averaging are enough to compensate for delays in corrections of active cases to avoid overestimating the hospitalizations ratio-if we have more official active cases than actual ones due to delays in correction.
The network dynamics model for the transmission of COVID-19, proposed by Zhu et al., showed a feasible approach to simulate the COVID-19 epidemic with different interventions, thus providing information and advice on how and when to open large-scale public facilities [14].
Next, we assume that the risk for ICU admission is also exponential and can be approached with the same linear model for the ratio of people over 60 years. Finally, we accept the discrepancies between model and actual cases as evidence for hospitals above capacity and not for evidence against the model.
Here, have some anecdotal evidence for limited capacity in hospitals in many regions of Bulgaria that were hit first in the delta wave, such as Targovishte and Burgas areas [15]. We also rely on the official policy of the Ministry of Health to transfer more cases into home treatment under supervision from general practitioners [16]. This is the core of our hypothesis-that the exponential risk of hospitalizations can be used to estimate the burden on hospitals and to provide evidence for exceeded capacity. Once again, it was confirmed that effective government communication is critical for public health crises [17], such as that observed in our country.
We show that there is a change after the alpha wave in the spring that leaves patients for home treatment and thus increases their mortality risk. The hospital system was left unprepared for the delta wave, leading to a lower percentage of people treated without delay. Instead of expanding the capacity for the expected large wave that came with two months delay, the system has shrunk it.

Discussion
The present methodology has several limitations connected to the data. First and foremost, the exponential risk model we used is obtained from data in the UK, where the numbers of tests per 100,000 are several times higher. The system of contact tracing and isolation is much more thorough. The data, released in the open portal from the Ministry of Health, indicate that new hospitalizations follow the shape of the exponential model but are almost double, which is expected because Bulgaria has a double case fatality ratio of 4%, instead of the average 2%, which indicates that we find half of the infected in comparison with many other countries. However, our modeling is done so that we do not depend on that discrepancy between the model and actual new hospitalizations-we assess the ratio of the active cases and active hospitalizations, in which the errors in the numerator and denominator cancel out. Secondly, we assume that 21 days of averaging are enough to compensate for delays in corrections of active cases to avoid overestimating the hospitalizations ratio-if we have more official active cases than actual ones due to delays in correction.
The network dynamics model for the transmission of COVID-19, proposed by Zhu et al., showed a feasible approach to simulate the COVID-19 epidemic with different interventions, thus providing information and advice on how and when to open large-scale public facilities [14].
Next, we assume that the risk for ICU admission is also exponential and can be approached with the same linear model for the ratio of people over 60 years. Finally, we accept the discrepancies between model and actual cases as evidence for hospitals above capacity and not for evidence against the model.
Here, have some anecdotal evidence for limited capacity in hospitals in many regions of Bulgaria that were hit first in the delta wave, such as Targovishte and Burgas areas [15]. We also rely on the official policy of the Ministry of Health to transfer more cases into home treatment under supervision from general practitioners [16]. This is the core of our hypothesis-that the exponential risk of hospitalizations can be used to estimate the burden on hospitals and to provide evidence for exceeded capacity. Once again, it was confirmed that effective government communication is critical for public health crises [17], such as that observed in our country.
We show that there is a change after the alpha wave in the spring that leaves patients for home treatment and thus increases their mortality risk. The hospital system was left unprepared for the delta wave, leading to a lower percentage of people treated without delay. Instead of expanding the capacity for the expected large wave that came with two months delay, the system has shrunk it.
The lack of hospital beds and late hospitalization in some patients with COVID-19 leads to the late application of adequate treatment. This leads to an increase in severe disease cases, long COVID-19, and mortality. However, the problem with the lack of hospital beds is more due to the lack of trained medical staff to treat patients. This situation often leads to the need for urgent reorganization, which does not always comply with the rules for preventing nosocomial infections. As a result, the risk of healthcare-associated infection for both staff and patients increase. This creates fear in many patients to seek medical attention in a timely manner. These may be patients with COVID-19 or other diseases. This will inevitably complicate their condition and increase the risk of death. For this reason, increased mortality is expected in other diseases as well.

Conclusions
We propose a methodology to estimate how many patients from the active cases should be in hospital based on a model for the hospitalization risk. We use the exponential nature of the risk to simplify the task and use the ratio of people over 60 years old to all people in cases. We infer the deviation from the expected trajectory when both hospitalized patients and the ratio are growing and cointegrated with the active cases' growth. We infer the age distribution of hospitalizations from that model and the age structure of cases.
Our work suggests that during the growth of delta cases in the summer, the policy of accepting patients changed and did not correspond to the exponential risk of hospitalizations, based on age. That implies the limited capacity of the hospital system and people left at home that should be accepted in the hospital. This is especially visible for intensive care beds when we have a decrease of occupied beds with the increase of the ratio of 60+ people in active cases. There is some uncertainty in estimating actual beds needed in our model. Still, its main goal is to show if there is a discrepancy between increases in active cases and hospitalized cases due to hospital system overburdening.