Modeling the early phase of the Belgian COVID-19 epidemic using a stochastic compartmental model and studying its implied future trajectories

2.1 Mathematical compartmental transmission model

We use an adapted version of an SEIR mathematical compartmental model to describe COVID-19 disease dynamics. In this model, individuals are susceptible to infection when in compartment S, and after an effective contact (between a susceptible and infectious individual) the susceptible individual moves to an exposed state E at age- and time-specific rate λ(t), referred to as the force of infection (with boldface notation representing a vector including age-specific rates). After a latent period, the individual becomes infectious and moves to a pre-symptomatic state I_presym at rate γ. Afterwards, individuals either develop symptoms (state I_mild) with probability 1 − p or remain completely free of symptoms (compartment I_asym, probability p). Asymptomatic cases recover at rate δ₁. Symptomatic infections are either very mild and such cases recover at rate (for an age-dependent fraction of these individuals, represented by the vector ϕ₀ = δ₂/(δ₂ + ψ), where and or they move to a state I_sev prior to requiring hospitalization (i.e., severe infection is defined as requiring hospitalization). When severely ill, implying hospitalization, individuals move to state I_hosp with probability ϕ₁, or become critically ill (I_ICU) with probability 1 −ϕ₁. Hospitalized and critically ill patients admitted to the Intensive Care Unit (ICU) recover at rate δ₃ and δ₄ with probabilities {1 − µ_hosp(a)} and {1 − µ_icu(a)}, respectively, where µ_hosp(a) = τ ₁/(τ ₁ + δ₃) and µ_icu(a) = τ ₂/(τ ₂ + δ₄) represent the age-specific case-fatality rates (i.e., probabilities of dying when severely ill and hospitalized on a general ward or admitted to ICU). Hospitalized and ICU patients die at rate τ ₁ or τ ₂ with probabilities µ_hosp(a) and µ_icu(a), respectively. A schematic overview of the compartmental model is given in Figure 1. Individuals in the red compartments are able to transmit the disease.

• Download figure
• Open in new tab

Figure 1:

Schematic overview of the flows of individuals in the compartmental model: Following SARS-CoV-2/COVID-19 infection susceptible individuals (S) move to an exposed state (E) and after a latent period individuals further progress to a pre-symptomatic state (I_presym) in which they can infect others. Consequently, individuals stay either completely symptom-free (I_asym) or develop mild symptoms (I_mild). Asymptomatic individuals will recover over time. Upon having mild symptoms, persons either recover (R) or require hospitalization (going from I_sev to I_hosp or I_icu) prior to recovery (R) or death (D).

The following set of ordinary differential equations describes the flows in the (deterministic version of the) proposed age-structured compartmental model: where, for example, S = (S₁(t), S₂(t), …, S_K(t))^T represents the vector of number of susceptible individuals in age group k = 1, …, K in the population at time t. A full account on the definition of the different compartments and the notation used for the number of individuals therein can be found in Table A1 in Appendix A. An overview of the different parameter definitions can be found in Table B3.

The proposed age-structured compartmental transmission model consists of 10 age classes, i.e., [0-10), [10-20), [20-30), [30-40), [40-50), [50-60), [60-70), [70-80), [80-90), [90, ∞) with the number of individuals in each age class obtained from Eurostat.

2.2 Social contact data and transmission rates

As mentioned previously, the infectious phase of COVID-19 disease is divided into two different states: a pre-symptomatic state occurring before the end of the incubation period, followed by a state in which individuals may either remain asymptomatic or develop (mild to severe) symptoms (see Figure 1). Transmission of the disease is governed by an age- and time-dependent force of infection. The age-specific force of infection in age group k = 1, …, K, denoted by λ(k, t), is defined as the instantaneous rate at which a susceptible person in age group k acquires infection at time t. Furthermore, the time-invariant transmission rate β(k, k ′) represents the average per capita rate at which an infectious individual in age group k makes an effective contact with a susceptible individual in age group k, per unit of time. Consequently, the force of infection is defined as where I(k ′, t) denotes the total number of infectious individuals in age group k ′at time t and β(k, k′) can be rendered as when relying on the so-called social contact hypothesis [6]. This hypothesis entails that c(k, k ′) are the per capita rates at which an individual in age group k makes contact with an individual in age group k ′, per unit of time, and q is a proportionality factor capturing contextual and host- and disease-specific characteristics such as susceptibility and infectiousness. The (K × K)-matrix C containing the elements c(k, k ′) is referred to as the social contact matrix describing mixing behaviour within and between different age groups in the population. Social contact rates c(k, k) are estimated based on social contact data from Flanders (Belgium) collected in 2010 [7, 8, 4, 9]. We hereby assume that contact rates for Flanders can be used for all regions in Belgium.

In this manuscript, we rely on social contact matrices C_sym and C_asym estimated for symptomatic and asymptomatic individuals, implying defining transmission rates for both symptomatic and asymptomatic cases, respectively (see Appendix C). Here, we assume that the relative infectiousness of symptomatic versus asymptomatic cases is equal to . The age-dependent force of infection is defined as: where λ(t) = (λ(1, t), λ(2, t), …, λ(K, t)) written as a matrix multiplication with boldface notation for vectors and matrices. Note that hospitalized individuals are assumed not to contribute to the transmission process because of isolation.

2.3 Discrete time stochastic epidemic model

The spread of the virus is hampered by reductions in the number of contacts and changes in the way contacts are made, either voluntarily or as a consequence of government intervention. These time- (and age-) dependent behavioural changes introduce substantial uncertainty in the further course of the outbreak and require stochastic model components to evaluate the effectiveness of the intervention strategies and to make future predictions in terms of, for example, new hospitalizations. Moreover, stochastic epidemic models allow to determine the probability of extinction based on multiple realizations of the model. Therefore, we amended the deterministic model hitherto described into a discrete time stochastic epidemic model to describe the transmission process under the mitigation strategies as highlighted hereabove.

Our chain binomial model, originally introduced by Bailey [10], is a so-called discrete-time stochastic alternative to the continuous-time deterministic model based on the health states and transitions presented in Figure 1. The chain binomial model assumes a stochastic version of an epidemic obtained through a succession of discrete generations of infected individuals in a probabilistic manner. Consider a time interval (t, t+h], where h represents the length between two consecutive time points at which we evaluate the model, here h = 1/24 day. Let us assume that there are S_t(k) susceptible individuals at time t in age group k, we expect newly exposed individuals at time t + h, i.e., where I_t is the total number of infected individuals at time t and p_t(k) represents the transmission probability conditional upon contact between a susceptible individual in age group k and an infected individual. The probability that a susceptible individual escapes infection (during a single contact with an infected individual) is equal to q_t(k) = 1 − p_t(k), hence, assuming all contacts to be equally infectious, the escape probability is in case the susceptible individual contacts m infectious individuals. In this setting, the probability of infection for a susceptible individual in age group k = 1, … K can be obtained as: The number of individuals in age group k leaving the exposed state (and entering the pre-symptomatic compartment) within the specified time interval is where 1/γ equals the mean length of the latency period. Probabilistic transitions in the other compartments are derived similarly, hence, a discretized age-structured stochastic model (with step size h = 1/24 days) is fully specified by and Predictions based on the stochastic discrete age-structured epidemic model will account for two sources of variability, namely (1) variability coming from the observational process reflected in uncertainty about the model parameters; and (2) variability introduced by the stochastic process. An overview of the fixed parameter values, sources (incl. literature), and distributional assumptions are listed in Table B3 of Appendix B.

2.4 Next generation matrix and basic reproduction number

The basic reproduction number R₀ for the proposed compartmental model can be obtained by means of the next-generation approach [11]. More specifically, the basic reproduction number is equal to the leading eigenvalue of the next generation matrix, i.e., R₀ is where M ΔV operates by multiplying the ith row of matrix M with the ith element of column vector V. Note that the vector N ≈ S(0) denotes the population age distribution (i.e., the number of individuals in each age group in the population). The time-dependent effective reproduction R_t is obtained by replacing N with the number of susceptible individuals S(t) at time t.

2.5 Intervention measures

Intervention measures mainly targeted the reduction of face-to-face contacts as an effective way of breaking the transmission chains of COVID-19 disease. These measures have led to significant alterations in social mixing patterns, hence, changing the trajectory of the dynamics of COVID-19. To assess the impact of these measures, we utilize social contact matrices derived using the online SOCRATES tool [4], developed for social contact data sharing and assessment of mitigation strategies, and including survey data for various locations, i.e, home, work, school, transport, leisure and other places. The imposed measures have changed the contacts made on these respective contact locations and altered disease transmission. Different social contact matrices are considered to describe the data prior to the lockdown measures and those quantifying contact patterns after the interventions taken. Different choices with regard to the reduction in social contacts are as outlined in Table 1 and Appendix C, and their performance in terms of model fit is compared using the Deviance Information Criterion (see Appendix C for specific details). Note that the choice of the intervention matrix quantifies the extent of social contact reductions, thereby determining the reduction in the effective reproduction number following the installment of stringent lockdown measures [12].

Compliance to the intervention measures taken is assumed to be gradual and is therefore modelled in a flexible way. More specifically, we consider a logistic compliance function where t_I is the time at which the interventions are initiated. The slope parameter associated with the compliance function (i.e., ) is estimated based on the available data.

2.6 Exit strategies

Following the intervention measures that the Belgian government imposed towards limiting the spread of COVID-19 disease, well-tailored exit strategies are needed in order to enable individuals individuals to resume their normal social life whilst protecting the health care system from unprecedented pressure leading to unnecessary loss of lives. Here, we explore and compare possible approaches in lifting imposed measures. The different aspects within the exit strategies are listed below:

• Progressive lifting of lockdown measures on key sectoral pillars of the economy that require physical presence for workers/staff while keeping non-essential service providers closed. This will entail progressively re-adjusting the social contact matrices made at work, during travel/transport, and contacts at other places.

• Gradual re-opening of schools. This will entail re-adjusting the contacts made at school. In line with the exit strategies adopted in Belgium, various partial re-opening dates are considered and their joint impact is explored.

• Opening social places like restaurants, retail stores and hotels. This will involve re-adjusting the social contact matrices for those contacts made during leisure activities, work and transportation as well as those made at other places.

Note that the aforementioned exit strategies cannot be looked at independently since, e.g., parents going back to work will have to rely on childcare/schools to take care of their children. Therefore, we will refer to exit scenarios rather than individual strategies in the remainder of the paper.

To assess the effectiveness of the individual exit strategies and combined scenarios, several comparisons will be made as follows: each exit scenario will be compared to a (baseline) situation without changes, and with each other. More details on the exit scenarios and the translation to the relative number of contacts compared to the pre-pandemic situation are outlined in Table 2. The different scenarios presented in this table give rise to a gradual relief of the lockdown measures taken, similar to the current strategy in Belgium:

Phase 1a – May 4: Although remote work remains the norm, business-to-business services and companies that are able to comply with physical distancing measures re-opened;

Phase 1b – May 11: Shops re-opened under strict requirements related to organisation of the work and restricting access to the store to avoid overcrowding;

Phase 2a – May 18: Schools partially re-opened (first phase - selected grades in primary and secondary schools);

Phase 2b – June 2: Schools partially re-opened further (second phase - pre-primary schools);

Phase 3 – June 8: Restaurants, bars, and cafes re-opened under strict measures including physical distancing and a limited number of customers;

These comparisons were mainly made on the basis of the daily number of new hospitalizations and admissions to the ICU. Furthermore, the implementations thereof in combination with the timing of holiday and school periods will be studied to explore whether rebound effects will occur, i.e., whether deconfinement results in subsequent COVID-19 waves.

2.7 Data and estimation procedure

In this section we briefly describe the parameter estimation procedure and different data sources that are considered to fit the models.

3.1 Baseline scenario accounting for mitigation strategies

Both the probability of experiencing an asymptomatic infection and the probability of having mild symptoms upon contracting COVID-19, i.e., p and ϕ₀ are assumed to be age-dependent. The latter is estimated using a prior distribution based on current literature (see Tables B3 and B4 for more details). The relative infectiousness of asymptomatic versus symptomatic individuals r_β is fixed at a value of 0.51 [19]. Other model parameters are either fixed or estimated based on the available data (see Table B3 and Appendix F).

The best fitting model included social contact matrices with an 80% reduction of the normal work and transportation contacts (α = 0.2 in Appendix C), with no school contacts and with 10% of the regular leisure contacts and contacts related to other activities. In Figure 2, we graphically depict 25 stochastic realizations of the hospital admissions and deaths since March 1 based on a thinned chain from the joint posterior distribution of the model parameters together with pointwise 95% credible intervals derived from stochastic realizations based on 5000 random draws from the joint posterior distribution of the model parameters. The figure clearly shows that the observed daily number of hospitalizations and deaths (black dots) are well described by the model. Furthermore, the estimated age-dependent daily numbers of new hospitalizations and deaths are graphically depicted in Figure 3 for the 10 age categories.

• Download figure
• Open in new tab

Figure 2:

Stochastic realizations of the compartmental model based on a thinned MCMC chain from the joint posterior distribution of the model parameters and relying on an ‘asymptomatic’ and ‘symptomatic’ social contact matrix composed of 20% of regular work and transportation contacts, no school contacts and 10% of leisure contacts and contacts related to other activities. Shaded areas represent 95% credible intervals. Reported daily number of hospitalizations and deaths are represented by black circles.

• Download figure
• Open in new tab

Figure 3:

Stochastic realizations of the compartmental model based on a thinned MCMC chain from the joint posterior distribution of the model parameters and relying on an ‘asymptomatic’ and ‘symptomatic’ social contact matrix composed of 20% of regular work and transportation contacts, no school contacts and 10% of leisure contacts and contacts related to other activities. Number of new hospitalizations (left upper and lower panels) and deaths (right upper and lower panels) are shown for all 10 age groups. Shaded areas represent 95% credible intervals. Reported daily number of hospitalizations and deaths are represented by circles.

3.2 Posterior distributions

In Table 3, we present summary measures for the posterior distributions of the most important (implicit) model parameters including the posterior mean, median, standard deviation and 95% credible intervals (CIs). An overview of posterior quantities for all model parameters is included in Appendix F).

The basic reproduction number R₀ at the start of the epidemic - prior to any government intervention - is estimated to be 2.900 with 95% CI (2.885, 2.918). On May 4, 2020 the effective reproduction number after the lockdown was estimated to be 0.738 (95% CI: 0.732, 0.744). The time-dependent effective reproduction number R_t is shown in Figure F3 of Appendix F. The posterior means for the parameters ω⁽ − 1)(k) are in line with estimates of the median duration between symptom onset and hospitalization [20]. Furthermore, the average length of the incubation and asymptomatic infectious period, ρ⁻¹ and are estimated to be 3.481 days (95% CI: 3.304, 3.642) and 6.284 days (95% CI: 5.959, 6.555), respectively, very similar to values reported in the literature (see Table B3). For individuals experiencing severe symptoms, the average length of the infectious period, constrained due to isolation in hospital, depends on the age-specific time between symptom onset and hospitalization. Full compliance to the intervention measures was obtained after approximately 6 days (see Appendix F).

Boxplots of the marginal posterior distributions of the probability of hospitalization are presented in Figure 4. The probability of hospitalization increases substantially with increasing age. A decrease in hospitalization probability is observed in age class [80, 90) after which it increases again for individuals of age 90+. However, the probability of hospitalization, as a proxy of disease severity, is likely time-dependent as well as biased for the oldest age groups due to differential referral policy in elderly homes or end of life choices in the last will of severely ill persons. More specifically, the general population essentially consists of two subpopulations, i.e., a nursing home and non-nursing home population. It is likely that in reality, the nursing home residents were more often exposed than elderly in the general population of the same age. One other explanation for the decrease in probability of hospitalization is that nursing home residents (which constitute an important part of the age group between 80 and 89 years of age) were less likely referred to the hospital (as would have been the case when they were not living in a nursing home), at least during the initial part of the epidemic, thereby implying a lower hospitalization probability in that age group. Moreover, persons in age class [80, 90) living in the non-nursing home population are believed to be in better condition than nursing home residents of that age, being more frail when suffering from more severe comorbidities. Hence, this could lower the probability of hospitalization further. Estimated mortality and infection fatality rates are shown in Figure F4 in Appendix F.

• Download figure
• Open in new tab

Figure 4:

Boxplots of the marginal posterior distributions of the probability of hospitalization by age group.

3.3 Estimated age- and time-dependent (sero)prevalence of COVID-19

In Figure 5, we show the estimated and observed seroprevalence on March 30, 2020 (left panel) and on April 20, 2020 (right panel) with 95% credible intervals in gray dashed lines and asymptotic 95% error bars for the weighted seroprevalence.

• Download figure
• Open in new tab

Figure 5:

Estimated age-dependent seroprevalence of COVID-19 with 95% credible interval on March 30, 2020 (left panel) and April 20, 2020 (right panel). Observed seroprevalences are depicted using red dots with 95% confidence intervals in blue. The confidence interval for the age group [0, 10) is wide due to the low number of individuals (n = 36).

The estimated age- and time-dependent prevalence of COVID-19 in the Belgian population is shown in Figure 6. Clearly, the estimated prevalence is higher in the oldest age groups ([80,90), 90+) which is in line with the observed seroprevalence, the latter providing a cross-sectional snapshot of the (delayed) build-up of seropositivity in the population upon infection (see Figure 5). The estimated overall weighted prevalence of COVID-19 is equal to 0.069 (95% CI: 0.064, 0.073) on May 4, 2020 (black dashed line with gray shaded area).

• Download figure
• Open in new tab

Figure 6:

Estimated time-dependent prevalence of COVID-19 in the different age groups and its weighted average (right panel; black dashed line - right y-axis) based on 5000 stochastic realizations given random draws from the joint posterior distribution of the model parameters with 95% credible intervals (shaded regions).

3.4 Exit strategies

We performed short-term predictions of the dynamics of COVID-19 in the Belgian population through the use of scenario analyses. In Figure 7, we show the impact of relaxing lockdown measures by increasing contacts made by people at different locations. More specifically, we present stochastic realizations of the model to predict the number of new hospitalizations under different scenarios with changes in contact behaviour as of May 4, 2020 (Phase 1a, Scenarios S1–S3) and May 11, 2020 (Phase 1b). From May 4 onward, we presume that work-related contacts will increase from 20% (baseline scenario - see Table 2) to 30% of the pre-pandemic number of contacts at work (or transmissibility is reduced to an extent equivalent with the assumed reduction in work-related contacts), and that the number of transport contacts and contacts during leisure and other activities will increase respectively from 20% to 30% and from 10% to 20% of their pre-pandemic values. On top of that, work, transport and leisure contacts stay constant (scenario S1) or increase to 40%, 40% and 30% (S2) or 50%, 50% and 40% (S3) of the pre-pandemic values, respectively, from May 11 onward. Moreover, a delay of one week is considered for each change in social contact behaviour (i.e., a full extent of all changes in behaviour is reached on May 18, 2020). A small to moderate increase in the contacts at work, transportation and during leisure (blue and purple scenarios) leads either to a complete reduction of hospitalizations or a constant number of new hospitalizations over time. Only the most extreme increases in contacts give rise to a resurgence of COVID-19 implying a second wave of COVID-19 infections (scenario S3 - orange lines).

• Download figure
• Open in new tab

Figure 7:

Impact of various exit strategies in terms of the number of work- and leisure-related contacts on the number of new hospitalizations in the absence of re-opening of schools.

A partial re-opening of schools as of May 18 (Phase 2a) is studied in detail in Figure 8. Work- and transport-related contacts and contacts during leisure and other activities increase as of May 4. School-related contacts are assumed to be 20% (S4 - blue lines), 40% (S5 - purple lines) or 60% (S6 - orange lines) of such contacts prior to the epidemic. This increase in school-related contacts is imposed between May 18, 2020 and July 1, 2020. The start of the summer holiday on July 1, 2020 implies a reduction of all school-related contact to 0%. A partial re-opening of schools in combination with a moderate increase in work, transportation and leisure activities leads to a small to moderate increase in the number of new hospitalizations after lockdown measures are relaxed.

• Download figure
• Open in new tab

Figure 8:

Impact of partial re-opening of schools on the number of new hospitalizations.

Finally, we investigate long-term predictions of subsequent COVID-19 waves for a selection of possible exit scenarios (Figure 9). In those scenarios, we mimic the timing of the Belgian exit strategy. More specifically, schools are partially re-opened on May 18 and June 2, 2020 yielding 20% of school contacts as of May, 18 (S7), an increase from 20% to 40% or 60% of school contacts between May 18 and June 2 for scenarios S8 and S9, respectively. Schools are closed during the vacation period starting from July 1, 2020 until August 31, 2020. We assume that contact behaviour at schools following partial re-opening on September 1, 2020 is equivalent to 60% of the pre-pandemic social contacts made at school. In the lower panel of Figure 9, scenarios S10–S12 show the impact of an increase of leisure contacts to 20%, 40% or 60% of pre-pandemic leisure contacts as of June 8, while assuming school-related contacts to be equal to 20% upon Phases 2a and 2b. Under the assumption of unadapted behaviour given a contact, we observe that due to an insufficient depletion of susceptibles during a second wave of COVID-19 infections (or a phase with a stable daily number of new hospitalizations) a large increase in the number of new hospitalizations will occur by the end of the year with a higher peak size if the one of the second wave (or the plateau level) was lower. The cumulative number of hospitalizations over time is presented in Appendix F. Moreover, leisure contacts are important in determining the peak size of the wave at the end of the year (lower panel of Figure 9). In Figure 10, boxplots of the estimated prevalence over time is shown for scenarios S7–S9 and age groups [0, 10), [30, 40), [60, 70), 90+. The largest increase in prevalence between May 1 (baseline) and December 1, 2020 is observed in the highest age category with an average increase ranging between 36.5% and 38.4% across different scenarios. In all age categories, the increase in prevalence is smallest for scenario S7 and highest for scenario S9.

• Download figure
• Open in new tab

Figure 9:

Long-term predictions of the impact of various exit strategies on the number of new hospitalizations.

• Download figure
• Open in new tab

Figure 10:

Predictions of the prevalence in exit scenarios S7–S9 for age group [0, 10) (top row), [30, 40) (row 2), [60, 70) (row 3) and 90+ (bottom row). Increments in prevalence compared to the prevalence on May 1, 2020 is added on top of the boxplots.

3.5 Validation of the model

Validation of the model is done based on (1) data on new hospitalizations and deaths following the relaxation of the lockdown measures, (2) serological survey data collected in a third round and (3) infection fatality derived from Belgian mortality data [21]. In Figure 11, we show stochastic realizations under the baseline scenario (without change in social contact behaviour after lockdown measures are relaxed) overlaid with new data points after May 4, 2020 (black solid circles). In general, the stochastic model describes the observed data well, even in the absence of changes in contact behaviour after intervention measures were relaxed. Although the observed number of hospitalizations tends to remain constant, thereby deviating from a further decrease thereof in the baseline scenario, no large differences between observed and predicted values are found. Following the gradual relief of the intervention measures, no resurgence of the disease is noticeable to date (end of June).

• Download figure
• Open in new tab

Figure 11:

Stochastic realizations based on the baseline scenario (without change in contact behaviour upon relaxing the stringent lockdown measures) for the number of new hospitalizations (left panel) and the number of new deaths (right panel) together with observed data points used for fitting (black open circles) and for validation (black solid circles).

Next, the estimated seroprevalence based on the model is related to the one obtained from a third cross-sectional serological survey. Data was collected in a third round, after the initiation of the gradual relaxation of the stringent measures, between May 18, 2020 and May 25, 2020. The overall weighted seroprevalence was estimated to be 6.87% (95% confidence interval: 5.89%, 8.01%) [5]. In our model, the posterior mean of the seroprevalence is 6.8% with 95% credible interval (6.4%, 7.2%) which is similar to the aforementioned values. Furthermore, the estimated infection fatality rates (IFRs) (see Appendix F) are in line with those reported by Molenberghs et al. (2020) [21]. These authors report an overall IFR of 0.43% (95% confidence interval: 0.30%, 0.62%) in the non-nursing home population whereas our model suggests a posterior mean of 0.507% (95% CI: 0.480%, 0.536%) which is nicely in line. Age-specific IFRs are presented in Appendix F.