Modeling Effects of Spatial Heterogeneities and Layered Exposure Interventions on the Spread of COVID-19 across New Jersey

Abstract

COVID-19 created an unprecedented global public health crisis during 2020–2021. The severity of the fast-spreading infection, combined with uncertainties regarding the physical and biological processes affecting transmission of SARS-CoV-2, posed enormous challenges to healthcare systems. Pandemic dynamics exhibited complex spatial heterogeneities across multiple scales, as local demographic, socioeconomic, behavioral and environmental factors were modulating population exposures and susceptibilities. Before effective pharmacological interventions became available, controlling exposures to SARS-CoV-2 was the only public health option for mitigating the disease; therefore, models quantifying the impacts of heterogeneities and alternative exposure interventions on COVID-19 outcomes became essential tools informing policy development. This study used a stochastic SEIR framework, modeling each of the 21 New Jersey counties, to capture important heterogeneities of COVID-19 outcomes across the State. The models were calibrated using confirmed daily deaths and SQMC optimization and subsequently applied in predictive and exploratory modes. The predictions achieved good agreement between modeled and reported death data; counterfactual analysis was performed to assess the effectiveness of layered interventions on reducing exposures to SARS-CoV-2 and thereby fatality of COVID-19. The modeling analysis of the reduction in exposures to SARS-CoV-2 achieved through concurrent social distancing and face-mask wearing estimated that 357 [IQR (290, 429)] deaths per 100,000 people were averted.

1. Introduction

1.1. Background: Characterizing Pathways of COVID-19 Transmission

COVID-19 presents the greatest public health challenge humanity has faced in over a century. In 2020 it became the third leading cause of deaths in the US behind heart disease and cancer [1,2] and in the Fall of 2021 still remains a major threat to the nation’s and the world’s health and well-being. In the initial phases of the pandemic, before effective preventive and therapeutic pharmacological interventions became available, controlling population exposures to the novel coronavirus SARS-CoV-2 was the only public health option available for mitigating the disease. However, developing rational and effective exposure intervention strategies, especially during the first half of 2020, was substantially hindered by complexities and uncertainties associated with the physical and biological processes affecting the spread of the disease. On one hand, it was recognized quite early in the pandemic that SARS-CoV-2 is transmitted not only by symptomatic but also by presymptomatic and asymptomatic individuals [3,4]. On the other hand, in the initial phase of the pandemic there were significant uncertainties and disagreements regarding the relative importance of the various routes and pathways involved in viral transmission [5]. Droplet spray was initially thought to be the primary means of SARS-CoV-2 transmission [6]. Viral respiratory droplets are produced when an infected person sneezes, coughs, or talks; these droplets can be directly deposited on the noses, mouths, and eyes of people nearby (usually at distances < 2 m) and can cause an infection: this is typically referred to as “direct contact of transmission”. Direct/indirect contact provides another potential pathway, where people are exposed by touching surfaces carrying the virus, such as an infected person’s hand and other fomites, and then touching their own eyes, noses, and mouths.

The importance of aerosol inhalation in COVID-19 transmission was not fully recognized initially. In April 2020, Morawska and Cao [7] pointed out the following: “Hand washing and maintaining social distance are the main measures recommended by the World Health Organization (WHO) to avoid contracting COVID-19. Unfortunately, these measures do not prevent infection by inhalation of small droplets exhaled by an infected person that can travel distances of meters or tens of meters in the air and carry their viral content.” Almost three months later, in July 2020, a follow-up article by Morawska and Milton [8], was co-signed by 239 scientists and started with the statement: “We appeal to the medical community and to the relevant national and international bodies to recognize the potential for airborne spread of coronavirus disease 2019 (COVID-19). There is significant potential for inhalation exposure to viruses in microscopic respiratory droplets (microdroplets) at short to medium distances (up to several meters, or room scale), and we are advocating for the use of preventive measures to mitigate this route of airborne transmission.” This commentary was widely publicized and eventually agencies, such as the World Health Organization, adopted exposure mitigation measures that incorporated face-mask wearing. Aerosols (e.g., droplets or particles with effective diameters less than 5–10 µm) can contain infectious virions and remain suspended in air for hours [9,10,11], resulting in them being transported and inhaled by individuals at distances substantially beyond 1–2 m from the infected persons [12]. Aerosol-based transmission of SARS-CoV-2 has been explored in several laboratory studies [13,14,15], while extensive work is ongoing to further elucidate the role of this transmission mode on the spread of the COVID-19 pandemic [16,17,18]. A wide spectrum of studies reveal multiple physical and biological process complexities associated with the aerosol transmission pathway, ranging from the role of indoor ventilation patterns [19] to different viral loads in aerosols for different variants of SARS-CoV-2 [20,21].

It is now widely accepted [5,22,23] that all of the above pathways of transmission must be controlled in order to reduce exposures to SARS-CoV-2 and minimize the spread of COVID-19. To achieve this, a hierarchy of “layered controls” aiming to contain the emission of virus-containing droplets and aerosols, curtail inter-individual contact, and frequent cleaning/disinfection were recommended to the public as exposure-reducing interventions [5]. Wearing a face mask provides “source control” that reduces both emissions and intake of virus-laden droplets and aerosols; for example, N95 masks are personal protective equipment (PPE) for both the wearer and nearby individuals, while both disposable surgical masks and reusable cloth masks have been widely used around the world as “community protective equipment” (CPE). Physical distancing and spacing mitigate the pathways of direct contact, prevent large expired droplets from the infected individuals reaching other persons, and also reduce the risk of exposures to viral aerosols. It should be noted that people in proximity to the infected individuals are expected to be exposed to higher concentrations of virus-laden aerosols than those that are more distant. Ventilation and filtration are indoor environmental control strategies that can reduce the risk of aerosol inhalation by diluting/removing the virus-containing droplets and aerosols in the air. Cleaning and disinfection practices eliminate or reduce the viral load on fomites, mitigating direct/indirect contact, and also reduce the risk of exposure to infectious aerosols that can potentially be resuspended from fomites into the air.

Governing authorities of countries, regional entities, and cities instituted combinations of public health interventions to reduce general population exposures to SARS-CoV-2, including closures of schools and non-essential businesses, social distancing, wearing face masks, travel restrictions, testing, tracing and isolation of exposed individuals, lockdown of specific areas, etc. [24,25,26]. The effectiveness of these policies has varied substantially across different regions and has been closely related to how communities have complied with them. In reality, both person-to-person infectious disease transmission and the resulting severity of viral infection are complex processes that are driven by multiple biological, demographic/behavioral, and environmental factors. Horton [27] correctly stated that COVID-19 is not a pandemic but a “syndemic” in order to emphasize the importance of numerous biological, social and environmental determinants of health (e.g., sex, age, obesity, medical histories, occupation, prior exposures to chemical, physical, biological and psychosocial stressors, etc.). For instance, older individuals [28,29,30], particularly those residing in nursing homes [31], and males with preexisting conditions (e.g., metabolic disorders, such as diabetes; chronic lung conditions, such as moderate to severe asthma; cardiac diseases with complications, etc.) [32] have been far more vulnerable to COVID-19 than healthy younger individuals. Workers employed in essential services (e.g., medical facilities, food and medical supply chains, transportation infrastructure, government operations, security, etc.) face significantly higher risks of exposure to the virus than those who can work from home [33]. This has resulted to a significantly disproportionate impact of the disease on communities of color [34,35].

1.2. Modeling COVID-19: Computational Approaches

Computational models provide essential tools for framing our understanding of the complexity of disease dynamics, including underlying exposure and transmission processes, and for developing quantitative assessments and evaluations of alternative strategies needed to manage the pandemic [36]. Mechanistic modeling of epidemic dynamics originated with the publications of Ross [37] in 1910 and of Kermack and McKendrick in 1927 [38]; their work established a discipline that, as it evolved over time, has utilized a wide spectrum of methods, including Ordinary Differential Equations (ODEs), Difference Equations (DEs), Partial Differential Equations (PDEs), Integrodifferential Equations (IDEs), Cellular Automata (CA), Agent-Based Models (ABMs), Network Models, etc. (see, e.g., [39,40,41,42]). Many of these methods have been applied to studies of COVID-19; CDC [43] has been reporting predictions of COVID-19 spread and outcomes from an ensemble of over 30 publicly available models, which can be broadly grouped into three main categories, i.e., compartmental models, individual and network-based models, and statistical and machine learning models. Classic compartmental models, such as the SEIR (Susceptible, Exposed, Infectious, Recovered Individuals) model, represent a standard and widely used method in infectious disease epidemiology [44,45]. A SEIR model employs systems of either deterministic or stochastic ODEs to describe the dynamics of an epidemic. The pathogen transmission process is formulated as population transition/movement between “states” or “compartments”, comprising a dynamic system that is driven by epidemiological, biological, environmental and other related parameters at the community level. Due to their mechanism-based nature, compartmental models provide flexible frameworks for capturing epidemic dynamics, detecting potential resurgences, and exploring the efficacy of mitigation measures. Individual and network-based models, particularly agent-based models, expand on compartmental SEIR models by attempting to capture the complex interactions and behavior patterns of agents representing individuals in the populations of concern [46,47]. This approach can provide a more detailed description of disease transmission via a bottom-up framework, and potentially account for population heterogeneities that emulate “emergent” phenomena such as superspreading and spread within clusters. However, a reliable agent-based model generally poses excessive requirements for data inputs and parameters that need to be informative at local scales. Statistical and machine learning models use data-driven techniques aiming to predict future conditions by fitting curves defined through specified functions using early data [48], by capturing historical trends and correlations of time series [49], and by “learning” complex patterns and relationships between adverse outcomes and various underlying factors [50]. Statistical learning is valuable for modeling short-term trajectories of an epidemic; however, this approach alone cannot produce robust and reliable COVID-19 predictions without the support of mechanistic models, such as SEIR approaches. Statistical modeling has limited “extrapolation” ability (i.e., long-term prediction performance) due to its data-driven foundations and evaluations are typically conducted only for a few days ahead. The limitations of statistical modeling are especially true when attempting to describe COVID-19 spread and outcomes, since the data used to build the statistical model reflect conditions that are drastically changing temporally (due to changes in people’s behavior and responses, availability of treatments, possible virus adaptations, etc.) further curtailing its predictive ability. In contrast, mechanistic models can narrow the modeling/searching space by utilizing epidemiological knowledge to improve prediction performance [51].

As mentioned earlier, in the context of SEIR formulations, epidemic modeling frameworks can be either deterministic or stochastic [39]. The substantial uncertainties associated with the processes affecting COVID-19 dynamics, especially during the initial phases of the pandemic, severely limit the applicability of deterministic frameworks to “real world” studies of the pandemic. Stochastic frameworks are more appropriate for such studies as they can account for the inherent randomness and variability of the system and can generate the confidence intervals for each state variable via repeated simulations, starting with an estimated set of parameters and initial conditions. Monte Carlo approaches are then employed to consider the uncertainties of different model parameters within “most probable” domains. Instead of assuming homogeneous uniform mixing—as is the case with deterministic modeling (typically resulting in a system of ODEs)—stochastic frameworks can directly describe transmission heterogeneities via a hierarchical structure of uncertain parameters and simulate representative fat-tailed distributions that involve superspreader events [52]. Technical analyses of the relationships between the two frameworks can be found in the literature [53,54,55,56].

The majority of models discussed above have been predicting COVID-19 trajectories for countries or large areas within countries or states [57,58,59,60,61], an approach that does not account for important heterogeneities within these large geographical areas. Kain et al. [62] used a stochastic SEIR model to capture the heterogeneities of disease transmission from five selected counties/cities across the US and concluded that interventions truncating the individual-level transmission rate distribution while partially relaxing social distancing can be effective in maintaining epidemic control. Thomas et al. [63] used a network-based model to demonstrate that heterogeneity in population distribution can have large impacts on local pandemic timing and severity, based on simulations for 19 US cities. It should be mentioned that, although PDE-based approaches were not included in the CDC ensemble of models [43], they provide systematic ways to account for the effects of spatial heterogeneities on epidemic dynamics, and such approaches have been utilized in studies of the spread of COVID-19 in Arizona, US [64] and the Lombardy Region of Italy [65,66,67,68]. However, there has been a lack of modeling work focusing on evaluating the effectiveness of interventions that reduce exposures at local regions, while also exploring how key variables (e.g., population density, compliance with wearing face masks, prior exposures to various stressors, age stratification, comorbidities in the population, etc.) could affect disease dynamics and contribute to local heterogeneities of COVID-19 spread and outcomes. For example, the higher population density of urban settings causes social distancing to be more difficult to achieve than in suburban or rural areas and as a result populations of highly urbanized counties are expected to experience greater potential exposures, consequently resulting in higher numbers of cases and deaths. The deaths, however, are also related to the availability and quality of health care and hospitals, and these may be better in urban than in rural settings.

1.3. New Jersey as a “Microcosm” of COVID-19 Spread Heterogeneities

New Jersey has been one of the earliest and largest hotspots for COVID-19 in the United States, ranking first in the nation in per capita death rates until September 2021 (with over 300 deaths per 100,000 residents [69]) when it was surpassed by Mississippi. As mentioned earlier, in addition to individual-level risk factors, such as age, sex and underlying medical conditions, it has been widely recognized that multiple demographic, environmental and socioeconomic factors are strongly correlated with the patterns and severity of COVID-19 spread. In particular, racial and ethnic minorities, economically disadvantaged populations, and environmental justice communities have been disproportionally impacted by COVID-19. For example, over 47% of individuals with confirmed positive tests in NJ during the first wave of the pandemic were Black and Hispanic, while these minorities constitute approximately 31% of the State’s population; furthermore, the NJ counties with the highest rates of COVID-19 deaths (Essex, Union, Passaic, Hudson, Bergen) have historically high levels of hazardous air pollutants, such as airborne diesel exhaust particles, with concentrations that rank above the 98th percentile of national values. In fact, New Jersey is the State with the highest population density in the nation, with approximately 1208 people/mile², while population densities of individual counties range from less than 200 people/mile² (Salem County) to more than 14,600 people/mile² (Hudson County). New Jersey is also one of the most ethnically and socioeconomically diverse States; furthermore, urban centers, suburban sprawl, an aging industrial infrastructure, active agriculture, extensive preservation areas and a densely populated coast, all in close proximity to each other, constitute a remarkably heterogeneous environmental landscape (Figure S1). Therefore, New Jersey is a veritable “microcosm” that is representative of much of what is happening across the entire contiguous US. As such, understanding the factors impacting the dynamics of COVID-19 across New Jersey can provide valuable insight and quantitative information that is directly applicable to a wide range of locations across the United States. It should be noted that the majority of casualties in New Jersey are associated with the first wave of the pandemic (March to September 2020), before effective pharmaceutical interventions became available and while uncertainties were persisting regarding both the significance of different exposure pathways and the proper medical protocols for treating hospitalized cases. This situation posed significant challenges with respect to both predicting trajectories of the disease across different regions of the State and evaluating the relevance and efficacy of intervention strategies for reducing exposures. Meeting these challenges required a flexible modeling framework, combining mechanistic description of stochastic epidemic dynamics with computationally efficient data science methods for parameter optimization.

In this article we present the development and deployment of a stochastic SEIR modeling platform that was designed to characterize and quantitatively assess the dynamics of the first wave of the pandemic for each individual county across New Jersey. To better quantify the effects of layered exposure interventions that are affected by spatial heterogeneities, we used state-of-the-art statistical inference tools and intelligent selection algorithms in a Sequential Quasi-Monte Carlo (SQMC) framework to reliably capture epidemic dynamics [70,71,72,73]. Model predictions, as well as other local information, for counties and municipalities across New Jersey, were made available to the public since the early days of the pandemic, via an online dashboard (https://ccl-eohsi.shinyapps.io/covid19_dashboard/ (accessed on 6 October 2021)). The findings of this work can provide guidance to local agencies for planning for future epidemics of infectious diseases, as it supports evaluation of mitigation strategies that can minimize potential exposures while allowing certain essential communal and economic activities.

2. Materials and Methods

2.1. Stochastic SEIR Model

A stochastic compartmental SEIR (Susceptible, Exposed, Infectious, Recovered Individuals) modeling system with time-varying probabilistic transmission parameters [45] was implemented in a Sequential Quasi-Monte Carlo framework to simulate local COVID-19 spread dynamics at the county level across New Jersey. The present model divides the residents of each county into ten compartments (Figure 1). Susceptible individuals (S) become exposed (E, infected but not yet infectious) through contacts with individuals from the infectious compartments comprised of presymptomatic, asymptomatic, and symptomatic individuals (I_P + I_A + I_S + I_M). Some of the exposed individuals become infectious but not yet symptomatic; these presymptomatic individuals (I_p) become either mildly symptomatic (I_M) or severely symptomatic (I_S). For simplicity, it is assumed that individuals with mild symptoms can recover (R) without medical treatment and hospitalization, while those with severe symptoms should be treated/hospitalized (T/H_D + T/H_R). The treated/hospitalized individuals then progress to either death (D) or recovery (R) with a treatment/hospitalization fatality rate that can change over time. Because not all infected individuals show symptoms, some individuals progress to the asymptomatic compartment (I_A) and to recovery (R) in the end. The above transmission dynamics constitute stochastic Markov Processes [74,75], where each compartment represents a discrete stochastic state. Transitions of individuals between compartments (equivalently, the stochastic states) are modeled probabilistically using binomial or multinomial distributions.

The mathematical relationships governing the dynamics of the stochastic SEIR model are listed in Equation (1): C_A, C_P, C_M and C_S denote the relative infectiousness of asymptomatic, presymptomatic, mildly symptomatic, and severely symptomatic individuals, respectively; λ_I, λ_P, λ_A, λ_S and λ_M are rate constants representing the reciprocal of infectious periods for states E, I_P, I_A, I_S and I_M, respectively; ρ_R and ρ_D are rate constants representing the reciprocal of duration from treatment/hospitalization to either recovery or death, respectively; α, μ and HFR denote the proportion of asymptomatic infections, symptomatic infections that do not require treatment/hospitalization and treatment/hospitalization fatality, respectively, where HFR is assumed to have a linear decreasing trend for the first wave, satisfying HFR = HFR₀ − k∙t. The time-varying transmission rate is denoted as β_t = β₀∙γ_t, where β₀ is a constant representing the basic transmission rate, and γ_t is the exposure control-based intervention scaling factor at time t. The binomial and multinomial distributions are denoted as ℬ(∙) and ℳ(∙), respectively. The differential operator is “d”, so for example, dSE denotes the number of individuals transitioning from the Susceptible compartment/state to the Exposed compartment/state per time unit.

(1)$\begin{array}{l}\mathrm{d}SE~\mathcal{B}\left(S,1-\exp \left(-{\beta }_t\frac{C_A{I}_A+C_P{I}_P+C_M{I}_M+C_S{I}_S}{N}\right)\right)\\ \left(\begin{array}{c}\mathrm{d}EE\\ \mathrm{d}E{I}_A\\ \mathrm{d}E{I}_P\end{array}\right)~\mathcal{M}\left(E,\left(\begin{array}{c}\exp (-{\lambda }_I\mathrm{d}t)\\ \left(1-\exp (-{\lambda }_I\mathrm{d}t\right))\alpha \\ \left(1-\exp (-{\lambda }_I\mathrm{d}t\right))\left(1-\alpha \right)\end{array}\right)\right)\\ \mathrm{d}{I}_{A}R~\mathcal{B}\left(I_A,1-\exp (-{\lambda }_A\mathrm{d}t)\right)\\ \left(\begin{array}{c}\mathrm{d}{I}_P{I}_P\\ \mathrm{d}{I}_P{I}_M\\ \mathrm{d}{I}_P{I}_S\end{array}\right)~\mathcal{M}\left(I_P,\left(\begin{array}{c}\exp (-{\lambda }_P\mathrm{d}t)\\ \left(1-\exp (-{\lambda }_P\mathrm{d}t\right))\mu \\ \left(1-\exp (-{\lambda }_P\mathrm{d}t\right))\left(1-\mu \right)\end{array}\right)\right)\\ \mathrm{d}{I}_{M}R~\mathcal{B}\left(I_M,1-\exp (-{\lambda }_M\mathrm{d}t)\right)\\ \left(\begin{array}{c}\mathrm{d}{I}_S{I}_S\\ \mathrm{d}{I}_S{H}_R\\ \mathrm{d}{I}_S{H}_D\end{array}\right)~\mathcal{M}\left(I_S,\left(\begin{array}{c}\exp (-{\lambda }_S\mathrm{d}t)\\ \left(1-\exp (-{\lambda }_S\mathrm{d}t\right))HFR\\ \left(1-\exp (-{\lambda }_S\mathrm{d}t\right))\left(1-HFR\right)\end{array}\right)\right)\\ \mathrm{d}{H}_{R}R~\mathcal{B}\left(H_R,1-\exp (-{\rho }_R\mathrm{d}t)\right)\\ \mathrm{d}{H}_{D}D~\mathcal{B}\left(H_D,1-\exp (-{\rho }_D\mathrm{d}t)\right)\\ \frac{\mathrm{d}S}{\mathrm{d}t}=-\mathrm{d}SE\\ \frac{\mathrm{d}E}{\mathrm{d}t}=\mathrm{d}SE-\mathrm{d}E{I}_A-\mathrm{d}E{I}_P\\ \frac{\mathrm{d}{I}_A}{\mathrm{d}t}=\mathrm{d}E{I}_A-\mathrm{d}{I}_{A}R\\ \frac{\mathrm{d}{I}_P}{\mathrm{d}t}=\mathrm{d}E{I}_P-\mathrm{d}{I}_P{I}_M-\mathrm{d}{I}_P{I}_S\\ \frac{\mathrm{d}{I}_M}{\mathrm{d}t}=\mathrm{d}{I}_P{I}_M-\mathrm{d}{I}_{M}R\\ \frac{\mathrm{d}{I}_S}{\mathrm{d}t}=\mathrm{d}{I}_P{I}_S-\mathrm{d}{I}_S{H}_R-\mathrm{d}{I}_S{H}_D\\ \frac{\mathrm{d}{H}_R}{\mathrm{d}t}=\mathrm{d}{I}_S{H}_R-\mathrm{d}{H}_{R}R\\ \frac{\mathrm{d}{H}_D}{\mathrm{d}t}=\mathrm{d}{I}_S{H}_D-\mathrm{d}{H}_{D}D\\ \frac{\mathrm{d}R}{\mathrm{d}t}=\mathrm{d}{I}_{A}R+\mathrm{d}{I}_{M}R+\mathrm{d}{H}_{R}R\\ \frac{\mathrm{d}D}{\mathrm{d}t}=\mathrm{d}{H}_{D}D\end{array}$

The biological and epidemic-relevant parameters that determine the compartmental transition probability and disease transmission dynamics are provided in the Supporting Information. These parameters were specified with either point (Table S1) or range estimates (Table S2), to account for various uncertainties associated with the dynamics of the first wave of the COVID-19 pandemic. Exposure control-based intervention strategies and public health policies are important factors affecting the transmission parameters, forcing them to vary over time. During early March, 2020, telecommuting or work-from-home arrangements were suggested/implemented for all businesses and non-profits that could operate in such a mode in the State. Following a series of mitigation measures (such as closures of educational facilities and non-essential services), the Governor issued a directive on 21 March, that residents should stay at home. When the Centers for Disease Control and Prevention (CDC) recommended wearing face masks to limit the spread of COVID-19 in early April 2020, the New Jersey Governor’s office implemented a face masks directive starting on April 11. In the model, we considered all of the above mitigation measures and policies. We divided the timeline into four intervention periods, i.e., non-intervention, work from home, shelter-in-place, and implementation of face masks directive, for estimating the effects of those measures on the spread of COVID-19 (Figure S4).

The treatment/hospitalization fatality rate (T/HFR) is an important metric for evaluating disease severity of COVID-19. A large T/HFR (e.g., 5–28%) was associated with the initial rise in the epidemic. After the initial period, as improved medical care protocols started becoming available, the T/HFR decreased [76]. Figure S5 depicts an approximate estimate of T/HFR in New Jersey, that decreased initially and then remained within a stable range. To account for the changing T/HFR that can affect COVID-19 death trajectories, we assumed a time-varying T/HFR beginning in mid-April (uncertain interval: 4/11–4/15). Note that April 14 is the date when the hospitalizations reached a peak; as such, we used a linear decreasing function, setting 1% as the lower bound of the estimated T/HFRs, to characterize such reduction and embedded it in the stochastic SEIR model.

2.2. Model Calibration

The stochastic SEIR model was calibrated using the reported daily numbers of new confirmed COVID-19 deaths for each NJ county, a metric that is substantially more reliable than other available data, such as the daily reported cases. The death data were based on reports from the New Jersey Department of Health (NJDOH) [77] and were compared with numbers reported in national repositories (NY Times, Atlantic’s COVID Tracking Project, Johns Hopkins). It should be noted that in late June 2020, New Jersey started reporting “probable” cases and deaths. A probable COVID-19 case/death is defined based on evidence from clinical, epidemiological or serological testing, or from vital records, but without a confirmatory laboratory RNA test. We excluded the probable deaths in implementing the model simulations in order to avoid discontinuities and fluctuations in the data and reduce related uncertainties. We adopted a SQMC framework for the optimization of the stochastic dynamic SEIR models, as such frameworks have been demonstrated to be more accurate and robust than other state/parameter estimation methods, especially for nonlinear and non-Gaussian dynamic systems [78]. The SQMC framework takes advantage of advanced filtering and sampling techniques (e.g., Bayesian Filtering and Sequential Monte Carlo) to “predict-and-correct” states recursively, in a manner consistent with the known dynamic model structure and the available observed data. In this study, an improved filtering algorithm [79], called Iterated Filtering [75], was used to optimize four parameters of the stochastic SEIR model, specifically a basic transmission rate ${\widehat{\beta }}_0$, a shelter-in-place transmission scaling factor ${\widehat{\gamma }}_{\mathrm{SIP}}$, a face masks directive transmission scaling factor ${\widehat{\gamma }}_{\mathrm{FMD}}$, and a treatment/hospitalization fatality rate slope k_T/HFR. The work from home transmission scaling factor ${\gamma }_{\mathrm{WFH}}$ was not estimated but was drawn randomly from a specified range (Table S2).

In order to consider the uncertainties associated with different model parameters, we generated 2000 parameter sets (low discrepancy Sobol sequences [45]) by “randomly” drawing samples from the specified parameter ranges. We constructed the model based on each Sobol sequence and calibrated the corresponding 2000 SEIR sub-models for each of the 21 counties of the State of New Jersey. The statistical inference and estimation algorithm has certain intrinsic limitations that in some cases may cause optimization failure (with non-convergent iteration), trapping into local optimum (with few initial values and particles), etc. To obtain reliable estimates for each calibration, we used 100 particle-filter iterations, 3000 particles and 5 particle filtering simulations [75]: these algorithm parameter settings ensure optimized performance while maintaining a practical computational load. We selected a time step of 4 h to ensure stable estimation for each state variable, and as such the daily predictions (number of deaths) are the sum of six estimates calculated in each day.

2.3. Simulation

The optimized parameters, along with the corresponding Sobol sequences, were combined to run multiple sub-simulations (2000) for each of the 21 NJ counties, where each sub-simulation was repeated multiple times (200) to account for stochasticity in state transitions. We considered three scenarios to assess the effects of mitigation measures on disease transmission: (a) Baseline, a counterfactual scenario using the estimated initial basic transmission rate ${\widehat{\beta }}_0$ throughout the epidemic; (b) Social Distancing without Face Masks, a counterfactual scenario using the estimated shelter-in-place transmission rate ${\widehat{\beta }}_0\cdot {\widehat{\gamma }}_{\mathrm{SIP}}$ after March 21; and (c) Social Distancing with Face Masks, a “real world” scenario using the estimated shelter-in-place transmission rate ${\widehat{\beta }}_0\cdot {\widehat{\gamma }}_{\mathrm{SIP}}$ from March 21 to mid-April, and using the face masks directive transmission rate ${\widehat{\beta }}_0\cdot {\widehat{\gamma }}_{\mathrm{FMD}}$ after mid-April. A time-varying treatment/hospitalization fatality rate T/HFR was also included for the latter two scenarios (Figure S4).

The predicted death trajectories are the result of aggregated sub-simulations performed using the sets of different calibrated models with corresponding Sobol sequences. It should be noted that not all sub-simulations are suitable for evaluation. The reasons are twofold: (a) using relatively wide parameter ranges causes some Sobol sequences to deviate substantially from the underlying true parameter values. Accordingly, the optimized parameters based on these “biased” Sobol sequences, as well as the “biased” sub-models, do not reflect the exact epidemic dynamics; (b) due to the intrinsic limitations of the optimization algorithm, the estimated parameters from some sub-models do not represent “global” optima. These “poor” estimates tend to fail to capture the exact epidemic dynamics, so the need arises to systematically select the “proper” Sobol sequences. Kain et al. [62] selected Sobol sequences using a likelihood metric; however, that metric cannot adequately reflect the performance of the calibration process and may become invalid or even misleading, because it is calculated based on certain ideal distributions. For this reason, we designed an efficient algorithm for selecting the “most suitable” Sobol sequences and associated sub-models: we first select the “qualified” Sobol sequences that correspond to good predictions of both the daily new deaths and the total deaths (e.g., below an error tolerance with 30% of the maximum); we then use expert knowledge and screening rules (if necessary) to target the “most suitable” Sobol sequences (~10) located within the most frequent/stable intervals (Figure S6). The intensive computation for optimization and exploration of the global optimum, as well as the efficient algorithm used to select the “most suitable” Sobol sequences are crucial for capturing robust and reliable epidemic dynamics across different areas.

3. Results

The spatial domain considered in this study covers the entire State of New Jersey (with 21 counties), where five clustered geographic regions were identified to better explore the spatial heterogeneities of disease transmission (Figure 2).

3.1. Spatiotemporal Analysis of COVID-19 Deaths across New Jersey

Table 1 provides a summary overview of the COVID-19 death rates during the first wave of the pandemic for each of the 21 New Jersey counties, as well as for the five geographic regions identified in the map of Figure 2. The death rate is defined as the number of the confirmed deaths per 100,000 residents during a specified time period. The peaks occurred mainly in mid-April 2020 for areas near the epicenter of New York City (e.g., Gateway Region), while there was a two to four weeks delay for those in Southern New Jersey (e.g., Delaware River Region, Southern Shore Region). Until August 2020, areas in southern New Jersey had death rates ~40% lower than those in the north. Figure 3a,b show the spatial distribution of the death rates calculated for April and August 2020, respectively. The local variations of the death rates in April 2020 are significantly associated with the proximity to the “epicenter”, that is, the closer a county is to New York City (NYC), the higher was the death rate. However, this pattern became more complex by August 2020: for instance, counties far from the NYC “epicenter” may have larger death rates than those counties that are nearer (e.g., Salem and Gloucester in Figure 3a,b), indicating that the effects from other factors, such as population density, policy compliance, etc., became more prominent as the COVID-19 spread progressed.

3.2. Comparison of Predicted and Reported Confirmed Deaths

Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 present the predicted trajectories of the daily new/total confirmed deaths (with median and 95% confidence interval) for the three scenarios considered, along with the reported numbers of the confirmed deaths for each county. The same scale was used for the plots for counties within the same geographic region, in order to better show between/within-region variations. Models were calibrated with the daily new confirmed deaths before 2 August 2020 and were then tested for 8 days and 52 days ahead, i.e., 3 August to 10 August to 30 September. As mentioned earlier, regularly updated model predictions were made available to the public throughout the first wave of the pandemic via the dashboard: https://ccl-eohsi.shinyapps.io/covid19_dashboard/ (accessed on 6 October 2021).

Three statistical measures, i.e., root mean squared error (RMSE), mean absolute percentage error (MEPE), and envelope bias (EB), were used to evaluate model performance (Table S3). There is very good agreement between the predicted and reported numbers of total deaths. The reported total confirmed deaths are all within the 95% confidence interval for each county (except for Mercer with EB = 1 day), with MAPE ranging from 9% to 24%. The reported daily new deaths show inter-day fluctuations that, as expected, cannot be captured by the averaged estimates of the model, but most of which are within the 95% confidence interval. The RMSE of the daily new deaths varies for different regions, with the highest in the Gateway Region (larger cardinal death number); the RMSE in the rising phase of the first wave of the epidemic was generally lower than that in the falling phase.

3.3. Realistic and Counterfactual Intervention Scenarios

Figure S7 shows the predicted median of the new/total confirmed deaths from the 21 New Jersey counties for each of the three scenarios, i.e., Baseline (scenario 1), Social Distancing without Face Masks (scenario 2), and Social Distancing with Face Masks (scenario 3). When considering the situation of practicing social distancing but without wearing face masks (scenario 2 vs. scenario 1), the model predicts that there is a 16 [IQR (14, 21)] days delay of peaks, with 179 [IQR (115, 212)] deaths averted per 100,000 people. The effectiveness of social distancing in “flattening the curve” results from the significant reduction of disease transmission/exposure, according to the estimates of transmission rates shown in Figure S8. Adopting the face masks recommendation (scenario 3 vs. scenario 2), is shown to have saved an additional 169 [IQR (116, 281)] lives per 100,000 people compared to just using social distancing, and to have resulted in a peak delay of three [IQR (−2, 5)] days compared to the baseline scenario. The concurrent presence of social distancing and wearing face masks produce synergistic effects that considerably reduce deaths while slightly delaying the peaks; see Table S4 for the evaluation of effects for each NJ county.

4. Discussion

4.1. Spatially Heterogeneous Transmission Rate

The basic transmission rate (β₀), a measure of the contacts per unit time (contact rate) multiplied by the likelihood of infection per contact, is a key component of the force of infection. Sy et al. [80] reported positive correlations between population density and the rates of transmission of SARS-CoV-2. However, this study was based on a statistical analysis (employing a linear mixed model) that cannot reveal the underlying mechanistic patterns. Our estimates for the transmission rates, optimized with the county-level SEIR models, reveal a nonlinear heterogeneous relationship to the population density (Figure 9). The basic transmission rate initially increases at lower population densities and then becomes saturated at higher densities (>2000 people/mile²). Such patterns reflect the combined effects of “density-dependent” (low density) and “frequency-dependent” (high density) mechanisms between contact rate and population density [81]. It should be pointed out that this trend is consistent with a mechanistic nonlinear scaling function (dashed line in Figure 9) that was derived from a spatial contact model which consider contacts of individuals within a population and is analogous to models that are widely used in kinetic theory [82]. Our estimation for the county-level transmission rates provides evidence on the significance of distinguishing density-dependent and frequency-dependent conditions, for better assessing the spatiotemporal patterns of the COVID-19 dynamic with partial differential equation models [83].

The transmission rate is also related to geographic location. Larger transmission rates were observed in locations closer to the epicenter of New York City (Figure S8a), with the highest rate occurring in the Gateway Region (1 to 1.3) and the lowest in the Southern Shore Region (0.4 to 0.7). Susceptible individuals living closer to the epicenter probably had a larger number of contacts with infectious individuals, due to both proximity to NYC and the higher population density of counties comprising this region. In practice, transmission rates are affected by different demographic (age, type of employment, income, education, ethnicity/race, etc.) and sociobehavioral factors, and the resulting spatially varying patterns are one of main reasons causing local heterogeneities in the spread and outcomes of COVID-19. Modeling at higher resolutions provides more flexible model structures which allow for capturing these local variations, and, moreover, improving the prediction performance at multiple scales. Figure 10 shows that the statewide death trajectory produced by combining outcomes from the 21 county-level SEIR models, considering heterogeneous transmission rates, is more reliable than the trajectory produced by a State-level SEIR model with a homogeneous transmission rate across the State.

4.2. Reproduction Number at County Level

The basic reproduction number (R₀) is a fundamental epidemiologic metric used to describe the contagiousness of an infectious agent in the absence of interventions [84]: It is defined as the average number of secondary cases arising from a typical primary case in a completely susceptible population, and calculated as the product of the basic transmission rate and the duration of contagiousness for an infected person. The R₀ for the present stochastic SEIR model is calculated via Equation (2):

(2)${\mathrm{R}}_0={\beta }_0\cdot \left[\alpha \cdot C_A\cdot \frac{1}{{\lambda }_A}+\left(1-\alpha \right)\cdot \mu \cdot \left(\frac{1}{{\lambda }_P}+\frac{1}{{\lambda }_M}\right)+\left(1-\alpha \right)\cdot \left(1-\mu \right)\cdot \left(\frac{1}{{\lambda }_P}+\frac{1}{{\lambda }_S}\right)\right]$

where the terms in the square bracket represent the approximated duration of contagiousness (to recovery) that “averages” the effects of asymptomatic, mildly symptomatic, and severely symptomatic transmission. Early studies [85,86,87] reported an average estimated R₀ for COVID-19 that was over 3 in multiple countries/regions with some estimates close to 6. Our estimates show significant variation across New Jersey, with R₀ ranging from 3 to 7 (Figure 3c). Such variation is related to the spatial heterogeneity of the estimated transmission rates across New Jersey, as discussed in the previous subsection, as well as to the local uncertainties of the duration of contagiousness that are reflected in the values of the biological/epidemic parameters of the model, such as the relative infectiousness of asymptomatic individuals, and the proportions of asymptomatic, severely symptomatic, etc. individuals within the population (Equation (2), Table S2). Caution should be taken when comparing with estimates from previous studies, since R₀ can vary considerably for different countries/regions [88]. It should also be pointed out that R₀ is a model-dependent parameter and a proper comparison of R₀ values should be based on the same model structure and assumptions [84].

The effective reproduction number (R_eff) is a time-varying epidemiologic metric, defined as the actual average number of secondary cases arising from a typical primary case in a population that may have individuals with immunity and may have intervention measures in place at a specific time point [89]. In the present study, we assumed that individuals in the recovered state/compartment can directly gain immunity. The R_eff for the present stochastic SEIR model is calculated through Equation (3):

(3)${\mathrm{R}}_{\mathrm{eff}}=frac\_S\cdot \gamma \cdot {\mathrm{R}}_0$

where frac_S is the fraction of the susceptible in time t; γ is the transmission scaling factor either specified in Table S2 or optimized via the SQMC framework, according to the type of the exposure control-based intervention strategies; R₀ is the basic reproduction number calculated in Equation (2). Previous studies mainly used statistical methods [90,91,92] to estimate R_eff in a fast/timely manner for evaluating the daily contagiousness of COVID-19. In this work, we used the stochastic SEIR model to calculate the R_eff as the product of R₀, the fraction of the susceptible population, and the transmission scaling factor given the interventions at time t (Equation (3)). The trend of our estimated R_eff (green lines in Figure S9) is consistent with that reported by NJDOH [90]: our county estimates of R_eff in May and June 2020 range from 0.18 to 0.86, compared to the estimates (0.6–0.9) at the state level by NJDOH. Figure S9 presents the time-varying R_eff for each of the 21 New Jersey counties for the three scenarios considered. The R_eff values for scenario 1 (red lines) decrease in a fast and smooth manner due to the rapid reduction of the susceptible population, and then converge to a small value (below 0.1). In scenario 2 (purple lines) and scenario 3 (green lines), the R_eff values decrease according to the strengths of the mitigation measures and the proportion of the susceptible population in each county.

4.3. Reducing COVID-19 Fatality by Layered Exposure-Relevant Interventions

Layered interventions, via interruption of multiple transmission routes and pathways, are efficient in reducing exposures to SARS-CoV-2. To control the epidemic, New Jersey implemented different levels of social distancing since March 2020: keeping six feet from others can prevent transmission from direct droplet spray contact; school closures, mass gatherings bans, and shelter-in-place directives also directly reduce the contact rate. The recommendation to wear face masks in mid-April further reduced the likelihood of infection per contact, by limiting the transfer of viral droplets and aerosols [93]. In essence, the goal of the combined interventions was to scale (reduce) the transmission rates of COVID-19, as expressed in the SEIR model. The modeling results show that social distancing, particularly the shelter-in-place order, considerably reduced the county-level transmission rates from 0.92 [IQR (0.61, 1.07)] (Figure 9 and Figure S8a) to 0.30 [IQR (0.15, 0.55)] (33% [IQR (24%, 51%)] reduction) (Figure S8b), and that the face masks directive further reduced the transmission rates to 0.20 [IQR (0.10, 0.41)] (65% [IQR (50%, 75%)] reduction) (Figure S8c). It should be noted that “social” (physical) distancing involves compliance with a series of government directives and guidelines mentioned above, though shelter-in-place is expected to play a dominant role in reducing the transmission rate before the implementation of face masks directive.

Population density not only affects the (basic) transmission rate, as discussed in the previous subsection, but also affects the effectiveness of the interventions/policies that scale the transmission rate. According to the counterfactual analysis for scenario 2, i.e., Social Distancing without Face Masks, counties where more deaths per capita were averted by social distancing are those with lower population density (Figure S10a). This can be explained by the fact that social distancing is usually more difficult to achieve in densely populated areas. By comparing scenario 2 and scenario 3 (Social Distancing with Face Masks), it is concluded that the face masks directive was more effective in areas with higher population density, i.e., the difference in averted deaths per capita (Figure S10b) between scenario 3 and scenario 2 is more negative, and the corresponding ratio (Figure S10c) is larger than 1. This could indicate that residents of more densely populated areas tend to be more concerned and exhibit higher compliance [94] with wearing face masks. Another possible explanation is that when social distancing can be accomplished in low population density areas, the effectiveness of a mask is lower, since emitted droplets and aerosols do not commonly reach a large number of individuals. The other issue is how people travel to work: in dense cities, there is substantially higher use of mass transit where masks are required and are effective, while in rural areas private automobiles are more likely used and masks would not offer additional benefits.

In addition to the transmission rate, the fatality rate is another crucial metric that determines death trajectories. We modeled the decline of treatment/hospitalization fatality rate (T/HFR) when the hospitalizations reached their peak in New Jersey. The estimates for the decline slope of the T/HFR curve is −0.07% [IQR (−0.09%, −0.05%)], consistent with the estimate of −0.06% reported by CEBM [76]. Such a decline (presumably due to improved treatment protocols) saved 38 [IQR (28, 60)] lives per 100,000 people, a much lower number than the number of lives saved by layered interventions, i.e., 357 [IQR (290, 429)]. This indicates that layered interventions, aiming to mitigate exposures to SARS-CoV-2, have been more effective in reducing mortality than medical treatments available during the first wave of the pandemic.

4.4. Limitations

There are certain limitations to the study presented here. First, this version of the stochastic SEIR model does not account explicitly for heterogeneities due to age disparities in susceptibility and transmission patterns [95]; incorporating age stratification should further reveal patterns in transmission/fatality rates (Figure 9 and Figure S10) with respect to population density and relevant demographic factors. Though we have shown the benefits of prediction at finer spatial scales (Figure 10), in the present work, transmission of COVID-19 was modeled independently for each county without considering the effects of inter-county interactions; incorporating relevant information (such as population mobility) could have further improved the assessments of spatial transmission [64,83,96,97]; however, it should be recognized that the data-driven parameter estimation and optimization implicitly captures the effects of inter-county differences such as “differentials” in community levels. Second, the model calibration process considers key elements of transmission dynamics for different intervention measures and timing, while for simplicity, it assumes fixed transmission rates after the implementation of those mitigation measures. The current model was calibrated with the number of the reported confirmed deaths (with less fluctuation and uncertainty) and evaluated using the same type of the data. Numbers of projected new/total cases can be generated but would be expected to be much higher than numbers of reported cases. Using reported case data in model calibration to predict and assess case trajectories for each of the 21 NJ counties would have required considering additional parameters (such as COVID-19 testing rates) associated with great uncertainty and variability, particularly in the early phase of the pandemic. Comparisons and analyses of the predicted case trajectories were therefore not a main focus of this study. Third, it is known that a large proportion of deaths in New Jersey and across the nation particularly in the beginning of the epidemic, were associated with Long Term Care (LTC) facilities. Due to the unavailability of relevant data, when the model was implemented in the first months of the pandemic, an “average” fatality rate that did not distinguish deaths that occurred at hospitals from those that occurred at nursing home was used. Given that spread of COVID-19 in LTC facilities exhibited more rapid patterns, it could have been described separately from the spread occurring in populations surrounding around these facilities [77]. However, the stochastic SEIR framework that is presented here can actually model the effects of “superspreading” events, such as those that occurred in LTC facilities, by probabilistically considering individual-level heterogeneities in transmission as the stochastic SEIR model can capture “superspreader” events in the tails of the skewed distributions. It should be noted that the current model assumes a single average transmission rate which may lead to overestimation of effectiveness of some interventions in LTC facilities (e.g., social distancing is difficult to achieve in LTC due to residents’ frailty and close living quarters). However, the estimated population-average county transmission rates are significantly associated with the number of LTC beds (Pearson coefficient 0.68) in the corresponding counties: these parameters can provide crucial information for further optimization that would consider additional effects of heterogeneities within each county. Overall, COVID-19 spread is a very complex process that is affected by multiple demographic, socioeconomic/behavioral, and environmental factors. Tiered modeling can be useful in extracting complex mechanisms consistent with alternative hypotheses, while it should be complemented by statistical analysis involving large numbers of exposure-relevant factors [98], a task that is currently an on-going next step of this work [99]. Better understanding of transmission dynamics and characterization of critical parameters/factors can help guide individual and government decision making aiming to mitigate the future trajectory of any epidemic.

5. Conclusions

This article presents a hybrid application of systems science, employing a mechanistic stochastic SEIR modeling framework, combined with computational data science, using a Sequential Quasi-Monte Carlo framework for optimal estimation of multiple model parameters, for quantifying the dynamics of COVID-19 within each of the 21 counties of New Jersey. This work applied stochastic epidemiological modeling to multiple heterogeneous jurisdictions within a single region, where the same mitigation policies had been imposed, to investigate the effectiveness and compliance of layered exposure-relevant interventions (such as social distancing and face mask use) on epidemic control. The modeling results in fact revealed significant local heterogeneities of COVID-19 dynamics across New Jersey; the captured spatiotemporal patterns were found to be related to a wide range of factors, such as population density, proximity to the epicenter of the epidemic, and compliance with wearing face masks.

We estimated transmission rates varying both in space and time and derived a nonlinear relationship with respect to population density that is consistent with mechanistic predictions analogous to those derived in kinetic theory. Population density was also found to be related to the effectiveness and compliance of policies that “scale” the epidemic transmission rate. The modeling results indicate that social distancing and wearing face masks were effective in reducing exposures to SARS-CoV-2, averting an estimated 357 deaths [IQR (290, 429)] per 100,000 people while slightly delaying the peak by three [IQR (−2, 5)] days for the first wave in New Jersey. Furthermore, before effective pharmacological interventions became available, gradual improvements of treatment protocols during the first wave saved lives (38 [IQR (28, 60)] per 100,000 people) but at a much lower rate than of the lives saved by layered exposure interventions. The present computational modeling framework is currently being expanded to the study of socioeconomic, behavioral and environmental factors on pandemic dynamics [99] and to quantification and assessments of pharmaceutical interventions, such as vaccination and loss of immunity, that are currently of public concern [100,101,102].

Author Contributions

Conceptualization, X.R. and P.G.G.; methodology, X.R. and P.G.G.; software, X.R. and P.G.G.; formal analysis, X.R.; writing—original draft preparation, X.R.; writing—review and editing, C.P.W. and P.G.G.; visualization, X.R.; project administration, P.G.G.; funding acquisition, P.G.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Center for Environmental Exposures and Disease at EOHSI (NIEHS Grant P30ES005022) through a Pilot Grant. Additional funding was provided by the State of New Jersey Department of Environmental Protection (Grant AQ05-011) and by the NJ Alliance for Clinical Translational Science (NIH Grant UL1TROO3017). The funding agencies have not reviewed and are not responsible for the contents and the conclusions of this manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

An online dashboard for this project is accessible at https://ccl-eohsi.shinyapps.io/covid19_dashboard/ (accessed on 6 October 2021). Data on COVID-19 cases and deaths are publicly available on the following URL: https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-counties.csv (accessed on 6 October 2021).

Acknowledgments

We gratefully acknowledge the Office of Advanced Research Computing (OARC) at Rutgers The State University of New Jersey for providing unlimited access to the Caliburn supercomputing cluster. We would also like to thank Zhongyuan Mi (Rutgers, GIS and dashboard support), Galen Collier (Rutgers, computing cluster support), and M. P. Kain et al. (Stanford) for making their SEIR model code (https://github.com/morgankain/COVID_interventions (accessed on 6 October 2021)) available for further development. In addition, we acknowledge multiple organizations sharing/exchanging information on COVID-19, including: CDC (COVID-19, WONDER), NJDOH, the Models of Infectious Disease Agent Study (MIDAS), the JHU CSSE dashboard, the COVID Tracking Project (The Atlantic), the NY Times GitHub repository; Christoph Schoenenberger; the Yu group (Berkeley) and the Mordecai group (Stanford).

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Block diagram of the stochastic compartmental SEIR model. Population “movements” between compartments (blocks representing stochastic states) follow binomial (B) or multinomial (M) distributions that are controlled by biological/epidemic parameters as marked on the arrows. S = susceptible, E = exposed, Ip = presymptomatic, IA = asymptomatic, Is = severely symptomatic, IM = mildly symptomatic, T/HD = treated/hospitalized to death, T/HR = treated/hospitalized to recovery, R = recovered, D = dead. β = transmission rate. Definitions of other parameters are listed in Table S1.

Figure 2. The 21 Counties and the five geographical regions identified based on similarities of exposure-relevant factors across the State of New Jersey.

Figure 3. Spread of the COVID-19 epidemic across New Jersey: (a) death rates (confirmed deaths per 100,000 people) by the end of April 2020, (b) death rates by the end of August 2020, and (c) basic reproduction number (R0) at the county level.

Figure 4. Predicted trajectories of new/total confirmed deaths from COVID-19 corresponding to three exposure scenarios for the counties in the Gateway Region of NJ, March to September 2020, compared with reported death data (vertical bars). The color-shaded areas depict the 95% Confidence Interval (CI) of the corresponding predictions.

Figure 5. Predicted trajectories of new/total confirmed deaths from COVID-19 corresponding to three exposure scenarios for the counties in the Shore Region of NJ, March to September 2020, compared with reported death data (vertical bars). The color-shaded areas depict the 95% Confidence Interval (CI) of the corresponding predictions.

Figure 6. Predicted trajectories of new/total confirmed deaths from COVID-19 corresponding to three exposure scenarios for the counties in the Southern Shore Region of NJ, March to September 2020, compared with reported death data (vertical bars). The color-shaded areas depict the 95% Confidence Interval (CI) of the corresponding predictions.

Figure 7. Predicted trajectories of new/total confirmed deaths from COVID-19 corresponding to three exposure scenarios for the counties in the Skyland Region of NJ, March to September 2020, compared with reported death data (vertical bars). The color-shaded areas depict the 95% Confidence Interval (CI) of the corresponding predictions.

Figure 8. Predicted trajectories of new/total confirmed deaths from COVID-19 corresponding to three exposure scenarios for the counties in the Delaware River Region of NJ, March to September 2020, compared with reported death data (vertical bars). The color-shaded areas depict the 95% Confidence Interval (CI) of the corresponding predictions.

Figure 9. Nonlinear relationship between the estimated county level transmission rates and county population density. The black dashed line is a fitted mechanistic nonlinear scaling function consistent with kinetic theory modeling as described in [82].

Figure 10. Comparison of projected trajectories of new/total deaths from COVID-19 calculated for different exposure scenarios with reported death data (vertical bars) for New Jersey using (a) summation of predictions from 21 county-level SEIR models, and (b) predictions from the statewide SEIR model. The color-shaded areas depict the 95% Confidence Interval (CI) of the corresponding predictions.

Table 1

COVID-19 confirmed deaths per 100,000 residents by August 31 2020 across New Jersey. The color bars in the time strip denote the weekly deaths per 100,000 people.

County/State	Total Deaths	Deaths per 100,000	[Image omitted. Please see PDF.]	Geographic Region	Total Deaths	Deaths per 100,000	[Image omitted. Please see PDF.]
New Jersey	14,068	158	[Image omitted. Please see PDF.]
Bergen	1768	189	[Image omitted. Please see PDF.]	Gateway Region	8458	197	[Image omitted. Please see PDF.]
Essex	1868	234	[Image omitted. Please see PDF.]
Hudson	1336	198	[Image omitted. Please see PDF.]
Middlesex	1218	147	[Image omitted. Please see PDF.]
Passaic	1098	218	[Image omitted. Please see PDF.]
Union	1170	210	[Image omitted. Please see PDF.]
Monmouth	763	123	[Image omitted. Please see PDF.]	Shore Region	1715	140	[Image omitted. Please see PDF.]
Ocean	952	158	[Image omitted. Please see PDF.]
Hunterdon	69	55	[Image omitted. Please see PDF.]	SkylandRegion	1559	130	[Image omitted. Please see PDF.]
Morris	684	138	[Image omitted. Please see PDF.]
Somerset	487	147	[Image omitted. Please see PDF.]
Sussex	159	113	[Image omitted. Please see PDF.]
Warren	160	151	[Image omitted. Please see PDF.]
Burlington	440	99	[Image omitted. Please see PDF.]	Delaware River Region	1866	111	[Image omitted. Please see PDF.]
Camden	541	107	[Image omitted. Please see PDF.]
Gloucester	211	72	[Image omitted. Please see PDF.]
Mercer	589	159	[Image omitted. Please see PDF.]
Salem	85	136	[Image omitted. Please see PDF.]
Atlantic	237	89	[Image omitted. Please see PDF.]	Southern Shore Region	470	92	[Image omitted. Please see PDF.]
Cape May	90	97	[Image omitted. Please see PDF.]
Cumberland	143	95	[Image omitted. Please see PDF.]

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/ijerph182211950/s1, Figure S1: heterogeneity and diversity in demographic and environmental attributes across NJ; Figure S2: NJ resident deaths, Figure S3: numbers of new confirmed COVID-19 cases at municipality level in NJ from March to September 2020, Figure S4: schematic diagram of the stochastic SEIR model with time-varying transmission rate, Figure S5: approximated treatment/hospitalization fatality rate in NJ, Figure S6: illustration of the selection process for the estimated treatment/hospitalization fatality rates, Figure S7: projected median of new/total deaths from COVID-19 for three exposure scenarios. Figure S8: spatial distribution of the median estimates for basic transmission parameter and transmission scaling factors, Figure S9: time-varying effective reproduction number, Figure S10: effectiveness of social distancing and wearing face masks with respect to population density in five categories. Figure S11: screenshots of different tabs of the dashboard for the COVID-19 project. Table S1: Biological and epidemic-relevant parameters specified with point estimates for New Jersey, Table S2: biological and epidemic-relevant parameters specified with range estimates for NJ, Table S3: statistical measures of agreement between modeled and reported confirmed COVID-19 death data, Table S4: counterfactual analysis: deaths averted per 100,000 people and peak date delayed due to social distancing and face masks use.

References

1. Ahmad, F.B.; Anderson, R.N. The leading causes of death in the US for 2020. JAMA; 2021; 325, pp. 1829-1830. [DOI: https://dx.doi.org/10.1001/jama.2021.5469] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33787821]

2. Ahmad, F.B.; Cisewski, J.A.; Miniño, A.; Anderson, R.N. Provisional mortality data—United States, 2020. Morb. Mortal. Wkly. Rep.; 2021; 70, pp. 519-522. [DOI: https://dx.doi.org/10.15585/mmwr.mm7014e1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33830988]

3. Bai, Y.; Yao, L.; Wei, T.; Tian, F.; Jin, D.-Y.; Chen, L.; Wang, M. Presumed asymptomatic carrier transmission of COVID-19. JAMA; 2020; 323, pp. 1406-1407. [DOI: https://dx.doi.org/10.1001/jama.2020.2565] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32083643]

4. Yu, P.; Zhu, J.; Zhang, Z.; Han, Y. A familial cluster of infection associated with the 2019 novel coronavirus indicating possible person-to-person transmission during the incubation period. J. Infect. Dis.; 2020; 221, pp. 1757-1761. [DOI: https://dx.doi.org/10.1093/infdis/jiaa077] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32067043]

5. Deziel, N.C.; Allen, J.G.; Scheepers, P.T.; Levy, J.I. The COVID-19 pandemic: A moment for exposure science. J. Expo. Sci. Environ. Epidemiol.; 2020; 30, pp. 591-593. [DOI: https://dx.doi.org/10.1038/s41370-020-0225-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32350383]

6. CDC. How COVID-19 Spreads. 2020; Available online: https://www.cdc.gov/coronavirus/2019-ncov/prevent-getting-sick/how-covid-spreads.html (accessed on 6 October 2020).

7. Morawska, L.; Cao, J. Airborne transmission of SARS-CoV-2: The world should face the reality. Environ. Int.; 2020; 139, 105730. [DOI: https://dx.doi.org/10.1016/j.envint.2020.105730] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32294574]

8. Morawska, L.; Milton, D.K. It is time to address airborne transmission of coronavirus disease 2019 (COVID-19). Clin. Infect. Dis.; 2020; 71, pp. 2311-2313. [DOI: https://dx.doi.org/10.1093/cid/ciaa939] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32628269]

9. Fears, A.C.; Klimstra, W.B.; Duprex, P.; Hartman, A.; Weaver, S.C.; Plante, K.S.; Mirchandani, D.; Plante, J.A.; Aguilar, P.V.; Fernández, D. Persistence of severe acute respiratory syndrome coronavirus 2 in aerosol suspensions. Emerg. Infect. Dis.; 2020; 26, 2168. [DOI: https://dx.doi.org/10.3201/eid2609.201806]

10. Fennelly, K.P. Particle sizes of infectious aerosols: Implications for infection control. Lancet Resp. Med.; 2020; 8, pp. 914-924. [DOI: https://dx.doi.org/10.1016/S2213-2600(20)30323-4]

11. Van Doremalen, N.; Bushmaker, T.; Morris, D.H.; Holbrook, M.G.; Gamble, A.; Williamson, B.N.; Tamin, A.; Harcourt, J.L.; Thornburg, N.J.; Gerber, S.I. Aerosol and surface stability of SARS-CoV-2 as compared with SARS-CoV-1. N. Engl. J. Med.; 2020; 382, pp. 1564-1567. [DOI: https://dx.doi.org/10.1056/NEJMc2004973]

12. CDC. Scientific Brief: SARS-CoV-2 and Potential Airborne Transmission. 2020; Available online: https://www.cdc.gov/coronavirus/2019-ncov/more/scientific-brief-sars-cov-2.html (accessed on 15 June 2020).

13. Kim, Y.-I.; Kim, S.-G.; Kim, S.-M.; Kim, E.-H.; Park, S.-J.; Yu, K.-M.; Chang, J.-H.; Kim, E.J.; Lee, S.; Casel, M.A.B. Infection and rapid transmission of SARS-CoV-2 in ferrets. Cell Host Microbe; 2020; 27, pp. 704-709. [DOI: https://dx.doi.org/10.1016/j.chom.2020.03.023] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32259477]

14. Richard, M.; Kok, A.; de Meulder, D.; Bestebroer, T.M.; Lamers, M.M.; Okba, N.M.; van Vlissingen, M.F.; Rockx, B.; Haagmans, B.L.; Koopmans, M.P. SARS-CoV-2 is transmitted via contact and via the air between ferrets. Nat. Commun.; 2020; 11, pp. 1-6. [DOI: https://dx.doi.org/10.1038/s41467-020-17367-2]

15. Shi, J.; Wen, Z.; Zhong, G.; Yang, H.; Wang, C.; Huang, B.; Liu, R.; He, X.; Shuai, L.; Sun, Z. Susceptibility of ferrets, cats, dogs, and other domesticated animals to SARS–coronavirus 2. Science; 2020; 368, pp. 1016-1020. [DOI: https://dx.doi.org/10.1126/science.abb7015]

16. Chia, P.Y.; Coleman, K.K.; Tan, Y.K.; Ong, S.W.X.; Gum, M.; Lau, S.K.; Lim, X.F.; Lim, A.S.; Sutjipto, S.; Lee, P.H. Detection of air and surface contamination by SARS-CoV-2 in hospital rooms of infected patients. Nat. Commun.; 2020; 11, pp. 1-7. [DOI: https://dx.doi.org/10.1038/s41467-020-16670-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32472043]

17. Lednicky, J.A.; Lauzardo, M.; Fan, Z.H.; Jutla, A.S.; Tilly, T.B.; Gangwar, M.; Usmani, M.; Shankar, S.N.; Mohamed, K.; Eiguren-Fernandez, A. Viable SARS-CoV-2 in the air of a hospital room with COVID-19 patients. MedrXiv; 2020; 4, pp. 1-20. [DOI: https://dx.doi.org/10.1016/j.ijid.2020.09.025]

18. Liu, Y.; Ning, Z.; Chen, Y.; Guo, M.; Liu, Y.; Gali, N.K.; Sun, L.; Duan, Y.; Cai, J.; Westerdahl, D. Aerodynamic analysis of SARS-CoV-2 in two Wuhan hospitals. Nat. Med.; 2020; 582, pp. 557-560. [DOI: https://dx.doi.org/10.1038/s41586-020-2271-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32340022]

19. Pei, G.; Taylor, M.; Rim, D. Human exposure to respiratory aerosols in a ventilated room: Effects of ventilation condition, emission mode, and social distancing. Sustain. Cities Soc.; 2021; 73, 103090. [DOI: https://dx.doi.org/10.1016/j.scs.2021.103090]

20. Adenaiye, O.O.; Lai, J.; de Mesquita, P.J.B.; Hong, F.H.; Youssefi, S.; German, J.R.; Tai, S.-H.S.; Albert, B.J.; Schanz, M.; Weston, S. Infectious SARS-CoV-2 in Exhaled Aerosols and Efficacy of Masks During Early Mild Infection. Clin. Infect. Dis.; 2021; [DOI: https://dx.doi.org/10.1093/cid/ciab797] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34519774]

21. Port, J.; Yinda, C.K.; Avanzato, V.; Schulz, J.; Holbrook, M.; van Doremalen, N.; Shaia, C.; Fischer, R.; Munster, V. Increased aerosol transmission for B. 1.1. 7 (alpha variant) over lineage A variant of SARS-CoV-2. bioRxiv; 2021; [DOI: https://dx.doi.org/10.1101/2021.07.26.453518]

22. Morawska, L.; Tang, J.W.; Bahnfleth, W.; Bluyssen, P.M.; Boerstra, A.; Buonanno, G.; Cao, J.; Dancer, S.; Floto, A.; Franchimon, F. How can airborne transmission of COVID-19 indoors be minimised?. Environ. Int.; 2020; 142, 105832. [DOI: https://dx.doi.org/10.1016/j.envint.2020.105832] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32521345]

23. Mahmood, A.; Eqan, M.; Pervez, S.; Alghamdi, H.A.; Tabinda, A.B.; Yasar, A.; Brindhadevi, K.; Pugazhendhi, A. COVID-19 and frequent use of hand sanitizers; human health and environmental hazards by exposure pathways. Sci. Total Environ.; 2020; 742, 140561. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2020.140561] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32623176]

24. Gostin, L.O.; Wiley, L.F. Governmental Public Health Powers During the COVID-19 Pandemic: Stay-at-home Orders, Business Closures, and Travel Restrictions. JAMA; 2020; 323, pp. 2137-2138. [DOI: https://dx.doi.org/10.1001/jama.2020.5460]

25. Lazzerini, M.; Putoto, G. COVID-19 in Italy: Momentous decisions and many uncertainties. Lancet Glob. Health; 2020; 8, pp. 641-642. [DOI: https://dx.doi.org/10.1016/S2214-109X(20)30110-8]

26. Pan, A.; Liu, L.; Wang, C.; Guo, H.; Hao, X.; Wang, Q.; Huang, J.; He, N.; Yu, H.; Lin, X. et al. Association of Public Health Interventions With the Epidemiology of the COVID-19 Outbreak in Wuhan, China. JAMA; 2020; 323, pp. 1915-1923. [DOI: https://dx.doi.org/10.1001/jama.2020.6130]

27. Horton, R. Offline: COVID-19 is not a pandemic. Lancet; 2020; 396, 874. [DOI: https://dx.doi.org/10.1016/S0140-6736(20)32000-6]

28. Bialek, S.; Boundy, E.; Bowen, V.; Chow, N.; Cohn, A.; Dowling, N.; Ellington, S.; Gierke, R.; Hall, A.; Macneil, J. et al. Severe Outcomes Among Patients with Coronavlrus Disease 2019 (COVID-19)—United States, February 12–March 16, 2020. MMWR Morb. Mortal. Wkly. Rep.; 2020; 69, pp. 343-346.

29. Onder, G.; Rezza, G.; Brusaferro, S. Case-Fatality Rate and Characteristics of Patients Dying in Relation to COVID-19 in Italy. JAMA; 2020; 323, pp. 1775-1776. [DOI: https://dx.doi.org/10.1001/jama.2020.4683]

30. Richardson, S.; Hirsch, J.S.; Narasimhan, M.; Crawford, J.M.; McGinn, T.; Davidson, K.W.; Barnaby, D.P.; Becker, L.B.; Chelico, J.D.; Cohen, S.L. Presenting characteristics, comorbidities, and outcomes among 5700 patients hospitalized with COVID-19 in the New York City area. JAMA; 2020; 323, pp. 2052-2059. [DOI: https://dx.doi.org/10.1001/jama.2020.6775]

31. Gandhi, M.; Yokoe, D.S.; Havlir, D.V. Asymptomatic Transmission, the Achilles’ Heel of Current Strategies to Control Covid-19. N. Engl. J. Med.; 2020; 382, pp. 2158-2160. [DOI: https://dx.doi.org/10.1056/NEJMe2009758]

32. Chen, T.; Wu, D.; Chen, H.; Yan, W.; Yang, D.; Chen, G.; Ma, K.; Xu, D.; Yu, H.; Wang, H. et al. Clinical characteristics of 113 deceased patients with coronavirus disease 2019: Retrospective study. BMJ; 2020; 368, pp. 1-12. [DOI: https://dx.doi.org/10.1136/bmj.m1091]

33. Koh, D. Occupational risks for COVID-19 infection. Occup. Med.; 2020; 70, pp. 3-5. [DOI: https://dx.doi.org/10.1093/occmed/kqaa036]

34. Chowkwanyun, M.; Reed, A.L., Jr. Racial health disparities and Covid-19—Caution and context. N. Engl. J. Med.; 2020; 383, pp. 201-203. [DOI: https://dx.doi.org/10.1056/NEJMp2012910]

35. Owen, W.F.; Carmona, R.; Pomeroy, C. Failing another national stress test on health disparities. JAMA; 2020; 323, 19. [DOI: https://dx.doi.org/10.1001/jama.2020.6547] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32293642]

36. Metcalf, C.J.E.; Morris, D.H.; Park, S.W. Mathematical models to guide pandemic response. Science; 2020; 369, pp. 368-369. [DOI: https://dx.doi.org/10.1126/science.abd1668] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32703861]

37. Ross, R. The Prevention of Malaria; John Murray: London, UK, 1910.

38. Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. Ser. A; 1927; 115, pp. 700-721.

39. Schneckenreither, G.; Popper, N.; Zauner, G.; Breitenecker, F. Modelling SIR-type epidemics by ODEs, PDEs, difference equations and cellular automata—A comparative study. Simul. Model. Pract. Theory; 2008; 16, pp. 1014-1023. [DOI: https://dx.doi.org/10.1016/j.simpat.2008.05.015]

40. Bacaër, N. A Short History of Mathematical Population Dynamics; Springer Science & Business Media: New York, NY, USA, 2011.

41. Chalub, F.A.; Souza, M.O. The SIR epidemic model from a PDE point of view. Math. Comput. Model.; 2011; 53, pp. 1568-1574. [DOI: https://dx.doi.org/10.1016/j.mcm.2010.05.036]

42. Zhuang, Q.; Wang, J. A spatial epidemic model with a moving boundary. Infect. Dis. Model.; 2021; 6, pp. 1046-1060. [DOI: https://dx.doi.org/10.1016/j.idm.2021.08.005]

43. CDC. COVID-19 Mathematical Modeling. Available online: https://www.cdc.gov/coronavirus/2019-ncov/covid-data/mathematical-modeling.html (accessed on 6 October 2021).

44. Giordano, G.; Blanchini, F.; Bruno, R.; Colaneri, P.; Di Filippo, A.; Di Matteo, A.; Colaneri, M. Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy. Nat. Med.; 2020; 26, pp. 855-860. [DOI: https://dx.doi.org/10.1038/s41591-020-0883-7]

45. Childs, M.L.; Kain, M.P.; Harris, M.J.; Kirk, D.; Couper, L.; Nova, N.; Delwel, I.; Ritchie, J.; Becker, A.D.; Mordecai, E.A. The impact of long-term non-pharmaceutical interventions on COVID-19 epidemic dynamics and control: The value and limitations of early models. Proc. R. Soc. B; 2021; 288, 20210811. [DOI: https://dx.doi.org/10.1098/rspb.2021.0811]

46. Bonabeau, E. Agent-based modeling: Methods and techniques for simulating human systems. Proc. Natl. Acad. Sci. USA; 2002; 99, pp. 7280-7287. [DOI: https://dx.doi.org/10.1073/pnas.082080899]

47. Grefenstette, J.J.; Brown, S.T.; Rosenfeld, R.; DePasse, J.; Stone, N.T.; Cooley, P.C.; Wheaton, W.D.; Fyshe, A.; Galloway, D.D.; Sriram, A. FRED (A Framework for Reconstructing Epidemic Dynamics): An open-source software system for modeling infectious diseases and control strategies using census-based populations. BMC Public Health; 2013; 13, pp. 1-14. [DOI: https://dx.doi.org/10.1186/1471-2458-13-940]

48. COVID, I.; Murray, C.J. Forecasting COVID-19 impact on hospital bed-days, ICU-days, ventilator-days and deaths by US state in the next 4 months. MedRxiv; 2020; pp. 1-26. Available online: https://www.medrxiv.org/content/10.1101/2020.1103.1127.20043752v20043751 (accessed on 6 October 2021).

49. Petropoulos, F.; Makridakis, S. Forecasting the novel coronavirus COVID-19. PLoS ONE; 2020; 15, e0231236. [DOI: https://dx.doi.org/10.1371/journal.pone.0231236] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32231392]

50. Altieri, N.; Barter, R.L.; Duncan, J.; Dwivedi, R.; Kumbier, K.; Li, X.; Netzorg, R.; Park, B.; Singh, C.; Tan, Y.S. Curating a COVID-19 data repository and forecasting county-level death counts in the United States. 2020; Available online: https://arxiv.org/abs/2005.07882 (accessed on 6 October 2021).

51. Gu, Y. COVID-Projections. Available online: https://covid19-projections.com/about/#historical-performance (accessed on 6 October 2021).

52. Wong, F.; Collins, J.J. Evidence that coronavirus superspreading is fat-tailed. Proc. Natl. Acad. Sci. USA; 2020; 117, pp. 29416-29418. [DOI: https://dx.doi.org/10.1073/pnas.2018490117] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33139561]

53. Allen, L.J.; Burgin, A.M. Comparison of deterministic and stochastic SIS and SIR models in discrete time. Math. Biosci.; 2000; 163, pp. 1-33. [DOI: https://dx.doi.org/10.1016/S0025-5564(99)00047-4]

54. Britton, T. Stochastic epidemic models: A survey. Math. Biosci.; 2010; 225, pp. 24-35. [DOI: https://dx.doi.org/10.1016/j.mbs.2010.01.006] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20102724]

55. West, R.W.; Thompson, J.R. Models for the simple epidemic. Math. Biosci.; 1997; 141, pp. 29-39. [DOI: https://dx.doi.org/10.1016/S0025-5564(96)00169-1]

56. Li, M.Y.; Muldowney, J.S. Global stability for the SEIR model in epidemiology. Math. Biosci.; 1995; 125, pp. 155-164. [DOI: https://dx.doi.org/10.1016/0025-5564(95)92756-5]

57. Chang, S.; Harding, N.; Cliff, O.; Prokopenko, M. Modelling transmission and control of the COVID-19 pandemic in Australia. arXiv; 2020; arXiv: 2003.10218[DOI: https://dx.doi.org/10.1038/s41467-020-19393-6]

58. Chatterjee, K.; Chatterjee, K.; Kumar, A.; Shankar, S. Healthcare impact of COVID-19 epidemic in India: A stochastic mathematical model. Med. J. Armed Forces India; 2020; 76, pp. 147-155. [DOI: https://dx.doi.org/10.1016/j.mjafi.2020.03.022] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32292232]

59. Ferguson, N.; Laydon, D.; Nedjati Gilani, G.; Imai, N.; Ainslie, K.; Baguelin, M.; Bhatia, S.; Boonyasiri, A.; Cucunuba Perez, Z.; Cuomo-Dannenburg, G. Report 9: Impact of non-pharmaceutical interventions (NPIs) to reduce COVID19 mortality and healthcare demand. Imp. Coll. Lond.; 2020; pp. 1-20.

60. Kucharski, A.J.; Russell, T.W.; Diamond, C.; Liu, Y.; Edmunds, J.; Funk, S.; Eggo, R.M.; Sun, F.; Jit, M.; Munday, J.D. et al. Early dynamics of transmission and control of COVID-19: A mathematical modelling study. Lancet Infect. Dis.; 2020; 20, pp. 553-558. [DOI: https://dx.doi.org/10.1016/S1473-3099(20)30144-4]

61. Matrajt, L.; Leung, T. Evaluating the Effectiveness of Social Distancing Interventions to Delay or Flatten the Epidemic Curve of Coronavirus Disease. Emerg. Infect. Dis.; 2020; 26, pp. 1740-1748. [DOI: https://dx.doi.org/10.3201/eid2608.201093] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32343222]

62. Kain, M.P.; Childs, M.L.; Becker, A.D.; Mordecai, E.A. Chopping the tail: How preventing superspreading can help to maintain COVID-19 control. Epidemics; 2021; 34, 100430. [DOI: https://dx.doi.org/10.1016/j.epidem.2020.100430] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33360871]

63. Thomas, L.J.; Huang, P.; Yin, F.; Luo, X.I.; Almquist, Z.W.; Hipp, J.R.; Butts, C.T. Spatial Heterogeneity Can Lead to Substantial Local Variations in COVID-19 Timing and Severity. arXiv; 2020; pp. 24180-24187. arXiv: 2005.09850[DOI: https://dx.doi.org/10.1073/pnas.2011656117] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32913057]

64. Wang, H.; Yamamoto, N. Using a partial differential equation with Google Mobility data to predict COVID-19 in Arizona. arXiv; 2020; arXiv: 2006.16928[DOI: https://dx.doi.org/10.3934/mbe.2020266]

65. Bertaglia, G.; Boscheri, W.; Dimarco, G.; Pareschi, L. Spatial spread of COVID-19 outbreak in Italy using multiscale kinetic transport equations with uncertainty. arXiv; 2021; arXiv: 2106.07262[DOI: https://dx.doi.org/10.3934/mbe.2021350]

66. Bertaglia, G.; Pareschi, L. Hyperbolic compartmental models for epidemic spread on networks with uncertain data: Application to the emergence of Covid-19 in Italy. arXiv; 2021; arXiv: 2105.14258[DOI: https://dx.doi.org/10.1142/S0218202521500548]

67. Viguerie, A.; Lorenzo, G.; Auricchio, F.; Baroli, D.; Hughes, T.J.; Patton, A.; Reali, A.; Yankeelov, T.E.; Veneziani, A. Simulating the spread of COVID-19 via a spatially-resolved susceptible–exposed–infected–recovered–deceased (SEIRD) model with heterogeneous diffusion. Appl. Math. Lett.; 2021; 111, 106617. [DOI: https://dx.doi.org/10.1016/j.aml.2020.106617]

68. Viguerie, A.; Veneziani, A.; Lorenzo, G.; Baroli, D.; Aretz-Nellesen, N.; Patton, A.; Yankeelov, T.E.; Reali, A.; Hughes, T.J.; Auricchio, F. Diffusion–reaction compartmental models formulated in a continuum mechanics framework: Application to COVID-19, mathematical analysis, and numerical study. Comput. Mech.; 2020; 66, pp. 1131-1152. [DOI: https://dx.doi.org/10.1007/s00466-020-01888-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32836602]

69. JHU (Johns Hopkins University). Death rates from coronavirus (COVID-19) in the United States as of November 6, 2020, by state. Available online: https://www.statista.com/statistics/1109011/coronavirus-covid19-death-rates-us-by-state/ (accessed on 6 October 2021).

70. Caflisch, R.E. Monte carlo and quasi-monte carlo methods. Acta Numer.; 1998; 1998, pp. 1-49. [DOI: https://dx.doi.org/10.1017/S0962492900002804]

71. Chopin, N.; Papaspiliopoulos, O. An Introduction to Sequential Monte Carlo; Springer: Berlin/Heidelberg, Germany, 2020.

72. Gerber, M.; Chopin, N. Sequential quasi-Monte Carlo. J. R. Statist. Soc. B; 2015; 77, pp. 509-579. [DOI: https://dx.doi.org/10.1111/rssb.12104]

73. Lemieux, C. Monte Carlo and Quasi-Monte Carlo Sampling; Springer Science & Business Media: New York, NY, USA, 2009.

74. Endo, A.; van Leeuwen, E.; Baguelin, M.J.E. Introduction to particle Markov-chain Monte Carlo for disease dynamics modellers. Epidemics; 2019; 29, 100363. [DOI: https://dx.doi.org/10.1016/j.epidem.2019.100363]

75. King, A.A.; Nguyen, D.; Ionides, E.L. Statistical inference for partially observed Markov processes via the R package pomp. 2015; Available online: https://arxiv.org/abs/1509.00503 (accessed on 6 October 2021).

76. CEBM (Center for Evidence-Based Medicine). Declining death rate from COVID-19 in hospitals in England. Available online: https://www.cebm.net/covid-19/declining-death-rate-from-covid-19-in-hospitals-in-england/ (accessed on 6 October 2021).

77. NJDOH. Has the State Used Data to Make Decisions and Slow the Spread of COVID-19?. Available online: https://covid19.nj.gov/faqs/nj-information/reopening-guidance-and-restrictions/has-the-state-used-data-to-make-decisions-and-slow-the-spread-of-covid-19 (accessed on 6 October 2021).

78. Safarishahrbijari, A.; Teyhouee, A.; Waldner, C.; Liu, J.; Osgood, N.D. Predictive accuracy of particle filtering in dynamic models supporting outbreak projections. BMC Infect. Dis.; 2017; 17, pp. 1-12. [DOI: https://dx.doi.org/10.1186/s12879-017-2726-9]

79. Ionides, E.L.; Nguyen, D.; Atchadé, Y.; Stoev, S.; King, A.A. Inference for dynamic and latent variable models via iterated, perturbed Bayes maps. Proc. Natl. Acad. Sci. USA; 2015; 112, pp. 719-724. [DOI: https://dx.doi.org/10.1073/pnas.1410597112]

80. Sy, K.T.L.; White, L.F.; Nichols, B.E. Population density and basic reproductive number of COVID-19 across United States counties. medRxiv; 2020; pp. 1-17. Available online: https://www.medrxiv.org/content/10.1101/2020.1106.1112.20130021v20130021 (accessed on 6 October 2021).

81. Anderson, R.M.; Anderson, B.; May, R.M. Infectious Diseases of Humans: Dynamics and Control; Oxford University Press: 1992; Available online: https://aasldpubs.onlinelibrary.wiley.com/doi/abs/10.1002/hep.1840150131 (accessed on 6 October 2021).

82. Hu, H.; Nigmatulina, K.; Eckhoff, P. The scaling of contact rates with population density for the infectious disease models. Math. Biosci.; 2013; 244, pp. 125-134. [DOI: https://dx.doi.org/10.1016/j.mbs.2013.04.013]

83. Grave, M.; Viguerie, A.; Barros, G.F.; Reali, A.; Coutinho, A.L. Assessing the spatio-temporal spread of COVID-19 via compartmental models with diffusion in Italy, USA, and Brazil. arXiv; 2021; arXiv: 2102.07208[DOI: https://dx.doi.org/10.1007/s11831-021-09627-1]

84. Delamater, P.L.; Street, E.J.; Leslie, T.F.; Yang, Y.T.; Jacobsen, K.H. Complexity of the basic reproduction number (R0). Emerg. Infect. Dis.; 2019; 25, pp. 1-4. [DOI: https://dx.doi.org/10.3201/eid2501.171901]

85. Liu, Y.; Gayle, A.A.; Wilder-Smith, A.; Rocklöv, J. The reproductive number of COVID-19 is higher compared to SARS coronavirus. J. Travel Med.; 2020; 27, pp. 1-4. [DOI: https://dx.doi.org/10.1093/jtm/taaa021] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32052846]

86. Sanche, S.; Lin, Y.T.; Xu, C.; Romero-Severson, E.; Hengartner, N.; Ke, R. High Contagiousness and Rapid Spread of Severe Acute Respiratory Syndrome Coronavirus 2. Emerg. Infect. Dis.; 2020; 26, pp. 1470-1477. [DOI: https://dx.doi.org/10.3201/eid2607.200282] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32255761]

87. Yuan, J.; Li, M.; Lv, G.; Lu, Z.K. Monitoring Transmissibility and Mortality of COVID-19 in Europe. Int. J. Infect. Dis.; 2020; 95, pp. 311-315. [DOI: https://dx.doi.org/10.1016/j.ijid.2020.03.050] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32234343]

88. Ridenhour, B.; Kowalik, J.M.; Shay, D.K. Unraveling R0: Considerations for public health applications. Am. J. Public Health; 2018; 108, pp. S445-S454. [DOI: https://dx.doi.org/10.2105/AJPH.2013.301704r]

89. Nishiura, H.; Chowell, G. The effective reproduction number as a prelude to statistical estimation of time-dependent epidemic trends. Mathematical and Statistical Estimation Approaches in Epidemiology; Springer: Berlin/Heidelberg, Germany, 2009; pp. 103-121.

90. NJDOH. New Jersey COVID-19 Dashboard. Available online: https://www.nj.gov/health/cd/topics/covid2019_dashboard.shtml (accessed on 6 October 2021).

91. Obadia, T.; Haneef, R.; Boëlle, P.-Y. The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks. BMC Med. Inform. Decis. Mak.; 2012; 12, pp. 1-9. [DOI: https://dx.doi.org/10.1186/1472-6947-12-147]

92. rt.live. New Jersey Effective Reproduction Rate (Rt). Available online: https://rt.live/us/NJ (accessed on 6 October 2021).

93. Zhang, R.; Li, Y.; Zhang, A.L.; Wang, Y.; Molina, M. Identifying airborne transmission as the dominant route for the spread of COVID-19. Proc. Natl. Acad. Sci. USA; 2020; 117, pp. 14857-14863. [DOI: https://dx.doi.org/10.1073/pnas.2009637117]

94. Fisher, K.A.; Barile, J.P.; Guerin, R.J.; Esschert, K.L.V.; Jeffers, A.; Tian, L.H.; Garcia-Williams, A.; Gurbaxani, B.; Thompson, W.W.; Prue, C.E. Factors associated with cloth face covering use among adults during the COVID-19 pandemic—United States, April and May 2020. MMWR Morb. Mortal. Wkly. Rep.; 2020; 69, pp. 933-937. [DOI: https://dx.doi.org/10.15585/mmwr.mm6928e3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32673303]

95. Davies, N.G.; Klepac, P.; Liu, Y.; Prem, K.; Jit, M.; Eggo, R.M. CMMID COVID-19 working group. Age-dependent effects in the transmission and control of COVID-19 epidemics. Nat. Med.; 2020; 26, pp. 1205-1211. [DOI: https://dx.doi.org/10.1038/s41591-020-0962-9]

96. Kraemer, M.U.; Yang, C.-H.; Gutierrez, B.; Wu, C.-H.; Klein, B.; Pigott, D.M. Open COVID-19 Data Working Group du Plessis, L.; Faria, N.R.; Li, R. The effect of human mobility and control measures on the COVID-19 epidemic in China. Science; 2020; 368, pp. 493-497. [DOI: https://dx.doi.org/10.1126/science.abb4218]

97. Pei, S.; Kandula, S.; Shaman, J. Differential effects of intervention timing on COVID-19 spread in the United States. Sci. Adv.; 2020; 6, eabd6370. [DOI: https://dx.doi.org/10.1126/sciadv.abd6370] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33158911]

98. Naughton, S.X.; Raval, U.; Harary, J.M.; Pasinetti, G.M. The role of the exposome in promoting resilience or susceptibility after SARS-CoV-2 infection. J. Expo. Sci. Environ. Epidemiol.; 2020; 30, pp. 776-777. [DOI: https://dx.doi.org/10.1038/s41370-020-0232-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32398780]

99. Georgopoulos, P.G.; Ren, X.; Mi, Z. Socioexposomics of COVID-19: The Case of New Jersey. Proceedings of the International Society of Exposure Science (ISES) Conference; Virtual Conference, 30 August–2 September 2021.

100. Dashtbali, M.; Mirzaie, M. A compartmental model that predicts the effect of social distancing and vaccination on controlling COVID-19. Sci. Rep.; 2021; 11, pp. 1-11. [DOI: https://dx.doi.org/10.1038/s41598-021-86873-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33854079]

101. Moghadas, S.M.; Vilches, T.N.; Zhang, K.; Wells, C.R.; Shoukat, A.; Singer, B.H.; Meyers, L.A.; Neuzil, K.M.; Langley, J.M.; Fitzpatrick, M.C. The impact of vaccination on COVID-19 outbreaks in the United States. medRxiv; 2021; [DOI: https://dx.doi.org/10.1101/2020.11.27.20240051]

102. Rachaniotis, N.P.; Dasaklis, T.K.; Fotopoulos, F.; Tinios, P. A two-phase stochastic dynamic model for COVID-19 mid-term policy recommendations in Greece: A pathway towards mass vaccination. Int. J. Environ. Res. Public Health; 2021; 18, 2497. [DOI: https://dx.doi.org/10.3390/ijerph18052497] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33802501]