The structure of coevolving infection networks
Disease awareness in infection dynamics can be modeled with adaptive contact networks whose rewiring rules reflect the attempt by susceptibles to avoid infectious contacts. Simulations of this type of models show an active phase with constant infected node density in which the interplay of disease dynamics and link rewiring prompts the convergence towards a well defined degree distribution, irrespective of the initial network topology. We develop a method to study this dynamic equilibrium and give an analytic description of the structure of the characteristic degree distributions and other network measures. The method applies to a broad class of systems and can be used to determine the steady-state topology of many other adaptive networks.
Stochastic analysis methods Stochastic modeling Structures and organization in complex systems
Research on the influence of network topology on the global dynamics of systems composed of many interacting units has been very active for the last ten years. More recently, adaptive and coevolving networks, that is, networks whose topology changes along with the dynamics, have begun to be considered on two standings. The first one is to explore the consequences of this interplay from the viewpoint of global dynamics . The second is to understand what sort of interactions may give rise to certain network architectures . Understanding the origin of degree distributions and correlations of real networks was first addressed by considering networks that evolve independently of the dynamics , and it is natural to extend this idea to networks whose evolution is determined by the dynamics they support.
Infection dynamics is a natural setting to investigate the coupled evolution of the interacting elements and the network of interactions, because disease awareness translates into susceptibles that try to evade infection by changing their contact patterns according to the disease status of their neighbors. The SIS (susceptible, infected, susceptible) and the SIR (susceptible, infected, recovered) models are the simplest models for the spread of epidemic infections . They are based on the idea that an infected individual infects susceptibles at a certain rate for a certain period, and then recovers and becomes susceptible again (if the disease does not confer immunity), or recovers and becomes immune. Stochastic versions of these models on networks can be mapped onto bond percolation .
Simple epidemic models like these exhibiting trivial global dynamics on static networks have been shown to display interesting behavior on adaptive networks, in which links between nodes may be created or removed according to rewiring rules dictated by disease awareness. In the network scheme of SIS dynamics, susceptible nodes (S-nodes, a fraction of the total number of nodes) are infected with rate along links connecting them to infected nodes (I-nodes, the remaining fraction ), whereas the latter recover with rate . Gross et al.  incorporated disease awareness into the conventional SIS model by introducing a topology-changing, yet link-number preserving, rewiring mechanism: S-nodes try to evade infection by retracting links from infected neighbors with rate , and rewiring them to randomly selected S-nodes. The description of the model in the framework of the standard pair approximation yields a complex phase diagram. In addition to a frozen phase (where the disease dies out) and to a stationary active phase (where is constant and , the number per node of links connecting nodes of type and , , is also constant), other active phases emerge that are absent in the SIS model without rewiring, namely an oscillatory phase and a bistable phase . Simulations in the stationary active phase show the states of the nodes and the links coevolving to produce and maintain a dynamic network topology characterized by well defined degree distributions not only for the global network, but also for the subsets of S- and I-nodes. This remarkable fact and the structure of the characteristic degree distribution were left mostly unexplained. In the wake of this first paper based on SIS infection dynamics, several contributions have addressed extensions and modifications  of the original model, focusing mostly on the global dynamics and the conditions for disease persistence. In a recent paper , steady-state degree distributions from individual-based simulations were shown to be very well approximated by the outcome of long-term numerical integration of a compartmental model for the evolution of the node status together with the number of its infected and susceptible neighbors. This model provides an alternative to simulations to produce approximately the observed network topology, but it is not amenable to analytic treatment. Overall, the description of the network topology underlying the steady states of these models has received comparatively little attention. Here we develop an analytic method to compute degree distributions and other measures that characterize this network topology.
2 The node-cycle model and its solution
We illustrate the method taking the original model  as an example of an adaptive network, and focusing on the stochastic process that each node and its links follow as the former undergoes status change from susceptible to infected, and gains (in the susceptible S-stage) or loses (in the infected I-stage) links. This process can be treated analytically, and the stationary degree distributions for the full node cycle (NC), from susceptible to infected and back again, can be obtained and compared with the degree distributions of the susceptible and infective subnetworks of the original network model. The NC approach is based on the ergodic properties of the network stationary state, which are a consequence of the random rewiring process that ensures mixing. Indeed, we observed in simulations of the steady state of the network model that recording ensemble statistics yields the same distributions as sampling a single node over a sufficiently long time.
In each stage of the NC we keep track of the numbers and of susceptible and infected neighbors of the node (its joint degree (x,y)). Due to infection and rewiring processes, these numbers change at certain rates, defining a time-homogeneous Markov process that can be described in each stage as a finite random walk on a degree grid spanned by and (Fig. 1). The time-continuous random walks performed in each stage are one-step processes in both coordinates. They are coupled by stage transition, which takes place at constant rate from the I-stage to the S-stage, and at rate from the S-stage to the I-stage. The probability of a walker in the S-stage, having started at coordinates , being at at time obeys the Master Equation
with the boundary conditions . The first term on the right hand side of Eq. (2) takes into account the swapping of an infected neighbor with a susceptible neighbor due either to rewiring or to recovery of the infected neighbor. The second term on the right represents the transition to the I-stage and an overall probability mass loss as in the active phase all susceptible nodes eventually become infected. The stochastic dynamics of the local random variables and is coupled to the global network dynamics through the total degree gain rate and through the force of infection that determines the infection rate of susceptible neighbors of the node. These two transitions are represented by the third and fourth terms on the right-hand side of Eq. (2). A similar Master Equation holds for the I-stage, where instead of the rate accounts for the force of infection on susceptible neighbors of an I-node. We assume the correspondence parameters , and to be constant and will assign values later on.
a linear first order PDE
Its solution, given that , is
where consist of linear combinations of exponential functions of time and .
The Taylor expansion of Eq. (2) up to arbitrary order is easy to compute, yielding . With the remaining probability mass at time being , the total occupation time at , having started at , reads as
Given that the rate of reinfection of a S-node with infected neighbours is , we calculate
to obtain the probability of, having started at , ending the walk in the S-stage at . In a similar way, a closed-form expression for the I-stage can be obtained, where because of recovery being coordinate-independent.
The elements of , , are the transition probabilities that map a random walker’s initial coordinate distribution (ICD) in stage , , onto the ICD of the subsequent stage according to
The ICD of the S-stage after a full node cycle changes according to
where the are given by the Chapman-Kolmogorov identity
Setting a cutoff for the maximum degree, the elements of encode a finite ergodic Markov chain, and therefore an arbitrary ICD in the S-stage will converge upon iteration of the node cycle to a unique stationary state which, according to the Perron-Frobenius theorem, can be found as the right eigenvector of associated with the eigenvalue . Once are known, the normalized total occupation times of coordinates in the stationary state can be computed in either stage as
where in particular
because of Eq. (4) and . In each stage , , the average number of infected and susceptible neighbors, and , and the average degree, , can be computed from the distributions .
The lifetime distributions in each stage of random walkers with initial coordinates and the generating functions are related through
with overall lifetime distributions
and survival functions in the steady state
The average duration of each stage in the stationary state is thus given by
It then follows from the solution of the analogue of Eq. (2) for the I-stage that as expected, because there random walkers switch to the S-stage at a coordinate-independent constant rate .
3 Correspondence with the network model and comparison with simulations
The stationary state of the NC process was solved analytically for any choice of the parameters , , and , , . These solutions encompass a much larger set of systems than the original network model of , because the correspondence parameters, tying network to NC dynamics, have yet to be determined. A first estimate for these parameters may be obtained from the equilibrium link number per node predicted by the pair approximation model of . Combining this approximation with the NC description yields , and similar expressions hold for the other correspondence parameters. The results for the subensemble degree distributions obtained through this approach reproduce the qualitative behaviour of the network simulations but the quantitative agreement is poor (results not shown). An alternative way of constraining the correspondence parameters comes out of embedding the NC description in a given network of fixed degree. On one hand, the correspondence parameters must be such that link depletion in the I-stage and link addition in the S-stage balance out for that degree. On the other hand, they should ensure that computing the numbers of links between nodes of type and does not depend on whether one sums over the infected neighbors of all the S-nodes or over the susceptible neighbors of all I-nodes. These consistency conditions will provide a self-contained determination of the correspondence parameters within the NC framework which is independent of the pair approximation and does not involve network simulations for parameter fitting.
For the NC process to represent the stationary dynamics of the network, so that the correspond to the network subensemble degree distributions and and are proportional to the fractions and of infected and susceptible nodes in the network, we shall therefore impose the following two conditions:
(i) The average node degree must be equal to the average degree that is chosen a priori in the network model,
where are the average degrees computed from Eq. (5).
(ii) The S- and I-stages represent nodes that are each other’s neighbors, and therefore for the mean number of susceptible and infected neighbors in the two stages, with the averages again computed from Eq. (5).
Computing the NC stationary state and choosing values of the correspondence parameters such that the consistency conditions (i) and (ii) hold, the NC description can be mapped onto the stationary state of a particular network model. However, the constraints set by (i) and (ii) may be impossible to satisfy exactly. Because of the mean field approximation involved in the assumption of constant and , the stochastic process we have dealt with does not capture the weak correlations present in the network. This translates into (i) and (ii) not being exactly fulfilled, and so the description of the stationary state given by the NC is only approximate. In order to compare the analytic results with simulations on networks, we have used the free parameters , and to minimize the squared distance of the NC output to the target values set by conditions (i) and (ii), bringing the stochastic model to represent as accurately as possible the state and degree evolution of a network node. We shall refer to the analytic stationary degree distributions obtained through the method outlined above for the optimal choice of , and as the NC degree distributions. Results for a high and a low rewiring regime in a network of average degree are presented in Fig. 2. For different choices of , and in the stationary active phase, the NC degree distributions as well as show very good quantitative agreement with those obtained from long Monte Carlo (MC) simulations of the network model of  and , as do cumulative node lifetime distributions as given by Eq. (2) (Fig. 3).
Our method predicts, through combinatorial handling of , steady-state densities of any star motif. With steady-state node densities being given through
the density of triplets with a central I-node connected to a S- and another I-node, normalized by the system size, is illustratively computed as
For the parameters of Fig. 2(a), this yields , with from MC simulations. The small discrepancies (see Fig. 2) between steady-state moment densities observed in MC simulations and those predicted by the NC method are another consequence of the constant and constant approximation. An improvement on this approximation would be to consider and to be linear functions of time chosen to represent status and degree correlations of the network model. For instance, should increase in time along the I-stage, because older (lower degree) I-nodes have susceptible neighbors that are younger than average S-nodes, and younger S-nodes in turn have more than the average number of infected neighbors. Although it would no longer be possible in this case to find a simple solution like Eq. (2) for the generating functions , the whole NC procedure could in principle be carried out numerically for any set of positive rate functions to yield an irreducible transition matrix , a unique stationary degree distribution for each stage, and improved approximations for all the moment densities that characterize the network’s steady state.
4 Discussion and conclusions
Apart from giving a complete description of the network’s steady state that is in excellent quantitative agreement with the results of MC simulations, the NC method also shows that convergence of an arbitrary ICD to a unique steady-state topology determined by the parameters , , and by the network’s fixed degree will take place once the rate functions associated with the correspondence parameters , and are defined. These may be simply chosen as the constants that optimize the consistency conditions (i) and (ii), or, as mentioned above, may be taken as functions of time with more free parameters, allowing a closer match of the NC model with the network-based system. In either case, encodes a time-homogeneous Markov chain, meaning that the method applies only to the network model’s stationary state and not to its transients. However, the correspondence parameters embody the mean-field approximations involved in the NC construction, and in the simplest case they are related only with two steady-state moment densities, and . If, as predicted in  and confirmed by simulations, in the stationary active phase these densities reach equilibrium independently of network structure, then the NC approach shows that this implies that the network topology, as described by various probability distributions, will reach equilibrium as well.
Hence the pairwise model of  and our NC approach should be considered complementary analytic frameworks to describe approximately adaptive networks. The former allows for a complete global treatment of system dynamics, at the cost of confining the level of description to low-order moment densities. The latter is a local model, in that it identifies active steady states and gives a comprehensive description of their equilibrium dynamics through various characteristic distributions.
In conclusion, taking the SIS network model with rewiring of Gross et al. as an example, we developed a self-sufficient method to derive analytically the steady-state subensemble degree distributions, and other measures that characterize the unique steady state in the simple endemic phase of the model. The only requirements of the construction are that we have a cyclic two-state system, that the state changes are due to node or contact processes with constant rates, and that the network dynamics is due to link processes also with constant rates. The requirement of constant rates may be relaxed at the expense of an analytic solution for the probability generating functions not being available in this case. As for the improved approximations briefly discussed in the previous section, dealing with non-linear transition rates would require the whole NC procedure to be carried out numerically to yield the stationary degree distribution for each stage. The requirement of a two-state system may also be relaxed to include cyclic multi-state systems such as SIRS or rock-paper-scissors dynamics . The method developed here can therefore serve as a tool to explore the relation between the node’s status change and rewiring rules and the resulting network structure in different contexts, ranging from infection to opinion dynamics.
Acknowledgements.The authors would like to thank Andrea Parisi for many helpful discussions. Financial support from the Portuguese Foundation for Science and Technology (FCT) under Contract POCTI/ISFL/2/261 is gratefully acknowledged. The first and second authors were also supported by FCT under Grant No. SFRH/BD/45179/2008 (S.W.) and CFTC-618-BII-02/08 (T.A.).
-  \NameP. Holme M. E. J. Newman \REVIEWPhys. Rev. E742006056108; \NameN. H. Fefferman K. L. Ng \REVIEWPhys. Rev. E762007031919; \NameS. Van Segbroeck, F. C. Santos, T. Lenaerts J. M. Pacheco \REVIEWPhys. Rev. Lett.1022009058105; \NameA. Szolnoki M. Perc \REVIEWEPL86200930007.
-  \NameM. G. Zimmermann, V. M. Eguíluz M. San Miguel \REVIEWPhys. Rev. E692004065102; \NameO. Gräser, C. Xu P. M. Hui \REVIEWEPL87200938003; \NameF. Vazquez, V. M. Eguíluz M. San Miguel \REVIEWPhys. Rev. Lett.1002008108702; \NameG. Iñiguez, J. Kertész, K. K. Kaski R. A. Barrio \REVIEWPhys. Rev. E802009066119.
-  \NameA.-L. Barabási R. Albert \REVIEWScience2861999509.
-  \NameJ. D. Murray\BookMathematical Biology\PublSpringer, New York \Year1993.
-  \NameP. Grassberger \REVIEWMath. Biosci.631983157; \NameM. E. J. Newman Phys. Rev. E662002016128; \NameL. M. Sander, C. P. Warren I. M. Sokolov \REVIEWPhysica A32520031.
-  \NameT. Gross, C. J. Dommar D’Lima B. Blasius \REVIEWPhys. Rev. Lett.962006208701.
-  \NameT. Gross I. G. Kevrekidis \REVIEWEPL82200838004.
-  \NameB. Bagnoli, P. Lió L. Sguanci \REVIEWPhys. Rev. E762007061904; \NameL. B. Shaw I. B. Schwartz \REVIEWPhys. Rev. E771008066101; \NameS. Risau-Gusman D. H. Zanette \REVIEWJ. Theor. Biol.257200952; \NameS. Funk, E. Gilad, C. Watkins V. A. A. Jansen \REVIEWProc. Natl. Acad. Sci. USA10620096872; \NameL. B. Shaw I. B. Schwartz \REVIEWPhys. Rev. E812010046120; \NameC. Lagorio, M. Dickison, F. Vazquez, L. A. Braunstein, P. A. Macri, M. V. Migueles, S. Havlin and H. E. Stanley \REVIEWPhys. Rev. E832011026102.
-  \NameV. Marceau, P. A. Noël, L. Hébert-Dufresne, A. Allard L. J. Dubé \REVIEWPhys. Rev. E822010036116.
-  \NameN. G. van Kampen\BookStochastic Processes in Physics and Chemistry\PublElsevier, Amsterdam \Year1981.
-  \NameD. T. Gillespie \REVIEWJ. Comput. Phys.221976403.
-  \NameG. Demirel, R. Prizak, P. N. Reddy T. Gross \ReviewEur. Phys. J. B\Year2011 \doihttp://dx.doi.org/10.1140/epjb/e2011-10844-4