Burstiness in activity-driven networks and the epidemic threshold

TL;DR

This paper analyzes how burstiness in activity-driven networks influences epidemic spreading, deriving an analytical epidemic threshold and revealing that burstiness lowers the threshold but has complex effects on prevalence depending on infectivity.

Contribution

It provides a closed-form analytical expression for the epidemic threshold considering burstiness in activity-driven networks, clarifying its dual effects on epidemic prevalence.

Findings

Burstiness lowers the epidemic threshold.

In low-infectivity systems, burstiness increases infection probability.

In high-infectivity systems, burstiness weakens epidemic spreading.

Abstract

We study the effect of heterogeneous temporal activations on epidemic spreading in temporal networks. We focus on the susceptible-infected-susceptible (SIS) model on activity-driven networks with burstiness. By using an activity-based mean-field approach, we derive a closed analytical form for the epidemic threshold for arbitrary activity and inter-event time distributions. We show that, as expected, burstiness lowers the epidemic threshold while its effect on prevalence is twofold. In low-infective systems burstiness raises the average infection probability, while it weakens epidemic spreading for high infectivity. Our results can help clarify the conflicting effects of burstiness reported in the literature. We also discuss the scaling properties at the transition, showing that they are not affected by burstiness.

Equations

$$\rho (a_i)=\nu a_m^{\nu }{a}_i^{-(\nu +1)}$$

$${\Psi }_{a_i}(\tau )=\alpha {\xi }_i^{\alpha }{\tau }^{-(\alpha +1)}$$

$${\partial }_t{Z}_a(t,t^{\prime })=-\mu Z_a(t,t^{\prime })+\lambda (1-Z_a(t,t^{\prime })){\int }_0^{\infty }d{a}^{\prime }\rho (a^{\prime })a^{\prime }{Q}_{a^{\prime }}(t)$$

$$Z_a(t^{\prime },t^{\prime })=Q_a(t^{\prime })+\lambda [1-Q_a(t^{\prime })]{\int }_0^{\infty }d{a}^{\prime }\rho (a^{\prime })P_{a^{\prime }}(t^{\prime })$$

$$Z_a(t,t^{\prime })=\left[Q_a(t^{\prime })+\lambda \overline{P}(1-Q_a(t^{\prime }))-\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}\right]e^{-(\mu +\lambda \overline{aQ})(t-t^{\prime })}+\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}$$

$$Q_a(t)$$

$$P_a(t)$$

$$Q_a(t)$$

$$\begin{array}{rl}\overline{P}& =h(\overline{P},\overline{aQ})=\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}+{\int }_0^{\infty }da\rho (a)a\frac{1-L_a(\mu +\lambda \overline{aQ})}{(\mu +\lambda \overline{aQ}{)}^2}\frac{\lambda \mu \overline{P}}{1-(1-\lambda \overline{P})L_a(\mu +\lambda \overline{aQ})}\\ \overline{aQ}& =f(\overline{P},\overline{aQ})={\int }_0^{\infty }da\rho (a)a\frac{\left[\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}+\left(\lambda \overline{P}-\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}\right)L_a(\mu +\lambda \overline{aQ})\right]}{1-(1-\lambda \overline{P})L_a(\mu +\lambda \overline{aQ})}\end{array}$$

$$\left[\begin{array}{c}{\overline{P}}_{n+1}\\ {\overline{aQ}}_{n+1}\end{array}\right]\sim J\left[\begin{array}{c}{\overline{P}}_n\\ {\overline{aQ}}_n\end{array}\right]=\left[\begin{array}{cc}\frac{\lambda ⟨a⟩}{\mu }& \frac{\lambda }{\mu }\\ \lambda {\int }_0^{\infty }da\rho (a)a\frac{L_a(\mu )}{1-L_a(\mu )}& \frac{\lambda ⟨a⟩}{\mu }\end{array}\right]\left[\begin{array}{c}{\overline{P}}_n\\ {\overline{aQ}}_n\end{array}\right]$$

$$\frac{\lambda ⟨a⟩}{\mu }+\lambda \sqrt{\frac{1}{\mu }{\int }_0^{\infty }da\rho (a)a\frac{L_a(\mu )}{1-L_a(\mu )}}<1.$$

$$\frac{\lambda ⟨a⟩}{\mu }=\frac{⟨a⟩}{⟨a⟩+\sqrt{⟨a^2⟩}}$$

$$\lambda r+\lambda \sqrt{\nu {\left(\frac{r(\nu -1)}{\nu }\right)}^{\nu }{\int }_{\frac{r(\nu -1)}{\nu }}^{\infty }dx\frac{1}{x^{\nu }}\frac{\alpha {\left(\frac{\alpha -1}{\alpha x}\right)}^{\alpha }\Gamma \left(-\alpha ,\frac{\alpha -1}{\alpha x}\right)}{1-\alpha {\left(\frac{\alpha -1}{\alpha x}\right)}^{\alpha }\Gamma \left(-\alpha ,\frac{\alpha -1}{\alpha x}\right)}}<1$$

$$P_d(t,L)\sim t^{-\delta -1}\mathcal{F}(t{L}^{-z})$$

$${\tau }_{life}(\Delta ,N)\sim \ln (N)\mathcal{G}(\Delta N^{\gamma })$$

$$\overline{P}\sim (r-r_C{)}^{\beta }$$

$$\begin{array}{rl}{\tilde{Q}}_a(s)={\int }_0^{\infty }dt{e}^{-st}{Q}_a(t)=& \frac{1}{s}\left[\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}+\left(\lambda \overline{P}-\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}\right)L_a(\mu +\lambda \overline{aQ})\right]\\ & +(1-\lambda \overline{P}){\tilde{Q}}_a(s){\int }_0^{\infty }d{\tau }^{\prime }{e}^{-(\mu +\lambda \overline{aQ}+s){\tau }^{\prime }}{\Psi }_a({\tau }^{\prime })+{\epsilon }_a\end{array}$$

$$\begin{array}{rl}{\tilde{P}}_a(s)={\int }_0^{\infty }dt{e}^{-st}{P}_a(t)=& \frac{1}{s}\left[\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}+a\left(\lambda \overline{P}-\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}\right)\frac{1-L_a(\mu +\lambda \overline{aQ})}{\mu +\lambda \overline{aQ}}\right]\\ & +a(1-\lambda \overline{P}){\tilde{Q}}_a(s){\int }_0^{\infty }d{\tau }^{\prime }{e}^{-(\mu +\lambda \overline{aQ}+s){\tau }^{\prime }}{\int }_{{\tau }^{\prime }}^{\infty }d{\tau }^{\prime \prime }{\Psi }_a({\tau }^{\prime \prime })\end{array}$$

$${\tilde{Q}}_a(s)=\frac{C_1}{s}+\mathcal{O}(\frac{1}{s})$$

$${\tilde{P}}_a(s)=\frac{C_3}{s}+\frac{C_1}{s}a(1-\lambda \overline{P})\frac{1-L_a(\mu +\lambda \overline{aQ})}{\mu +\lambda \overline{aQ}}+\mathcal{O}(\frac{1}{s})$$

$$C_1=\frac{\left[\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}+\left(\lambda \overline{P}-\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}\right)L_a(\mu +\lambda \overline{aQ})\right]}{1-(1-\lambda \overline{P})L_a(\mu +\lambda \overline{aQ})},L_a(s)={\int }_0^{\infty }d\tau e^{-s\tau }{\Psi }_a(\tau )$$

$$C_3=\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}+a\left(\lambda \overline{P}-\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}\right)\frac{1-L_a(\mu +\lambda \overline{aQ})}{\mu +\lambda \overline{aQ}}$$

$$Q_a(t)\underset{t\to \infty }{\to }{C}_1=Q_a^0,P_a(t)\underset{t\to \infty }{\to }{C}_4=P_a^0$$

$$C_4=C_3+C_{1}a(1-\lambda \overline{P})\frac{1-L_a(\mu +\lambda \overline{aQ})}{\mu +\lambda \overline{aQ}}=C_3+a(1-\lambda \overline{P})\frac{1-L_a(\mu +\lambda \overline{aQ})}{\mu +\lambda \overline{aQ}}{Q}_a^0$$

$$P_a^0=\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}+a\frac{1-L_a(\mu +\lambda \overline{aQ})}{(\mu +\lambda \overline{aQ}{)}^2}\frac{\lambda \mu \overline{P}}{1-(1-\lambda \overline{P})L_a(\mu +\lambda \overline{aQ})}$$

$$Q_a^0=\frac{\left[\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}+\left(\lambda \overline{P}-\frac{\lambda \overline{aQ}}{\mu +\lambda \overline{aQ}}\right)L_a(\mu +\lambda \overline{aQ})\right]}{1-(1-\lambda \overline{P})L_a(\mu +\lambda \overline{aQ})}$$

$$P_a^0(\mu +\lambda \overline{aQ}{)}^2=\lambda \overline{aQ}(\mu +\lambda \overline{aQ})+\lambda \mu \overline{P}\frac{a(1-L_a(\mu +\lambda \overline{aQ}))}{1-(1-\lambda \overline{P})L_a(\mu +\lambda \overline{aQ})}$$

$$Q_a^0(\mu +\lambda \overline{aQ})=\lambda \overline{aQ}+\frac{\mu \lambda \overline{P}{L}_a(\mu +\lambda \overline{aQ})}{1-(1-\lambda \overline{P})L_a(\mu +\lambda \overline{aQ})}$$

$$-\mu \overline{P}+\lambda B\overline{P}-\frac{{\lambda }^2}{\mu }{\overline{aQ}}^2(\overline{P}-1)+\lambda \overline{aQ}-2\lambda \overline{aQ}\overline{P}=0$$

$$\frac{\lambda }{\mu }{\overline{aQ}}^2+\left(1-\frac{\lambda ⟨a⟩}{\mu }\right)\overline{aQ}-\lambda H\overline{P}=0$$

Full text

, , and

Keywords: Network dynamics, Epidemic modelling, Phase transitions into absorbing states, Critical phenomena of socio-economic systems.

1 Introduction

The heterogeneous distribution of times between two consecutive actions – often denoted as burstiness – is a typical signature of time-resolved records of human activities [1, 2]. Mobile phone calls, social contacts and many other complex real-world dynamics exhibit patterns of enhanced activity within short time periods, followed by long lags of inactivity. Actually, the frequency of events shows large variability and strong temporal fluctuations [1, 3, 4, 5] and the event dynamics strongly differs from a Poisson process, where the rate at which events occur is constant. In “bursty” processes, the probability distribution for the inter-event time between two successive actions is not an exponential as in Poisson processes but, typically, it features a fat tail, easily measured from large datasets [1, 2].

Recent studies have addressed the problem of the origin of this highly heterogeneous behavior, arising from the decision-based queuing process that humans apply in performing tasks and allocating their priority [5, 6]. Besides, a large amount of work has been devoted to clarify the effects of intermittent patterns on the temporal structures of interactions, as described by temporal networks [7]. Bursty dynamics can indeed influence the structure of links on the local and on the global scale [8, 9, 10, 11, 12]. More importantly, heterogenous temporal patterns in the evolution of time-varying networks can affect in a non-trivial way dynamical processes. Analytical arguments [13, 14, 15] and numerical approaches [16, 17] have indeed shown that burstiness can significantly modify processes mediated by interactions, such as random walks, epidemics, information diffusion, consensus formation, percolation. Among these, epidemic spreading is one of the most representative and widely applicable example [18]. Analytical results on epidemics in temporal network have focused on the epidemic threshold [19, 20], in general pushed to lower values by bursty effects. In other works, burstiness is introduced on a static network as a non-Markovian effect in infection (or recovery) processes [21, 22, 23, 24, 25].

Interestingly, when comparing bursty with Poisson dynamics in epidemics, conflicting observations have been reported. On the one hand, numerical results [14, 17, 16] and modeling techniques [15, 26, 24, 27] provided strong evidence of a slowing down for the late time spreading in the presence of burstiness, while opposite effects are reported in the early-time dynamics [28]. In [14] and [16], burstiness is found to slow down spreading in early times, while other works observe the opposite effect [26, 17].

To understand such conflicting observations and keep track at the same time of the temporal evolution of interactions, in this paper we focus on the Susceptible-Infected-Susceptible (SIS) process in the presence of burstiness on activity-driven networks [29, 30, 31]. In activity-driven networks, the propensity to engage an interaction is moved from the links to nodes ($N$), by assigning to each node $i$ its activity potential $a_{i}$ measuring the typical number of activations (link formations) per time performed by agent $i$. We then model the bursty behavior by a fat tailed distribution of inter-event times between two successive activations of the same node. Using an activity based mean-field approach, which is here exact in the large $N$ limit, we find a closed analytical form for the epidemic threshold for general activity and inter-event time distributions including, as a particular case, bursty dynamics. The analytical results are in excellent agreement with extensive numerical simulations. We show that, as expected, burstiness lowers the epidemic threshold, while its effect on prevalence is twofold: burstiness tends to raise the average stationary infection probability in low infective systems, therefore strengthening the epidemic, while it weakens the epidemic in high infective systems, lowering the prevalence. This result can help to clarify the conflicting effects of burstiness reported in the literature. We also analyze in details the scaling behavior at the transition, showing that the critical exponents are mean-field exponents of the directed percolation (DP) universality class and that they are unchanged in the presence of burstiness.

The paper is organized as follows. In Section 2 we describe the activity-driven network with bursty dynamics and we define the epidemic process unfolding on top of it. In Section 3 we present the analytical approach and we derive an exact expression for the epidemic threshold, comparing the result with extensive numerical simulations. In Section 4 we analyze the effects of burstiness on disease spreading, focusing on the epidemic threshold as well as on prevalence, and we discuss the scaling properties. Finally, Section 5 contains a summary of our results and some perspectives of our work.

2 Epidemics on activity-driven networks

2.1 Activity-driven temporal networks

We introduce the activity-driven network by associating to each node $i=1,\dots,N$ an activity $a_{i}$ drawn from a distribution $\rho(a_{i})$. The activation rate of node $i$ is defined by its inter-event time distribution $\Psi_{a_{i}}(\tau)$ [8, 9, 10, 11, 12]: this sets the statistic of time intervals between consecutive activations of node $i$, so that the average inter-event time equals the inverse of the activity of node $i$, i.e. $\langle\tau\rangle_{i}=a_{i}^{-1}$. The network evolves according to a Gillespie-like algorithm [32, 22, 33]. We initialize the system at $t=0$, drawing from $\Psi_{a_{i}}(t)$ the first activation time, $t_{i}$, of each node. Then the agent $i$ with the lowest $t_{i}$ becomes active and create $m$ links connecting to $m$ randomly chosen nodes. We fix the next activation time $t_{i}$ for node $i$ as $t_{i}+\tau$ with $\tau$ drawn from $\Psi_{a_{i}}(\tau)$. Finally all links are removed and the process is iterated, activating the node with smallest activation time $t_{i}$.

Empirical measures on real datasets provide a strong evidence of heterogeneous activity patterns, with $\rho(a_{i})$ being a broad, heterogeneous function with large fluctuations [34, 29, 35]. We model such a heterogeneity in the activity distribution as a power law:

[TABLE]

for $a_{i}>a_{m}$ and $\rho(a_{i})=0$ for $a_{i}<a_{m}$, where $a_{m}$ is a lower cut-off. However, as we show in the next sections, our results hold for any $\rho(a_{i})$ and therefore for any activity-driven network.

The inter-event time distribution sets the dynamics of activation for each site. For example, an exponential $\Psi_{a_{i}}(\tau)=a_{i}e^{-a_{i}\tau}$ leads to a Markovian Poisson process as in standard continuous time activation dynamics. When the inter-event time distribution of a single node is heterogeneous, with trains of activation events with short time separation, followed by long inactivity, we take into account the bursty behavior choosing a power-law (Pareto) distribution of inter-event times with a cut-off $\xi_{i}$:

[TABLE]

for $\tau>\xi_{i}$ and $\Psi_{a_{i}}(\tau)=0$ for $\tau<\xi_{i}$. More specifically; for $\alpha>1$, we define $\xi_{i}=\frac{\alpha-1}{\alpha a_{i}}$ so that $\langle\tau\rangle_{i}=a_{i}^{-1}$, and the distribution of $\xi_{i}$ is $\Phi(\xi_{i})\sim\rho(\frac{\alpha-1}{\alpha\xi_{i}})\xi_{i}^{-2}$. Doing so, we can effectively compare systems with the same average inter-event time, fixed by the inverse of the node activity, but with different burstiness distributions $\Psi_{a_{i}}(\tau)$. In the case $0<\alpha<1$, $\langle\tau\rangle_{i}$ diverges and the activity is not properly defined. We can nevertheless introduce the inter-event time distribution according to Eq. (2). As in [9, 8, 11] we set $\xi_{i}=a_{i}^{-1}$, with the cut off $\xi_{i}$ distributed among the nodes according to $\Phi(\xi_{i})=\xi_{i}^{-2}\rho(\xi_{i}^{-1})$. However, in this case $a_{i}$ does not correspond to the activity of the node, i.e. the inverse of average activation time, and non-trivial aging effects may be present [11]. In particular, for the SIS model, since node activation becomes always more unlikely as time increases, we will show that epidemics do not spread in the system.

2.2 The SIS epidemic model on activity-driven networks

We consider the epidemic Susceptible-Infected-Susceptible (SIS) model. Each node can be infected ($I$) or susceptible ($S$-healthy) and the dynamics is based on two elementary processes: contagion, when a susceptible individual, coming in contact with an infected person, has a probability $\lambda$ to contract the infection $S+I\xrightarrow[]{\lambda}2I$; spontaneous recovering, when an infected turns susceptible with a rate $\mu$: $I\xrightarrow[]{\mu}S$. The SIS model on activity-driven networks is defined by a Gillespie-like algorithm as follows [29, 31]):

Up to a relaxation time, $t_{0}$, the network evolves in the absorbing state where all nodes are susceptible, so that the node activation dynamics relaxes to its equilibrium. 2. 2.

At time $t=t_{0}$ we start the epidemics setting the initial condition, for which we divide the population into a configuration of susceptible (S) and infected (I). 3. 3.

Agent $i$ with the lowest activation time $t_{i}$ activates. Infected nodes at time $t$ are recovered at $t_{i}$ with probability $1-e^{-\mu(t_{i}-t)}$, i.e. the recovery is a Poisson process. Then we set $t=t_{i}$. 4. 4.

The active node $i$ generates $m$ links with other randomly-chosen nodes. If the nodes involved in the contacts (links) are both infected ($I$) or susceptible ($S$) nothing happens, otherwise a contagion can occur with probability $\lambda$: $S+I\xrightarrow[]{\lambda}2I$. Here we consider $m=1$, with no loss of generality. 5. 5.

A new inter-event time for the active node $i$ is drawn from $\Psi_{a_{i}}(\tau)$, hence the new activation time $t_{i}$ would be $t_{i}+\tau$. Then all its links are removed and the process is iterated from point 3.

As well known, the SIS model features a phase transition in the class of directed percolation (DP) between an absorbing phase, where the system evolves towards a state with the whole population susceptible (healthy), and an endemic phase, with the epidemics spreading persistently in the population [18, 36]. We consider a control parameter $r$ (e.g. $r=\frac{\lambda}{\mu}$ on static networks, $r=\frac{\lambda\langle a\rangle}{\mu}$ on activity-driven networks with Poisson dynamics [29, 31]), so that the phase transition occurs at a critical value $r_{C}$, the epidemic threshold: if $r<r_{C}$ the system is in the absorbing phase, while if $r>r_{C}$ the system is in the endemic phase. This is the quantity that we will investigate in the next section.

3 The epidemic threshold

3.1 Analytical results

We consider the SIS model on activity-driven networks with arbitrary $\rho(a)$ and $\Psi_{a}(\tau)$. If the average inter-event time $\langle\tau\rangle$ is finite (i.e. $\alpha>1$ considering the power-law distribution of Eq. (2)), nodes can be identified according to their activity, so we divide the population into classes with the same $a_{i}=a$ and study the probability that a node belonging to a class $a$ is infected at time $t$. This corresponds to an activity-based mean-field approach [29, 37], which turns out to be exact for the present system given that connections are continuously reshuffled at each time step, destroying local correlations.

At each contact event occurring at time $t$, we consider the following three quantities to fully characterize the system state:

• $Q_{a}(t)$: the probability for a node in the class $a$, that activates at time $t$, to be infected right before $t$;

• $P_{a}(t)$: the probability for a non-active node, with activity $a$, to be infected at $t$;

• $Z_{a}(t,t^{\prime})$: the probability for a node with activity $a$ to be infected at time $t$, knowing that its last activation occurred at time $t^{\prime}<t$.

The Master Equation for the evolution of $Z_{a}(t,t^{\prime})$ can be written as:

[TABLE]

where first term accounts for the recovery process with rate $\mu$, while the second term describes (for $\alpha>1$) the contagion process; actually, this is the product of the probability to be susceptible at time $t$, $(1-Z_{a}(t,t^{\prime}))$, times the probability to be infected by an active node, which, for $\alpha>1$, can be expressed as $\lambda\int_{0}^{\infty}da^{\prime}\rho(a^{\prime})a^{\prime}Q_{a^{\prime}}(t)$ 111 The probability for node $i$ to infect a node in the interval $[t,t+dt]$ is $\lambda n_{i}(t)Q_{i}(t)/N$, where $n_{i}(t)$ is the number of activations of node $i$ in $dt$. For $\alpha>1$ we have $n_{i}(t)=dt/\langle\tau\rangle_{i}=a_{i}dt$, where $\langle\tau\rangle_{i}$ is the mean time between two activations. Thus, the probability to be contacted by any infected node in $dt$ is $N^{-1}\sum_{i=1}^{N}a_{i}\,dt\,Q_{i}(t)=dt\int da^{\prime}\rho(a^{\prime})a^{\prime}Q_{a^{\prime}}(t)$. Notice that this expression –together with the whole analytic approach– does not apply for $0<\alpha<1$, since the average inter-event time diverges and $a_{i}dt$ is not the probability of activation of node $i$ in the interval $dt$. We observe that in Eq. (3) we only include passive contagion events because, by definition of $Z_{a}(t,t^{\prime})$, the node does not activate between times $t^{\prime}$ and $t$.

Eq. (3) needs to be solved using the initial condition at time $t^{\prime}$:

[TABLE]

where, in analogy with the equation above, the first term corresponds to the case where the active node with activity $a$ is already infected at the contact time $t^{\prime}$; the second term is the probability $(1-Q_{a}(t^{\prime}))$ that the node is susceptible before $t^{\prime}$, then it connects to an infected node with probability $\int_{0}^{\infty}da^{\prime}\rho(a^{\prime})P_{a^{\prime}}(t^{\prime})$, and it is infected with probability $\lambda$.

Eq. (3) with initial condition (4) can be solved in the stationary condition, in the hypothesis that the system reaches an asymptotic steady state $(P_{a}^{0},Q_{a}^{0})$ with $(P_{a}(t),Q_{a}(t))\xrightarrow[t\rightarrow\infty]{}(P_{a}^{0},Q_{a}^{0})$. In this regime, also the node activation dynamics needs to reach the steady state to avoid aging effect in the network evolution. For $\alpha>1$ we expect this condition to be fulfilled if the relaxation time is much larger than the average activation time ($t_{0}\gg\langle a\rangle^{-1}$).

In the stationary regime, one can consider the averages on the activities as independent of time $\int_{0}^{\infty}da\rho(a)aQ_{a}(t)\sim\int_{0}^{\infty}da\rho(a)aQ_{a}^{0}=\overline{aQ}$ and $\int_{0}^{\infty}da\,\rho(a)\,P_{a}(t)\sim\int_{0}^{\infty}da\rho(a)P_{a}^{0}=\overline{P}$ to obtain:

[TABLE]

It is now possible to write explicitly $Q_{a}(t)$ and $P_{a}(t)$ as functions of $Z_{a}(t,t^{\prime})$ using the relations:

[TABLE]

and plugging the solution (5) into Eqs. (6), one obtains:

[TABLE]

where $L_{a}(\mu)=\int_{0}^{\infty}d\tau\,e^{-\mu\tau}\Psi_{a}(\tau)$ is the Laplace transform of $\Psi_{a}(\tau)$. In order to determine the epidemic threshold, one needs to study the stability of the absorbing state solution $(P_{a}(t),Q_{a}(t))=(0,0)\,\forall\,a$ of the Eqs. (7).

Solving Eqs. (7) for $t\rightarrow\infty$, two self-consistency equations for $\overline{P}$, $\overline{aQ}$ are obtained:

[TABLE]

(see A for a detailed derivation). Eqs. (8) can be considered as a two dimensional map: $(\overline{P}_{n+1},\overline{aQ}_{n+1})={\cal M}(\overline{P}_{n},\overline{aQ}_{n})$. The fixed point $(\overline{P}_{n},\overline{aQ}_{n})=(0,0)$ corresponds to the solution $(P_{a}(t),Q_{a}(t))=(0,0)\,\forall\,a$ in Eq. (7). Then, the stability of the absorbing state can be studied by linearizing ${\cal M}(\overline{P}_{n},\overline{aQ}_{n})$ in $(\overline{P}_{n},\overline{aQ}_{n})=(0,0)$:

[TABLE]

The absorbing state is stable if all the eigenvalues of $J$ are smaller than unity, i.e.:

[TABLE]

Eq. (9) provides an analytical form of the epidemic threshold for the SIS model on an activity-driven network ($m=1$) with arbitrary $\rho(a)$ and $\Psi_{a}(\tau)$. As we will show, numerical simulations agree with the analytical estimate, suggesting that the activity based mean-field approach and the stationary hypothesis are exact for this model.

In particular, for an exponential distribution of inter-event times $\Psi_{a}(\tau)=ae^{-a\tau}$ and an arbitrary activity distribution $\rho(a)$, one can get the epidemic threshold from Eq.(9):

[TABLE]

which is the known result for the SIS model on activity-driven networks with standard non-bursty dynamics [29].

In the general case, the parameters $\lambda$, $\mu$ and $\langle a\rangle$ are coupled in a non-trivial way, since Eq. (9) depends on the whole functional form of the activity and inter-event time distributions, so that the control parameter of the phase transition is not easily identified. Hereafter, we will use the adimensional control parameter $r=\frac{\langle a\rangle}{\mu}$. Notice that, unlike the exponential case, the critical value depends in a non-trivial way on $\lambda$.

For a bursty distribution of inter-event times $\Psi_{a}(\tau)\sim\tau^{-(\alpha+1)}$ and the power law activity distribution $\rho(a)\sim a^{-(\nu+1)}$ defined by Eqs. (2) and (1) respectively, Eq. (9) becomes:

[TABLE]

where $\Gamma(s,x)=\int_{x}^{\infty}t^{s-1}e^{-t}dt$ is the upper incomplete Gamma function. In this case, the epidemic threshold, $r_{C}$, can be obtained by solving numerically Eq.(11) for $r$.

3.2 Numerical simulations

We verify with numerical simulations the analytic result (9), checking the validity of our hypotheses. Henceforth, we fix $\lambda=1$, with no loss of generality. We first let the network evolve for a time $t_{0}$ so that network dynamics is relaxed to the equilibrium of the activation dynamics and it has no memory of the initial condition at $t=0$ (when all sites are considered to be active). For $\alpha>1$ this means that $t_{0}$ needs to be much larger than the average activation time $t_{0}\gg\langle a\rangle^{-1}$.

The epidemic model evolves following the steps described in Section 2. The epidemic threshold is numerically computed with the lifetime (or lifespan) method [38, 39]. The lifetime $\tau_{life}$ is the average time that the system takes to reach the absorbing state when infecting at $t_{0}$ a single node of the system, discarding in the average all the endemic realizations which remain infected for an infinite time. In simulations, we consider as endemic a realization whose coverage (number of distinct, ever infected nodes) reaches $N/2$. Observe that fluctuations can cause finite-time recovering even in the endemic phase so that the lifetime as a function of $r$ displays a peak at the critical point $r=r_{C}$ behaving like an effective susceptibility as discussed in [39]. Also, it has been shown that the most effective initial condition to obtain a good estimate of the threshold on an inhomogeneous network is to infect at $t_{0}$ only the most connected site[38]; this is the strategy we adopt here, infecting the most active one.

Figure 1 shows the lifetime as a function of the control parameter for different network sizes $N$, setting all nodes with same activity (i.e. $\rho(a)=\delta(a-a_{m})$). Left and right panels consider an exponential $\Psi_{a}(\tau)\sim e^{-a\tau}$ and a power-law $\Psi_{a}(\tau)\sim\tau^{-(\alpha+1)}$ distribution, respectively. In both cases, for sufficiently large $N$, the numerical data agree very well with the analytical prediction of the threshold.

Figure 2 shows the lifetime as a function of $r$ for a network where $\rho(a)\sim a^{-(\nu+1)}$ and $\Psi_{a}(\tau)\sim\tau^{-(\alpha+1)}$ are power-law distributed according to Eqs. (1) and (2) respectively. Several values of the exponent $\alpha$ have been considered. When $0<\alpha<1$, as we discussed above, analytic results do not apply and we can introduce the inter-event time distribution $\Psi_{\xi_{i}}(\tau)\sim\xi_{i}^{\alpha}\tau^{-(\alpha+1)}$ with a broadly distributed cut-off $\xi_{i}$ (as stated in Section 2). Simulations reveal that the lifetime does not feature a peak since the epidemics always reaches the absorbing state. In particular, the average recovering time $\mu^{-1}$ is finite while node activation becomes always more unlikely as time increases, implying that, eventually, all nodes recover: an endemic or active phase does not exist in this case.

Finally, in Figure 3 we compare numerically and analytically obtained epidemic thresholds in the case of bursty dynamics for two different activity distributions $\rho(a)$, illustrating the excellent agreement between theory and numerics.

4 Burstiness and epidemic spreading

4.1 Epidemic thresholds for various distributions

As explained above, for a generic distribution of inter-event times and activities, we use $r=\frac{\langle a\rangle}{\mu}$ as a control parameter for the phase transition, setting $\lambda=1$. In particular, for a system with $\Psi_{a}(\tau)$ and $\rho(a)$ given by Eqs. (2) and (1) respectively, the critical threshold $r_{C}$ is obtained solving Eq. (11); results are shown in Figure 4. The behavior of $r_{C}$ is qualitatively the same for all $\nu$: it vanishes for $\alpha\rightarrow 1^{+}$ and it saturates to a constant depending on $\nu$ for $\alpha\rightarrow\infty$. For $\alpha<1$, as explained above, the epidemics always ends up in the absorbing state.

As expected, a broader activity distribution leads to a lower epidemic threshold and favors the propagation of the epidemics. In particular, for large $\nu$ the activity distribution tends to a $\delta$-function while in the limit $\nu\to 2$, i.e. when the activity fluctuations diverge, the epidemic threshold vanishes.

Let us now show that this is a general property for any inter-event distribution. Since $\int_{0}^{\infty}\tau\Psi_{a}(\tau)d\tau=a^{-1}$, $\Psi_{a}(\tau)$ tends to a $\delta$-function peaked in $\tau=0$ when $a\to\infty$. Hence, it is possible to expand $e^{-\mu\tau}$ around $\tau=0$ when evaluating $L_{a}(\mu)$ for $a\to\infty$, in this way: $L_{a}(\mu)\xrightarrow[a\to\infty]{}1$ and $1-L_{a}(\mu)=\int_{0}^{\infty}(1-e^{-\mu\tau})\Psi_{a}(\tau)d\tau\xrightarrow[a\to\infty]{}\int_{0}^{\infty}\mu\tau\Psi_{a}(\tau)d\tau=\mu/a$. Plugging these expression into Eq. (9) we obtain that the integral over the activities diverges for $\nu<2$ due to the slow decay at large $a$. Hence, Eq. (9) cannot be satisfied for $\nu<2$ and the absorbing state is never stable. The behavior for $\nu\to 2$ and $\alpha\to 1$ confirms that, keeping the average activity fixed, fluctuations and burstiness lower the epidemic threshold.

In Figure 5, we show the epidemic threshold for a power-law activity distribution Eq. (1) with $\nu=2.5$ and compare several inter-event time distributions $\Psi_{a}(\tau)$. In the left panel the continuous (blue) line is the epidemic threshold $r_{C}$ for inter-event times distributed according to the power law distribution with a cut-off in Eq. (2), Pareto type I distribution. The plot shows that for large $\alpha$ the threshold becomes even higher than in the exponential case, i.e. the dotted (orange) horizontal line. Indeed, for $\alpha\to\infty$, the power law distribution $\Psi_{a}(\tau)$ tends to a $\delta$-function and the node dynamics is periodic [16], hence more regular than the Poisson dynamics (see the green dash-dotted line for comparison in Figure 5 - left panel). The dashed (red) curve in the left panel represents the epidemic threshold $r_{C}$ for the Lomax (or Pareto type II) distribution $\Psi_{a}(\tau)=\alpha\xi_{a}^{\alpha}(\tau+\xi_{a})^{-(\alpha+1)}$ with $\tau>0$: here $\xi_{a}=(\alpha-1)a^{-1}$ so that $\langle\tau\rangle_{a}=a^{-1}$, i.e. we introduce a different lower cut-off. In this case, for $\alpha\to\infty$ we have that $\Psi_{a}(\tau)\to ae^{-a\tau}$, so in this limit the epidemic threshold tends to the value of the Poisson process (dashed orange line). The Lomax distribution is more heterogeneous than an exponential distribution, therefore the epidemic threshold is smaller for all $\alpha$.

Another interesting choice for the inter-event time –usually considered in epidemic spreading– is the Weibull distribution $\Psi_{a}(\tau)=\beta\delta_{a}^{-\beta}\tau^{\beta-1}e^{-(\tau/\delta_{a})^{\beta}}$ with $\tau>0$ [24, 23]: here $\delta_{a}=(a\Gamma(1+1/\beta))^{-1}$ so that $\langle\tau\rangle_{a}=a^{-1}$. The Weibull distribution interpolates between a broad distribution for $\beta\to 0^{+}$ (it tends to the Zipf distribution [24]) and a $\delta$-function for $\beta\to\infty$. Moreover, for $\beta=1$ it reduces to an exponential distribution. Notice that for any $\beta>0$ for large $\tau$ it decays faster than any power-law. In the right panel of Figure 5 we plot $r_{C}$ for the Weibull distribution with a solid (red) line. The threshold vanishes when $\beta\to 0$, it increases with $\beta$ and for $\beta=1$, $r_{C}$ coincides with that of the exponential case (dotted orange line). For $\beta>1$ the threshold increases and for $\beta\to\infty$ it converges to the threshold of a periodic dynamics (dash-dotted green line).

In general, from Figure 5, we see that the value of $r_{C}$ depends on the whole shape of the inter-event time distribution, displaying, however, some general features. For broad distributions, i.e. small $\alpha$ and $\beta$, the threshold is small and it vanishes for $\alpha\to 1$ and $\beta\to 0$. Then, the threshold increases with $\alpha$ and $\beta$ and it reaches a maximum value for a periodic dynamics i.e. if $\Psi_{a}(\tau)$ tends to a $\delta$-function. From here on, we consider a power-law distribution of inter-event times $\Psi_{a}(\tau)$ as defined in Eq. (2).

4.2 Effects of burstiness on prevalence

Within our analytic approach, beside the epidemic threshold value, a numerical solution of Eq. (8) provides us with the value of $\overline{P}$ in the stationary state. This probability, called prevalence, is the order parameter of the phase transition: it is zero below the epidemic threshold ($r<r_{C}$) and is different from zero above it ($r>r_{C}$). In Figure 6 we plot the prevalence as a function of the control parameter $r=\frac{\langle a\rangle}{\mu}$ for different values of the burstiness exponent $\alpha$, when $\rho(a)\sim\delta(a-a_{m})$ (left panel) and when $\rho(a)\sim a^{-(\nu+1)}$ as in Eq. (1) (right panel).

Plots confirm that burstiness lowers the epidemic threshold as discussed above. However, they also illustrate that the global effect of burstiness ($\alpha\rightarrow 1$) is twofold: it raises the prevalence in low infective systems ($r<1$) strengthening the epidemic, while it weakens the epidemic in high infective systems ($r>1$), lowering the prevalence. This behavior helps clarifying the conflicting effects of burstiness reported in the literature. At $\lambda=1$ (see Fig. 6), the switch from a regime where the burstiness strengthens the epidemics to the regime where epidemics is weakened occurs at $r=1$ and $\overline{P}=1/2$: in B we show that, regardless the shape of $\Psi_{a}(\tau)$ and $\rho(a)$, this property holds when $r=1/\lambda$. One can argue that the twofold effect of burstiness on prevalence is purely due to fluctuations. In this perspective, in C we show that the same behavior can be induced by activity fluctuation in systems without burstiness.

4.3 Critical scaling behavior

The study of so-called spreading experiments – in which a single seed of infection is placed into an otherwise quiescent network – is a standard tool to locate the phase transition between the active and the absorbing state. In such experiments one usually computes, after the initial seed is placed, the probability $P_{d}(t)$ that a run reaches the absorbing state in a time interval $[t,t+dt]$, and the total number of infected nodes $N(t)$ at a given time $t$. These quantities scale as power laws at criticality, with exponents $\delta$ and $\eta$, respectively, defined by $P_{d}(t)\sim t^{-\delta-1}$ and $N(t)\sim t^{\eta}$ [40]. In particular, in the DP universality class the mean-field exponents take the values $\delta=1$ and $\eta=0$ [40]. These are typically observed in high dimensional lattices, dense networks and fully connected systems. In Figure 7, for $\rho(a)\sim\delta(a-a_{m})$, we verify that the exponent $\delta$ in our activity driven model model is consistent with the mean field value $\delta\simeq 1$ both for exponential and power-law inter-event time distributions. A similar result can also be obtained for the exponent $\eta\simeq 0$.

In finite size systems, the probability of reaching the absorbing state is:

[TABLE]

where ${\cal{F}}(tL^{-z})$ is an unspecified scaling function, controlling the finite-size effects, $z$ is the dynamical exponent and $L$ is –in general– a measure of the linear size of the system [40]. The case of a lattice of size $N$ we get $N\sim L^{d}$, where $d$ is the dimension. Thus, in general, the average lifetime is $\tau_{life}=\langle t\rangle\sim\int dt\,t^{-\delta-1}\,t\,{\cal{F}}(tN^{-z/d})\sim{N}^{z(1-\delta)/d}$. Since in the mean-field case $\delta=1$, the integral results into a logarithm and one expects $\tau_{life}$ to grow logarithmically with network size $\tau_{life}\sim\log(N)$. Similarly, as $\tau_{life}$ is peaked around a central value, the position of the peak, $\tau_{life}^{max}$, scales in the same way.

As shown in Figure 8 this is, indeed, the case in numerical simulations for either exponential (Poisson) or power-law distributed inter-event time distributions keeping $\rho(a)\sim\delta(a-a_{m})$ (see figure caption).

In the light of the above, let us adopt a finite-size scaling approach – similar to that proposed in [39] – to analyze lifespans (such as those shown in Figure 1) as a function of $N$ and the distance $\Delta$ to the critical point: $\Delta\equiv(r-r_{c})/r_{c}$. Since right at the critical point ($\Delta=0$) lifespans grows as $\ln(N)$, slightly away from the epidemic threshold we get the finite-size scaling relation:

[TABLE]

where $\gamma$ is expected to take its mean-field value $1/2$ 222Actually, the characteristic length-scale scales in mean field as $L\sim\Delta^{-1/2}$ and $N\sim L^{d}$, where in mean-field $d=d_{C}=4$ is the critical dimension, so that $N\sim\Delta^{-2}$ and, thus $\gamma=1/2$ [40]. Finally, in Figure 9 we verify the scaling hypothesis of Eq. (13) both in the exponential and in the power law case, showing a rather good curve collapse when rescaling (dividing) $\tau_{life}$ with $\ln(N)$ and $\Delta$ with a factor $N^{\gamma=1/2}$.

We also measure the critical exponent $\beta$, that describes the vanishing of the prevalence $\overline{P}$ at the epidemic threshold:

[TABLE]

for $r>r_{C}$. In C we show that $\beta=1$ for the Poisson process, regardless of $\rho(a)$. Figure 10 shows that also in the bursty case $\beta\simeq 1$ for every $\alpha$ and activity distribution $\rho(a)$.

Summarizing, despite the non-trivial temporal dynamics, critical exponents suggest that burstiness does not change the universality class of the absorbing-active phase transition for the SIS model on activity-driven network, when the average inter-event time is finite ($\alpha>1$). The system presents the same mean-field like scaling properties for non-Poissonian (bursty) and Poissonian dynamics, regardless of the exponent $\alpha$, and shows the mean-field exponents of DP universality class. A similar equivalence is obtained in [21], where burstiness is introduced on static networks as a non-Poissonian effect in infection process.

5 Conclusions

Our results provide a coherent understanding of the interplay between the bursty dynamics of an activity-driven network and the epidemic process unfolding on it. Focusing on the SIS model, by means of an activity-based mean-field approach we analytically derive a closed form for the epidemic threshold on an activity-driven network: the threshold depends on the whole functional form of the activity and inter-event time distributions. The mean-field approach turns out to be exact in this system, as connections are continuously reshuffled destroying local correlations. We show that the analytical results are in excellent agreement with extensive numerical simulations.

The model captures two key aspects driving the evolution of social networks: burstiness and heterogeneity in human activities. As is well known, heterogeneity produces a decrease in the epidemic threshold. We show that burstiness also reduces the epidemic threshold, while the effect on prevalence is more subtle: in low-infective systems burstiness raises the prevalence, while it weakens the epidemics in high infective systems. We showed that this behavior can be attributed to fluctuations effects. Finally, despite the non-trivial temporal dynamics, the values of the scaling exponents suggest that burstiness does not change the universality class of the absorbing transition.

Our results can help to clarify the effect of burstiness on dynamical process in temporal network and, in turn, to design efficient immunization strategies [41] and control protocols, based on temporal structure of the human interactions. The model can be extended to different epidemic processes and the analytic results for the thresholds should hold also in the case of the SIR model, while a different behavior is expected in the active phase [37]. We can also include in the model further features of realistic social networks, i.e. the presence of communities [42], adaptive behavior [43, 44] and memory [37], these combined with burstiness can produce non-trivial effects [8].

Appendix A Evaluation of the epidemic threshold

Let us consider the system of equations (7). We set an initial condition, adding $\varepsilon_{a}\delta(t)$ to $Q_{a}(t)$: i.e. we impose that at $t=0$ all nodes are active and a node with activity $a$ is infected with probability $\varepsilon_{a}$. For $\widetilde{Q}_{a}(s)$ and $\widetilde{P}_{a}(s)$, the Laplace transforms of $Q_{a}(t)$ and $P_{a}(t)$, we obtain:

[TABLE]

[TABLE]

We substitute $\widetilde{Q}_{a}(s)$ in $\widetilde{P}_{a}(s)$ and expand both for $s\sim 0$, since for the calculation of the epidemic threshold we only need the asymptotic behavior ($t\rightarrow\infty$). We obtain at first order in $\mathcal{O}(\frac{1}{s})$:

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

Going back to the temporal domain, in the asymptotic limit we obtain:

[TABLE]

where

[TABLE]

In this way we obtain two self-consistency equations for $Q_{a}^{0}$ and $P_{a}^{0}$:

[TABLE]

[TABLE]

since $\overline{P}=\int_{0}^{\infty}da\,\rho(a)\,P_{a}^{0}$ and $\overline{aQ}=\int_{0}^{\infty}da\,a\,\rho(a)\,Q_{a}^{0}$. Integrating the above expression over $a$ we get the self consistence equations described in Eq. (8).

Appendix B Prevalence at $r=1/\lambda$

We can write Eqs. (19) and (20) for $P_{a}^{0}$ and $Q_{a}^{0}$ as:

[TABLE]

[TABLE]

We multiply Eq. (21) by $\rho(a)$ and Eq. (22) by $\rho(a)a/\mu$. Then we integrate over the activity obtaining for $\overline{aQ}$ and $\overline{P}$:

[TABLE]

[TABLE]

where $B=\int_{0}^{\infty}da\,\rho(a)\frac{a(1-L_{a}(\mu+\lambda\overline{aQ}))}{1-(1-\lambda\overline{P})L_{a}(\mu+\lambda\overline{aQ})}$, $H=\int_{0}^{\infty}da\,\rho(a)\frac{aL_{a}(\mu+\lambda\overline{aQ})}{1-(1-\lambda\overline{P})L_{a}(\mu+\lambda\overline{aQ})}$. Setting $\langle a\rangle/\mu=r=1/\lambda$, Eq. (24) gives:

[TABLE]

Plugging Eq. (25) in Eq. (23) and setting $\mu=\lambda\langle a\rangle$, we obtain:

[TABLE]

Since

[TABLE]

then Eq. (26) becomes:

[TABLE]

Since $H>0$, $\overline{aQ}\geq 0$ and $\overline{P}\geq 0$, Eq. (28) has only two solutions: the absorbing state $(\overline{P},\overline{aQ})=(0,0)$ which is unstable because from Eq. (9) we get $r_{C}<1/\lambda=r$; and $\overline{P}=1/2$ which represents an endemic state, stable at $r=1/\lambda$.

In this way we proved that $\overline{P}=1/2$ if $r=1/\lambda$ regardless of the inter-event time distribution $\Psi_{a}(\tau)$ and of the activity distribution $\rho(a)$. Therefore, as observed in Fig. 6, the prevalences $\overline{P}$ calculated for different $\alpha$ cross at $r=1/\lambda$ and so there is a switch at $r=1/\lambda$ in the effect of burstiness on epidemic spreading.

Appendix C Exponential inter-event time distributions

In this section we discuss two properties of the model with exponential inter-event time distribution. First we calculate exactly the critical exponent $\beta=1$. In this case the probability $P_{a}(t)$ is governed by the standard activity based mean-field differential equation:

[TABLE]

Imposing the steady state condition $\partial_{t}P_{a}(t)=0$ we get the stationary value $P_{a}^{0}$ of the $P_{a}(t)$:

[TABLE]

Averaging Eq. (30) over $a$ we get $\overline{aP}$ as a function of $\overline{P}$:

[TABLE]

Plugging Eq. (31) in Eq. (30) we obtain for $P_{a}^{0}$:

[TABLE]

Multiplying both members for $\rho(a)$ and integrating over the activity we obtain a self-consistency equation for $\overline{P}$:

[TABLE]

Eq. (33) admits the solution $\overline{P}=0$ representing the absorbing state. Other possible solutions can be obtained by solving for $\overline{P}$ the equation:

[TABLE]

To evaluate $\beta$ we assume $r\sim r_{C}$ and we can expand Eq. (34) near the critical point, for $\overline{P}\sim 0$. We get:

[TABLE]

Integrating over the activity and keeping only the leading orders in $\overline{P}$ we get:

[TABLE]

where

[TABLE]

and

[TABLE]

and for the exponential case

[TABLE]

Since $\langle a^{2}\rangle>\langle a\rangle^{2}$ we have that $D(r)>0$. Moreover close to the critical point:

[TABLE]

Hence for $r\sim r_{C}$ Eq. (36) has two solutions: $\overline{P}=0$ i.e. the absorbing state, and for $r>r_{C}$

[TABLE]

which means that $\beta=1$.

We now consider the behavior of the prevalence $\overline{P}(r)$ in the active phase as a function of the activity fluctuations, i.e. we vary the exponent $\nu$ characterizing the distribution $\rho(a)$ defined by Eq. (1). We notice that, according to Eq. (39), the threshold $r_{C}$ increases by increasing $\nu$, hence at small $\overline{P}$ close to the epidemic thresholds the prevalence grows with $\nu$. Then the results of B show that at $r=1/\lambda>r_{C}$, for any values of $\nu$, $\overline{P}=1/2$, i.e. the prevalence is independent of activity fluctuations. Finally, for $r>1/\lambda$ the opposite behavior is expected i.e. the system presents a larger prevalence at large $\nu$. The overall behavior is illustrated in Figure 11: large activity fluctuations and large burstiness give rise to a similar effect on the prevalence i.e. an increase at small $r$ and a decrease at large $r$.

Acknowledgements

R.B. acknowledges funding from the INFN BIOPHYS project. M.A.M. acknowledges the Spanish Ministry of Science as well as the Agencia Española de Investigación (AEI) for financial support under grant FIS2017-84256-P (FEDER funds). We warmly thank V. Buendia for very useful comments.