key: cord-0573432-ozcrolrg
authors: Caccioli, Fabio; Martino, Daniele De
title: Epidemic oscillations induced by social network control: the discontinuous case
date: 2020-12-04
journal: nan
DOI: nan
sha: 58dfb2a9b5509cd3290590c9cc60f865e7753669
doc_id: 573432
cord_uid: ozcrolrg

Epidemic spreading can be suppressed by the introduction of containment measures such as social distancing and lock downs. Yet, when such measures are relaxed, new epidemic waves and infection cycles may occur. Here we explore this issue in compartmentalized epidemic models on graphs in presence of a feedback between the infection state of the population and the structure of its social network for the case of discontinuous control. We show that in random graphs the effect of containment measures is simply captured by a renormalization of the effective infection rate that accounts for the change in the branching ratio of the network. In our simple setting, a piece-wise mean-field approximations can be used to derive analytical formulae for the number of epidemic waves and their length. A variant of the model with imperfect information is used to model data of the recent covid-19 epidemics in the Basque Country and Lombardy, where we estimate the extent of social network disruption during lock downs and characterize the dynamical attractors.

The onset of oscillations in a system as a consequence of feedback has been highlighted since the inception of control theory [1, 2] . Traditionally, feedback-induced oscillations have been studied in engineering artificial tools like thermostats and steering devices [3] . More recently, research focused also on its application to natural systems [4] , in particular homeostasis and its disruption in biological systems, a classical example being glycemic control and diabetes in human metabolism [5] .

Feedback-induced oscillations are currently emerging as governments are trying to control the evolution of the covid-19 pandemic crisis with containment measures such as social distancing, lock downs and quarantine.

The modeling of containement measures in compartmentalized epidemic models [6] is thus under the focus of intense research [7, 8] . It has been very recently rigorously demonstrated that compartimentalized epidemic models display oscillations in presence of feedback between infection rate and infection states [9] ,and that in general a feedback between order and control parameters in large interacting systems subject to phase transitions triggers self oscillations [10, 11] , where an Andronov-Hopf bifurcation takes over the usual phase transition.

As infection and recovery rates are changed, epidemic models on networks display out-of-equilibrium phase transitions between a phase where a disease is prevented from spreading and a phase where a macroscopic finite fraction of the population becomes infected [12] .

In this article, we will study the SIS and SIR models in a full microscopic settings on random networks in presence of a feedback that changes the structure of the underlying social network, and we will show that such feedback triggers self-oscillations along the theory proposed in [10] , where suitably defined connectivity properties play the role of the control parameter. In order to mimic the occurrence of lock downs, we will focus on a a simple discontinuous feedback control, where a certain fraction of links is deleted if the fraction of infections exceeds a given threshold 2. The same links are then be reinstated once the fraction of infections has been reduced below a second threshold value 1 < 2.

Oscillations in epidemic spreading have been studied mainly from the point of view of seasonal effects that act as an external driving forces, while in the case studied here oscillations are autonomously driven by an internal feedback. The resulting models are described by time-independent equations and parameters, and such oscillations can be considered as emerging self-oscillations [13, 14, 15] .

The article is organized as follows: In the first section we define the model and illustrate its behavior with results from numerical simulations on an instance of a real social network. In the second section we then study the model on mean-field uncorrelated networks, where we will show that the overall effect of lock downs on the dynamics is captured by a renormalization of the effective infection rate through a change of the network branching ratio.

This finding is then exploited in the third section, where we analyze simple piece-wise well-mixed models with point transformation techniques, leading to analytical formulae for the number of waves and their length. In the last section, we consider the realistic case of imperfect information on the infection state, and we infer parameters from data on the current evolution of the covid-19 pandemic, for which we estimate the extent of social network disruption and characterize its limit cycle attractors.

We consider compartmentalized epidemic models on random networks, specifically the SIS and SIR models (see [16] for a review), where individual agents are represented as the nodes of a social network and can be in different states: Infected, Susceptible and Recovered. Infected individuals recover with a Poissonian rate , which is a parameter of the model, becoming either susceptible (in the SIS model) or recovered (in the SIR model), and they infect neighboring susceptible nodes with a Poissonian rate , which is the second parameter of the model. In order to model containment measures and their relaxation, we consider a feedback between the network structure and the infected state and its history in the following terms:

• Starting from a state with few infections (whose relative number we will indicate with ) that fastly spread, if the spreading overcomes a certain threshold > 2 a central authority decides to disrupt the network structure by randomly removing a macroscopic fraction of the links. This will eventually revert back the spreading.

• Starting from a regressing infected state in a disrupted network, when the infection state is reverted to an acceptably low value < 1, the network structure is restored back to its initial conditions. → : In a dense social network few infected highly contagious people give rise to an epidemic spreading.

→ : during the epidemic outbreak a centralized authority decides for containment measures by severing the network.

→ Under confinement, the epidemic regresses to few cases.

→ Once the epidemic is supposedly under control containment measures are withdrawn and the social network is restored. Red dots represent infected individuals, while white dots represent susceptible ones. Light grey links in panels C and D represent links that have been removed because of the containment measure.

For the SIS model, this will eventually lead to an infection cycle as il-lustrated in Fig. 1 , where we show results of simulations on a school friendship network reconstructed in [17] (number of nodes = 134). In Fig. 2 , we show instead a simulation of the SIR model on the same network with and without feedback control, thus illustrating the effect of enforcing containment measures. The control successfully reduces the spreading of the infection, but for this to occur a series of lock downs have to be put in place. We will consider in the next section the SIS and SIR models for the case of large uncorrelated random networks, where it is possible to characterize analytically general features of the dynamics.

We consider here the case of large annealed uncorrelated random graphs with degree distribution ( ). At odds with static networks, in annealed networks we assume that links are randomly rewired over a faster time scale of the spreading process, while the assumption of uncorrelated networks implies there are no correlations between the degrees of neighboring nodes. The locally treee-like structure of these networks makes it possible to make analytical progress in the study of dynamical processes taking place on them, since it allows to recur to well-controlled approximations for the factorization of probability states. In particular, we will consider here the heterogeneous mean-field approximation, where nodes are grouped in classes according to their degree.

Let us start from the SIS model. In absence of feedback, the rate equation for the fraction of infected individuals of degree can be written as follows [18] 

where Θ( ) = ∑︀ ( ) ⟨ ⟩ ( ) is the probability that a randomly selected neighbor of a node of degree is infected, and we denote by ⟨ ⟩ the average degree of the network.

Here we consider the case in which, when the fraction of infected individuals ( ) = ∑︀ ( ) ( ) exceeds a given threshold 2, a lock down measure is implemented that removes a fraction of links, which are then reinstated once the condition ( ) < 1 is satisfied. The equations of the model in presence of this feedback mechanism can therefore be written in terms of a state-dependent infection rate as followṡ

wherẽ

We note that we can express˜as a function of the fraction of infected population and its derivative because we are considering deterministic rate equations. A more general representation for the discrete stochastic case would require the introduction of a binary state variable to denote the occurrence or absence of a lock down.

In Figure 3 , we compare the result of numerical simulations with the numerical solution of equation (2) for Erdős-Rényi and scale-free random networks. In both cases we clearly see the emergence of oscillations due to the feedback. In Figure 4 , we visualize the feedback-induced oscillations by means of a phase portrait, where we plot the fraction of infected individuals vs. the fraction of new positives. The figure clearly shows the emergence of a limit cycle as an attractor of the dynamic. We also note from both Figures 3 and 4 that the dynamics on scale-free networks display bigger sample-to-sample fluctuations than that on Erdős-Rényi networks. The same feedback mechanism can be considered for the SIR model as well. If we now define as the fraction of recovered nodes with degree , the rate equations that describe the dynamic of the SIR model are

where, as before, Θ( ) = ∑︀ ( ) ⟨ ⟩ ( ) 1 and˜( ,˙) is given as before by equation (3).

In figure 5 we show the evolution over time of the fraction of infected and recovered individuals for the case of Erdős-Rényi and scale-free networks. We see that the introduction of the feedback can lead to oscillations corresponding to multiple infection waves. Clearly, in contrast to the case of the SIS model, these oscillations will eventually come to an end once a large enough fraction of the population has been infected. The number of infection waves depends on the parameters of the model. In the next section we provide an analytical estimation for a well-mixed population in the limit when is close to 1.

The analysis reported in the previous section, where the effect of network embedding and the feedback enacting on it is captured by a simple renormalization of the bare infection rate, suggests to study the model on well-mixed populations to get analytical insights. In this section , and will indicate the fraction of infected, susceptible and recovered individuals in the total population, respectively. The feedback law mimicking containment, with parameters 1 > 0 and 2 > 1, is given bỹ

For this form of the function˜( ,˙), mixed epidemic models can be readily analytically solved piece-wise in each sector, where the infection rate is constant, and the solutions can be joined at the boundaries (see Andronov [2] ).

For instance, in the mixed approximation the SIS model with no feedback and infection rate has the simple solution ( + = 1)

In presence of the feedback, the piece-wise constructed solution shows that for 1 > > 0, 2 < 1 − / 1 the dynamics settles into a limit cycle, and the periods of the quiescent epidemic spreading ( 1) and of the recovery under lock downs ( 2) are given by the following analytical formulae

For a swift and resolute population lock down, we can approximate 2 << 1 − / 1, 0 << and obtain for the total duration of an epidemic wave

where 0 = 1/ . For instance, from the values 0 ∼ 3, , 2/ 1 ∼ 100 and 1/ ∼ 2 weeks we can calculate ∼ 3 months. For the SIR model, the solution in each interval -starting from initial conditions , , at time -reads

The total number of lock downs can be worked out analytically by joining solutions piece-wisely, and a first order expansion in ( 0/ 1, 0/ ) (see the appendix) gives the formula

where

and

For instance for the values 0 ∼ 3, 2 − 1 ∼ 0.1 ± 0.05 we get * ∼ 3 ± 2.

In the next section we will illustrate data modeling applications of our framework.

In this section, we illustrate our framework in the context of modeling epidemic data of the covid-19 infection in 2020 in the Italian region of Lombardy and the Spanish region of the Basque Country. Data includes daily reports of new infections and active cases 2 , plus a single prevalence estimate. For the purpose of data analysis, we consider a variant of the model with imperfect information by splitting the total number of infections into detected and undetected cases, whose numbers we denote by and respectively. We assume the existence of a detection process by which undetected infected individuals are spotted with rate and then put in isolation (which is equivalent to removing them). The mean-field rate equations are as followṡ

where˜( ,˙) is the piece-wise constant function defined in equation (6) (parametrized by 1, 0, 1, 2), and analogous equations for and/or , depending on whether we consider the SIS and/or the SIR model. The inference of the model parameters has been performed by approximately solving Bayes equations under the hypothesis of Gaussian noise by means of Monte Carlo methods and dynamical system numerical simulations (see the appendix). In both panels, the solid blue line refers to the model inferred maximizing the likelihood, red dots refer to daily data and green stars refer to weekly data.

In Fig. 6 we show the time series of new daily cases for the two analyzed cases against the model with maximum likelihood parameters, the latter showing a clear cusp peak and piece-wise exponential trends corresponding to the lock down event. The model suggests that a closed trajectory must be observed in the plane of new and active daily cases. This is shown in Fig. 7 , where we report data smoothened by a 7-days moving average as well. This closed trajectory in the ideal case of a pure SIS model with feedback would be a limit cycle attractor of the dynamics.

Connectivity plays a crucial role in the definition of the parameters that control the collective behavior of a system. This finding has striking consequences, like the absence of an epidemic threshold in epidemic spreading models defined on scale-free networks [19] , and it suggests that it is possible to control the spreading by acting on the network of social interactions.

In this article we have shown that feedback control at the level of the social network in epidemic models triggers self-oscillations along the theory proposed in [10] .

We have investigated self-oscillations induced by a simple discontinuous feedback control mimicking lock down events in classical compartmentalized epidemic models (SIS and SIR) on networks.

On random graphs, for Erdős-Rényi as well as scale-free networks with naive populations, we have shown that the effect of lock downs simply amounts at renormalizing the effective infection rate to account for the reduction in the network branching ratio.

This led to simple piece-wise mean field approximations that we solved analytically by means of transformation point methods, recovering formulae for the number of waves and their extent in terms of the model parameters. These formulae can be in principle tested against data, once a certain amount of evidence accumulates on the number of lock downs and their lengths, prevalence estimates, and basic infection numbers region by region.

A problem related to data collection during the covid-19 epidemic outbreak was the fact that many positives were undetected. In order to bring the model to data, we have therefore extended it by assuming the existence of a fraction of undetected positives, who can then be detected at a given rate (for instance through testing). We have applied our extended framework to analyze data from the first epidemic wave of covid-19 in Lombardy and Basque country, where parameters have been inferred leading to a characterization of the dynamical attractors in the phase space. Apart from applications to predictive modeling -which would require more extensive data analysis [20] and methods of system identification--we do point out here briefly some potentially interesting theoretical problems stemming from this work.

First, the issue of optimal scheduling [21] in the control of the social network, leading to continuous feedback and potentially smoother oscillations. Current qualitative evidence from the second epidemic wave of covid-19 seems indeed to show in some regions smoother trends and os-cillations around the phase transition point ( ∼ 1), due to attempts of finer control like partial restrictions and selected closures taken in due course. Recently proposed analytical frameworks [7] can be very useful in this respect.

Second, within the framework proposed here, where epidemics can be regarded as self-oscillators, it comes naturally the question of coupling and synchronization [22, 23] of epidemic waves running on different networks that are weakly connected, e.g. by migration processes.

Finally, another interesting issue concerns the impact of periodically external drive on oscillators: analogously to well-known forced double well oscillator [24] , the combined effect of seasonal changes and feedback could potentially lead to chaotic oscillations in strange attractors, an aspect that adds to the problem of predictability of such systems, and that we leave for future investigations.

In this section we derive formula (14) (15) (16) . In the SIR model at fixed infection rate (as in a simple model without feedback, or in a given interval for the piece-wise feedback model) the peak value of the infected fraction is (when˙= 0 and = / )

This value is not achieved if it is greater than the one triggering the lock down, i.e. when > 2. Thus we will assume as halting condition that ≤ 2, since in this case no lock down takes place and the system proceeds towards herd immunity. We will now work out a series for fraction of susceptible individuals at the various stages of the epidemic waves, exploiting the piece-wise analytical solutions in each interval.

Suppose we are at the beginning of a wave = 1 with given susceptible fraction = ,−, the system (with = 1) will evolve towards 2 and a given ,+ that satisfies

where Δ = 2 − 1. Then we have the lock down = 0, and the system will evolve towards 1 with a given +1,− that satisfies

These equations define a series eventually halting when ≤ 2.

Vanishing 0

Suppose 0 = 0. In this case +1,− = ,+ ≡ and we have

If we start from 0 ≈ 1, from the halting condition we find that the number * of lockdowns is

where we denote by ⌊ ⌋ the integer part of .

A first order expansion in 0 leads to

Then defining = ,− − / 1 log ,−

we have the recursion relation

that can be bounded by ( 0) terms (given that 0 ≤ ,− ≤ 1)

and the halting criterion leads to

In Figure 8 we show the agreement between numerical simulations on an Erdős-Rényi network and the analytical prediction given by equation (26) for 0 = 0. 

Here we consider the task of fitting against epidemic data a model including imperfect information on the state of the system. This is an instance of a system identification problem [25] , which we solved along the following lines: We consider the time series of observed new daily and active cases ( , ) ( = 1 . . . is the temporal index in days, starting from the 1st of March, = = 155 for Lombardy, and = = 100 for the Basque country), and we assume it as coming from an instance of the model plus a noise term

where we highlighted the dependence of the model trajectory by the dynamical parameters and boundary values. We assume shot-noise of the form

which for large numbers we assume to be distributed normally. Upon assuming an uniform prior, we have the following formula for the loglikelihood of the parameters

From the Bayes formula, the posterior probability distribution of parameters (℘) ∝ −ℒ(℘) has been sampled by a Metropolis Montecarlo rule, where the evaluation of ℒ(℘) has been done by numerical integration of the model equations. More explicitly we have been following the following flowchart:

• Start from some value of the parameters ℘0: a warm start has been provided by fitting the curve of new daily cases alone in linear approximation ( << ).

• Propose a change for the parameters ℘ → ℘ +1: we used independent geometrical random walks of stepsize 10 −3 .

• Numerically integrate the model equations with the new proposed parameters to evaluate their log-likelihood ℒ(℘ +1). We used the standard Verlet algorithm.

• Accept the proposed new parameters with probability min(1, exp ℒ(℘ +1) − ℒ(℘ )) (Metropolis rule), otherwise keep the old parameters.

This defines a series that asymptotically uniformly samples the posterior probability for the parameters, whose peak values have been used for the results showed in Fig. 6 and 7 . Finally, we do point out that when we numerically integrate the model equations we rejected solutions that i) do not include at least one lock down event and ii) do not agree with prevalence estimate from serological data. The latter gives a lower bound for the total number of infected individuals at a certain time, ≥ that can be reformulated as an inequality between the model parameters (in particular , 1 and 0) as follows:

can decomposed in detected and undetected cases = + , where = ∑︀ Δ + is given by the data, while for the latter we have 

In the first and second lines we have used the model hypothesis, in the third we have decomposed the sum in terms of the increasing and decreasing part of the wave, and finally in the last line is the fraction of the detected infections during the decreasing part of the epidemic wave. We have finally the inequality

Given the simplicity of the model employed -in particular with respect to the hypothesis of a constant detection rate -we obtain a fairly high value of the 2 ∼ 20, with a concomitant acceptably low average relative error of ∼ 20% across data points that can be considered apt for a qualitative description of the data. We report in the following table the maximum likelihood inferred parameters with their standard deviation.

Region 0 = 1/ 1 − = 0/ 1 (day −1 ) Basque country 2.9 ± 0.1 1.3 ± 0.4 · 10 −2 3 ± 1 · 10 −3 Lombardy 5.2 ± 0.2 4 ± 1 · 10 −2 8.4 ± 0.6 · 10 −3 We do point out an anomalously high 0 for Lombardy, due to a very low average recovery rate of approximately 1 month (while the one from Basque country is around two weeks in agreement with WHO estimates). This is apparent upon looking at the much slower decay of the infection curves during the lock down, and it is probably due to a biased oversampling of critically ill cases whose average recovery is typically longer.

Aleksandr Adol'fovich Vitt, and Semen Emmanuilovich Khaikin. Theory of Oscillators

Feedback systems: an introduction for scientists and engineers

Oscillation patterns in negative feedback loops

Metabolic regulation: a human perspective

Epidemics with containment measures

Covid-19 and flattening the curve: a feedback control perspective

Active control and sustained oscillations in actsis epidemic dynamics

Feedback-induced self oscillations in large interacting systems subjected to phase transitions

Oscillations in feedbackdriven systems: Thermodynamics and noise

Non-equilibrium critical phenomena and phase transitions into absorbing states

Self-sustained vibrations in volcanic areas extracted by independent component analysis: a review and new results

Experimental study of self-oscillations of the trachea-larynx tract by laser doppler vibrometry

Self-oscillation

Contact patterns in a high school: a comparison between data collected using wearable sensors, contact diaries and friendship surveys

Epidemic processes in complex networks

Epidemic spreading in scale-free networks

The evolving worldwide dynamic state of the covid-19 outbreak

Optimal control theory: an introduction. Courier Corporation

Coupled self-oscillating systems: Theory and applications

Filippo Miele, and Timoteo Carletti. Ginzburg-landau approximation for self-sustained oscillators weakly coupled on complex directed graphs

Nonlinear oscillations, dynamical systems, and bifurcations of vector fields

System identification: an introduction