A Mathematical Model for COVID-19 with Variable Transmissibility and Hospitalizations: A Case Study in Paraguay

: Forecasting the dynamics of the number of cases with coronavirus disease 2019 (COVID-19) in a given population is a challenging task due to behavioural changes which occur over short periods. Planning of hospital resources and containment measures in the near term require a scenario analysis and the use of predictive models to gain insight into possible outcomes for each scenario. In this paper, we present the SEIR-H epidemiological model for the spread dynamics in a given population and the impact of COVID-19 in the local health system. It was developed as an extension of the classic SEIR model to account for required hospital resources and behavioural changes of the population in response to containment measures. Time-varying parameters such as transmissibility are estimated using Bayesian methods, based on the database of reported cases with a moving time-window strategy. The assessment of the model offers reasonable results with estimated parameters and simulations, reﬂecting the observed dynamics in Paraguay. The proposed model can be used to simulate future scenarios and possible effects of containment strategies, to guide the public institution response based on the available resources in the local health system.


Introduction
Recently, the use of mathematical models in public health decision-making has been increasing, especially to forecast the spread of epidemics [1]. Timely country-based decisions have become critical during the coronavirus disease 2019 (COVID-19) pandemic due to the high contagiousness and the non-negligible proportion of severe cases [2,3].
After a year of intensive research, several types of vaccines for COVID-19 were developed and approved, but are not widely available in most of the countries [4]. Thus, the use of mathematical models in forecasting the course of this pandemic and the required hospital resources is still relevant in order to enforce measures to maintain the situation at levels that can be sustained by a local health system.
A vast amount of works use mathematical models to study the dynamics of the spread of COVID-19 in different contexts [5][6][7][8][9][10][11][12][13][14]. The dynamics of epidemic transmission can be represented by compartmental models based on a system of ordinary differential equations [15]. In these models, the population in a region is divided into groups or compartments such as susceptible, exposed, infectious, recovered, among others, where other groups can be added and/or removed depending on the dynamics or mechanism to be modelled.
Among the simplest and most traditional models for the dynamics of transmission of COVID-19, we find the Susceptible-Infectious-Recovered (SIR) model and, the Susceptible-Exposed-Infectious-Recovered (SEIR) model. Cooper et al. [5] have presented modifications of the SIR model by considering variable the susceptible population, and Ellison [6] have added heterogeneity to the population. Further modifications of the model can be seen in the work of Crokidakis [7], who has included a "quarantined" group in order to analyse the effect of isolation of newly confirmed cases. Additionally, He et al. [8], based on the SEIR model, proposed to divide the infected group into two groups, being with and without interventions, and have also incorporated the "quarantined" and hospitalised groups.
On the other hand, Arino and Portet [11] have considered the latently infected and asymptomatic infectious groups in order to take into account the pre-symptomatic or asymptomatic transmission, which are characteristics of the COVID-19 [9,10]. Additionally, they have subdivided each group into stages to incorporate the Erlang distribution of residence times in each group [16]. There are more elaborate models, such as, SIDARTHE proposed in [12] or θ-SEIHQRD presented in [13], where more compartments are included. In metapopulation models, the spread dynamic of disease is split into age groups or regions, and these groups are interconnected to each other by additional interaction or mobility terms. For instance, Arenas et al. [14] applied a metapopulation model with age-stratified and mobility-based for the spread of COVID- 19. Although there have been published several works proposing mathematical models for the spread dynamics of COVID-19, there are still challenges in adapting particular models in specific situations. These challenges are related to: • the amount of available data and the quality of such data, that vary from one country to another, and • the containment strategies applied in each country.
In particular, highly detailed models are composed of many parameters, which have to be determined independently, using different groups of data. When the available data is limited, then the models with many parameters are over-fitted, which increases the prediction error [17]. However, the simple SIR or SEIR models with constant parameters along the epidemic cycle, fail to represent the dynamics when different containment strategies are applied along the time.
Thus, we developed a mathematical model, which is as simple as possible, but at the same time, applicable to several situations. Our model, named SEIR-H, is based on the simple SEIR model, where R stands for reported instead of recovered. The reported group is composed of individuals who have been diagnosed as positive confirmed cases. Also, in order to answer the specific question about required hospital resources, our model tracks the number of hospitalisations and admissions to an Intensive Care Unit (ICU) of COVID-19 patients. Some parameters, such as the latent period or the average residence times in each compartment or group, are maintained constants as in traditional SEIR models. A key feature that distinguishes our model is to allow the transmissibility and proportions of individuals who need hospital beds or ICU, to vary daily with a moving time-window strategy.
Most of the mathematical models developed in the context of the COVID-19 pandemic were applied to large countries or those with early outbreaks, with a large amount of available data. However, in the context of limited data availability in low-to-middle-income countries, specific models and strategies are needed to overcome this difficulty. In particular, the mathematical model developed in this work is based on limited data available from Paraguay, where except around the capital city, there is no regular testing policy in the rest of the country, nor a good record of region-based mobility. The data used in order to calibrate the model was provided by the Ministry of Public Health and Social Welfare (MSPyBS: Ministerio de Salud Pública y Bienestar Social) and Health Monitoring Centre in Paraguay (DGVS: Dirección General de Vigilancia de la Salud) (http://dgvs.mspbs.gov.py/, accessed on 20 May 2021).
The spread of COVID-19 in Paraguay is characterised by a strongly enforced initial lockdown, which kept the amount of new positive cases low for several months. In this sense, it is worth discussing the effect of early lockdown in the dynamics of spread, as well as, how the time-varying parameters can explain different scenarios observed in Paraguay. The feature of time-varying parameters allows a user to draw future scenarios and use our model as a prediction tool to guide decision-making to contain outbreaks and plan hospital resources accordingly.
The work is organised as follows: the proposed model is presented in Section 2 including the spread dynamics in Section 2.1.1 and the dynamics of hospital resources in Section 2.1.3, as well as the parameter estimations in Section 2.3; Section 3 shows the results obtained and the forecast capability of the model, and it is followed by the discussions in Section 4, the relevance of the present work in Section 5, and final remarks in Section 6.

Materials and Methods
In this section, a description of the SEIR-H model is given in Section 2.1. Section 2.2 gives a brief overview of the restriction policies applied in Paraguay. Finally, the methodology used for parameter estimation is presented in Section 2.3.

SEIR-H Model
The SEIR-H model is a compartmental model. For the COVID-19 spread dynamics, the SEIR-H model is nearly identical to the classic SEIR model, with the difference that in the proposed model, the group R stands for the reported individuals, who are diagnosed as confirmed positive cases. This reported group is split mainly into three groups: (i) H: hospitalised group, composed of individuals who have severe symptoms, and thus are hospitalised, (ii) U: a group composed of individuals who are admitted to ICU due to the worsening of the symptoms, and (iii) F: death. Thus, the SEIR-H model is separated into two parts representing two dynamics: (i) the dynamics of spread/transmission of the disease, and (ii) the dynamics of the use of hospital beds. In this model, the dynamics of hospital resources do not affect the spread dynamics of the disease, i.e., the compartments S, E, I, R are not influenced by H, U and F. This assumption is reasonable in the sense that the disease transmissions from the hospitalised COVID-19 patients are negligible compared to the transmissions occurring in the community outside the hospitals. Figure 1 schematically illustrates the proposed SEIR-H model. Both parts of the model are highlighted: (i) the spread dynamics with yellow background, (ii) the dynamics of reported cases in green background. The parameters in blue are associated with the proportions of individuals moving from one compartment to the next, and the parameters in orange are associated with the reciprocal of the average length of stay in each compartment. In the next subsections, both parts are discussed in detail by presenting the systems of differential equations along with the definitions of the parameters. The compartment O in Figure 1, corresponds to those who are recovered from COVID-19. In this work, we do not track the count of recovered individuals, because the data of those recovered are unreliable. Thus, the equation for this compartment is omitted. Figure 1. Schematic diagrams of the SEIR-H model: the region in yellow indicates the transmission dynamics given by a SEIR with external reported cases, and the region in green represents the dynamics related to hospital resources. A more detailed version of the green region is used to show that the R compartment contains the dynamics related to the health system, represented mainly by H, U and F.

Dynamics of The Spread
The spread dynamics of the disease, highlighted in yellow in Figure 1, is similar to the SEIR model [15], and is given by: In this model, the Equations (1)-(4), respectively, show the dynamics of the compartments of: the susceptible S, the exposed E, the infectious I and the reported R. The term R im represents a rate of reported cases from the travellers, which was relevant at the beginning of the dynamics of COVID-19 in Paraguay. The susceptible group S represents those who have not yet contracted the disease and could contract it in the future.
The exposed group E represents those who were exposed to the infection and have been infected, but are in a latent period not yet capable of infecting others. The latent period refers to the time between exposure to the virus and the moment when the individual becomes contagious and is one of the epidemiological characteristics of COVID-19 [8,[18][19][20][21] While in the classic SEIR model, compartment E could include those in the incubation period [22], in this work, however, we only include in E those in the latent period. It is crucial to establish this distinction between the latent period and the incubation period, since COVID-19 is characterised by a pre-symptomatic infection [9,10,23]. Furthermore, the inclusion of the E-compartment can resemble the Erlang distribution in the infectious period [11].
The infectious compartment I represents individuals who become contagious. The key feature which distinguishes compartments I from E, is the onset of infectivity. Individuals in I may transmit the disease to susceptibles in S, even without presenting symptoms.
The compartment R is composed of those who are reported with the confirmation of COVID-19 positive diagnosis. In the classic SEIR model, the R compartment is associated with those who are recovered or removed from the dynamics after a certain time. In this work, the R compartment also has the meaning of removed, given the assumption that those who are reported are also isolated and removed from the transmission dynamics. However, since R will be fitted to the officially positive confirmations, we rename it as the reported compartment. The population N is maintained constant throughout the time, equal to the initial population. This condition differs from the actual population that varies due to the reported cases from the travellers R im , and also because of the individuals who are travelling abroad (which are not considered in the proposed model).
Three parameters appear in this model: (i) β(t) is the rate at which susceptible individuals acquire the disease, also known as transmissibility, which is time-dependent in this work, (ii) α is the rate at which an exposed individual becomes infectious, also defined as the reciprocal of the average latent period, and (iii) γ is the rate at which an infectious individual is reported, expressed as the reciprocal of the average time an individual remains contagious to other individuals. Table 1 shows the parameters of the SEIR-H model, with its descriptions and values, and whether it is considered to be constant or time-dependent.

Reproduction Number
The dynamics of disease transmission can be expressed using a single parameter called the reproduction number R 0 . This number can be defined as the expected number of secondary cases that arise after introducing an infectious individual into a susceptible population. This reproduction number is useful to define a threshold value to be maintained for a certain time, for the epidemic to persist (R 0 > 1), or be extinguished (R 0 < 1) [15,27,28].
The reproduction number is affected by various biological, socio-behavioural and environmental factors [27]. The numerical value depends on the estimation method, and the specific value obtained from the mathematical models based on differential equations should be understood as a threshold rather than as the number of secondary cases [28]. The numerical value can be considered a function of three parameters: the duration of the infectious period, the average probability of transmission of the virus in each contact and the average rate of contacts [27].
The reproduction number of the classic SEIR model can be expressed in terms of the parameters β and γ, as follows: As the spread of disease is evolving, the number of susceptible individuals declines monotonically. Thus, the number in which effectively works as an epidemic threshold is: also known as the effective reproduction number. This expression is obtained by adding the right-hand sides of the Equations (2) and (3) and then making it equal to zero. From these equations, it can be observed that the value of the reproduction number R 0 or R e (t) can be altered by interventions on parameters β and γ. Containment and social distancing measures can affect the parameter β, by reducing the average rate of contacts. However, hygiene recommendations and proper use of face masks can also affect β by reducing the average probability of virus transmission in each contact. While the parameter γ depends on the disease itself, it can be altered by effective tracing, testing and isolation of infected individuals.

Dynamics of Required Hospital Resources
The dynamics of reported cases, highlighted in green in Figure 1, is mainly focused on the number of reported individuals in need of general hospitalisations and Intensive Care Units (ICUs). Hence, it is important to analyse the portion of infected individuals that develop severe cases of the disease as well as the fatality rate.
A compartmental model is considered: (i) the compartment H includes general hospitalisation of COVID-19 patients, (ii) the compartment U includes all those reported individuals requiring ICU, and (iii) the compartment of the deceased F. In this work, we use the letter F for the deceased because it is the first letter of the translation in Spanish of the deceased (Fallecido).
The equations that describe the dynamics in hospitals can be expressed as follows: The dynamics of the portion of infected individuals who develop severe cases and need to be hospitalised, is represented in Equation (7), by using a parameter λ I→H . A portion of those who are hospitalised, given by the parameter λ H→U , require intensive care. The average length of stay at a normal bed before requiring an intensive care unit, is given by 1/δ H→U . However, hospitalised individuals are dying with portion λ H→F with average staying in a hospital of 1/δ H→F . The rest of the portion λ H→O = 1 − λ H→U − λ H→F is recovered with average staying in hospital of 1/δ H→O .
The Equation (8) represents the dynamics of the occupation of ICU. The admission in ICU is restricted to individuals with severe cases who are initially hospitalised in normal beds. The portion of dying individuals from the ICU is λ U→F , with an average stay in ICU of 1/ϕ U→F . The recovered portion is λ U→O = 1 − λ U→F , with average stay in ICU of 1/ϕ U→O .
The last compartment corresponds to accumulated death, given by the Equation (9). Here, it is included those who die from those hospitalised and also those who were in ICU. In addition, a certain portion given by λ I→F dies without admission in the hospitals.

Containment Measures in Paraguay
The first confirmed COVID-19 positive case in Paraguay was reported on 7 March 2020 (First case of the new coronavirus in Paraguay: https://www.mspbs.gov.py/portal/20535/ primer-caso-del-nuevo-coronavirus-en-el-paraguay.html, accessed on 20 May 2021). Before the COVID-19 pandemic has reached Paraguay, the country was emerging from one of the worst dengue fever epidemics (The worst Dengue epidemic has been contained: https://www. mspbs.gov.py/portal/20525/la-peor-epidemia-de-dengue-esta-contenida.html, accessed on 20 May 2021). Thus, the government had adopted very strong containment policies in order to gain some time to strengthen the health system (Gs. 530 billion in supplies and 2700 health contracts: https://www.mspbs.gov.py/portal/20576/g-530-mil-millones-en-insumosy-2700-contratos-para-salud.html, accessed on 20 May 2021; Personal protective equipment delivered to health regions: https://www.mspbs.gov.py/portal/20621/entregan-equipos-deproteccion-personal-a-regiones-sanitarias.html, accessed on 20 May 2021). This lockdown lasted at least two months, after which the government adopted a system of gradual lifting of the containment measures with enumerated phases: from the most rigorous confinement at Phase 0 to the more relaxed at Phase 4. After that, the policy named "New normality" is applied, which is the most relaxed measure with specific restrictions in certain areas depending on its situation. Table 2 summarise the different phases of quarantine (while the term quarantine is used for mobility restrictions on exposed individuals, in Paraguay, this term is widely adopted to refer to the various degrees of lockdown measures) in Paraguay. Local and temporary restrictions that differ from the main phases can be consulted in Appendix A.

Estimation of Model Parameters
This section presents the methodology to estimate the parameters of the proposed model. For this purpose, we use data provided by the Ministry of Public Health and Social Welfare, and the Health Monitoring Centre. The available data are daily records of: diagnosed with positive confirmed cases or daily reported cases D R (t j ), excluding the count of confirmed cases from travellers, 2.
diagnosed with positive confirmed cases in the isolated groups who were travellers from abroad R im (t j ), 3. occupation in general bed in the hospitals by the patients of COVID-19 D H (t j ), 4. occupation in ICU by patients of COVID-19 D U (t j ), 5. confirmed daily deaths by COVID-19 D F (t j ). These data in a specific day are assumed to be independent of each other. These data are from t 1 = 7 March 2020 until t j end = 22 May 2021. From the data provided by the Health Monitoring Centre, it was determined: (i) the average time of stay in general beds in the hospitals by the COVID-19 patients, and (ii) the average time of ICU stay by COVID-19 patients.
The estimation of the model parameters is mainly performed in two steps: (i) the estimation of the transmissibility β of the disease for the spread dynamics using the Equations (1)-(4), and (ii) the different portions λ * , with * ∈ {I → H, I → F, H → U, H → F, U → F}, for bed occupations in hospitals using the Equations (7)-(9) in addition with the Equations (1)-(4). The estimation of transmissibility β is performed from the date of the first confirmed positive case in Paraguay, at t 1 . Due to the initial strong containment, the amount of hospital bed occupations by COVID-19 patients remained very low for several months.
Because of this conditions, and due to the high intermittency of the available data, there were difficulties in estimating the portions λ * , thus these values are estimated only as of t j h = 5 October 2020, when the "new normality" policy measures are applied. Consequently, the states of H, U and F, which are associated with these proportions are simulated from t j h .
The time-varying parameters are determined using the moving time-window. In this work, we fixed the size of the time-window to W = 14 days. In the classic SEIR model, the time-window size would correspond to the entire epidemic cycle, keeping the transmissibility β invariant over the cycle. This hypothesis is a reasonable assumption in seasonal flu epidemics, because most of the population do not change their behaviour during the season. However, this condition is not the case with the COVID-19 pandemic, because behavioural changes of the population are continuous and vary mainly in response to official containment measures. Thus, in order to capture the relevant dynamics of transmission, a mechanism is needed to approach a time-varying transmissibility.
The strategy of using a moving time-window in the parameter estimation resembles the moving average for time-series data, i.e., a smaller time-window size will better fit daily noisy data with higher oscillations, and on the other hand, a larger time-window size will offer smoother results and reveal trends. Thus, parameters obtained using timewindows with size W = 14 days can be understood as an average during this time span. This size was chosen mainly because of two reasons: (i) according to the officials at the Health Monitoring Centre in Paraguay (DGVS), the delay from the day when a contagion occurs until the detection is normally at least 10 days, and (ii) the strategy of containment measures in Paraguay is usually dictated in a two-week time span.
Algorithm 1 shows schematically the procedure for parameter estimation. The required information is the set of parameters which are maintained constant (see Table 1), the data D R (t j ), D H (t j ), D U (t j ), D F (t j ), R im (t j ), from j = 1 to j = j end , population N and size of time-window W. The results are the transmissibility β(t j ) from j = 0 to j = j end − W, and the proportions λ * (t j ) from j = j h to j = j end − W.

. Estimation of Initial Condition
The population of Paraguay is estimated at N = 7,252,671 according to the National Institute of Statistics (INE: Instituto Nacional de Estadística https://www.ine.gov.py/default. php?publicacion=2, accessed on 20 May 2021). The initial time is set on t 0 = 6 March 2020, one day before the day when the first case is reported. On this day, the reported group is set as R(t 0 ) = 0, and the initial conditions for: the exposed e 0 = E(t 0 ), the infectious i 0 = I(t 0 ) and the transmissibility β(t 0 ) are estimated using D R 0 = {D R (t 0 + k) | k = 1, . . . , W}, available data of daily reported or confirmed cases in the first W = 14 days (Algorithm 1, line 1). The initial susceptible population is Some parameters that do not depend on the particularities of each country or region but on the disease itself are set according to the values presented in the literature. The average of the latent period is 1/α = 3 days [18][19][20][21]24,25] and the average of the infectious period is 1/γ = 7 days [18,20,25,26] (see Table 1). Although the infectious period can be reduced by early tracing and isolation (quarantine) strategies for infectious persons, after the movement restriction period and throughout much of the pandemic in Paraguay, early tracing and isolation of infectious persons was difficult.
The estimation is performed using the Bayesian method, which can be expressed as: where P(e 0 , i 0 , β 0 , σ R 0 |D R 0 ) is the posterior probability distribution of the parameters {e 0 , i 0 , β 0 , σ R 0 } given the observation of the data D R 0 , P(D R 0 |e 0 , i 0 , β 0 , σ R 0 ) is the probability of observing the data D R 0 given specific values of the parameters, which is also known as likelihood function, and P(e 0 , i 0 , β 0 , σ R 0 ) is the prior probability distribution of the parameters. We assumed that the observed data D R 0 are normal-distributed around the daily reported cases dR e 0 ,i 0 ,β 0 from the simulations results using the parameter set {e 0 , i 0 , β 0 } with standard deviation σ R 0 . Thus, the likelihood function can be expressed as: where given a set of values of the parameters {e 0 , i 0 , β 0 }, the simulated daily reported cases dR e 0 ,i 0 ,β 0 are obtained by: where R e 0 ,i 0 ,β 0 is the accumulated count of reported cases from the simulation results. This standard deviation σ R 0 is restricted to the interval [0,100]. The priors are assumed as restricted uniform distributions, with a search range of [0, 100] for both initial exposed and infectious populations, e 0 and i 0 , and [0, 0.8] for the transmissibility β (see Table 3).
The library stan [31] in the platform R [32] is used for estimating the posterior distributions of the parameters. In the estimation, stan uses the Hamiltonian Monte Carlo (HCM) algorithm of Markov chain Monte Carlo method (MCMC) [33]. The log probability l p (the logarithm of the right-hand side of the Equation (10)) is evaluated internally in stan, and the parameter set from the posterior distribution samplings that gives the maximum log probability is chosen for the optimum values: where N iter is the number of samplings of the MCMC method. This number of samplings with the amount of Markov chains are maintained large enough to ensure the convergence of the parameters obtained: 30 Markov-chains with 5000 iterations for each chain including, 1000 warm-up iterations, and the average proposal acceptance probability adapt_delta [31] is set to 0.94. These optimum values of the parameters are used to advance one day, and to obtain the state of the variables: (S(t 0 + 1), E(t 0 + 1), I(t 0 + 1), R(t 0 + 1)). To solve the system of differential equations, stan uses the fourth and fifth-order Runge-Kutta method.
As the parameters of transmissibility β are time-dependent, it can be expressed as: where β opt j is a specific value among the MCMC samples that gives the optimal l p at time t j using the time-window [t j , t j + W], and χ is the indicator function defined as: Similarly, the proportions λ * can be defined as: where λ opt * j are a MCMC sample which gives the optimal l p at time t j .

Estimation of Transmissibility
The estimation of the transmissibility of the j-th day, β j at t j = t 0 + j, from the initial day t 0 is described (Algorithm 1, line 3). The estimated transmissibility of a day before, β opt j−1 , is used in order to obtain the state of the variables at time t j : (S(t j ), E(t j ), I(t j ), R(t j )). This state is adopted as the initial condition for estimating the transmissibility β j by performing the Bayesian method as follows: The interpretation of the different terms of this equation is similar to the Equation (10), but in a time-window of W starting on t j : The observed data is D R j = {D R (t j + k)|k = 1, . . . , W}, and these data are assumed to be normal-distributed around the daily reported cases dR β j from the simulation results with standard deviation σ R j (restricted to the interval [0, 10,000]). The likelihood function can be expressed as: where, the daily reported cases dR β j are obtained from: where R β j is the accumulated count of reported cases from the simulation using the parameter β j . The priors for β j are also restricted to uniform distributions with a search range of [0, 0.8] (see Table 3). The value for β j is estimated as the one that gives the maximum log probability (the logarithm of the right-hand side of the Equation (16)) performed using MCMC: β opt j = {β j m |max m∈[1,N iter ] l p }. 30 Markov-chains with 3000 iterations for each chain, including 1000 warm-up iterations are used, and the average proposal acceptance probability adapt_delta [31] is set initially to 0.94. This average proposal acceptance probability is dynamically augmented by 0.01 when divergent transitions after warm-up are observed. In most cases, this value is remained in 0.94, and there are a few cases when it is augmented up to 0.97.

Estimation of the Parameters for the Dynamics in the Hospitals
As mentioned above, the different proportions λ * , with * ∈ {I → H, I → F, H → U, H → F, U → F} for the dynamics in the hospitals, are estimated as of t j h = 5 October 2020 (Algorithm 1, line 5). In order to estimate the proportions λ * (t j h ), the initial condition for the SEIR part is obtained using the optimal estimated transmissibility at time t j h − 1 (β opt j h −1 ), and the initial state of the compartments H, U and F are set using the 7-days average (from t j h − 3 to t j h + 3) of the data of hospitalisation (normal bed), ICU and daily death, respectively.
For the estimation of the proportions λ * at time t j = t 0 + j (j > j h ), the state (S(t j ), E(t j ), I(t j ), R(t j ), H(t j ), U(t j ), F(t j )), obtained using the optimal values of the estimated parameters β opt j−1 and λ opt * j−1 at day before, is adopted as initial condition. The parameters associated with the average lengths of stay in the hospitals are set to be constant. These times are estimated from individualised case record data provided by the DVGS (see Table 1 for the specific values).
The estimation of the proportions λ * is performed using also the Bayesian method, maintaining constant the parameter β opt j obtained beforehand, and it is expressed as: The observed data are:  (19) is a product of three likelihoods: P(D H j |λ * j , σ H j ), P(D U j |λ * j , σ U j ), and P(D F j |λ * j , σ F j ). The likelihood functions can be expressed as: where, H λ * j and U λ * j are the hospitalised (in normal bed) and in ICU from the simulation using {β opt j , λ * j }. The daily deaths from the simulation dF λ * j are obtained from: where F λ * j is the accumulated count of death, obtained from the simulation using {β opt j , λ * j }. Based on previous data, before 5 October 2020, priors are assumed as restricted uniform distributions, with the search ranges for: and particularly for λ U→F it is considered a normal distribution with mean µ = 0.4, standard deviation σ = 0.1, and search range [0.2, 1] (see Table 3). The values for λ opt * j are selected those which give the maximum log probability (the logarithm of the right-hand side of the Equation (19)) performed by stan. 60 Markov-chains with 3000 iterations for each chain including 1000 warm-up iterations are used, and the average proposal acceptance probability adapt_delta [31] is set initially to 0.94. Similar to the transmissibility estimation, this average proposal acceptance probability is dynamically augmented by 0.01 when divergent transitions after warm-up are observed. However, this value remains at 0.94 at almost all times, and there are a few cases when it is augmented up to 0.95.

Results
The results obtained from the parameter estimations are presented in Section 3.1. Section 3.2 shows the model assessment together with the projection error horizon considering nearly constant parameters. Initially, the optimal transmissibility is about 0.3. This transmissibility gives the reproduction number about 2, which is actually an estimated value presented in the first studies conducted on COVID-19 pandemic [34].

Estimated Parameters
A few days after the appearance of the first reported case, Paraguay has implemented a strict pandemic containment regime with nationwide mobility restrictions for at least 2 months (see Table 2  To further analyse these variations in transmissibility between May and June, Figure 3 overlays the β parameter with the simulated values for the exposed E and infectious I compartments, as well as the local confirmed data, i.e., excluding imported cases R im . The distinction between local and imported cases is important, as in the model considered, imported cases R im are mainly associated with those who are already isolated from the outset in shelters and do not participate in the transmission dynamics. Figure 3 shows how the amounts of exposed and infectious also decrease to values very close to zero. This can be interpreted as a resetting of the model dynamics around Transmissibility E from simulation using optimal opt I from simulation using optimal opt Daily reported from data Optimal opt Value of for R 0 = 1 From Figure 3, it can also be observed that the exposed E and infectious I decrease if β is lower than the threshold and increase as it becomes larger. Furthermore, the delay of at least ten days normally observed between the occurrence of contagion and the detection can be clearly noticed here when beta adopts higher values and its effects are seen several days later with the increase of daily reported.
From August onwards, there is a sustained increase in the number of cases reported daily, but the various containment measures implemented have allowed that the transmissibility remains around the value that corresponds to a reproduction number equal to one, at least until December. The particularity of December as a summer festive month and the increase in the number and frequency of social gatherings, reveal an increase in the transmissibility in the last week of December, as shown in Figure 2.
In general, we observe that the simulations of the proposed model in its transmissibility dynamics reflect the different circumstances of COVID-19 in Paraguay throughout the year 2020. Figure 4 shows the evolution of the estimated λ opt * (t) proportions defined in Section 2.3.3 as of 5 October 2020, from the beginning of the new normality phase as described in Table 2. According to data provided by the Ministry of Public Health and Social Welfare, by 31 December 2020, 163 out of 461 available intensive care units (ICU) for respiratory diseases were occupied by patients diagnosed with COVID-19. As can be observed from the data points in Figure 5, between January and March 2021, while the number of hospitalisations was nearly constant the number of patients in ICU at least doubled. This trend is reflected in the parameter λ H→U growing from 0.15 to 0.25 in this period.

Estimated Parameters for Hospital Resources
Although the total ICU capacity to receive COVID-19 patients was extended to a total of around 600 units by May 2021, the rapidly growing number of hospitalisations between March and May 2021 overwhelmed the health system, leading to increased waiting lines to access an ICU and a sudden increase of daily deaths. This trend is coherent with the decrease of λ H→U due to the unavailability of ICU and increase of λ H→F . The increase of λ I→H from 0.08 to 0.13 could be either attributed to a more aggressive COVID-19 variant or to the limited testing regime, which leads to underreported cases. In general, λ U→F tends to stay around 0.4.

Model Assessment
In this section, the ability of the model to fit the available data is presented, as well as an analysis of the error in the projections under the assumption of nearly constant model parameters. Figure 5 shows the state trajectories of the SEIR-H model along with reported data points. The trajectories are obtained by solving the system of differential Equations (1)-(4) and (7)-(9), with β opt and λ opt * (t) ( * = I → H, I → F, H → U, H → F, U → F) parameters obtained from the parameter estimation explained in Section 2.3. The trajectories shown (from the top to bottom) are the daily reported dR β opt ,λ opt * , normal bed occupation in hospitals H β opt ,λ opt * , occupation in ICU U β opt ,λ opt * , and daily death dF β opt ,λ opt * . As happened with the proportions λ * , the simulated trajectories of hospitalised, ICU and daily death are plotted starting on 5 October 2020, which is when the "new normality" measures were enforced. The simulation trajectories stop on 1 April 2021 and the posterior distributions estimated for that day are used to advance trajectories for an additional month with the percentile 95 credibility intervals for the trajectories.  From Figure 5, we observe that the trajectories from the simulation can follow very well the data. Actually, there is a lot of variability in the data, which is due to the limited testing regime and delays in reporting of deaths. On the weekends, the amount of testing and also the number of reports of the deaths are reduced. Thus, higher variability is observed in daily reported and daily deaths. The relative mean deviation of the data can be defined using a moving standard deviation and moving average of length 7 days. The moving standard deviation of the data D at time t j is expressed as:

Comparisons
where movave D (t j ) is the moving average of the data at time t j , which can be defined as: The relative mean deviation of the data is expressed as: For the hospitalised and ICU, the relative mean deviation are rel_dev D H ≈ 0.039 and rel_dev D U ≈ 0.036, respectively, meanwhile, for the daily reported and daily deaths are rel_dev D R ≈ 0.413 and rel_dev D F ≈ 0.204, respectively. Data from j * = 1 to t j end − W is used in the computation of rel_dev D R , and from j * = j h to t j end − W for the rest.
In order to assess the model, the deviation of the relative error between the data and the simulation results is computed. For this task, the relative error between the data and the simulation results at time t j is defined as: where X is just a representation of any simulation results, whether it is daily reported, hospitalised, ICU or daily death. Thus, the deviation of the relative error is expressed as: For the hospitalised and ICU the deviation of relative error are dev_rel_error H ≈ 0.047 and dev_rel_error U ≈ 0.050, respectively, meanwhile, for the daily reported and daily deaths are dev_rel_error R ≈ 0.668 and dev_rel_error F ≈ 0.234, respectively (see Table 4). These values are comparable to the relative deviation of the data, and this means that the methodology adopted in this work, by using the moving time-window to estimate the parameters, allows us to obtain the time-varying parameters that give the simulation trajectories almost similar to the moving average of the data.

Projection Error under the Assumption of Nearly Constant Parameter
This section presents a statistical analysis of the behaviour of the error as a function of the projection horizon. An issue criticised to epidemiological models during the pandemic was the error in the projections. These errors can have their origin in different factors, e.g., the difficulty in detecting positive cases due to a limited testing rate, delay and difficulties in confirming deaths related to COVID-19, the assumptions that social behaviour will be homogeneous and stationary, among others. These sources of error are dynamic, therefore they force a continuous update of the model parameters and limit the projection horizon for which the parameters estimated at a given moment will be appropriate. Furthermore, from the perspective of a decision-maker, it is important to be informed about the limitations of a modelling framework to anticipate scenarios.
We estimate the error in the projection under the assumption of nearly constant values for the parameters. For the projection from the time t j+W , a Sigmoid function is used, where the range is from the parameter value at t j to an average in 2W = 28 days of the parameter, with a variation width of about W/2 = 7 days starting from t j+W . For the transmissibility parameter, the mentioned Sigmoid function can be expressed as: whereβ is obtained by averaging the values of transmissibility for the last 2W = 28 days: from t j−W to t j+W − 1. Using the available data until t j+W , β can be determined until t j , i.e., β(t j ) gives the better fit of data D R j . Thus, the average transmissibilityβ for the last 2W days is composed of W days with β varying daily from t j−W to t j−1 , and β is constant with value β(t j ) from t j . This average can be expressed as: .
The Sigmoid function is chosen because, in one hand, the parameter value β(t j ) fit better to data for the last W days, but in the other hand,β gives a representative transmissibility for longer-term if there is no change in the containment policies or behaviour on the population scale. The Sigmoid function is just a simple way to join smoothly these two values of the parameters.
Similar functions are used for the parameters of proportions. Thus, the projections are obtained using these parameters, and solving the system of differential equations. The results obtained are named as dR j , H j , U j , and dF j .
The relative error (Equation (27)) is iteratively computed using the results obtained from the simulations dR j , H j , U j , and dF j , with j such that, t j h = 5 October 2020 until 1 March 2021. Thus, vectors of relative error are obtained, which are denoted as rel_error X j (t j+k ) for k = 1, 2, . . . , W + W p : composed of W = 14 days of training and W p = 30 days of test or projection. These vectors are visualised using box plots for each k = 1, 2, . . . , W + W p . In the analysis, the days are remapped in order that the projection days could start from 1. Thus, the day n is related to k by n = k − 14, where k = 1, 2, . . . , W + W p . Figure 6 illustrates the statistical behaviour of the relative error corresponding to the observation day n. As can be seen in the box plots, the most probable error in all projections tends to show an average of 0. However, this average may increase in function of the projection days for deaths, evidencing an underestimation tendency in terms of precision (observed higher than expected). The red lines in the figures indicate a 20% threshold of error. It was evident that up to W p =30 days of projection (without considering the W = 14 of training), the error is more frequently within this threshold, with a more frequent tendency of underestimation. Also note that as the projection horizon increases, the margin of error observed increases by the interquartile range (IQR), whiskers Q {1,3} ± 1.5(IQR), and outliers, that is, the accuracy is compromised. Consequently, it is evident that both the precision and accuracy of the model are compromised with the increase in projection days, thus increasing the uncertainty in the estimate over time. This result reinforces the need for periodic updates in the estimation of the model parameters for the projections.

Discussion
The estimated parameters helped us to interpret and analyse the spread dynamics of COVID-19 in different circumstances, which was particular in Paraguay. The relative error between the simulation results using the time-varying parameters and the available data is comparable to the relative deviation of the data using the moving time-window. This means that the simulation trajectories represent very well the variability in the population scale. The obtained transmissibility β changes depending on the epidemic containment policies applied. And, the different proportions λ are also changing due to the containment policies, but primarily on the availability of the beds and ICU in the hospitals. As the number of severe cases of COVID-19 patients becomes larger than the available beds and ICU in the hospitals, the proportions of deaths are also increasing.
The statistical analysis of the behaviour of the error showed reasonable results [35], with a tendency to underestimate the observed data. This problem is somewhat inevitable considering the assumption of nearly constant parameters. As a decision-making tool, several scenarios can be drawn, for example: (i) plateau: by considering the transmissibility that stabilises daily reported infections as approximately a constant, (ii) the average value of transmissibility over the last month, (iii) 10% decrease in transmissibility when stronger quarantine measures are enforced, or (iv) 10% increase in transmissibility when quarantine measures are relaxed. All these scenarios can be based on the recently estimated parameters, which give additional advantages in the sense that, in one hand, the transmissibility of the disease can change over time due to the different strains of the viruses, and, on the other hand, after several containment measures, a portion of the population respecting the measures can vary compared to previous phases.
The early lockdown in Paraguay allowed for a slow increase in the number of COVID-19 patients, and thus, severe cases could be treated adequately by the health system until the beginning of 2021. Most of the time in 2020, the daily cases and deaths reported followed a plateau, but starting February 2021, cases increased substantially and the death proportions are becoming larger due to the collapse of the health system.

Relevance of the Work
Despite its simplicity, the SEIR-H model proposed in this work offers results relevant for the evolution of cases for COVID-19 in Paraguay. A visualisation tool for these results based on data from Paraguay is available at http://epymodel.uaa.edu.py (accessed on 12 October 2021). The Ministry of Public Health and Social Welfare and its Health Monitoring Centre in Paraguay use and reference this tool for analysis and decision-making (http: //dgvs.mspbs.gov.py/page/#1387, accessed on 20 May 2021). Currently, results from this simple model are used to evaluate future scenarios and plan for resources allocation in the public health budget.
Results and projections obtained using the proposed SEIR-H model, exhibits similar statistical performance compared with other public platforms such as the COVID-19 Projections from IHME (Institute for Health Metrics and Evaluation) (https://covid19. healthdata.org/, accessed on 20 May 2021) [35] which serves as a reference for multiple governments and agencies. A sample comparison of the SEIR-H and IHME model results is presented in Appendix B.

Conclusions and Future Work
In this work, we proposed the SEIR-H epidemiological model with time-varying parameters in the form of step functions. This model is based on the classical SEIR model maintaining its simplicity in the spread dynamics of the disease. However, the SEIR-H model contains the dynamics of the occupations of normal beds and ICU in hospitals in order to guide decisions based on available hospital resources. This model was designed to have as few parameters as possible to be estimated, while still representing the dynamics observed in Paraguay.
We summarise the contributions of this work as follows: • Integration of health system dynamics into the SEIR model, without losing the simplicity of the SEIR model. • The methodology of parameter estimation for the time-varying parameters using a moving time-window, which allows us to determine the variability of the parameters according to the spread dynamics and social behaviour. • The inclusion of the reported cases that were travellers from abroad in the inflow of reported cases, which is important at the beginning of an epidemic. • Model assessment and statistical analysis of the error as a function of the projection horizon.
• Discussion of spread dynamics in Paraguay using the estimated parameters and trajectories obtained.
The proposed model has two main limitations. The first limitation is the assumption that all cases end up being reported ignoring the under-reporting of cases, which can be quite significant [36]. The second limitation is to ignore the limits of hospital capacity, both in terms of available care units for medium and severe cases. This issue is especially important for projections of deaths when the health system is overburdened [37]. Future work can address these aspects by adding a model for under-reported cases and by limiting the capacity of the health system.
Although, the SEIR-H model is fitted to the case of Paraguay, this model can be useful for any country or region to have a macroscopic overview of the dynamics of an epidemic like COVID-19, especially in contexts of limited availability of data detailing the heterogeneity of the population.

Conflicts of Interest:
The authors declare no conflict of interest. Figure A1 gives a geographical and demographic contextualisation of the detailed restrictions. Figure A1. The population of Paraguay

Appendix B. Performance Comparison
In this appendix, we present a comparison between the proposed SEIR-H model results with the COVID-19 daily death projections from the IHME platform. Daily death reports are considered the most complete, in the sense that it is less subject to underreporting. Figure A2 illustrates a comparative analysis of the statistical performance for daily deaths forecasts (at a time interval of 41 days) between these two models.
The expected and observed forecasts of the proposed model adjust to a line with a slope 0.79. However, the expected and observed daily deaths from the forecasts of the IHME platform (COVID-19 current projection) show a slope 0.78. Therefore, both projections are similar, which can be confirmed by the last graph where the fit of the expected values from the IHME and the expected values from the SEIR-H model shows a line with unit slop. Thus, this analysis shows a remarkable statistical correspondence between both models.