Most of the epidemic models divide the target population into a certain number of compartments, consisting of individuals with identical statuses concerning a given disease. The foundations of the entire approach to epidemiology based on compartmental models were laid by public health physicians in the early 1900s. One of the first applications of the compartmental model was made by R. Ross, who demonstrated the dynamics of the transmission of malaria between mosquitoes and humans and consequently was awarded the Nobel Prize in Medicine in 1902 (7, 8). Since then, compartmental models are still widely used to simulate the spread of a variety of infections (9).

One of the most popular extensions of the SIR model is the SEIR model (10), a traditional method used to simulate infectious disease that incubates inside the hosts for a while before the hosts become infectious. The SEIR model considers the incubation period by introducing a new compartment E (Exposed) to the compartmental system. This model and its modifications were already adapted to simulate the COVID-19 virus in many countries (11–13). In this work, we adopt a widely-used modification of the model—SEIRD (14) with the death compartment D. More recently, a SEIRD model with relaxed parameters has also been proposed to consider the rapidly changing social scenario arising from the period of the COVID-19 (15). Likewise, ameliorating compartment models according to the scenario would be a promising direction to study.

The Markov chain Monte Carlo (MCMC) is a large class of sampling algorithms widely used for probabilistic problems. MCMC was first introduced in 1953 as a new method to simulate the distribution of states for the system of idealized molecules (16). However, the application of the algorithm did not limit itself to the physics field. It was later adapted and generalized by Hastings (17) to focus on statistical problems, opening its application to a wide range of domains. Due to its ability to handle complex types of analyses, the MCMC approach was widely used in finance (18, 19), communication (20, 21), computational biology (22), linguistics (23, 24), and other fields with probabilistic settings. By no surprise, these methods are widely popular for estimating effects in complex epidemiological analyses as well (25–27). For example, Cauchemez et al. (28) has shown how to model influenza transmission using the Bayesian MCMC approach, and lots of variations of MCMC methods were used to infer features of an Ebola virus and analyze its transmission mechanism (29, 30). Recent reports are also benefited from the Bayesian MCMC methods to infer COVID-19 virus transmission dynamics. Zhou et al. (31) implemented such inference based on a probabilistic compartmental model using daily confirmed COVID-19 cases and applied it to six states of the United States.

MCMC algorithms are also successfully applied to change-point models. The objective is to detect the abrupt property changes lying behind the time-series data (32). Recent work showed that MCMC algorithms with Bayesian parameter inference could be used to detect change-points in COVID-19 spread using SIR and SEIR epidemiological models of South Africa (33). According to their results, South Africa experienced two change-points: the first at the time of the national lockdown and the second after the massive screening and testing program. Dehning et al. (34) adopted a similar approach to the case study of coronavirus spread in Germany by utilizing SIR models with MCMC sampling for detecting change-points in effective growth rate that correlates well with the times of publicly announced interventions. Following these examples, we applied the change point model with an extended SEIRD epidemiology model to identify where policies affected the COVID-19 virus transmission rate in 10 countries with different interventions.

A wide range of work was done to estimate the efficiencies of the policies imposed by different countries to prevent COVID-19. Many of them were focused on the individual country cases considering their unique demographic features (35–37), while other reports compared many countries by the independent effects of a single category of policy (38–40). For instance, Iwata et al. (41) used Bayesian method analysis. They did not reveal the effectiveness of school closures that occurred in Japan in mitigating the risk of coronavirus infection in the nation. Another recent work by Sharov et al. (42) used a modified SIR model to compare the effectiveness of lockdown measures introduced during the coronavirus pandemic in 13 European countries, comparing them to two baseline countries (Sweden and Iceland) that did not implement the lockdown policies. For evaluation, this work used the herd immunity level and time of formation to indicate the effectiveness of lockdown measures (42). According to Sharov's results, lockdown and no-lockdown modes of containment led to roughly similar results.

There are also reports considering multiple policies across the globe (43, 44). One example is work by Flaxman et al. (45), which investigates effects of applied non-pharmaceutical interventions (NPIs) across 11 European countries for the period from the start of the COVID-19 epidemics. According to their results, major non-pharmaceutical interventions, like lockdowns, have had a significant effect on reducing the transmission of the virus. However, a subsequent study by Haug et al. (46), which assessed the efficacy of 6,068 NPIs across 226 countries and gave a detailed analysis of the country-specific “what-if” scenarios, showed different results. They analyzed the impact of government interventions on the effective reproduction number R_t by combining several analytical approaches. By utilizing statistical, inference, and artificial intelligence tools, they concluded that combinations of some less disruptive and less costly NPIs could be as effective as more expensive and harsh ones like national lockdowns. Brauner et al. (47) came to the same conclusion by analyzing 41 countries during the first wave of the pandemic. According to their study, less harsh NPIs can be more effective in mitigating COVID-19 transmission than more strict stay-at-home orders (47). Singh et al. (48) exploited the spatial and temporal variation in the introduction and lifting of non-pharmaceutical interventions (NPIs) across counties using a staggered difference-in-differences (DID) approach. They compared US counties with NPIs in place (treated) with counties that do not have NPIs in place (control) before and after implementation. Enabled by datasets with rich population characteristics, they stratified the datasets into several groups and analyzed the impact of implementing and lifting NPIs by the population groups they target. However, as we will discuss in Section 5, it takes a certain amount of lagging time to see the effects of NPI implementations. More meaningful analysis can be obtained with DID by considering the delay.

However, one of the significant limitations of recent studies is that none of them perform a comprehensive analysis considering the economic factors that affect the efficiencies of the policies. There were some reports regarding the economic cost of the pandemic situation across the globe (49). For example, McKibbin and Fernando (50) simulated a global economic model to explore seven scenarios, which differ in the proportion of the population who become infected or dead. According to their estimations, in a scenario where COVID-19 develops into a global pandemic, the cost of lost economic output begins to escalate into trillions of dollars (51). However, they do not include the effect of policy interventions in their simulations. To address this limitation, we propose another method of cost estimations by involving policy effects in Section 6.

The most basic compartmental model, namely the SIR model, uses three compartments of Susceptible (S), Infectious (I), and Recovered (R). Each individual can move from a compartment to another compartment, resembling the progress of the disease. We could use S, I, and R to denote the number of individuals in their respective compartments. For the COVID-19 case, there is an incubation period in which people are infected but not yet infectious. Hence, we adapted the epidemiological SEIRD model (14) for our simulations, extending the SIR model with the E compartment of exposed individuals and the D compartment for deaths.

The below differential equations (Equations 1–5) describe the transition between compartments (14).

where the effective reproduction number (R_e) is the expected number of people that each infected individual can transmit the virus to during the outbreak, the basic reproduction number (R₀) is the natural reproduction number when there is no intervention, the incubation time (tE=1σ) is the average time in which an individual is exposed but not yet infectious, the recovery time (tI=1γ) is the average time after which an the infected case become concluded (recovered/dead), the case fatality proportion μ is the proportion of fatal cases among all concluded cases, and the waning time (tR=1α) is the time that recovered individuals retain immunity. Since we deal with the initial stages of the pandemic, we assume people carry immunity to the disease upon recovery (α = 0) and the population stays constant over time and is equal to N.

Substantial amounts of COVID-19 cases are not reported due to testing availability and testing strategy. Hence, the exact total number of COVID-19 cases is unknown and typically not uniquely determined from the number of confirmed cases (52–54). Korolev (53) emphasizes that neglecting unreported cases leads to biased parameter estimation, so it is important for our model to address this issue.

To distinguish observable cases among all possible cases, we use a parameter called response rate ρ, which is the probability of a case being reported. At each timestamp (in our case, each day), if the SEIRD model estimates the number of new cases (transition from S to E) to be S2E and the number of recovered cases (transition from I to R) to be I2R, our model, the number of reported confirmed cases and recovered cases will be corrected by the response rate:

where Binomial(n, p) is the discrete probability distribution of the number of successes in a sequence of n experiments. Here, a reported case can be understood as a success, and each case can be understood as an experiment. Response rate ρ is inferred together as a parameter for our SEIRD model. Since it varies widely across countries, the estimated value is used internally within the country and is not generalized for others.

To implement the probabilistic models, we used the probabilistic programming language Pyro (55). For this particular inference task, we adopted Pyro's Epidemiology framework (56) for scaling up our experiments with a restricted class of stochastic discrete-time discrete-count compartmental models. This framework uses the Markov Chain Monte Carlo (MCMC) algorithm to fit the SEIRD model to infer COVID-19-related parameters: reproduction number R₀, recovery time, incubation time, transmission rate, and mortality rate.

We ran the model to estimate the parameters of the Swedish data before April 1st, 2020. We chose this early stage of the COVID-19 pandemic since Sweden did not impose any strict policies and aimed to achieve herd immunity (62). We assumed that the virus transmission rate was unaffected by any interventions, so we used Sweden as a baseline case and perform the experiments to infer unaffected COVID-19 parameters. To run our probabilistic model, we set the prior according to the estimations of the World Health Organization (63). It was reported that mild cases typically recovered within two weeks, the incubation period was on average 5–6 days, and R₀ was typically around 2. Mortality and recovery rates differed depending on the region and stage of the virus spread, but the case-fatality rate was roughly 2.5%. Prior of the response rate ρ is given as Beta(10, 10), which favors the initial value around 0.5, then converges in the range of 0 to 1.

The obtained posterior virus-related statistics for Swedish data are shown in Table 1. For more accurate results, we ran the model six times and reported the averaged values. The results are reasonable enough to use in our further simulations.

There are several works that explain and compare the transmission of viruses under various situations by using the slope of the log-transformed instances. Both the rate of change of the log-transformed case incidence and the instantaneous reproduction number, R_e, are considered important for the investigation of the virus (64). For example, Caspi et al. (65) modeled replication rate (RR) as the slope of the logarithmic curve of confirmed cases to compare the coronavirus spread in different climates. Another study by Gebski et al. (66) observed changes in the slopes of log-transformed incidents of Staphylococcus aureus (MRSA) infections in hospitals to evaluate the success of interventions. In this work, we adopt similar strategy by evaluating effectiveness of the interventions by comparing the slopes of the confirmed cases before and after the policy establishment.

We will first explain how to measure the policy strength by referring to the SEIRD compartmental model (see Section 4.1) with some simple assumptions and new terms we describe below. The transmission rate β is the number of susceptible individuals that an infected individual can infect in a day, which is calculated as β = R_eγ. At the very beginning, almost everyone in our setting is in the Susceptible compartment, so we can assume that S is equal to the total population size N, or S = N. With this approximation, Equation (1) can be rewritten as follows:

Additionally, since the incubation period is much shorter than the recovery time, we can ignore the E, R, and D compartments at the initial stage of the simulation. Thus, we can approximate the total population size as N = S+I and Equation (3) becomes:

The size of the Infected compartment rises exponentially with the rate w = β (that is, on day t, the number of infected cases is calculated as e^βt). Due to the exponential nature, it is appropriate to investigate the case counts using a log scale. In the log scale, the exponential spread is represented as a linear line, with the transmission rate of the virus β represented by the slope w.

Consider an example in Figure 4. In the beginning, the graph is a steep line with the slope w₁ (representing a rapid, exponential spread). After a corresponding policy is applied, the graph bends and becomes less steep with the slope w₂≪w₁ (slower spread). Therefore, the graph roughly consists of two lines of different slopes, w₁ and w₂, with a separation point in-between, which we call a change-point (the black dotted vertical line in the graph). The slope w₁ and w₂ represent the transmission rates before and after the change-point. Since w₁ and w₂ represent the transmission rates before and after the policy takes its effect, we can define the strength of the target intervention in terms of their ratio:

Given the incubation period, we expect the policy will show effect after around 2–4 weeks after the policy establishment. In the next subsection, we will introduce a probabilistic programming approach to find the change-point when the policy takes effect.

As was mentioned in the related work (see Section 2.2), probabilistic models were successfully used to detect change-points in transmission rates of coronavirus. In the present subsection, we describe the probabilistic programming approach to detect change-points as well as estimating policy strength. Algorithm 2 sketches how the change-point detection model based on the probabilistic programming operates. From the case trajectory of the COVID-19 pandemic, we used the HMC algorithm with NUTS (58) to estimate the change-point τ and slopes w₁ and w₂.

We investigate major initial interventions applied by several countries to mitigate the virus spread. For more accurate results, we chose nine countries presented in Table 2 that strongly imposed corresponding policies, assuming that they were applied to the fullest extend. By investigating the countries that applied a policy most stringently, we find a meaningful upper bound for each policy's efficacy. This upper bound is helpful for policymakers to determine the most appropriate intensity of the policy (more details in Section 6). In this experiment, we focused on five main policy categories:

By using the probabilistic programming model described in the previous section, we detected the amount of time the policy needed to take effect after establishment, the change-point, and policy strength for each case. Since we focused on the initial stage of the pandemic, the time frame we simulated is 3–6 months from the first recorded case, including the date that the policy is enforced and its effect could be seen. In all experiments, results have converged to the values consistent with our priors. The posteriors also fit well with the actual data (Figures 5–12). The summary of policy efficiencies is shown in Table 2. In the next section, we discuss the results in more detail by investigating each policy separately by country.

To simulate the infection and fatality cases, we used the SEIRD (see Equation 5). We follow the differential equations Equations (1)–(6) and the virus and policy statistic derived in Sections 4, 5. However, the fatality rate will not stay constant as we considered the hospital capacity.

Apart from the parameters related to virus statistics and policy efficiency, we also need the input of the economic effect of the policies as well as the hospital statistics (hospital capacity, percentage of cases requiring hospital admission, and the death rates with and without hospital treatment). Our model gives flexibility to each country to input its own parameters into the model.

We illustrate the operation of our model on an imaginary country with a population of 1 million. In addition to the parameters we inferred from previous results, we applied some additional assumptions:

We also estimate the economic and human capital cost for each policy:

These estimations are reasonably set based on the following facts:

We ran the model without any policy and with the lockdown applied from day 30 to 60. As you can see from Figure 14, applying lockdown for 1 month simply postpones the virus spread. Another problem is that it significantly hits the country's economy, so it cannot be applied for a long time. Thus, even though lockdown is estimated to have the highest efficiency of 0.96, it might not be the best policy to apply. So further experiments are required to identify how, when, and for how long the policies should be applied.

Finally, we want to devise the best initial response to the virus. Using the inferred statistics, we conducted experiments on the policies and performed simulations to develop the optimal policy with minimal loss (both economic loss and life loss). We designed an imaginary country with a population of 1 million and a GDP of $30,000 per capita. The country had a population of 1 million, and the simulation spanned three months. We assume that policy-makers revised the policy every month, and a policy is applied exhaustively, partially (50% efficacy), or not applied at all. Policies could be applied together. The goal was to minimize the cost.

Due to the challenge of separating the policy effects, our study has some limitations within which our findings need to be interpreted carefully. First, we have investigated major policies applied in relatively wealthy countries. For economic loss estimations, we utilized policies' costs reported from various developed countries using different sources. Depending on the parameters adjusted for a particular country, the results might be different. Second, the estimated efficiency of each policy in Section 5 is measured on the most successful cases. Although it provides a meaningful upper bound for each policy's efficiency, we cannot simply assume that every country can achieve this maximum efficiency. Therefore, the “best” initial response in Section 6.4 should be understood under the context that every policy can be feasibly applied to its fullest extent. Third, future work can utilize our model and extend it by considering the confounding effects of other interventions or changes in the susceptible population size. Such a model can be used to estimate the efficiencies of the policies applied in the later stages of the pandemic.

Given the unique socioeconomic state, each country has its own feasible stringency and the price tag for each measure. Many factors contribute to this variability, such as public acceptance, political climate, and government priority. To cope with that, we provide a flexible end-to-end pipeline, which can be tailored to each country's specific needs. Decision-makers can adjust the corresponding parameters or apply their country's cost estimation to adapt to their own situation. They can also exclude infeasible policy settings in their country when running the simulation. The resilience to control model's parameters allows countries to see how the pandemic will play out under different scenarios and build their own strategies based on the model's output.