Abstract
As a response to the pandemic caused by SARSCov-2 virus, on 15 March, 2020, the Republic of Serbia introduced comprehensive anti-epidemic measures to curb COVID 19. After a slowdown in the epidemic, on 6 May, 2020, the regulatory authorities decided to relax the implemented measures. However, the epidemiological situation soon worsened again. As of 15 October, 2020, a total of 35,454 cases of SARSCov-2 infection have been reported in Serbia, including 770 deaths caused by COVID19. In order to better understand the epidemic dynamics and predict possible outcomes, we have developed a mathematical model SEIRDS (S-susceptible, E-exposed, I-infected, R-recovered, D-dead due to COVID19 infection, S-susceptible). When developing the model, we took into account the differences between different population strata, which can impact the disease dynamics and outcome. The model can be used to simulate various scenarios of the implemented intervention measures and calculate possible epidemic outcomes, including the necessary hospital capacities. Considering promising results regarding the development of a vaccine against COVID19, the model is enabled to simulate vaccination among different population strata. The findings from various simulation scenarios have shown that, with implementation of strict measures of contact reduction, it is possible to control COVID19 and reduce number of deaths. The findings also show that limiting effective contacts within the most susceptible population strata merits a special attention. However, the findings also show that the disease has a potential to remain in the population for a long time, likely with a seasonal pattern. If a vaccine, with efficacy equal or higher than 65%, becomes available it could help to significantly slow down or completely stop circulation of the virus in human population. The effects of vaccination depend primarily on: 1. Efficacy of available vaccine(s), 2. Prioritization of the population categories for vaccination, and 3. Overall vaccination coverage of the population, assuming that the vaccine(s) develop solid immunity in vaccinated individuals. With expected basic reproduction number of Ro=2.46 and vaccine efficacy of 68%, an 87%- coverage would be sufficient to stop the virus circulation.
1. Introduction
On 11 March, 2020, the World Health Organization characterised the disease caused by the novel SARS Cov-2 virus as a pandemic [1]. Initial epidemic outbreak in China spread outside the Wuhan area, and subsequently on a global scale. On 6 March, 2020, the first case of the novel corona virus infection was reported in the Republic of Serbia. Taking into consideration the escalation of the disease and limited effects of the initially implemented measures, the state of emergency was declared throughout the country on 15 March, 2020. Comprehensive anti-epidemic measures (e.g. lockdown of entire country) were introduced in the entire country [2]. Due to the absence of specific pharmaceutical intervention, Serbia, like other countries, implemented an anti-epidemic strategy based on social distancing, school and university closures, reduced number of workers present in the workplaces, closure of places of worship for public religious services, reduced working hours of cafés and restaurants, avoiding mass gatherings, events, sports games, tracing and identification of infected people’s contacts, etc. After a slowdown in the epidemic, as shown in the relevant officially published data, the regulatory authorities decided to relax the introduced measures on 6 May, 2020. However, the epidemiological situation soon worsened again, resulting in the reinstatement of some measures, as well as the introduction of new measures [2]. Although the return of extensive measures has yielded favourable results, the further development of the epidemic is not clear.
For these reasons, mathematical modelling has a crucial role in understanding the epidemic and predicting possible outcomes. Modelling is a particularly useful tool for devising strategies for combating the epidemic, capacity planning, and selection of efficient measures, especially in the absence of specific pharmaceutical treatments [3, 4, 5]. Mathematical modelling based on differential equations dates back to the first half of 20th century. In 1927, Kermack and McKendrick developed the basic model of disease transmission consisting of three compartments: susceptible (S), infected (I) and removed (R). The model was based on a connected system of nonlinear differential equations as a special case of the general epidemiological model [6, 7]. Subsequent models, became more complex and adapted to the needs of research [5].
Since the outbreak of COVID19, many published papers have dealt with the implementation of mathematical modelling and prediction of possible outcomes of COVID19 epidemics. Most of these research efforts have been based on the implementation of the SIR (susceptible-infected-removed) model. Many of the other models provide a clear picture of dynamics of COVID19 spreading, including the overloading of the relevant health systems. For example, Ferguson et al., developed one of the first models for COVID19 simulation, which was, among other things, used to plan the health care resources [8]. Wu et al., developed the SEIR model to examine the dynamics of SARS-Cov-2 transmission from person to person. This model was also used to calculate the basic reproduction number Ro, which we use in this paper as one of the key parameters [9]. The classical SIR model assumes that there is homogenous mixing of infected and susceptible populations and that the total population is constant and does not change over time. Moreover, according to the classical SIR model, there is a monotonous decline in susceptible population towards zero [10]. However, such assumptions are not objective in the case of COVID-19 spreading and they are the basic problem in the modelling of this pandemic. In reality, human population fluctuates constantly [11]. In order to account for such fluctuation, and better understand the COVID 19 epidemic in the Republic of Serbia, we have employed mathematical modelling of the epidemic using the available data on the characteristics of the disease, such as incubation period, latent period, recovery period, severity of clinical signs, and mortality rate caused by COVID19. Unlike the classic SIR model, SEIRDS (S-susceptible, E-exposed, I-infected, R-recovered, D-dead due to COVID19 infection, S-susceptible) epidemic model, developed for this research, simulates the spreading of COVID 19 in an open population. Taking into account that the population is constantly changing and that various measures are applied for different strata or subgroups of the population (such as preschool children, children attending primary school, high school students, students, employees, the unemployed and retirees), as well as changes in the intensity of applied measures, we have proposed the use of a model that takes these circumstances into account. Based on input disease parameters taken from scientific literature and specific data related to Serbia, this model simulates daily disease occurrence, including the number of hospitalised patients and cases which require intensive care. The model also predicts the expected number of deaths, as well as hospital capacities necessary to accommodate the patients. It provides a possibility to simulate different scenarios of disease control and intervention measures. Considering the expectations of successful development of the vaccine against COVID 19 in the near future, we added an option to model vaccination of different strata of the population as a set of disease control strategies.
2. Methodology
This section presents the research methodology and the proposed model, which was used to predict the further dynamics of the epidemic in Serbia. We also presented the data that were used to model the epidemic, a simulated strategy to combat COVID 19, and a sensitivity analysis.
2.1. SEIRDS mathematical model
Classical SEIRDS model divides the population into compartments, i.e. groups, and follows the disease dynamics at all times. The population is divided into the following compartments: the portion of the population susceptible to the infection is denoted by S, those latently infected with SARSCov-2 (exposed to) are denoted by E, the infected individuals who are able to spread the disease are denoted by I, the ones recovered from the infection are denoted by R, and those who died due to disease with D. Assuming that individuals mix homogenously, the force of infection λ (the rate at which susceptible persons are infected) is related to per capita contact rate β. Also, the risk of infection is closely related to the number of infectious individuals in the population It. It depends on the number of infectious individuals (It) and how frequently they make contacts with other persons. In a situation of homogenous mixing among the population, the force of infection λ can be express ae follows:
The change of rates in every compartment per unit time in SEIRDS model is presented in the following series of differential equations:
where f is rate of onset of infections expressed as the reciprocal of the latent infection period, r is the rate at which infectious individuals are recovered, δ is the rate at which infectious individuals die from COVID 19 infection and ω is rate of waning of immunity. The total population at any particular interval of time t is:
where parameters b and m are per capita daily birth and death rates unrelated to COVID19.
However, considering that implemented anti-epidemic measures against COVID 19 do not have an identical impact on population’s age subgroups and that COVID19 pathogenesis varies in different age subgroups, we propose the use of multi-compartment version of standard SEIRDS model. In this model susceptible population was further stratified within the compartment S according to age and occupations. Grouping into various strata was done according to the real age structure of the Republic of Serbia population as follows: pre-school children (Sps), elementary school children (Ses), high school children (Shs), college students (Scs), unemployed population (Sua), employed population (Sea), and elderly/retired (Sr) (Table 1). To simulate the epidemic progression through different population strata-subgroups, we used appropriate, stratum-specific, model parameters and factors of effective contact reduction (anti-epidemic intervention measures - ρ), which were adapted to the relevant population groups: lockdown of the entire country, closures of pre-school, school and college sessions closures, reduced number of workers allowed in the workplaces, work from home, restrictions of mobility of elderly people, etc. During the simulation, we monitored the effects of various levels of contact reduction, ranging from 25% to 75%, taking into account the realistic possibilities of maintaining a minimum work process, functioning of the society and feasibility of such measures.
Structure of different population strata in the Republic of Serbia [23]
Given that intervention measures, applied in response to the emergence of COVID19, are not the same for all population strata, homogeneous mixing can be expected only within same population stratum. The rate of effective contacts β, after the application of intervention measures, is no longer identical at the level of all strata of the population. Effective contacts are limited by the intensity and types of measures applied and are identical only when it comes to individuals within the same population strata. Furthermore, persons in different population strata become infected at different rates depending on how frequently they interact with other persons in their own subgroup and other subgroups. If we assume that force of infection differs between different strata of population, the equation for the force of infection would be as follows:
where λi(t) is force of infection in the ith population strata, βij is the rate at which susceptible persons in the ith population strata and infectious persons in jth population strata come into effective contact per unite of time, and Ij(t) is the number of infectious persons in jth population strata. Also, in this model the number of recovered and dead is conditioned with different ages and genders.
Now our model will be expressed as follows:
In this model susceptible, exposed, infectious, recovered, deaths and total population are:
2.2. Determining the herd immunity threshold and control of CIVID 19 by vaccination policy
Considering the undergoing worldwide efforts to develop a vaccine against COVID 19 and promising results, we extended the model to simulate and analysed the effects of a hypothetical vaccination on the epidemic dynamics, and to estimate the extent of coverage of vaccination which can interrupt the chain of infection. The control of COVID 19 by vaccination means targeting the entire susceptible population with mass vaccination until critical herd immunity achieved. In such situation there is a “race” between the exponential growth of epidemic and mass vaccination. The level of herd immunity threshold (HIT) can be calculated by the following formula:
and the critical vaccination coverage required to achieved herd immunity can be obtained by multiplying herd immunity threshold with reciprocal value of vaccine efficacy, ve:
Most people infected with SARS-CoV-2 develop an immune response followed by the development of specific antibodies between 10 and 21 days after getting infected [12]. Specific IgM and IgG antibodies against SARS-CoV-2 develop 6 to 15 days after the onset of the disease [13-17]. According to some studies, the presence of antibodies has been confirmed in less than 40% of the patients within 1 week after the onset of the disease, whereas percentage reaches 100% of subjects 15 days after the onset of disease [18]. Although duration of the immune response against CVOVID 19 is still unknown, comparing with other coronaviruses, immunity wane within 12 to 52 weeks after the first symptoms appear [19], while in the case of SARS-CoV-1 infection, the presence of IgG antibodies was confirmed in 90% and 50% of infected patients, respectively, over two and three years, respectively [29]. Based on these findings, we assumed that in the event of the development of a successful vaccine, immunity against the SARRSCov-2 virus could last for a year, as well as after recovery after a natural infection.
For the purpose of modelling, additional compartment to the model was added denoted with V(t), in which there are vaccinated persons who have successfully developed protective immunity after vaccination. The vaccination parameter, υ, is the daily rate of vaccination of susceptible population and it represents the proportion of susceptible population immunized per unite time. The critical daily rate of vaccination, υc, is ѱc = (b+ω)(Ro -1), required to interrupt the infection [5]. The basic reproduction number under the vaccination is Rop = (1-p)Ro. The proportion of effectively protected persons, p, is conditioned by parameters the vaccine efficacy, ve. This parameter represents a proportion of person who successfully developed protective immunity after vaccination, whereas total number of actively protected individuals in time t is V(t) = number of vaccinated x ve [5]. In this compartment the daily rate of waning of immunity at which immunity of vaccinated population fades out is ω, and it is reciprocal to the period of lasting of immunity. Vaccinated persons, after losing their immunity, become sensitive again and removed to the compartment S. The change of rate in this compartment per unit time is as follow:
The compartment S(t) is slightly modified as follow:
The other of compartments of SERIDS model remain unchanged.
2.3. Model parametrisation
In proposed model, β is per capita effective contact rate at which specific persons come into effective contact per unit of time. An effective contact is defined as a contact sufficient to cause disease transmission [10, 20, 21]. We calculated the parameter β by using the formula:
where Ro is a basic reproduction number of the disease, i.e. the average number of newly infected people with COVID 19 (secondary infection cases), infected by one infectious individual in a totally susceptible population, N is total population, and TR is the average duration of infectious period [10, 20, 21]. The R values of 2.46 and 3.1 are adopted from the relevant literature. The Ro values were based on the data obtained during the initial phase of the epidemic in Italy [22]. Since the implemented measures and disease transmission were simulated through various population strata, we corrected the β parameter with a relevant, stratum-specific contact reduction factor ρi. In this way, we obtained the per capita contact rate specific for each separate stratum based on following formula. βi = β(1-ρi). The values of the ρi factor in different population strata ranged from 0.25 to 0.75 (effective contact reduction ranged from 25% to 75%).
Parameters such as daily birth and death rates were calculated based on the data published by the Office of Statistics of the Republic of Serbia, and data published by the World Bank regarding the life expectancy in the Republic of Serbia [23, 24]. The latter study reported that the life expectancy in Serbia was 79.06 years in 2017 [23]. By using this figure, we expressed the daily death rate as a value reciprocal to life expectancy m = 0.000036. We calculated the daily per capita birth rate of b= 0.000025 based on the figure of 9.2 births in the Republic of Serbia per 1000 people in a year. These estimates were needed to realistically simulate fluctuations of the total population. To simplify the calculations, we assumed that the general morality rate m is applicable for all population strata.
The infectivity rate, i.e. the rate of transfer from compartment E to I, was derived as a value reciprocal to the COVID 19’s average latent period. The data on the average duration of latent infection (f-1) and the average period during which an infected person is shading the SARSCov-2 virus (TR) were adopted from the relevant literature as f-1 = 4.6 days and TR = 6.5 days, respectively [8]. The same study is used as a source of the data: on the percentage of hospitalised patients and those whose therapy requires intensive care, used for prediction of required hospital capacities, as shown in Table 1 [8].
Parameters such as δ and r are related to the infectious fatality rate (IFR) for COVID 19 and average times taken from infection to death (TD) or recovery (TR). These parameters were calculated using the following formulas:
The IFRs, shown in Table 3, were taken from literature and compared with local IFR value which was calculated based on officially registered deaths published by the health system of the Republic of Serbia [2]. The Calculation of local IFR is presented in section 2.4. Population data, (e.g. total population, age structure, and stratification) are presented in Tables 1 and 2. A summary of all model parameters is given in Table 3.
Age structure of the population of the Republic of Serbia and expected percentage of hospitalised patients, patients in intensive care, and death rate caused by COVID19.
SEIRD model parameters
2.4. Setting disease control scenarios
Five different scenarios were developed for simulating the COVID19 epidemic control based on non-pharmaceutical interventions. SC1 implies a base-case scenario where the epidemic spreads in susceptible population without any anti-epidemic measures being implemented. In the other scenarios, the extent of contacts was reduced, for each population stratum separately, according to objective possibilities and measures which were implemented during the actual epidemic in the Republic of Serbia. Scenarios are presented in table 4.
Description of different simulated non-pharmaceutical intervention scenarios
Additional four scenarios of control of COVID 19 were simulated based on vaccination policy. We assumed that vaccine efficacy was 50%, 68%, 80%, and 85%. The initial conditions assumed that all other anti-epidemic measures are excluded from the model and replaced with mass vaccination. Indicators of epidemic dynamics were monitored, such as: CI, hospitalized patients, patient in intensive care units and deaths.
2.5. Model sensitivity analysis and calibration
Considering the world experience with detection of COVID19 cases, as well as the unreliable data on COVID19 infections which are currently available worldwide, model calibration is very challenging, and can result in obtaining inaccurate values for the parameters [26]. This is especially due to the facts that a significant percentage of the infected individuals do not exhibit any symptoms. The other issue is small percentage of tested population [26].
As part of the national seroepidemiological study, 1,006 subjects were tested in Serbia from May 11th to June 25th, 2020, for the purpose of estimating the extent of COVID-19 infection. According to the published data, seroconversion was confirmed in 6,4% of the subjects. On the other hand, a total of 13,372 cases of the infection were reported, which means that those who were infected constitute around 0.19% of the overall population. However, it is our opinion that the data on reported deaths caused by COVID19 infection is more reliable for use in model calibration, e.g. infection fatality rate. Alex et al. reached a similar conclusion when simulating COVID19 by using the SEIRD model with heterogeneous diffusion [26]. When we compared the data recorded during the beginning of the epidemic in Serbia with the results obtained during the simulation, such as the initial doubling time, the two data series matched well. However, later, the obtained results did not match well the officially registered data on the number of infected, especially after the beginning of the implementation of measures in Serbia. We attribute these differences to the methodology by which official authorities register cases of infection, and collect the data. However, on October 15th, 2020 official authorities have published new epidemiological estimates of COVID19 population infection, which is 20% of the total population, which approximately coincides with the results obtained from our simulation, (e.g. 1,396,521 infected cases obtained from simulation compared to 1,120,460 cases estimated by the Serbian official authorities).
The parameter that determines the number of deaths is the IFR. It is the number of persons who die of the COVID 19 among all infected individuals regardless of whether the infected show symptoms of the disease or not. As with many diseases, IFR is not always equivalent to the number of reported deaths caused by COVID 19. This is because a significant number of deaths, although caused by COVID 19, will not be recognized as deaths caused by COVID 19 [27]. Also, there are many asymptomatic cases of infection which are never detected [28, 29, 30].
Due the fact that there is a lag in time between when people are infected and when they die, patients who die on any given day were infected much earlier, and thus the denominator of the mortality rate should be the total number of patients infected at the same time as those who died [27]. David et al. estimated mortality rate by dividing of deaths on a given date by the total number of persons confirmed as COVID 19 cases 14 days before [27]. It is based on the assumption that maximum incubation period is14 days [34]. However, if we take into account that the number of registered cases of COVID 19 infection is usually significantly lower than the actual number, assuming that the data on deaths are accurate, the real IFR value is significantly lower than the calculated value [26]. If we apply this to the situation in Serbia, the daily value of IFR on July 10th, 2020, when the largest number of deaths was registered in one day, was 9.33%, considering that 18 deaths were registered on July 10th and 14 days earlier 193 confirmed cases of COVID 19 infection. The raw values of IFR for the period between March 6th and August 10th were as follows: median of IFR = 2.11%, and average value of IFR = 7.15% bounded in interval 4.17%-10.13%. When we compared these values with those published by the WHO, CDC and other authors [29, 30, 31] we concluded that they differ significantly. Considering these findings, the IFR values adopted in the model (for various population groups and genders) were primarily taken from the literature, with a remark that the selection of IFR values was based on preliminary comparison of the overall Serbian IFR with similar IFRs found in the literature, taking into account the registered deaths and most probable number of infected individuals [8]. To make this possible, the first step was to correct the local raw IFR value mentioned above. Based on real data we first calculated the population at risk of dying from COVID 19 infection for each individual day since the outbreak, ending on August 10th, 2020. The number at risk on a given day should correspond to the number of deaths from COVID 19 infection, considering the lag period from infection to death. For this calculation, we used the data on the number of deaths Dt in Serbia registered on daily bases [2]. We hypothesized that the distribution of time periods tn from the moment of COVID 19 infection to death follows the lognormal distribution defined by the parameters m = 26.8 and σ = 12.4 days [32]. Taking into account the published data on the expected percentages of asymptomatic cases in the population, we assumed that the registered number of cases of infection is 20% of the total number of actually infected [29, 33, 34]. Based on the formula: by reverse, we calculated the population at risk of dying from COVID 19 infection for each individual day, where mt is the probability that the time between infection and death is t days and follows the lognormal distribution (m = 26.8, σ = 12.4) [35]. After that, the daily IFR values were calculated according to the formula IFR(t) = Dt /Ir(t) [35]. Based on IFR values calculated in this way, we made descriptive statistics and obtained the mean value of IFR = 0.70%, bounded in the interval 0.46-0.94% and a median of 0.19%. It is important to note here that this value corresponds with the COVID 19 IFR values found in Eastern European countries and Spain [35, 36]. Taking these findings into account, we decided to take the IFR values specific to certain population strata recorded in Spain as the most appropriate for our case [36]. The adopted IFR values are listed in Table 3.
The model sensitivity analysis was conducted by changing the most sensitive model parameters: β, f, r. The values of these parameters were increased by 5%, 10% and 25% relative to the base scenario and changes in output indicators were observed.
3. Results
3.1. Predicting the number of sick, hospitalized patients and deaths caused by COVID 19 in the absence of any intervention measures
After the simulation, the model demonstrated that with Ro = 2.46, and without the implementation of any anti-epidemic measures, the initial doubling time of the infection was five days. The epidemic wave would peak 158 days after the outbreak, and it would yield 157,424 infected individuals in a day. Afterwards, the infection rate would decline for 180 days, eventually reaching the daily incidence of 225 newly infected, after which the next epidemic wave would start. The second wave would peak 569 days after the onset of epidemic and yielding 36,922 infected individuals in a day. The third epidemic wave would peak 914 days later, with 18,791 infected individuals in a day. The true cumulative incidence in the first year of the epidemic would be 6,633,834 infected people with SARSCov-2 virus, while the apparent cumulative incidence would be 1,326,767 infected. A total of 5,525 patients would die due to COVID19 consequences. With Ro = 3.1, the following results were obtained: the initial doubling time of the infection was five days, true cumulative incidence 7,035,433, apparent cumulative incidence 1,407,089, and the total deaths of 5,868. Fig. 1, panels a) and b) show daily variations in the number of susceptible, latently infected, infected and recovered patients, at basic reproduction numbers of Ro=2.46 and Ro=3, respectively. Panels c) and d) of the same figure show daily fluctuations in susceptible, recovered and net reproduction rates Rn for Ro=2.46 and Ro=3, respectively. Panels e) and f) of Fig. 1 show daily variations in Rn, true and apparent disease incidences at basic reproduction numbers of Ro=2.46 and Ro=3. The shaded area corresponds to the period when the daily number of new COVID19 infected individuals increasing, and therefore all of the following hold: Rn>1, proportion susceptible >1/Ro and the proportion of population that is recovered (immune) is below the herd immunity threshold. Fig. 2 panels a) and b) show a prediction of necessary hospital capacities. Panels c) and d) of Fig. 2 show predicted numbers of sick and dead due to COVID19 at Ro=2.46 and Ro=3 and age structure of hospitalized patients and deaths.
Model prediction of latently infected, diseased, recovered and daily fluctuations of Rn.
Model prediction of required hospital capacities under the assumption of different intervention measures.
3.2. Predicting the number of sick, hospitalized patients and deaths caused by COVID 19 in the conditions of application of restrictive anti-epidemic measures
When the spread of COVID19 epidemic through totally susceptible population in the Republic of Serbia is simulated, under an assumption of only incidental movement among the population, basic reproduction number of Ro = 2.46, and with lock-down of entire country, a significant slowdown of the epidemic was observed. Initial infection doubling time was 11 days. The peak of the epidemic wave was reached 269 days after the epidemic onset and it was 25,754 infected in one day. In the first year of the epidemic 3,283,691 people would get infected and 3,391 people would die. When the basic reproduction number was increased to Ro=3.1, while keeping the same conditions of locking down the entire country, results changed significantly. The initial infection doubling time was 9 days, the epidemic wave peak was reached after 164 days and it yield 65,327 infected people in one day. In the first year of the epidemic 5,021,561 people would get infected and 4,923 would die. Table 6 provides the overview of epidemic indicators obtained from the simulations of all five scenarios with Ro = 2.46 and Ro = 3.1. Figures 3. and 4 provide a comparative overview of results of all five simulated scenarios. Panels a) and b) of Fig. 3 show the values of cumulative incidences on a daily basis for Ro=2.46 and Ro=3.1, respectively. Panels b) and c) of Fig. 3 show the expected total number of hospitalised patients and patients in intensive care units (ICU) for Ro=2.46 and Ro=3.1, respectively. Fig. 4 provides overview of required hospital capacities e.g. hospital bed occupancy and the occupancy of beds in ICU on daily bases for Ro=2.46 and Ro=3.1, respectively. The results show that after applying various measures to slow down the circulation of SARSCov-2, the number of newly infected people, hospitalized patients and the occupancy of hospital capacities are the lowest in the situations where rigorous anti-epidemic measures are applied to all population strata (Scenario 2 in Table 6 and in Fig. 3 and Fig. 4). Openings of pre-school and elementary school’s facilities leads to a visible jump in the number of infected and hospitalized in all strata. This finding clearly shows that children, although least susceptible to developing more severe clinical pictures, are important when transmitting SARSCov-2 (Scenario 3 in Table 6 and Fig. 3). Opening of the high schools and colleges also leads to a visible increase in the number of newly infected and hospitalized patients, including an increase in the number of deaths (Scenario 5 in Table 6 and Fig. 3). Without the application of any intervention measures, the greatest burden on the health system could be expected 196 days from the beginning of the epidemic at Ro=2.46 or 144 days at Ro=3.1. Depending on the Ro value used in simulation, it would be necessary to provide 12,740 (19,380) hospital beds and an additional 8,486 (12,355) in intensive care units. In the case of Scenario 2, there is a significant slowdown in the epidemic, so that after 621 (238) days at the moment the highest occupancy of hospital facilities would be to provide 461 (4,181) beds in COVID hospitals and 358 (3,120) beds in intensive care units (Fig. 4).
Results of different simulated scenarios (Ro = 2.46 and Ro = 3.1). The data refers to the period of 365 days from epidemic onset.
Model prediction of expected number of hospitalised patient and patient in intensive care.
Model prediction of required hospital capacities needed to treat patients with COVID 19.
3.3. COVID19 simulation and disease control by implementing a hypothetical vaccine
Assuming that the disease is spreading at the basic reproduction number of Ro=2.46, the herd immunity threshold (when the disease can be expected to slow down and the chain of infection is expected to be broken) would be 59.35%, while at Ro=3.1, it would be 67.74%. Depending on the efficacy of the potential vaccine, the required vaccination coverage should be 87% (ve =68%), 74.19% (ve =80%), 69.82% (ve=85%) and 95.41% (ve=71%), 84.68% (ve=80%), 79.70% (ve=85%) for Ro=2.46 and Ro=3.1, respectively. Fig. 5 shows different scenarios of COVID19 control strategies based solely on vaccination.
Results of simulated COVID19 control based solely on vaccination, scenarios 6-9 (Ro=2.46 and Ro=3.1).
3.4. Results of the model sensitivity analysis
The results of the sensitivity analysis are presented in table 9. The table shows increased values of input parameters and the percentage of the parameter increase relative to the basic scenario, as well as the values of output results obtained after the simulation of the changed scenario.
Results of model sensitivity analysis.
4. Conclusions and discussion
For the needs of this research, we augmented the classic deterministic model by adding the compartments of vaccinated and latently infected subjects. By adding birth and death rates, we enabled daily fluctuations of the overall simulated population, which brought us closer to the real conditions in which the disease is transmitted. When we assumed that the recovered lose immunity over time, we obtained dynamic oscillations of epidemic waves through susceptible population. When we simulated different disease control scenarios of COVID19 epidemic based on non-pharmaceutical intervention measures, scenario number 2 proved to be the most effective approach to the disease control, because it implemented the most comprehensive anti-epidemic measures (entire country lock down). However, the basic problem of this approach is feasibility and possibility to maintain the measures in the long term.
The model demonstrated that students, children and younger school-age generations have an important role in transmitting COVID19, especially if they come into contact with a more vulnerable population. The model showed that, in the case of returning school children of all ages to schools, an increase of 32.98% of the estimated deaths and 57.50% of the number of infected is possible, when compared to the conditions before opening of the schools (Scenarios 2 and 4). However, most dead and seriously diseased people are found in the older population. This is particularly important when planning intervention measures, especially when deciding on which restrictions to be lifted and how (opening schools, students’ return to faculties etc.). The model demonstrated that COVID19 has a potential to spread rapidly and linger in population. Due to a large number of the infected and duration of the disease, there are significant needs for hospital capacities, especially in the conditions when the disease is suppressed by implementing partly relaxed anti-epidemic measures, or in the case of absence of any measures. Without the application of any intervention measures, at the moment of the greatest load, depending of actual Ro, the health system would have to provide 12,740 (Ro=2.46) hospital beds for the care of the patients and an additional 8,486 in intensive care units. On the other hand, in the case of the application of the strictest anti-epidemic measures, the needs decrease to only 461 beds in COVID hospitals and additional 358 beds in intensive care units. In the case of continued implementation of current measures, which are significantly less intense than the measures applied at the beginning of the epidemic (Scenario 4), it is necessary to provide 3,508 beds in COVID hospitals and 2,617 beds in intensive care units in the entire country. The model also shows that the needs for hospital capacities decline with the ending of the first epidemic wave, since daily incidence decreases and during the second and third waves it never reaches the initial peaks, but these needs still remain substantial. For example, in the case of Scenario 1, at the top of the second epidemic wave, it is necessary to provide 3,643 beds in COVID hospitals and 2,575 in intensive care units, which makes 30% and 25% of the required capacities of the first wave of the epidemic.
Based on the cyclical patterns of the epidemic waves and duration of simulated epidemics, the model demonstrated that the disease has a potential to linger in population and that it will most probably have a seasonal pattern. Therefore, potential vaccines can have an enormous potential and significance for COVID19 control. Depending on the efficacy of future vaccines, the disease can be stopped and curbed almost solely by implementing the measure of vaccination. However, the necessary conditions for these predictions and expectations are the efficacy of potential vaccines and the ability of a health systems to implement vaccination to a satisfactory extent and rapidly, especially with regards to the most sensitive categories of population. Depending on the Ro, a vaccine that would have an efficacy ≥ 68-71% could stop the pandemic and break the chain of infection. However, even vaccines with lower efficacy could be useful as they would significantly reduce the number of cases and deaths, especially if used in combination with the other disease control measures. The required minimal vaccination coverage should be 87% (ve =68%), 74.19% (ve =80%), 69.82% (ve=85%) and 95.41% (ve=71%), 84.68% (ve=80%), 79.70% (ve=85%) for Ro=2.46 and Ro=3.1, respectively. The minimum daily vaccination rate should be 0.47% for vaccines with an efficiency of 85%, and 0.59% for vaccines with an efficiency of 68%.
Based on the obtained results, we can conclude that at this point, without the application of specific pharmaceutical products, COVID19 suppression is highly dependent on the basic reproduction number (Ro), and that more intensive contacts and relaxed measures can result in a dramatic spread of the virus. The choice of intervention measures depends on the feasibility of their implementation and their efficacy in different social contexts. COVID19 will likely have to be suppressed in this way for a certain period of time. This will most probably last until sufficient quantities of a reliable and effective vaccine are available, and thereafter until optimal vaccination coverage is achieved.
Data Availability
The data used to support the findings of this study are included within the article and can be used without restriction.
Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.