Self-adapting infectious dynamics on random networks

Complex adaptive systems¹ often arise when interacting nodes (or agents) adapt their dynamics (or strategies) due to external influences. This occurs, for example, not only in socio-economic systems, such as epidemic and information spreading,^2,3 markets,⁴ evolution of languages,⁵ or evolutionary games with environmental feedback⁶ but also in automation problems including machine learning⁷ or control theory.⁸ In standard network models, adaptive dynamics is often implemented by creation, deletion, or rewiring of edges between pairs of nodes,^9,10 some aspects of which have been extended to higher-order networks.^11–13 From a more general perspective, adaptivity is characterized by a set of strategies that are selected either by the system itself or by individual agents in the system, each leading to potentially different dynamical rules and future evolution.

The adaption of strategies may be characterized by a change of parameters. This occurs, for example, in piecewise deterministic Markov processes (PDMPs),^14–17 where a specific parameter of the system evolves according to an additional, independent stochastic process. Assuming that the parameter space represents a possible set of strategies, this is equivalent to a system where some decision maker (randomly) changes the strategy. If this parameter crosses a bifurcation, the stability of the system changes abruptly, typically leading to an entirely different dynamical evolution. On the other hand, if parameter changes are influenced by the history of the system and occur with some temporal delay, usually delay equations¹⁸ are considered. Examples are biomedical models of cancer evolution,^19,20 population dynamics,²¹ or machine learning.²² It is reasonable to assume that in real world systems, both effects occur: delayed adaptation depending on the history of the state and piecewise deterministic evolution coupled to sudden strategy changes.¹⁰ 

Particularly interesting are systems with policy makers who adapt either the rules of the system or their own behavior according to some, usually complex, evaluation mechanism. For example, restrictions and lockdown laws have been enforced and lifted during the Sars-Cov-2 pandemics depending on the recent course of the pandemic^23–25 and are regularly changed, e.g., due to the emergence of new variants or available vaccinations. Therefore, it is reasonable to ask how basic epidemic models behave in the context of adaptive policies. This is related to recent models for the influence of individual risk perception on the epidemic spreading of Sars-Cov-2,²⁶ where direct feedback mechanisms are considered.

Epidemic modeling has a long history,³ the most famous examples being the SIS and SIR models.^27,28 In these compartmental models, the population is divided into susceptible (⁠ $S$⁠) and infected individuals (⁠ $I$⁠), which after recovery become either susceptible again or recover (⁠ $R$⁠), for SIS and SIR systems, respectively. These models are often considered on contact networks, where existing links determine possible infections between individuals (see Ref. 3, and references therein). Many extensions have been studied, from additional compartments,²⁹ to adaptive rewiring or (de)activation of edges,^30–35 adaptive force of infection,^36–38 and including vaccination^39–43 or awareness control of behavioral changes.⁴⁴ Depending on the specific model and parameters, the system state typically either approaches a stable disease-free or endemic equilibrium, oscillates, or shows more complex dynamics. Apparently, such adaptive models often become analytically intractable, even if the high-dimensional network dynamics is reduced to the level of differential equations by appropriate closure techniques.

In this paper, we consider a framework for history dependent strategy adaption with piecewise deterministic dynamics. This is applied to simple models of adaptive epidemics for which an analytical treatment is possible. In particular, we consider a threshold-based activation and deactivation of edges within a fixed contact network for the examples of SIS and SIR systems. We observe and describe a large parameter region with stable oscillations, as well as bifurcations toward stable disease-free and stable endemic states for the case of SIS. We derive analytic expressions for the period which are validated with numerical simulations for random uniform networks but are unsuitable to describe the dynamics in scale-free networks. In summary, our model provides one of the simplest settings to generate oscillations in epidemic dynamics under global non-Markovian coupling.

The paper is structured as follows. In Sec. II, we formalize a model of strategy adaption. Section III recapitulates SIS and SIR models on networks. In Sec. IV, we apply the adaptive framework to SIS and SIR epidemics on networks. A summary and outlook is presented in Sec. V.

There are many ways to include adaptivity into dynamical systems, e.g., by rewiring rules^10,30 on the network level or by state-dependent parameter changes³⁷ on the level of differential equations. In contrast, here we study a system with a strategy space $S$ such that the dynamics depends on the currently active strategy $S i∈ S$⁠. For example, in a system with infectious dynamics, the strategies $S i$ could refer to different combinations of reduced numbers of average contacts and/or increased measures against transmission. Alternatively, such a strategy could impose the emergence (or deletion) of additional compartments $X∈{ Q, E, I ~,…}$ into the model, e.g., enforcing quarantine, separately counting exposed but not infectious agents or introducing new variants of the disease.

It is natural to assume that the chosen strategy depends on some observable function $g:Γ→ R n$⁠, which evaluates the state $x∈Γ$ of the system, where $Γ$ denotes the state space. This is, in particular, relevant if the strategies are chosen by some decision maker(s) based on their evaluation metric. Note that in real social systems, often many different metrics are used to evaluate the same state, e.g., incidence, hospitalization rate, and vaccination rate are all possible options to evaluate the state of an epidemic. For the sake of simplicity, we consider a one-dimensional metric, $n=1$⁠, $g(x)∈ R$⁠, in the following.

We define an adaptive dynamical system with discrete strategies as follows. Let $Φ$ be a map from the strategy space $S$ to a flow on $Γ$⁠, i.e., $Φ S:=Φ(S)$ is a mapping $Φ S:T×Γ→Γ$ for all $S∈ S$ with $Φ S(0,x)=x$ and $Φ S[ t 2, Φ S( t 1,x)]= Φ S( t 1+ t 2,x)$ for all $t 1, t 2∈T$ and $x∈Γ$⁠. An adaptive dynamical system is then defined as the tuple $(T,Γ,Φ, S,A)$⁠, such that for each $S∈ S$ the triple $(T,Γ,Φ(S))$ is a dynamical system and where the strategy $S t$ evolves according to the adaption function $A$ as defined above.

We emphasize that the time evolution of the system depends on the current strategy $S t$ and simultaneously the current strategy depends on the time evolution of the system. This leads to a nontrivial feedback mechanism and possibly complex and interesting dynamics. In order to introduce such an adaptive mechanism to some dynamical system, it is, therefore, necessary to specify a set of strategies $S$ and their implications on the internal dynamics, some observable $g$ and strategy function $J$⁠, and how the system adapts its dynamics to changes in $J$ through the function $A$⁠.

This framework connects piecewise deterministic Markov processes and delay equations: If the switching of strategies $A$ is governed by a Markov process the full system is equivalent to PDMPs (see, e.g., Refs. 14–17). On the other hand, if the selected strategy depends on some parameter evaluated at a delayed time $τ$ in the past, i.e., the integration kernel is $ρ(t)=δ(t−τ)$⁠, the resulting system can be reduced to delay equations.⁴⁵ In the following, we apply this adaptive framework to epidemiological models on networks.

We consider an epidemic system on a network $G=( N, E)$ with transmission rate $τ$ and recovery rate $γ$⁠. The basic network structure of $G$ is assumed to be fixed, i.e., the set of edges $E$ does not change. Adaptive activation and deactivation of edges occurs within this framework by assigning a weight $w: E→{0,1}$ to each edge $e∈ E$⁠, which determines if the edge is active, $w(e)=1$⁠, or not $w(e)=0$⁠. More generally, one could assume transmission rates $τ(e)=w(e)τ$⁠, which are relevant in the context of different preventive measures, such as quarantine^48,49 or social distancing.⁵⁰ In order to make this system adaptive, we consider the simplest case with only two strategies, $S={ S 0, S −}$⁠, which determine the number of active edges. Here, the null-strategy $S 0$ corresponds to the case where all (existing) edges in the network are active, $w(e)=1∀e∈ E$⁠. The strategy $S −$ corresponds to a lockdown strategy, where each edge is deactivated with probability $p cut$⁠, i.e., $P[w(e)↦0]= p cut$⁠. This means that the average degree in the adaptive network becomes time dependent with $k −= k 0(1− p cut)$⁠.

For the adaptive mechanism $A$⁠, we consider prevalence and incidence as suitable observable functions, i.e., $g 1(t, x I)= x I$ and $g 2(t, x I)=τ n S I/N≈τk(t) x I x S$⁠. Furthermore, the kernel in Eq. (1) is chosen such that $g$ is averaged over fixed time intervals $[t−Δt,t]$ for some time span $Δt$⁠, i.e., $ρ(t)= 1 Δ t[Θ(t)−Θ(t+Δt)]$ with Heaviside function $Θ$⁠. This leads to $J 1 , 2(t)= ∫ t − Δ t t g 1 , 2[ t ′, x I( t ′)] d t ′$⁠. We emphasize that such an average has been a commonly used measure during the Corona pandemic, e.g., in terms of the so-called $7$-day incidence (see, e.g., Ref. 23). Let us emphasize that $J$ measures the seriousness of a pandemic situation, where more infections or larger prevalence manifest in larger values of $J$⁠.

Note that there are numerous possibilities for defining the adaption function, which can either deterministically or also stochastically select the future strategy. In the latter case, the system contains two distinct sources of stochasticity, one from the dynamics and one from the adaption mechanism, which adds another layer of complexity. Here, we focus on the deterministic case.

We emphasize that our proposed model differs quite substantially from recent stochastic SIRS models,^15–17 which are piecewise deterministic Markov processes (PDMPs). In these systems, the switching of parameters is guided by an additional random stochastic process, which is independent of the (deterministic or stochastic) dynamics. Here, in contrast, the dynamical evolution of the system feeds back into the adaptive mechanism. If the underlying dynamics is stochastic, the adaption of strategies is guided by the (stochastic) feedback process implemented with the adaptive function. On the other hand, in a purely deterministic setting, the proposed system can be seen as an extension of the phase-space with the set of strategies, $X× S$⁠, where switching between different branches occurs deterministically according to the adaptive function $A$⁠.

For the SIS epidemics on a network, one expects an initial growth of the prevalence $x I$ until $J$ exceeds $ξ +$⁠. At this point, the social distancing strategy $S −$ is applied by removing a proportion of $p cut$ edges from the system [see illustration in Fig. 2(a)]. If sufficiently many edges are removed, the prevalence decreases and $J$ becomes smaller than $ξ −$ such that the null-strategy $S 0$ is reapplied. Within this regime, we expect periodic behavior.

In the following, we first derive different dynamical regimes in the adaptive pairwise closed SIS model and second compare this to numerical results obtained from simulation on networks. Recall that for each strategy, $S 0$ and $S −$⁠, there is one stable equilibrium $x ∗(β)$ with $x ∗(β)= x 0=0$ for $β=kτ/γ<1$ and $x ∗(β)= x en(β)=1− β − 1$ for $β>1$ [see the red curve in Fig. 2(a)]. In the following, we assume arbitrary but fixed thresholds $0< ξ −< ξ +<1$ and specify how the dynamics depends on the choice of $( β 0, β −)$⁠. For an illustration, see Fig. 2(b). For simplicity, we also assume that $Δt=0$⁠, such that $J(t)=g( x I(t))$⁠. First, if $β 0<1$ the state of the system converges to the stable equilibrium $x 0=0$ and the strategy remains in $S 0$ for all times. Second, for $β 0>1$ the state of the system approaches the endemic state $x en( β 0)$⁠. The strategy switches, if $g( x I)≥ ξ +$⁠, which implies the limit $x en( β 0)=1− β 0 − 1= g − 1( ξ +)$⁠. In particular, for $1< β 0< 1 1 − g − 1 ( ξ + )$⁠, the threshold $ξ +$ is too large and the state of the system converges to $x en( β 0)$⁠. On the other hand, for $β 0> 1 1 − g − 1 ( ξ + )$⁠, the strategy will eventually switch to $S −$ for some $t>0$⁠. After switching, the state converges to the stable equilibrium $x ∗( β −)$⁠. Similar considerations as above lead to the following limits depending on the threshold $ξ −$⁠. If $β −> 1 1 − g − 1 ( ξ − )$⁠, the state approaches the stable endemic equilibrium $x en( β −)$⁠, without reaching the lower threshold, thus remaining in strategy $S −$ for all times. If $β −< 1 1 − g − 1 ( ξ − )$⁠, the stable equilibrium is below the threshold such that after some finite time, the system switches back to $S 0$ again. Hence, in this regime, the dynamics is periodic. Conversely, if $β 0>1$ and $β −$ are fixed, it is possible to specify limits for the thresholds $ξ ±$⁠, within which we expect periodic dynamics. In particular, one obtains $ξ +≤g(1−1/ β 0)$ and $ξ −≥g( x ∗( β −))$⁠. With $β −= ( 1 − p cut ) k 0 τ γ$⁠, these conditions similarly imply a lower bound for the cutting probability $p cut$ at given $ξ −$⁠, given by $p cut≥1− γ [ 1 − g − 1 ( ξ − ) ] k 0 τ$⁠. In Fig. 3, we illustrate the time dependence of the relative prevalence $x I$ for the adaptive SIS epidemics with strategy function $J 2$ on two different network types, networks from the random Erdös–Renyi (ER or random uniform) ensemble $G(N,p)$^51,52 and scale-free Barabási–Albert (BA) networks.⁵³ We emphasize that these network types are chosen as commonly known examples to illustrate our approach. The parameters are $γ=0.25$⁠, $k 0=50$⁠, $β 0=2$⁠, and $p cut=0.8$ with thresholds $ξ +=0.025$ and $ξ −=0.005$⁠, i.e., the lockdown strategy is enforced when each infected individual infects on average $2.5%$ of the population (over the past time-window $Δt=5$⁠) and it ends below $0.5%$⁠. These parameters correspond to the periodic regime in Fig. 2(b).

For both network types, the fraction of infected individuals $x I$ oscillates see top panels in Figs. 3(a) and 3(b). The corresponding incidence function also oscillates between the predefined thresholds $ξ ±$⁠, shown as the blue curve and red dotted lines in the bottom panels. Note that the discontinuities of the incidence function are caused by the sudden (de)activation of edges, which immediately changes the number of $S I$-edges and thereby the possible number of new infections. For comparison the analytic result of the pairwise closed moment system is shown as a dashed black curve, see the  Appendix. This agrees qualitatively very well with the simulation for the ER network over the entire scale of our simulation runs. In contrast, the initial spreading of the epidemics in the BA network is much steeper, showing also larger maximal values of $x I$⁠. Consequently, the frequency of fluctuation for these scale-free networks is increased compared to the pairwise closed system.

In Fig. 4, we show $T$ vs the upper threshold $ξ +$ for two different values of $β 0∈{1.6,2}$ and compare network simulations (green errorbars) to the analytic result (black curve). The upper threshold is chosen in the interval $ξ −≤ ξ +< x en( β 0)$⁠. Note that the expected period becomes zero at some $ξ +$ smaller than $ξ −$ (see Fig. 9), e.g., for $Δt=0$⁠, this occurs at $ξ += ξ −$⁠. On the other hand, if $ξ +→ x en( β 0)$ the expected period goes to infinity for the analytic model. This is also observed in the network simulations: for random ER networks, the analytic expression of the pairwise closed model fits very well for both $β 0$ [see (a) and (c)]. For the scale-free BA network and $β 0=1.6$⁠, the observed period is systematically too large for smaller $ξ +$ and too small for larger $ξ +$⁠.

Intuitively, this is explained as follows. The prevalence increases initially much faster in scale-free networks due to nodes with very large degrees (superspreaders), thereby overshooting the threshold significantly stronger [see Fig. 3(b)]. Thus, the initial prevalence in the second half of the cycle is increased compared to the random networks, thereby increasing the period significantly. This effect becomes negligible if the upper threshold gets closer to the endemic equilibrium, since here the prevalence typically flattens for all network types and the overshooting effect vanishes. Instead, the initially much faster speed of infection leads to shorter periods in this regime.

The strong dependence on the topology is revealed when considering larger values of $β 0=2$ and leaving the other parameters unchanged. Note that with $p cut=0.6$⁠, this still implies $β −=0.8$ significantly smaller than one. In the scale-free network, the observed periods are up to four times longer than predicted by the moment equations [Fig. 4(d)]. This means that with the same number of deleted edges in scale-free networks, the progression of the epidemic spreading is slowed down more efficiently than in random networks. Note that these longer periods are mainly caused by a much flatter slope during the lockdown phase $S −$⁠. This corresponds to an effective cutting probability of $p cut≈0.49$ (gray line).

This shows that the network topology leads to observable differences in adaptive network systems. We believe that such topology induced effects on oscillations can also be observed in more complex models.

The basic reproductive number $R 0$ is commonly defined as the expected number of new infections caused by one infected individual in a fully susceptible community.^46,56 For the pairwise closed SIS model, this gives the ratio of the total rate of infection to the total rate of recovery, $R 0=β=kτ/γ$⁠.⁴⁶ If for infectious disease $R 0>1$⁠, an epidemic outbreak occurs, then for $R 0<1$⁠, the disease vanishes.

In Fig. 5, we illustrate the reproductive ratio $R 0$ as a function of time, comparing numerical values (blue) to the simplest estimation based on $β$ (black). The numerically observed $R 0$ oscillates in a similar pattern around $R 0≈1$⁠, however, fluctuating on faster scales and typically below $R 0= β 0=2$⁠. We assume that the latter is caused by infectious clusters with a small number of susceptibles, leading to much fewer infections for large prevalence.

We emphasize that such fluctuations of the reproductive number around $R 0≈1$ are also observed in many different countries during the Corona pandemic, e.g., for evolution in Germany, see Fig. 6. This indicates that in real epidemics there is some form of self-organized criticality around $R 0=1$⁠, i.e., there is intrinsic adaptive self-tuning of the system toward the critical bifurcation value. This can be explained due to the conflicting goals of reducing the number of infections and increasing the freedom of individuals. In particular, our simple model with adaptive strategies already mimics the complex adaptive dynamics due to stricter lockdown measures and more general cautious behavior at higher incidence levels. It is, therefore, reasonable to expect similar observations in more sophisticated models, e.g., when each agent chooses from a set of strategies and adapts in a more complex manner to the development of an epidemic. This supports that epidemic models with self-limiting feedback mechanisms should be viewed through the lens of self-organized criticality.⁵⁷ 

We further apply the proposed adaptive strategy approach to the SIR epidemics on networks.⁴⁶ In the SIR model recovered individuals do not contribute to the epidemic spreading anymore. Similar to the SIS model there is a regime where initial transient oscillations are expected. When a significant proportion of individuals is recovered, $R$⁠, the epidemic dies out after a certain number of cycles.

This is illustrated in Fig. 7, showing the time dependence of relative prevalence and incidence for the SIR dynamics on two different network types, similar to Fig. 3. The considered thresholds are $ξ +=0.01$ and $ξ −=0.002$⁠. Note that for real epidemics, these thresholds should be chosen according to the capacity of the health system, which typically can sustain only a very small fraction of infected individuals. Furthermore, choosing larger thresholds in the simulations leads to many infected and recovered individuals already before strategy adaption takes place. For both network types, we observe the expected initial oscillatory dynamics in Fig. 7, similar to the SIS system. We observe that Erdös–Rényi networks show comparable results to the pairwise closed solution (black dashed line), even though individual simulations differ substantially. In contrast, the epidemic spreading in the scale-free network takes place much faster. This is intuitively explained with the power law degree distribution and the existence of a small number of highly connected nodes. First, if one of these nodes is infected, it triggers a lot of subsequent infections. Second, after its recovery, it efficiently blocks the epidemic from spreading.

We investigate numerically how the times of strategy adaption depend on the probability to cut edges during the lockdown phase $p cut( S −)$⁠. This is illustrated in Fig. 8 for random ER networks and compared to numerical solutions of the first order moment SIR system with strategy adaption (black lines). For both network types and small $p cut≤0.5$⁠, the strategies switch only twice, from $S 0$ to $S −$ (violet colors) and back (green colors). In this regime, the corresponding $β −≥1$ such that the epidemic spreading slows down but does not revert. Note that for $p cut=0$⁠, no edges are removed and time evolution is equivalent to the SIR dynamics on the original networks. For increasing $p cut$⁠, the number of adaptions increases and we find very good agreement between simulations on the ER network and the solution of the closed moment equation (see Fig. 8). The regime with $p cut≳0.8$ is characterized by initial transient oscillations similar to SIS epidemics. We observe that the strategy adaptions in scale-free networks occur much earlier than in the closed moment equations (not shown), which is expected due to the initial faster epidemic spreading in these networks.

In this paper, a self-adaptive mechanism with a finite set of strategies has first been formalized, which leads to a coupled system of piecewise deterministic dynamics with history dependent strategy switching.

This adaptive framework is applied to epidemic models on networks. For the SIS system, we observe a stable regime of strategy induced oscillations. Their period depends on the network type and the degree distribution. Based on the pairwise closed model, an analytic prediction for the period is derived, which fits very well for Erdös–Rényi networks and also approximates the period observed in scale-free networks. We believe that our model identifies one of the simplest and reduced mechanisms and how global observation of the incidence can lead to oscillations in epidemic dynamics.

We emphasize that the proposed adaptive mechanism deterministically depends on the state of the system and its history. One promising generalization is to consider a stochastic process for the adaptive mechanism, with state (and history) dependent transition rates between the strategies. This adaptive mechanism could be coupled with deterministic ODE models (like the pairwise closed SIS) or with stochastic models (like the SIS network model). The latter implies two distinct sources of random behavior, which could lead to interesting phenomena.

It would be further interesting to investigate in more detail the observation of the fluctuating reproductive rate around its critical value $R 0=1$⁠. There is a connection of this observation to self-organized criticality, which has been observed in a wide variety of adaptive/co-evolutionary network models.^58–60 We conjecture that having an observable controlling a switching mechanism that entails lowering or increasing the infection numbers is one simplest mechanism to obtain self-organized criticality in epidemic dynamics. One could now also investigate power law distributions of epidemic event sizes in various models that could further confirm this conjecture. Yet, for our SIS-based model, the underlying bifurcation-theoretic mechanism is the switching around a transcritical bifurcation, where one can calculate mathematically that power laws emerge close to the bifurcation point.^17,61 Hence, finding evidence for self-organized criticality in very complex epidemic models and/or long-time data sets with multiple outbreaks are the most challenging aspects.

K.C. and C.K. thank the VolkswagenStiftung for support via the grant “Self-Dynamics of Self-Adapting Networks” within a Lichtenberg Professorship awarded to C.K. We also thank two anonymous referees, who have helped to improve the exposition of the manuscript via their comments.

The authors have no conflicts to disclose.

Konstantin Clauß: Investigation (lead); Software (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). Christian Kuehn: Conceptualization (equal); Funding acquisition (lead); Investigation (supporting); Methodology (lead); Supervision (lead); Writing – review & editing (lead).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

In the following, we derive an analytic expression for the expected period in the oscillatory regime for the SIS model with threshold adaption of strategies. For a simplified analysis the pairwise closure, of the system is considered, Eq. (3), $x I ′=−γ x I+τk(1− x I) x I$⁠. The average node degree takes one of two values, $k∈{ k 0, k −}$⁠, referring to the strategies $S 0$ and $S −$⁠, respectively. If in $S −$ the probability to cut links is $p cut$⁠, then $k −= k 0(1− p cut)$⁠. The ratio of the expected rate of infection to the rate of recovery, $β(k)=τk/γ$⁠, also takes one of two values $β∈{ β 0, β −}$ with $β −=β( k −)= β 0(1− p cut)$⁠.

Note that a similar calculation for $J$ being the average incidence over past $Δt$ leads to the condition $ξ += τ k 0 Δ t ∫ t − Δ t tx( t ′)[1−x( t ′)] d t ′$ which becomes difficult to solve analytically.

In Fig. 9, we illustrate the period $T$ as a function of $ξ +$ for different values of $p cut$ within the periodic regime. For smaller $p cut$⁠, the system converges toward an endemic state and for the larger one to the disease-free state (see Fig. 2). Increasing (a) $Δt=0$ to (b) $Δt=5$ mainly influences the regime with small $ξ +$⁠. For $Δt=0$⁠, lower and upper boundaries of the periodic regime in terms of $ξ +$ are given by $max{ ξ −, x en( β −)}$ and $x en( β 0)$⁠, respectively, indicated by gray dashed lines.