key: cord-010705-lez9dcao
authors: de Oliveira, Marcelo M.; Alves, Sidiney G.; Ferreira, Silvio C.
title: Dynamical correlations and pairwise theory for the symbiotic contact process on networks
date: 2019-11-04
journal: nan
DOI: 10.1103/physreve.100.052302
sha: 
doc_id: 10705
cord_uid: lez9dcao

The two-species symbiotic contact process (2SCP) is a stochastic process in which each vertex of a graph may be vacant or host at most one individual of each species. Vertices with both species have a reduced death rate, representing a symbiotic interaction, while the dynamics evolves according to the standard (single species) contact process rules otherwise. We investigate the role of dynamical correlations on the 2SCP on homogeneous and heterogeneous networks using pairwise mean-field theory. This approach is compared with the ordinary one-site theory and stochastic simulations. We show that our approach significantly outperforms the one-site theory. In particular, the stationary state of the 2SCP model on random regular networks is very accurately reproduced by the pairwise mean-field, even for relatively small values of vertex degree, where expressive deviations of the standard mean-field are observed. The pairwise approach is also able to capture the transition points accurately for heterogeneous networks and provides rich phase diagrams with transitions not predicted by the one-site method. Our theoretical results are corroborated by extensive numerical simulations.

Discontinuous transitions [1] , observed in distinct classes of cooperative systems, have led to a renewed burst of interest within different contexts, such as social interactions [2, 3] , coinfections [4] [5] [6] , synchronization [1, 7] , and percolation [1, 7] , to cite only a few fundamental processes. Many of these investigations have focused on phenomena occurring at the top of complex networks [8] , which constitute basic substrates for describing interacting patterns of complex systems [9] [10] [11] . While in a continuous transition the order parameter varies continuously from zero, in the discontinuous case it suddenly jumps to a finite value at the transition point, involving a macroscopic portion of the system, and this transition is commonly referred to as catastrophic or abrupt [12] .

Coinfection epidemics [4] , i.e., when a host can be infected simultaneously by two distinct diseases, can result in coexisting thresholds when competitive interactions are considered [13, 14] . Nevertheless, cooperative or synergistic interactions result in richer phase diagrams, which may include discontinuous phase transitions, and they have been a topic of intense research [5, 6, [14] [15] [16] [17] [18] . Cooperation and competition have been investigated in multispecies interacting models in distinct ecological contexts [19] [20] [21] [22] [23] [24] . In particular, symbiotic interactions were recently investigated in a two-species contact process (2SCP) [24] . In the standard, * mmdeoliveira@ufsj.edu.br † sidiney@ufsj.edu.br ‡ silviojr@ufv.br single-species contact process (CP) [25] , individuals lay on the vertices of a graph, which was originally a lattice, but the process was later extended to networks [26, 27] . In the standard CP model, each individual can reproduce asexually at a rate λ, with its offspring occupying one empty nearest neighbor randomly chosen, or it can die at a rate μ (which can be taken as μ = 1 without loss of generality). In the 2SCP [24] , two species, A and B, inhabit the same graph. The symbiotic interaction is modeled via a reduced death rate, μ < 1, at sites doubly occupied (by one individual of each species). With the exception of this interaction, the two populations evolve independently according to the standard CP.

Apart from its interest as an elementary model of symbiosis, the 2SCP is fundamentally interesting for the study of nonequilibrium phase transitions. Extinction represents an absorbing state, i.e., a frozen state without fluctuations of populations [28] . On regular lattices, it was found that the 2SCP exhibits a continuous phase transition in one and two dimensions [24] . However, the transition becomes discontinuous in the regime of strong symbiosis if diffusion is introduced [29] . The 2SCP was recently investigated in complete graphs and random regular networks [30] , and it was conjectured that the nature of its transition changes at the upper critical dimension, from continuous to discontinuous. The phase diagram determining the regions of the 2SCP space parameter μ versus λ was determined in the ordinary, onesite mean-field level [24, 29, 30] , which neglects all dynamical correlations with the assumption of statistical independence among the states of individuals even if they are nearest neighbors. This assumption is a strong approximation that can make the theory inaccurate even for high-dimensional systems 1 such as complex networks [32] , where the smallworld property [8] resembles a mean-field (fully connected) regime. Dynamical correlations can be introduced using pairwise approximations [28, 33] where the mean-field equations for pairs of connected vertices are considered. Despite being inaccurate in low dimensions [28] , pairwise correlations greatly improve the theoretical results in the case of dynamical processes on complex networks [32, [34] [35] [36] [37] .

In the present work, we investigate the role of dynamical correlations in the 2SCP using a homogeneous pairwise meanfield (PMF) theory. The theoretical predictions are compared with stochastic simulations on different random networks, including homogeneous, Poissonian, and scale-free degree distributions [38] . We observe that the PMF theory substantially improves the ordinary mean-field in all investigated cases, being very accurate for determining the transition points and phase diagrams. Moreover, the PMF phase diagrams are more complex, showing discontinuous transitions with either symbiosis parameter μ or rate infection λ fixed in contrast with the ordinary mean-field where only the former one can happen.

The remainder of this paper is organized as follows. In Sec. II we present the definitions of the model and networks used as substrates for simulations. The mean-field theories are developed in Sec. III and Appendix while the simulation methods used in the present work are described in Sec. IV. Section V is devoted to a comparison between mean-field theories and simulation outcomes. Finally, our conclusions and prospects are drawn in Sec. VI.

We investigate the symmetric 2SCP [24] , in which the involved rates are the same for both species. The model is defined for a networked substrate as follows. Each vertex of the network can hold at most one individual of each species A and B. So, the state S of a vertex i can be empty, represented by S i = ×, occupied by only one individual of type A (S i = •) or B (S i = •), or occupied by one individual of each species, represented by S i = ••. The reproduction of new individuals of a given species happens independently of each other according to the standard CP rules [25, 28] . A vertex i with one individual of type A, S i = • or S i = ••, creates an offspring of type A at a randomly selected nearest neighbor j at a rate λ if j does not have an individual of type A (i.e., if S j = × or S j = •). Same rules and rates are used for the replication of a type B. Any individual in a vertex i spontaneously dies with rate 1 if they are alone (S i = • or •) or with rate μ otherwise (S i = ••). If μ = 1, the two processes evolve independently, but for μ < 1 they interact symbiotically since the annihilation rates are reduced at vertices where both species are present. The stationary 2SCP can evolve asymptotically to a fully active state, where both species coexist, or to an inactive (absorbing) phase in which both species are extinct. There are also two partly active states, where only one of species A or B is extinct, that can be reached for specific initial conditions (the dynamics in these partly active states is beyond our interest since it is that of a standard single-species CP).

In this work, we use random graph models as substrates on which the 2SCP takes place. The number of neighbors of a vertex i, the vertex degree, is denoted by k i . We consider both homogeneous and heterogeneous networks. In random regular (RR) networks [39] , all vertices have the same degree, k i = k, and connections are performed at random following the configuration model [40] avoiding both multiple and selfconnections. In the Erdös-Renyi (ER) model [8] , each pair of vertices is connected with probability p. When the size of the graph N → ∞, its degree distribution is a Poissonian with a finite mean k = pN. Both RR and ER models belong to the class of homogeneous degree networks, where large deviations from the average value do not occur [31] , whereas only RR has a strictly homogeneous degree distribution. The third substrate is the Barabási-Albert (BA) model [38] , in which the networks are generated through a preferential attachment mechanism [38] starting from a fully connected set with m 0 + 1 vertices and vertices with m edges being added one at a time. We used m = m 0 such that the average degree k = m can be fixed accordingly for comparison with the homogeneous case. The BA degree distribution follows asymptotically a power law (PL), p k ∼ k −γ , with γ = 3 and is representative of a highly heterogeneous network possessing a heavy-tailed degree distribution.

We consider a homogeneous mean-field approximation where all vertices are assumed to have the same degree k i = k. To apply this theory to the heterogeneous cases of ER and BA, we use the same equations of the homogeneous theory replacing k by k . This strategy has been used to compare the transition points of the CP obtained in a homogeneous PMF theory, given by λ c = k/(k − 1) [28] , with the simulations on heterogeneous networks using λ c = k /( k − 1) [27, 41] .

Let [S] be the probability that a vertex is in state S, and let [S, S ] be the joint probability that a vertex is in the state S and its nearest neighbor is in the state S . Symmetries in the rates with respect to distinct species imply that 

Following this approach, one identifies three independent one-site variables, [•], [••], and [×], related by the closure relation

The dynamical equations for these variables are readily derived as

and

Equations (2)-(4) are exact but not closed since the evolution of the one-site probabilities depends on pairs. Cutting the correlations at a vertex level with the approximation [S, S ] = [S][S ], we obtain the ordinary mean-field equations [24] 

and

Note that if one species is extinct, the above system reduces to the mean-field theory for the one-species CP,

with a transition point at λ = 1.

The stationary solution of Eqs. (5)- (7) is given by [24] [

and

For μ 1/2, [•] grows continuously from zero at λ = 1, marking the latter value as the transition point. The activity grows linearly,

For μ < 1/2, however, the expression is already positive for λ = √ 4μ(1 − μ) < 1, and there is a discontinuous transition at this point.

Further improvement in including dynamical correlations is given by a PMF theory, in which the system is described in terms of pairwise variables [S, S ]. The dynamical equations of pairs will depend on triplets [S, S , S ] in which a pair S, S of nearest neighbors is connected to a vertex in the state S through the vertex in the state S . Here, we neglect that S and S can also be connected, i.e., we assume the network does not form triangles having thus a negligible clustering coefficient [42] We now proceed with the standard pairwise approximation [33, 37] [S,

and apply the generic closure relation

to obtain a set of seven dynamical equations for the independent pairwise variables, given by Eqs. (A8)-(A15) in Appendix. These equations, in addition to the one-site dynamical equations (2)-(4), build a closed system. It is worth mentioning that pairwise approximations have been presented for one-dimensional chains and square lattices in Ref. [29] while our approach is valid for generic homogeneous graphs (conditioned to have a very low clustering coefficient).

Our simulations of the 2SCP on networks were implemented using an optimized Gillespie algorithm [43] , in which we maintain two lists, one of singly and another of doubly occupied vertices. Let N 1 and N 2 denote, respectively, the numbers of such nodes. The total rate of (attempted) transitions is λN 1 + 2λN 2 + N 1 + 2μN 2 ≡ ( t ) −1 , where t is the average time increment associated with a given step of the simulation.

At each time step, we randomly choose among the events: (i) creation attempt by an isolated individual, with probability λN 1 t; (ii) creation attempt by an individual at a doubly occupied node, with probability 2λN 2 t; (iii) death of an isolated individual, with probability N 1 t; (iv) death of an individual at a doubly occupied node, with probability 2μN 2 .

Once the event type is selected, a vertex i is chosen at random from the respective list. The creation attempt occurs at a vertex j randomly selected among the nearest neighbors of i. If j is already occupied by an individual of the species to be created, no change of state is implemented, and the simulation continues to the next step. If node i is doubly occupied, the species of its offspring in a creation event is chosen to be A or B with equal probability. In the same way, in an annihilation event at a doubly occupied node, we choose the species to be removed at random. Time is incremented by t.

In the simulations, we sample the quasistationary (QS) distribution of active nodes employing the QS simulation method [44] . This scheme consists in replacing the absorbing state, every time the system attempts to visit it, with an active configuration randomly taken from the history of the simulation. This procedure optimizes the numerical simulations restricting the dynamics of the process to active states, and it is a powerful tool in analyzing continuous absorbing state phase transitions in networks [27, 37, [45] [46] [47] . For a recent analysis on such simulation methods applied to networks, see Ref. [48] . In this work, we performed QS simulations for systems of sizes up to N = 10 6 nodes, with each run lasting at least 10 6 time units. Averages are taken in the QS regime, after discarding an initial transient which depends on the system size and symbiosis strength used.

Finally, mean-field analysis was performed via a numerical integration of the respective dynamical systems using the fourth-order Runge-Kutta method with a time step t = 10 −5 . The steady state is computed after a relaxation of t > 10 6 .

The quasistationary densities of sites occupied by a single ρ 1 = ρ • + ρ • and by two species ρ 2 = ρ •• are compared with mean-field approximations in Fig. 1 for a RR network with connectivity degree k = 6, using distinct values of the symbiotic strength parameter μ. This dynamics can exhibit bistability with both active and absorbing stable stationary states [24] . So, we first analyze an initial condition having all vertices in a doubly occupied state ρ 2 = 1 such that this threshold represents the loss of global stability of the absorbing state, marking the lower spinodal point. We observe that the ordinary mean-field theory, although qualitatively predicting the discontinuous or continuous nature of the transition, is not quantitatively accurate as indicated by the dashed curves in Fig. 1 . However, a far better result is obtained with the PMF theory, represented by the solid curves, presenting an excellent agreement between theory and simulations, for both near and above the transition point regimes irrespective of the μ values. Increasing the connectivity k reduces the transition value of λ as can be seen in Table I . Actually, the threshold will converge to the one-site mean-field value for k 1 since this limit corresponds to the fully connected graph for which the one-site theory is exact in the thermodynamical limit. Figure 2 compares the results from simulations and meanfield theories for ER networks with k = 6 using the substitution of k by the mean connectivity k . As observed in RR networks, the one-site approximation reproduces qualitatively well the nature of transition but underestimates the transition point. The PMF accurately matches the transition point and density of singly occupied vertices in the active phase of simulations but underestimates the density of doubly occupied vertices and consequently the overall density of active vertices ρ = ρ 1 + ρ 2 .

The role of heterogeneity is further investigated in Fig. 3 where the results on BA networks with k = 10 are shown. The PMF theory outperforms the ordinary one but still underestimates the location of the phase transition for small values of μ, where the phase transition is discontinuous while yields a very good agreement for the transition point for the continuous cases with μ > 0.5 despite the high heterogeneity of the BA networks, where hubs are expected to play some important role. This finding in the continuous case is in agreement with one-species CP on scale-free networks [27, 41] .

The transition points for the investigated networks (RR, ER, and BA) with different average degrees are presented in Table I . We observe that the PMF theory exhibits a higher accuracy for large average connectivity, as expected, since in the limit k 1 it corresponds to the fully connected graph where mean-field theories become exact. Another important observation is that the PMF performs worst for BA networks when μ is small while the accuracy tends to be reduced for larger μ in RR networks. In BA networks, where hubs are present, small μ will enhance activity localized in the neighborhood of hubs, which, in turn, are not sufficiently mixed for the regime of scale-free networks with γ = 3 [49] . More precisely, the star subgraph containing a hub and its nearest neighbors can stay active in isolation for long periods through a feedback mechanism where the hub activates its neighbors, which in turn reactivate the hub recurrently [45, 50, 51] . Since localization is opposed to the homogeneous mixing hypothesis of the mean-field methods, larger deviations from the theory are indeed expected in this regime. Another intriguing feature of the data of Table I is that the PMF thresholds are nearer to the simulations for the slightly heterogeneous ER than from the homogeneous RR networks, which becomes more evident for lower k and higher μ. A similar phenomenon was also observed for the one-species CP on networks [37] and associated with the approximation k ≈ k , which would not be observed in a degree-based mean-field theory [9] .

Up to this point, we have addressed the global loss of the stability of the absorbing state that, in the case of discontinuous transitions, involves the transition from an absorbing to a bistable stationary state where both the absorbing phase and active phases can be stable depending on the initial condition. We now also consider the loss of global stability of the active phase, i.e., the transition from bistable to active phases. For this aim, we consider an initial condition very close to the absorbing state with a very low density of doubly occupied vertices (ρ 1 = 0 and ρ •• 1) in conjunction with the fully occupied state ρ •• = 1 used previously. We can therefore compute hysteresis curves where the phase coexistence can active FIG. 4 . Hysteresis analysis for 2SCP on RR networks with k = 6 and μ = 0.2. Solid curves: PMF theory (arrows indicate the hysteresis flow). Dashed curves: Ordinary mean-field theory. Symbols: Results from simulations for the RR network with N = 10 5 vertices using different initial conditions. be derived as in Fig. 4 where one sees absorbing, bistable, and active regions.

Phase diagrams in the parameter space (λ, μ) obtained via PMF for different connectivities are shown in Fig. 5 . The diagrams exhibit three connected regions corresponding to globally inactive (dotted region), globally active (empty region), and bistable (hashed region). An interesting difference with respect to the ordinary mean-field is a transition from a globally stable active phase to a bistable dynamics obtained either by fixing the infection rate λ and varying the symbiotic parameter μ or using fixed μ and changing λ, whereas in an ordinary mean-field theory such a transition occurs only for fixed μ; see, for example, Fig. 3 of Ref. [30] , which resembles very much the limit of large degree shown in Fig. 5(c) . Indeed, the region with double transitions (inactive to bistable and bistable to active) at fixed λ, occurring in the small μ region, shrinks as the connectivity increases toward a fully connected graph limit as shown in Figs. 5(b) and 5(c). Triple points shown in the phase diagrams depend on connectivity, departing from (λ t , μ t ) = (1, 1/2) for the complete graph (one-site mean-field) [24] to λ t = k k −1 and μ t = 0.501(1), 0.545(5), and 0.595(5) for k = 10 3 , 20, and 6, respectively. It is worthwhile to comment that phase diagrams with qualitatively similar spinoidals found in 2SCP with our PMF theory have been reported for other models with abrupt transitions on networks [52] [53] [54] , with the difference that they were obtained by a first-order mean-field theory, whereas in the 2SCP model we needed to go up to a second-order pairwise theory. Figure 5 (a) also presents the phase diagram obtained via simulations on different graphs with average degree k = 6. One observes a remarkably good match between simulations and theory for RR and ER (homogeneous) networks. The quantitative agreement is less satisfactory for the BA model, but it is still qualitatively correct. The dynamical correlations in these infinite dimensional systems have proven to lead not just to a quantitatively better theory but also to new behaviors not predicted by the ordinary mean-field theory. 

Stochastic interacting systems with cooperative or synergistic couplings can be used for modeling different phenomena, such as coinfections of pathological agents, social relations, and ecological interactions, among others.

One key issue of these models involves the nonequilibrium phase transitions among active and absorbing phases that can emerge in some systems, especially the abrupt ones, where macroscopic order parameters change discontinuously. Moreover, understanding the behavior of these transitions taking place at the top of complex networks has become very relevant since most of these dynamical processes indeed involve interactions mediated by networked systems. In the present work, we investigate a simple model in which two species lying on networks interact symbiotically, i.e., the symbiotic contact process [24] . We used both stochastic simulations and mean-field pairwise approximations, the latter reckoning dynamical correlations not considered in previous analytical (one-site mean-field) studies of this model.

The PMF theory outperforms the ordinary one and, more importantly, predicts features not captured by the latter, in agreement with the simulations on distinct types of networks with homogeneous, Poissonian, and scale-free degree distributions. In the case of homogeneous networks, the PMF theory is accurate at describing very well both the dynamics near and above the transitions (supercritical region). Although it disregards the heterogeneity of the networks, our theory also gives very good estimates of the transition points in the case of Poissonian and scale-free networks, but it cannot quantitatively capture the prevalence in the highly active regime.

We also investigated the phase diagrams exhibiting globally absorbing, bistable, and globally active regions. The PMF theory, beyond being quantitatively more accurate than the ordinary one, provides a much richer phase diagram where discontinuous transitions, involving bistable phases can be obtained by either varying the infection rate λ or the symbiosis parameter μ independently, whereas this transition happens only at fixed μ in the one-site approach. These PMF results are confirmed by extensive numerical simulations. Our findings provide an important example of the role played by dynamical correlations in cooperative processes on networked substrates, which is commonly not considered in other theoretical approaches. Pair approximations such as the one developed in this work can be used to improve the accuracy of phase diagrams observed in other related works [52, 55] where first-order mean-field theory captures the qualitative behavior but analytical spinoidals deviate from simulations when the average connectivity is reduced. We expect that our work will stimulate future analysis of dynamical correlations on other cooperative processes.

The theories presented in this paper do not include the network heterogeneity and, especially, the degree distribution. As a prospect for forthcoming analysis, one could use degree-based [9] or individual-based [56] mean-field theories to tackle the role of heterogeneity.

Explosive phenomena in complex networks

Statistical physics of human cooperation

Coevolution spreading in complex networks

Avalanche outbreaks emerging in cooperative contagions

Phase transitions in cooperative coinfections: Simulation results for networks and lattices

Explosive transitions in complex networks' structure and dynamics: Percolation and synchronization

Statistical mechanics of complex networks

Epidemic processes in complex networks

Statistical physics of social dynamics

A simple rule for the evolution of cooperation on graphs and social networks

Threshold Effects for Two Pathogens Spreading on a Network

Interacting epidemics and coinfection on contact networks

Outbreaks of coinfections: The critical role of cooperativity

First-order phase transitions in outbreaks of co-infectious diseases and the extended general epidemic process

Mutually cooperative epidemics on power-law networks

Explosive spreading on complex networks: The role of synergy

The paradox of enrichment in phytoplankton by induced competitive interactions

Environmental Versus Demographic Variability in Two-Species Predator-Prey Models

Parasites on parasites: Coupled fluctuations in stacked contact processes

Incorporating spatial correlations into multispecies mean-field models

A simple population theory for mutualism by the use of lattice gas model

Symbiotic two-species contact process

Contact interactions on a lattice

Non-Mean-Field Behavior of the Contact Process on Scale-Free Networks

Quasistationary simulations of the contact process on quenched networks

Nonequilibrium Phase Transitions in Lattice Models, Collection Alea-Saclay: Monographs and Texts in Statistical Physics

Phase diagram of the symbiotic two-species contact process

Symbiotic contact process: Phase transitions, hysteresis cycles, and bistability

Network Science

High-Accuracy Approximation of Binary-State Dynamics on Networks

Mean-field (n, m)-cluster approximation for lattice models

The effects of local spatial structure on epidemiological invasions

Modeling dynamic and network heterogeneities in the spread of sexually transmitted diseases

Binary-State Dynamics on Complex Networks: Pair Approximation and Beyond

Heterogeneous pair-approximation for the contact process on complex networks

Emergence of scaling in random networks

Critical behavior of the contact process on small-world networks

Generation of uncorrelated random scale-free networks

Griffiths Phases on Complex Networks

Collective dynamics of 'smallworld' networks

Optimized Gillespie algorithms for the simulation of Markovian epidemic processes on large and heterogeneous networks

How to simulate the quasistationary state

Collective versus hub activation of epidemic phases on networks

Disease localization in multilayer networks

Analytical Computation of the Epidemic Threshold on Temporal Networks

Sampling methods for the quasistationary regime of epidemic processes on regular and complex networks

Universal scaling of distances in complex networks

Contact processes on random graphs with power law degree distributions have critical value 0

Nature of the Epidemic Threshold for the Susceptible-InFected-Susceptible Dynamics in Networks

Spontaneous recovery in dynamical networks

Critical behaviors in contagion dynamics

Failurerecovery model with competition between failures in complex networks: A dynamical approach

Failure and recovery in dynamical networks

Epidemic phase transition of the SIS type in networks

The authors acknowledge the financial support of Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG).

In the pair approximation, the seven two-site variables to be considered are [ . All other combinations are equivalent to some of them, due to symmetry. Considering that every vertex has k nearest-neighbors, after some algebra we obtain the following set of dynamical equations in Eqs. (A1) to (A7), and find that they can be approximated as