# Activity driven model threshold

So I am trying to figure out a threshold for an activity driven model. That is described in this paper: https://www.nature.com/articles/srep00469

That is I don't see how the researchers got R0 to be between 0.2 and 0.5 (page 6). If they are using gamma=-2.1 and epsilon=10^-3 and threshold is defined in equation 4, my calculations give me threshold of 2.020396043515874.

Now if the power law distribution (activity) is a random number i.e. this: http://mathworld.wolfram.com/RandomNumber.html

this yields average activity of 0.005489691793543477 using this equation:

((1-gamma)/(2-gamma))*((1-(epsilon**(2-gamma)))/(1-(epsilon**(1-gamma))))


so then following equation 4, the threshold is 2.020396043515874, but that's not what Figure 4 (B) is showing so I'm very confused.

I would really appreciate if someone could clear this out for me.

• can you fill in a few more steps of your calculation? Where does your average activity formula come from? I'm getting 0.59 following the mean and variance formulas from en.wikipedia.org/wiki/Pareto_distribution - but in any case your question will be much more useful with more details – Ben Bolker Mar 31 '18 at 20:05
• @BenBolker Average activity formula comes from $\int_{\epsilon}^{1}af(x)da$. At each time $t$ a node will be activate with probability $a$, which is between $\epsilon = 10^{-3}$ and 1 and each active nodes will get m links. We then infect 1% of node, i.e. 1% of 10,000 and run SIR model on that. – Terry BRETT Apr 2 '18 at 16:04

I think you just get some calculations wrong, maybe for the second moment $<a^2>$ of the activity distribution. All the calculations are just simple algebra and definite integrals of powers. For a power-law distribution $F(a)\propto a^{-\gamma}$ with $a\in[\epsilon,1]$, in order to have $\int^1_\epsilon F(a)\:da=1$ i.e. the correct normalization, we have: $$F(a)=\frac{\gamma-1}{\epsilon^{1-\gamma}-1}a^{-\gamma}$$ The moments of the distributions $<a^n>=\int^1_\epsilon a^nF(a)\:da$ are: $$<a^n>=\frac{\gamma-1}{\gamma-1-n}\frac{\epsilon^{1+n-\gamma}-1}{\epsilon^{1-\gamma}-1} \quad if \quad n\ne\gamma-1$$ The case $n=\gamma-1$ have a logarithm and it is not relevant in our case. For $\gamma=2.1$, $\epsilon=.001$, $n=1$ and $n=2$ we get: $$<a>= \frac{\gamma-1}{\gamma-2}\frac{\epsilon^{2-\gamma}-1}{\epsilon^{1-\gamma}-1}= 0.00548969179354348\\ <a^2>=\frac{\gamma-1}{\gamma-3}\frac{\epsilon^{3-\gamma}-1}{\epsilon^{1-\gamma}-1}=0.00061164650162922$$ With these numbers I get: $$R_0=\frac{2<a>}{<a>+\sqrt{<a^2>}}=0.36330095782279$$ which is between $0.2$ and $0.5$. I'm not the best at explaining things, so let me know if there are any doubts! Goodluck!

• How will the threshold change on a directed graph please? – Terry BRETT Jun 22 '18 at 14:56
• How directed? In the sense that infection can occur only if the active vertex is infected ? – Uollano Jun 23 '18 at 18:15
• Directed in a way that an infected node will be able to infect another, but not vice versa.So if node A is infected, and it has a direction to node B, A->B, the infection can travel from A to B, but not back (if A would recover). At the moment the graph wasn't directed, which made it bidirectional really A<->B, but now I want to introduce that direction in which infection can travel, which would change the threshold – Terry BRETT Jun 24 '18 at 14:34
• If I understood well you question, you should get R_0=2 . Can you run simulations of such a system? – Uollano Jul 13 '18 at 13:17
• I get R_0 = 2 as well, does that value change for you if you have different values of m? – Terry BRETT Aug 8 '18 at 15:39

(A little too long for a comment, and maybe useful)

Assuming that a is supposed to represent the degree of a node, and using the formulas from Wikipedia for the Pareto distribution, I have

• mean (which I'm assuming is <a>) = gamma*epsilon/(gamma-1)

and

• variance (equal to <a^2>-<a>^2, or <a^2>= mean^2 + var) = epsilon^2*gamma/((gamma-1)^2*(gamma-2))

Translating into R code (rather than doing the algebra):

pmean <- function(gamma=2.8,epsilon=0.001) {
gamma*epsilon/(gamma-1)
}
pvar <- function(gamma=2.8,epsilon=0.001) {
epsilon^2*gamma/((gamma-1)^2*(gamma-2))
}
prms <- function(gamma=2.8,epsilon=0.001) {
sqrt(pvar(gamma=gamma,epsilon=epsilon)+ pmean(gamma=gamma,epsilon=epsilon)^2)
}


When I use the value of gamma=2.1 in the figure caption and feed it into equation (4):

2*pmean(gamma=2.1)/(pmean(gamma=2.1)+prms(gamma=2.1))


I get 0.588 (still not right, but closer ... ??)

• $a$ is the activity of a node at time $t$ if a node is active, it will connect with $m$ other nodes (in this case 2). The activity is a power law from the link I've got above. – Terry BRETT Apr 2 '18 at 15:57