Inferring the effect of interventions on COVID-19 transmission networks

Countries around the world implement nonpharmaceutical interventions (NPIs) to mitigate the spread of COVID-19. Design of efficient NPIs requires identification of the structure of the disease transmission network. We here identify the key parameters of the COVID-19 transmission network for time periods before, during, and after the application of strict NPIs for the first wave of COVID-19 infections in Germany combining Bayesian parameter inference with an agent-based epidemiological model. We assume a Watts–Strogatz small-world network which allows to distinguish contacts within clustered cliques and unclustered, random contacts in the population, which have been shown to be crucial in sustaining the epidemic. In contrast to other works, which use coarse-grained network structures from anonymized data, like cell phone data, we consider the contacts of individual agents explicitly. We show that NPIs drastically reduced random contacts in the transmission network, increased network clustering, and resulted in a previously unappreciated transition from an exponential to a constant regime of new cases. In this regime, the disease spreads like a wave with a finite wave speed that depends on the number of contacts in a nonlinear fashion, which we can predict by mean field theory.


February 26 -March 15
The early phase of the epidemic is characterized by a low number of cumulative infections. We can therefore directly use the absolute numbers of new infections as input for our agent-based model, as we are always far away from the epidemic threshold. We choose broad, uniform priors for all the parameters that can be found in table S1. Note that the probability for random connections varies on several orders of magnitude p ∈ [10 −6 , 10 0 ] and we therefore infer this parameter on a logarithmic scale. We use n ABC = 100 and obtain an effective sample size of n eff ≈ 46.

March 16 -June 6
In the time period from March 16 to June 6 Germany recorded 177652 cases in total. This means that it becomes computationally unfeasible to replicate the population directly in our model without noticing a strong effect of the removed (immune) agents. Therefore, we scale down the total number of infections to our system size and compare the relative number of cases per 300,000 people instead. Assuming our hypothesis that the NPIs lead to a strongly clustered transmission network holds, we expect a large number of unconnected communities in Germany in that time period, which we represent as distinct model instances. We use n ABC = 200 and obtain an effective sample size of n eff ≈ 105. Our priors for this period can be found in table S2.

June 7 -September 15
To infer the parameters of the system during this time period we first sample parameters from the posterior distribution which we obtained for the previous time period, and let the system evolve for 81 days (corresponding to the time period from March 16 to June 6). Next, we change the infection probability p I and assign a new transmission network based on a new set of parameters p, k. For these parameters, we choose the same prior distributions as for the previous time period, see Table S3. We use n ABC = 200 and obtain an effective sample size of n eff ≈ 164. 2 Parameter scan

SEIR model
To investigate the disease dynamics in the small-world network, we perform a parameter scan. We vary the network parameters p, k while keeping the rest of the parameters fixed at n = 10 5 , p I = 0.02, n E (0) = 0, n I (0) = 10. Our choice for p I during the parameter scan is motivated by reports of the COVID-19 individual-level secondary attack rate (SAR) in the household of 17 %. Inverting Eq 3 we obtain We vary the number of contacts from k = 2, 4, . . . , 24 and sample the probability for random contacts in eleven equally-spaced steps on the log scale from log 10 p = −5, . . . , 0. As initial condition, 10 random agents are set to the infectious state and the system is simulated until there are no more exposed and infectious agents. We repeat this process five times per parameter combination (p, k). As output we determine the peak of simultaneously infectious people and the cumulative infection curves

SIR model
To compare our analytical prediction of the wave speed in the highly clustered network (see below), we define a simplified SIR model. The difference to the SEIR model is that there is no exposed state, and that the waiting times for the transition from the infectious to the removed state are drawn from an exponential distribution with mean τ I = 10d. For the parameter scan, we initialze the system of n = 10 5 agents in a ring-like topology (no random links) with a single infectious agent and let it evolve for 365 time steps. We repeat this process 20 times per parameter k. We then record the cumulative infection curves see Fig. 5a. We calculate the linear growth rate from the cumulative infections as where we neglect the initial exponential growth by skipping t min = τ I = 10 time steps. We also determine the maximum time t max,k until the epidemic dies out for each parameter k so that after t max,k time steps, the cumulative number of infections did not increase in any realization of the system 3 Mathematical analysis

Epidemic threshold
We can calculate an upper bound for the number of contacts k c by demanding R 0 = 1, i. e. a single infectious person in a network of susceptible people will effect on average one other person. Using Eq 3 (Main Text) we can calculate the expected number of infections as For p I = 0.02 and τ I ≈ 10 we obtain k c ≈ 5.5.

Derivation of differential-equation approximation
In order to predict the linear growth of infections in the highly clustered small-world network, we consider a simplified variant of our original model, where we neglect the exposed state and assume that the progression times are distributed exponentially (SIR model, section 2.2 in Supplementary Information). Then we describe the state of an agent j in our model by a set of three Boolean stochastic variables s j (t), i j (t), r j (t) = 0, 1, where s j + i j + r j = 1 and s j = 1 indicates that the agent is susceptible, i j = 1 means he is infectious and if r j = 1 he is removed. We represent the event "agent j becomes infectious at time t" by the Boolean stochastic variable α j (t) with P (α j (t) = 1) where I j (t) := m∈Nj i m (t) is the number of infectious agents in the neighborhood N j of agent j. Similarly, the event "agent j is removed at time t" is given by the stochastic variable β j (t) with P (β j (t) = 1) = p R = 1/τ I .
During one time step the state of all agents changes as We want to calculate the expected value of the state variable under a mean-field approximation, i. e. we replace the expected value of a function f (X) of any random variable X by the function evaluated at the expected value of the random variable, f (X) ≈ f ( X ). In particular, this also means that we neglect any correlations between the state variables of neighboring nodes. We denote the expected values of the state variables as σ j (t) := s j (t) , ι j (t) := i j (t) , ρ j := r j (t) . For the expected number of infectious neighbors we obtain where N is the total number of nodes in the network and p is the probability of a random link (see Model definition).
Next, we also introduce a small time step length τ > 0 and continuous timet = τ t along with the transition rates κ I := p I /τ, κ R := p R /τ . In the following we only consider continuous time and drop the hat for better readability. For τ → 0 we can approximate the expected value of α j (t) by a Taylor approximation around p I = 0 as P (α j (t) = 1) = 1 − (1 − κ I τ ) Ij (t) ≈ κ I τ I j (t). (S13) We can further simplify our system by considering the regime N → ∞, and introducing the spatial step length ∆x > 0 so that N ∆x = L = const. We can then replace the state probabilities σ j (t), ι j (t), ρ j (t) by the probability densities σ(x = j∆x, t), ι(x = j∆x, t), ρ(x = j∆x, t). This allows us to approximate the expected number of infectious agents in the neighborhood by expanding the terms ι(x + m∆x, t) ≈ ι(x, t) + m∆x∂ x ι(x, t) + m 2 ∆x 2 /2∂ xx ι(x, t) to obtain wherek := (k/2 + 1)(k + 1)/12. Approximating all time-dependent functions by their Taylor approximation up to second order f (t + τ ) ≈ f (t) + τ f (t) + τ 2 /2f (t) and rearranging terms we obtain the following set of non-linear PDEs for the expected value of the agents' states We introduce the constant D k := ∆x 2k τ and can now identify two separate time scales: A fast time scale, with terms ∝ 1 τ which characterizes the disease progression, and a slow timescale with terms ∝ D k that describes the wave of infections in the network. Finally, with I(t) := 1/L L 0 ι(x, t) dx as the total number of infectious people, we obtain the following set of nonlocal PDEs In the regime τ → 0 we expect the the system to quickly reach a local steady state, so that the fast time scale terms vanish. We can use this to obtain an approximate value for the wave speed of infections.

Calculation of the infection wave speed in the highly clustered network
To calculate the wave speed of infections we consider the regime p = 0, so that we can neglect the nonlocal coupling by I(t) in Eq S18. As described above, if we let lim τ, ∆x → 0 with finite D k , we expect the system to always be in a local steady state, corresponding to the fast time scale terms being 0. Then, we arrive at the following equation for the density of infected people ∂ tt ι(x, t) = κ I kD k σ(x, t)∂ xx ι(x, t).
If we consider a point far away from the wave front, we have σ(x, t) ≈ 1 (everyone is susceptible), and Eq S21 reduces to the wave equation ∂ tt ι(x, t) = κ I kD k ∂ xx ι(x, t) with the wave speed c = √ κ I kD k ∝ k √ k. Figure S1. Dynamics in the agent-based SEIR model. Susceptible agents can become exposed, if they are linked to infectious agents. Exposed agents become infectious after the waiting time τE, and infectious agents are removed after τI . All nodes are updated simultaneously at every time step t.