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:


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.

  • $\begingroup$ 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 $\endgroup$ – Ben Bolker Mar 31 '18 at 20:05
  • $\begingroup$ @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. $\endgroup$ – 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!

  • $\begingroup$ How will the threshold change on a directed graph please? $\endgroup$ – Terry BRETT Jun 22 '18 at 14:56
  • $\begingroup$ How directed? In the sense that infection can occur only if the active vertex is infected ? $\endgroup$ – Uollano Jun 23 '18 at 18:15
  • $\begingroup$ 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 $\endgroup$ – Terry BRETT Jun 24 '18 at 14:34
  • $\begingroup$ If I understood well you question, you should get R_0=2 . Can you run simulations of such a system? $\endgroup$ – Uollano Jul 13 '18 at 13:17
  • $\begingroup$ I get R_0 = 2 as well, does that value change for you if you have different values of m? $\endgroup$ – 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)


  • 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) {
pvar <- function(gamma=2.8,epsilon=0.001) {
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):


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

  • $\begingroup$ $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. $\endgroup$ – Terry BRETT Apr 2 '18 at 15:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.