††thanks: Joint first author††thanks: Joint first author

Also at ]Department of Computational Medicine, University of California, Los
Angeles, CA

# Epidemic Management and Control Through Risk-Dependent Individual Contact
Interventions

Tapio Schneider  California Institute of Technology, Pasadena, CA  Oliver R.
A. Dunbar  California Institute of Technology, Pasadena, CA  Jinlong Wu
California Institute of Technology, Pasadena, CA  Lucas Böttcher
Computational Social Science, Frankfurt School of Finance and Management,
Frankfurt a. M., Germany  [  Dmitry Burov  California Institute of Technology,
Pasadena, CA  Alfredo Garbuno-Iñigo  Departamento de Estadística, Instituto
Tecnológico Autónomo de México, Ciudad de México, México  Gregory L. Wagner
Massachusetts Institute of Technology, Cambridge, MA  Sen Pei  Department of
Environmental Health Sciences, Mailman School of Public Health  Chiara Daraio
California Institute of Technology, Pasadena, CA  Raffaele Ferrari
Massachusetts Institute of Technology, Cambridge, MA  Jeffrey Shaman
Department of Environmental Health Sciences, Mailman School of Public Health

###### Abstract

Abstract

Testing, contact tracing, and isolation (TTI) is an epidemic management and
control approach that is difficult to implement at scale. Here we demonstrate
a scalable improvement to TTI that uses data assimilation (DA) on a contact
network to learn about individual risks of infection. Network DA exploits
diverse sources of health data together with proximity data from mobile
devices. In simulations of the early COVID-19 epidemic in New York City,
network DA identifies up to a factor 2 more infections than contact tracing
when harnessing the same diagnostic test data. Targeting contact interventions
with network DA reduces deaths by up to a factor 4 relative to TTI, provided
compliance reaches around 75%. Network DA can be implemented by expanding the
backend of existing exposure notification apps, thus greatly enhancing their
capabilities. Implemented at scale, it has the potential to precisely and
effectively control the ongoing or future epidemics while minimizing economic
disruption.

Popular Summary

During the ongoing COVID-19 pandemic, exposure notification apps have been
developed to scale manual contact tracing to larger groups of users. The apps
use proximity data from mobile devices to automate finding and notifying
direct contacts of an infection source. They suffer from lack of trust as they
provide users only with rare and binary alerts. Here we present network data
assimilation (DA) as a new digital approach to epidemic management and
control. Network DA uses the same data as exposure notification apps but does
so much more effectively and provides frequently updated individual risk
assessments to users.

Network DA is based on automated learning about individuals’ risk of
infectiousness from crowd-sourced health data and proximity data. The data are
aggregated with models of disease transmission and progression to produce
statistical assessments of users’ risks. In an extensive simulation study of
the COVID-19 epidemic in New York City (NYC), we show that network DA with
regular diagnostic testing could have achieved epidemic control with fewer
than half the deaths that actually occurred during NYC’s lockdown, all while
typically having only 5–10% of the population in isolation at any given time.

Implemented at scale, network DA thus has the potential to effectively control
the ongoing or future epidemics while minimizing the economic and social
disruption of lockdowns.

††preprint: Arxiv preprint

##  I Introduction

Until a majority of the global population has reached immunity against
continuously evolving virus variants through vaccination or infection, the
ongoing COVID-19 pandemic, as well as future pandemics, will need to be fought
with non-pharmaceutical interventions (NPIs). They include social distancing,
mask usage, and restrictions of mass gatherings. But NPIs such as lockdowns
come at catastrophic costs to individuals, economies, and societies, with
disproportionate burdens carried by disadvantaged groups [29, 14]. Even if
imposed only intermittently and regionally, lockdowns are an inefficient means
of epidemic management and control: they isolate much of the population,
although even at extreme epidemic peaks, only a small fraction is infectious
[42, 28]. If individuals who are at high risk of being infectious could be
identified before they infect others, control measures could be made more
efficient by targeting them to this high-risk group.

Testing and contact tracing have been discussed and partly implemented as
strategies to identify individuals who are at high risk of being infectious
[23, 19, 30]: testing determines who is infectious, contact tracing identifies
those who may have been exposed through contact with an infectious individual,
and this high-risk group is then isolated. However, controlling the COVID-19
epidemic by testing, contact tracing, and isolation (TTI) has been complicated
by frequent asymptomatic and presymptomatic transmission, which support silent
spread, and a short serial interval, the period between the onset of any
symptoms in infector and infectee [23, 33, 40, 22]. Even in ideal scenarios,
contact tracing needs to identify upward of 75% of infections to achieve
epidemic control [19, 40]. Quickly diagnosing such a large fraction of
infections and manually identifying exposed individuals requires testing and a
contact tracing workforce at a scale that has been challenging to realize in
most countries [2, 59].

To scale up contact tracing without a massive expansion of the workforce,
exposure notification apps have been developed. They rely on proximity data
from smartphones or other mobile devices to identify close contacts between
users [6, 12]. If an individual user is identified as being infectious, prior
close contacts are notified and can then self-isolate. The exposure
notification is deterministic (a user is only notified when potentially
exposed), and it only uses nearest-neighbor information on the network of
close contacts among users. Exposure notification apps have not seen
widespread use, in part perhaps because of early implementation difficulties
and privacy concerns but likely also because they do not provide users with
information except in the rare case when they receive an exposure
notifications [32]. Nonetheless, where they have been used, these apps have
helped prevent spread of infections [61].

Here we present a new and much more effective way of exploiting the same
information on which exposure notification apps rely. Unlike these apps,
however, this method provides users with continuously updated information
about their individual risks. The core idea is to learn about individual risks
of exposure and infectiousness by propagating crowdsourced information about
infection risks over a dynamic contact network assembled from proximity data
from mobile devices. Instead of the deterministic assessments of exposure
notification apps, our approach exploits data from diverse sources
probabilistically. Various types of information, including their
uncertainties, can be harnessed. For example:

  * •

Diagnostic tests, including sensitive but slow molecular tests, less sensitive
but rapid antigen tests, or pooled diagnostic tests [53].

  * •

Serological tests, which indicate a reduced probability of susceptibility when
antibodies specific to SARS-CoV-2 (or the causative agent of another targeted
disease) are detected.

  * •

Self-reported clinical symptoms, elevated body temperature readings, or other
wearable sensor data, which can indicate an elevated probability of
infectiousness and virus transmission [10, 44].

Quantification of individual risks is achieved by assimilating data into a
model of virus transmission and disease progression defined on the dynamic
contact network. For decision making, periodically updated individual risks of
having been exposed or of being infectious take the place of the deterministic
assessments in exposure notification apps. The probabilistic network approach
propagates data farther along the contact network than contact tracing. It
harnesses more information than contact tracing, both by being able to include
diverse data sources with their uncertainties and by exploiting information
inherent in the network structure itself: an individual with many contacts
generally is at greater risk of having been exposed than an individual with
fewer contacts [35, 39].

The network and the information it contains are dynamically updated in
periodic data assimilation (DA) cycles. These cycles resemble the daily DA
cycles that weather forecasting centers use operationally [8]. The
quantitative information that is provided by the risk assessment platform we
outline in what follows can be used in similar ways as weather forecasts: to
inform personal decisions by users based on their desire to avoid risk (in the
weather forecasting analogy, staying home rather than going on a mountain in
the face of a likely downpour) and to inform public policy when aggregate risk
measures indicate that wider mandates are necessary (analogous to evacuating a
city to protect lives and avoid overwhelming public health and social
infrastructures when a hurricane is likely to make landfall).

##  II Network Data Assimilation

Our point of departure is a variant of the widely used
susceptible–exposed–infectious–resistant (or recovered) (SEIR) model of
epidemiology, extended through inclusion of hospitalized (H) and deceased (D)
compartments to an SEIHRD model [7]. While compartmental epidemiological
models have traditionally been applied on the level of aggregated individuals
(e.g., the population of a city or country) [9], here we apply the SEIHRD on
an individual level on a time-dependent contact network [39, 27]. Each
individual is represented by a node on the network; time-dependent edges
between the nodes are established by close contacts between individuals, as
recorded by proximity data from mobile devices. Virus transmission can occur
during close contacts from infectious or hospitalized nodes to susceptible
nodes, which thereupon become exposed. The probability of transmission
increases with contact duration, and the transmission rate can vary from node
to node and with time, for example, to reflect time-varying transmission rates
resulting from virus mutations or a reduced transmission rate when masks are
worn. From being exposed, nodes progress to becoming infectious, and later
they may progress to requiring hospitalization, recover, or die.

At any time $t$, each node $i$ is in one of the six health and vital states
$S_{i}(t)$, $E_{i}(t)$ etc. of the SEIHRD model (see Methods). Network DA
learns about the probabilities $\langle S_{i}(t)\rangle$, $\langle
E_{i}(t)\rangle$, etc. of finding an individual node $i$ at time $t$ in each
of the different states. We adopt a sequential Bayesian learning approach that
propagates an ensemble of individual probabilities $\langle S_{i}(t)\rangle$,
$\langle E_{i}(t)\rangle$, etc. across the network and periodically updates
them and the SEIHRD model parameters with new data [4, 51, 41, 33]. Data
falling within a DA window of length $\Delta$ (typically $\Delta\approx
1~{}\mathrm{day}$) are incorporated into the model by adjusting the ensemble
to minimize the misfit to the data in the window. An interval $\Delta$ later,
the updating procedure is repeated (see Methods for details). Such DA cycles
and the underlying algorithms are used daily in weather forecasting to
estimate up to $10^{9}$ variables characterizing the state of the atmosphere;
they easily scale to network epidemiology models with millions of nodes or
more. Essentially all types of data and their error characteristics can be
assimilated with this approach, even data that are less sensitive to
infectiousness, such as readings of heart rates [45] or body temperatures [10,
44] (Fig. 1).

Figure 1: Schematic of the personalized risk assessment platform. Proximity-
tracking data from mobile devices is used to assemble a contact network, in
which nodes represent individuals, and edges represent close contacts between
individuals. An epidemiological model defined on the contact network is then
fused with diverse health data, including diagnostic tests, hospitalization
status, and possibly data such as body temperature readings. The model spreads
risk of infectiousness from a positive individual (red) to others, taking into
account knowledge about disease progression, the time and duration of
contacts, and the use of personal protective equipment (PPE), among other
factors. The result of the network DA is an assessment of individual risks,
for example, of being infectious, which then can be used to target contact
interventions.

##  III Synthetic Network for Proof-of-Concept

To demonstrate the methods, we construct a large synthetic contact network
with $N=97,942$ nodes and about 1 million connections among them. The network
has typical characteristics of a human contact network. It has a time-
dependent contact rate minimum at night and a maximum midday, and it has a
connectivity (degree) distribution similar to human networks: there are many
individuals with few connections and a few highly connected individuals who
are more likely to become superspreaders [16] (Fig. S5). The network also
contains a block representing hospitals, where hospitalized patients are
connected to healthcare workers, who in turn are connected to the community in
the rest of the network. However, the present network omits some features of
human contact networks, such as household, workplace, and school structures.

As surrogates for the real world, we use stochastic simulations of the
epidemiology model for the state variables $S_{i}(t)$, $E_{i}(t)$, etc. on the
network. To obtain simulations resembling the real world, we randomly assign
ages to nodes according to the age distribution of New York City (NYC), and we
use those ages to model age-dependent disease progression according to current
knowledge about COVID-19 (Tables 2 and 3). With these surrogate worlds, we
explore how individual risk assessment and epidemic management and control can
be achieved in what-if scenarios.

##  IV Results

###  IV.1 Lockdown and world avoided

Figure 2: Evolution of an outbreak in surrogate-world simulations with a
lockdown (blue) and without (orange). The left column shows infection rates
and the right column death rates. Upper row for cumulative counts and lower
row for daily counts, smoothed with a 7-day moving average filter. Red bars
represent confirmed and probable COVID-19 deaths and confirmed infection rates
for New York City [1], with the red axis labels on the right for confirmed
infection counts. Solid lines indicate the corresponding counts in the
simulations, with the black axis labels on the left for infections on the
network. (The axes for infections in the simulations are stretched by a factor
10 relative to the axes for confirmed NYC infections, reflecting the
undercount of total infections by confirmed cases [62].) The light lines show
20 simulations that only differ by random seeds, with the thicker lines
indicating the ensemble mean. The average contact rate for all nodes is
reduced by 58% from March 25, 2020 onward to mimic the lockdown effect (blue
solid line).

As an illustrative example, we simulate the evolution of an epidemic that,
when scaled up from our network size to the NYC population of 8.3 million,
resembles the early evolution of the COVID-19 epidemic in NYC in 2020 (Fig.
2). If the contact rate on the network remains unchanged in the simulations,
the number of infections and deaths rises from early March into April, with
daily deaths reaching a peak of around 21 per 100,000 population in the second
half of April.

However, this world was avoided by a lockdown, which became mandatory in NYC
from March 22 onward. In its wake, the number of daily new cases and deaths
began to decline from mid-April onward (Fig. 2). We can reproduce similar
behavior in the stochastic simulations by reducing the average contact rate of
all nodes by 58% starting March 25 (Fig. 2). The infection rates in the
stochastic simulations exceed the number of confirmed cases in NYC because the
latter undercount actual infections [62]. However, the death rates in the
stochastic simulations are close to the NYC death rate (Fig. 2), and so are
the hospitalization rates (Fig. S7).

The similarity of the synthetic and observed epidemics shows that the
synthetic epidemics can be used to illustrate the concepts of the risk
assessment platform we envisage.

###  IV.2 Accuracy of individual risk assessment

To demonstrate the accuracy of individual risk assessments, we assume the
network DA platform has $\tilde{N}\leq N$ users who exchange proximity data
with each other, with 25% to 100% of the population in the user base (i.e.,
$0.25\leq\tilde{N}/N\leq 1$). For subsets of the $\tilde{N}$ users, we
assimilate results of simulated rapid diagnostic tests from the corresponding
nodes in the surrogate-world simulation. A fraction $f=5\%$, $10\%$, or $25\%$
of the user base is assumed to be tested daily, with results available on the
day of test administration; that is, every user is tested on average every 20,
10, or 4 days (Methods).

We first illustrate network DA in the worst-case scenario of the free-running
synthetic epidemic, in which contact patterns do not change. DA begins on
March 5. We show results for April 9, near the epidemic peak, when about 7% of
the population are infectious, and for April 30, when new infections are
waning (Fig. 2). (In this free-running epidemic, the maximum prevalence of
infectiousness is considerably higher than in the lockdown simulations, in
which prevalence peaks at 1.5%–2%.)

Figure 3: Receiver operating characteristic (ROC) curves for classification as
possibly infectious. ROC curves trace out the true positive rate (TPR) vs the
predicted positive fraction (PPF) as the classification threshold is varied.
TPR and PPF are given relative to the user base size $\tilde{N}$. Green shades
of the ROC curves from lighter to darker correspond to increasing diagnostic
test rates $f$. Left column for April 9; right column for April 30. (a, b) For
the ideal user base of $\tilde{N}/N=100\%$. For comparison, the filled circles
are for a test-only scenario when only users with positive diagnostic tests
are classified as positive; the open circles are for a contact-tracing
scenario in which additionally prior close contacts of users with positive
diagnostic tests are classified as positive. Also shown is a sensors-only
scenario in which 75% of the user base is assumed to provide daily body
temperature readings. (c, d) For user bases of 25%, 50%, and 75%, with the
same test rates $f$ in shades of green as in (a, b). The black dashed line
represents a random classifier.

We classify users $i$ as possibly infectious (“positive”) when the estimated
probability of infectiousness $\langle I_{i}(t)\rangle$ exceeds a threshold
$c_{I}$. The true positive rate (TPR) at which users are correctly classified
as positive (i.e., they are indeed infectious in the stochastic surrogate-
world simulation) increases as the classification threshold $c_{I}$ decreases;
at the same time, the positive predicted fraction (PPF)—the rate at which
users overall are classified as positive, whether correctly or
incorrectly—also increases because the false positive rate (FPR) increases.
Receiver operating curves (ROC) trace out these competing changes in TPR and
PPF (or FPR) as the classification threshold $c_{I}$ is varied (Fig. 3).
Choosing a classification threshold $c_{I}$ means finding a trade-off between
wanting a high TPR while keeping PPF and FPR low. In the ideal albeit
unrealistic scenario when the user base encompasses the whole population
($\tilde{N}/N=100\%$), TPRs for April 9 range from 19% to 47% for a PPF of 8%
and test rates from $f=5\%$ to $25\%$; for April 30, TPRs range from 27% to
59% (Fig. 3a, b). That is, the classification results improve as the network
model learns about the evolution of the epidemic. The classification results
are insensitive to the user base coverage $\tilde{N}/N$ (Fig. 3c, d), as well
as to the user base topology; for example, classification performance is not
substantially affected whether the user base consists of neighborhoods in the
total population network (Fig. 3) or of random nodes (Fig. S8).

To put these results in context, consider the following two traditional
approaches:

  * •

If only users with positive diagnostic tests are classified as positive, TPRs
reach 4%, 9%, and 22% for test rates $f=5\%$, $10\%$, and $25\%$,
respectively, with PPFs below 1.8% on April 9 and below 0.6% on April 30 (Fig.
3a, b, solid circles). This test-only TPR is close to but slightly smaller
than $f$ because the test sensitivity is less than 100%. Classification by
network DA can achieve substantially higher TPRs than testing alone, at the
expense of increased but still modest PPFs.

  * •

Contact tracing and exposure notification apps classify as positive users with
positive diagnostic tests, plus their potentially exposed nearest neighbors on
the network. If individuals are classified as positive if, over the 10 days
preceding the diagnosis, they had at least one contact of more than 15 minutes
length with a user who had a positive diagnostic test, the contact-tracing
TPRs range from 11% to 45% for test rates from $f=5\%$ to $25\%$, with PPFs
between 5% and 23% (Fig. 3a, b, open circles). Network DA exploits the same
data as contact tracing but achieves substantially higher TPRs at the same
PPF. For example, at the contact-tracing PPF, network DA achieves about a 40%
higher TPR than contact tracing for April 9, and about a factor 2 higher TPR
for April 30 (Fig. 3a, b, vertical lines above open circles).

Network DA can also be used to assess quantitatively to what extent lower-
fidelity data can improve classification. As an example, we conducted a set of
experiments in which 75% of the users were assumed to report body temperatures
daily—for example, with wearable sensors [44]—with infectiousness indicated by
elevated temperature readings with 20% sensitivity [10]. Such temperature
readings improve the classification when no or very few diagnostic tests are
available, but they do not provide a substantial benefit when $f=5\%$ of the
user base or more can be tested daily (Fig. 3a, b). Nonetheless, if widely
adopted, temperature sensors can provide a modest benefit when diagnostic
testing capacity is low [44].

The results show that network DA allows identification of a large fraction of
infectious individuals, provided widespread testing is available. For a
sufficiently large user base coverage $\tilde{N}/N$, this capability can be
used to tailor individualized contact interventions for epidemic management
and control.

###  IV.3 Risk-tailored contact interventions

The individual risk assessments can be used to prompt those who are classified
as possibly infectious for contact interventions. As an illustrative example
of such individual contact interventions, we assume that users of the app
self-isolate by reducing their contact rate with others by 91%, to an average
of 4 contacts per day, during the time when they are classified as positive
and 5 subsequent days; all others in the population, whether app users or not,
do not change their behavior. For comparison, we present TTI scenarios with
the same contact rate reduction but continuing over 14 days after diagnosis or
identification as possibly exposed through contact with an infectious
individual; for a direct comparison with the network DA, TTI compliance is
assumed to be confined to the user base. We focus the presentation on the end-
member testing regimes with test rates $f=5\%$ and $25\%$. As classification
threshold, we choose $c_{I}=1\%$, resulting in $\mathrm{TPR}\gtrsim 40\%$ and
$\mathrm{PPF}\lesssim 9\%$ when contact interventions commence.

In the idealized case with full user base coverage ($\tilde{N}/N=100\%$), the
epidemic is more strongly suppressed with the network DA interventions than in
the lockdown scenario, with 50–70% fewer cumulative deaths (Fig. 4). However,
whereas in the lockdown scenario the entire population has reduced contacts,
with network DA only a small fraction of the population self-isolates. The
self-isolation fraction has an initial peak of 15–17% for about a week and
then falls quickly to 5–10%, with damped relaxation oscillations over several
weeks in the case with lower test rates ($f=5\%$); 50% of those who isolate do
so for 7 days or less, and 90% for 14 days or less. That is, in this idealized
case, risk-tailored self-isolation achieves effective epidemic control with
isolation of only a small fraction of the population at a time. For
comparison, TTI with 100% compliance also achieves epidemic control; however,
at the test rate $f=5\%$, cumulative deaths are 3 times higher than with
network DA because TTI misses more infections than network DA. At the test
rate $f=25\%$, the cumulative death rate with TTI is comparable to or lower
than with network DA, but at the expense of a 2–5 times higher isolated
fraction of the population. Whereas the performance of TTI is strongly test-
rate dependent, that of network DA is less sensitive to test rate, and it is
always more efficient to TTI.

Figure 4: Comparison of different contact intervention scenarios for full user
base with $\tilde{N}/N=100\%$. Shown are the free-running epidemic without
interventions (orange) and the lockdown scenario (blue) from Fig. 2, the
results of network DA and isolation of positive individuals for test rates
$f=5\%$ and $25\%$ (greens), and the results of TTI with the same test rates
as for network DA (purples).

In the more realistic case with $\tilde{N}/N=75\%$ user base coverage, risk-
tailored self-isolation still achieves epidemic control at both test rates
(Fig. 5), with a cumulative death rate similar to the 100% user base. The
fraction of the population in isolation again peaks at just over 15% initially
and then drops to 5–10%. As before, TTI with 75% compliance and with the
higher test rate ($f=25\%$) also achieves epidemic control, but with a higher
isolated fraction of the population; at the lower test rate ($f=5\%$), TTI
results in an about four times higher cumulative death rate than isolation
tailored by network DA, which additionally isolates fewer individuals.

Figure 5: Comparison of different contact intervention scenarios for a user
base with $\tilde{N}/N=75\%$. Plotting conventions as in Fig. 4. TTI here is
confined to the same user base as network DA, implying 75% compliance.

With a user base coverage of $\tilde{N}/N=50\%$, classification remains
accurate and isolation tailored by network DA can still achieve epidemic
control and can remain more effective than a lockdown in preventing infections
and deaths (Fig. S10). The initial fraction of the population in isolation
increases to around 30%, and then drops again to between 5–10%. However, this
means that initially, the majority of the user base (50% of population) is in
isolation, which, in a case of perverse incentives, effectively puts the user
base, but not others, in a lockdown. TTI with 50% compliance fails to control
the epidemic for $f=5\%$ but still achieves control at $f=25\%$, albeit with a
higher isolated population fraction than with network DA.

These results for reduced user bases are for sub-networks consisting of
neighborhoods in the overall population network. Results for user bases
consisting of nodes selected at random from the overall population are
qualitatively similar for $\tilde{N}/N=75\%$, albeit with an adjusted
classification threshold and a higher fraction of the population in isolation
(Figs. S9). For $\tilde{N}/N=50\%$ with a random user base, network DA, while
still being able to identify a large fraction of infectious individuals in the
user base (Fig. S8), ceases to be effective for epidemic control (Fig. S11);
however, it is possible its performance may be improved with data-adaptive
classification thresholds.

##  V Discussion

We have demonstrated the concept of a platform for individual health risk
assessment, which exploits the same proximity data from mobile devices that
exposure notification apps rely upon but is substantially better at
identifying infectious individuals. It achieves these gains by assimilating
crowdsourced data from diverse sources into an epidemiological model defined
on a contact network. Network DA provides informative and actionable risk
assessments for individuals, even when only a modest fraction of the
population uses the app necessary to obtain proximity data. When the user base
is sufficiently large (covering around 75% of the population), the platform
can be used to tailor interventions that are more efficient for epidemic
management and control than lockdowns or TTI. For example, with a user base
covering 75% of the population and users tested every 20 days, simulations for
NYC showed that risk-tailored self-isolation achieves epidemic control with
63% fewer deaths than during NYC’s lockdown, with typically only 5–10% of the
population in isolation at any given time. This risk-tailored isolation
approach is also more effective in preventing infections and deaths than TTI
that uses the same network and diagnostic testing data.

The platform can be realized with existing technology by expanding the backend
of exposure notification apps. High-precision proximity data are now available
through Bluetooth protocols [6], and lower-precision location data from mobile
devices have been exploited commercially for some time. To be effective, the
platform requires that users opt-in to provide proximity data and other
crowdsourced data, such as test results and reports of clinical symptoms. The
more detailed data users make available, the more accurate and detailed risk
profiles can be produced in return.

The data would have to be transferred temporarily to a central computing
facility for data assimilation, posing challenges in protecting users’
privacy. Because of the central processing requirement, it is somewhat more
difficult to harden the platform against malicious exploitation than it is for
exposure notification apps, which only require central data exchange when
there is direct evidence of an infection [24]. Nonetheless, the data need not
be stored beyond a data assimilation window that is at most a few days long.
Additionally, the platform requires only anonymized proximity data but not
absolute location data, and it does not rely on humans in the loop, reducing
risks of malicious exploitation. There may be ways to harden the platform
itself and the data exchange with users against privacy breaches [57].

The network DA platform provides obvious benefits in managing and controlling
epidemics, for example, in reducing the need for lockdowns while preventing
infections and deaths, and in providing users tools to manage their personal
risks. It provides a scalable and more accurate alternative to the digital or
manual TTI programs developed by many governments, including that of NYC [37].
The effectiveness of such programs has been modelled [23, 19, 30], but their
impact in practice is only beginning to be elucidated [61]. Given that many
TTI programs are voluntary, and documentation of contacts is subjective, it
will be important to compare both the control and cost effectiveness of manual
trace programs with the more objective and automated network DA approach
presented here.

In addition to its health impacts, the COVID-19 pandemic has exacted an
enormous economic toll on countries throughout the world [29, 14]. There is a
continuing need to identify approaches that precisely and effectively control
epidemics while minimizing economic disruption. With sufficient uptake and
testing, the platform described here provides a means for achieving these dual
aims.

## Acknowledgments

We thank Tobias Bischoff, Mason Porter, and Andrew Stuart for helpful
discussions. This research was supported by Eric and Wendy Schmidt and Schmidt
Futures (T.S., O.R.A.D., J.W., D.B.); Swiss National Science Foundation
(P2EZP2_191888), National Institutes of Health (R01HL146552), and Army
Research Office (W911NF-18-1-0345) (L.B.); National Science Foundation
(DMS-2027369) (S.P., J.S.) and the Morris-Singer Foundation (J.S.).

## References

  * [1] NYC Coronavirus Disease 2019 (COVID-19) Data.  https://github.com/nychealth/coronavirus-data, 2020.
  * [2] D. Allen, S. Block, J. Cohen, P. Eckersley, M. Eifler, L. Gostin, D. Goux, D. Gruener, V. Hart, Z. Hitzig, J. Klein, J. Langford, T. Nordhaus, M. Rosenthal, R. Sethi, et al., Roadmap to pandemic resilience: Massive scale testing, tracing, and supported isolation (TTSI) as the path to pandemic resilience for a free society, tech. rep., Edmond J. Safra Center for Ethics at Harvard University, April 20 2020.
  * [3] F. Altarelli, A. Braunstein, L. Dall’Asta, A. Lage-Castellanos, and R. Zecchina, Bayesian inference of epidemics on networks via belief propagation, Phys. Rev. Lett., 112 (2014).
  * [4] J. L. Anderson, An ensemble adjustment Kalman filter for data assimilation, Mon. Wea. Rev., 129 (2001), pp. 2884–2903.
  * [5] , Localization and sampling error correction in ensemble Kalman filter data assimilation, Mon. Wea. Rev., 140 (2012), pp. 2359–2371.
  * [6] Apple/Google, Privacy-preserving contact tracing, 2020.  https://www.apple.com/covid19/contacttracing.
  * [7] A. Arenas, W. Cota, J. Gómez-Gardeñes, S. Gómez, C. Granell, J. T. Matamalas, D. Soriano-Paños, and B. Steinegger, Modeling the spatiotemporal epidemic spreading of COVID-19 and the impact of mobility and social distancing interventions, Phys. Rev. X, 10 (2020), p. 041055.
  * [8] P. Bauer, A. Thorpe, and G. Brunet, The quiet revolution of numerical weather prediction, Nature, 525 (2015), pp. 47–55.
  * [9] A. L. Bertozzi, E. Franco, G. Mohler, M. B. Short, and D. Sledge, The challenges of modeling and forecasting the spread of COVID-19, Proc. Natl. Acad. Sci., (2020).
  * [10] M. Bielecki, G. A. G. Crameri, P. Schlagenhauf, T. W. Buehrer, and J. W. Deuel, Body temperature screening to identify SARS-CoV-2 infected young adult travellers is ineffective, Travel Med. Infect. Dis., 37 (2020), p. 101832.
  * [11] C. Brown, A. Noulas, C. Mascolo, and V. Blondel, A place-focused model for social networks in cities, in 2013 International Conference on Social Computing, IEEE, 2013, pp. 75–80.
  * [12] G. Cencetti, G. Santin, A. Longa, E. Pigani, A. Barrat, C. Cattuto, S. Lehmann, M. Salathé, and B. Lepri, Digital proximity tracing on empirical contact networks for pandemic control, Nature Comm., 12 (2021), p. 1655.
  * [13] Centers for Disease Control and Prevention, Contact tracing for COVID-19.  https://www.cdc.gov/coronavirus/2019-ncov/php/contact-tracing/contact-tracing-plan/contact-tracing.html, 2021\.
  * [14] S. Chang, E. Pierson, P. W. Koh, J. Gerardin, B. Redbird, D. Grusky, and J. Leskovec, Mobility network models of COVID-19 explain inequities and inform reopening, Nature, 589 (2021), pp. 82–87.
  * [15] A. Duval, T. Obadia, L. Martinet, P.-Y. Boëlle, E. Fleury, D. Guillemot, L. Opatowski, and L. Temime, Measuring dynamic social contacts in a rehabilitation hospital: effect of wards, patient and staff characteristics, Sci. Rep., 8 (2018), pp. 1–11.
  * [16] A. Endo, Q. J. Leclerc, G. M. Knight, G. F. Medley, K. E. Atkins, S. Funk, and A. J. Kucharski, Implication of backward contact tracing in the presence of overdispersed transmission in COVID-19 outbreaks, Wellcome Open Res., 5 (2021), p. 239.
  * [17] N. M. Ferguson, D. A. Cummings, S. Cauchemez, C. Fraser, S. Riley, A. Meeyai, S. Iamsirithaworn, and D. S. Burke, Strategies for containing an emerging influenza pandemic in Southeast Asia, Nature, 437 (2005), pp. 209–214.
  * [18] N. M. Ferguson, D. A. T. Cummings, C. Fraser, J. C. Cajka, P. C. Cooley, and D. S. Burke, Strategies for mitigating an influenza pandemic, Nature, 442 (2006), pp. 448–452.
  * [19] L. Ferretti, C. Wymant, M. Kendall, L. Zhao, A. Nurtay, L. Abeler-Dörner, M. Parker, D. Bonsall, and C. Fraser, Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing, Science, 368 (2020), p. eabb6936.
  * [20] D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, J. Phys. Chem., 81 (1977), pp. 2340–2361.
  * [21] J. P. Gleeson, S. Melnik, J. A. Ward, M. A. Porter, and P. J. Mucha, Accuracy of mean-field theory for dynamics on real-world networks, Phys. Rev. E, 85 (2012), p. 026106.
  * [22] X. He, E. H. Lau, P. Wu, X. Deng, J. Wang, X. Hao, Y. C. Lau, J. Y. Wong, Y. Guan, X. Tan, et al., Temporal dynamics in viral shedding and transmissibility of COVID-19, Nature Medicine, 26 (2020), pp. 672–675.
  * [23] J. Hellewell, S. Abbott, A. Gimma, N. I. Bosse, C. I. Jarvis, T. W. Russell, J. D. Munday, A. J. Kucharski, W. J. Edmunds, F. Sun, et al., Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts, Lancet Global Health, 8 (2020), pp. E488–E496.
  * [24] R. Hinch, W. Probert, A. Nurtay, M. Kendall, C. Wymant, M. Hall, K. Lythgoe, A. B. Cruz, L. Zhao, A. Stewart, L. Ferretti, M. Parker, A. Meroueh, B. Mathias, S. Stevenson, D. Montero, J. Warren, N. K. Mather, A. Finkelstein, L. Abeler-Dörner, D. Bonsall, and C. Fraser, Effective configurations of a digital contact tracing app: A report to NHSX, tech. rep., April 2020.
  * [25] P. L. Houtekamer and F. Zhang, Review of the ensemble Kalman filter for atmospheric data assimilation, Mon. Wea. Rev., 144 (2016), pp. 4489–4532.
  * [26] B. Karrer and M. E. Newman, Stochastic block models and community structure in networks, Phys. Rev. E, 83 (2011), p. 016107.
  * [27] I. Z. Kiss, J. C. Miller, and P. L. Simon, Mathematics of Epidemics on Networks, Springer, 2017.
  * [28] S. M. Kissler, C. Tedijanto, E. Goldstein, Y. H. Grad, and M. Lipsitch, Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period, Science, 368 (2020), pp. 860–868.
  * [29] M. König and A. Winkler, COVID-19: Lockdowns, fatality rates and GDP growth, Intereconomics, 56 (2021), pp. 32–39.
  * [30] A. J. Kucharski, P. Klepac, A. J. K. Conlan, S. M. Kissler, M. L. Tang, H. Fry, J. R. Gog, and W. J. Edmunds, Effectiveness of isolation, testing, contact tracing, and physical distancing on reducing transmission of SARS-CoV-2 in different settings: a mathematical modelling study, Lancet Infect. Dis., 20 (2020), pp. P1151–1160.
  * [31] S. A. Lauer, K. H. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. Meredith, A. S. Azman, N. G. Reich, and J. Lessler, The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application, Ann. Intern. Med., 172 (2020), pp. 577–582.
  * [32] D. Lewis, Contact-tracing apps help reduce COVID infections, data suggest, Nature, 591 (2021), pp. 18–19.
  * [33] R. Li, S. Pei, B. Chen, Y. Song, T. Zhang, W. Yang, and J. Shaman, Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2), Science, 368 (2020), pp. 489–493.
  * [34] Q.-H. Liu, M. Ajelli, A. Aleta, S. Merler, Y. Moreno, and A. Vespignani, Measurability of the epidemic reproduction number in data-driven contact networks, Proc. Natl. Acad. Sci., 115 (2018), pp. 12680–12685.
  * [35] L. A. Meyers, B. Pourbohloul, M. E. Newman, D. M. Skowronski, and R. C. Brunham, Network theory and SARS: predicting outbreak diversity, J. Theor. Biol., 232 (2005), pp. 71–81.
  * [36] J. Mossong, N. Hens, M. Jit, P. Beutels, K. Auranen, R. Mikolajczyk, M. Massari, S. Salmaso, G. S. Tomba, J. Wallinga, J. Heijne, et al., Social contacts and mixing patterns relevant to the spread of infectious diseases, PLoS Medicine, 5 (2008), p. e74.
  * [37] New York City, Test and trace corps.  https://www1.nyc.gov/site/coronavirus/get-tested/test-trace-corps.page, 2021.
  * [38] M. E. Newman, Spread of epidemic disease on networks, Phys. Rev. E, 66 (2002), p. 016128.
  * [39] R. Pastor-Satorras, C. Castellano, P. V. Mieghem, and A. Vespignani, Epidemic processes in complex networks, Rev. Mod. Phys., 87 (2015), pp. 925–979.
  * [40] C. M. Peak, R. Kahn, Y. H. Grad, L. M. Childs, R. Li, M. Lipsitch, and C. O. Buckee, Individual quarantine versus active monitoring of contacts for the mitigation of COVID-19: a modelling study, Lancet Infect. Dis., 20 (2020), pp. 1025–1033.
  * [41] S. Pei, S. Kandula, W. Yang, and J. Shaman, Forecasting the spatial transmission of influenza in the United States, Proc. Natl. Acad. Sci., 115 (2018), pp. 2752–2757.
  * [42] S. Pei, S. Kandulaa, and J. Shaman, Differential effects of intervention timing on COVID-19 spread in the United States, Sci. Adv., 6 (2020), p. eabd6370.
  * [43] T. P. Peixoto, The graph-tool Python library.  doi:10.6084/m9.figshare.1164194, 2014.
  * [44] G. Quer, J. M. Radin, M. Gadaleta, K. Baca-Motes, L. Ariniello, E. Ramos, V. Kheterpal, E. J. Topol, and S. R. Steinhubl, Wearable sensor data and self-reported symptoms for COVID-19 detection, Nature Medicine, 27 (2021), pp. 73–77.
  * [45] J. M. Radin, N. E. Wineinger, E. J. Topol, and S. R. Steinhubl, Harnessing wearable device data to improve state-level real-time surveillance of influenza-like illness in the USA: a population-based study, Lancet Digital Health, 2 (2020), pp. e85–e93.
  * [46] S. Richardson, J. S. Hirsch, M. Narasimhan, J. M. Crawford, T. McGinn, K. W. Davidson, and the Northwell COVID-19 Research Consortium, Presenting characteristics, comorbidities, and outcomes among 5700 patients hospitalized with COVID-19 in the New York City area, JAMA, 323 (2020), pp. 2052–2059.
  * [47] E. S. Rosenberg, J. M. Tesoriero, E. M. Rosenthal, R. Chung, M. A. Barranco, L. M. Styer, M. M. Parker, S.-Y. J. Leung, J. Morne, D. Greene, et al., Cumulative incidence and diagnosis of SARS-CoV-2 infection in New York, Ann. Epidemiology, 48 (2020), pp. 23–29.
  * [48] M. Salathé, M. Kazandjiev, J. W. Lee, P. Levis, M. W. Feldman, and J. H. Jones, A high-resolution human contact network for infectious disease transmission, Proc. Natl. Acad. Sci., 107 (2010), pp. 22020–22025.
  * [49] H. Salje, C. T. Kiem, N. Lefrancq, N. Courtejoie, P. Bosetti, J. Paireau, A. Andronico, N. Hoze, J. Richet, C.-L. Dubost, et al., Estimating the burden of SARS-CoV-2 in France, Science, 368 (2020), p. eabc3517.
  * [50] N. Sethuraman, S. S. Jeremiah, and A. Ryo, Interpreting diagnostic tests for SARS-CoV-2, JAMA, 323 (2020), pp. 2249–2251.
  * [51] J. Shaman and A. Karspeck, Forecasting seasonal outbreaks of influenza, Proc. Natl. Acad. Sci., 109 (2012), pp. 20425–20430.
  * [52] K. J. Sharkey, Deterministic epidemiological models at the individual level, J. Math. Biology, 57 (2008), pp. 311–331.
  * [53] N. Shental, S. Levy, S. Skorniakov, V. Wuvshet, Y. Shemer-Avni, A. Porgador, and T. Hertz, Efficient high-throughput SARS-CoV-2 testing to detect asymptomatic carriers, Science Adv., 6 (2020), p. eabc5961.
  * [54] Sociopatterns, DATASET: High school dynamic contact networks.  http://www.sociopatterns.org/publications/contact-patterns-among-high-school-students/, 2014\.
  * [55] U.S. Census, Age demographics, New York, NY.  https://datausa.io/profile/geo/new-york-ny#demographics, 2018.
  * [56] U.S. Food and Drug Administration, EUA authorized serology test performance.  https://www.fda.gov/medical-devices/emergency-situations-medical-devices/eua-authorized-serology-test-performance, 2020\.
  * [57] M. B. van der Waal, C. dos S. Ribeiro, M. Ma, G. B. Haringhuizen, E. Claassen, and L. H. M. van de Burgwal, Blockchain-facilitated sharing to advance outbreak R&D, Science, 368 (2020), pp. 719–721.
  * [58] W. Wang, Y. Xu, R. Gao, R. Lu, K. Han, G. Wu, and W. Tan, Detection of SARS-CoV-2 in different types of clinical specimens, JAMA, 323 (2020), pp. 1843–1844.
  * [59] C. Watson, A. Cicero, J. Blumenstock, and M. Fraser, A national plan to enable comprehensive COVID-19 case finding and contact tracing in the US, tech. rep., Johns Hopkins Bloomberg School of Public Health, Center for Health Security, Baltimore, MD, April 2020.
  * [60] Z. Wu and J. M. McGoogan, Characteristics of and important lessons from the coronavirus disease 2019 (COVID-19) outbreak in China: Summary of a report of 72314 cases from the Chinese Center for Disease Control and Prevention, JAMA, 323 (2020), pp. 1239–1242.
  * [61] C. Wymant, L. Ferretti, D. Tsallis, M. Charalambides, L. Abeler-Dörner, D. Bonsall, R. Hinch, M. Kendall, L. Milsom, M. Ayres, C. Holmes, M. Briers, and C. Fraser, The epidemiological impact of the NHS COVID-19 app, Nature, (2021).
  * [62] W. Yang, S. Kandula, M. Huynh, S. K. Greene, G. Van Wye, W. Li, H. T. Chan, E. McGibbon, A. Yeung, D. Olson, A. Fine, and J. Shaman, Estimating the infection-fatality risk of SARS-CoV-2 in New York City during the spring 2020 pandemic wave: a model-based analysis, Lancet Infect. Dis., 21 (2021), pp. 203–212.

## Appendix: Methods

### SEIHRD Model on a Contact Network

We consider a population of $N$ individuals $i$ (with $i=1,\dots,N$). At any
time $t$, an individual $i$ is in exactly one of six possible health and vital
states:

  1. 1.

Susceptible $S_{i}(t)$ when they can get infected with the virus;

  2. 2.

Exposed $E_{i}(t)$ when infected with the virus but not yet infectious;

  3. 3.

Infectious $I_{i}(t)$ when shedding the virus (with or without clinical
symptoms) but not hospitalized;

  4. 4.

Hospitalized $H_{i}(t)$ when hospitalized with active disease, in which case
individuals are assumed to be shedding the virus;

  5. 5.

Resistant $R_{i}(t)$ when immune to the disease through either vaccination or
immunity conferred by a prior infection (we assume lifelong resistance for now
but can relax that assumption if evidence becomes available that immunity is
temporary);

  6. 6.

Deceased $D_{i}(t)$.

We take $S_{i}(t)$, $E_{i}(t)$, $I_{i}(t)$, $H_{i}(t)$, $R_{i}(t)$, and
$D_{i}(t)$ to be Bernoulli random variables that depend on time $t$ and take
only the values $0$ and $1$. That is, $S_{i}(t)=1$ when individual $i$ is
susceptible at time $t$, and otherwise $S_{i}(t)=0$ (and analogously for the
other variables). Because the six SEIHRD states enumerate all health and vital
states of individuals in this model, we have

$S_{i}(t)+E_{i}(t)+I_{i}(t)+H_{i}(t)+R_{i}(t)+D_{i}(t)=1\,.$ (1)

Therefore, there are only 5 independent states.

In the network epidemiology model, a close contact between individuals $i$ and
$j$ establishes a temporary network edge with weight $w_{ji}(t)=1$ for the
duration $\tau$ of the contact; outside the contact period, $w_{ji}(t)=0$.
Transmission along the temporary edges from one node to another and
transitions between health and vital states within each node are modeled as
independent Poisson processes [35, 17, 39, 27]. Each process is characterized
by a rate that may vary from node to node and may depend on external variables
such as age, sex, and medical risk factors (see Fig. S1 for a schematic).

We make the following assumptions about the transmission rate and the
parameters characterizing transition rates between SEIHRD states, including
prior distributions used in the network model for DA:

  * •

_Transmission rate_ : During the contact period between an infectious or
hospitalized individual ($I_{j}(t)=1$ or $H_{j}(t)=1$) and a susceptible
individual ($S_{i}(t)=1$), virus can be transmitted across the edge between
nodes $j$ and $i$.When transmission occurs, the susceptible node $i$ becomes
exposed and switches state to $E_{i}(t)=1$. During the contact period in which
an edge is active ($w_{ji}(t)=1$), we assume the transmission rate to a
susceptible node with $S_{i}(t)=1$ from an infectious node with $I_{j}(t)=1$
is $\kappa_{ji}^{I}=a_{ji}(t)\beta$, and that from a hospitalized node with
$H_{j}(t)=1$ is $\kappa_{ji}^{H}=a_{ji}^{\prime}(t)\beta$. The parameter
$\beta$ is a transmission rate across active edges, which we set to a global
constant in the stochastic surrogate-world simulations and learn on a nodal
basis in the model used for DA; $a_{ji}(t)$ and $a_{ji}^{\prime}(t)$ are time-
dependent transmission modifiers that can be adjusted to incorporate
additional information that may be available, for example, user-supplied
information that individual $i$ is using PPE at time $t$. In our proof-of-
concept simulations, we use $a_{ji}(t)=0.1$ within hospitals and $a_{ji}=1$
otherwise. A typical value for the transmission rate of respiratory viruses is
around $\beta=0.5~{}\mathrm{h}^{-1}=12~{}\mathrm{day}^{-1}$ [48].

Because we model transmission as a Poisson process, the probability that
transmission occurs during contact increases with the duration of the contact
period $\tau$, e.g., for an infectious node as [38]

$T_{ji}(\tau)=1-e^{-\kappa_{ji}^{I}\tau}\,.$

(This holds provided the contact period $\tau$ is short relative to the
duration of infectiousness, so that the infectiousness status of a node does
not change during contact.)

In the model used for DA, we do not assume perfect knowledge of the
transmission rate; instead, we learn a partial transmission rate $\beta_{i}$
for each node $i$, and compute transmission rates from node $j$ to node $i$ as
the averages $\kappa_{ji}^{I}=0.5a_{ji}(t)(\beta_{i}+\beta_{j})$ and
$\kappa_{ji}^{H}=0.5a_{ji}^{\prime}(t)(\beta_{i}+\beta_{j})$. We assume
independent normal priors for $\beta_{i}$ for each node, with a mean of
$12~{}\mathrm{day}^{-1}$ and a standard deviation of $3~{}\mathrm{day}^{-1}$.
We truncate these distributions to
$[1~{}\mathrm{day}^{-1},20~{}\mathrm{day}^{-1}]$, though in practice these
bounds are rarely reached.

  * •

_Latent period_ : Exposed nodes with $E_{i}(t)=1$ transition to being
infectious with $I_{i}(t)=1$ at the rate $\sigma_{i}$, which is the inverse of
the latent period: the time it takes for an exposed individual to become
infectious. For COVID-19, the latent period lies between about 2 days and
about 12 days [31, 33, 28]. We take the latent period $\sigma_{i}^{-1}$ to be
fixed for each node $i$ but heterogeneous across nodes. In the model used for
DA, we represent it as $\sigma_{i}^{-1}=1~{}\mathrm{day}+l_{i}$, where $l_{i}$
has a gamma prior distribution with shape parameter $k=1.35$ and scale
parameter $\theta=2~{}\mathrm{day}$; hence, the minimum latent period is
$1~{}\mathrm{day}$, and its prior mean value is 3.7 days
($1~{}\mathrm{day}+k\theta$).

  * •

_Duration of infectiousness in community_ : Infectious nodes with $I_{i}(t)=1$
transition to resistant, hospitalized, or deceased at the rate $\gamma_{i}$,
which is the inverse of the duration of infectiousness in the community (i.e.,
outside hospitals). For COVID-19, the median duration of infectiousness is
around 3.5 days [33], but its distribution has a long tail, for example, from
individuals with serious or critical disease progression [22]. Like
$\sigma_{i}$, we take $\gamma_{i}$ to be fixed for each node $i$ but
heterogeneous across nodes. In the model used for DA, we model the duration of
infectiousness as $\gamma_{i}^{-1}=1~{}\mathrm{day}+g_{i}$, where $g_{i}$ has
gamma prior distribution with shape parameter $k=1.1$ and scale parameter
$\theta=2~{}\mathrm{days}$; hence, the minimum duration of infectiousness is
$1~{}\mathrm{day}$, and its prior mean value is $3.2~{}\mathrm{days}$ [33,
22].

  * •

_Duration of hospitalization_ : Hospitalized nodes with $H_{i}(t)=1$
transition to resistant or deceased at the rate $\gamma_{i}^{\prime}$, which
is the inverse of the duration of hospitalization. As before, we take
$\gamma_{i}^{\prime}$ to be fixed for each node $i$ but heterogeneous across
nodes. In the model used for DA, we model the duration of hospitalization as
$\gamma_{i}^{\prime-1}=1~{}\mathrm{day}+g_{i}^{\prime}$, where
$g_{i}^{\prime}$ has a gamma prior distribution with shape parameter $k=1.0$
and scale parameter $\theta=4~{}\mathrm{days}$; hence, the minimum duration of
hospitalization is $1~{}\mathrm{day}$, and its prior mean value is
$5~{}\mathrm{days}$. We assume hospitalized nodes are infectious. (If there is
evidence that a hospitalized patient no longer sheds the virus, this can be
taken into account by setting the transmission rate modifier $a_{ji}(t)$ from
the corresponding node to zero; however, we are not considering this situation
in our proof-of-concept.)

  * •

_Hospitalization rate_ : We assume a fraction $h_{i}$ of infectious nodes with
$I_{i}(t)=1$ requires hospitalization after becoming infectious. More
precisely, we assume that infectious nodes transition to becoming hospitalized
at the rate $h_{i}\gamma_{i}$. This implies that, over a period $\Delta t$
that is short relative to the duration of infectiousness $\gamma_{i}^{-1}$,
the probability of transitioning from being infectious to hospitalized,
relative to the total probability of leaving the infectious state, is

$\frac{1-e^{-h_{i}\gamma_{i}\Delta t}}{1-e^{-\gamma_{i}\Delta t}}\approx
h_{i}\quad\text{for}\quad\gamma_{i}\Delta t\ll 1\,.$

We take $h_{i}$ to be fixed for each node $i$ but heterogeneous across nodes;
it generally depends on age and other risk factors [7, 60]. We model the age
dependence in the stochastic surrogate-world simulations according to clinical
data as described below (Table 3), and we assume the same parameters in the
model used for DA.

  * •

_Mortality rate_ : We assume a fraction $d_{i}$ of infectious nodes with
$I_{i}(t)=1$ and a fraction $d_{i}^{\prime}$ of hospitalized nodes with
$H_{i}(t)=1$ die. More precisely, we assume infectious nodes die at the rate
$d_{i}\gamma_{i}$, and hospitalized nodes die at the rate
$d_{i}^{\prime}\gamma_{i}^{\prime}$. Both $d_{i}$ and $d_{i}^{\prime}$ are
fixed for each node but are heterogeneous across nodes, depending on age and
other risk factors [7, 60]. Both in the stochastic surrogate-world simulation
and in the model used for DA, we assume the same age-dependent mortality rates
(Table 3).

  * •

_Resistance_ : For now, we assume resistance to be lifelong, so that an
individual who becomes resistant remains so indefinitely and does not return
to being susceptible. This assumption can be relaxed by allowing transitions
back to the susceptible state if resistance is not permanent.

The health and vital states and transition rates define a Markov chain for the
individual-level SEIHRD states. The SEIHRD Markov chain on a contact network
can be simulated directly with kinetic Monte Carlo methods [20], as in
previous studies [17, 18, 48, 34]. We use kinetic Monte Carlo simulations both
to benchmark a model for the SEIHRD probabilities and to provide a surrogate
for the real world in our proof-of-concept simulations.

### Reduced Master Equations

We are principally interested in the individual SEIHRD probabilities, which
are the expected values $\langle S_{i}(t)\rangle$, $\langle E_{i}(t)\rangle$,
etc. associated with the Bernoulli random variables for the states. That is,
$\langle S_{i}(t)\rangle$ is the probability that individual $i$ is
susceptible at time $t$.

These probabilities could be obtained as averages over an ensemble of kinetic
Monte Carlo simulations; however, it is more computationally efficient to
solve reduced master equations for the probabilities directly. The equations
are [52, 39]

$\displaystyle\langle\dot{S}_{i}\rangle$
$\displaystyle=-\left[\zeta_{i}+k_{i}^{x}\langle w_{i}\rangle
P(t)\eta_{i}\right]\langle S_{i}\rangle\,,$ (2a)
$\displaystyle\langle\dot{E}_{i}\rangle$
$\displaystyle=\left[\zeta_{i}+k_{i}^{x}\langle w_{i}\rangle
P(t)\eta_{i}\right]\langle S_{i}\rangle-\sigma_{i}\langle E_{i}\rangle\,,$
(2b) $\displaystyle\langle\dot{I}_{i}\rangle$ $\displaystyle=\sigma_{i}\langle
E_{i}\rangle-\gamma_{i}\langle I_{i}\rangle\,,$ (2c)
$\displaystyle\langle\dot{H}_{i}\rangle$ $\displaystyle=h_{i}\gamma_{i}\langle
I_{i}\rangle-\gamma_{i}^{\prime}\langle H_{i}\rangle\,,$ (2d)
$\displaystyle\langle\dot{R}_{i}\rangle$
$\displaystyle=(1-h_{i}-d_{i})\gamma_{i}\langle
I_{i}\rangle+(1-d_{i}^{\prime})\gamma^{\prime}_{i}\langle H_{i}\rangle\,,$
(2e) $\displaystyle\langle\dot{D}_{i}\rangle$
$\displaystyle=d_{i}\gamma_{i}\langle
I_{i}\rangle+d_{i}^{\prime}\gamma^{\prime}_{i}\langle H_{i}\rangle\,,$ (2f)
where $\zeta_{i}(t)=\frac{\sum_{j=1}^{N}w_{ji}(t)(\kappa_{ji}^{I}\langle
S_{i}(t)I_{j}(t)\rangle+\kappa_{ji}^{H}\langle
S_{i}(t)H_{j}(t)\rangle)}{\langle S_{i}\rangle}$ (2g)

is the total infectious pressure on node $i$ from within the network formed by
the $\tilde{N}$ users. The infectious pressure represents the possibility of
transmission to node $i$ from all network nodes that are at least temporarily
connected with node $i$. Additionally, we have included an exogenous infection
rate $\eta_{i}$. This allows for infection from outside the network of
$\tilde{N}$ users when the master equation network represents only a subset of
a larger network with $N$ nodes, and so transmission can occur from
unaccounted nodes. The exogenous infection rate $\eta_{i}$ is scaled by the
number of external neighbors $k_{i}^{x}$ of node $i$ that are not part of the
user network, by the probability $\langle w_{i}\rangle$ of an edge of node $i$
being active, and by the time-dependent prevalence of infectiousness $P(t)$,
estimated from the network of $\tilde{N}$ users as described below in (11).
Thus, a user surrounded by other users will have no exogenous infection rate,
while users with many external neighbors will have a larger exogenous
infection rate. The probability of exogenous infection increases with the
prevalence of infectiousness within the user base, taken as a proxy of
prevalence outside the user base.

We integrate these ordinary differential equations with a Runge-Kutta-Fehlberg
4(5) scheme, with an adaptive timestep of maximum length 3 hours. The weights
$w_{ji}(t)$ vary on shorter timescales. This is taken into account in the
numerical integration by averaging $w_{ji}(t)$ over the length of a time step,
rather than evaluating $w_{ji}(t)$ at discrete intervals.

### Closure of Reduced Master Equations

The master equations (2) for the probabilities are not closed because they
depend on the joint probabilities $\langle S_{i}(t)I_{j}(t)\rangle$ and
$\langle S_{i}(t)H_{j}(t)\rangle$ in the infectious pressure (Eq. 2g). Various
approaches to closing this term have been proposed [52, 39, 27]. Our approach
is to estimate it from the ensemble used in the DA cycle, as follows.

The joint-event probability $\langle S_{i}(t)I_{j}(t)\rangle$ and the marginal
probabilities $\langle S_{i}(t)\rangle$ and $\langle I_{j}(t)\rangle$ in the
master equations (Eq. 2) are related through the ratio

$\mathcal{C}\bigl{[}S_{i}(t),I_{j}(t)\bigr{]}=\frac{\langle
S_{i}(t)I_{j}(t)\rangle}{\langle S_{i}(t)\rangle\langle I_{j}(t)\rangle},$ (3)

which is the rescaled joint probability of $S_{i}(t)$ and $I_{j}(t)$. We
estimate the rescaled joint probability
$\mathcal{C}\bigl{[}S_{i}(t),I_{j}(t)\bigr{]}$ by its ensemble analogue

$\mathcal{C}_{M}\bigr{[}S_{i}(t),I_{j}(t)\bigr{]}=\frac{\overline{\langle
S_{i}(t)\rangle\langle I_{j}(t)\rangle}}{\overline{\langle
S_{i}(t)\rangle}~{}\overline{\langle I_{j}(t)\rangle}},$ (4)

where $\overline{(\cdot)}=M^{-1}\sum_{m}(\cdot)$ denotes the mean over the
ensemble (with $m=1,\dots,M$ labeling ensemble members). Thus, we approximate
the joint probability in the infectious pressure (Eq. 2g) as

$\langle S_{i}(t)I_{j}(t)\rangle=\langle S_{i}(t)\rangle\,\langle
I_{j}(t)\rangle\,\mathcal{C}_{M}\bigl{[}S_{i}(t),I_{j}(t)\bigr{]}.$ (5)

With this empirical approximation, we obtain a closed-form expression for the
second moment $\langle S_{i}^{m}(t)I_{j}^{m}(t)\rangle$ for each ensemble
member $m$, which we use in the reduced master equations. The second moment
$\langle S_{i}^{m}(t)H_{j}^{m}(t)\rangle$ follows analogously. If
$\mathcal{C}_{M}\bigl{[}S_{i}(t),I_{j}(t)\bigr{]}=1$ and
$\mathcal{C}_{M}\bigl{[}S_{i}(t),H_{j}(t)\bigr{]}=1$, this reduces to the
mean-field approximation that is commonly made in epidemiological models [52,
39] and that often is accurate on real-world networks [21].

We verified this closure against direct kinetic Monte Carlo simulations of the
SEIHRD model on the synthetic network described below, in the free-running NYC
simulation without lockdown (Fig. 2). The closure has similar performance as
the mean-field approximation and adequately, albeit not perfectly, captures
the stochastic network dynamics (Figs. S2 and S3). The closure correction
coefficients ((4)) concentrate close to the value of 1 (Fig. S4), which
explains the similar performance to the mean-field approximation.

### Data Assimilation Algorithm

For DA, we use a version of the ensemble adjustment Kalman filter (EAKF) [4],
which has previously been used with epidemiological models [51, 41, 33, 42].
EAKF treats an ensemble of $M$ model parameters and states $\langle
S^{m}_{i}(t)\rangle$, $\langle E^{m}_{i}(t)\rangle$, etc. ($m=1,\dots,M$) from
a previous DA cycle as a prior and then linearly updates the ensemble of model
parameters and states to obtain an approximate Bayesian posterior given the
available data. Unlike other algorithms for computing Bayesian posteriors on
networks [3], it makes no assumptions about the network structure, and it
scales well to high-dimensional problems [4].

To learn about parameters and the states of nodes on the network, we use a
scheme based on iterating forward passes of the master equations over a time
window $\Delta$, with EAKF updates between each pass; a similar scheme has
been used in epidemiology models before [51, 41, 33, 42]. In this way, we
effectively use EAKF as a smoother, harnessing all available data in a DA
window $(t_{f}-\Delta,t_{f})$. There are two parts to the DA procedure:

  1. 1.

Update stage: An EAKF update step is performed to assimilate all data
available for the window $(t_{f}-\Delta,t_{f})$, using the previous ensemble
model run as prior. The mismatch between the simulated ensemble trajectories
and the data is used to update the combined ensemble of parameters and states
at the initial time $t_{f}-\Delta$.

  2. 2.

Forecast stage: The updated ensemble of states $\langle
S^{m}_{i}(t_{f}-\Delta)\rangle$, $\langle E^{m}_{i}(t_{f}-\Delta)\rangle$,
etc., with the updated model parameters, is integrated forward up to time
$t_{f}$, to serve as prior for the next update cycle.

EAKF relies on linear updates and assumes Gaussian error statistics. However,
the forward equations (2) are nonlinear. As a result, the EAKF update does not
always conserve total probability, in the sense that SEIHRD probabilities for
each node will not always sum to 1. We therefore augment the state with an
additional total probability conservation variable, with observation equal to
the target probability sum 1. The Gaussian assumption is also at odds with
probabilities in $[0,1]$. We have experimented with approaches of transforming
variables to an unbounded space, leading to total probability conservation
becoming highly nonlinear. We found it to be more robust to work in the
original space where total probability conservation is a linear constraint.
This approach does, however, violate Gaussianity assumptions about the
ensemble when we reinforce the probability bounds by clipping the states
$\langle S^{m}_{i}(t_{f}-\Delta)\rangle$, $\langle
E^{m}_{i}(t_{f}-\Delta)\rangle$, etc. to $[0,1]$.

We assume data errors to be uncorrelated, so that their error covariance
matrix is diagonal (see below for how we specify error variances in the proof-
of-concept simulations). Prior information on parameters and states is
specified in EAKF through the initial condition of the ensemble. We draw the
parameters of the ensemble from the above-specified prior distributions, and
we initialize the state by seeding each ensemble member with a fraction of
(possibly different) infectious nodes, the rest being susceptible. The initial
fraction of infectious nodes is drawn from a beta distribution with shape
parameters $\alpha=0.0016$ and $\beta=1$ (not to be confused with the
transmission rate $\beta$). The mean fraction (here, $0.16\%$) of initially
infected nodes agrees with the fraction of initially infected nodes in the
stochastic surrogate-world simulations.

To account for the multi-fidelity nature of the assimilated data, we perform
EAKF in multiple passes. This allows for better conditioned data covariance
matrices and for different hyperparameter choices for the different types of
data. We perform the following passes to assimilate data from the lowest to
the highest fidelity:

  * •

In a first EAKF pass, we update parameters and states at $t_{f}-\Delta$ using
the poorest fidelity data (e.g., temperature sensor data), followed by a
forecast over $(t_{f}-\Delta,t_{f})$;

  * •

In a second EAKF pass, we update parameters and states at $t_{f}-\Delta$ using
moderate-accuracy diagnostic test data, followed by another forecast over
$(t_{f}-\Delta,t_{f})$;

  * •

In a final EAKF pass, we update parameters and states at $t_{f}-\Delta$ using
data about hospitalization and death status with small error variances,
followed by a final forecast over $(t_{f}-\Delta,t_{f})$.

There are three well-established challenges that ensemble-based filters must
tackle when assimilating a number of parameters/states that is large relative
to the ensemble size [25]: overestimation of long-range covariances,
underestimation of variances, and ensemble collapse.

  1. 1.

To prevent spurious long-range covariances, we localize the effect of
observations on states within a single node [5, 25]. That is, direct updates
of a nodal state are only due to observations at that node during the DA
window. This also provides large computational savings because EAKF updates
may be performed sequentially node-by-node, in any order.

  2. 2.

To prevent underestimation of variances by the finite-size ensemble, which can
lead to discounting of data points [4], and to ensure well-posedness of the
matrix inversions, we use regularization of the ensemble covariance matrix
$\Sigma$. If $\Lambda_{\min}$ and $\Lambda_{\max}$ denote the minimum and
maximum eigenvalues of $\Sigma$, we replace $\Sigma$ in the EAKF algorithm
with with $\Sigma+\max(\delta(\Lambda_{\max}-\Lambda_{\min}),\delta_{\min})I$.
We choose $\delta=5/M$ to assimilate diagnostic test data, and $\delta=1/M$ to
assimilate hospitalization/death status; $\delta_{\min}$ is taken to be the
mean observational noise standard deviation of the update.

  3. 3.

To prevent ensemble collapse, we add a hybrid inflation to an assimilated
state with a map $x\mapsto a(x-\bar{x})+\bar{x}+N(0,b\bar{x})$, where
$\bar{x}$ is the ensemble mean state and $N(0,b\bar{x})$ is Gaussian noise
with mean zero and standard deviation $b\bar{x}$ [25]. We take $a=3.0$ and
$b=0.1$.

Because of the binary nature of the hospitalization and death data, we do not
update these states directly; doing so can lead to shocks in the system
dynamics. We only update the SEIR states $\langle
S_{i}(t_{f}-\Delta)\rangle,\langle E_{i}(t_{f}-\Delta)\rangle,\langle
I_{i}(t_{f}-\Delta)\rangle,\langle R_{i}(t_{f}-\Delta)\rangle$ at the
beginning of a DA window $t_{f}-\Delta\leq t\leq t_{f}$ when assimilating
hospitalization and death data that fall within the DA window.

### Synthetic Network for Proof-of-Concept

We generate a synthetic time-dependent contact network in two steps:

  1. 1.

We generate a static degree-corrected stochastic block model (SBM) [26, 43],
consisting of $N$ nodes in three groups. The three groups represent (a)
hospitalized patients, (b) healthcare workers with contacts both within
hospitals and in the community, and (c) the community of all remaining
individuals (e.g., people in an urban environment). Hospital beds in group (a)
are filled when infected nodes become hospitalized; we assume an infinite
supply of hospital beds. Healthcare workers in group (b) make up 5% of all
nodes, and the remaining 95% of nodes constitute group (c).

We describe connections within groups (a) and (b) with an Erdős–Rényi model
and use mean degrees of 5 in group (a) and 10 in group (b), based on a social-
contact analyses in a hospital setting [15]. Hospitalized patients in group
(a) can interact only with each other and with the healthcare workers in group
(b). To model the interactions between groups (a) and (b), we set the
corresponding mean degrees per node to 5 for edges connecting the groups.

We parameterize the contacts among nodes in the community group (c) with a
power-law degree correction (corresponding to what is sometimes called a
“scale-free” degree distribution) and set the exponent to $2.5$, the mean
degree to $\hat{k}=10$, and the maximum degree to $100$ [11]; Fig. S5 shows
the degree distribution. The community (c) as a group only interacts with
healthcare workers (b), and we set the corresponding mean degree to 5.

  2. 2.

To model time-dependence of the network, we make the edges of the static SBM
network created in the first step time-dependent by switching them on and off.
That is, neighbors of all nodes remain fixed in time, but their connections
are activated or deactivated with time. We account for day/night cycles in the
edge weights $w_{ji}(t)$, but we ignore, e.g., weekly cycles. We generate a
diurnal cycle that replicates some properties of observed time-dependent human
contact networks [54]: The fraction of active edges is small at night and in
the early morning hours, reaches a maximum around noon, and approaches small
values again in the evening.

To model the time-dependence of $w_{ji}(t)$, we use a birth-death process
commonly used in queueing theory. The birth-death process is a Markov chain in
which arrivals (edge activations) are inhomogeneous Poisson processes with a
diurnally varying mean rate $A_{ji}(t)$; contact durations are exponentially
distributed with a mean contact duration $\tau$ (i.e., a mean rate parameter
$\mu=\tau^{-1}$). We choose a mean duration of $\tau=2~{}\mathrm{min}$ (hence
$\mu=720~{}\mathrm{day^{-1}}$), based on high-resolution human contact data
[48]. We model the mean edge activation rate as

$A_{ji}(t)=\frac{1}{\hat{k}}\max\Bigl{\\{}\min(\lambda_{j,\min},\lambda_{i,\min}),~{}\min(\lambda_{j,\max},\lambda_{i,\max})\left[1-\cos^{4}\left(\frac{\pi
t}{1~{}\mathrm{day}}\right)\right]^{4}\Bigr{\\}}.$ (6)

Here, $t=0$ starts at midnight, and $\hat{k}$ is the mean degree of the
network in the community group (c), so that $\hat{k}A_{ji}$, when averaged
over edges, is an average contact rate per node. The diurnally averaged edge
activation rate then is

$\bar{A}_{ji}=\frac{1}{1~{}\mathrm{day}}\int_{0}^{1~{}\mathrm{day}}A_{ji}(t)\,\mathrm{d}t.$
(7)

For the minimum and maximum contact rates per node, $\lambda_{i,\min}$ and
$\lambda_{i,\max}$, we choose the default values
$\lambda_{\min}=4~{}\mathrm{day}^{-1}$ and
$\lambda_{\max}=84~{}\mathrm{day}^{-1}$. If the default contact rates apply
for all nodes, this gives for the community group (c) a mean contact rate per
node of

$\hat{k}\bar{A}_{ji}\approx 37.7~{}\mathrm{day}^{-1};$ (8)

this is about a factor 3–4 larger than typical human contact rates as assessed
by self-reports [36], consistent with the fact that we also take fleeting
contacts into account that would likely not be self-reported. The minimum and
maximum contact rates $\lambda_{i,\min}$ and $\lambda_{i,\max}$ for a node $i$
are the principal parameters we vary to explore the effect of contact
interventions. Reducing $\lambda_{i,\min}$ and $\lambda_{i,\max}$ for a node
reduces the fraction of time edges connecting to node $i$ are active.

The time dependencies of all edges $w_{ji}(t)$ are treated as independent. We
update the time-dependence of each edge at midnight every simulated day,
running independent birth-death processes with parameters $A_{ji}(t)$ and
$\mu$ for the next day.

If a node becomes hospitalized, it is deactivated at its previous location in
the network and transferred to the hospital group (a). Hospitalized nodes are
assumed to be infectious. (This assumption may later be relaxed to model the
situation that a patient is no longer infectious but may still be hospitalized
with ongoing disease.)

### Selection of Subnetwork for User Base

To select a subnetwork for a user base with $\tilde{N}<N$ users, we construct
subgraphs with two different topologies. First, a neighbor-based subgraph is
constructed from a randomly selected seed user, by adding all neighbors of
this user to the subgraph, then in a greedy fashion adding all neighbors of
each member of this new subgraph, and so on, until a desired user population
is reached. Second, a random subgraph is constructed by randomly choosing
nodes from the full network. Figure S6 illustrates the different user base
topologies, and Table 1 summarizes their characteristics.

Type | Population | Interior | Exterior connectivity
---|---|---|---
75% neighbor | 73,456 | 21,499 | 1.9
50% neighbor | 48,971 | 6,301 | 5.2
25% neighbor | 22,381 | 2,107 | 10.0
75% random | 73,353 | 7,061 | 3.1
50% random | 48,371 | 550 | 6.3
25% random | 24,482 | 33 | 9.3
Table 1: Details of the different user bases. The percentage represents
approximately $\tilde{N}/N$, for the user population $\tilde{N}$. The interior
defines how many users are completely surrounded by other users. The exterior
connectivity gives the average number of exterior nodes connected to a node
inside the user base.

### Surrogate World Simulation

To generate surrogate worlds with which to test the DA algorithm and
interventions, we simulate epidemics on the synthetic network stochastically
with kinetic Monte Carlo methods [20]. For these stochastic simulations (but
not for the model used for DA), we choose mean transmission and transition
rates between SEIHRD states that are homogeneous across nodes, except for
hospitalization and mortality rates that depend on age. The mean rates we use
are based on current knowledge about COVID-19 (Table 2). We simulate 20
epidemics that only differ in their random seed. They are initialized on March
5, 2020, with 0.16% of nodes randomly assigned to be infectious.

Parameter | Value
---|---
$\beta$, $\kappa^{I}$ | $12~{}\mathrm{day}^{-1}$
$\kappa^{H}$ | $0.1\beta$
$\sigma$ | $(3.7~{}\mathrm{day})^{-1}$
$\gamma$ | $(3.2~{}\mathrm{day})^{-1}$
$\gamma^{\prime}$ | $(5~{}\mathrm{day})^{-1}$
$\lambda_{\max}$ | $84~{}\mathrm{day}^{-1}$
$\lambda_{\min}$ | $4~{}\mathrm{day}^{-1}$
Table 2: Mean transmission and transition rates and maximum/minimum contact
rates for the surrogate-world simulations with the stochastic SEIHRD model
(Fig. S1). The mean rates are taken to be the same for all nodes; hence, the
nodal indices are suppressed. The latent period $\sigma^{-1}$ and duration of
infectiousness in the community $\gamma^{-1}$ are approximated from those in
refs. [33] and [22]; the duration of hospitalization ${\gamma^{\prime}}^{-1}$
is from ref. [46], and the transmission rate $\beta$ is fit to be consistent
with data for respiratory viruses [48] and to roughly reproduce NYC data.

The dependencies of the hospitalization rate $h_{i}=h(a)$ and mortality rates
$d_{i}=d(a)$ and $d^{\prime}_{i}=d^{\prime}(a)$ on age $a$ are estimates based
on recent data (Table 3).

Age group $a$ (yrs) | $f(a)$ | $h(a)$ | $d(a)$ | $d^{\prime}(a)$
---|---|---|---|---
0–17 | 20.7% | 0.2% | 0.0001% | 1.9%
18–44 | 40.0% | 1.0% | 0.001% | 7.3%
45–64 | 24.5% | 4.0% | 0.1% | 19.3%
65–74 | 8.3% | 7.6% | 0.7% | 32.7%
$\geq 75$ | 6.5% | 16.0% | 1.5% | 51.2%
Table 3: Age-dependent mean hospitalization and mortality rates in the
surrogate-world simulation. The share $f(a)$ of the population in each age
group $a$ is taken from U.S. Census data [55]. The age-dependent death rate in
hospitals $d^{\prime}$ is obtained from cumulative hospitalization and death
rates in NYC by June 1, 2020 [1], under the assumption that $90\%$ of deaths
occurred in hospitals. Age-dependent hospitalization rates $h(a)$ and
mortality rates $d(a)$ in the community (outside hospitals) are difficult to
obtain directly from NYC data because of an age-dependent undercount of
infections [62]. We choose hospitalization rates $h(a)$ that approximate data
from France [49], adjusting the rates so that the overall hospitalization rate
is $\sum_{a}f(a)h(a)\approx 3.1\%$, which is NYC’s overall hospitalization
rate if one assumes a cumulative COVID-19 incidence rate of 23% [47], together
with NYC’s actual hospitalization count (52,333 on June 1, 2020) and
population (8.34 million) [1]. Similarly, the mortality rate in the community
$d(a)$ is chosen such that the overall infection fatality rate is
$\sum_{a}f(a)\bigl{[}d(a)+h(a)d^{\prime}(a)\bigr{]}\approx 1.1\%$, which is
NYC’s overall infection fatality rate if one considers the same cumulative
incidence of 23% and the confirmed and probable cumulative death count from
COVID-19 (21,607 by June 1, 2020).

To model age-dependent disease progression, we randomly assign ages to network
nodes in the community group (c) according to the age distribution for NYC, as
given by U.S. Census data [55]. Assumptions about the age dependence of
disease progression and other model parameters are immaterial for the risk
assessment platform—the platform can learn such parameters and parametric
functions from data; we merely use them to create a reasonably realistic
surrogate of the real world. We assign ages to nodes in the healthcare worker
group (b) according to the age distribution among working-age adults (21–65
years old). Initially, there are no hospitalized nodes (i.e., group (a) is
empty).

The network has 97,942 nodes. We choose the global mean transmission rate
$\beta$ so that our simulations are qualitatively aligned with the evolution
of the COVID-19 epidemic in NYC [1]. We find that a global value of
$\beta=12~{}\mathrm{day^{-1}}$ can qualitatively reproduce the observed rate
of infections and can quantitatively reproduce the rate of hospitalizations
and deaths during the initial phases of the epidemic in NYC.

### Synthetic Data

We sample synthetic data from the stochastic surrogate-world simulation on the
network with $N$ nodes and assimilate data for a possible subset of
$\tilde{N}\leq N$ users in the reduced master equation model. We consider the
following data and error rates:

  * •

A positive virus test for node $i$ is taken to imply

$\langle I_{i}(t)\rangle=\mathrm{PPV}$ (9)

at the time the test sample is taken. The positive predictive value (PPV) is
calculated as

$\mathrm{PPV}=\frac{\mathrm{sensitivity}\times
P(t)}{\mathrm{sensitivity}\times
P(t)+(1-\mathrm{specificity})\times(1-P(t))},$ (10)

where we take the sensitivity of the test to be 80% and the specificity to be
99%, which we use as an approximation of the currently imprecisely known
actual sensitivities and specificities [50, 58]. As an estimate of the
infectiousness prevalence $P(t)$ in the population, we use the average of the
infectiousness probabilities both over the network of size $\tilde{N}$ and
over the ensemble of size $M$,

$P(t)=\max\left(\frac{1}{\tilde{N}M}\sum_{m=1}^{M}\sum_{i=1}^{\tilde{N}}\langle
I_{i}^{m}(t)\rangle,\frac{1}{\tilde{N}}\right).$ (11)

The cutoff of ${\tilde{N}}^{-1}$ is included to guard against erroneously
assuming prevalence to be zero because of subsampling on the reduced network.
For the DA, we assume an error rate of $1-\mathrm{PPV}$ for a positive test
result.

  * •

A negative virus test for node $i$ is similarly taken to imply

$\langle I_{i}(t)\rangle=\mathrm{FOR}$ (12)

at the time the test sample is taken. The false omission rate (FOR) is
calculated as

$\mathrm{FOR}=\frac{(1-\mathrm{sensitivity})\times
P(t)}{(1-\mathrm{sensitivity})\times
P(t)+\mathrm{specificity}\times(1-P(t))},$

with the same sensitivity, specificity, and prevalence as for a positive virus
test. For the DA, we assume an error rate equal to $\mathrm{FOR}$ for a
negative test result.

  * •

To assimilate low-fidelity data such as those from temperature sensors, we
assume $\langle I_{i}(t)\rangle=\mathrm{PPV}$ as for a positive virus test
when they indicate infectiousness (e.g., when a temperature reading is
elevated). However, we use a sensitivity of 20% and specificity of 98% to
reflect the lack of sensitivity of temperature sensors in detecting COVID-19
infection [10]. For the DA, we assume an error rate equal to $1-\mathrm{PPV}$,
analogous to a positive virus test.

  * •

Data about hospitalization with COVID-19 imply that $H_{i}(t)=1$ for the
duration of hospitalization. We assume the hospitalization status of all users
to be known with certainty, that is, we assimilate the hospitalization status
$H_{i}(t)=0$ or $H_{i}(t)=1$ for all users; however, we only update the SEIR
probabilities at the beginning of a DA window $t_{f}-\Delta\leq t\leq t_{f}$
with hospitalization data.

  * •

Death implies $D_{i}(t)=1$. We assume the vital status of all users to be
known with certainty, that is, we assimilate $D_{i}(t)=0$ or $D_{i}(t)=1$ for
all users; as for hospitalization, however, we only update the SEIR
probabilities at the beginning of a DA window $t_{f}-\Delta\leq t\leq t_{f}$
with death data.

  * •

For completeness, we state that a positive serological test for SARS-CoV-2 for
node $i$ would be taken to imply

$\langle R_{i}(t)\rangle=\mathrm{PPV},$

with the positive predictive value calculated from (10). Typical values for
sensitivity would be 90% and for specificity 95% [56], and the prevalence of
resistance can be estimated from the resistance probabilities on the reduced
master equation network. We would again assume an error rate equal to
$1-\mathrm{PPV}$. However, we did not assimilate simulated serological tests
in our proof-of-concept because currently achievable serological test rates
are low.

To model data errors, we randomly corrupt the synthetic data sampled from the
surrogate world network with the false positive and false negative rates
implied by the sensitivity (false negative rate = $1-\mathrm{sensitivity}$)
and specificity (false positive rate = $1-\mathrm{specificity}$).

### Testing Strategy

Because of the low prevalence of disease, we use a simple testing strategy
that randomly tests a given budget of nodes once per day. It is possible that
more intelligent targeting strategies for tests, based on the estimated nodal
states and the network structure, would be more effective in finding new
cases. Our framework can provide a testbed for developing such strategies.

#### Parameter Learning

In addition to assimilating probabilities of SEIHRD states, we can in
principle learn about parameters in the reduced master equation model (2), for
example:

  * •

Individual partial and time-dependent transmission rates $\beta_{i}$;

  * •

Individual inverse latent periods $\sigma_{i}$;

  * •

Individual inverse durations of infectiousness $\gamma_{i}$ and
hospitalization $\gamma_{i}^{\prime}$;

  * •

Individual hospitalization rates $h_{i}$ and mortality rates $d_{i}$ and
$d_{i}^{\prime}$;

  * •

Exogenous infection rates $\eta_{i}$.

We have not fully explored the efficacy of learning about all these different
parameters from data. For now, we include only the partial transmission rates
$\beta_{i}$, the inverse latent periods $\sigma_{i}$, and the durations of
infectiousness $\gamma_{i}$ and hospitalization $\gamma^{\prime}_{i}$ in the
DA, all with the priors stated above. Figure S12 shows the prior distributions
of the parameters at the beginning of the epidemic. and the posterior
distributions as the epidemic evolves and the network model learns about the
parameters. The results show that the DA does not refine the prior estimates
of the parameters.

The hospitalization rates $h_{i}$ and mortality rates $d_{i}$ and
$d_{i}^{\prime}$ are fixed at the same values as in the stochastic surrogate-
world simulation (Table 3). We assume the exogenous infection rates $\eta_{i}$
to be equal to the partial transmission rates $\beta_{i}$, and we estimate the
probability of an edge of node $i$ being active as

$\langle w_{i}\rangle=\frac{\bar{A}_{i}}{\mu+\bar{A}_{i}},$ (13)

where $\bar{A}_{i}$ is the diurnally averaged edge activation rate,

$\bar{A}_{i}=\frac{1}{1~{}\mathrm{day}}\int_{0}^{1~{}\mathrm{day}}A_{i}(t)\,\mathrm{d}t,$
(14)

with

$A_{i}(t)=\frac{1}{\hat{k}}\max\left\\{\lambda_{i,\min},\lambda_{i,\max}\left[1-\cos^{4}\left(\frac{\pi
t}{1~{}\mathrm{day}}\right)\right]^{4}\right\\}.$ (15)

With our parameters, this is
$\bar{A}_{i}=\lambda_{i,\min}/\hat{k}=0.4~{}\mathrm{day^{-1}}$ for isolated
nodes and $\bar{A}_{i}=3.77~{}\mathrm{day^{-1}}$ otherwise. For a stationary
birth-death process, this estimate for $\langle w_{i}\rangle$ is the
stationary probability of an edge being active; it approximates the diurnally
averaged probability in the case of the birth-death process with diurnally
varying edge activation rates. Through this probability $\langle
w_{i}\rangle$, the exogenous infectious pressure depends on the isolation
status of a node.

#### Classification of Disease Status

Nodes $i$ in the community group (c) are classified as possibly infectious
($\mathcal{I}_{i}=1$) or not ($\mathcal{I}_{i}=0$) according to

$\mathcal{I}_{i}=\begin{cases}1&\text{if }\overline{\langle
I_{i}^{m}\rangle}>c_{I},\\\ 0&\text{otherwise}.\end{cases}$ (16)

Here, $c_{I}$ is a classification threshold, which can be determined from
receiver operating characteristic (ROC) curves as some optimum tradeoff
between wanting to achieve high true positive rates while keeping false
positive rates modest. The ROC curves we use are adapted to the setting in
which we are primarily interested in the fraction of users that is classified
as possibly infectious (and thus may be asked to self-isolate). We plot the
true positive rate (TPR, nodes with $\mathcal{I}_{i}=1$ for which $I_{i}=1$ in
the stochastic simulation) against the predicted positive fraction (PPF,
fraction of nodes with $\mathcal{I}_{i}=1$ in the user base of size
$\tilde{N}$). ROC curves are traced out by lowering the classification
threshold $c_{I}$, thereby increasing both TPR and PPF.

#### Contact Interventions

We implement two types of intervention scenarios in our test cases. In the
first, a lockdown scenario (Fig. 2), we set $\lambda_{i,\max}$ for all nodes
in the community group (c) to $33~{}\mathrm{day^{-1}}$. This amounts to a
reduction of the mean contact rate ((8)) in group (c) by 58%. In the second, a
time-limited isolation intervention, we reduce the contact rates of targeted
high-risk nodes by setting
$\lambda_{i,\max}=\lambda_{i,\min}=4~{}\mathrm{day^{-1}}$; thus, these high-
risk nodes are assumed to self-isolate, with only 4 contacts per day on
average, corresponding to a reduction of their average contact rate by $91\%$.

The duration of contact reduction takes three possible values. In the lockdown
scenario (Fig. 2), all nodes have contact reduction from the inception of the
lockdown until the end of the simulation. In the TTI scenario, self-isolating
nodes have contact reduction for 14 days, in accordance with current CDC
guidelines [13], after which contact rates are reset to the original values.
For the network DA scenario, self-isolating nodes have contact reduction until
they are classified as negative ($\overline{\langle I_{i}^{m}\rangle}\leq
c_{I}$) for a period of 5 consecutive days, after which contact rates are
reset to their original values. This corresponds to reinstatement of original
contact rates for 50%, 90% and $>$99% of isolated nodes within 7, 14 and 21
days, respectively.

## Supplementary Figures

Figure S1: Schematic of SEIHRD model [7]. Infected and hospitalized nodes
infect susceptible nodes at rates $\kappa^{I}$ and $\kappa^{H}$, respectively.
After being infected, susceptible nodes become exposed. Exposed nodes become
infectious at rate $\sigma$. Infected nodes may get hospitalized at rate
$h\gamma$, die at rate $d\gamma$, or become resistant at rate $(1-h-d)\gamma$.
Once hospitalized, nodes either become resistant at rate
$(1-d^{\prime})\gamma^{\prime}$ or die at rate $d^{\prime}\gamma^{\prime}$.

Figure S2: Overall epidemic dynamics from SEIHRD model using mean-field
approximation.

Figure S3: Overall epidemic dynamics from SEIHRD model using mean-field
approximation with ensemble correction.

Figure S4: Histograms of correction coefficients (top row)
$\mathcal{C}_{M}\bigl{[}S_{i}(t),I_{j}(t)\bigr{]}$ and (bottom row)
$\mathcal{C}_{M}\bigl{[}S_{i}(t),H_{j}(t)\bigr{]}$ at different times during
the simulated epidemic.

Figure S5: Distribution of degrees $k$ in synthetic contact network with
97,942 nodes.

Figure S6: Illustration of different user base topologies. (a) Neighbor-based
user base, constructed by iteratively adding neighborhoods. (b) Random
subnetwork of users. Red nodes and edges are part of a user base, grey nodes
and edges of the overall population. The shown networks have 982 nodes and
5,916 edges. Both user bases contain 5% of all nodes.

Figure S7: Hospitalization rates in surrogate-world simulation with a lockdown
(blue) and without (orange). The left panel shows cumulative hospitalizations
and the right panel daily hospitalizations per 100,000 population for the same
simulations as those in Fig. 2. Red bars represent COVID-19-related
hospitalization rates for New York City [1]. As in Fig. 2, the simulation data
are smoothed with a 7-day moving average filter.

Figure S8: Receiver operating characteristic (ROC) curves for classification
as possibly infectious. As in Fig 3, but for subnetworks with randomly
selected nodes rather than for subnetworks with a neighborhood topology.

Figure S9: Comparison of different contact intervention scenarios for random
user base with $\tilde{N}/N=75\%$. As in Fig. 5, but with a subnetwork with
randomly selected nodes and with a classification threshold $c_{I}=0.25\%$.

Figure S10: Comparison of different contact intervention scenarios for
neighborhood user base with $\tilde{N}/N=50\%$ and with a classification
threshold $c_{I}=0.5\%$.

Figure S11: Comparison of different contact intervention scenarios for random
user base with $\tilde{N}/N=50\%$ and with a lower classification threshold
$c_{I}=0.25\%$. As in Fig. S10, but with a subnetwork with randomly selected
nodes.

Figure S12: Distribution of ensemble averaged model parameters across nodes as
a function of time during the epidemic. The shaded regions contain 50%, 80%
and 90% of the distribution. The dashed line represents the true parameters in
the stochastic simulation. During the first 8 days, no DA is performed, and
the parameter distributions are the prior distributions.