key: cord-125979-2c2agvex
authors: Mata, Ang'elica S.
title: An overview of epidemic models with phase transitions to absorbing states running on top of complex networks
date: 2020-10-05
journal: nan
DOI: nan
sha: 
doc_id: 125979
cord_uid: 2c2agvex

Dynamical systems running on the top of complex networks has been extensively investigated for decades. But this topic still remains among the most relevant issues in complex network theory due to its range of applicability. The contact process (CP) and the susceptible-infected-susceptible (SIS) model are used quite often to describe epidemic dynamics. Despite their simplicity, these models are robust to predict the kernel of real situations. In this work, we review concisely both processes that are well-known and very applied examples of models that exhibit absorbing-state phase transitions. In the epidemic scenario, individuals can be infected or susceptible. A phase transition between a disease-free (absorbing) state and an active stationary phase (where a fraction of the population is infected) are separated by an epidemic threshold. For the SIS model, the central issue is to determine this epidemic threshold on heterogeneous networks. For the CP model, the main interest is to relate critical exponents with statistical properties of the network.

The present paper briefly reviews the modeling and theory of non-equilibrium dynamical systems on networks. A key class of non-equilibrium process are those that exhibit absorbing states, i.e. states from which the dynamics cannot escape once it has fallen onto them. A relevant feature of systems that presents absorbing states is a non-equilibrium phase transitions among an active state, in which the activity lasts forever in the thermodynamic limit, and an absorbing state, in which activity is absent [1, 2] .

The same type of transition occurs in epidemic spreading processes [3] since a fully healthy state is absorbing in the above sense, provided that spontaneous birth of infected individuals is not allowed. The susceptibleinfected-susceptible (SIS) [4] model and the contact process (CP) [5] represent some of the simplest epidemic models possessing an absorbing-state phase transition.

Lattice systems that exhibit such absorbing state phase transitions have universal features, determined by conservation laws and symmetries which allow to group them in a same universality class 1 [6] . The most robust class of absorbing state phase transitions is the directed percolation that was originally introduced as a model for directed random connectivity [7] . Both CP and SIS models are interacting particle systems involving self-annihilation and catalytic creation of particles that presents an absorbing-phase transition and thus belong to directed percolation class. The SIS dynamics is indeed the most studied model to describe epidemic spreading on networks. Although the CP was initially thought as a toy model for epidemics, lately it has been widely used as a generic reaction-diffusion model to study phase transition with absorbing states.

Other epidemic model that also presents an absorbing phase transition is the susceptible-infected-recoveredsusceptible model (SIRS) [4] . It is an extension of the standard SIS model, allowing a temporary immunity of nodes. Both SIS and SIRS models are equivalent from the mean-field theory perspective, but the mechanism of immunization changes the behavior of the epidemic dynamics depending on the heterogeneity of the network structure. The susceptible-infected-recovered (SIR) model is another example of epidemic models with permanent immunity, it means that a recovered node can no longer return to the susceptible compartment, so the system present many absorbing states since each configuration that have only susceptible and recovered nodes is absorbing [8, 9] .

In the face of this context, we reviewed the SIS and CP models as examples to investigate absorbing phase transitions in complex networks. We firstly explain, in section II, some basic concepts related to complex networks required to understand the main idea of this paper. Then we describe both epidemic models in section III and, in the section IV, we present distinct theoretical approaches devised for them. In section V, we described some commonly used simulation techniques to analyze both models numerically. For the SIS model, the central issue is to determine an epidemic threshold separating an absorbing, disease-free state from an active phase on heterogeneous networks [10] [11] [12] [13] [14] [15] [16] [17] [18] . While for the CP model, most of the interest is to relate the critical exponents with statistical properties of the network, in particular the degree distribution [19] [20] [21] [22] [23] [24] . In sections VI and VII we present a discuss about these points related to CP and SIS models, respectively. Finally, in section VIII we draw our final comments.

Network analysis is a powerful tool that provide us a fruitful framework to describe phenomena related to real-world complex systems. Here we will describe just some features of complex networks that it will be used throughout the paper. We will also present the uncorrelated configuration model (UCM) [25] , the substrate that will be used to model the dynamics of the epidemic process on networks.

We can represent a network by means of an adjacency matrix A. A graph of N vertices has a N × N adjacency matrix. The edges can be represented by the elements A ij of this matrix such that [26] A ij = 1, if the vertices i and j are connected 0, otherwise,

for a undirected and unweighted graph. In this case, the adjacency matrix is symmetric, it means A ij = A ji . A relevant information gives from the adjacency matrix is the degree k i of a vertex i defined as the number of links that the vertex i has, i.e., the number of nearest neighbors of the vertex i. The degree of the vertices can be written as [27] 

When it concerns to very large systems a suitable description can be done by means of statistical measures as the degree distribution P (k). The degree distribution provides the probability that a vertex chosen at random has k edges [28, 29] . The average degree is an information that can be extracted from P (k) and it is given by the average value of k over the network, it means,

Similarly, it can be useful to generalize and calculate the n-th moment of the degree distribution [27] k n = k k n P (k).

(4)

We can classify networks according to their degree distribution. The basic classes are homogeneous and heterogeneous networks. The first ones exhibit a fast decaying tail, as for example, a Poisson distribution. Here the average degree value corresponds to the typical value in the system. Heterogeneous networks exhibit heavy tail that can be approximated by a power-law decay, P (k) ∼ k −γ . In this kind of network, the vertices often have a small degree, but there is a non-negligible probability of finding nodes with very large degree values thus, depending on γ, the average degree does not represent any characteristic value of the distribution [27] .

Many network model have been created in order to describe real systems. The advantage of using a model is to reduce the complexity of the real world to a level that one can be treated, for example, from the perspective of mean-field approach. In this context, uncorrelated random graphs are important from a numerical point of view, since we can test the behavior of dynamical systems whose theoretical solution is available only in the absence of correlations. For this propose, Catanzaro and collaborators [25] presented an algorithm to generate uncorrelated random networks with power law degree distributions, called uncorrelated configuration model (UCM) as described below.

To construct this networks we started with a set of N disconnected vertices. Each node i is signed with a number k i of stubs, where k i is a random variable with distribution P (k) ∼ k −γ under the restrictions k 0 ≤ k i ≤ N 1/2 and i k i even. It means that no vertex can have either a degree smaller than the minimum degree k 0 or larger than the cutoff k c = N 1/2 . The network is constructed by randomly choosing two stubs and connecting them to form links, avoiding both self and multiple connections [25] .

It is possible to show that [30] , to avoid correlations in the absence of multiples and self-connections, the random network must have a structural cutoff scaling at most as k c (N ) ∼ N 1/2 . As said previously, this algorithm is very useful in order to check the accuracy of many analytical solutions of dynamical process on networks. Because of that, it was chosen as a substrate for implementing the dynamics of the SIS and CP models in this review work.

In the SIS epidemic model, each vertex i of the network can be only one of two states: infected or susceptible. Let us assume the most general case where a vertex i becomes spontaneously healthy at rate µ i , and transmits the infection to each one of its k i neighbors at rate λ i . For classical SIS, one has µ i = µ and λ i = λ for every vertex [9] .

As in the SIS model, vertices in the CP model can be infected or susceptible, which in reaction-diffusion system's jargon are called occupied and empty, respectively. The spontaneous cure process is exactly the same as in the SIS model: infected vertices become susceptible at rate µ i = µ. However, the infection is different. An infected vertex tries to transmit the infection to a randomly chosen neighbor at rate λ, implying that the transmission rate of vertex i is λ i = λ/k i , where k i is the number of neighbors of the i-th node. This reduces drastically the infective power of very connected vertices in comparison with the SIS dynamics.

Both SIS and CP dynamics exhibit a phase transition between a disease-free (absorbing) state and an active stationary phase where a fraction of the population is infected. Originally, these regions are separated by an epidemic threshold λ c [8, 9] . The density of infected nodes ρ is the standard order parameter that describes this phase transition, as shown in Figure 1 . However, for a finite system the unique true stationary state is the absorbing state, even above the critical point, due to dynamical fluctuations. To overcome the difficulty to study the active state of finite systems some simulation strategies were proposed in the literature [2, 18, 31, 32] , as we will show in section V C. The usual behavior of the density of infected nodes ρ in function of the control parameter λ, in a epidemic model as SIS or CP in the thermodynamic limit. The value λc is the epidemic threshold that separates an absorbing state to an active phase with ρ > 0.

When we study epidemic processes running on the top of heterogeneous networks a more complex behavior can emerge. Indeed, the accurate theoretical understanding of epidemic models running on the top of complex networks rates among the hottest issues in the physics community [10, 11, 14, [16] [17] [18] [19] [20] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] . Much effort has been devoted to understand the criticality of the absorbing state phase transitions observed in CP [19] [20] [21] [22] 36] and SIS [10, 11, 14, 16-18, 35, 37-42] models, mainly based on perturbative approaches (first order in ρ) around the transition point [10, 11, 16, 17, 19, 34] . In the following we will present the basic mathematical approaches for the epidemic dynamics.

The prediction of disease evolution can be conceptualized within a variety of mathematical approaches. All theories aim at understanding the properties of epidemics in the equilibrium or long term steady state, the existence of a non-zero density of infected individuals, the presence or absence of a threshold, etc. We will start with the simplest case namely the homogeneous meanfield theory, and thereafter we will review other more sophisticated mathematical approaches.

The simplest theory of epidemic spreading assumes that the population can be divided into different compartments according to the stage of the disease (for example, susceptible and infected in both SIS and CP models) and within each compartment, individuals (vertices in the complex networks' jargon) are assumed to be identical and have approximately the same number of neighbors (edges), k ≈ k . The idea is to write a time evolution equation for the number of infected individuals I(t), or equivalently the corresponding density ρ(t) = I(t)/N , where N is the total number of individuals. For example, the equation describing the evolution of the SIS model is [43] :

considering uniform infection and cure rates (λ i = λ and µ i = µ, ∀ i = 1, 2, · · · N ). The first term on the righthand side of the Eq. (5) refers to the spontaneous healing and the second one to the infection process, that is proportional to the spreading rate λ k , the density of the susceptible vertices is 1 − ρ(t) that may become infected, and the density of infected nodes ρ(t) in contact with any susceptible individual. Note that the evolution of the SIS model is completely described by Eq. (5), since the density of susceptible individuals is S(t)/N = 1−ρ(t).

The mean-field character of this equation comes from the fact that the correlations among different nodes were neglected. Thus, the probability that one infected vertex is connected to a susceptible one is approached as

Near the phase transition between an absorbing state and an active stationary phase we can assume that the number of infected nodes is small ρ(t)

1. In this regime, we can use a linear approximation neglecting all ρ 2 terms 2 . So the Eq. (5) becomes 3 :

The solution is ρ(t) ∼ e −(1−λ k )t , implying that ρ = 0 is an stable fixed point for 1 − λ k > 0. Thus, one obtains,

Here, λ c is the epidemic threshold such that for any infection rate above this value the epidemic lasts forever [43] . Similar analysis can be done for the CP model. In this case, the homogeneous mean-field equation read as

since the transmission rate of each node is λ/ k . Performing the same linear stability analysis in the steadystate, one obtains λ c = 1.

In this framework, one considers that the connectivity patterns among individuals are homogeneous, neglecting the highly heterogeneous structure of the contact network inherent to real systems [28] . Many biological, social and technological systems are characterized by heavy tailed distributions of the number of contacts k of an individual (the vertex degree), characterized by a power law degree distribution, P (k) ∼ k −γ . For such systems the homogeneity hypothesis is severely violated [27, 28, 44] . Complex networks are, in fact, a framework where the heterogeneity of the contacts can be naturally afforded [28] . Indeed, this heterogeneity plays the main role in determining the epidemic threshold. To take it into account other approaches have been proposed, as we showed in the next sections. The major aim is to understand how the epidemic spreading can be strongly influenced by the topology of networks.

An important and frequently used mean-field approach to describe epidemic dynamics on heterogeneous networks is the quenched mean-field (QMF) theory [10] , that explicitly takes into account the actual connectivity of the network through its adjacency matrix. The central idea is to write the evolution equation for the probability ρ i (t) that a certain node i is infected. For the SIS model the dynamical equation for this probability takes the form [10] :

where A ij is the adjacency matrix. The first term on the right-hand side considers nodes becoming healthy spontaneously while the second one considers the event that the node i is healthy and gets the infection via a neighbor node. Performing a linear stability analysis around the trivial fixed point ρ i = 0, one has

where the Jacobian matrix is

δ ij being the Kronecker delta symbol. The transition occurs when the fixed point loses stability or, equivalently, when the largest eigenvalue of the Jacobian matrix is Υ m = 0 [45] . The largest eigenvalue of L ij is given by

Since A ij is a real non-negative symmetric matrix, the Perron-Frobenius theorem states that one of its eigenvalues is positive and greater than, in absolute value, all other eigenvalues, and its corresponding eigenvector has positive components. So, one obtains the epidemic threshold of the SIS model in a QMF approach [10] :

where Λ m is the largest eigenvalue of the adjacency matrix. For the CP dynamic, the Eq. (9) becomes,

Performing the same linear stability analysis around the trivial fixed point, as was done for the SIS model, one obtains

where the Jacobian matrix is given by

Once again the transition point is defined when the absorbing phase becomes unstable or, equivalently, when the largest eigenvalue of L ij is null [45] . The largest eigenvalue of L ij is given by

Now, supported by the Perron-Frobenius theorem [44] , we conclude that the largest eigenvalue of C ij is Λ m = 1 resulting in the transition point λ c = 1, as obtained in a homogeneous approximation.

Returning to the SIS model, the equation (12) can be complemented with the results of Chung et. al. [46] who calculated the largest eigenvalue of adjacency matrix of networks with a power law degree distributions as

where k c is the degree of the most connected node and k n is the n-th moment of the degree distribution. Since k c grows as a function of the network size for any γ and k 2 diverges for 2 < γ < 3, the central result of equation (16) is: Λ m diverges for enlarging networks with power law degree distributions even when k 2 remains finite [46] . Therefore, the epidemic thresholds scale as [47] 

which vanishes for any power-law degree distribution. The reasons for this difference of λ c predicts by QMF approach, for γ larger or smaller than 5/2 are explained in ref. [47] . In processes allowing endemic steady-states, the activation mechanisms depend on the degree of heterogeneity of the network. For γ > 5/2 the hub sustains activity and propagates it to the rest of the system while for γ < 5/2 the innermost network core collectively turns into the active state maintaining it globally [47] [48] [49] . However, the behavior of the SIS model on random networks with power-law degree distribution can be much more complex than previously thought. We will discuss these different possible scenarios in section VII. Although mean-field theories are a simplified description of models, it is expected that they correctly predicts the behavior of dynamical process on networks, due to its small-world property. However, dynamical correlations are not taken into account since the states of a node and its neighbors are considered independent. One can consider dynamical correlations by means of a pairapproximation [2] in which the dynamic of an individual is explicitly influenced by its nearest neighbors as we showed in the section IV D 1.

In the degree-based theories, called heterogeneous mean-field (HMF) theory, dynamical quantities, as the density of infected individuals, depend only of the vertex degree and do not of their specific location in the network. Actually, the HMF theory can be obtained from the QMF one performing a coarse-graining where vertices are grouped according to their degrees. To take into account the effect of the degree heterogeneity we have to consider the relative density ρ k (t) of infected nodes with a given degree k, i.e., the probability that a node with k links is infected. Again using the SIS model as an example, the dynamical mean-field rate equation describing the system can thus be written as [9] :

The first term on the right-hand side considers nodes becoming healthy at unitary rate. The second term considers the event that a node with k links is healthy and gets the infection via a nearest neighbor. The probability of this event is proportional to the infection rate λ, the number of connections k and the probability that any neighbor vertex is infected P (k |k)ρ k . The linearization of Eq. (18) gives

where L kk = −δ kk +λkP (k |k). Therefore, the epidemic threshold is

where Λ m is the largest value of C kk = kP (k |k).

It is difficult to find the exact solution for Λ m for a general form of P (k |k). But it is possible to extract the value of the epidemic threshold. In the case of uncorrelated networks, P (k |k) = k P (k )/ k and C kk = k kP (k )/ k . So, it is easy to check that v k = k is an eigenvector with eigenvalue k 2 / k that, according to Perron-Frobenius theorem, is the largest. Thus, we obtain the epidemic threshold:

Equation (21) has strong implications since several real networks have a power law degree distribution P (k) ∼ k −γ with exponents in the range 2 < γ < 3 [28] . For these distributions, the second moment k 2 diverges in the limit of infinite sizes implying a vanishing threshold for the SIS model or, equivalently, the epidemic prevalence for any finite infection rate. Both theories HMF and QMF predict vanishing thresholds for γ < 3 despite of different scaling for 5/2 < γ < 3. However, HMF predicts a finite threshold for networks with γ > 3 unlike the QMF theory that still predicts a vanishing threshold [47] .

For the CP model, the heterogeneous mean-field equation, analogous to Eq. (18), can be written as

Again we assume degree-uncorrelated networks and a simple linear stability analysis shows the presence of a phase transition, located at the value λ c = 1, as found in homogeneous mean-field theory, independent of the degree distribution and degree correlations. According to simulations [19, [21] [22] [23] [24] 36 ] the transition point does not quantitatively reproduced the predictions of both approaches, HMF and QMF. However, the advantage of HMF over the QMF theory, is that we can analytically obtain the critical exponents for the dynamical model and compare with numerical results. Indeed, some works have shown that the contact process running on the top of highly heterogeneous random networks is well-described by the heterogeneous mean-field theory [22, 23] . However, some important aspects such as the threshold and strong corrections to the finitesize scaling observed in simulations are not clarified in this theory. We summarized the intense scientific discussion [19, 20, 22-24, 50, 51] about this subject in section VI.

Improvement of both HMF and QMF theories including dynamical correlations by means of a pairapproximation do not change qualitatively the results, however they promote a quantitative refinement, as we showed hereinafter.

In the paper [16] , the authors investigated the role of dynamical correlations on the dynamic of the SIS epidemic model on different substrates. We can start rewriting Eq. 9 as follows:

where A ij is the adjacency matrix and φ ij represents the probability of a pair of nodes i and j of being in the state φ ij = [S i I j ], this means the node i is susceptible and its neighbor j is infected. When we uses a simple approximation, this joint probability was factorized:

The first three terms represents the processes related to the pair of neighbors i and j, spontaneous annihilation reactions

The other terms represent processes related to the interaction with the other neighbors of i and j, this means another vertex l that can infect i or j: (23) and (24) cannot be solved due to the triplets. In turn, the dynamical equations for triplets will depend on quadruplets, and so on. Then we have to break these correlations in some point. To obtain an pair-approximation solution, we can apply the standard pair-approximation [1, 52] :

After a few steps (details are in reference [16] ), similar to what we did in the one-vertex mean-field approximaation as to perform a linear stability analysis around the fixed point ρ i = [S i I j ] = [I i I j ] = 0 and a quasi-static approximation for for t → ∞, dρ i /dt ≈ 0 and dφ ij /dt ≈ 0, in Eqs. (23) and (24), we can find the Jacobian matrix

As in the case of one-vertex QMF theory, the critical point is obtained when the largest eigenvalue of L ij is nul. Analytical solution for simple networks as random regular networks, star an wheel graphs can be directly obtained from Eq. 26. For power law arbitrary random networks, as UCM model, the largest eigenvalue of Eq. 26 can be numerically determined. In reference [16] , Mata and Ferreira showed that the thresholds obtained in pair and one-vertex QMF theories have the same scaling with the system size but the pair QMF theory is quantitatively much more accurate than the one-vertex theory when compared with simulations. In reference [38] the authors also studied the impact of dynamic correlations on the SIS dynamics on static networks. Performing the same analysis for the contac process, one obtains the Jacobian matrix:

.

The critical point is also obtained when its largest eigenvalue is null. A general analytical expression is not available for large random power-law degree networks but, in principle, it can be obtained numerically for any kind of network. However, for simple graphs as a random regular network, the transition point can be obtained after some algebraic manipulations and it is given by:

where m is the degree of all nodes of the network. That is the same value yield by the simple homogeneous pairapproximation [2] .

For the SIS model, the equation for the probability that a vertex with degree k is occupied takes the form

where the conditional probability P (k |k), which gives the probability that a vertex of degree k is connected to a vertex of degree k , weighs the connectivity between compartments of degrees k and k and φ kk represents a pair of nodes with degree k and k respectively, in the state [S k I k ]. The dynamical equation for φ kk is

The one-vertex mean-field equation [Eq. (18) ] is obtained factoring the joint probability φ kk ≈ (1 − ρ k )ρ k in Eq. (29) . The factor k − 1 preceding the first summation in Eq. (30) is due to the k neighbors of middle vertex except the link of the pair [0 k 0 k ] (similarly for k − 1 preceding the second summation).

Following the same line of reasoning as the previous calculations, we now approximate the triplets in Eq. (30) with the pair-approximation of Eq. (30), performing a linearization around the fixed point ρ k ≈ 0 and φ kk ≈ 0, and performing a quasi-static approximation for t → ∞, in which dρ k /dt ≈ 0 and dφ kk /dt ≈ 0.

Considering uncorrelated random networks, we obtain the Jacobian L kk

with δ kk being the Kronecker delta symbol. Again, the absorbing state is unstable when the largest eigenvalue of L kk is positive. Therefore, the critical point is obtained when the largest eigenvalue of the Jacobian matrix is null, thus obtaining:

This threshold coincides with that of the susceptibleinfected-recovered (SIR) model in a one-vertex HMF theory [27] . This results was also proposed in Ref. [18] using heuristic arguments. They argued that, dynamical correlations, that are neglected in a one-vertex HMF theory, account for the fact that infected nodes have higher probability to be still infected. This means that, in the next step, this node can be considered immunized (recovered for a while). So, a better upper bound for the spreading of the disease is given by the HMF theory of the SIR model, that is exactly the threshold given by Eq. (32) .

In a similar way, we can perform the same analysis for the CP model and we obtain the Jacobian:

Assuming that the network does not present correlated degree, we have P (k |k) = k P (k )/ k . So, the u k = k is an eigenvector of C kk with eigenvalue

Since C kk > 0 is irreducible and u k > 0, the Perron-Frobenius theorem [44] warranties that Λ is the largest eigenvalue of C kk . The critical point is given by −1+Λ = 0, it means we have to solve the transcendent equation

numerically to obatin the transition point for any kind of uncorrelated degree network. Using a random regular network as an example, this means, P (k) = δ k,m , where m is the degree of all nodes of the network, we obtain again:

that is the same of the homogeneous and quenched pair mean-fiel theory.

In reference [53] the authors have also determined the critical exponents in the pair HMF approach for the CP model. For the infinite size limit the exponents are the same as the one-vertex theory, as expected. However, the finite-size corrections to the scaling obtained in the pair HMF theory allowed a remarkable agreement with simulations for all degree exponents (2.0 ≤ γ ≤ 3.5) investigated, suppressing a deviation observed for low degree exponents in the one-vertex HMF theory [23] .

Numerical simulations is an essential tool to predict the accuracy of mean-field approaches in the study of epidemic processes on complex networks. Although this tool is widely used, strict implementations of epidemic processes on networks with high heterogeneity on degree distribution are not simple. In reference [54] , the authors showed that, depending on the network properties, the threshold of the SIS model can be altered when occur modifications in the SIS dynamics, even preserving the basic properties of spontaneous healing and potencial of infection of each vertex growing unlimitedly with its degree.

The classical algorithm to model continuous-time Markov processes is known as Gillespie algorithm [55] . In this recipe, we associated each dynamical transition (infection and healing for the SIS or CP model, for example) with a Poisson process, this mean, independent spontaneous processes. At each change of state, we have to update a list containing all possible spontaneous processes. However, for very large networks, this is computationally unfeasible. So, we used an optimized Gillespie algorithm proposed by Cota and Ferreira [56] . In the following sections, we summarized these optimized algorithms for the SIS and CP dynamics. For details of how to implement optimized algorithms of continuous-time Markovian processes based on Gillespie algorithm see reference [56] .

The SIS dynamics in a network of size N can be simulated in a very simple way: Select a vertex at random with equal chance. If the selected vertex i is infected we turn it to susceptible with probability

where k max is maximal number of connections, n i = j A ij σ j is the number of infected nearest neighbors of the vertex i, and σ j = 1 corresponds to infected node and σ j = 0 otherwise. Here, we are considering the simplest case of λ i = λ and µ i = µ for all vertices. If the selected vertex i is susceptible, it becomes infected with probability

This algorithm is accurate and can used for any generic SIS dynamics. However, if one is interested in regions close to the threshold where the great majority of the vertices are susceptible, the algorithm is very inefficient since changes happens only in the neighborhood of infected vertices. Therefore, we can use a more efficient strategy based on the previous algorithm. This strategy requires to keep and constantly update a list P with the positions of all infected vertices where changes will take place. The list update is simple. The position of a new infected is added at the end of the list. When a infected vertex becomes susceptible, the last entry of the list is moved to the index of the cured vertex.

The total rate that a infected vertex becomes susceptible in the whole network is R = µN i , where N i is the number of infected vertices. Analogously, the total rate that one susceptible vertex is infected is given by J = λN e , where N e is the number of vertices emanating from infected nodes. So, the SIS dynamics can be simulated according to the algorithm proposed by Ferreira et al. [12] as follows: the step is incremented by ∆t = 1/(R + J). With probability p = R/(R + J) an infected vertex i is selected randomly and turns it to susceptible. With complementary probability q = J/(R+J) an infected vertex is selected at random and accepted with probability proportional to its degree. In the infection attempt, a neighbor of the selected vertex is randomly chosen and if susceptible, it is infected. Otherwise nothing happens and simulations run to the next time step.

The CP dynamics can also be efficiently simulated if a list of occupied vertex P is used analogously to the SIS algorithm. The total rate of cure is also given by R = µN i . The total creation rate is J = λN i [2] . An infected vertex i is selected with equal chance. With probability p = R/(J + R) = µ/(µ + λ) it is cured. With probability q = J/(J + R) = λ/(µ + λ) one of the k i neighbors of i is selected and, if susceptible, is infected. The time is incremented by ∆t = 1/[N i (λ + µ)]. Notice that infected vertices are selected independently of their degrees and with probability n i /k i reach an already infected neighbor.

Although equivalent for strictly homogeneous graphs (k i ≡ k ∀ i), the SIS and the CP models are very different for heterogeneous substrates. The universality class of CP and SIS is the same in homogeneous lattices. Both models belong to the directed percolation universality class [6] . Nevertheless, in complex networks, heterogeneity affects both models and, at the heterogeneous meanfield approach, they have different critical exponents (see discussion in Ref. [57] ).

In both SIS and CP models, the control parameter of the dynamics is the infection rate λ. In the thermodynamic limit, above a critical value λ c (epidemic threshold), a finite fraction of the population is infected. However, for λ < λ c , the epidemic can not survive and the dynamics goes to an absorbing state where everyone is susceptible. Nevertheless, for a finite system the unique achievable stationary state is the absorbing state, even above the critical point, because of dynamical fluctuations. Some simulation methods were proposed in the literature to solve the issue of study the active state in finite systems , as we will see in the next section.

Mean-field approaches and field theory renormalization are the key analytical tools to investigate dynamical processes with absorbing-state phase transition [58] . While the former is only valid above the upper critical dimension, application of the latter in physical dimensions is hindered by large technical difficulties. For this reason, most of our knowledge about the properties of absorbing-state phase transitions is based in the computer simulation of different representative models.

The numerical analysis of these computer data also represents a challenge mainly because of finite size effects. In finite systems, any realization of the dynamics reaches the absorbing state sooner or later, even above the critical point, due to dynamic fluctuations. This difficulty was traditionally overcome by starting with a finite initial density of active sites and averaging only over surviving samples, i.e., realizations which have not yet fallen into the absorbing state [2, 59] . Analyzing the quasistationary state defined by surviving averages, the critical point and various critical exponents can be performed by studying the decay of the survival average of different observables as a function of the system size. However averaging over surviving runs is computationally so inefficient since surviving configurations are increasingly rare at long times.

A much more efficient strategy is provided by the quasistationary (QS) method, proposed by de Oliveira and Dickman [32, 60] , in which every time the system is on the verge of fall in an absorbing state, it jumps to an active configuration previously stored during the simulation. This can be computationally implemented by saving, and constantly updating a sample of the states already visited. The update is done by periodically replacing one of these configurations by the current one.

The characteristic relaxation time is always short for epidemics on random networks due to the very small average shortest path [44] . The averaging time, on the other hand, must be large enough to guaranty that epidemics over the whole network was suitably averaged. It means that very long times are required for very low QS density (sub-critical phase) whereas relatively short times are sufficient for high density states [31, 32] .

Both equilibrium and non-equilibrium critical phenomena are hallmarked by simultaneous diverging correlation length and time which microscopically reflect divergence of the spatial and temporal fluctuations [1] , respectively. Even tough a diverging correlation length has little sense on complex networks due to the small-world property [61] , the diverging fluctuation concept is still applicable. We used different criteria to determine the thresholds, relied on the fluctuations or singularities of the order parameters.

The QS probabilityP n , that does not depend on the initial condition, is defined as the probability that the system has n occupied vertices in the QS regime, is computed during the averaging time and basic QS quantities, as lifespan and density of infected vertices, are derived fromP n . Indeed, we have that ρ = 1 N nP n and τ = 1/P 1 [32] , where τ is the lifespan of the epidemic.

Thus, thresholds for finite networks can be determined using the modified susceptibility [12] 

that does exhibit a pronounced divergence at the transi-tion point for SIS [12, 14, 16] and contact process [53, 62] models on networks. The choice of the alternative definition, Eq. (39), instead of the standard susceptibilitỹ χ = N ( ρ 2 − ρ 2 ) [6] is due to the peculiarities of dynamical processes on complex networks. In a finite system of size N , χ shows a diverging peak at λ = λ QS p (N ), providing a finite size approximation of the critical point. In the thermodynamic limit, λ QS p (N ) approaches the true critical point with the scaling form [63] 

as we can see for the SIS model in Figure 2 (a). The network was generated with the uncorrelated configuration model [25] , where vertex degree is selected from a power-law distribution P (k) ∼ k −γ with a lower bound k 0 = 3.

In the context of epidemic modeling on complex networks [64] , Boguñá et. al. [18] proposed another strategy which considers the lifespan of spreading simulations starting from a single infected site as a tool to determine the position of the critical point. Each realization of the dynamical process is characterized by its lifespan and its coverage C, where latter is defined as the fraction of different sites which have been occupied at least once during the realization. In the thermodynamic limit realizations can have either finite or infinite, according to whether they proceed below or above the critical point. Endemic realizations have an infinite lifetime and their coverage is equal to 1. Finite realizations have instead a finite lifetime and a coverage vanishingly small in the limit of diverging size.

In finite systems this distinction is blurred, since any realization is bound to end, reaching the absorbing state, although this can occur over long temporal scales. In practice, a realization is assumed as active whenever its coverage reaches a predefined threshold value 4 C th , which was generally takes equal to C th = 0.5. Realizations ending before value C = C th is reached are considered to be finite.

In this method the role of the order parameter is played by the probability Prob(λ − λ c , N ) that a run is longterm, while the role of susceptibility is played by the average lifetime of finite realizations τ . For small values of λ all realizations are finite and have a very short duration τ . As λ grows the average duration of finite realizations increases, but for very large λ almost all realizations are long-term, only very short realizations remaining finite. For this reason τ exhibits a peak for a value λ p (N ) depending on N and converging to λ c in the thermodynamic limit. We can then use the average lifespan to determine numerically the critical point, as shown in Figure 2 (b). As we observe in Figure 3 , the critical points as a function of network size obtained by both methods are in very well agreement. In reference [67] , Sander and collaborators sumarize alternative methods for work around the problem related to the absorbing phase in simulations. They mentioned the reflecting boundary condition [68] that consists basically in avoiding the absorbing state by reverting the system to the configuration that it was immediately before visit the absorbing state. Other strategy is to use a uniform external field that creates particles spontaneously at a given rate which disappears in the thermodynamic limit [69] . In addition, one can use the hub reactivation method on heterogeneous networks. If the system reaches the absorbing state, the dynamics starts again with the most connected vertex of the network infected.

In Ref. [24] , Castellano and Pastor-Satorras derived the HMF theory for the CP dynamic in the limit of infinite network size. They obtained the following scalinḡ

At the transition point λ = λ c ,

Also the relaxation time scales as

These exponents are also obtained using a pair HMF approximation as shown in reference [53] . It is not possible to check this predictions with numerical analysis because of the finite-size effects. A comparison became possible using the finite-size scaling (FSS) ansatz [59] , adapted to the network topology, and they previously concluded that CP dynamics on networks was not described by the HMF approximation. However, it was assumed in Ref. [24] that heterogeneous networks follow the same FSS known for regular lattices [2] . Indeed, the FSS on networks is more complicated than previously assumed. The behavior of the CP on networks of finite size depends not only on the number of vertices N but also on the moments of the degree distribution [19] . This implies that, for scale-free networks, the scaling around the critical point depends explicitly on how the largest degree k c diverges with the system size N . Such dependence introduces very strong corrections to scaling. However, when such corrections are suitably taken into account, they showed that the CP on heterogeneous networks agrees with the predictions of HMF theory with good accuracy [22] . Ferreira and collaborators [23] started from Eq. (22) , and they considered, in addition, uncorrelated networks with P (k|k ) = kP (k)/k, to obtain the equation for the overall density ρ = k ρ k P (k):

A mean field theory for the FSS can be obtained using the strategy proposed by Castellano and Pastor-Satorras [19] , in which the motion equation is mapped in a one-step process, in the limit of very low densities, with transition rates

where W (n, m) represents the transitions from a state with m infected vertices to another state with n infected vertices. In the stationary state, dρ(t)/dt = 0, the Eq. (44) read as [22] ρ k = λkρ/ k 1 + λkρ/ k .

Close to the criticality, when the density at long times is sufficiently small such thatρk c 1, Eq. (46) becomes ρ k λkρ/ k . Substituting this result in Eq. (45) , one finds that the first-order approximation for the one-step processes is

where g = k 2 / k 2 . The master equation for a standard one-step process is [70] dP n dt = m W (n, m)P m (t) − m W (m, n)P n (t). (48) Substituting the rates (47), we find dP n dt = (n + 1)P n+1 + u n−1 P n−1 − (n + u n )P n (49) with u n = λn(1 − ng). Since the probability for the process not to end up in the absorbing state up to time t, named survival probability, is given by P s (t) = n≥1 P n (t), we can define the quasistationary (QS) dis-tributionP n as [2] P n = lim

withP 0 ≡ 0 and normalized condition n≥1P n = 1 (see more details in section V C). The solutions of the equation (49) have already been exhaustively investigated, then we merely report the results of Ref. [22] where it was found that the critical QS distribution for large systems has the following scaling form 5

where f (x) is as scaling function with the following prop-

The critical quasistationary density scale as

Similarly, the characteristic time scales as

For a network with degree exponent γ and a cutoff scaling with the system size as k c ∼ N 1/ω , where ω = max[2, γ − 1] for uncorrelated networks with power law degree distribution [25] , the factor g scales for asymptotically large systems as g ∼ k 3−γ c for γ < 3 and g ∼ const. for γ > 3. The result is a scaling law ρ ∼ N −ν and τ ∼ Nα where the exponentsν andα are given bŷ

(54) In Ref. [23] , Ferreira et al. investigated the CP on heterogeneous networks with power-law degree distribution by performing quasistationary simulations, and concluded that heterogeneous mean-field theory correctly describes the critical behavior of the contact process on quenched networks. However, some important questions remained unanswered. The transition point λ c = 1 predicted by this theory does not capture the dependence on the degree distribution observed in simulations. Subleading corrections to the finite-size scaling, undetected by the one-vertex HMF theory, are quantitatively relevant for the analysis of highly heterogeneous networks (γ → 2), for which deviations from the theoretical finitesize scaling exponents were reported [23] .

The HMF theory assumes that the number of connections of a vertex is the quantity relevant to determine its state and neglects all dynamical correlations. But in reference [53] , the authors present a pair HMF approximation, the simplest way to explicitly consider dynamical correlations, for the CP on heterogeneous networks. Despite they found the same critical exponents obtained in the simple HMF approximation, the corrections of the finite-size scaling were better, supporting that degree based theory estimate correctly the scaling exponents of the contact process on scale-free networks.

As we saw in the previous sections, distinct theoretical approaches were devised for the SIS and CP models to determine an epidemic threshold λ c separating an absorbing, disease-free state from an active phase [10-12, 14-16, 18, 71] . The quenched mean-field (QMF) theory [71] explicitly includes the entire structure of the network through its adjacency matrix while the heterogeneous mean-field (HMF) theory [9, 72] performs a coarsegraining of the network grouping vertices accordingly their degree. However, for the SIS model, both theories predicts different thresholds. The HMF theory predicts a vanishing threshold for the range 2 < γ ≤ 3 while a finite threshold is expected for γ > 3. Conversely, the QMF theory states a threshold inversely proportional to the largest eigenvalue of the adjacency matrix, implying that the threshold vanishes for any value of γ [10] . Regardless, Goltsev et al. [11] proposed that QMF theory predicts the threshold for an endemic phase, in which a finite fraction of the network is infected, only if the principal eigenvector of adjacency matrix is delocalized. In the case of a localized principal eigenvector, that usually happens for large random networks with γ > 3 [73] , the epidemic threshold is associated to the eigenvalue of the first delocalized eigenvector. For γ < 3, there exists a consensus for SIS thresholds: both HMF and QMF are equivalent and accurate for γ < 2.5 while QMF works better for 2.5 < γ < 3 [12, 16] .

Lee et. al. [14] proposed that for a range λ QM F c < λ < λ c with a nonzero λ c , the hubs in a random network become infected generating epidemic activity in their neighborhoods. This activity has a characteristic lifespan τ (k, λ) depending on the degree k and the infection rate λ. On networks where almost all hubs are directly connected the activity can be spread among them if the lifespan τ (k, λ) is large enough. Then, above λ QM F c , the network is able to sustain an endemic state due to the mutual reinfection of connected hubs.

However, when hubs are not directly connected, the reinfection mechanism does not work and high-degree vertices produce independent active domains. These independent domains were classified as rare-regions, in which activity can last for very long periods increasing exponentially with the domain size [74] . This means, usually we have two distinct states: λ > λ c corresponds to a supercritical phase where the system is globally active and λ < λ c corresponds to an absorbing inactive state. However, the SIS dynamics running on top of power-law networks presents a region in which λ is smaller than the epidemic threshold -but greater than a certain value be-low which the epidemic actually ends -where the activity survives for very long times. This results in a slow dynamics known as Griffiths phase [36, 75, 76] . The sizes of these active domains increase for increasing λ leading to the overlap among them and, finally, to an endemic phase for λ > λ c . In the thermodynamic limit these regions vanish because they decreases as soon as the network size increases [17, 73, 77] .

This anomalous behavior in the subcritical phase was also investigated in reference [37] . The authors used extensive simulations to show that the SIS model running on the top of power-law networks with γ > 3 can exhibit multiple peaks in the susceptibility curve that are associated with large gaps in the degree distribution among the few most highly connected nodes, which permits the formation of these independent domains of activity. However, if the number of hubs is large, as occurs for networks with γ < 3, the domains are directly connected and the activation of hubs implies in the activation of the whole network.

The arguments presented by the authors of reference [37] are in agreement with the scnario investigated in refs. [11, 14] that leads to the conclusion that the threshold to an endemic phase is finite in random networks with a power law degree distribution for γ > 3. Inspired in the appealing arguments of Lee et al. [14] , Boguñá, Castellano and Pastor-Satorras [18] reconsidered the problem and proposed a semi-analytical approach taking into account a long-range reinfection mechanism and found a vanishing epidemic threshold for γ > 3.

As reported by Lee et. al. [14] , when the hubs on a network are directly connected, the activity can be spread throughout the network even in the limit λ → 0. However, when higher degrees nodes are distant from each other these hubs are able to sustain local active domains around them and only with a nonzero λ c , the endemic state is reached. Nevertheless, Boguñá, Castellano and Pastor-Satorras [18] revisited the problem taking into account long range dynamic correlations in a coarse-grained time scale. As explained in Ref. [18] , their approach states that a directed connection between hubs is not a necessary condition for reinfected them. There is a possibility of "long-range" reinfection since the network has a small-world property. In this approach, it was concluded that the epidemic threshold vanishes for random networks with P (k) decaying slower than exponentially, in particular, P (k) ∼ k −γ with any γ. It was rigorously proved by Chatterjee and Durret [42] in the thermodynamic limit, for networks with degree distribution P (k) ∼ k −γ with γ > 3.

Afterward, Mountford and collaborators [40] expanded the result found in reference [42] including the range 2 < γ ≤ 3 of the degree exponent. They also analysed the behavior of the density of infected nodes in function of λ close to the epidemic threshold and they predicted analytical exponents which were also found in the numerical results of reference [41] .

Recently, Castellano and Pastor-Satorras [39] enlight-ened the understanding of the SIS dynamics elaborated a mathematical formulation of the mutual reinfection process, using the cumulative merging percolation (CMP) process proposed by Menard and Singh [78] . In this process, each node can be considered active with a certain probability. Inactive nodes play just as a bridge between active ones. A initial cluster with size 1 contains only an active node but, in a interactive process, two clusters can colapse into one if the criterion of topological distance between them is satisfied. So the CMP process creates a cluster composed by a set of active nodes that were aggregated due to iteration of merging events. Such nodes are part of the same connected component of the underlying network. The insight of Castellano and Pastor-Satorras [39] were classified the hubs able to sustain the epidemic as these active nodes of the CMP process. Therefore they could relate this process to the reactivation of hubs that are not directly connected. Consequently, they observed that the presence of a CMP giant component is related to an endemic stationary state. In their paper, they showed that the epidemic threshold does not behavior as QMF prediction but it vanished more slowly, with an exponent that decreases as soon as γ increases. The dependence of the epidemic threshold with the network size that they found is in agreement with the asymptotic scaling found analytically by Mountford and collaborators [40] and recently by Huang and Durret [79] .

In this paper, we have reviewed the main features of the SIS and CP models, which have been very used to describe epidemic dynamics on complex networks. Throughout this review we described both models, presented distinct theoretical approaches devised for them, their advantages and disadvantages, and the main differences among them. We also exposed simulation techniques to analyze both models numerically. Since both models are examples to investigate absorbing phase transitions in complex networks, we presented the main simulation strategies to overcome the difficulty to study the active state of finite networks. Finally we reported the central issue for each model and we summarized the difficulties and discussions that came up in the literature related to these issues. For SIS model, problems related to determine the epidemic threshold on heterogeneous networks. For CP model, concerns related to the critical exponents and the degree distribution of the network.

Although there are some books and articles reviewing such models, the main idea of this manuscript is to provide, in summary, an overview of the main points of this subject: the theories and simulation techniques, as also the main concerns about the investigation of epidemic models with phase transitions to absorbing states running on top of complex networks.

The progress in this area grows incredibly fast and it is not possible to discuss all recent results. But we try to mention just a few of many research lines. In the last decades, epidemic models have also been studied in hypergraphs [80] , temporal networks [81] [82] [83] , metapopulations [84, 85] and also multiplex subtrates [86, 87] . They analyzed, among other issues, epidemic spreading with awareness, social contagion, measures of epidemic control and how patterns of mobility affects the transmission of the disease. There are also studies about the impact of infectious period or recovery rates on epidemic spreading [88] [89] [90] , spectral properties of epidemic in correlated networks [41] , the speed of disease spreading [91] . It is also important to mention the challenges in modelling spreading diseases related to public health, global transmission [92, 93] as well livestock and vector-borne diseases [94, 95] . The relevance of studying epidemic models is also evident when faced with alarming situations such as the recent pandemic of COVID-19 caused by the new coronavirus [96] [97] [98] . These references and the others cited throughout the manuscript provide accurate studies for readers who wish to go deep into the subject.

Non-equilibrium phase transitions

Nonequilibrium phase transitions in lattice models

Mathematical epidemiology of infectious diseases: model building, analysis and interpretation

Infectious diseases in humans

Nonequilibrium phase transition: Absorbing Phase Transitions

Complex Webs in Nature and Technology

Dynamical Processes on Complex Networks

The mathematical theory of infectious diseases and its applications

Networks: An Introduction

Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers

Proc. Natl. Acad. Sci. USA

Critical Dynamics

of Current Physics-Sources and Comments

Monte Carlo simulation in statistical physics: an introduction

Stochastic Processes in Physics and Chemistry

This work was partially supported by CAPES, FAPEMIG and CPNq. The author thanks the financial support from FAPEMIG (Grant No. APQ-02482-18) and CNPq (Grant No. 423185/ 2018-7). The author also wishes to express her deepest gratitude for Silvio C. Ferreira for reading this overview and providing relevant suggestions.