Epidemic management and control through risk-dependent individual contact interventions

Testing, contact tracing, and isolation (TTI) is an epidemic management and control approach that is difficult to implement at scale because it relies on manual tracing of contacts. Exposure notification apps have been developed to digitally scale up TTI by harnessing contact data obtained from mobile devices; however, exposure notification apps provide users only with limited binary information when they have been directly exposed to a known infection source. Here we demonstrate a scalable improvement to TTI and exposure notification apps that uses data assimilation (DA) on a contact network. Network DA exploits diverse sources of health data together with the proximity data from mobile devices that exposure notification apps rely upon. It provides users with continuously assessed individual risks of exposure and infection, which can form the basis for targeting individual contact interventions. Simulations of the early COVID-19 epidemic in New York City are used to establish proof-of-concept. In the simulations, network DA identifies up to a factor 2 more infections than contact tracing when both harness the same contact data and diagnostic test data. This remains true even when only a relatively small fraction of the population uses network DA. When a sufficiently large fraction of the population (≳ 75%) uses network DA and complies with individual contact interventions, targeting contact interventions with network DA reduces deaths by up to a factor 4 relative to TTI. Network DA can be implemented by expanding the computational backend of existing exposure notification apps, thus greatly enhancing their capabilities. Implemented at scale, it has the potential to precisely and effectively control future epidemics while minimizing economic disruption.


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 and future epidemics will need to be fought with non-pharmaceutical interventions (NPIs) [1,2]. 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 [3,4]. 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 [5,6]. 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 [7][8][9]: 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 [7,[10][11][12]. Even in ideal scenarios, contact tracing needs to identify upward of 75% of infections to achieve epidemic control [8,11]. 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 [13,14].
To scale up the contact tracing component of TTI 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 May 10, 2022 3/66 users [15,16]. 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 also because they do not provide users with information except in the rare case when they receive an exposure notification [17]. Nonetheless, where they have been used, these apps have helped prevent the spread of infections [18].
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 assessments of 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 [19].
• 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 [20,21].
Quantification of individual risks is achieved by assimilating data into a model of virus transmission and disease progression defined on the dynamic contact network assembled from proximity data. 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 May 10, 2022 4/66 propagates data farther along the contact network than contact tracing, consistent with models of disease progression and rates of virus transmission. 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 [22,23], and such contact rates are available from the proximity data from mobile devices.
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 [24]. 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).

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 [25]. Compartmental epidemiological models have traditionally been applied on the level of aggregated individuals (e.g., the population of a city or country) [26]; here we follow more recent work and apply the SEIHRD on an individual level on a time-dependent contact network [23,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 May 10, 2022 5/66 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  [36], 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 infections by confirmed cases [37].) The light lines show 20 simulations that only differ by random seeds, with the thicker lines indicating the ensemble mean; thus, they give an indication of sampling variability. The average contact rate for all nodes is reduced by 58% from March 25, 2020 onward to mimic the lockdown effect (blue solid line).
COVID-19 (Tables 2 and 3). While we assign ages to nodes randomly, the realism of the model could be improved with age-stratified contact patterns [34,35]. With the resulting surrogate worlds of contact patterns and disease progression, we explore how individual risk assessment and epidemic management and control can be achieved in what-if scenarios.
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, presumably because the latter undercount actual infections [37]. However, the death rates in the stochastic simulations are close to the NYC death rate (Fig. 2). The hospitalization rates in the simulation also track the actual hospitalization rates closely (Fig. S8).
Thus, the simulated epidemics on the synthetic network reproduce statistics similar to the actual early epidemic in NYC, with realistic parameter choices for transmission rates and disease progression (Methods). Notwithstanding the simplifications of the network structure, this points to the qualitative adequacy of the synthetic network epidemiology model as a testbed for network DA, which makes no a priori assumptions about the structure of the network.

Accuracy of individual risk assessment
To demonstrate the accuracy of individual risk assessments, we assume the network DA platform hasÑ ≤ N users who exchange proximity data with each other, with 25% to 100% of the population in the user base (i.e., 0.25 ≤Ñ /N ≤ 1). In an idealization, the contact patterns of individuals within the user base are assumed to be known completely; the contact patterns of individuals outside the user base are assumed unknown. We also assume the number of external contacts of individuals in the user base to be known, for example, from proximity-sensing devices that can also detect devices of non-users. ( (Fig. 3). Choosing a classification threshold c I means finding a trade-off between wanting a high TPR while keeping FPR and hence PPF low.
In the ideal albeit unrealistic scenario when the user base encompasses the whole population (Ñ /N = 100%), TPRs for April 9 are 12%, 19%, and 47% for a PPF of 8% and test rates from f = 1%, 5%, and 25%. Later in the epidemic, for April 30, TPRs are 13%, 27%, and 59% ( Fig. 3a, b). That is, the classification results improve as the  and 100% of the total population, even though the scenarios with more limited user bases only use contact information for the users, not for non-users (Fig. 3c, d). The results are also insensitive to the user base topology (Fig. S7): classification performance is not substantially affected whether the user base consists of neighborhoods in the total population network (Fig. 3) or of randomly selected nodes (Fig. S9).
To put these results in context, compare them with the following two traditional approaches: • If only users with positive diagnostic tests are classified as positive, TPRs reach 0.8%, 4%, and 22% for test rates f = 1%, 5%, and 25%, respectively, with PPFs 0.07%, 0.4%, and 1.8% on April 9 and corresponding TPRs of 1%, 4%, and 23% with PPFs 0.02%, 0.1%, and 0.7% 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 much higher TPRs than testing alone, especially at low test rates, 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, following standard contact tracing protocols, 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 so-obtained contact-tracing TPRs for April 9 are 2.4%, 11%, and 45% for test rates from f = 1%, 5% and 25%, with PPFs of 1%, 5% and 23% ( Fig. 3a 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 [21]-with infectiousness indicated by elevated temperature readings with 20% sensitivity [20]. Such temperature readings improve the classification when no or few (f = 1%) diagnostic tests are available; however, 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 [21].
The results show that network DA allows identification of a large fraction of infectious individuals, provided widespread testing is available. The improved identification of infectious individuals over traditional methods is insensitive to the fraction of the population covered by the user base, to the user base topology, and to stochastic variability of the epidemic. Network DA extends classification beyond the nearest network neighbors on which contact tracing and exposure notification apps focus. This gives it an advantage especially when testing capacity is limited.
The capability of network DA to identify infectious individuals can be used to tailor individualized contact interventions for epidemic management and control. For epidemic management and control to be effective, however, it is important not only that the classification accuracy is high but also that the user base coverage is sufficiently large so that a large fraction of infectious individuals can be identified in the population, rather than just within the user base.

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 days thereafter; all others in the population, whether app users or not, do not change their behavior. As a baseline for comparison, we present TTI scenarios with the same contact rate reduction but May 10, 2022 13/66 continuing over 14 days after diagnosis or identification as possibly exposed through contact with an infectious individual. For this baseline TTI scenario, an individual is classified as exposed if over the preceding 10 days, they had at least one contact lasting more than 15 minutes with an individual who had a positive diagnostic test; that is, the contact trace stage of this baseline TTI emulates techniques used in exposure notification apps, relying on the same data as those available for network DA in our synthetic examples. For a direct and fair comparison with network DA, TTI compliance is assumed to be confined to the user base. We use uniform testing regimes with test rates f = 1%, 5%, and 25% within the user base. As classification threshold, we choose a fixed threshold c I = 1%, resulting in TPR 40% and PPF 9% when contact interventions commence. Choosing the classification threshold c I adaptively, in response to current prevalence of infectiousness in the population, may further improve the results.
In the idealized but unrealistic case with full user base coverage (Ñ /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. Network DA does not squash daily infections to zero, because the classification threshold c I was chosen as a compromise between wanting a reliable classification with a high TPR while avoiding isolation of a too large fraction of the population with a too high PPF (Fig. 3). For comparison, TTI with 100% compliance does not achieve epidemic control at a test rate f = 1%; at a 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 than TTI.
In the somewhat more realizable case withÑ /N = 75% user base coverage, we simulate a demanding scenario in which testing and contact interventions are confined to the user base; no contact information among non-users is harnessed, and non-users maintain their contact patterns without isolation. In this case, risk-tailored self-isolation still achieves epidemic control at all test rates of f = 1%, 5%, and 25% within the user base ( Fig. 5), and attains 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 highest test rates (f = 25%) also achieves epidemic control, but with a higher isolated fraction of the population. At the 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. TTI fails to achieve epidemic control at a test rate f = 1%.
With a further reduced user base coverage ofÑ /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. 4. TTI here is confined to the same user base as network DA, implying 75% compliance.
( Fig. S11). 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 the population) is in isolation, which creates perverse incentives: it effectively puts the user base, but not others, in a lockdown. TTI with 50% compliance fails to control the epidemic for test rates below f = 5% but still achieves some control at f = 25%, albeit with a higher isolated population fraction than with network DA.
For the yet smaller user base coverage ofÑ /N = 25%, classification remains accurate (Fig. 3); however, here the dominance of non-users within the population, who do not isolate, rules out epidemic control (Fig. S13). As with any epidemic management measure, control cannot be achieved with low compliance rates.
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Ñ /N = 75%, albeit with an adjusted classification threshold and a higher fraction of the population in isolation (Fig. S10). ForÑ /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. S9  In scenarios in which the user base and/or test rates are too small to achieve epidemic control, there is still a pronounced reduction in the cumulative death rate of users relative to the general non-user population (Fig. 6). For test rates f = 1%, 5%, and 25% per day within theÑ /N = 25% user base consisting of neighborhoods in the overall population network, the cumulative death rate is respectively 29%, 48%, and 42% lower than the death rate among non-users. Additionally, although the contact interventions are confined to the user base, the death rate in the non-user population is still reduced by about 50% compared with the no-intervention scenario (Fig. 2). For ã N /N = 25% user base consisting of nodes selected at random, the results are qualitatively similar: Death rates among users relative to non-users are reduced by 47%, 52%, and 56% for respective test rates f = 1%, 5%, and 25% (Fig. S17).
That is, risk-tailored isolation on the basis of network DA generally outperforms TTI as an epidemic management and control approach when both are presented with the same contact and test data. Even when it does not achieve epidemic control because of low compliance rates, it still offers advantages to users in terms of reduced death rates.

Discussion
We have demonstrated a platform concept 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. The accuracy of the risk assessments is largely independent of the fraction of the population using the platform and of the user base topology; it improves with increasing diagnostic test rates, as should be expected.
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 more effective at preventing infections and deaths than a TTI approach that uses the same contact and diagnostic test data. Our experiments were solely based on self-isolation among app users, without considering other public health interventions. As a result, 75% coverage may be a conservative estimate. In reality, multiple non-pharmaceutical interventions will likely be employed simultaneously at the population level, which may reduce the user coverage required to achieve epidemic control.
We have produced a modular codebase that allows for exploration and benchmarking of tools to manage and control epidemics in a synthetic setting. The platform has a relatively low barrier to widespread implementation. It can be realized by expanding the computational backend of existing exposure notification apps.
High-precision proximity data are now available through Bluetooth protocols [15], and lower-precision location data from mobile devices have been exploited commercially for some time. Statistical techniques may be required to optimize the reconstruction of contact networks from such proximity data in practical implementations with imperfect knowledge of contact patterns [38]. To be effective, the platform requires that users 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. Uptake rates of exposure notification apps have already reached up to 75% in some urban areas, as in our simulated scenarios (e.g., more than 90% of Singapore's population over 6 years of age [39] is using an exposure notification app). Uptake rates on a national scale so far have been more modest (e.g., a third of the UK population [18]), in part, for example, because of rural-urban digital divides but also, probably, because of the limited information provided by current exposure notification apps. However, smartphone usage rates worldwide are around 50% and continue to grow rapidly [40]; thus, widespread use of network DA in future epidemics will become technically possible. And while routine surveillance test rates in much of the world are still low, more widespread surveillance testing on the scale of major cities or regions at this point is feasible; for example, NYC currently is already testing up to roughly 2.5% of its population daily [41]. Our conclusions provide further evidence of the benefits of widespread testing, especially when that is combined with network DA to spread the test information over dynamic contact networks assembled from proximity data.
Challenges to widespread and successful adoption of a network DA platform center around equity, compliance, and privacy questions. Smartphone use is not equitably distributed within the population, and there are disincentives (e.g., unavailability of sick leave) to comply with individual contact interventions. Conversely, classification of users as "low risk" may encourage risky and counterproductive behavior. It is also privacy. The network DA platform requires data to be transferred temporarily to a central computing facility for data assimilation [42]. This makes the platform more difficult to harden against malicious exploitation than exposure notification apps, which only require central data exchange when there is direct evidence of an infection [43].
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 [44].
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 alternative to manual TTI programs, and a backend that delivers more accurate and actionable information than current digital TTI and exposure notification programs developed by many governments [39,45]. The effectiveness of such programs has been modelled [7][8][9], but their impact in practice is only beginning to be elucidated [18]. Given that many TTI programs are voluntary, and documentation of contacts in manual programs is subjective, it will be important to compare both the control and cost effectiveness of manual and digital 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 [3,4]. 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.

SEIHRD model on a contact network
We consider a population of N individuals i (with i = 1, . . . , N ). At any time t, an individual i is in exactly one of six possible health and vital states: 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 τ 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 [22,23,27,46]. 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 κ I ji = a ji (t)β, and that from a hospitalized node with H j (t) = 1 is κ H ji = a ji (t)β. The parameter β 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 (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, to reflect the rarity of SARS-CoV-2 transmission in hospitals in which PPE is worn [47].
(In reality, however, depending on the types of PPE and adherence to hygiene protocols, the degree of transmissibility reduction may vary substantially among hospitals [47].) A typical value for the transmission rate of respiratory viruses is around β = 0.5 h −1 = 12 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 τ , e.g., for an infectious node as [49] T ji (τ ) = 1 − e −κ I ji τ .
(This holds provided the contact period τ is short relative to the duration of infectiousness, so that the infectiousness status of a node does not change during In the model used for DA, we do not assume perfect knowledge of the transmission rate; instead, we learn a partial transmission rate β i for each node i, and compute transmission rates from node j to node i as the averages κ I ji = 0.5a ji (t)(β i + β j ) and κ H ji = 0.5a ji (t)(β i + β j ). We assume independent normal priors for β i for each node, with a mean of 12 day −1 and a standard deviation of 3 day −1 . We truncate these distributions to [1 day −1 , 20 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 σ 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 [6,10,50]. We take the latent period σ −1 i to be fixed for each node i but heterogeneous across nodes. In the model used for DA, we represent it as σ −1 i = 1 day + l i , where l i has a gamma prior distribution with shape parameter k = 1.35 and scale parameter θ = 2 day; hence, the minimum latent period is 1 day, and its prior mean value is 3.7 days (1 day + kθ).
• Duration of infectiousness in community: Infectious nodes with I i (t) = 1 transition to resistant, hospitalized, or deceased at the rate γ 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 [10], but its distribution has a long tail, for example, from individuals with serious or critical disease progression [12]. Like σ i , we take γ 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 γ −1 i = 1 day + g i , where g i has gamma prior distribution with shape parameter k = 1.1 and scale parameter θ = 2 days; hence, the minimum duration of infectiousness is 1 day, and its prior mean value is 3.2 days [10,12].
• Duration of hospitalization: Hospitalized nodes with H i (t) = 1 transition to resistant or deceased at the rate γ i , which is the inverse of the duration of hospitalization. As before, we take γ i to be fixed for each node i but heterogeneous across nodes. In the model used for DA, we model the duration of May 10, 2022 23/66 hospitalization as γ −1 i = 1 day + g i , where g i has a gamma prior distribution with shape parameter k = 1.0 and scale parameter θ = 4 days; hence, the minimum duration of hospitalization is 1 day, and its prior mean value is 5 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 γ i . This implies that, over a period ∆t that is short relative to the duration of infectiousness γ −1 i , the probability of transitioning from being infectious to hospitalized, relative to the total probability of leaving the infectious state, is We take h i to be fixed for each node i but heterogeneous across nodes; it generally depends on age and other risk factors [25,51]. 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 of hospitalized nodes with H i (t) = 1 die. More precisely, we assume infectious nodes die at the rate d i γ i , and hospitalized nodes die at the rate d i γ i .
Both d i and d i are fixed for each node but are heterogeneous across nodes, depending on age and other risk factors [25,51]. Both in the stochastic surrogate-world simulation and in the model used for DA, we assume the same age-dependent mortality rates (Table 3). 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 [52], as in previous studies [46,48,53,54]. 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 S i (t) , E i (t) , etc. associated with the Bernoulli random variables for the states. That is, S i (t) 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 [23,55] where is the total infectious pressure on node i from within the network formed by theÑ node i that are not part of the user network, by the probability w i of an edge of node i being active, and by the time-dependent prevalence of infectiousness P (t), estimated from the network ofÑ users as described below in (11). The probability of exogenous infection then increases with the prevalence of infectiousness P (t) within the user base, which is taken as a proxy of prevalence outside the user base. In an idealization that may not be achievable in practice, we take the number of external neighbors k x i as given from the network structure. In practice, the number of external neighbors can be estimated through use of the same proximity technologies (e.g., Bluetooth) on which exposure notification apps rely, which allow the sensing of other nearby mobile devices, even if they do not participate in the proximity sensing and exposure notification protocol. While this is unlikely to yield perfect knowledge about the number of external neighbors, it may be combined with statistical approximations [38]. The net effect of these assumptions and approximations is that 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. We have confirmed in simulations that exact knowledge of the number of neighbors can be replaced by statistical knowledge; for example, replacing the node-dependent k x i by the user-network average for all nodes (external connectivity in Table 1) yields similar results (Figs. S15, S16).
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 S i (t)I j (t) and S i (t)H j (t) in the infectious pressure (Eq. 2g).
Various approaches to closing this term have been proposed [23,27,55]. Our approach is to estimate it from the ensemble used in the DA cycle, as follows.
The joint-event probability S i (t)I j (t) and the marginal probabilities S i (t) and May 10, 2022 26/66 I j (t) in the master equations (Eq. 2) are related through the ratio which is the rescaled joint probability of S i (t) and I j (t). We estimate the rescaled joint probability C S i (t), I j (t) by its ensemble analogue where (·) = M −1 m (·) denotes the mean over the ensemble (with m = 1, . . . , M labeling ensemble members). Thus, we approximate the joint probability in the infectious pressure (Eq. 2g) as With this empirical approximation, we obtain a closed-form expression for the second moment S m i (t)I m j (t) for each ensemble member m, which we use in the reduced master equations. The second moment S m i (t)H m j (t) follows analogously. If C M S i (t), I j (t) = 1 and C M S i (t), H j (t) = 1, this reduces to the mean-field approximation that is commonly made in epidemiological models [23,55] and that often is accurate on real-world networks [56].
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) [28], which has previously been used with epidemiological models [5,10,29,30]. EAKF treats an  [57], it makes no assumptions about the network structure, and it scales well to high-dimensional problems [28].
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 ∆, with EAKF updates between each pass; a similar scheme has been used in epidemiology models before [5,10,29,30]. In this way, we effectively use EAKF as a smoother, harnessing all available data in a DA window (t f − ∆, t f ). There are two parts to the DA procedure: 1. Update stage: An EAKF update step is performed to assimilate all data available for the window (t f − ∆, 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 − ∆. 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 S m 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 α = 0.0016 and β = 1 (not to be confused with the transmission rate β). 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 − ∆ using the poorest fidelity data (e.g., temperature sensor data), followed by a forecast over • In a second EAKF pass, we update parameters and states at t f − ∆ using moderate-accuracy diagnostic test data, followed by another forecast over • In a final EAKF pass, we update parameters and states at t f − ∆ using data about hospitalization and death status with small error variances, followed by a final forecast over (t f − ∆, 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 [58]: overestimation of long-range covariances, underestimation of variances, and ensemble collapse.
1. To prevent spurious long-range covariances, we localize the effect of observations on states within a single node [58,59]. That is, direct updates of a nodal state are only due to observations at that node during the DA window. This also provides May 10, 2022 29/66 large computational savings because EAKF updates may be performed sequentially node-by-node, in any order.
2. To prevent underestimation of variances by the finite-size ensemble, which can lead to discounting of data points [28], and to ensure well-posedness of the matrix inversions, we use regularization of the ensemble covariance matrix Σ. If Λ min and Λ max denote the minimum and maximum eigenvalues of Σ, we replace Σ in the EAKF algorithm with with Σ + max(δ(Λ max − Λ min ), δ min )I. We choose δ = 5/M to assimilate diagnostic test data, and δ = 1/M to assimilate hospitalization/death status; δ min is taken to be the mean observational noise standard deviation of the update.
3. To prevent ensemble collapse, we add a hybrid inflation to an assimilated state with a map x → a(x −x) +x + N (0, bx), wherex is the ensemble mean state and N (0, bx) is Gaussian noise with mean zero and standard deviation bx [58]. 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 the beginning of a DA window t f − ∆ ≤ t ≤ 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. We generate a static degree-corrected stochastic block model (SBM) [60,61], 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. As pointed out in [63], when groups are ignored, degree distributions associated with social interactions are well-described by a negative binomial distribution, which, for example, has also been used to describe degree distributions associated with sexual-contact networks [64]. In the presence of groups, however, degree distributions of social-interaction networks have been found to exhibit a power-law tail with an exponent of about 2.5 [63]. In accordance with the results presented in [63,65], we therefore describe parts of the synthetic contact network by a stochastic block model with power-law degree correction with exponent 2.5, mean degreek = 10, and maximum degree 100; Figure 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.
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 varying mean rate A ji (t); contact durations are exponentially distributed with a mean contact duration τ (i.e., a mean rate parameter µ = τ −1 ). We choose a mean duration of τ = 2 min (hence µ = 720 day −1 ), based on high-resolution human contact data [48]. We model the mean edge activation rate as Here, t = 0 starts at midnight, andk is the mean degree of the network in the community group (c), so thatkA ji , when averaged over edges, is an average contact rate per node. The diurnally averaged edge activation rate then is For the minimum and maximum contact rates per node, λ i,min and λ i,max , we choose the default values λ min = 4 day −1 and λ max = 84 day −1 . If the default contact rates apply for all nodes, this gives for the community group (c) a mean contact rate per node ofkĀ ji ≈ 37.7 day −1 ; this is about a factor 3-4 larger than typical human contact rates as assessed by self-reports [67], 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 λ i,min and λ i,max for a node i are the principal parameters we vary to explore the effect of contact interventions. Reducing λ i,min and λ i,max for a node reduces the fraction of time edges connecting to node i are active. The contact rate and total contact duration over the network for five simulated days are displayed in Figure S6.
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 µ 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

Selection of subnetwork for user base
To select a subnetwork for a user base withÑ < 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 S7 illustrates the different user base topologies, and Table 1 summarizes their characteristics.
From  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 σ −1 and duration of infectiousness in the community γ −1 are approximated from those in refs. [10] and [12]; the duration of hospitalization γ −1 is from ref. [68], and the transmission rate β is fit to be consistent with data for respiratory viruses [48] and to roughly reproduce NYC data.

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 [52]. 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.
The dependencies of the hospitalization rate h i = h(a) and mortality rates d i = d(a) and d i = d (a) on age a are estimates based on recent data (  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 [69]. The age-dependent death rate in hospitals d is obtained from cumulative hospitalization and death rates in NYC by June 1, 2020 [36], 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 [37]. We choose hospitalization rates h(a) that approximate data from France [70], adjusting the rates so that the overall hospitalization rate is a f (a)h(a) ≈ 3.1%, which is NYC's overall hospitalization rate if one assumes a cumulative COVID-19 incidence rate of 23% [71], together with NYC's actual hospitalization count (52,333 on June 1, 2020) and population (8.34 million) [36]. Similarly, the mortality rate in the community d(a) is chosen such that the overall infection fatality rate is a f (a) d(a) + h(a)d (a) ≈ 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). epidemic in NYC [36]. We find that a global value of β = 12 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Ñ ≤ 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 I i (t) = PPV (9) at the time the test sample is taken. The positive predictive value (PPV) is calculated as 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 [72,73]. 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Ñ and over the ensemble of size M , The cutoff ofÑ −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 − PPV for a positive test result.
• A negative virus test for node i is similarly taken to imply at the time the test sample is taken. The false omission rate (FOR) is calculated with the same sensitivity, specificity, and prevalence as for a positive virus test.
For the DA, we assume an error rate equal to FOR for a negative test result.
• To assimilate low-fidelity data such as those from temperature sensors, we assume I i (t) = 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 [20]. For the DA, we assume an error rate equal to beginning of a DA window t f − ∆ ≤ t ≤ 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 − ∆ ≤ t ≤ 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 with the positive predictive value calculated from (10). Typical values for sensitivity would be 90% and for specificity 95% [74], 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 − 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 − sensitivity) and specificity (false positive rate = 1 − specificity).
Testing strategy We use a simple testing strategy that randomly tests a given budget of nodes once per day. Our framework provides a testbed for different strategies.
We found random testing consistently outperformed three other simple strategies: (i) concentrating the test budget on near-neighbors of positively tested nodes; (ii) continuous testing of a fixed subset of the population; and (iii) testing nodes with high predicted infectiousness values. We attribute this to the low prevalence of disease.
However, it is possible that more effective testing strategies can be discovered that exploit estimated nodal states, the network structure, and the intervention strategy. In a real-world scenario, systematic biases in testing (e.g., testing biased toward certain workplaces or educational institutions) may also affect quantitative details of our results.

May 10, 2022
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 β i ; • Individual inverse latent periods σ i ; • Individual inverse durations of infectiousness γ i and hospitalization γ i ; • Individual hospitalization rates h i and mortality rates d i and d i ; • Exogenous infection rates η i .
We have not fully explored the efficacy of learning about the different parameters from data. For now, we include only the partial transmission rates β i , the inverse latent periods σ i , and the durations of infectiousness γ i and hospitalization γ i in the DA, all with the priors stated above. Figure S18 shows the prior distributions of the parameters at the beginning of the epidemic, as well as 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. When priors were not initially centered on the true values, they remained biased during the simulation. Further investigations focusing, for example, on learning statistical averages of parameters rather than individual node-per-node parameters would be beneficial.
The hospitalization rates h i and mortality rates d i and d i are fixed at the same values as in the stochastic surrogate-world simulation (Table 3). We assume the exogenous infection rates η i to be equal to the partial transmission rates β i , and we estimate the probability of an edge of node i being active as whereĀ i is the diurnally averaged edge activation rate, 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 I i = 1 for which I i = 1 in the stochastic simulation) against the predicted positive fraction (PPF, fraction of nodes with I i = 1 in the user base of sizeÑ ). ROC curves are traced out by lowering the classification threshold c I , thereby increasing both TPR and PPF.
Classification in TTI For the TTI scenarios, we assume the dynamic contact network among users is known, as in the network DA scenarios, and we assume instantaneous tracing. When a node i in the community group (c) is tested positive, it is classified as infectious; all nodes that have had at least one 15-minute contact with node i within the preceding 10 days are classified as exposed. All infectious and exposed community nodes are immediately isolated. This TTI scheme mimics the methods of typical exposure notification apps; although it is idealized and overestimates TTI May 10, 2022 39/66 performance achievable in practical settings [8,75], it provides a fair baseline for comparison with network DA.

Contact Interventions
We implement two types of intervention scenarios in our test cases. In the first, a lockdown scenario (Fig. 2), we set λ i,max for all nodes in the community group (c) to 33 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 λ i,max = λ i,min = 4 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 [76], This corresponds to reinstatement of original contact rates for 50%, 90% and >99% of isolated nodes within 7, 14 and 21 days, respectively.