A country-specific COVID-19 model

Model motivation

The novel coronavirus (COVID-19) is a highly contagious and deadly airborne disease. It has claimed millions of lives worldwide since 2020, caused untold economic damages in every country, and is considered as the worst global public health crisis since the Spanish Flu. While each pandemic is inextricably characterized by the virus lineage and pathophysiology, COVID-19 is uniquely marked not only by the significant roles played by public and private mitigation strategies, accelerated developments in vaccination and therapeutics, but also the prevalent social and political beliefs in each country.

The challenge to cast all the causes and effects in a mathematically coherent framework is formidable. It is widely recognized that a credible pandemic model is one that captures the virological heuristics, that reflects the heterogeneity of immunity in a given population, and that describes and interprets empirical data with a demonstrable capability to forecast. Additional properties attractive to modelers encompass robustness and adaptability when applied to different geographies and over time, balanced between sufficient complexity and efficient model calibration.

The classical Susceptible-Infectious-Removed (SIR) compartmental framework pioneered by Kermack and McKendrick (1927) has been a foundational model of choice for communicable diseases. It provides a formulation with which the disease transmission mechanism is expressed in a mathematical form. It has been widely used for evaluation, prediction and infection prevention purposes. Extensive applications of SIR are found in the literature for chickenpox (Giraldo and Palacio 2008), measles (Bjørnstad, Finkenstädt, and Grenfell 2002) and MERS (Saleh et al. 2017), as discussed in Brauer (2005). Propelled by advanced computing and biomedical technology, numerous extensions of SIR made by creating additional compartments augmented with testable features have led to significant progress in the science of epidemiology, contributed to a richer and broader understanding across behavioral and demographical boundaries.

Specifically, at time t, the SIR model places the population of interest into one of three non-intersecting compartments, namely:

Under the construct of SIR, S(t) represents the population who are at risk of becoming infected, I(t) signifies the source and prevalence of the pathogen in circulation, and R(t) measures the population immune to the infection, no longer contagious, both healed or deceased are counted (Table 1).

The compartmental transmission mechanism described by SIR is diagrammed in Figure 1, with the model parameters defined in Table 2:

The SIR model is constructed as a system of 3-dimensional Ordinary Differential Equations (ODE), with the core discrete-time expression by Equation (1):

Figure 1:

Epidemiological flow of a pandemic under SIR.

The implicit assumptions made by SIR, as discussed by Brauer (2005), are threefold: (1) an average member of the population makes contact sufficient to transmit infection to β others per unit time; (2) the infected leave the infective class at a rate of γ per unit time, equivalent to an average of infection duration of 1 γ in units of time; (3) there is no entry into or departure from the population, except possibly through death from the disease. A crucial feature shared by the types of pathogens suitable of SIR modeling, as pointed out by Blackwood and Childs (2018), is that the recovered cohort is awarded with life-long immunity.

A non-linear numeric solution of this closed ODE system can be found once the parameters and the initial conditions are known. The world governed by SIR is the one in which the infected population I(t) increases precipitously in an outbreak, then reaches a plateau and eventually decreases at a much slower speed, as the immunity builds and the susceptible population S(t) is being depleted. By design, the SIR model is a self-regulating autonomous system as the drift terms ∆S(t) + ∆I(t) + ∆R(t) are netted to zero intrinsically at all times. While SIR is capable of tracking an endemic from the outset to reaching its finality, the inter-compartment transitions are however deterministic, with the evolution in I(t) features a single peak over the entirety of an epidemic lifecycle.

When applied to COVID-19 pandemic, these well-recognized SIR properties are found to be unrealistic and counterfactual, as several peaks and valleys have already been chronicled since 2020 worldwide. Attempts in fitting empirical data to the basic SIR model only lead to poor parameter estimates, offering no temporal or comparative value for surveillance purpose, and certainly an unlikely candidate for forecasting. The dynamism and realism required by the natural phenomena in COVID-19 pandemic give rise to the proposed stochastic model, allowing a probabilistic approach to describe the variability in the infectivity, resolution and fatality through time and across a variety of demographics.

On a consistent basis, available COVID-19 pandemic data is limited to the confirmed total (cumulative) as well as new cases of Infection and Death at a daily interval. Essentially we are left with all S(t), I(t), R(t) three quantities unknown and in need of a model estimation. It should be remarked that the fatality counts, for individuals deceased from infection, is known empirically although not present in Equation (1) as a standalone compartment.

In the cannon of epidemiological study, the Basic Reproduction Number represents the number of secondary infections that will be produced by a single initially infected case in a population of susceptible individuals. Put more precisely, it is a ratio between new infection cases and removed cases within a given time interval. In a SIR setting, the Basic Reproduction Number, represented by the symbol ρ hence after,^[1] is defined as:

In practice, the Basic Reproduction Number is often viewed as a precursor – signaling the degree to which the virus replicates in a population and urgent medical care is needed near-term. When ρ > 1, it signifies the virus infects individuals more rapidly than the number of individuals resolving their infection status, thus the total infection and death case counts are expected to both rise as a result. Conversely, when ρ < 1, the infected population is expected to decline in numbers along with the follow-on deaths. For individuals, medical professionals as well as public health officials, assuming ρ a constant and knowable only ex-post is of little value, even if calibrated from empirical data averaging over a lengthy period of time. Instead, allowing the Basic Reproduction Number to be instantaneous and vary over time as a state variable, is not only preferable but mandatory in order to be deployable in real-time.

The Mortality rate from COVID-19 is of the utmost importance both medically and economically, as discussed in Atkeson, Kopecky, and Zha (2020). Judged by the clinical information and evidence gathered across countries, the determinants for fatality are thought to range from patients’ age and underlying medical conditions, to hospital locality and the availabilities of ventilators and ICU beds, not to mention the great variety of COVID-19 variants. A time-varying dynamic Mortality Rate process, replacing the static or deterministic projection, is a critical feature in our proposal.

By modeling the infection reproduction and fatality rate as random processes, we should be able to achieve our dual-objectives: (1) synchronized monitoring of the disease transmissibility and its virology, critical to medical treatments and developments of therapeutics; (2) continuously evaluating the intervention effects from social and medical counter measures. This information is instantly distributable and consumable to the general public, useful for regulating self-protective behavior, and forming effective and risk-based social distancing measures with respect to remote working, home schooling, in-store dining, or safe transportation.

The remainder of the paper is organized as follows: In Section 2 we propose the eSIR model and discuss the epidemiological and model parameters in a structured setting. In Section 3 we describe the operationalized method for a general country-specific model calibration. Section 4 illustrates the procedures used in treating the raw empirical data and presents the main results accompanied by our observations and comments. In Section 5 we demonstrate the usefulness of the stochastic modeling choice by conducting historical simulation and forecasting. Section 6 concludes and discusses future research ideas.

The proposed e(nhanced)-SIRD model

Making use of a set of Stochastic Differential Equations (SDE) should lead to a probabilistic framework better suited for modeling COVID-19 and its intermittent transmissibility. Largely inspired by Perla et al. (2021) and Atkeson, Kopecky, and Zha (2020), we show that our proposed methodology admits a significant model improvement in reproducing empirical data and forecasting, without compromising the expected transparency and interpretability associated with the classical SIR model.

The 2-factor eSIRD model hereby proposed extends the deterministic structural SIR framework by modifying the Basic Reproduction Number ρ as a time-varying random variable ρ(t) which is referred to as the Reproduction Number, and by adding the fatality compartment D(t) formally which is jointly determined by the size of the infected population I(t) and another random factor – the case Mortality Rate m(t). In summary, the epidemiological schema ascribed by our eSIR model is presented in Figure 2.

Figure 2:

Epidemiological flow of a pandemic under eSIR.

Where all the state variables and associated model parameters are defined collectively in Table 3 for completeness:

It should be noted that, with the introduction of D(t), the new definition of compartment R(t) is narrowed to contain the recovered individuals only, modified from Table 2 whereas the removed compartment hosts both the recovered and the deceased. It should also be highlighted that this group does not include people with vaccine-induced immunity. Likewise, I(t) group contains only the people who are contagious and actively shedding virus, but not caused by vaccinations. The resolution rate γ, however, remains unchanged from Table 2, a constant parameter representing the outflow rate from the infected I(t) compartment irrespective of the receiving compartment being recovered R(t) or deceased D(t).

The flowchart of Figure 2 prescribes and delineates the stages of infection in a given population: individuals in the susceptible compartment, S(t), are first infected and moved into the infected individuals (tested positive) compartment I(t), at a rate of ρ t γ . Subsequently, the infected individuals resolve their status, at a rate of γ, either by death from COVID infection with a net rate of m(t)γ and moved into the fatality compartment D(t), or by recovering with a healing rate of 1 − m t γ and moved into the recovered compartment R(t), in each time step.

Ostensibly the population in the three compartments of interest, namely S(t), R(t) and D(t), experiences only single directional (in or out) flows, while the dynamic of the Infected compartment, I(t), is in-fact multi-directional, increasing from new infections out from S(t) and reducing to R(t) or D(t) concurrently. Same as with SIR, since the existence of COVID-19 is relatively short, the assumption of a closed demography is therefore justified, i.e. the natural births and deaths unrelated to the pandemic are not considered.

Finally, with variables and parameters defined in Table 3, we proceed to write the proposed eSIR model formally as 2 sets of inter-connected differential systems as Equations (3) and (4):

where ρ(t) denotes the instantaneous Reproduction Number and m(t) the instantaneous case Mortality Rate. Furthermore, they are each assumed to evolve following a well-known CIR stochastic process per (Cox, Ingersoll, and Ross 1985):

The CIR parameters required by Equation (4) are:

CIR model was introduced first in a seminal paper CIR (1985) and has been prominently applied in the field of financial mathematics, in particular for modeling interest rates. Stationary, positivity, local volatility and mean-reverting are the most distinguishing features when a CIR process is characterized (Table 4). This can be seen in our proposal where the diffusion coefficients are specified as σ ρ ρ t , σ m m t corresponding to strictly positive valuations of { ρ t , m t , as the infection reproduction and death rate shall never be negative in epidemiological terms. The long term mean, i.e. the latent asymptotic value ρ ̄ = lim t → ∞ ρ ( t ) or m ̄ = lim t → ∞ m ( t ) , has a well-understood role in a mean-reversion process and gives a well-defined asymptotic distribution. This is once again consistent with the repeated outbreaks already observed for COVID-19, and the biostatistical verity that is, with sufficient time, all pandemics do subside and eventually extinct.

The sources of risk or randomness that generate this 2-factor stochastic system in Equation (4) are the Brownian motions {W_ρ(t), W_m(t)} drawn from independent^[2] normal distributions N(0,dt) with zero mean and a dt-sized variance. In our approach, these native innovations are first transmitted through the state variables { ρ t , m t and then propagated to all four SIRD compartments governed by Equation (3), with the case counts of new or cumulative infections and deaths as the ultimate observables.

It is worth noting that reinfections from COVID-19, while not common, do occur. Repeat cases have been sporadically observed and reported, as the gains in herd immunity from natural infection or inoculation are partially offset by the countervailing factors such as low vaccine rates and viral mutations in a continuum. Without a clear clinical definition in reinfection, without a universal contact tracing and persistent screening test targeting each variant, an epidemiological model is unable to accurately measure and discern these effects judiciously. Absent of reliable and specific data and low prevalence, viewed from a global perspective, we disallow COVID-19 reinfections for modeling purpose to achieve our research objectives. As can be seen in the following sections, this simplification is also mathematically necessary to solve a highly nonlinear multivariate stochastic system empirically.

By construction, the epidemiological identity is preserved by our eSIRD model. With the rapid progression of COVID-19, it suffices to assume an individual must be in one of four said compartments at any time t. With this identity, namely N = S(t) + I(t) + R(t) + D(t), the population infection status follows a closed system with a net zero migration rate:

It is important to note that our model assumes homogeneity in a given country. This is in contrast to Britton (2020) and Lipton and de Prado (2020), whereas all intra-country heterogeneities such as locality, age or pre-existing condition cohorts, social activities are accounted for. Furthermore, our country-specific homogenous model does not attempt to make any causal inferences or distinguish the evolving COVID-19 virulent nature or the heterogeneous immunity. Although the waves of spread have been incrementally slowed down by varied containment policies, such as travel restriction, face-masking, social distancing, a precise measure of their effectiveness is hampered by inconsistencies both in policies and the actual adherences. It is widely acknowledged that behavior modification is a function of the infection prevalence, thus its standalone impacts are difficult to identify, or found to be statistically insignificant as argued by Atkeson (2021) and Last (2021).

However, the realized transmissibility, as a confluence of these multifaceted and interlinked longitudinal and cross-sectional factors, is ultimately captured effectively by the Reproduction Number ρ(t) following a reduced-form approach in our proposal. The stochastic characteristics prescribed by our model afford a more accurate measurement of pandemic progression transparent to all stakeholders – the general population and the state actors. It allows closer inspections and real-time reading of the interactions between the pathology of viral infection and the uneven human responses notwithstanding the widely available vaccination programs since late 2020. Under this parsimonious setting, we expect that the Reproduction Number vary substantially through the successions of infections, irrespective of endogenous and exogenous causes.

Historical calibration

Modeling epidemics has proved to be a challenging task for researchers. While it is widely known that COVID-19 case data is incomplete and highly imperfect, due largely to the lack of nation-wide rigorous testing and virology monitoring on a continuous basis, the credibility of a historically calibrated model is still primarily judged by the plausibility of epidemiological parameters. This subtlety is consequential, a topic to which we will return in our subsequent discussions.

While COVID-19 is considered to be moderately mutable, within the coronavirus family such as SARS and MERS, there are several invariant properties to which an empirical calibration should adhere. Constrained by the well-defined biomedical definitions of the infection resolution rate γ, and the asymptotic values ( ρ ̄ and m ̄ ) of the CIR processes, valuation with realism and reasonableness should be an overriding priority in determining the optimal model choice and the optimal calibration method.

According to studies from the World Health Organization (WHO) and Johns Hopkins University (JHU), it is a medical fact that the infection resolves within 1–2 weeks for mild symptoms, and up to 6 weeks for the seriously ill. Exploiting the biological information should lead to a superior interpretability and yield a substantial simplification needed for model calibration. We accordingly declare γ as a constant biomedical value. In other words, the COVID-19 infection resolution rate γ is assumed holistically given with a value ranging between 7 and 14 days, independent from model specifications or numerical calibrations.

With γ eliminated from the estimation, the remainder of required eSIR model parameters listed in Table 4 are gathered and collectively referred to as θ:

This global infectious statistics is centrally tracked by government information portals, led by the WHO^[3] and the Centers for Disease Control and Prevention (CDC^[4]), or by privately-owned data aggregators and research institutions such as JHU.^[5] The most followed data are the reported new cases of Infection ∆ C (t), and the daily reported new cases of Death Δ D t : = d D ( t ) .

Following the framework of Equation (3), and using ‘hat’ henceforward to symbolize the estimated latent variables given θ via historical calibration (discussed by the following Section), we can express mathematically the equivalence between the observables and the latent model variables as:

As stated earlier, reinfection disallowance implies that the rate of change in the Susceptible compartment S(t) is strictly negative, and that of the Deceased compartment D(t) is strictly positive. The min and max functions in Equation (7) can therefore be safely removed. In addition, d S ̂ t | θ is replaced by d S t in Equation (7) with the associated direct equivalence, considered as given and independent from parameter assessment.

We further assume that, for COVID-19 pandemic, the initial conditions for all compartments is uniformly disease free across all nations when t = 0. The terminal state and the epidemic equilibrium however, are yet to be determined.

Results and discussions

Since the declaration from WHO identifying COVID-19 as a public health emergency of international concern in January 2020, numerous governmental and healthcare organizations have been mobilized to track the pandemic status. For our country-specific analysis, the raw COVID-19 data is sourced from the Humanitarian Data Exchange (HDX),^[8] an open platform for sharing data across crises and organizations. It records the daily new and cumulative numbers of confirmed infection and death in a machine-readable format, for each country across the global. This dataset has been cited by a number of studies of COVID-19 pandemic, Last (2021) and Abbasimehr, Paki, and Bahrini (2021) are counted among the notables.

Data preparation

Given the difficulties in screening and uneven testing capabilities, intermittently changing public health policies and heterogeneous beliefs in every population, the COVID-19 data quality is significantly tainted with delays and extensive inaccuracies. In the presence of widespread noise, we choose to preprocess the raw data and use the smoothened time series instead as model inputs for calibration. Utilizing a Gaussian kernel smoothing algorithm, we can show in Figure 3 a visual comparison of raw data and their smoothened counterparts to illustrate the effects and the operational necessity:

Figure 3:

Left: COVID-19 cases | Right: COVID-19 deaths. Empirical data (original label) vs. Smoothened data (labeled with prefix ‘s’) of USA.

We then derive the solution Reproduction Number ρ ̂ t | γ and Mortality Rate m ̂ t | γ governed by Equation (10) and graphed in Figure 4:

Figure 4:

(Purple) Reproduction number ρ ̂ t | γ and (Black) mortality rate m ̂ t | γ – USA. Left panel assumes 7 days infection resolution | Right panel assumes 14 days infection resolution.

The estimated Reproduction Number ρ ̂ ( t | γ ) exceeded 10 during the initial acute onset of COVID-19. Shortly after Q2 of 2021, we observed the Reproduction Number dipping below 1 for the first time, resulting from the success of effective vaccine intervention and the efforts and resources invested by the Federal and State-sponsored inoculation programs to the entire US population. This is encouraging as having the Reproduction Number sub 1 is the first indicator of a subsiding infection trend in the population. Since the summer of 2021, however, the Reproduction Number has been fluctuating and moving upward slightly from being over 1 to above 2. There have been several subsequent spikes, with the timing coinciding with the emergence of new variants,^[9] holiday seasons and the varied prevalence of unvaccinated people in certain communities. From Alpha, Beta, to more recent Delta and Omicron variants, the unpredictability of COVID-19 mutating path is largely responsible for the recent uptick in transmissibility thus the increasing Reproduction Number.

In a global struggle against COVID-19, the death cases are systematically tracked and the Mortality Rate is studiously monitored to assess the ultimate vulnerability in a population. Unfortunately, the Mortality Rate is also a hard quantity to evaluate with precision. The initial high fatality numbers observed in Figure 4 reflect the rapid succession of COVID-19 outbreaks around the world, which have since moderated substantially. Most epidemiologists and healthcare providers believe that the underestimation of Total Confirmed Infection Cases is widespread, persistent and significant. This is plausible since very limited testing and contact tracing are conducted, developed and emerging market economies notwithstanding. As a result, the empirically estimated Mortality Rate as a ratio of T o t a l C o n fi r m e d D e a t h s T o t a l C o n fi r m e d I n f e c t i o n C a s e s would likely have been overestimated. OurWorldinData^[10] has reported the Mortality Rate to highly vary across time and geography, from 14% in Brazil, 6% in US in April-May 2020 to generally below 3% for many countries since 2021. By a different estimation, by taking in a smaller and self-contained sample such as the Diamond Princess cruise-ship outbreak, the Mortality Rate was found to be at around 1–1.5%. We believe our systemically generated m(t) is sufficiently validated by the referenced empirical results.

Guided by the low sensitivity to either 7 day or 14 day infection resolution time range, as demonstrated in Figure 4, we proceed with the subsequent analysis using 14 days fixed period as the COVID-19 resolution parameter γ for all countries.^[11]

The onset effects

To evaluate the performance of the regression based calibration, we present Figure 5 in which the residuals of Equation (14), namely ε t ρ and ε t m , are plotted against the timeline and Figure 6 their density distributions. While the onset effects are not statistically insignificant upon inspection, it can also be said that the residuals are concentrated in a small fraction of the duration of the still on-going COVID-19 pandemic.

Figure 5:

Time series of residuals (Left ε t ρ |Right ε t m ) of USA.

Figure 6:

Density functions of residuals (Left ε t ρ |Right ε t m ) of USA.

Early March through July of 2020 was the period during which the world first discovered COVID-19. These single peaks found on Figure 5 correspond to the devastating harm inflicted on an unsuspecting population, marked by mass hospitalizations and shocking number of fatalities. It was only after the acute stage that some stationarity began to incrementally emerge. Helped by voluntarily learned behavior modifications or forced by government mandates, and rapid developments in pharmaceutical intervention measures, the initial unimpeded outbreaks have evolved to become cascading infection surges. We believe these counter-interactive effects have created the necessary condition under which the Brownian motions became admissible. The inherent uncertainties in a protracted pandemic are captured by random walks featured by a stable long-run convergence once a global-scaled equilibrium is reached. However, as informed by the residual analysis, the onset effects of COVID-19 may have negatively impacted the calibration and model efficacy.

As an attempt to bring in the onset effects to the generation of a time-varying Reproduction Number for COVD-19, authors from Linka et al. (2020) propose a hyperbolic tangent function, in a form of ρ t = R 0 − 1 2 1 + tanh t − t * T R 0 − ρ t with t* representing the adaptation time or the Onset period and R₀ the assumed initial value of the Reproduction Number.

Separately, we have contemplated a time-varying drift term to each of ρ(t) and m(t) processes. Instead of the probabilistic expressions as in Equation (4), let us suppose two deterministic ODE processes as in Equation (15):

(15) d ρ ̃ t = φ ρ t d t = a ρ + b ρ t e λ ρ t d t   d m ̃ t = φ m t d t = a m + b m t e λ m t d t  

As an example, Figure 7 shows the aforementioned hypothetical onset effects after 2 sets of {a, b, λ} parameters calibrated based on ρ ̂ t | γ and m ̂ t | γ as inputs. Needless to say, the non-linear functional forms of φ_ρ(t) and φ_m(t) can be more general, and the differential order is not necessarily restricted to 1.

Figure 7:

Left ρ ̂ t | γ (black) and ρ ̃ t | a ρ , b ρ , λ ρ (red) | Right m ̂ ( t | γ ) (black) and m ̃ t | a m , b m , λ m (red) – USA.

In a fully integrated model with the onset and saturated random dynamics both captured, conceptually, the targeted random processes of state variable ρ(t) and m(t) should then be thought off as a pair of SDE processes of CIR overlaid with φ ρ t and φ m t drifts. The stochastic Equation (4) conceivably could be revised to:

(4*) d ρ t = φ ρ t + κ ρ ρ ̄ − ρ t d t + σ ρ ρ t d W ρ t   d m t = φ m t + κ m m ̄ − m t d t + σ m m t d W m t  

Compared to our proposed model as in Equation (4), adding onset effects has several nontrivial and unintended consequences. With the expansion of parameters, undoubtedly the calibration of Equation (4*) would require a robust numerical optimization algorithm of considerable power, within Method (2), (3) or (4) most likely with reasons discussed in Section 3.2. Short of a rapid convergence required by all countries, we find that the canonical trade-off between a superior model specification and an efficient calibratability is a delicate one, a subject certainly worthy of further study.

We therefore conclude that our eSIRD model, operationalized in its original simplified form, is the optimal proposal, fitting for our intended purpose at the current stage of COVID-19 pandemic, useful to produce reliable and actionable insights.

Simulation

This robust and efficient regression-based calibration is capable of generating realistic scenarios reflecting the real world COVID-19 empirical data. With the parameter set θ calibrated, we are able to simulate CIR^[15] processes of { ρ ̂ t | γ , m ̂ t | γ }, then track the resulting case counts of model simulated { S ̂ t | γ , θ , , R ̂ ( t | γ , θ ) ) }. A single path simulation of the modeled evolutions is graphed in Figure 8 for illustration purpose.

Figure 8:

Simulated single path {S(t), I(t), R(t)} solution – USA.

The simulated SIR state variables are then transformed into the observables, i.e. the cumulative cases of infections C ̂ (t) and deaths D ̂ ( t ) . Figure 9 compares the empirical observable with the sampled single path simulation results in the US. Reported 64 million cumulative infection cases (in purple), and over 800,000 deaths (in black) as of January 2022 are in close proximity to the simulated model results.

Figure 9:

Empirical { C (t), D ( t ) } vs. simulated { C ̂ (t), D ̂ ( t ) } – USA.

One of the benefits of the SDE model is to allow looking ahead and gaining perspectives well into the future through Monte Carlo simulation. As depicted in Figure 10, our eSIRD model can be used to harness future shocks within the framework therefore forecast both infection and fatality.

Figure 10:

Empirical { C (t), D ( t ) } (solid lines) vs. experimental forecasting (dotted lines).

Conclusion

In this paper, we apply an enhanced version of the stochastic epidemic model to the current ongoing COVID-19 pandemic. Our model retains its original compartmentality in a SIR structural form, fitted by two latent random variables of the Reproduction Numbers and the Mortality rate to be parameterized in a reduced-form fashion. Illuminated by a medical fact that the resolution time for infected individuals is nominally invariant, we are able to estimate the Reproduction Numbers and Mortality Rate closed-form and in real-time. This simple, effective and interpretable approach is directly and immediately useful for the general public and policy makers who are responsible for health protocols, as well as economic decisions.

Uniquely, our model is capable of capturing the evolution of COVID-19 pathogen and human behavioral and medical responses in each country, corresponding to changing government policy, risk tolerance and the availability of vaccination and medical treatments. The same methodology can be easily applied to regions or cities as a tool for a more granular public health risk management.

The calibration results and discussions presented in the paper are based on in-sample statistical regression. Partitioning the time series into training and forecasting sets could shed additional light into the choice of model specification and method for calibration. Compare to the strong performance in modeling the transmissibility, the stochastic process of Mortality Rate may require improvements. Further work, by incorporating the onset effects while preserving the mean-reverting features observed in the real-world, and effectively and efficiently calibrating a combined process could further improve the model predictability and tractability.