key: cord-0644874-t4qx73xr
authors: Herbrich, Ralf; Rastogi, Rajeev; Vollgraf, Roland
title: CRISP: A Probabilistic Model for Individual-Level COVID-19 Infection Risk Estimation Based on Contact Data
date: 2020-06-09
journal: nan
DOI: nan
sha: 92ea0ee554abc07965ed9fc114a441d9e8bae73c
doc_id: 644874
cord_uid: t4qx73xr

We present CRISP (COVID-19 Risk Score Prediction), a probabilistic graphical model for COVID-19 infection spread through a population based on the SEIR model where we assume access to (1) mutual contacts between pairs of individuals across time across various channels (e.g., Bluetooth contact traces), as well as (2) test outcomes at given times for infection, exposure and immunity tests. Our micro-level model keeps track of the infection state for each individual at every point in time, ranging from susceptible, exposed, infectious to recovered. We develop a Monte Carlo EM algorithm to infer contact-channel specific infection transmission probabilities. Our algorithm uses Gibbs sampling to draw samples of the latent infection status of each individual over the entire time period of analysis, given the latent infection status of all contacts and test outcome data. Experimental results with simulated data demonstrate our CRISP model can be parametrized by the reproduction factor $R_0$ and exhibits population-level infectiousness and recovery time series similar to those of the classical SEIR model. However, due to the individual contact data, this model allows fine grained control and inference for a wide range of COVID-19 mitigation and suppression policy measures. Moreover, the algorithm is able to support efficient testing in a test-trace-isolate approach to contain COVID-19 infection spread. To the best of our knowledge, this is the first model with efficient inference for COVID-19 infection spread based on individual-level contact data; most epidemic models are macro-level models that reason over entire populations. The implementation of CRISP is available in Python and C++ at https://github.com/zalandoresearch/CRISP.

The COVID-19 pandemic has spread rapidly around the world, with the number of infections and deaths steadily growing. Most governments around the world have been completely unprepared to deal with the COVID-19 outbreak, which UN Secretary-General Antonio Guterres has referred to as humanity's worst crisis since World War II. While governments around the world had plans in place in the event of a pandemic, the peculiarities of COVID-19 (e.g., delayed onset of symptoms, asymptomatic transmission) have challenged these preparations. Governments have reacted by implementing measures such as nationwide lock-downs, that require people to stay inside their homes, enforcing social distancing and therefore breaking the COVID-19 infection chain. However, a blunt mechanism such as a lock-down (over an extended period) can cause severe damage to the economy, and so, there is a need to find alternative measures to slow down or stop the spread without incremental effects in other areas of society. These alternatives have to be built in a solid foundation such as widespread testing and the isolation of infected (or potentially infected) individuals via contact-tracing.

Contact tracing technologies [20, 19] have shown promise in tracking the spread of the disease across the population. These mobile apps capture social contact information between users such as contact duration, distance, etc. using Bluetooth signals on devices. The fine-grained contact data of individuals collected by the apps can enable:

• Individual risk score prediction. The contact data, combined with information about COVID-19 positive test cases, can be used to predict the likelihood of infection for each individual. The individual risk scores can be leveraged by governments and organizations to prioritize testing as well as to identify individuals that need to enter isolation/quarantine. • Hotspot detection. Tracing technologies can help authorities identify areas with a high density of contacts and/or individuals with high infection risk. This can allow policymakers to make more effective decisions, for example, by imposing highly restrictive measures such as lock-downs, shelter-at-home, or school closures only in COVID-19 hotspots while allowing activities to remain closer to normal in unaffected areas. • Insights about infection spread. Contact tracing can provide insights into the relative importance of different modalities of disease transmission (e.g., through intermediate surfaces vs individual contact), risk of infection transmission based on contact characteristics such as duration and distance, most likely locations (e.g., schools, work, malls) for the spread of disease, and "super spreaders" who come in close proximity with a large number of individuals and so must be frequently tested for infection.

To achieve the above-mentioned benefits, we need to devise new models and inference algorithms for analyzing contact tracing data. This is because existing epidemics models [3, 15, 14, 6, 5 ] focus on estimating population-level statistics such as percentage of the population infected, number of days for the epidemic to peak, etc. as opposed to the infection state of each individual in the population. Other models [16, 17] that use ML-based inference techniques assume complete knowledge of the infection state of each individual at each time instant. However, in the COVID-19 scenario, (1) the infection status of individuals is not known until they are tested, and (2) infectious time of individuals are unknown since individuals may infect others while asymptomatic. Finally, governments are using contact tracing data [20, 19] to identify and test individuals who have come in direct contact with COVID-19 positive test cases. However, the fact that asymptomatic individuals may have infected a large number of people prior to displaying symptoms and being tested, delays the detection of these newly infected individuals by only using contact tracing apps.

Our main contributions can be summarized as follows:

• We propose CRISP (COVID-19 RIsk Score Prediction), a probabilistic graphical model for COVID-19 infection spread through diverse contacts channels between individuals. Our model uses latent variables to represent the epidemiological states of individuals based on the SEIR model [12] at different points in time, and captures both the transitions between states as well as test outcomes. 

We classify related work into four broad categories: epidemic models, Machine Learning (ML) based inference of model parameters, influence maximization in social networks, and contact tracing apps.

In recent years, there has been research on modeling individual dynamics of epidemics [3, 15, 14] . However, this work typically resorts to mean-field theory to model virus spread over a network, and thus does not characterize the dynamic infectious state of each individual over time.

Ferguson et al. [6, 5] use a compartmental transmission model to simulate the spread of influenza across a population, and analyze interventions such as antiviral prophylaxis and social distancing to halt a pandemic. The authors use a stochastic model of individuals co-located in households that are randomly distributed across a geographical region, and infection risk from 3 sources -household, place and random contacts in the community. The infection transmission rates for the 3 sources and recovery rates are based on analysis of historical data. In contrast, we leverage real individual contact tracing data and outcomes of tests on individuals to infer the infection transmission rate for each contact and the likelihood of infection for each individual.

Lorch et al. [10] propose a spatiotemporal epidemic model that uses marked temporal processes to represent the epidemiological condition of each individual (based on a variation of the SEIR compartment models), individual mobility patterns, test outcomes, and testing and contact tracing strategies. The authors design an efficient sampling algorithm for the model using Monte Carlo roll-outs that is able to predict the spread of COVID-19 under different testing & tracing strategies, social distancing measures, and business restrictions, given contact histories of individuals. They use Bayesian optimization techniques to infer model parameters (e.g. infection transmission rate) that minimize the difference between the real positive COVID-19 cases and those in the Monte-Carlo simulations. In addition, they demonstrate the efficacy of their model using real COVID-19 data and mobility patterns of Tübingen, Germany. Our Monte Carlo EM inference algorithm for model parameters is computationally much more efficient than the Bayesian optimization techniques employed in [10] .

In [16] , the authors consider the problem of inferring latent social networks based on network diffusion or disease propagation data. Given the times when nodes become infected, but not who infected them, the authors identify the optimal network that best explains the observed data. The authors present a maximum likelihood approach based on convex optimization with a L 1 -like penalty term (that encourages sparsity) to estimate the conditional probability of infection transmission between every node pair. A key difference from our work is that [16] assumes complete knowledge of infected nodes and infection times. In contrast, in the COVID-19 scenario, (1) the infection status of nodes is not known until they are tested, and (2) infection times of nodes are unknown since nodes may not show symptoms even though they are infected (and infecting others).

Warriyar et al. [17] introduce a novel R statistical software package EpiILM for simulating infectious disease spread, and carrying out Bayesian MCMC-based statistical inference for spatial and/or (contact) network-based models in the Deardon et al. [4] individual-level modelling framework. In individual-level models (ILMs), the epidemiological state of each individual (e.g., susceptible or infected) is assumed to be perfectly known at each time instant, which makes it relatively straightforward to estimate model parameters such as infection transmission probabilities (as a function of covariates) using maximum likelihood estimation or Bayesian inference using Metropolis-Hastings MCMC. However, in the COVID-19 scenario, epidemiological states of individuals are hidden until they are tested, and this complicates Bayesian inference in our probabilistic model setting. Note that this model has no cycles as we assume the infection status z u,t only depends on variables z v,t before time step t, t < t. However, due to the "memory" that the state z u,t = E and z u,t = I have, we require edges into the entire past of an infection trace.

The Influence Maximization problem aims to select k users in a social network that maximize influence spread, and was first modeled as an algorithmic problem by Kempe et al. [8] . [9] presents a comprehensive survey of different diffusion models that capture the information diffusion process and approximation algorithms to maximize influence. The papers assume that diffusion model parameters such as influence probabilities are given and focus on selecting k users to maximize influence spread. In contrast, our paper focuses on the problem of estimating model parameters related to infection transmission probabilities for each contact, given social contact information between users and COVID-19 test results for users.

[11] addresses the problem of finding the "backbone" of an influence network. It employs network sparsification to preserve only the links that play an important role in the propagation of information. [7] considers the problem of estimating influence probabilities between users in a social graph. Given a social graph and a log of actions by users, the Maximum Likelihood Estimator (MLE) of influence probability of node u on node v is simply the fraction of actions performed by u that are also performed by v. Unlike [7] , in our setting, the infection status and times of nodes are latent, and need to be inferred by our algorithms.

To combat the spread of COVID-19, governments have launched contact tracing apps [20, 19] that use Bluetooth signals on mobile phones to track contacts between users. Users who have come in direct contact with COVID-19 positive test cases are considered to be at high risk of infection, and subject to tests and quarantine actions. However, a key problem with this approach is that COVID-19 infected users are typically tested only after they show symptoms, and typically, infected users show symptoms 5-6 days post infection. These asymptomatic users may have infected a large number of users over multiple hops prior to displaying symptoms and being tested. This delays detection of infected users using contact tracing apps, and limits their effectiveness to proactively test and isolate infected users to contain the spread of COVID-19. In contrast, our probabilistic modeling algorithm CRISP predicts the likelihood of a user getting infected with COVID-19 through a chain of social contacts involving asymptomatic users, and identifies infected users early, even though they may be multiple hops from a user who has tested positive for COVID-19 and even before they begin to show symptoms. Our inference algorithm also learns infection transmission probabilities for each contact channel.

Our CRISP model is an SEIR model at the level of every individual (see [12] for an introduction). Note that we consider discrete time steps t implicitly assumed to be at the level of single days. We assume that we are given the following two datasets for a given set S of individuals:

Here we assume that the feature vector x i describes the overall contact between u i and v i via the number x i,j of mutual contacts over channel j (e.g., Bluetooth encounters, queuing together, sharing public transportation). We assume D contact to be symmetric so that

We model the T discrete time steps of infection status for each individual. Our model has |S| × T many latent variables Z := {z u,t } u∈S,t=1,...,T ∈ {S, E, I, R} |S|×T that represent the four stages of infection 2 :

• z u,t = S: individual u has not been infected and is susceptible, • z u,t = E: individual u is infected but not contagious, • z u,t = I: individual u is infected and is contagious, • z u,t = R: individual u has recovered and is not contagious. Let us use the notation Z u,t := {z u,1 , . . . , z u,t } to denote the set of latent states z u,t of individual u up to and including time t and Z t := u Z u,t . In addition to the latent infection status of all individuals at each time step, we also model

Then, our graphical model G := (V, E) between the variables V := Z O has the following edges:

1. E time = u E u and E u := {(z u,t , z u,t ) t<t }. All edges between the latent infection states of a single individual u. The edges E u will be used to describe the probability of P(z u,t |z u,t−1 , . . . , z u,1 ) and describe the full time series of being susceptible, exposed, infectious and then recovered. 2. E contact := (u,v,t,x)∈Dcontact {(z u,t , z v,t+1 )}. All edges between two individuals u and v who had a contact at time t. 3. E test := (u,t,o)∈Dtest {(z u,t , o u,t )}. All edges between a test outcome at time t and the corresponding infection status z u,t of individual u. These edges will be used to describe the probabilities P(o u,t |z u,t ) of a test outcome given the infection status of u at that same time.

The full edge set E is the union of these three edge types: Figure 3 shows an example graphical model with these three edge components.

In order to define the joint probability distribution, note that all edges E time and E contact are pointing forward in time. Thus, all {z u,t+1 } u∈S are conditionally independent of each other given all the past states Z t . Also, all edges E test have the property that a test outcome of individual u at time t only depends on the infection status of u at t, z u,t . The joint probability distribution is given by

where

Since we are using an SEIR model, the only non-zero probabilities P(z u,t+1 |Z t ) are the transitions S → S, S → E, E → E, E → I, I → I, I → R, R → R.

2 Note that we are assuming that recovered individuals are immune until time step T .

Infection Model In order to define f, we assume that an infection occurs from exogenous influences with a fixed probability p 0 ∈ [0, 1] or with probability of p j ∈ [0, 1] for every instance of a contact through the contact channel j if the contact was already in the state I. Thus, the probability that no infection occurred at time t equals

Infection Status Model In order to define g and h, let us assume we have a point density function

for d E time steps (and similarly for the duration of the infectiousness). Examples of functions q E and q I are the probability mass functions of the binomial, negative-binominal or geometric distributions. However, for the case of COVID-19, we will use discrete probabilities established from analysis of the population in [1] and [18] . Moreover, let

be the conditional probability (according to q) that the duration is exactly n time steps given that the duration is at least n time steps. Then,

Note that the first argument to both g and h is the number of E and I states up to and including time t in the state sequence Z u,t .

Test Outcome Model Finally, we need to define the probability of a test outcome o given the infection status z u,t of individual u at time t. Since there are two types of mistakes of a test, we use

We assume both 0 < α 1 and 0 < β 1. It is easy to implement more sophisticated test accuracy models here, in particular to distinguish between different infection states. Also, we can easily model α and β which are dependent on how many days an individual has been in state I; this change would not affect the block-Gibbs sampling scheme in Section 4 in an adverse way.

Prior Model In order to complete the description of the full probabilistic model, we have to specify P(Z 1 ). For simplicity, we assume these probabilities to be a delta-peak at state S, that is, P(z u,1 = S) = 1 for all u ∈ S.

For inference in the aforementioned model we are interested in computing the infection risk score of every individual u at every time step t given the test outcomes O available as well as the hyperparameters θ := (p 0 , p 1 , . . . , p J ) which cannot be set by knowledge of the diseases: the J parameters p j represent the probabilities of COVID-19 infection transmission through the contact channel j and p 0 captures the probability that an infection occurs at any time-step from exogenous influences.

In order to estimate θ, we will maximize the log-likelihood of the data O, that is

where the second decomposition explicitly contains the posterior P (Z T |O). However, this posterior is not analytically tractable and therefore we will approximate it by performing block-Gibbs sampling of an infection trace z u := (z u,1 , . . . , z u,T ) of individual u keeping all other infection traces {z v:v =u } fixed. This requires a computationally efficient procedure to sample from the conditional probability distribution P (z u | {z v:v =u } , O, θ) which we describe in the following subsection.

Since we assume that the total number of days of the model, T , is not large 3 , we will enumerate all possible sequences of infection traces z u and compute the un-normalized probability of

for all terms that depend on elements of the trace z u in order to re-normalize and draw from this distribution. Also, as our model is an SEIR model, we know that each sample z u can be uniquely represented by a triple ω There are three groups of factors that (might) involve z u in the (un-normalized) conditional probability

The Factors A(z u ) In order to derive a compact representation of A(z u ), we assume that it can be written in terms of

Since the infectious status, z v,t = I of other individuals v that had contact with u only affect u in the susceptible state, we can derive l 0 (t 0 ) and l infected from (4) by collecting the f terms (see (5)):

where p u,t := (v,u,t,x)∈Dcontact:zv,t=I j (1 − p j ) x j . Note that l 0 (t 0 ) is the density function of the geometric distribution. Similarly, given (4) and (7) we can derive l E (d E ) as

A similar derivation shows that l I (d I ) = q I (d I ) which proves that the computational complexity of computing A(z u ) has been reduced to one factor for each contact during the S states of u and three additional factors corresponding to the compact representation ω for z u . The number of factors do not directly scale up with T .

Factors B(z u ) In order to derive a compact representation of B(z u ), we note that only the cases of z v,t = S potentially contain the value of z u,t for v = u (see the value range of the function g and h in (4)). In fact, looking at (5) it becomes evident that it requires z u,t = I. Thus, B(z u ) is defined by

where C S (u, t) and C E (u, t) are the individuals that u met at time t who were susceptible and have either stayed susceptible or got exposed, respectively:

In this subsection, we describe how we can speed up the block-Gibbs sampling step for drawing a sample infection trace z u for individual u.

A key observation is that the term C(z u ) = Contacts into u The remaining term l infected in A(z u ) (see (14) ) captures the contribution due to contacts into individual u from other (infectious) individuals v who are in state z v,t = I prior to u herself getting infected at t 0 . As a result, l infected depends on the values of other infection traces {z v:v =u } and needs to be recomputed during each block-Gibbs sampling step to draw sample z u . Let us define

where p u,t := (v,u,t,x)∈Dcontact:zv,t=I j (1 − p j ) x j (see also (17) ). Then, we have

Thus, at the start of each block-Gibbs sampling step, we pre-compute (22) and (23) for each time step t, and then use (24) to compute l infected for a particular infection trace z u . Note that this involves only t 0 multiplications (or additions in the log-domain).

Contacts out from u We next turn our attention to computing B(z u ) for each infection trace z u that captures the contribution due to contacts out from u. We introduce two states Z I t and Z ¬I t which are identical to Z t except for the value of infection state z u,t which is I is

be the inner terms in (21). Note that for all contacts (u, v, t, x) ∈ D contact , the terms B(z u , t, Z I t ) and B(z u , t, Z ¬I ) differ only in the factor

but not in f(v, t, Z ¬I t ) since u is infectious at this time t in Z I t but not in Z ¬I t . Also, note that values of z u,t for t < t do not affect B(z u , t, Z t ). We can then obtain B(z u ) for each infection trace z u value as

where Constant is the product of B(z u , t, Z ¬I t ) over all time steps t and can be ignored due to normalization of the sampling distribution. Again, the ratio B(z u , t, Z I t )/B(z u , t, Z ¬I t ) can be precomputed for all time steps t at the start of the block-Gibbs sampling step for z u , and then used to compute B(z u ) for each infection trace z u as in (26).

Putting it all together Note that the quantities l infected (t), l infected (t) and B(z u , t, Z I t )/B(z u , t, Z ¬I t ) only depend on the infection status of individual u at time t because the infection traces z v of all other individuals v are fixed when we are drawing a block-Gibbs sample for z u . Thus, the (un-normalized) conditional probability for each z u value is obtained by taking the product of A(z u ), B(z u ) and C(z u ), which in turn are computed efficiently as described above from pre-computed values of l 0 (t 0 ), l E (d E ), l I (d I ) and C(z u ) at the start of the algorithm, and l infected (t), l infected (t) and B(z u , t, Z I t )/B(z u , t, Z ¬I t ) for all time steps t at the start of the block-Gibbs sampling step.

Additional implementation optimizations We use two additional ideas to accelerate the implementation of the block-Gibbs sampling algorithm:

• We never materialize the infection trace z u because it is uniquely described by the triple ω = (t 0 , d E , d I ); each value z u,t can be computed by no more than three comparisons of t with t 0 , t 0 + d E and t 0 + d E + d I . Thus, the whole state of the latent variable model is represented by 3 × |S| integers. • We carry out all computations of probabilities in the log-domain so all functions become sums and products instead of products and powers.

Algorithm 1: Block-Gibbs sampling algorithm for CRISP model /* Initialization */ 1 Initialize each z u = S · 1 /* Precomputations independent of contact data */ 2 forall (t 0 , d E , d I ) ∈ N + × N + × N + such that t 0 + d E + d I ≤ T do 3 Pre-compute l 0 (t 0 ), l E (d E ) and l I (d I ) 4 Construct the sequence z u with t 0 states S, d E states E, d I states I and

Pre-compute C(z u ) according to (13) Pre-compute l infected (t) using (22) and l infected (t) using (23) 10 Pre-compute ratio B(z u , t,

Compute log(A(z u )) = log l 0 (t 0 ) + log l E (d E ) + log l I (d I ) + log(l infected ) using (24) 14 Compute log(B(z u )) using (26) 15 Set l t 0 ,d E ,d I = log(A(z u )) + log(B(z u )) + log(C(z u )) /* Block-Gibbs sampling step */ 16

Set z u with (S,E,I,R) states corresponding to (t * 0 , d * E , d * I ) 18 return Z i = Z 19 until convergence Block Glibbs Sampling Algorithm Algorithm 1 is block-Gibbs sampling algorithm for sampling Z i T from our CRISP model. It cycles through (random) individuals u, sampling the vector of latent variables z u from the conditional distribution P(z u |{z v:v =u }, O, θ) until convergence. We can use the samples Z 1 T , . . . , Z m T drawn by this algorithm to compute the infection risk score for an individual u at time t by taking the fraction of samples Z i T in which the latent infection state z u,t ∈ {E, I}.

In order to estimate the hyper-parameters θ of the CRISP model, would like to find θ * that maximizes the log-likelihood log (12) . However, since this is intractable, we propose to use the Monte Carlo Expectation-Maximization (EM) algorithm [2] . We will use EM to refine θ in successive iterations. Let θ old be the value of θ computed in the previous iteration. Then, in the E step of the current iteration, we will estimate the expected complete-data log-likelihood

We will use the block-Gibbs sampling procedure described in Algorithm 1 to approximate the posterior distribution P(Z T |O, θ old ) over the latent infection status of individuals u. If the samples drawn from the posterior P(Z T |O, θ old ) are Z 1 T , . . . , Z m T , then in the M step, we will compute θ that maximizes the expected complete-data log-likelihood

where z i u,t+1 is the infection state for individual u at time t + 1 in sample Z i T . If t i 0 denotes the number of initial S states in the sample infection trace z i u , we note that by virtue of (4) only the first t i 0 terms depend on θ which reduces the above maximization term to only

We use stochastic gradient descent to compute the θ values that maximize the above expression. We also note that for numerical stability, we re-parameterize p j via w j as p j = exp(w j )/(1 + exp(w j )) which allows for an unconstrained optimization over w.

We can extend the block-Gibbs sampling algorithm in CRISP to a federated learning setting [13] where local contact and test outcome data for an individual u are utilized to compute the block-Gibbs sample z u on the individual's mobile device without ever needing to be shared with anyone else. This has two benefits: (1) We distribute the block-Gibbs sampling algorithm across hundreds of millions of mobile devices in the world and thereby utilize their distributed computational power, and (2) Contact and test outcome data for an individual are stored only on the individual's mobile device and not shared with other mobile devices--this preserves a user's privacy. In the federated setting, the contact data is never centralized-instead for each individual u, her device executes the block-Gibbs sampling step to draw sample z u only using the locally available contacts and test outcome data for u, as well as additional "minimal statistics" sent to u by the devices of its past contacts. In the following two paragraphs, we explain how to compute the factors A(z u ), B(z u ) and C(z u ) in (13) in a federated setting (see Algorithm 2 for the pseudo-code which runs on every mobile device).

Factors A(z u ) and C(z u ) A key observation is that the terms l 0 (t 0 ), l E (d E ), l I (d I ) in A(z u ) as well as the factor C(z u ) can all be computed locally on the device with the contact and test outcome information available on the device. In order to compute the remaining term l infected in A(z u ), we only require information on the individuals v who had a contact with u at each time step t and the infection status z v,t of v at the time of the contact. Individual u's mobile device already has the contact information for u; thus all that is required to compute l infected are the current infection traces z v for all individuals v who have had contacts with u. Since each infection trace is uniquely characterized by a (t 0 , d I , d E ) triple, we require the mobile devices of all individuals v who have had a contact with u to send u's device the (t 0 , d I , d E ) triple corresponding to z v .

Factor B(z u ) In order to compute B(z u ) as defined in (21), we require the term f(v, t, Z t ) for each individual v who has had a contact with u at time t and whose infection state z v,t = S. Let f −u (v, t, Z t ) be defined as in (5) over all contacts of v at time t except for individual u. Then, the device for each individual v who has had a contact with u at time t and whose infection state z v,t = S sends to u's device the quantity f −u (v, t, Z t ) computed based on v's view of the infection traces of its contacts. These terms are used by u's device to compute B(z u ) as defined in (21).

In this section, we present two types of experimental evaluations:

1. Population Level COVID-19 Infection Spread. In the first set of experiments, we will demonstrate that CRISP is capable of modelling infection spread across an entire population. We will relate our individual-level parameters θ to more classical measures of infection spread such as reproduction factor R 0 and demonstrate that the structure of the contact patterns allow more fine grained control of the infection spread which can be used for alternative containment measures of the COVID-19 pandemic.

2. Test and Quarantine Efficacy of CRISP Model. In the second set of experiments, we will assess the test and quarantine efficacy of the CRISP model by comparing the population health after 5 months under three testing and quarantining policies: (1) symptom-based, (2) contact-tracing-based, and (3) CRISP model-based.

Algorithm 2: Federated block-Gibbs Sampling algorithm for CRISP model /* Initialization */ 1 Initialize z u = S · 1 2 Initialize N v = (∞, 0, 0, ∅) for all v // Stores minimal statistics from contacts 3 /* Precomputations independent of contact data */ 4 repeat /* Update minimal statistic received in the incoming queue */

Pre-compute C(z u ) for this sequence according to (13) /* Precomputations dependent on contact data */ 10 forall time steps t do 11 Pre-compute l infected (t) using (22) and l infected (t) using (23) and

Compute log(A(z u )) = log l 0 (t 0 ) + log l E (d E ) + log l I (d I ) + log(l infected ) using (24) 16 Compute log(B(z u )) using (26) 17 Set

Set z u with (S,E,I,R) states corresponding to (t * 0 , d * E , d * I ) /* Send minimal statistic to all contacts */ 20 forall v in past contact list do

In all these experiments, we use the parameters α = 0.001 and β = 0.01 in (9) and match the distribution q E and q I of exposure and infectiousness duration to the empirical distributions provided in the medical literature [1, 18] . This is both used in the generation of the simulated test outcome data as well as for the CRISP inference algorithms as these parameters are publicly known. We will use the notation q I for the expectation of the empirical distributions q I .

In order to simulate realistic epidemiological spread, we need to translate a reproduction factor R 0 at t = 0 into contact data. By definition, R 0 is the average number of individuals that an infected person will infect over the entire period of being infectious. Thus, for a reproduction factor R 0 and a contact channel j with transmission probability p j ∈ [0, 1], we need to generate C(R 0 , p j ) := R 0 /(q I · p j ) many connections on average for all individuals in each time step. Conversely, for any process that generates η j connections to unique and distinct individuals over channel j in each time step, the effective R 0 over contact channel j with 100% transmission probability equals q I j η j . The actual number of contacts is drawn form a binomial distribution with n = S − 1 and a rate p = C(R 0 ,p j ) 2(S−1) . Note that the rate is one half of the target contact rate because all contacts are symmetrically mirrored.

In order to assess if the CRISP model is able to provide realistic population-level statistics for COVID-19 infection spread, we simulate a population of |S| = 10, 000 individuals over a period of 274 days (9 months). We single out an individual u for whom we set p 0 = 1 so that she will get infected with probability 100% at t = 1 ("patient 0"); for all other people we assume a p 0 = 10 −6 to model a miniscule chance of infection spread from exogenous sources. We assume a single contact channel with a 1% chance of transmission, p 1 = 0.01. We simulate five scenarios:

• No Mitigation. Since R 0 of COVID-19 is estimated to be 2.5, at any time t we generate C(2.5, p 1 ) random connections for every individual at every time step.

• Social Distancing After 60 Days. Intuitively, the "locality" of the contact patterns should play a role in the infection spread of COVID-19: if an individual is in contact with a broad range of other individuals, the spread should be faster than if unique number of people in contact over time is small. In order to demonstrate that this concept has indeed an effect on the infection spread, we performed an additional simulation where we kept the unique number of people that every individual meets in every time step at C(2.5, p 1 ) but introduced the concept of "social bubbles" where all individuals form groups of 20 who have a large number of interactions with each other (i.e., equivalent to C(2, p 1 ) but only rare interactions with people from other bubbles equivalent to C(0.5, p 1 ) (see Figure 3 for a picture of the contact matrix with random connections and with "social bubbles").

• Mitigation After 60 Days. For t ≤ 60, we generate C(2.5, p 1 ) random connections for every individual at every time step. Afterwards, we assume that mitigation measures are taken which reduce the reproduction rate to 1.0. Thus, we generate only C(1.0, p 1 ) random connections for every individual at every time step t > 60.

• Suppression After 60 Days. For t ≤ 60, we generate C(2.5, p 1 ) connections for every individual at every time step. Afterwards, we assume that lock-down measures are taken to suppress the pandemic which reduce the reproduction rate to 0.5.

• Suppression After 60 Days and Release of Lock-down after 120 Days. This scenario is similar to the previous scenario but we assume that due to very low infection numbers, the lock-down is released after 60 days. Thus, we generate C(2.5, p 1 ) random connections for every individual for t > 120. In Figure 2 , we show the plot of u P(z u,t = z) over t = 1, . . . , 274 days for z ∈ {E, I, R} (orange = E, red = I, blue = R) from 100 forward samples of the CRISP model for these scenarios. As one can see, with no mitigation there is a high peak around day t * = 180 and eventually herd-immunity is achieved at 85% of infected population. Even though the number of unique contacts in each time step is the same, "social bubbles" flatten the curve, thus slowing down the infection but growth rates of infected people are still super linear until large parts of the population had been in contact with the disease. Note that a similar mitigation policy is currently used in Belgium. In case of mitigation to R 0 = 1.0, growth rates are pushed to sub linear but the pandemic is still continuously going on after 9 months. Not surprisingly, suppression is most effective at bringing the infections back to nearly 0% after 120 days. However, if the lock-down is lifted after 120 days, a second wave of infections will cause an exponential increase in infectiousness after only two weeks (dashed lines). Note that all these effects were computable by simply forward sampling our individual-level CRISP model.

In order to assess the test and quarantine efficacy of the CRISP model, we consider a population of |S| = 1, 000 individuals for 150 days (5 months) with a uniformly random contact pattern of C(2.5, 0.025) = 5.03 contacts on average per individual and day. We simulate the actual infection spread by applying the following sequence in each time step (i.e., day): At the beginning of each time (bottom-right) Lock-down at day 60 and reduction of R 0 to 0.5 (solid lines). In dashed lines we show the effect of a subsequent re-opening of a subsequent contact rate increase to R 0 of 2.5 starting at day 120.

step, we query the testing-and-quarantining policy for a list of individuals which need to be tested and need to be in quarantine during this step (this will only be done after t * = 30 to simulate an undetected initial outbreak). Each policy is constrained to select no more than 10 test candidates per day (1% of the total population). Given the quarantined individuals on that day, we remove contacts from and to the quarantined individuals for that day and then use the CRISP forward model (4) and CRISP test outcome model (9) to draw one sample of the next simulated infection state of every individual as well as the actual test outcomes of the requested test candidates. If the infection state of an individual changes from E to I in this sampling step, we assume that with 50% probability, the individual generates symptoms. Finally, at the end of the time step, the testing-and-quarantining strategy is revealed the test outcomes as well as the list of symptomatic individuals (again, provided t ≥ t * ). We single out an individual u for whom we set p 0 = 1 so that she will get infected with probability 100% at t = 1 ("patient 0"); for all other people we assume a p 0 = 10 −4 to model a small chance of infection spread from exogenous sources.

1. Symptom-Based Policy. For every time step t ≥ t * , we will request testing for up to 10 symptomatic individuals from the previous time step. For all individuals with a positive test outcome on the previous day, we will institute a quarantine for ρ time steps where ρ ranges from 2 to 21 days in our evaluation.

2. Contact-Tracing Policy. For every time step t ≥ t * , we will request testing for up to 10 symptomatic individuals from the previous time step. If there are less than 10 symptomatic individuals, then we will request the remaining tests for individuals in quarantine sorted in descending order of the number of contacts they have had in the past 7 days with people who have tested positive. For every individual with a positive test outcome, we will not only quarantine her but also all the contacts she had in the past 7 days for ρ time steps where ρ ranges from 2 to 21 days in our evaluation; for every individual with a negative test outcome, we will remove her from quarantine.

3. CRISP Model-Based Policy. For every time step t ≥ t * , we will use block-Gibbs sampling of 100 infection traces z u to estimate P(z u,t ) for every individual u at the current time step t based on the contacts and test outcomes prior to time step t. We will request testing for up to 10 symptomatic individuals from the previous time step. If there are less than 10 symptomatic individuals, then we will request the remaining tests for individuals (who have not tested positive before) in descending order ofP(z u,t = I). We will quarantine any individual who is not yet quarantined but whose estimated probabilityP(z u,t ∈ {E, I}) exceeds a given policy threshold τ EI ; we will release an individual from quarantine once their estimated probabilityP(z u,t ∈ {S, R}) exceeds a given policy threshold τ SR . Note that we increase p 0 in the block-Gibbs sampling by a factor of 10 to account for "patient 0".

In order to gauge the efficacy of each policy, we measure two quantities at the end of the simulation (t = 150): (1) Percentage of population that got infected during the 150 days, and (2) total number of days that individuals were quarantined (e.g., if a policy locks down for the entire 150 days, this would result in 150,000 quarantine days). Varying the policy parameters ρ, τ EI and τ SR results in curves on the two dimensions of infection percentage and quarantine days. The closer a curve is to the origin, the more effective is the policy in terms of "health" (infection) and "economic" (quarantining) cost.

In Figure 4 , we plot curves for the three policies with ρ ∈ {2, 7, 14, 21}, τ EI ∈ {0.2, 0.3, 0.4, 0.5}, τ SR = 0.9. For comparison, we also show the two extreme points corresponding to "no mitigation" (i.e., zero quarantine days but the largest infection percentage of 90%) and "full lock-down" (i.e., largest quarantine days of 120,000 and near-zero infection percentage). All three curves exhibit a negative slope where a higher percentage of quarantine days corresponds to a more effective mitigation of infection spread. Of the three policies, our CRISP-based policy achieves the best performance in terms of the smallest number of quarantine days for a given infection percentage (i.e., Pareto frontier). This is because our CRISP model is able to accurately identify infectious users (even though they may be asymptomatic) and test/quarantine them proactively-this helps to prevent infection from spreading across the population while at the same time quarantining fewer individuals with a high likelihood of getting infected. In contrast, the symptom-based policy only tests individuals with symptoms and then quarantines the individuals who have tested positive. As a result, since 50% of the infected individuals are asymptomatic, they never get tested and quarantined, thus resulting in a spread of infection to 60% of the population. Similarly, the contact-tracing policy, by isolating all contacts of positive tested individuals (many of whom may have low likelihoods of getting infected), is able to achieve the absolute smallest infection percentage but at the cost of massive quarantining (30% of the population). Figure 5 shows a visualization of these effects in one of the simulation runs for ρ = 14, τ EI = 0.3, and τ SR = 0.9.

In this paper, we proposed a probabilistic graphical model for COVID-19 infection spread through individual contacts that captures the epidemiological state of each individual based on the SEIR model. We developed a computationally efficient block-Gibbs sampling-based algorithm to infer In the bottom plots, we show a stacked bar chart of quarantined individuals per day grouped by actual infection status. While the number of quarantined individuals for the symptom-based policy is small, the infection spread is not contained and the quarantining keeps growing exponentially. In contrast, the contact-tracing policy effectively suppresses infection spread while regularly quarantining more than 25% of the population. The CRISP model-based policy is initially picking a large number of individuals for quarantining but is then able to keep it at a low-level, in particular of susceptible individuals.

the COVID-19 infection risk score of all individuals at any time, given test outcome and mutual contact information between individuals. An efficient C++-based Python implementation of our inference algorithm is available at https://github.com/zalandoresearch/CRISP. Through experiments with simulated data, we showed that the CRISP model is able to model macro-level characteristics of the COVID-19 infection at county level (≈ 10, 000 individuals) and effectively mitigate COVID-19 spread by pro-actively quarantining and testing individuals with high risk of infections.

As part of future work, we would like to further accelerate our inference procedure using other approximation techniques such as Variational Bayes and Expectation Propagation [2] . Our inference algorithm can also be speeded up by exploiting the parallelism inherent in our block-Gibbs Sampling algorithm. For example, it is possible to concurrently sample infection traces of two individuals with no contacts in common. It is also known that the hyper-parameters of the SEIR model vary with demographic attributes such as age, socio-economic status, or location (see, for example [10] who present a location-varying infection spread model). We would like to extend our model with group-level hyper-parameters to account for this variation. We would also like to explore the causal impact of mitigation or suppression policy measures (e.g., school closures, shop closures, small group gatherings) on COVID-19 infection spread when using contact-level data. Finally, we would like to consider more sophisticated models of COVID-19 transmission through different modalities, and contacts with varying duration and distance characteristics.

Incubation period of 2019 novel coronavirus (2019-nCoV) infections among travellers from Wuhan, China

Pattern Recognition and Machine Learning

Epidemic thresholds in real networks

Inference for individual-level models of infectious diseases in large populations

Strategies for mitigating an influence pandemic

Strategies for containing an emerging influenza pandemic in southeast asia

Learning influence probabilities in social networks

Maximizing the spread of influence through a social network

Influence maximization on social graphs

A spatiotemporal epidemic model to quantify the effects of contact tracing, testing, and containment

Aristides Gionis, and Antti Ukkonen. Sparsification of influence networks

Infectious diseases of humans: dynamics and control

Communication-efficient learning of deep networks from decentralized data

-homogeneous virus spread in networks

Virus spread in networks

On the convexity of latent social network inference

Individual-level modeling of infectious disease data: Epiilm

Clinical presentation and virological assessment of hospitalized cases of coronavirus disease 2019 in a travel-associated transmission cluster

Acknowledgments We would like to thank Sebastian Munoz, Christopher Gandrud, Jasvinder Kandola and Peter Herbrich for all their valuable input and feedback on earlier drafts of this paper. We are also indebted to Christoph Thöns for his support with Python and C++ programming.