key: cord-0551838-9y1wrqtf
authors: Gopalakrishnan, Vishrawas; Navalekar, Sayali; Ding, Pan; Hooley, Ryan; Miller, Jacob; Srinivasan, Raman; Deshpande, Ajay; Liu, Xuan; Bianco, Simone; Kaufman, James H.
title: Adaptive Epidemic Forecasting and Community Risk Evaluation of COVID-19
date: 2021-06-03
journal: nan
DOI: nan
sha: 22c2380cce89849a7ff85f9845f25593f0773954
doc_id: 551838
cord_uid: 9y1wrqtf

Pandemic control measures like lock-down, restrictions on restaurants and gatherings, social-distancing have shown to be effective in curtailing the spread of COVID-19. However, their sustained enforcement has negative economic effects. To craft strategies and policies that reduce the hardship on the people and the economy while being effective against the pandemic, authorities need to understand the disease dynamics at the right geo-spatial granularity. Considering factors like the hospitals' ability to handle the fluctuating demands, evaluating various reopening scenarios, and accurate forecasting of cases are vital to decision making. Towards this end, we present a flexible end-to-end solution that seamlessly integrates public health data with tertiary client data to accurately estimate the risk of reopening a community. At its core lies a state-of-the-art prediction model that auto-captures changing trends in transmission and mobility. Benchmarking against various published baselines confirm the superiority of our forecasting algorithm. Combined with the ability to extend to multiple client-specific requirements and perform deductive reasoning through counter-factual analysis, this solution provides actionable insights to multiple client domains ranging from government to educational institutions, hospitals, and commercial establishments.

Due to the COVID-19 pandemic, it is estimated that half of the world's 3.3 billion global workforces were/are at risk of losing their livelihoods [11] and in the U.S itself, there has been a -11% change in the number of low-wage category jobs, and overall a -6% decrease across categories. Similar disruptions are reported across different work sectors with AHA estimating losses of $202.6 billion for America's hospitals and health systems for the period of March-1 2020 to June-30 2020 [18] and over 1,300 colleges and universities across all 50 states canceled in-person classes affecting the quality of education as well as leading to financial losses [25] . To alleviate these issues, studies have recommended that governments and industries should, as soon as possible, try to resume operations even in a limited capacity and then adapt their operations based on the evolving situation [24] .

Towards this goal, our solution -part of IBM's Watson Works, provides timely information on the pandemic's current and future state in a region. As a solution framework, with an objective to cater to different specific client domains, it provides actionable insights that support planning safe reopening of offices, crafting school schedules, meeting hospitalization/ICU demands, planning for clinical vaccine trials, etc. Additionally, it has also helped local officials and state government test hypotheses and crafting their intervention strategies for the COVID-19 pandemic. Several efforts like [4, 6, 16, 27, 29] provide forecasts of COVID-19. However, the superiority of our method lies in its ability to provide a solution to:

• Model at a hyper-local level.

• Translate the current and predicted case load numbers to community reopening risk metrics in accordance with national/regional laws. • Easily adapt to business specific tasks and objectives like hospital/ICU demand projections or identifying appropriate location for vaccine trials. • Scale globally through distributed processing.

• Interact and simulate what-if scenarios.

• Incorporate tertiary signals apart from cases and deaths data to better estimate the future trends. As shown in Section 4, our system consistently performs better in both short and long-term projections, enabling users to take suitable corrective steps. The "base" prediction model is easily enhanced to handle user-specific tasks by plugging in custom analytics. For instance, this solution helped a major hospital in south Florida to accurately estimate hospital bed and ICU demands; thus, allowing them to plan/reschedule elective procedures and surgeries, which otherwise could have got cancelled. It also helped a pharmaceutical company to identify countries for their vaccine's clinical trials. Apart from customizing some of the known methodologies for data orchestration, the deployed solution also makes a significant contribution to core epidemiology by providing a template for models to incorporate important tertiary signals without increasing the parameter search space. Further, it also demonstrates a novel way of allowing parameters (e.g., transmission rate) to assume multiple values over the course of the pandemic to better reflect changes in the pandemic, including interventions, policy enactments, and their time-varying efficacy. We also report learnings and insights regarding variation in the relative importance of different input signals over time, the impact of geographic level (county vs. state) on the prediction, etc. The solution begins by first connecting to 'The Weather Channel' (TWC) API and consuming daily published COVID-19 cases, deaths, and other pertinent statistics. TWC scrapes data from individual federal/state/county government websites and stores data, including individual cases, test positivity rate, hospitalization data, etc. After performing audits on data quality, TWC pushes data into their database for downstream consumption.

The next step in the pipeline is Data Curation wherein raw COVID-19 cases and death numbers are pre-processed and merged with tertiary data, including mobility and other client-specific data such as individual hospital and ICU demand data. The latter is optional as the Client Specific Analysis module can be disabled. As shown in Figure 1 , pre-processing steps include adjusting spatial granularity, handling reporting irregularities, data smoothening, and inflection point detection.

The pre-processed data is then passed to the Prediction Module which provides case and death projections with upper and lower interval bounds. The current case and death data, along with the projections from the model, are provided to the Client Independent Analytics module to calculate various epidemiological metrics and derive a score between 1-6 (lower is better) reflecting the community transmission risk in a region. At the same time, the model prediction output along with supplemental client data drives Client Specific Analytics insights like hospitalization/ICU demand or simulating the effect of an event (holiday travel, sports event) or intervention (government restrictions) on the pandemic's transmission rate.

The entire process is wrapped in an Apache NiFi data flow pipeline that invokes relevant modules and manages the data flow. The pipeline is triggered once a day after the latest COVID-19 data is populated by TWC. The 'Prediction Module' is triggered once every three days as the projections rarely deviate a lot in such a short period. However, the 'Client Independent Analytics', which includes computing community risk, is triggered every day to provide up-to-date statistics like doubling period, 14-day case trends, 1-week, and 2-week ahead community risk levels, etc. The 'Client Specific Analytics' is run as per client requirements.

The following sections focus on a pipeline configured for U.S.A data, although the solution itself is applicable globally with relevant changes to account for the availability of tertiary data.

This section details various algorithms and steps used in the different modules shown in Figure 1 by first introducing the motivation/source of the problem, the solution intuition, and lastly, by providing specific operational examples.

This step entails two important tasks: (i) Denoising the COVID-19 data and (ii) Identifying approximate timestamps of the inflection points in disease dynamics.

Geo-Spatial Noise: Any pandemic is defined by its transmission rate -the average rate at which individuals pass on the infection. This, in turn, is determined by the population contact or mixing rate in the region. Assuming a similar contact rate across regions will lead to an incorrect susceptible population. For instance, counties of New York City-Newark-Jersey City, NY-NJ-PA Metropolitan Statistical Area have a higher mixing rate than rural counties of Alaska or North Dakota, which are sparsely populated and selfcontained. Thus, the first step is to normalize regions/counties Additionally, it is often found that in an area where there is a major hospital center, reported cases around that area are skewed as patients from nearby counties are reported at the major hospital. For example, COVID cases that occurred in Bronx County, N.Y. were reported as originating in Westchester County, N.Y. Towards this end, our solution uses Louvain's method of graph clustering [3] to group adjacent counties together. The algorithm uses the census bureau's residence to workplace commute data 2011-2015 1 to create an adjacency matrix where each node is a county and edges represent the strength of commute. Note that we symmetrize the adjacency matrix to reflect the return journey as well. Using this approach, we were able to group 3,132 counties into super-groups of 685 county clusters. There are about 4.59 counties per cluster, and 19% of the total clusters are singletons. The intersection of these county clusters with statistical regions like Metropolitan Statistical Areas is high. Further, the clustering method allows us to cover the entire USA as opposed to the limited commercially active zones covered by the census bureau's statistical area maps. Additionally, it provides the flexibility to impose state border conditions on the adjacency matrix; thus, supporting the roll-up of statistics to a state-level when required. Figure 2 shows the county cluster maps of Alaska and New York.

Reporting Noise: Often, public health data sources avoid making retrospective changes to correct reporting errors. As an example, Figure 3 shows that the cumulative incidence data for Bristol County, RI, (FIPS code 44001) is not monotonically increasing. The reason for such aberrations can be numerous, ranging from sources incorrectly assigning new cases to a hospital (instead of the patient's home address) to changes in the way cases are counted. To resolve these issues, we employ Isotonic regression [5] defined by:

Isotonic regression involves finding a weighted least-squares fit ∈ R to a vector ∈ R with weights vector ∈ R subject to a set of non-contradictory constraints of the kind ≤ . Figure 3 also shows the fit of Isotonic regression on the noisy data.

While the Isotonic regression ensures the cumulative number of cases is monotonically increasing, it does not safeguard against sudden spikes that arise periodically due to a backlog in lab tests 1 www.census.gov/data/tables/2015/demo/metro-micro/commuting-flows-2015.html reported (e.g., on the weekend) or correction in previously underreported numbers. Since we model both daily and cumulative incidences, these fluctuations would affect the model performance.

Towards that end, we use an adaptive-degree polynomial filter (ADPF) [2] to smoothen the daily cases. It is shown that that ADPF performs nearly as well as the optimally chosen fixed-degree Savitzky-Golay filter and outperforms sub-optimally chosen Savitzky-Golay filters [2] . Note that we only calculate cumulative deaths as the daily death numbers are small, and stochasticity makes daily optimization problematic for deaths.

As mentioned earlier, transmission rate defines the dynamics of the pandemic, and this varies over time. The change in the transmission rate (and other disease parameters) primarily stems from Non-Pharmaceutical Interventions (NPIs) and changes in human behavior. NPIs include social distancing practices, stay-at-home orders, night-time movement restrictions, restricted public gatherings, etc. Unlike pharmaceutical interventions like vaccines that reduce the susceptible population by moving individuals from to a 'removed' compartment (e.g., ), NPIs aim to lower the contact rate between infectious and susceptible individuals.

While there are efforts to capture the impact and effectiveness of NPIs in curtailing COVID-19 cases [22, 23, 26] , there are challenges in using them in a live system. Notably, there are region-specific factors like long weekends or festivities, which affect the disease parameters that are not usually part of the database. There is also an overhead of maintaining an up-to-date list of all interventions. Furthermore, transmission rate changes may happen due to human behavior as well. Lastly, there is often a time delay between the implementation (or relaxation) of an NPI and its effect on the observed cases and deaths. This time delay may be unknown until changes in the cases and deaths are observed. Our decision to auto-detect changes in the transmission rate stems from this last observation. Namely, the entire history of the pandemic in a region is split into multiple timepieces, where each piece represents the observed variable of parameters in that time range. Note that in this work, we focus only on transmission rate changes due to NPIs as there is not yet enough data on vaccination.

To identify changes in model parameters, we use Kneedle algorithm [21] on the COVID-19 case data to detect inflection points. The algorithm is run on smoothened daily case data. If, despite smoothening, the algorithm detects multiple inflection points in close vicinity (within 7 days), they are merged. We also take into consideration the 'significance' of each inflection point, wherein we not only measure the change in the slope of the curve around the inflection point but also ensure that if a 'knee' is detected, the subsequent values over two weeks are decreasing (vice-versa for an 'elbow'). Figure 4 illustrates the auto-detected intervention for county-level data using our algorithm and superimposes the local government imposed NPI.

These inflection points provide an approximate suggestion to our model (Section 3.2.3) regarding times where it might have to (re)adjust parameters to get an optimal fit. Note that this approach is just one of the two parts to dynamically adapt time sensitive parameters like transmission rate. This part handles scenarios where the effect of an intervention (or relaxation) is already visible. The scenario where the effect is anticipated is handled by incorporating mobility as a tertiary source (refer Section 3.2.2).

We use a compartment structure -a generic framework of models extensively used in epidemiology [1, 10, 13, 28] , to model COVID-19 cases and deaths. A compartment model includes the transition/flow of people across different stages of a disease outbreak. At the minimum, a compartment flow model contains two states -representing a group of susceptible people in a population and representing infectious people. Depending on the type of micro-parasitic infection, epidemiologists add compartments based on the organism's life cycle, patient state, and transmission mode. It is worth noting that the ecosystem can span beyond humans; for instance, in the case of Avian Influenza there is an interplay between humans and birds (and mosquitoes) with different sets of parameters. In the case of COVID-19, we would need at minimum two additional compartments to reflect people who recovered ( ) and died ( ). It is to be noted that the compartment should be thought to be as 'Removed' compartment; i.e., if the compartment is initialized with say 20% of the population in the region, it would then amount to saying that 20% of the people would never get affected -indirectly accounting for a concept called 'herd immunity'. So depending on the biological property of the virus, the 'herd immunity' ratio is set -for COVID-19, we initialize with 20% of the total population [7] . In our model, we add two additional 'infectious' compartments to account for Asymptomats, including unreported infectious people, ( ) and Pre-symptomats ( ). Like the normal flu, it is reported [12] that many people infected by COVID-19 do not exhibit any symptoms. However, they continue to shed viruses and contribute to the force of infection. Conversely, some infected people with mild symptoms may not visit a doctor to be counted in the COVID-19 infected list. The compartment represents these patients.

There is an incubation time for an individual to begin exhibiting the symptoms. Even during the time when they are 'pre-symptomatic', such people may shed the Sars-CoV-2 virus and contribute to the force of infection. It is necessary to account for this 'pre-symptomatic' period in our model. It is important to note that someone being marked COVID-19 positive corresponds to the transition from → . The time it takes for a person to move from → corresponds to symptom appearance rate ( ). Normally that represents just the incubation time, but in this case, it also accounts for latency in testing and reporting. Figure 5 represents our compartment model. The compartment holds the susceptible population after accounting for the percentage of 'removed' population used to initialize the compartment. From the compartment, an infected person can be either in or in , which is determined by two things: (a) Likelihood of the person getting infected -represented by the transmission rate ( ), and (b) Likelihood that patient ends up being tested if not showing any symptoms -represented by the case reporting rate ( ). As shown in Figure 5 , we also introduced a compartment , which does not contribute to the force of infection. The compartment accounts for patients who become worse and require extra attention and possibly hospitalization. We assume the people in are isolated in hospital and hence do not contribute to the force of infection. As this base model uses only reported incidence and death data for fitting, acts as a delay compartment to account for the time delay ( ) between infection and death. People can recover from , , ; albeit with different rates , , . People enter the death compartment at a rate of . COVID-19 being similar to influenza, also has an 'immunity loss rate' ( ), denoting the rate at which a person previously infected may lose their immunity and become susceptible again. We set it to 1/(10 months) [14] .

Mobility: Figure 5 also shows an isolation compartment ( ). The and compartments have a cyclic connection. The idea is to stage as a transient state to simulate people who are not mixing with the population. No restriction on population mixing is an assumption made by standard compartment models. However, with the stay-at-home order and decreased mobility/interaction, this assumption is violated. Continuing that assumption would lead to faster depletion of the susceptible population and overestimation of the asymptomats (note that we fit only with respect to cases and deaths, and so the compartment is indirectly determined by the testing/case reporting rate ). To account for the population 'hidden' from virus exposure and scenarios where effects of relaxation of an intervention are yet to reflect on observed cases (recalling discussions from Section 3.1.2), we set the transition rate from ↔ to be equal to rate of change in the daily mobility. This value is obtained from tertiary sources like Apple 2 , Google 3 , or Descartes Lab 4 . If the normalized difference between two days is positive, it means the mobility is on the rise and the flow from → is triggered, else from → . Thus, at any point in time, the transition between and is always unidirectional.

The advantage of this model is that we are able to include mobility data without any additional tunable parameter(s). Note that in theory, the isolation compartment can be attached to and with exits to and , respectively. However, this will lead to overcomplication of the model and an additional parameter → .

Here, we provide the mathematical description of the compartmental flow illustrated in Figure 5 . These differential equations represent changes in the state of population (each compartment) at any time instance. As the time scale is continuous, we need to add interpolation functions on any of our tertiary sources that directly control the flow -e.g., changes in the mobility that controls ↔ . Let be a function fitted on the mobility, then the equation from ↔ can be described as:

The fitted function can also be used to extrapolate mobility data to get likely mobility values in the future. Now, the compartmental flow rates can be defined as: 

Note that the in the above equations are time-varying parameters as explained in Section 3.1.2 -i.e., there will be one for each of the inflection-detected time intervals, and selecting which beta to modify is based on the time stamp of the ordinary differential equation (ODE) solver. Recall that the approach in Section 3.1.2 only provides an initial estimate for the time of the inflection. Since the ODE solver iterates over time steps, it can find the 'optimal' from the rough estimate to fit the piece-wise curve. Through this, we are able to mimic the natural evolution of the disease, and unlike a regression algorithm, this does not violate any epidemiological principle. We also set ( → , ), ( → ), ( → ) to be a time-sensitive parameter similar to .

We use an LSODA algorithm [8, 9, 20] to solve the ODEs. LSODA automatically selects between non-stiff (Adams) and stiff (BDF) methods. It uses the non-stiff method initially and dynamically monitors data to decide which method to use. The parameter estimation is done by Levenberg-Marquardt algorithm [15, 17] ; selected mainly for its simplicity and speed. The error function to minimize includes daily case NRMSE, cumulative case NRMSE, and cumulative death NRMSE. Because the entire system of equations receives continuous external 'shocks' via the mobility values, the solver takes additional time to converge. However, as discussed in Section 4.1.2, including this external data provides valuable information and improves performance.

Our solution provides a set of analytics regardless of which client use-case it is serving. These typically include epidemiological metrics like 0 , doubling time, 14-day rolling average, etc. Amongst these, 0 is one of great importance. It provides a likelihood of how much an infected individual is likely to infect others. It is directly connected to the transmission rate ( ) and recovery time ( ). However, it is also model specific. Solving Equations 6 to Equation 11 , under steady-state, we derive 0 :

where the , , and in the above equation refer to the first fitted value of , , and . Recall that we fit multiple mimicking the different transmission rates over the course of the pandemic in a region. If we calculate at each of those different times, then using the fraction of susceptibles during each of those times, we get −effective ( ). is used to denote the latest . Together with doubling time (Equation 13), these two statistics provide an accurate sense of the pandemic situation on the ground. 

3.3.1 Community Risk Evaluation. The community risk evaluation uses a rule-based approach to convert the epidemiological metrics into risk scores for a region between 1-6, with 1 being the safest. The features that determine the score include recent incidence trends, the case prevalence in the population, and local government 3 ), , , 2: Output: Current and 3 weeks projected community risk scores 3: is_strict_decr = 3 < 2 ∧ 2 < 1 4: is_strict_incr = 3 > 2 ∧ 2 > 1 5: if 3 , 2 , 1 < then 6: if 3 , 2 < then 

As mentioned above, our base model could be easily extended for client-specific use cases. In one of the client engagements, we augmented the base model with hospitalization data to predict both hospitalization and ICU demand. Figure 6 shows the model extension. The input to this extension is the predicted incidence from our base model. Note that we do not use the predicted deaths from the base model as-is, because in this case, we need to infer the mortality rate from the ICU training data. The hospitalization, ICU, and death rates are learnt using the same infrastructure used in solving the base model. The equations can be derived similarly to our base model described in Section 3.2.3. As shown in Section 4.2.3, extending to simulate what-if scenarios are also straightforward, as one can adjust the learnt parameters and run the model to solve the ODEs in a 'scoring' mode. 

Having described the relevant algorithms and modules, we now present a detailed evaluation analysis of the results of the various modules in the system. Our experiments can be broken down into two broad categories wherein we first measure the prediction performance of our base/core system, and then through various user story scenarios, explain the performance, usefulness, and dexterity of the various analytic modules.

Base Module Evaluation Scheme: We use Mean Absolute Percentage Errors (MAPE) to quantitatively measure the prediction performance. We compare our results with the baselines presented in CDC 5 . To avoid bias in measuring performance, we measure the accuracy at both the falling phase (end of first wave~Aug-3 ) and the rising curve (start of the second wave~Oct-26 ℎ ). Note that CDC did not report case prediction before Aug-1 . We also perform a time series ablation study to measure the importance of mobility data and then wrap up the base module analysis by presenting geo-spatial and running time analysis.

Analytic Module Evaluation Scheme: We consider real use-cases to characterize the performance of the 'Client Independent Analytics' and 'Client Specific Analytics' modules. In one of the client engagements pertaining to the state of Rhode-Island, we highlight how the compartment model helps track the disease progression in the region and compare the validity of the estimated community risk scores with the situation on the ground. Through our support for a major hospital in Tampa, FL, we evaluate the hospital and ICU demand prediction module. As part of our engagement, we also simulate what-if scenarios in response to the Thanksgiving and holiday season, enumerating the effectiveness of various counter-measures. Table 1 shows the cumulative weekly incidence MAPEs calculated from the end of training dates (August-3 and October-26 ℎ ). As we can see, the 'IBM' system consistently comes in top-3, and for the week-1 and week-1+2, it is at the top. The numbers in the bracket indicate the MAPE, and we can see that for the week-1, the difference is substantial. Note that for October, the rest of the top-5 baselines do not make it into top-5 for the remaining weeks. It is also to be noted that some of the methodologies like Oliver Wyman started publishing their results in September and hence do not appear in the August comparison. Table 2 shows the cumulative deaths MAPE. In the interest of space, we report only the October-26 ℎ value. We observe a similar trend like Table 1 where our model consistently appears in top-2.

We also noted that some reported methodologies predict only cases or deaths (e.g., YYG). Based on Table 1 and 2, IBM and Oliver Wyman are comparable models with IBM doing particularly well in the 1 week. There are two important points to remember while interpreting these tables. First, death numbers are substantially smaller than incidence, so the sensitivity of MAPE to incorrect predictions is higher. Second, as the MAPE is calculated based on weekly values, a small deviation in the reported ground truth numbers at the end of an evaluation week may adversely impact the MAPE calculation itself.

Week-1

Week-1+2

Week-1+2+3

Week-1

Week-1+2 Table 3 compares the MAPE of the model with and without mobility for multiple runs. While there is an improvement across the board, we note a 11% improvement in the pandemic's early days. Figure 7 shows how the mobility data is able to adjust the trajectory more accurately. While some studies have already suggested that mobility is a leading indicator for COVID-19 cases [19] , in this work, we further show how to incorporate that information in practice and its usefulness in not only improving the accuracy but also understanding the population mixing and the susceptible pool. We also observed that the contribution of mobility data towards prediction has waned over time. We suspect this may be due to increased awareness, mask adoption, and increasing herd immunity (i.e., as the population develops antibodies from infection or vaccination). Thus, it is our learning that using mobility data, if available, is a must to get accurate predictions during the early stages of a pandemic. It provides a more precise estimate of susceptibles and the magnitude of the potential second wave.

The choice of geo-unit for analysis is important. As shown in Table 4 , modeling at a state level as opposed to aggregating the projections from county-cluster leads to sub-par MAPE. This is because the local intricacies of population mixing are lost when the entire state is considered a well-mixed unit. Thus, crafting an accurate policy at a higher geo-level requires hyper-local modeling. A similar conclusion also holds true when using the state projections to estimate county clusters using area population to drill down. In general, we recommend modeling at county-clusters granularity during the pandemic's initial stages to avoid data sparsity issues. Once the counties within the cluster enter the exponential growth phase (vs. sporadic or intermittent case reports), modeling at the county level will prove beneficial, albeit with a trade-off on the running time.

The running time of the pipeline is determined by the choice of parameter optimization algorithm and the external "shocks" introduced through tertiary sources like mobility data. These shocks cause a sudden change in the equilibrium of the system. Furthermore, with the increase in the number of timestamps (and consequently the number of time-sensitive parameters), the effect is non-linear. Table 5 shows the running time of finding the best parameters for 20 initializer combinations. The time-sensitive parameters also make our ODEs stiff at certain places and our experiments with different differential equation solvers, including Runge-Kutta, BDF, and LSODA, showed that LSODA gave the best performance.

We also evaluated different parameter estimation algorithms. Table 6 shows the efficiency of the different algorithms with respect to Levenberg-Marquardt algorithm, determined to be the best choice. Note that we ensured convergence of approximately the same training error when measuring the running time efficiency.

The tri-state region of Rhode Island, Connecticut, and Massachusetts share a common border and worked in tandem to enact restrictions in March-April to control the pandemic. However, as shown in Figure 8 , their end date of restrictions are not similar. In the majority of the cases, RI had relaxed their restriction 2-3 weeks earlier.

Based on the caseloads and predictions at the time, our solution estimated that in mid-June early-July, CT and MA were safe to To validate these numbers, we examine the model's compartmental values around the time of June reopening. Figure 8 shows the population normalized fitted compartment values of (susceptible), (isolation), and + (total infectious) from March-1 to October-20 ℎ . It also superimposes the end date of the restrictions. Comparing the prevalence (normalized infection load) across states at the time of reopening, we can observe that RI's prevalence was high, and it stayed relatively high throughout the summer. Positivity rates also remained high (note that this also depends on state policy on who gets tested). Looking at the normalized fraction in and , we see that all states are comparable. Thus, the higher numbers in RI (from August onwards) have to do with the high sustained prevalence, and this was in part because of early relaxation of restrictions (15 days to 1 month in advance). This not only validates that the system correctly estimated the risk scores for RI but also shows that CT would most likely have a second wave earlier than MA again. As expected, the risk score of RI further worsened to 6, while CT and MA eventually slipped to risk level 6 from 3 on October-27 ℎ and October-29 ℎ , respectively. Figure 9 shows the fitting and the prediction of hospitalization and ICU beds for Hillsborough County, FL. The figure is for an engagement between April to November 2020, where we provided projections at a regular cadence (indicated in the legend). The model accuracy predicted a peak incidence on 7/20/2020 and peak ICU on 7/26/2020. The client used these projections to chart a hospital capacity forecast and take steps to handle the surge. The blue zone above and below the total ICU prediction represents the MAPE (5.9% on average). Based on the projections, the hospital decided not to cancel other elective procedures and surgeries in the time period shown.

Counterfactual analysis is important in simulating what-if scenarios. Clients, especially local authorities like hospitals or government, would like to understand the effects of future events (including seasonal holidays) on the state of the pandemic. At the same time, it also allows them to plan and take remedial steps. For instance, Figure 10 shows the base scenario wherein we show the prediction based on mobility and training as of October-26 ℎ , 2020 for Hillsborough County, FL. With the Thanksgiving holidays around the end of November, local government officials may be interested in knowing ways to reduce the spread of COVID-19. As mentioned earlier, the transmission of a virus is primarily based on population mixing, and one potent way to curb the population mixing is to enforce lock-downs. However, lock-downs have a detrimental effect on the economy, and so one needs to thoroughly understand the implication of different degrees of enforcement on the spread so that one can take the most appropriate step without overreacting. Figure 10 shows the change in daily COVID-19 cases with varying mobility. Note that we do not change any other parameters obtained from the base scenario. The most potent intervention for this state of the epidemic appears to be a 10% reduction in the mobility values observed on October-26 ℎ . However, that would require a significant level of curtailment. By comparing the historical case trends, the 7% decrement scenario matches structurally to the events observed in July. A 5% reduction achieves the effect of stayat-home; however, it leads to a slower decline in the subsequent epidemic wave (gray curve). The figure also reveals a minimum threshold level for mobility curtailment of >~2%. In the case of 2% reduction, although the mobility value has decreased, there is still an insufficient movement of people into the compartment to significantly reduce future incidence over the base scenario. Only for higher values does mobility reduction impact future cases. This form of what-if analysis allows the authorities to assess multiple possible interventions and provides insight into their relative efficacies based on measurements of the effects of past NPI's. 

The solution uses Python for modeling and JAVA for data orchestration. We built a Docker container and deployed it in IBM Cloud Kubernetes Cluster so that each geo-unit can be processed by an image independently. Apache NiFi orchestrates the entire data flow. There are two NiFi pipelines, one for the core epidemiological module, which is run once in three days. The second is for the 'Client Independent Analytics', a.k.a community risk pipeline, which is run nightly. A manifest configuration file containing the list of geo-ids and any list of hyper-parameters like model compartment initializer values is used to spin one new pod on the cluster for the data flow of a geo-unit. This step is throttled to 18 geo-units at a time due to cluster performance limitations. Note that each pod itself is highly parallelized to solve multiple compartment initializer combinations concurrently. On completion of a pod, the NiFi pipeline stores the output of the epidemiology run in an IBM COS bucket. The Community Risk NiFi pipeline, which is independently triggered nightly, collects the results from the COS bucket and then triggers the process to download the latest case data and merge it with the epidemiological predictions. After computing the statistics and risk scores, it finally pushes the output to DB2 On Cloud through a Data Access Layer (DAL) for clients to consume. The 'Client Specific Analytics' has a separate private pipeline that uses these outputs with any additional client data. Figure 11 shows the final UI presentation to the end-user.

In the near future, we intend to extend this solution with additional signals like vaccines, demography, regional comorbidity factors, etc. We strongly believe that the solution and findings of this paper would render epidemiologists and decision-makers better equipped under similar conditions in the future.

Age-stratified discrete compartment model of the COVID-19 epidemic with application to Switzerland

Smoothing and differentiation by an adaptive-degree polynomial filter

Fast unfolding of communities in large networks

Comparing ensemble approaches for short-term probabilistic COVID-19 forecasts in the US

Isotonic median regression: a linear programming approach

Forecasting COVID-19 impact on hospital bed-days, ICU-days, ventilator-days and deaths by US state in the next 4 months

The COVID-19 herd immunity threshold is not low: A re-analysis of European data from spring of 2020

LSODA, ordinary differential equation solver for stiff or non-stiff system

ODEPACK, a systematized collection of ODE solvers

Modeling the dynamics of dengue fever

Impact of COVID-19 on people's livelihoods, their health and our food systems

SARS-CoV-2 Transmission From People Without COVID-19 Symptoms

A contribution to the mathematical theory of epidemics

Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period

A method for the solution of certain non-linear problems in least squares

Nikolaos Trichakis, Thomas Trikalinos, Mohammad Fazel Zarandi, and Dimitris Bertsimas. 2020. Overview of DELPHI Model V3 -COVIDAnalytics

An algorithm for least-squares estimation of nonlinear parameters

Preserving elective surgeries in the COVID-19 pandemic and the future

Mobility trends provide a leading indicator of changes in SARS-CoV-2 transmission

Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations

Finding a" kneedle" in a haystack: Detecting knee points in system behavior

Predicting the impact of asymptomatic transmission, non-pharmaceutical intervention and testing on the spread of COVID19 COVID19. medRxiv

Creating and applying SIR modified compartmental model for calculation of COVID-19 lockdown efficiency

The impact of the COVID-19 pandemic on firm performance

Higher Education Responses to Coronavirus (COVID-19)

WNTRAC: Artificial Intelligence Assisted Tracking of Non-pharmaceutical Interventions Implemented Worldwide for COVID-19

COVID-19 Pandemic Navigator Core Model: Overview of methodology and use cases

Modeling, state estimation, and optimal control for the US COVID-19 outbreak

Epidemic Model Guided Machine Learning for COVID-19 Forecasts in the United States