A probabilistic epidemiological model for infectious diseases: The case of COVID-19 at global-level.

This study has developed a probabilistic epidemiological model a few weeks after the World Health Organization declared COVID-19 a pandemic (based on the little data available at that time). The aim was to assess relative risks for future scenarios and evaluate the effectiveness of different management actions for 1 year ahead. We quantified, categorized, and ranked the risks for scenarios such as business as usual, and moderate and strong mitigation. We estimated that, in the absence of interventions, COVID-19 would have a 100% risk of explosion (i.e., more than 25% infections in the world population) and 34% (2.6 billion) of the world population would have been infected until the end of simulation. We analyzed the suitability of model scenarios by comparing actual values against estimated values for the first 6 weeks of the simulation period. The results proved to be more suitable with a business-as-usual scenario in Asia and moderate mitigation in the other continents. If everything went on like this, we would have 55% risk of explosion and 22% (1.7 billion) of the world population would have been infected. Strong mitigation actions in all continents could reduce these numbers to, 7% and 3% (223 million), respectively. Although the results were based on the data available in March 2020, both the model and probabilistic approach proved to be practicable and could be a basis for risk assessment in future pandemic episodes with unknown virus, especially in the early stages, when data and literature are scarce.

INTRODUCTION

The World Health Organization (WHO) has declared the coronavirus disease (COVID‐19), which is caused by the severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), a pandemic (Livingston et al., 2020a). Indeed, it has been confirmed that there are infected people in all the six continents, that is, Europe (EU), North America (NA) (Canada, the United States, and Mexico), South America (SA) (all other countries in America including Central America), Asia (AS), Africa (AF), and Oceania (OC). It is a worldwide threat, and all countries and continents must come together as one against this common enemy.

For the purposes of this work, risk is a combined measure of the probability and consequence of an undesired event (Duarte & Droguett, 2016; Inaba & Sekine, 2004; Zhao & Milner, 2008). Quantitative risk assessment (QRA) is the formal and systematic approach to quantify the risks associated with a system (Duarte et al., 2019; Forbes et al., 2011, 2010; Harwell et al., 2012). Here, the system is an infectious disease spreading across the world and the undesired event is a high number of infections within a period of time. More specifically, we focus on the quantification and assessment of what we call the “risk of explosion” (i.e., the number of infected people increasing over 25% of the world population within 12 months).

Many strategies have been discussed and implemented to control the spread of the virus until a vaccine is developed, licensed, and manufactured on a global scale. These actions include restrictions on travel and business/studies/social activities (from now on, the term business will be used to refer to all three types of activity), social isolation (for this work, this is equivalent to “stay at home” measures), vertical isolation (i.e., when it affects only the elderly and groups at risk), and using therapeutics and new medical tools to reduce fatality rates (from now on, the term “medical tools” will be adopted and this includes the use of therapeutics).

Most decisions have been taken based on the subjective opinions of epidemiologists (Kupferschmidt & Cohen, 2020; Sohrabi et al., 2020; Toms & Petrie, 2020) or projections of deterministic models (Choi & Ki, 2020; Ferguson, Laydon et al., 2020; Ferguson, Walker et al., 2020; Khan & Atangana, 2020; T. M. Chen et al., 2020; Yang et al., 2021), but neither the risks associated to these decisions nor the uncertainties in the opinions and estimates have been quantified. In this sense, a recent study highlights the importance of acknowledging uncertainty as a central risk component to properly characterize and communicate risk (Aven & Bouder, 2020). Globally, policymakers are demanding tools to guide them on how to prioritize resources for designing control strategies. Indeed, they require objective answers for questions such as:

How many people may die globally, and how many may be infected if we decide on a hypothetical mitigation strategy (say A, B, or C) for the next year?

By what amount and for how long is the social isolation of nonelderly and elderly people necessary to reduce risk to a low or negligible level?

To what extent and for how long are restrictions on intercontinental flights necessary to reduce risk to a low or negligible level?

Which continents are at risk in the future? Which one deserves the most effort to control the disease? What is the order of prioritization?

Generally, some simplifying assumptions are necessary to model the dynamics of a disease, such as those involved in the transmission of SARS‐CoV‐2. For example, traditional approaches for modeling COVID‐19 are based on deterministic models that often rely on average data, and thus only provide expected results. These neither propagate the variability and uncertainty of data nor consider environmental stochasticity (i.e., the unpredictable natural fluctuation in vital rates; Fujiwara & Takada, 2017). As deterministic models for SARS‐CoV‐2 (Choi & Ki, 2020; Ferguson, Laydon et al., 2020; Ferguson, Walker et al., 2020; Khan & Atangana, 2020; T. M. Chen et al., 2020; Yang et al., 2021) lack stochasticity in parameters and uncertainty in results, potentially misleading conclusions may lead to imprudent decisions, which in turn, might lead to a much greater number of lives being threatened and lost.

In fact, one probabilistic model for SARS‐CoV‐2 was used for assessing the risk of outbreaks outside China (Boldog et al., 2020). This model made predictions from January 23rd until March 15th, 2020, and evaluated control measures (e.g., travel restrictions from China) for countries at risk. Severe limitations of this study include: (i) it was conducted at the beginning of February 2020 when, supposedly, most infected people were still in China, and thus it could simulate the spread of the disease only from China; indeed, the disease was not yet even considered a pandemic (ii) it does not separate age groups; (iii) it does not quantify the number of deaths, but only cases of infection; (iv) it does not consider social isolation as a control measure.

Our model is somewhat similar to the susceptible–infected–recovered (SIR) type models (e.g., Weiss, 2013). However, instead of using an existing SIR model and adapting it to the purposes of risk assessment at the global level, our model was built from its origin using RAMAS® Metapop v. 6.0 software (Akçakaya & Root, 2013), which is a computational tool for constructing population‐based solutions (most often ecological models, e.g., those by Andersen et al., 2017; Bino et al., 2020; McCusker et al., 2017), but also epidemiological approaches (e.g., Duarte et al., 2014) that rely on probabilistic simulation via the Monte Carlo method (Kalos & Whitlock, 2008). Thus, besides being stochastic, building an original model gave us the flexibility to include new features such as each SIR state is subdivided into age classes (nonelderly and elderly) with different probabilities of fatality and infection; the spatial structure of the infected population at the global level over time with potential for dispersal among continents; and density‐dependence (DD) effect to account for reductions in infection rates as the susceptible population density decreases.

In 2021, some new deterministic and stochastic models have been proposed. Regarding the deterministic ones (Carcione et al., 2020; Kumar et al., 2021; Lawal & Vincent, 2021; Mandal et al., 2021; Piccirillo, 2021; Varghese et al., 2021), the main difference, when compared to past studies (Choi & Ki, 2020; Ferguson, Laydon et al., 2020; Ferguson, Walker et al., 2020; Khan & Atangana, 2020; T. M. Chen et al., 2020; Yang et al., 2021), is the compartmental susceptible–exposed–infective–recovered (SEIR) approach (for more details see Bartlett, 1957). Besides being deterministic, these are not age‐structured, as we propose here. Regarding stochastic models offered in 2021 (Abrams et al., 2021; Barreiro et al., 2021; Engbert et al., 2020; Siqueira et al., 2021), all of them focused on a single specific country or region.

Another characteristic of our model is that it represents COVID‐19 globally. All the models mentioned above described COVID‐19 dynamics in a single country or city. To the best of the authors’ knowledge, Ferguson, Walker et al. (2020) provided the only model that predicts the global impact of COVID‐19 and evaluates strategies for mitigation and suppression. Indeed, Ferguson, Walker et al. (2020) deterministically predicted the number of infections and deaths in the world in seven different regions for 250 days ahead since March 23, 2020 for varying scenarios. However, Ferguson, Walker et al.’s approach (2020) cannot be adopted for the purposes of risk assessment because they used only single‐point estimates for the number of infections and deaths. Thus, such an approach does not tackle the inherent uncertainty of these predictions (EPA, 1998).

On the other hand, probabilistic models can consider uncertainty in parameters and give risk results in terms of probabilities. Using Monte Carlo simulations, for example, allows us to obtain a set of results for infections/deaths so that a nonparametric probability distribution can be drawn, which accommodates uncertainty, that is, there is a probability associated with each estimate of infections/deaths at time t. With such an approach, a model can successfully quantify risks to measure the likelihood of undesired consequences (infections/deaths), thus, supporting decisionmakers to understand the likelihood of the outcomes of a given measure and making informed decisions.

Therefore, this paper aims to develop an epidemiological model for COVID‐19 at the world‐population level that overcomes the drawbacks mentioned above and is tailored for QRA of hypothetical mitigation scenarios as a decision support tool for more effective planning. We set out to answer the questions mentioned above. To the best of our knowledge, this paper conducts the first probabilistic epidemiological modeling of COVID‐19 at the worldwide population level. It is worth mentioning that parameters estimates that we used to feed our model were based on the evidence available by March 2020, when this article was written. More sophisticated knowledge has been acquired since then. Therefore, our model does not attempt to make precise predictions, but instead it is only descriptive in a sense that it is a tool for describing the dynamics of COVID‐19 under predefined scenarios (different conditions of social isolation, travel restrictions, medical devices) to evaluate the role of such conditions and produce meaningful conclusions that can be used to steer public health decisions. Hence, the model is meaningful for decisions made under uncertainty, but due care must be taken when interpreting results.

The remainder of this paper is organized as follows. First, we present the structure of the model and assumptions, which is flexible in parameterization and can be used to simulate varying scenarios. Next, we discuss the materials and methods used to conduct a risk assessment employing this model. We then present results for representative scenarios, validate the model, and discuss its advantages and limitations. Finally, some conclusions are drawn, and suggestions are made for future research.

THE STRUCTURE OF THE MODEL AND ASSUMPTIONS

Our proposed model is probabilistic and, thus, provides meaningful information to decisionmakers because it allows for: (i) the assessment of uncertainty by specifying lower and upper bounds in the results; (ii) modeling the spatial dynamics of infected people across six continents; (iii) the quantification of the risk of explosion (i.e., the number of infected people increases to over 25% of the world population after 12 months). The 25% threshold for the risk of explosion was defined based on the most severe pandemic in the 20th century, namely the Spanish Flu, which lasted from January 1918 to December 1920 (caused by the A(H1N1) virus). It is estimated to have infected 500 million people, approximately 25% of the world population at that time, and caused 20–50 million deaths (Spreeuwenberg et al., 2018; World Health Organization, 2005).

These tasks may be performed for scenarios regarding different containment measures, thus assessing their effectiveness in terms of risk. Using this approach, it is also possible to identify the continents where SARS‐CoV‐2 might persist, and hence this may help target public control strategies to reduce human infections in those areas.

Data from the literature and public databases have been gathered to estimate the model parameters meeting the requirements of the approach. Due to the lack of access to private data, the estimation of some parameters (e.g., the number of flights arriving every day in each country) was based on simplifying assumptions.

Figure 1 shows a simplified schematic representation of our model. We separate the world human population into three states: susceptible, infected, and recovered. The population was divided into this three groups for the model estimates to assimilate to COVID‐19 data (Duarte et al., 2020) (i.e., daily number of cases/infected and recovered, whereas susceptible is everything else) and, thus, simplify the parameter estimation, as well as the understating of results and risk communication to the general public. Note that the primary structure of our model follows the SIR modeling approach, however, there are secondary structures that are integrated into our model and make it unique (e.g., age range structure, dispersal among continents, DD effect), as will be detailed below.

FIGURE 1

Simplified schematic representation of Covid‐19 dynamics in the human population

Indeed, we subdivided each state into two age groups: the nonelderly (<60 years) and the elderly (≥60 years). The 60 years threshold was based on Bonanad et al. (2020), who analyzed more than half million of COVID‐19 patients from different countries. The authors highlighted the determinant effect of age on mortality rates with the relevant threshold on > 60 years, and concluded that > 60 years patients should be prioritized in the implementation of preventive measures. Although a model with more population categories would be more realistic, it would become more complex and less tractable (Pastorok et al., 2002), since the number of parameters would grow exponentially with the number of age classes. Hence, the parameterization would be more challenging and simulations would last longer.

When a susceptible individual gets infected, whether they are nonelderly or elderly, they stay in this state for some time and then becomes either recovered or dies (represented by the “Deaths Counter” box). Although, after recovery, a person can get infected again with a greater than zero probability (Lan et al., 2020), by the time this paper was written, there was almost no data on reinfection. Thus, we do not consider reinfection in this work, that is, we assumed that the probability of a recovered individual being infected again is 0. The transition between states is governed by random variables that follow probability density functions (PDFs) that vary over time, that is, the parameters of the PDFs are time dependent.

The world population is divided into six continents (subpopulations) with potential for dispersal. The dispersal rates from one continent to another also follow PDFs. We group countries into continents to keep the model and communication of risk more straightforward. A global model that represents dispersal among all countries could become intractable, resulting in challenging risk communication to authorities and the public.

To sum up, the structure of the proposed model can be tailored to incorporate many realistic and case‐specific features, such as: (i) the spatial structure of the infected population at the global level over time with the potential for dispersal among continents; (ii) the population structured by age range (nonelderly and elderly) with different probabilities of fatality and infection; (iii) control measures such as business restrictions and social isolation (reducing exposure), medical tools targeted at decreasing death rates, and travel restrictions between continents. Tables I and II, present descriptions of the variables and parameters respectively.

TABLE 1

Definition of the model variables

Variable	Symbol	Description
Number of susceptible nonelderly in continent i at time t	N2i(t)	Assessment endpoint described as a minimum, average, and maximum values, with a 95% confidence interval
Number of infected elderlies in continent i at time t	N5i(t)	Assessment endpoint described as a minimum, average, and maximum values, with a 95% confidence interval
Continent‐specific frequency of infection	μ^i	Number of expected new cases of infection generated by one infected person in each continent per week
The standard deviation of the frequency rate	σRi	Continent‐specific standard deviation from the frequency of infection
Exposure level for nonelderly (s=2)and elderly (s=5) for each SCN k.	Esk	Accounts for the reduction in the exposure due to a SCN k
Duration of social isolation for nonelderly (s=2)and elderly (s=5) for each SCN k.	Ssk	Accounts for the duration time of reduction in the exposure due to a SCN k
Ceiling for each continent i	K^i	Total initial population for continent i
Travel restriction for each SCN	T_r	Proportion of flights reduced as a measure to lower the spread of the infection ()
Fatality rate for (s=2) and (s=5)	d_i	Proportion of individuals that die from the infection each week

Note: [s = 1] Nonelderly susceptible, [s = 2] nonelderly infected, [s = 3] nonelderly recovered, [s = 4] elderly susceptible, [s = 5] elderly infected, and [s = 6] elderly recovered.

TABLE 2

Definition of the model parameters

Parameter	Symbol	Assumptions (data source or rationale)	Mean	Standard deviation
Age‐specific exposure	β_su	The probability of the virus being transmitted among the nonelderly is higher than among the elderly (Exposure assessment section and Table 3).	β_22=1 β_25=β_52=0.5 β_55=0.25
Time to recover	T_rec	Most individuals take two weeks to recover (37).	2
Permanence rate in s=1 and s=4	a_11;a_44	The susceptible population is much larger than the infected population, so there is a slight decrease in the susceptible population as more people get infected (user/author input).	0.999	0.001
Infection rate for s=2 and s=5	a22i(t);a55i(t)	Directly proportional to the frequency of infection and corrected by the age‐specific exposure, recovery rate and exposure level in each SCN.	E2k(β_22μ^i−a_32); E5k(β_55μ^i−a_65)	E5k(β_22σRi−σ_32); E5k(β_55σRi−σ_65)
Infection rate from s=5 to s=2; from s=2 to s=5	a25i(t);a52i(t)	Same as a22i(t);a55i(t).	E2k(β_25μ^i−a_32); E5k(β_52μ^i−a_65)	E2k(β_25σRi−σ_32); E5k(β_52σRi−σ_65)
The recovery rate from s=2 to s=3; from s=5 to s=6	a_32;a_65	Assessment of frequency section.	1T_rec(1−α_2); 1T_rec(1−α_5)
Permanence rate in s=3 and s=6	a_33;a_66	The probability of a recovered individual being reinfected is zero (Lan et al., 2020).	1.0
Fatality rate for s=2 and s=5	α_2;α_5	Description of SCNs section	f^k∗d_2; f^k∗d_5
Reduction in fatality‐rate due to medical tools	f^k	Reduction in the fatality rate for each SCN k (user/author input)	{f^SCN−1=f^SCN−3=1f^SCN−1=0.5f^SCN−2=0.25
Exposure level for nonelderly (s=2)and elderly (s=5) for each SCN k.	Esk	Accounts for the reduction in the exposure due to a SCN k (user/author input)	See exposure assessment section
Duration of social isolation for nonelderly (s=2)and elderly (s=5) for each SCN k.	Ssk	Accounts for the duration time of reduction in the exposure due to a SCN k (user/author input)	See exposure assessment section.
Dispersal rate of individuals among continents	m_ij	Parameterization of the model and Initial Conditions section	(M_ijN^i)×T_r
Threshold for infected population explosion	I_exp	Explosion threshold (25% of the world population infected)	1, 947, 531, 610

Note: [s = 1] Nonelderly susceptible, [s = 2] nonelderly infected, [s = 3] nonelderly recovered, [s = 4] elderly susceptible, [s = 5], elderly infected and, [s = 6] elderly recovered.

and

;

The model describes a population consisting of six subpopulations in the continents EU, NA, LA, AS, AF and OC, with potential for dispersal from one patch to another. The structure of the state of each population is [s = 1] nonelderly susceptible, [s = 2] nonelderly infected, [s = 3] nonelderly recovered, [s = 4] elderly susceptible, [s = 5] elderly infected, and [s = 6] elderly recovered, where s is the state. Dispersal is state‐specific, which means we can simulate different dispersal rates for the nonelderly and elderly.

Let denote the number of infected people in the continent i at time t, and , the total population in i at t. Then, the model forward projects the number of infected people, , (i.e., number of nonelderly infected; the number of elderly infected ) for all i, for 52 weeks (1 year) from March 17th, 2020 to March 16th, 2021.

In population ecology, the DD is the modification in the influence of any factor that affects population growth as the population density changes (Akçakaya et al., 1999; Burgman et al., 1993). In this paper, we did not aim at evaluating the growth in the world population. Therefore, population size is held constant (there are no inclusions of new susceptible individuals), and then the DD type is considered to be the “Ceiling.” The basis for accounting for DD was the total subpopulation of each continent (susceptible + infected + recovered). However, DD was modeled to affect only the number of infected. Under this type of DD, the infected population exponentially grows until the total subpopulation (susceptible + infected + recovered) reaches the ceiling (e.g., when the number of susceptible individuals is so low and recovered is so high that the number of infections stops growing). Then, it remains at that level until there is a decline in the total population (e.g., a random fluctuation or emigration) that takes it below the Ceiling (Akçakaya et al., 1999; Akçakaya & Root, 2013).

Note that we do not consider inclusions of new susceptible individuals (babies being born), but there are inclusions of the newly infected and, then, of recently recovered individuals, so DD matters because it limits the number of infected individuals; otherwise, this number would grow exponentially and infinitely. The Ceiling is continent‐specific, denoted by , and it is defined to be the total current subpopulation in that continent. Given that, the following algorithm represents one replication for stochastically simulating the epidemiological model. For each iteration, we repeat the next steps for all continents i: where is the transition rate from state u to state s and is the permanence rate in state u (u s ∈ {1,2,3,4,5,6}), and α2 and α5 are the mortality of infected nonelderly and elderly individuals, respectively. Note that some transition rates (i.e., a 11, a 22, ) are random variables that follow PDFs with parameters that are constant over time; therefore, a value is randomly sampled from the associated PDF for an iteration and kept constant for the entire 52 time‐steps of this iteration. Other transition rates (i.e., and ) are nonparametric stochastic processes because they are both random and dependent on the interaction among individuals; therefore, their PDFs change over time, and a value is sampled from the associated PDF at t. where is a matrix comprising the dispersal rates () of individuals from continent j to continent i.

To estimate the continent‐specific stage numbers:

To update estimates to account for DD: .

To update estimates of to account for the dispersal of individuals between continents by adding the number of entries and subtracting the number of exits for each subpopulation:

RAMAS® Metapop v.6.0 software (Akçakaya & Root, 2013) was adopted for running the simulations with 10, 000 replications. A review of 22 Monte Carlo simulation studies showed that the minimum recommended number of replications was always less than 10,000 and, generally, 7500–8000 replications are sufficient to produce stable results (Mundform et al., 2011). Although 8000 replications would be sufficient in our study, the only (quite insignificant) advantage is that simulations would last a few seconds less in comparison to 10,000. For example, the simulation of this model took approximately 2 min 11 s and 2 min 49 s with 8000 and 10,000 replications, respectively, in a decent computer (i.e., Intel® Core™ i7‐8550U and an 8 GB 1.80 GHz RAM memory). Yet, we follow the RAMAS® manual recommendation (Akçakaya & Root, 2013, p. 96): “use the maximum number of replications unless you are running a test simulation or making a demonstration”. Figure 2 shows the CI limits of the risk curves (± percent risk of explosion) as a function of the number of replications. Both axes are in logarithmic scale. Note that for 10,000 replications, the 95% CI limits is almost null (i.e., ±0.004), in other words, if our model estimates zero risk of explosion, the actual risk may not be higher than 0.4% (based on Kolmogorov–Smirnov test; see, Sokal & Rohlf, 1981, p. 721). We share all the model files in RAMAS format (Duarte et al., 2020).

FIGURE 2

Log‐log plot of 95% confidence interval limits of risk curves as a funcion of the number of replications (± percent risk of explosion)

MATERIALS AND DATA

Quantitative risk analysis approach

We used the structure of the model above presented to conduct a QRA for COVID‐19 by using the following steps (Duarte et al., 2019): (i) to characterize the problem; (ii) to describe the scenarios (SCNs); (iii) to assess exposure; (iv) to assess frequency; (v) to parameterize the model and initial conditions; and (vi) to quantify and categorize the risks. Duarte et al. (2014) and Siqueira et al. (2021) already applied this methodology to, respectively, run QRA for schistosomiasis disease and COVID‐19 in Brazil.

QRA is closely linked to risk communication (i.e., the effective transfer of technical information regarding possible risks to non‐technical audiences) (Teaf & Kuperberg, 2004). The way risk is conceptualized and described could be significant for how the authorities judge the magnitude of the risk, communicate it to the public, and conclude what to do (Aven & Bouder, 2020). It is almost useless to quantify risks if they cannot be well understood by the nontechnical audience, as policymakers are not necessarily experts on the language of probability. Thus, quantified risk can be transformed into classes that are much easier to interpret.

Therefore, risk categories have been used in all fields of QRA to make risk communication more straightforward (e.g., industrial QRA (CPR18E, 2005), ecological QRA (IUCN, 2018), microbial QRA for water safety management (Haas et al., 2014), microbial QRA of schistosomiasis (Duarte et al., 2014). However, to the best of the authors’ knowledge, there is not yet consensus in the literature on categorizing the quantified risk of the explosion of a pandemic at the population level. Thus, we here propose four groups, and the correct understanding of these is vital for a correct interpretation of the results:

CRITICAL RISK (CR): probability of explosion within 6 weeks > 50%.

HIGH RISK (HI): probability of explosion within 12 weeks > 20%.

CONSIDERABLE RISK (CO): probability of explosion within 52 weeks > 10%.

NEGLIGIBLE RISK (NE): probability of explosion within 52 weeks < 10%.

The method for assigning the above categories is as follows. Quantitative risk has three dimensions: probability, undesired consequence, and time frame (Duarte et al., 2019; IUCN, 2018); then, we established bounds for these three dimensions to form a risk category. In our case, the undesired consequence is the explosion of the disease (more than 25% infected people in the world at t, similarly to the Spanish Flu); this threshold is the same for all categories. Regarding the probability dimension, the bounds are the same as in the red list categories of the International Union for the Conservation of Nature (IUCN, 2018) (i.e., > 50%, > 20%, > 10%, and < 10% for CR, HI, CO, and NE respectively). Concerning the time dimension, it is also based on the IUCN categories (i.e., 10, 20, 100, and 100 years for CR, HI, CO, and NE, respectively), but adapted to the time horizon in which we make the forecast (i.e., 52 weeks). Thus, we have 6, 12, 52, and 52 weeks for CR, HI, CO, and NE respectively.

Note that the proposed categories do not consider the probability of massive deaths as is common for industrial QRA (e.g., (CPR18E, 2005). Our categories seek to indicate the risk of overloading the health system, which is associated with infections, high numbers of sick people, and substantial societal and economic costs. Conversely, risk categories based on deaths could neglect very infectious diseases with low death rates, even though the health system would be overloaded. Thus, we preferred to consolidate our undesired consequence in terms of infections, as these categories can serve as a basis for epidemiological modeling in future pandemics.

Characterizing the problem

The problem consists of quantitatively assessing the risks of SARS‐CoV‐2 to provide health managers worldwide with objective answers about the dynamics of the disease under several control strategies. Because of the arguments in the previous paragraph, we chose as assessment endpoint the number of infected people. Nevertheless, our model can also calculate the number of deaths. This evaluation is based on a probabilistic method that provides risk results as a PDF for infections and deaths over time, along with an average value and a confidence interval.

This assessment is intended to be conservative. Thus, whenever different sources provided distinct parameters’ estimates for the PDF that governs a transition rate, , we have considered the most conservative values. The outputs are as follows: (i) estimate of the infected subpopulation size over time for each continent and in the world for 52 weeks; (ii) projection of the accumulated number of deaths in the world over time for 52 weeks; (iii) risk curves of explosion; (iv) time to explosion; (v) risk categorization; and a (vi) comparison of these results for all scenarios defined in the next section.

Data regarding the number of infected people for each day, from January 1st, 2020 until March 23rd 2020, for each country, were gathered from the public database managed by the Johns Hopkins University Center for Systems Science and Engineering (JHU CSSE, 2020). The processed data (grouped by continents) are available in the Appendix, Table A1. The primary sources of information used for a general qualitative understanding of the SARS‐CoV‐2 dynamics were (P. Chen et al., 2020; Choi & Ki, 2020; Ferguson, Laydon et al., 2020; Ferguson, Walker et al., 2020; Khan & Atangana, 2020; Lin et al., 2020; Livingston et al., 2020a; Pang et al., 2020). The specific source of information used for estimating each parameter of the model is presented in Table 2.

TABLE A1

Record of infected per day in each continent of interest. Adapted from: JHU CSSE (2020)

Date	Europe	North America	South America	Africa	Asia	Oceania
1/22/20	0	1	0	0	554	0
1/23/20	0	1	0	0	652	0
1/24/20	2	2	0	0	937	0
1/25/20	3	2	0	0	1429	0
1/26/20	3	6	0	0	2105	4
1/27/20	4	6	0	0	2912	5
1/28/20	8	7	0	0	5558	5
1/29/20	10	7	0	0	6143	6
1/30/20	10	7	0	0	8208	9
1/31/20	16	11	0	0	9891	9
2/1/20	20	12	0	0	11993	12
2/2/20	22	12	0	0	16740	12
2/3/20	24	15	0	0	19829	12
2/4/20	25	15	0	0	23838	13
2/5/20	25	16	0	0	27580	13
2/6/20	25	16	0	0	30761	14
2/7/20	28	18	0	0	34268	15
2/8/20	33	18	0	0	36992	15
2/9/20	34	18	0	0	40017	15
2/10/20	39	18	0	0	42553	15
2/11/20	41	19	0	0	44590	15
2/12/20	42	19	0	0	44968	15
2/13/20	42	20	0	0	60114	15
2/14/20	42	20	0	1	66587	15
2/15/20	43	20	0	1	68664	15
2/16/20	43	20	0	1	70788	15
2/17/20	43	21	0	1	72722	15
2/18/20	43	21	0	1	74512	15
2/19/20	43	21	0	1	74936	15
2/20/20	43	21	0	1	75481	15
2/21/20	60	24	0	1	76083	19
2/22/20	102	24	0	1	77794	22
2/23/20	195	24	0	1	78030	22
2/24/20	273	61	0	1	78518	22
2/25/20	373	62	0	2	79257	22
2/26/20	527	68	1	2	80057	22
2/27/20	789	71	1	2	81148	23
2/28/20	1061	75	1	4	82219	24
2/29/20	1420	92	2	4	83718	26
3/1/20	2120	103	9	5	85316	28
3/2/20	2610	130	9	9	86693	31
3/3/20	3194	153	12	12	88559	40
3/4/20	4119	187	17	21	89794	55
3/5/20	5483	259	23	24	91071	58
3/6/20	7108	317	36	43	93120	64
3/7/20	9150	462	44	43	94858	68
3/8/20	11526	589	73	86	96071	81
3/9/20	13912	667	84	94	96939	96
3/10/20	16693	1045	105	106	98138	112
3/11/20	21184	1397	147	122	99903	133
3/12/20	21912	1792	182	138	101207	133
3/13/20	33073	2384	354	176	103065	205
3/14/20	40113	2951	439	254	104982	256
3/15/20	47100	3792	502	320	106913	305
3/16/20	55715	5100	731	410	108513	386
3/17/20	76870	7166	1045	528	110493	464
3/18/20	90528	8751	1162	652	112517	588
3/19/20	108928	14891	1652	841	114974	710
3/20/20	129446	20621	2268	1042	117245	832
3/21/20	150950	27530	3013	1250	119944	1125
3/22/20	169466	35733	4146	1511	123004	1383
3/23/20	169334	35798	4164	1568	123045	1383

Description of scenarios

It is pretty intricate to predict or assess all the potential events (e.g., meteorological and environmental conditions, numerous control strategies, various novel medical tools, changes in hygiene and cleaning habits, transportation restrictions in all modes, and events like virus mutation) that might occur in the future and influence SARS‐CoV‐2 transmission. Thus, our model does not aim to be precisely predictive, only descriptive.

In this context, we defined three scenarios (SCNs), and then we compared them to a benchmark (SCN‐0) so that we can evaluate changes in SARS‐CoV‐2 dynamics (and the reduced/added risk) caused by each scenario. The benchmark scenario (SCN‐0) is defined as follows:

Do‐nothing plan;

Business‐as‐usual;

No social isolation neither for the non‐elderly nor the elderly (100% exposure);

No travel restrictions;

No new medical tools to reduce death rates.

To assess the efficiency of integrated containment strategies without excessive information to represent the many possible combinations, we evaluated three scenarios, which we believed to be representative for the next year:

Moderate mitigation (SCN‐1):

○Social isolation of the nonelderly (assumed 18% exposure) and the elderly (considered 9% exposure) for the next 2 and 6 weeks, respectively. After that, exposures go back to 100%;

○75% reduction in intercontinental flights in the next 6 weeks, and afterward flights are back with a 25% reduction in the usual volume;

○Medical tools targeted at reducing fatality rates by 50%;

○Ceteris paribus SCN‐0.

Strong mitigation (SCN‐2):

○Social isolation of the nonelderly (assumed 18% exposure) and the elderly (considered 9% exposure) for the next 7 and 17 weeks, respectively. After that, exposures go back to 100%;

○75% reduction in intercontinental flights in the next 12 weeks, and after that, flights are back with a 25% reduction in the usual volume;

○Medical tools targeted at reducing fatality rates by 75%;

○Ceteris paribus SCN‐0.

Also, there has been much discussion about the so‐called vertical isolation, which affects only the elderly and groups at risk. The strategy is criticized by international health organizations (Schuchmann et al., 2020). Then, we analyzed one more scenario to evaluate the effectiveness of this strategy.

Vertical isolation (SCN‐3):

○Business as usual for all the nonelderly (100% exposure);

○The elderly were entirely isolated (0% exposure) during 1 year;

○100% travel restriction for the elderly;

○Ceteris paribus SCN‐0.

With SCN‐3, our aim is only to assess the strategy's impact on risk, not the associated psychological harm that it could cause due to 1‐year isolation. The rationale for the exposure levels in each scenario will be further explained in the next section.

Exposure assessment

Human exposure to SARS‐CoV‐2 mainly occurs when people leave home in their day‐to‐day life (Gatto et al., 2020). In SCN‐0, we assumed the nonelderly are out for business 46 h/week and the elderly for 23 h/week. Thus, we set this level of exposure for our benchmark as 100% for the nonelderly and 50% exposure for the elderly. Next, we make changes in the exposure level to represent social isolation in mitigation scenarios. Table 3 shows a summary of the exposure assessment.

TABLE 3

Summary of the exposure assessment for each scenario (SCN)

	Estimated time out of home for any activity (h/week)		Exposure level (EsK)		Duration of social isolation (weeks) (Ssk)
Scenarios	Nonelderly	Elderly	Nonelderly	Elderly	Nonelderly	Elderly
SCN‐0	46	23	100%	50%	0	0
SCN‐1	8	2	18%	9%	2	7
SCN‐2	8	2	18%	9%	7	17
SCN‐3	46	0	100%	0%	0	52

Frequency assessment

When exposed to infected individuals, a susceptible individual may get infected. The frequency of infection per week can be estimated by processing the data on the daily number of confirmed cases in each country provided by (JHU CSSE, 2020), resulting in Table A1 (Appendix). To that end, we grouped countries into continents. Then, we calculated the weekly frequency of infection in each continent for every day D (i.e., the number of infected people at divided by the number of infected people at t). From this sample of values, we calculated the mean and standard deviation in each continent and checked if there were outliers outside a 99.7% CI (). If there were outliers, we calculated μ and σ and checked for outliers again. We repeated this process for each continent until there were no more outliers. Table IV presents the frequency of infection per patch, which corresponds to a do‐nothing scenario (SCN‐0).

TABLE 4

Frequency of infection per week (mean and standard deviation) for each continent

	Frequency of infection
Subpopulation	Mean	SD
Europe	4.3781	2.7887
North America	3.5772	2.6315
South America	10.2265	4.9092
Asia	1.2812	0.3117
Africa	5.3783	3.1823
Oceania	2.2280	1.1916

After an individual becomes infected, they may either die or recover. For the fatality rate ( (i.e., the rate at which infected individuals may die per week), we used as reference a study that presented the fatality rates in the world per age group deciles, that is, , until 80 + years old (Yanping, 2020). Then, we grouped fatality rates into the age classes of interest (<60 and ≥60 years old) and, then, estimated age‐specific fatality rates of the infected nonelderly and elderly: 0.006216 and 0.067575 per week, respectively. These rates mean that, on average, 0.62 and 6.76% of the total number of infected nonelderly and elderly populations die per week.

The recovery rates (i.e., the transition rate per week from infected to recovered) can be estimated based on the incubation and transmission periods of those who develop symptoms. A study (Huang et al., 2020) suggested that transmission of SARS‐CoV‐2 also occurs during the incubation period. Thus, recovery time can be considered as the sum of incubation and transmission periods. According to Lauer et al. (2020), most individuals take 14 days to recover under conservative assumptions. Thus, after 2 weeks, it is improbable that an infected individual would still be in the transmission period. This assumption is also under, and well supported by, the recommendation of the U.S. Centers for Disease Control and Prevention for the period of active monitoring of infected people (14 days = 2 weeks) (WhiteHouse, 2020). Thus, we estimated the mean recovery rate as , for the nonelderly individuals (note that the same can be done for the elderly, by using a 55 and α5).

Parameterizing the model and initial conditions

Tables 1 and 2 summarize the variables, parameters, and initial conditions of the model. Some studies in the literature have already estimated parameters governing the dynamics of SARS‐CoV‐2, which we used in the proposed model: the fatality rate per age class (Yanping, 2020); the mean incubation and transmission period (Lauer et al., 2020); the time taken to recover (Lauer et al., 2020); and the proportion of the nonelderly and elderly infected (Population Pyramid, 2019). Other parameters were estimated from raw data (Appendix, Table A1): the weekly frequency of infection and deaths for each continent. On the other hand, due to the distinctive characteristic of this disease, there was still a lack of scientific information, and then two parameters of the model were estimated via conservative educated opinions of the authors, that is, the permanence rate in state two (nonelderly susceptible); and the permanence rate in state four (elderly susceptible) (see Table 2 for the rationale and assumptions regarding these parameters). The estimates for these two parameters may be improved when more information becomes available.

Some parameters were estimated using a mean value and others a mean and standard deviation (mean and σ columns in Table 2). To make the latter uncertain, we consider that they follow a normal distribution. One can make good use of a Gaussian approach in the vital rates of biological models because there is a reasonable reason for random values not to be too far away from average, that is, physical limitations prevent substantial deviations and natural forces from equilibrium that brings vital rates back to their average values (DeGroot & Schervish, 2018; Montgomery & Runger, 2010). For probabilistic simulation, RAMAS® converts the normal distribution parameters into the corresponding Lognormal counterpart, which avoids bias resulting from truncation because all parameters are greater than zero.

We straightforwardly model dispersal. Indeed, we did not have access to private data of all flight arrivals and departures per country. Public information only shows the total number of arrivals per country and continent in 2017 and 2018 (World Tourism Organization, 2016). From this dataset, we took the average number of arrivals per year. We estimated the average number of arrivals per week by dividing it by 52 (the number of weeks in a year) for each continent. Thus, we model dispersal rates per week from continent j to continent i as , where is the number of arrivals in the continent i from the continent j per week, is the total subpopulation of the continent j and is the travel restriction management measure that can vary from 0 to 1. As a result, we had dispersal rates in the order of 10−3, 10−4, and 10−5. We built a dispersal matrix for SCN‐0 (Appendix, Table A2), where T(t) = 1; is dependent on each scenario (see Description of scenarios section) and varies over time as can be seen in the Appendix, Figure A1.

TABLE A2

Dispersal matrix between continents. Each element in the dispersal matrix means the proportion of the population of continent j (column) that travels to continent i (line) per week

	AS	EU	SA	NA	OC	AF
AS		0.00026	0.000365	0.000399	0.000208	0.000256
EU	0.000355		0.000456	0.000996	0.000260	0.000256
SA	0.000025	0.00026		0.000399	0.000104	0.000034
NA	0.000209	0.00091	0.000456		0.000260	0.000342
OC	0.000008	0.00010	0.000046	0.000199		0.000017
AF	0.000021	0.00016	0.000091	0.000266	0.000156

FIGURE A1

Flight restriction over time, Tr(t), for each scenario: SCN‐0 (business as usual), SCN‐1 (moderate mitigation), SCN‐2 (strong mitigation) and SCN‐3 (vertical isolation)

Although the current proportion of infected individuals is very low (less than 0.1%), we estimate the initial number of susceptible individuals by subtracting the infected individuals for each continent in each age class from the total population. The proportion of age classes was estimated based on available data regarding the age pyramid in each continent (Population Pyramid, 2019). The initial number of fatalities and recovered people were assumed to be zero because the proportion of fatalities was still low, and data about individuals recovered was scarce and difficult to collect.

RESULTS

Here, we present the main risk results of each scenario and make a comparison between them. SCN‐0 works as a baseline for comparing and quantifying the risk reduction caused by mitigation strategies. The results for the population were given in the following structure: average, ± SD, maximum, and minimum. Then, we present these outcomes as boxplots (Figure 3) to show the global infected population at the final time‐step of the SCN‐0 simulation, which helps identify the most likely continents where SARS‐CoV‐2 might increase. Based on the results, in the absence of interventions, AF is expected to be the continent with the most significant number of infected people shortly, followed by EU, SA, NA, OC, and AS.

FIGURE 3

Boxplots for the number of infections in each continent (in millions) at the final time‐step (after 52 weeks) for a business‐as‐usual scenario (SCN‐0). It presents the percentage of infections from the total subpopulation in each continent

Figure 4 illustrates the efficiency of integrated strategies (social isolation + flight restrictions + medical tools) for disease control: (A) a projection of the world infected population; (B) a projection of the accumulated number of deaths in the world; and (C) time to explosion (i.e., the cumulative probability distribution for the time taken for the percentage of infected people in the world to exceed 25%) for each scenario.

FIGURE 4

(A) Number of infected individuals in the world over time; (B) Death toll in the world over time; (C) Time to explosion (CR = Critical Risk; HI = High Risk; CO = Considerable Risk; NE = Negligible Risk). SCN‐0 (business as usual); SCN‐1 (moderate mitigation); SCN‐2 (strong mitigation; SCN‐3 (vertical isolation)

In Figure 4(A and B), results for each scenario are presented as mean values. For example, for SCN‐2m, the expected infected population in the world was estimated to be around 223 million, and the expected cumulative number of deaths is expected to be about 1.5 million in 52 weeks. Based on the results in Figure 4(C), it is possible to categorize the risks associated with each scenario (see Table 5): high (HI) for SCN‐0, High (HI) for SCN‐3, Considerable (CO) for SCN‐1 and Negligible (NE) for SCN‐2. Note that SCN‐2 is the only scenario in which the world infected population is estimated to be below the order of billions and the risk category is NE. Figure 4(C) also shows that, compared to SCN‐0, strong mitigation (SCN‐2) dramatically increases the time to explosion. In contrast, for moderate mitigation (SCN‐1) and vertical isolation (SCN‐3), it is still within one year, although the time to explosion is increased.

TABLE 5

Summary of the outputs for each SCN

	World infected population	Total death toll	Risk of explosion	Time to explosion	Risk category
SCN‐0 (benchmark)	Fluctuates between 3 and 3.6 billion	Fluctuates between 84 and 88 million	100%	8.6 weeks	HI
SCN‐1 (moderate mitigation)	Fluctuates between 1.45 and 3.46 billion	Fluctuates between 19.5 and 46.5 million	47.30%	Tends to infinity	CO
SCN‐2 (strong mitigation)	Fluctuates between 223 and 910 million	Fluctuates between 1.5 and 6.3 million	7.28%	Tends to infinity	NE
SCN‐3 (vertical isolation plan)	Fluctuates between 2.47 and 4.26 billion	Fluctuates between 15.4 and 26.6 million	100%	11.5 weeks	HI
SCN‐0 in AS; SCN‐1 in the other continents	Fluctuates between 1.67 and 2.71 billion	Fluctuates between 22.5 and 29 million	54.52%	19.8 weeks	CO
SCN‐2 in AF; SCN‐1 in the other continents	Fluctuates between 1.34 and 2.23 billion	Fluctuates between 37.2 and 62.1 million	50.60%	23.4 weeks	CO
SCN‐2 in AF and EU; SCN‐1 in the other continents	Fluctuates between 735 million and 2.03 billion	Fluctuates between 20.5 and 56.6 million	27.88%	Tends to infinity	CO

Note: HI = High Risk; CO = Considerable Risk; NE = Negligible Risk.

Other significant results are as follows: SCN‐1 (moderate mitigation) and SCN‐2 (strong mitigation) cause the risk of explosion to be reduced by 52.7% and 92.72% respectively compared to SCN‐0; the vertical isolation plan alone (SCN‐3) does not significantly reduce this risk when compared to SCN‐0, and so it is not helpful to maintain the prevalence rate below 25%.

To suggest a scenario with minimal social isolation for the nonelderly (hence, expected minimal impact on economy) that would provide NE risk, we carried out a sensitivity analysis of gradual decreases in the duration of social isolation for the nonelderly in SCN‐2, (Figure 5). This result showed that the explosion risk falls into the CO category for 6 weeks or less, suggesting that 7 weeks is the minimum duration of social isolation for the non‐elderly.

FIGURE 5

Sensitivity analysis for the non‐elderly isolation time

To account for continent‐focused strategies, we simulated: business‐as‐usual in the most populated continent (AS) and moderate mitigation in the rest, strong mitigation in the continent predicted to be the most hard‐hit (AF) and moderate mitigation in the rest, and strong mitigation in the two continents predicted to be the most hard‐hit (AF and EU). Results showed significant reductions in the explosion risks, although all of them were still within the CO risk region (see Table 5).

Scenario suitability analysis

Our model does not aim to be predictive, but the objective is to describe the importance of mitigation strategies in reducing risks. Thus, results are strongly dependent on the assumptions made in scenarios before simulation (Description of scenarios section), that is, on scenario‐dependent parameters , and (Table 2). These are flexible author inputs that can be easily changed to simulate alternative model scenarios.

The model scenarios (Materials and data section) were built and parameterized, simulated and their results generated (Results section) using data and information until March 23, 2020, which were quite scarce at that time. Not only the infections and deaths history was very small, but also the information available in the literature to assist in the establishment of the hypotheses for scenarios. Even so, with only 12 weeks of data, the model was able to generate useful and fast results to assist public managers in planning the duration and intensity of mitigation strategies.

After a few weeks have passed, it was possible to analyze the suitability of scenarios to reality. To that end, we compared the actual values of infections and deaths against results of a scenario that we chose as a reasonable fit for reality throughout the first six weeks of the simulation period, i.e., business as usual (SCN‐0) for AS and moderate mitigation (SCN‐1) for all other continents. This choice was based on many news from March 17th to April 28th (Bird et al., 2020; Machado, 2020; The Inquirer, 2020) that said Asian countries have reopened their businesses. We use boxplots with a 99% confidence interval (Potter, 2006) to analyze the estimated uncertainty interval and the root mean square error (RMSE) (Barnston, 1992) as a metric to compare the actual values against the expected simulated ones. The comparison is shown in Fig. 6A and Fig. 6B.

FIGURE 6

Scenario suitability analysis for the first six weeks after generation of results. The solid line represents the actual values and the boxplot the estimates for a business as usual SCN‐0 in Asia and moderate mitigation SCN‐1 in all the other continents: (A) number of infections (in millions); (B) number of deaths (in millions)

In Fig. 6A, one can see that the actual number of infections is within the 99% confidence interval, showing a good suitability when considering the propagation of uncertainty in the results. However, if one disregards the uncertainty of results and looks only to the expected values (horizontal lines in the boxes) one has a in the same order of magnitude as the existing infections in the period, indicating a relatively large deviation. Also, one can see that the expected values were all underestimated when compared to reality, indicating that our assumptions for the level (and duration) of social isolation were too optimistic, that is, , and (note: although our model was thought to be conservative, in the early stages of the pandemics, when this article was written, there was lack of information to estimate how long social isolation would take; by this time, we thought that we were being conservative in the assumed exposure levels and timing).

In Figure 6(B), the actual values of deaths are higher than the third quartile of the predicted boxplot for all six time steps, thereby indicating that our assumption of medical tools that would reduce α2 and α5 by 50% may not correspond to reality, at least for this particular 6‐week period. This divergence can also be justified by , which is almost half the confirmed deaths at the end of the first 6 weeks.

Furthermore, we compared the actual values of infections in the first 6 weeks with the results for each scenario applied to all continents in a general manner. We observed that SCN‐1 (moderate mitigation) in all continents was the one closest to what has happened in the world in general (no different continent‐specific scenarios). Figure 7 shows, for each continent, the comparison of the number of infections at week 6 (April 28th, 2020), where the black dots are the actual values, and the boxplots on the left of each graph are the predicted results for SCN‐1 in all continents. Note that there is a good suitability for all continents.

FIGURE 7

Scenario suitability analysis at week 6 (April 28th, 2020) for each continent, assuming: SCN‐1 in all continents (left boxplot); and SCN‐0 only in Asia and SCN‐1 in all the other continents (right boxplot). The actual values of infections are represented as dots and predicted values as boxplots

Discussion

In this section, we first discuss the advantages and then the limitations of using this model.

Advantages

This model was developed and the results were generated in the early stages of the COVID‐19 pandemics. Data and references used to feed our model was that available from February 10, 2020 until March 23, 2020. Even in a situation of data scarcity, the model made estimates for March 16, 2021 in the order of 109 and 108 for infections, and 107 and 106 for deaths, depending on the simulated scenario. The actual numbers on this date were in the order of 108 for infections and 106 for deaths (WHO, 2021), that is, within our range of estimates. We considered this an advantage, as the results could actually communicate the magnitude of the problem to society and public managers, and how much this magnitude could be reduced by mitigation strategies. As new pandemics caused by unknown virus emerges, a similar context will likely happen (i.e., lack of data and information about a totally new infectious disease). Specially in such situations, this study could be a useful reference for epidemiological modeling and risk assessment. Specific emphasis has been placed on the model flexibility, so that the same model structure and step‐by‐step approach could be used not only for COVID‐19, but also in future pandemics of unknown infectious diseases.

Our model proved to have great potential to be genuinely informative for decision‐making and not just one more deterministic prediction for managers to follow without understanding all the uncertainty around the data. Although the data were still very imprecise, our model propagates uncertainty in the results and gives answers regarding the distribution of consequences associated with probabilities. A “single‐point estimate” for the discretized time to explosion (e.g., in weeks), T, was calculated for every Monte Carlo run. After many Monte Carlo replications (e.g., 10,000), we had a set of “single‐point estimates” for the time to explosion and the number of occurrences of a “single‐point estimate”. Thus, we could calculate and present the probability of occurrence of each “single‐point estimate” (e.g., P(T) = number of occurrences of T/10,000). Then, for each time t, it was possible to calculate the probabilities of all T lower than t, which results in the Cumulative Distribution Function (CDF) for the time to explosion, that is, ). In summary, means the probability that the explosion will occur at or before t, and it was plotted in a graph (Figure 4C).

The great advantage of such an approach over deterministic analyses is that results show not only what could happen but how likely each outcome is. For example, another study (Wells et al., 2020) used Monte Carlo simulations to provide a risk graph that shows the cumulative probability over time (days) of exporting at least a single infected case from mainland China via international travel. Although the objective and scope of that study was different from ours, it also shows how probabilistic models and Monte Carlo simulation can provide results that incorporate the indelible uncertainty in the dynamics of COVID‐19.

We used 6 weeks (March 17th to April 28th, 2020) to analyze the suitability of scenarios by comparing actual against estimated infections and measuring the RMSE, and the most suitable scenario was SCN‐0 in AS and SCN‐1 in all other countries (Figure 6). Also, for this scenario, the continent‐specific estimates versus the actual number of infections (Figure 7) showed good adequacy for all continents. Note that we presented results for several scenarios so that one of them will probably correspond to reality. Our model cannot make precise predictions of what will exactly happen in the future. Any model that tries to do that most likely will miss some information because decisions are taken every day and change the future. In this sense, a great advantage of our model is to provide not only single‐point estimates (expected values) for infections and deaths, but also a confidence interval as a measure of uncertainty. Moreover, it allows for fast simulation (i.e., less than 3 min per scenario in a decent computer), so it can be used as a tool to estimate the impact of decisions before they are taken so that authorities of the most representative countries in each continent may be warned of the risks of their decisions to world health.

Finally, we proposed criteria to categorize the quantified risks of a pandemic at the population level, which can be helpful not only for COVID‐19, but also as a reference for classifying risks of any pandemic in the future. To the best of the authors’ knowledge, there have been no suggestions in the literature so far on how to categorize such risks. The rationale behind the categorization was explained in the QRA approach section. Readers must ensure that they fully understand it to be confident that they can correctly interpret the results of the risk category.

Limitations

The following limitations were mostly due to the scarcity of both data and peer‐reviewed publications about COVID‐19 at the time this study was conducted, as will be bellow detailed. We remind readers that our model was developed and the results were generated in the early stages of the pandemic.

This epidemiological model has been developed to represent the whole world population. However, a major simplification was to subdivide the spatial structure of the world into continents only, without detailing countries, regions, states, and cities. In fact, COVID‐19 affects various countries across all continents at different rates (of infections and deaths). Also, within a country it affects differently (e.g., hot spots). A model that captures this spatial variability of variables and parameters within continents and countries would be more realistic, but also more complex and less tractable.

In the present application of our model, we considered that the probability of a recovered individual being infected again is 0. Yet, there is evidence that recovered individuals may become reinfected (Lan et al., 2020). However, it was still not sure whether such individuals have really become reinfected or were not infected before, but the test result was a false positive. However, our model is flexible in parameterization, and this probability can be easily changed to a value greater than zero when more information becomes available.

The scenario suitability analysis section showed that the expected death toll is considerably lower than in reality (Figure 6B). SCN‐1 assumes that new medical tools would be developed and reduce fatality rates by 50%. We acknowledge that this assumption was too optimistic and can be improved in future applications, as more information regarding new medicines/therapeutics is available.

We assumed that it is improbable that an infected individual would still be in the transmission period after 2 weeks of infection. Although new evidence indicates that the transmission period has a significant chance of being greater than 2 weeks (Livingston et al., 2020b), this was a valid assumption when the model was developed. For future applications of the model, this can be improved by making changes to the standard deviation of the time to recover parameter, , which this paper assumed to be zero (Table 2).

There was significant uncertainty and lack of clarity regarding how data about the number of confirmed cases in each country have been collected. This uncertainty can be seen in the high values of standard deviations for the frequency of infection per week, especially in the poorest continents (AF and SA). Thus, the estimates for the order of the hardest‐hit continents (Figure 3) should be treated with care. This ranking should be thought only of an initial guide for prioritizing resources among continents. Nevertheless, unlike other models, our approach could propagate uncertainty in the global results by using probabilistic language expressed in boxplots (Figure 3).

DD was modeled in such a way that the frequency of infection per week, , was assumed to be constant over time until the number of susceptible + infected + recovered individuals reached a ceiling. Then, it remains at that level until a decline in the population (e.g., a random fluctuation or emigration) takes it below the Ceiling. This approach was a conservative and straightforward way of limiting the growth in the number of infected people. It would be more realistic (although less conservative, in the sense that it will decrease risks) if the frequency of infection, , gradually reduces as the number of recovered individuals rises. This conservative approach can be improved in future studies by modeling DD as Scramble or Contest‐type (Akçakaya et al., 1999). Such an approach would benefit the model since the current ceiling DD type makes the infection progress faster than natural after loosening the social isolation. A model with contest–type of DD was proposed in an epidemiological model used to assess risks of COVID‐19 to the public health system in Brazil (Siqueira et al., 2021), and then this approach can be adapted to any other country or even the world.

This paper did not consider the risks of mitigation strategies to the global economy. We quantified and categorized microbial risks only as a measure of the probability of massive infections and deaths. However, social isolation and business shutdown have caused the income of many people to plunge. There are studies (Mortensen et al., 2016; Rehnberg & Fritzell, 2016; Wolfson et al., 1999) that show a relationship between population income and mortality (e.g., the higher a person's income, the better they can eat and take care of their health and the lower their mortality), and then it is possible to estimate this link and integrate it into our model. A proposal for future studies would be to feature in our model this characteristic.

Our model was built and simulated using a paid software called RAMAS® (Akçakaya & Root, 2013). Although we share all the model files (Duarte et al., 2020), it is only helpful for those with the RAMAS® license. We acknowledge this impairs the ease of reproducing the model. There is a strong movement toward reproducibility in science, especially near‐term ecological forecasting (Anderson et al., 2020; Dietze et al., 2018; White et al., 2019). Therefore, a second proposal for a future line of research is to build and simulate the model in an open scriptable software so that other researchers can easily reproduce it.

Other proposals for future studies include capturing the spatial variability of variables and parameters within continents and countries; conducting a sensitivity analysis to identify essential control measures; and undertaking a long‐term QRA of COVID‐19 to evaluate the effectiveness of alternative vaccine types and mass vaccination programs.

CONCLUSIONS

This study has presented a probabilistic epidemiological model, which was developed during the first few weeks after WHO declared the COVID‐19 pandemic, when data and literature about the disease was still scarce. By means of this model, we have quantified, assessed, categorized, and ranked the risks related to varying mitigation scenarios for the future and provided the results so that authorities can make informed decisions regarding the consequences of these risks. Global risks were assessed, and continent‐specific risks were helpful as an initial guide for prioritizing resources among continents.

The adequacy of scenarios was analyzed by comparing the results with actual values in the six weeks from March 17th to April 28th 2020. The outcomes showed that the simulated number of infections (for a business‐as‐usual scenario in Asia and a moderate mitigation scenario in the rest of the world) were underestimated. Our claim is that the main cause for this deviation is not in the model structure, but on the exposure assumptions for such scenarios, which were too optimistic.

Moreover, the estimated number of deaths has been lower than reality, mainly because we were too confident that medical tools would be developed and reduce fatality rates. Another reason for underestimated estimates was the DD type used in this model (i.e., Ceiling). Contest‐type of DD would be more adequate and provide a more realistic evolution for the number of cases.

The main advantage of using this model compared to others is that it is probabilistic by nature, so it provides results that incorporate the indelible uncertainty in a pandemic. The main limitation is the great uncertainty in the input data, which is also a hurdle for all other models in the literature that use public sources on the daily number of confirmed cases per country. Nevertheless, unlike other approaches, our model can inform managers about where there is uncertainty in the results to understand the risks arising from their decisions.

The main next challenges for this model are: to change Ceiling into Contest‐type of DD; to capture the variability of parameters and variables within continents and countries, while still maintaining an viable level of model tractability; to include the effect of varying COVID‐19 candidate vaccines and mass vaccination programs and, then, to conduct a longer‐term simulation (e.g., 5–10 years) to evaluate the effectiveness of such vaccines and vaccination programs; and to build and simulate the model in an open scriptable software. We consider that others can use our model, approach and risk categories in future pandemics. Eventually, this study could be an important reference in epidemiological modelling, especially under data scarcity and great uncertainty.