Simulation and prediction of spread of COVID-19 in The Republic of Serbia by SEIRDS model of disease transmission.

As a response to the pandemic caused by SARS-Cov-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 7 February, 2021, a total of 406,352 cases of SARSCov-2 infection have been reported in Serbia, 4,112 deaths caused by COVID-19. In order to better understand the epidemic dynamics and predict possible outcomes, we have developed an adaptive mathematical model SEAIHRDS (S-susceptible, E-exposed, A-asymptomatic, I-infected, H-hospitalized, R-recovered, D-dead due to COVID-19 infection, S-susceptible). 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 COVID-19, the model is extended 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 COVID-19 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.

Introduction

On 11 March 2020, the World Health Organization characterised the disease caused by the novel SARS-Cov-2 virus as a pandemic (World Health Organization (WHO) 2020). The 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 coronavirus 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 (https://covid19.rs 2020).

Due to the absence of specific pharmaceutical intervention, Serbia, like other countries, implemented an anti-epidemic strategy based on physical 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 (https://covid19.rs 2020). Although the return of extensive measures has yielded favourable results, the further development of the epidemic is not clear.

For these reasons, mathematical modeling has a crucial role in understanding the epidemic and predicting possible outcomes. modeling 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 (Wangping et al., 2020; Kucharski et al., 2020; Keeling and Rohani, 2021). Mathematical modeling 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 (Kermack and McKendrick, 1927; Barnes and Fulford, 2021). Subsequent models, became more complex and adapted to the needs of research (Keeling and Rohani, 2021).

Since the outbreak of COVID-19, many published papers have dealt with the implementation of mathematical modeling and prediction of possible outcomes of COVID-19 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 COVID-19 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 (Imperial College COVID-19 Response Team 2020). 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 R, which we use in this paper as one of the key parameters (Wu and Googan, 2020). 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 (Vynnycky and Richard, 2010). However, such assumptions are not objective in the case of COVID-19 spreading and they are the basic problem in the modeling of this pandemic. In reality, the human population fluctuates constantly (2021). In order to account for such fluctuation, and better understand the COVID-19 epidemic in the Republic of Serbia, we have employed mathematical modeling 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 COVID-19.

Unlike the classic SIR model, the SEAIHRDS (S-susceptible, E-exposed, A-asymptomatic, I-infected, H-hospitalized, R-recovered, d-dead due to COVID-19 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 hospitalized 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.

Methodology

This section presents the research methodology and the proposed model, which were 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.

SEAIHRDS 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 I. It depends on the number of infectious individuals (I) 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 as follows:

The change of rates in every compartment per unit time in SEIRDS model is presented in the following series of differential equations: where ƒ is rate of onset of infectiousness 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 rate and death rate unrelated to COVID-19.

However, considering that implemented anti-epidemic measures against COVID-19 do not have an identical impact on the population's age subgroups and that COVID-19 pathogenesis varies in different age subgroups, we propose the use of multi-compartment version of standard SEIRDS model. The model, named SEAIHRDS, monitors the dynamics of following compartments: susceptible individuals (S), latently infected with SARSCov-2 (E) (exposed to/presymptomatic), asymptomatic infectious individuals (A), infectious individuals with symptoms/clinically ill (I), hospitalized patients (H), recovered individuals (R), and those who died due to disease (D). In this model the 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 (S), elementary school children (S), high school children (S), college students (S cs), unemployed population (S), employed population (S), and elderly/retired (S) (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 20% to 75%, taking into account the realistic possibilities of maintaining a minimum work process, functioning of the society and feasibility of such measures.

Table 1

Structure of different population strata in the Republic of Serbia (Chan et al., 2013).

Stratum	Population	Percentage of total population
Younger than 7 years	356,377	5.10%
Elementary school	550,527	7.88%
Secondary school	249,455	3.57%
Students	241,698	3.46%
Employed	2197,065	31.46%
Pensioners	1715,152	24.56%
Others	1672,330	23.95%

Given that intervention measures, applied in response to the emergence of COVID-19, 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 λ is force of infection in the i population strata, β is the rate at which susceptible persons in the i population strata and infectious persons in j population strata come into effective contact per unit of time, and I is the number of infectious persons in j 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: where α is the proportion of asymptomatic cases, η accounts for the relative infectiousness of asymptomatic carriers (in comparison to symptomatic carriers), r is the rate at which infectious individuals whit symptoms are recovered, γ is the rate at which asymptomatic infectious individuals are recovered, is rate at which infectious individuals are hospitalized, ε is the rate at which hospitalized patients are recovered, δ is the rate at which infected individuals die from COVID-19 infection and ω is rate of waning of immunity (Supplementary material).

Determining the herd immunity threshold and control of COVID-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 achieve herd immunity can be obtained by multiplying herd immunity threshold with reciprocal value of vaccine efficacy, v

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 (European Centre for Disease Prevention and Control 2020). Specific IgM and IgG antibodies against SARS-CoV-2 develop 6 to 15 days after the onset of the disease (Woelfel et al., 2020; Liu et al., 2020; Long et al., 2021; J Zhao et al., 2020; Okba et al., 2020). 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 (J Zhao et al., 2020). Although duration of the immune response against CVOVID 19 is still unknown, comparing with other coronaviruses, where immunity wane within 12 to 52 weeks after the first symptoms appear (Kellam and Barclay, 2020; Kissler et al., 2020), 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 (The World Bank 2021), we assumed that durable immunity against COVID-19 is possible (Kissler et al., 2020; Callow et al., 1990). Immunity to HCoVOC43 and HCoV-HKU1 appears to wane appreciably within 1 year (Callow et al., 1990; Dan et al., 2021), whereas SARS-CoV-1 infection can induce longer-lasting immunity (Chan et al., 2013). S. F. Lumley et. in a study conducted on 1246 persons recovered from COVID-19 found no symptomatic re-infections over 6 months (Lumley et al., 2021).

Based on these findings, and fact that SARS-Cov-2 virus is also beta coronavirus, we assumed that in the event of the development of a successful vaccine, immunity against the SARS-Cov-2 virus could last for a year, as well as after recovery after a natural infection. In a study conducted to determine the dynamics of SARS-Cov-2 transmission in the post pandemic period, Kissler et al. applied a similar approach in defining the possible length of the immunity period (Kissler et al., 2020).

For the purpose of modeling a control strategy based on vaccination, 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, ѱ, is ѱ required to interrupt the infection (Keeling and Rohani, 2021). The basic reproduction number under the vaccination is R The proportion of effectively protected persons, p, is conditioned by parameter the vaccine efficacy, v. 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 v (Keeling and Rohani, 2021). 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 follows:

The compartment S(t) is slightly modified as follows:

The other of compartments of SERIDS model remain unchanged.

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 (Vynnycky and Richard, 2010; Lumley et al., 2021; Thrusfield, 2007). 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 (Vynnycky and Richard, 2010; Chowell et al., 2009; Thrusfield, 2007). 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 (Marco D'Arienzo and Coniglio, 2020). 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 ρ. In this way, we obtained the per capita contact rate specific for each separate stratum based on following formula

The values of the ρ 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 (Demografski pregled 2017; The World Bank 2021). The latter study reported that the life expectancy in Serbia was 79.06 years in 2017 (The World Bank 2021). 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 (ƒ −1) and the average period during which an infected person is shading the SARS-Cov-2 virus (T) were adopted from the relevant literature as ƒ −1 = 3.5 days (Byrne et al., 2020; Shujuan Ma et al., 2020) and T = 9.3 days, respectively (Byrne et al., 2020). Also the data on the percentage of hospitalized patients and those whose therapy requires intensive care, used for prediction of required hospital capacities, as shown in Table 2, were taken from literature (Imperial College COVID-19 Response Team 2020; Verity et al., 2020; Byrne et al., 2020).

Table 2

Age structure of the population of the Republic of Serbia and expected percentage of hospitalized patients, patients in intensive care, and death rate caused by COVID-19.

Population age groups	*Population	Percentage of total population	**Expected% of hospitalized patients (σ)	***Expected% of patients whose treatment requires intensive care	****Infection fatality rate IFR (male/female)
0–9	458,199	6.56%	0.00%	5.00%	0.04%;0.01%
10–19	445,481	6.38%	0.04%	5.00%	0.00%;0.02%
20–29	1028,226	14.73%	1.04%	5.00%	0.00%;0.01%
30–39	951,615	13.63%	3.43%	5.00%	0.00%;0.05%
40–49	968,854	13.88%	4.25%	6.30%	0.08%;0.04%
50–59	963,229	13.79%	8.16%	12.20%	0.33%;0.20%
60–69	815,244	11.68%	11.80%	27.40%	1.62%;0.62%
70–79	696,045	9.97%	16.60%	43.20%	6.11%;2.68%
80-	655,711	9.39%	18.40%	70.90%	16.40%;6.49%

*[28], **[32], ***[8], ***[33].

Parameters such as δ and r are related to the infectious fatality rate (IFR) for COVID-19 and average times taken from onset of symptoms to death (T) or recovery (T). These parameters were calculated using the following formulas:

The IFRs, shown in Table 2 , 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 .

Table 3

SEAIHRDS model parameters.

Input parameters	Mark	Value	Source
Population	Nt0	6982,604	(Demografski pregled 2017)
Initial number of cases	It0	1	Assumed
Initially immune	Imm0	0	Assumed
Basic reproduction number	R0	2.46 (3.1)	(Marco D'Arienzo and Coniglio, 2020)
Effective contact rate	Ce	0.38	Estimated
Per capita contact rate	β	0.0000000378	Estimated
Daily infection rate (transfer E→I)	Ƒ	0.294118	Estimated
Recovery rate of symptomatic cases	R	0.107527	Estimated
Daily rate of waning of immunity	Ω	0.002739726	Estimated
Per capita birth rate	B	0.000025205	(Demografski pregled 2017)
Per capita death rate unrelated to COVID-19	M	0.000036006	Estimated
Life expectancy in years	L	76.09	(The World Bank 2021)
Duration of latent infection in days	ƒ−1	3.5	(Byrne et al., 2020)
Duration of infectious period in days (clinical cases)	TR	9.3	(Byrne et al., 2020)
Duration of immunity in days	Imm	365.00	Assumed
Incubation period in days	Inc	5.8	(McAloon et al., 2020)
Time period (day)	Dt	1.00	–
Average treatment duration in hospital	Days	15,9	(Boëlle et al., 2020)
Average time spent in intensive care	Days	27	(Boëlle et al., 2020)
Recovery rate of hospitalized cases	ε	0.062893	Estimated
Average times taken from onset of symptoms to death	Days	17	(Linton et al., 2020)
Infectious period for asymptomatic cases	Days	7.25	(Byrne et al., 2020; S Ma et al., 2020)
Recovery rate of asymptomatic cases	γ	0.137931	Estimated
Expected percentage of asymptomatic cases	–	30%	(Johansson et al., 2021)
Infectiousness of asymptomatic cases in relation to symptomatic cases
	Η	75%	(Johansson et al., 2021)

Setting disease control scenarios

Five different scenarios were developed for simulating the COVID-19 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 .

Table 4

Description of different simulated non-pharmaceutical intervention scenarios.

Mark	Scenario	Scenario description
SC1	Base-case scenario	The population relies on development of herd immunity. No anti-epidemic measures are implemented.
SC2	Lock down of the entire country	Pre-schools, schools, and colleges are fully closed – reduction in contacts at educational institutions by 75%; reduction in contacts in workplaces by 50%; reduction in contacts of the elderly (older than 65) by 50% at R0 = 2.46 or by 65% at R0 = 3.10; physical distancing of the unemployed population and in public places – reduction in contacts by 45% at R0 = 2.46) or by 55% at R0 = 3.10.
SC3	Partial lockdown of the country - I	Elementary and pre-school educational institutions are open. High schools and colleges are closed. Reduction in contacts by 75% at colleges and high-schools; reduction in contacts in workplaces by 50%; reduction in contacts of the elderly (older than 65) by 60% at R0 = 2.46 or by 65% at R0 = 3.10; social distancing of the unemployed population and in public places – reduction in contacts by 40% at R0 = 2.46 or by 55% at R0 = 3.10.
SC4	Partial lockdown of the country - II	Colleges are closed – reduction in contacts by 75%; reduction in contacts in the workplace by 50%; reduction in contacts of the elderly (older than 65) by 60% at R0 = 2.46 or by 65% at R0 = 3.10; social distancing of the unemployed population and in public places – reduction in contacts by 40% at R0 = 2.46 or by 55% at R0 = 3.10.
SC5	Partial lockdown of the country - III	Reduction in contacts in the workplace by 50%; reduction in contacts of the elderly (older than 65) by 60% at R0 = 2.46 or by 65% at R0 = 3.10; social distancing of the unemployed population and in public places – reduction in contacts by 40% at R0 = 2.46 or by 55% at R0 = 3.10.

The timing of the simulation of anti-epidemic measures, i.e. the reduction of individual contacts corresponds to the actual date when the implementation of measures in the real epidemic in Serbia began (March 15, 2020). Considering that it is not realistic to expect the desired level of reduction of physical contacts to be achieved in one day, in each scenario, a period of gradual introduction of measures was simulated (lag period of 7 days). Since each individual scenario was simulated at R0 of 2.46 and 3.10, different contact reduction rates were applied accordingly.

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.

Model sensitivity analysis and calibration

Considering the world experience with detection of COVID-19 cases, as well as the unreliable data on COVID-19 infections which are currently available worldwide, model calibration is very challenging, and can result in obtaining inaccurate values for the parameters (Viguerie et al., 2021). 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 (Viguerie et al., 2021).

As part of the national seroepidemiological study, 1006 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 COVID-19 infection is more reliable for use in model calibration, e.g. infection fatality rate. Alex et al. reached a similar conclusion when simulating COVID-19 by using the SEIRD model with heterogeneous diffusion (Viguerie et al., 2021). 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.

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 (Baud et al., 2020). Also, there are many asymptomatic cases of infection which are never detected (Manski and Molinari, 2021; WHO 2020; CDC 2020; BMJ 2020).

However, according to new findings, the overall estimate of the proportion of people who become infected with SARSCoV-2 and remain asymptomatic throughout infection was 20% (95% confidence interval) with a prediction interval of 3%–67% in 79 studies that addressed this question (Buitrago-Garcia et al., 2020; Meyerowitz et al., 2020). Michael A. Johansson et al. reported that 30% of individuals with infection never develop symptoms and are 75% as infectious as those who do develop symptoms, and concluded that persons with infection who never develop symptoms may account for approximately 24% of all transmission (Johansson et al., 2021).

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 (Baud et al., 2020). 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 (The World Bank 2021; Baud et al., 2020). It is based on the assumption that maximum incubation period is14 days (Marco D'Arienzo and Coniglio, 2020). 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 (Viguerie et al., 2021). 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 (WHO 2020; CDC 2020; CEBM- 2020) 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 (Pastor-Barriuso et al., 2020). 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 D in Serbia registered on daily bases (https://covid19.rs 2020). We hypothesized that the distribution of time periods t from the moment of COVID19 infection to death follows the lognormal distribution defined by the parameters m = 26.8 and σ = 12.4 days (CEBM- 2020; Wood, 2021).

Based on the formula:by reverse, we calculated the population at risk of dying from COVID-19 infection for each individual day, where m is the probability that the time between infection and death is t days and follows the lognormal distribution (m = 26.8, σ = 12.4) (CEBM- 2020; Wood, 2021). After that, the daily IFR values ​​were calculated according to the formula IFR(t) = Dt /I (Ghisolfi et al., 2020). 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 (Pastor-Barriuso et al., 2020; Ghisolfi et al., 2020). 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 (Pastor-Barriuso et al., 2020). The adopted IFR values are listed in Table 3.

In this section, we used sensitivity analysis to estimate the amount of change on outcomes when varying the input values used in the model. Sensitivity analysis is one of the methods most frequently used for the evaluation of disease spread models (Martınez-Lo´pez et al., 2012; Karsten et al., 2005). A sensitivity analysis is carried out to characterize the impact of uncertainty of input parameters on model outputs. Sensitivity analysis consists of assessing the impact that changes in input parameters have on model outcomes. We evaluate two aspects of the model: the global behavior of the model when perturbing a group of key parameters together, and the impact of changes when perturbing the key parameters used in the model individually.

The model sensitivity analysis was conducted by changing the most sensitive model parameters: 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.

Finaly, validation of the model was performed by comparing the historical data of the real epidemic in Serbia with the data obtained with the SEAIHRDS model. For validation purposes, the current epidemic of COVID-19 was simulated, along with actual anti-epidemic measures. Taking into account the risks described above related to the accuracy of historical data from the COVID-19 epidemic in Serbia and significantly higher confidence in the accuracy of data related to the number of deaths from COVIDA-19 compared to data on the number of infected, for the purpose of model validation, only data on the cumulative number of deaths were used. The reason for this assumption is that there is still uncertainty about the proportion of the infected population that is not reported due to a mild form of the disease or the patients are asymptomatic.

In statistical analysis the coefficient of determination, R2, was used to check goodness of fit of SEAIHRDS model with COVID-19 data recorded during the real epidemic. The regression coefficient compares predicted values (y) against actual data (x). To address model uncertainties, bias, mean absolute error (MAE), mean square error (MSE), the root mean square error (RMSE), maximum deviation (MaxDev) and normalized root mean square error (NRMSE) were also estimated (Supplementary material).

Results

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 predicted that with Ro = 2.46, and without the implementation of any anti-epidemic measures, the initial doubling time of the infection could be five days. The epidemic wave could peak 219 days after the outbreak, and it could yield 99,819 infected individuals in a day. Afterwards, the infection rate could decline for 215 days, eventually reaching the daily incidence of 492 newly infected, after which the next epidemic wave could start. The second wave could peak 706 days after the onset of epidemic and yielding 25,232 infected individuals in a day. The third epidemic wave could peak 429 days later, with 13,709 infected individuals in a day. The true cumulative incidence in the first year of the epidemic could be 6229,144 infected people with SARS-Cov-2 virus, while the apparent cumulative incidence could be 4360,401 infected. A total of 20,894 patients could die due to COVID-19 consequences.

With Ro = 3.1, the following results were obtained: the initial doubling time of the infection was five days, true cumulative incidence 7133,221, apparent cumulative incidence 4993,254, and the total deaths of 23,951. 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.1, 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.1. The shaded area corresponds to the period when the daily number of new COVID-19 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 COVID-19 at Ro=2.46 and Ro=3.1 and age structure of hospitalized patients and deaths.

Fig. 1

Model prediction of latently infected, asymptomatic infectious individuals, infected, recovered and daily fluctuations of Rn.

Fig. 2

Model prediction of required hospital capacities under the assumption of different intervention measures.

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 COVID-19 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 6 days. The peak of the epidemic wave could occur 702 days after the epidemic onset, when there could be 2848 infected in one day. In the first year of the epidemic 308,581 people could be infected and 1031 people could die.

When the basic reproductive number was increased to R0 = 3.1, and for certain segments of the population the contact reduction factor ρ increased compared to the values from the scenario with R0 = 2.46 (reduction of contacts in public places and contacts of persons over 65 years by 55% and 65%, respectively), the results changed significantly. The model predicted that with R0 = 3.10, the initial infection doubling time could be 5 days, the peak of epidemic wave could occur after 246 days and it yield 8110 infected people in one day. In the first year of the epidemic 2219,251 people could be infected and 7194 could die.

Table 5 provides the overview of epidemic indicators obtained from the simulations of all five scenarios with Ro = 2.46 and Ro = 3.1. Figs. 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 hospitalized 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 SARS-Cov-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 5 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 SARS-Cov-2 (Scenario 3 in Table 5 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 5 and Fig. 3). Without the application of any intervention measures, the greatest burden on the health system could be expected 228 days from the beginning of the epidemic at Ro=2.46 or 168 days at Ro=3.1.

Table 5

Results of different simulated scenarios (Ro = 2.46 and Ro = 3.1). The data refers to the period of 365 days from epidemic onset.

Scenario mark	SC1		SC2		SC3		SC4		SC5
Basic reproductive number	Ro=2.64	Ro=3.1	Ro=2.64	Ro=3.1	Ro=2.64	Ro=3.1	Ro=2.64	Ro=3.1	Ro=2.64	Ro=3.1
Cumulative incidence (CI)	6229,144	7133,221	308,581	2219,251	1286,227	2419,079	1375,416	2489,197	1408,262	2514,936
Apparent CI	4360,401	4993,254	216,007	1553,476	900,359	1693,355	962,791	1742,438	985,783	1760,456
Overall hospitalized	278,781	320,567	13,970	98,337	56,631	105,930	60,319	108,554	61,685	109,538
Overall in intensive care	85,633	98,333	4260	29,753	17,085	32,058	18,202	32,855	18,612	33,147
Overall deaths	20,894	23,951	1031	7194	4113	7754	4383	7948	4483	8018

Fig. 3

Model prediction of expected number of hospitalized patient and patient in intensive care.

Fig. 4

Model prediction of required hospital capacities needed to treat patients with COVID-19.

Depending on the Ro value used in simulation, it could be necessary to provide 42,351 (61,739) hospital beds and an additional 20,173 (28,740) in intensive care units. In the case of Scenario 2, there is a significant slowdown in the epidemic. According to the predictions obtained by the simulation of Scenario 2, after 713 (259) days at the moment of maximum occupancy of hospital capacities, it might be necessary to provide 1387 (3940) beds in COVID-19 hospitals and 784 (2162) beds in intensive care units (Fig. 4).

COVID-19 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% (v =68%), 74.19% (v =80%), 69.82% (v=85%), 65.94% (v=90%) and 95.41% (v=71%), 84.68% (v=80%), 79.70% (v=85%), 75.27% (v=90%) for Ro=2.46 and Ro=3.1, respectively. The different ways of including vaccination in the SEAIHRDS model are detailed in the supplementary material. Fig. 5 shows different scenarios of COVID-19 control strategies based solely on vaccination.

Fig. 5

Results of simulated COVID-19 control based solely on vaccination, scenarios 6–9 (Ro=2.46 and Ro=3.1).

Results of the model sensitivity analysis

The results of the sensitivity analysis are presented in Tables 6 and 7 . The tables show 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.

Table 6

Results of the model sensitivity analysis of individual parameters used in the model.

Model parameter mark	Change relative to the basic scenario	CI	Deaths (Dth)	Change in CI relative to the basic scenario	Change in Dth relative to the basic scenario
R0	5%	6398,486	21,470	2.72%	2.75%
R0	10%	6553,978	21,998	5.21%	5.28%
R0	25%	6982,604	23,782	12.10%	13.82%
ƒ−1	5%	6224,056	19,885	0.08%	4.83%
ƒ−1	10%	6219,502	18,968	0.15%	9.22%
ƒ−1	25%	6208,481	16,666	0.33%	20.24%
R	5%	6068,299	20,351	2.58%	2.60%
R	10%	5908,223	19,811	5.15%	5.18%
R	25%	5425,721	18,182	12.90%	12.98%

Table 7

Results of the model sensitivity analysis of group of parameters using a perturbation up to 25%.

Model parameter mark	Change relative to the basic scenario	CI	Deaths (Dth)	Change in CI relative to the basic scenario	Change in Dth relative to the basic scenario
R0, ƒ−−1, r	5%	6242,433	19,948	−0.21%	4.53%
R0, ƒ−1, r	10%	6257,685	19,093	−0.46%	8.62%
R0, ƒ−1, r	25%	6321,339	16,984	−1.48%	18.71%
R0, β, ƒ−1, r, b, m, γ, ε, η, ω, τ̅d	25%	6147,725	13,203	1.31%	36.81%

Results of the model validation

The model was validated on the most recent historical data (from January 1, 2021 to February 1, 2021). As shown in Table 8 . all measures of the prediction quality of deceased due to COVID-19 are low. The average difference between the actual number of people who died of COVID-19 and predicted one is only 2.05% and the maximum deviation between the predicted and actual number will not exceed 4.82% with a probability of 95%.

Table 8

Measures of the prediction quality.

	MAE	%Error	MSE	RMSE	Normalized MAE	Normalized MSE	Max Deviation
Deceased	73	2.04%	7212.20	84.92	2.05%	0.06%	4.82%

Based on Fig. 6 and Table 9 , we can conclude that the SEAIHRDS model fits historical data quite well. Pearson's r and coefficient of determination (R2) have shown strong the linear relationship between real deceased and the number of deaths predicted by SEAIHRDS model.

Fig. 6

The observed number of deceased individuals (blue), number of deceased individuals modeled with SEAIHRDS model (orange), and predicted number of deceased individuals (green) by model and corrected with real data (95% confidence interval between dotted lines).

Table 9

Regression statistics.

Multiple R	R Square	Adjusted R Square	Standard Error	Observations
0.981678289	0.963692263	0.963124954	86.24143528	66

More information on validation results are detailed in the supplementary material.

Conclusions and discussion

For the needs of this research, we augmented the classic deterministic model by adding the compartments of vaccinated, asymptomatic, hospitalized 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.

The input values for the parameters used to simulate the COVID-19 were obtained either from literature review or were calculated on the basis of data have taken from the literature and other official sources. Whereas some of these inputs are well documented; other input values are either not so well documented or based on potential subjective opinions (i.e., expert opinion, historical data from epidemic etc.). In any case, well documented or not input values have potential to impact results and, therefore, should be carefully evaluated. The results of conducted sensitivity analysis show that the sources of uncertainty are different for each output considered and it is necessary to consider multiple output variables for a proper assessment of the model. The most influential parameter is r than R0 and ƒ .

The model was tested on Serbian COVID-19 statistic data and obtained validation results allow us to conclude that the proposed model has good prediction ability and performance. However, although we obtained satisfactory results during the validation of the model, worth noticing also that some of the model parameters were estimated based on available data that might be less precise due to the difficulty of being measured. That could be the reason why the values of some parameters e.g. recovery rate, contact rate, daily infection rate, that are estimated in hospitals may differ from those acquired by this study. For that reasons SEAIHRDS model can be used for the long-term rough predictions of the epidemic. Obtained long-term predictions reflect the general dynamic of the outbreak and are especially useful for the healthcare system workers and governmental officials.

When we simulated different disease control scenarios of the COVID-19 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 the feasibility and practicability to maintain the measures in the long term.

The model predicted that students, children and younger school-age generations have an important role in transmitting COVID-19, 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 10.48% in the estimated deaths and 12.16% 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 predicted that COVID-19 has a potential to spread rapidly and linger in population. Due to a large number of the infected persons 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. According to the prediction, without the application of any intervention measures, at the moment of the greatest load, depending of actual Ro, the health system should provide 42,351 (Ro=2.46) hospital beds for the care of the patients and an additional 20,733 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 1387 beds in COVID-19 hospitals and additional 784 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, it is necessary to provide 3537 beds in COVID-19 hospitals and 1945 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 11,845 beds in COVID-19 hospitals and 6378 in intensive care units, which makes 27.96% and 30,76% 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 predicted 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 COVID-19 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 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. The required minimal vaccination coverage should be 87% (v =68%), 74.19% (v =80%), 69.82% (v=85%) and 95.41% (v=71%), 84.68% (v=80%), 79.70% (v=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, COVID-19 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. COVID-19 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.

Funding: This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

CRediT authorship contribution statement

Slavoljub Stanojevic: Conceptualization, Methodology, Formal analysis. Mirza Ponjavic: Validation. Slobodan Stanojevic: Investigation. Aleksandar Stevanovic: Validation. Sonja Radojicic: Supervision, Resources.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.