Online human aggregation under pressure moves beyond preferential attachment
Abstract
There is a significant amount of online human activity which is either clandestine or illicit in nature, and hence where individuals operate under fear of exposure or capture. Yet there is little theoretical understanding of what models best describe the resulting dynamics. Here we address this gap, by analyzing the evolutionary dynamics of the supporters behind the 95 proISIS online communities (i.e. selforganized social media groups) that appeared recently on a global social media site. We show that although they do not follow a conventional (i.e. sizebased) preferential attachment (PA) model, their dynamical evolution can be explained by a new variant that we introduce here, which we refer to as active attraction model (AA). This AA model takes into account the locality and group heterogeneity which undoubtedly feature in humans’ online behavior under pressure, but which are not contained in conventional PA models. The AA model captures both groupspecific and macroscopic observations over all size ranges – as opposed to just the tail for large groups or groups’ initial growth – suggesting that heterogeneity and locality play a crucial role in the dynamics of online extremist support. We derive approximate expressions for the group size distributions in two simple systems that involve simultaneously the mechanisms of group joining (governed by either PA or AA), group leaving, and account banning, and show how these processes influence the group size distributions. We believe this work will serve in helping understand a broad spectrum of online human activities which are either clandestine or illicit in nature, and hence where individuals operate under fear of exposure or capture.
I Introduction
Just as the Internet can be used for good, it also serves as an ideal vehicle for the more clandestine and illicit side of human activity. For example, extremist entities such as ISIS (known as Islamic State) make ample use of the Internet for spreading their message and propaganda materials, recruiting young people, and soliciting funds. One particular social media platform VKontakte (www.vk.com), which has million global users and is almost identical to Facebook, became the primary online social media source for ISIS propaganda and recruiting Johnson et al. (2016). Unlike on Facebook where proISIS activity is almost immediately eliminated, support on VKontakte develops around online groups (i.e. selforganized communities) which are akin to Facebook groups that support particular everyday topics such as a sports team. These online proISIS groups may not only organize premeditated attacks, but also incite decentralized lonewolf attacks Johnson et al. (2016). Hence, there is a compelling need to investigate the dynamics of such online groups, especially the earlygrowing ones whose group sizes are relatively small.
Here we analyze the evolutionary dynamics of the supporters behind the 95 proISIS online communities (i.e. selforganized social media groups) that appeared recently on a global social media site. Though focussed on ISIS support, our model, analysis and findings should serve in helping understand a broad spectrum of online human activities which are either clandestine or illicit in nature, and hence where individuals operate under fear of exposure or capture. Indeed, outside of extremism and online illicit or clandestine activity there is already increasing interest in understanding how communities, users, or groups attract new followers/members, and develop over time Backstrom et al. (2006); Mislove et al. (2007); Palla et al. (2007); Leskovec et al. (2009, 2008a); Traud et al. (2012); Backstrom et al. (2008); Berger and Perez (2016); Kairam et al. (2012). Many previous studies highlight the role played by a “nonlocal” preferential attachment (PA) effect in the formation of groups or network clusters Albert and Barabási (2002); Barabâsi et al. (2002); Backstrom et al. (2006); Mislove et al. (2007); Leskovec et al. (2009); Kumar et al. (2010). Unfortunately, we find that conventional PA models cannot well explain the unusual rapid and heterogeneous growths at the early growing stage of the online proISIS groups observed from empirical data (see Sec. II), possibly due to one or more of the following reasons:

The online proISIS supporters and the selforganized groups that they form online, are under pressure. The members in such a group are discussing an extreme topic and supporting an illegal terrorist organization. These individuals have to cohabit the same online space as opposing individuals and entities (e.g. the online organization called Anonymous) as well as government agencies, all of whom are not only trying to defuse the narrative of the extremist social media groups but are also possibly trying to track down the identities of particular group members. As a result, the extremist supporters are under continual pressure and likely act online in ways that help them maintain a more hidden profile. This challenges significantly the typical PA model assumption of implicit knowledge of all group sizes across the whole population.

Extremism discussed in online proISIS groups is a specialist niche topic. It seems less likely that people would be drawn to it simply because it is popular among others. Therefore, the probability of attaching should not only depend on node connectivity, but also should incorporate individual heterogeneity and locality – just as readers of an article in a physics journal are likely to be physicists, and within that subpopulation, readers of an article on network science are more likely to be from that network science community. These are often overlooked factors that can play an important role in group formation and dynamics Leskovec et al. (2008b).

The evolution of online groups is affected by moderators who have the right to ban a group. The conventional PA models fail to consider the consequence of such moderators. Their special role and powers suggest that their presence and activity can influence heavily the group evolution Kelley et al. (2011). While a full theory that includes them would require multispecies analysis that we cannot yet provide, it suffices to say that their presence in the system is another reason that can move the dynamics beyond PA.
These issues need to be addressed in order to get a more complete picture of human activity online. Here we go partway toward trying to address these issues. In so doing, our work also goes partway toward addressing the following wishlist in the network literature: (1) Many previous studies are based on particular definitions and techniques for detecting the groups. It would be very useful to have a more general framework for discussing what a group is Leskovec et al. (2009, 2008a); Yang and Leskovec (2015); Palla et al. (2007). (2) Although many previous studies focus either on the observational aspects Backstrom et al. (2006); Mislove et al. (2007); Backstrom et al. (2008); Berger and Perez (2016); Kairam et al. (2012) or theoretical modeling Sorensen et al. (1987); Gueron and Levin (1995); Leskovec et al. (2008a), knowledge about how network theories agree with the observation on both the microscopic and the macroscopic scales is rare. For instance, a large portion of studies focus on the global statistical properties such as the scaling behavior of the group/cluster size distribution and the evolution of the globally averaged quantities Albert and Barabási (2002); Backstrom et al. (2006); Mislove et al. (2007); Leskovec et al. (2009), yet little is known about the evolution at the groupspecific level. This includes the study in Ref. Johnson et al. (2016) which focuses on the mechanisms producing the tail in the size distribution at larger groups sizes. More detailed groupspecific studies for any group size are needed since global statistics could be misleading, e.g. due to the temporal variation of the global populationBarabâsi et al. (2002).(3) Although there are studies on how microscopic node behavior would reproduce the observed macroscopic statistical properties Leskovec et al. (2008b), knowledge of how individual behaviors contribute to the evolutionary property of a single group is missing.
Specifically, this paper proposes a simple growth model, namely the active attraction (AA) model, that goes beyond PA and takes into account the locality and heterogeneity of online activity. We show that this active attraction model (AA) captures both groupspecific and macroscopic observations over all size ranges – as opposed to just the tail for large groups Johnson et al. (2016) or groups’ initial growth and development. Our findings suggest that heterogeneity and locality play a crucial role in the dynamics of online extremist support. We also derive approximate expressions for the group size distributions in two simple systems that involve simultaneously the mechanisms of group joining (governed by either PA or AA), group leaving, and account banning, and show how these processes influence the group size distributions.
The outline of the paper is as follows. In Sec. II we introduce the dataset and the findings. In Sec. IV we model the system by a PA mechanism, and check if the groupspecific and macroscopic observations can be reproduced by the model. In Sec. IV we introduce the AA model, and show how the groupspecific and macroscopic observations are reproduced by the model. In Sec. V, we derive the analytic expressions for the group size distributions of two simple systems involving simultaneously the group joining (governed by either AA or PA), leaving, and banning processes. The conclusion is given in Sec. VI.
Ii Empirical data and analysis
Our dataset is assembled using the same methodology as Ref. Johnson et al. (2016). In contrast to Ref. Johnson et al. (2016) that focuses on the late stage when groups have become very large, we focus here on the early growing stage when the group sizes are relatively small. The dataset used in this work comprises the 95 groups identified as being proISIS Johnson et al. (2016) whose dates of first appearance are within our observational period of 320 days. These provide us with detailed information about the evolution of the group memberships with a high temporal resolution up to one day. Specifically, there are three main processes involved in the group evolution: the group joining and leaving events that may occur every day during a group’s lifetime, and the banning of a group by the moderators. For each group, the dataset provides information about the size, the number of joining and leaving events on each day, as well as the first appearance date which we take as being the first day on which the group has at least have one member, and the banning date if the group gets banned within the period of observation. This banning can be identified by an abrupt drop in the group size to zero.
We started by examining the evolution of the group sizes. According to PA, the size of an earlygrowing group should be small. However, we found that the group sizes can become very large within a few tens of days, causing the sharkfin shaped growth. This sharkfin shape is defined by the distinctive concave shape during growth (Fig. 1(a)). We then studied the correlation between the group joining and leaving rates on day and the group size on that day (Fig. 1(b)). For a given group, the joining/leaving rate on day is estimated by applying a linear fitting to the cumulative number of joining/leaving events during the 5 days around day . We carried this out for all groups and all days – except for the first and last 2 days, since we need 5 days around day to do the fitting. We checked that our findings are insensitive to the precise value of this time window. We found that the groupjoining data points are highly dispersed when the group size is small (we consider a group is small if its size is less than following Ref. Johnson et al. (2016)). This suggests that nonPA rules may apply during the early growing stage. Groupleaving data points can be well fit by a line whose slope is found to be from linear regression (Fig. 1(b)), indicating that leaving a group is more like an independent personal decision than the act of joining a group.
However, this is still insufficient to exclude PA as the governing mechanism for the group joining process during the early growth stage, for the following two reasons: First, a rapid growth of the (global) total number of group members or followers may still result in decelerating growth during the early growth stage Barabâsi et al. (2002). Second, it could be that the temporal fluctuations in the total number of followers caused the dispersion of data points observed in Fig. 1(b). Therefore in the next section we attempt to test out how well PA reproduces the observations, by assuming that PA governs the groupjoining process.
Iii Preferential Attachment (PA) model
iii.1 Size evolution of a single group
We first model the size evolution of a single group. According to PA, the increase of the group size on day (before the group gets banned) is given by
(1) 
where and describe the group joining and leaving rates, respectively, and is the total number of followers/members in the system on day (shown in Fig. 2(a)(b)). We define as the group’s appearance date, i.e., the first day that the group’s size is nonzero. Given and initial size , is estimated to be approximately (Fig. 1(b)). Since can be directly obtained from the data, we can iteratively estimate its size on all future days. Hence the curve fitting problem is to find the optimal and that minimize the Pearson’s :
(2) 
where is the day when the group gets banned, is the group size on day estimated using the iterative expression (i.e. Eq. 1), and is its corresponding observed one. The minimization can be carried out by conventional multivariable optimization algorithms. Note that Eq. 1 is only valid during the early growth stage when the saturation effect (i.e. the effect of finite population) can be ignored. Hence to do the fitting for each group in the dataset, we only used the sizes of the first 20 days when the group size is nonzero.
We find that the best estimates of are spread across a broad range, but the fittings for most of the groups are still poor (e.g. see Fig. 2(a) where the fittings barely reproduce the concave shapes of the empirical profiles). This indicates that the empirical growth rate of is far from sufficient to bring about the observed sharkfin shapes. Also against the PA is the fact that, rigorously speaking, conventional PA models implicitly assume a constant for all groups.
iii.2 Stochastic simulation of PA
We simulated the group growth using PA as follows. We set all the parameters (including the total number of groups, the creation and banning date, the total number of new joining events on each day, etc.) to be the same as the data, except that

we redistribute the new group joining events observed on each day to all the nonbanned groups following the PA rule (i.e. the probability that a user joins a group on day is proportional to the size of the group on day , );

we use the constant group leaving rate () estimated from Fig. 1(b);

for a group that appears for the first time on day , we manually assign a small initial size to it, e.g. 1. We also tested other larger values, but the main results are the same.
The exact steps in the simulation are as follows. On day , we first detect from the dataset which groups exist (i.e. have at least one member) and denote them by a set, . We also obtain directly from the dataset, the total number of new joining events (). Next we redistribute the new joining events to the alive groups by drawing a sample from the multinomial distribution, , where is a vector of probabilities whose values sum up to 1 and are proportional to the sizes of the groups on the previous day, By so doing, we ensure that the number of new joining events for a group is always proportional to its current size, and the total number of new joining events is the same as the data. Finally for each group, we consider the group leaving events by subtracting a value from its previous size. The value is drawn from the Binomial distribution , where is the size of the group on the previous day.
Based on the simulation results, we are also able to compare the probability density function (PDF) of the group sizes in the simulation to that in the empirical observation. The PDFs for both cases are obtained as follows. First, we record the sizes of all the groups on each day from the 80100th day to a list . We choose these days when doing the statistics because they correspond to the mature stage of the whole system, and hence have the maximum total number of groups. We also tested other days around and reduced the number of days, but the results are similar. Then we use the Fit module in the wellknown Python package powerlawAlstott et al. (2014) to obtain the PDF.
We find from this PA simulation that (1) the growths of most of the groups that are created at a later point in time, are effectively suppressed by a couple of groups that were created earlier (Fig. 1(c)). (2) The groupjoining data points in Fig. 1(d) are much less dispersed than those in the empirical observation (Fig. 1(b)). (3) The PAbased size distribution doesn’t agree well with the data. In short, neither the groupspecific nor the macroscopic observations can be well reproduced by a PA model.
Iv Active Attraction (AA) model
iv.1 Size evolution of a single group
Here we introduce the AA model. The AA starts by making two reasonable assumptions:

Public group events (e.g. a rapid sharing of an interesting post, or an invitation letter simultaneously sent to a large population, etc.) can attract users to join quickly even when the group size is still small.

Groups are heterogeneous in that a given group may not be easily accessible to all users, or it may not be of equal interest to all users.
The first assumption can introduce some nonPA effects, while the second stresses the locality and heterogeneity in usergroup interactions. With these two assumptions, one can imagine observing a stairlike growth in group size (e.g. Fig. 2(c)) if several major public group events occur in turn. This motivates our new AA model beyond PA, and its name which reflects a more specific usergroup active attraction. In particular, based on the fact that the groups in our dataset are under pressure and the topic itself will only appeal to an extreme fraction of the population of online users, a group may be accessible or of interest to only a portion of the whole population. This portion is determined by a saturation level denoted by . For simplicity, we assume changes linearly with time, i.e. , where is the initial size, and is the rate of change. Next, we define a joining rate ( constant in time, but could be different for each group), a leaving rate (, as mentioned earlier and found to be approximately for all groups), and an average (over all groups) groupbanning rate , which is only used in the simulation. With these settings, the evolution of the size of a single group before banning (hence no dependence) in the continuous limit is given by
(3) 
and . We further define as the effective joining rate, and as the effective saturation level. Let be , and be . Then the differential equation can be rewritten as , whose solution is given by
(4) 
which contains 4 free parameters (i.e., , , and ), and is the expression used for the fittings by the AA model shown in Fig. 2(a)(b). The fitting is done using a conventional multivariable optimization algorithm. Specifically, we used the curve_fit function with bounds set properly in the wellknown Python package SciPyJones et al. (01).
We find that Eq. 4 fits most of the sharkfin group growth profiles well, not only in terms of their early growth stage (e.g. Fig. 2(a)), but also over their entire lifetime in many cases (e.g. Fig. 2(b)). In addition, we find that there exists a high heterogeneity in the fitting parameters. That is, , and ( can be either positive or negative) all spread over a broad range – and indeed, they scale roughly like a powerlaw with the powerlaw exponent , except which is not very sensitive and is always just a couple of days earlier than the observed firstappearing day. We also tried to fix the three parameters (i.e., , and ) to a reasonable constant value (e.g., the median or mean value from previous fittings that allowed them to change), but found it resulted in poorer fittings. Hence, for the empirical data, the heterogeneity indeed exists in all these three parameters. Since the distributions of these three parameters are not the focus of this work, we show only the distribution of in Fig. 4, which is used in the simulations.
iv.2 Stochastic simulation of AA
We now show how we performed the simulation of the AA model. Similar to the simulation of the PA model, we controlled all the parameters to be the same as the dataset, except that (1) we adopt a constant groupleaving rate of ; (2) we redistribute the new joining events on each day among all the alive groups on that day according to the AA rule (i.e. the probability that a user joins a group on day is proportional to , where and are the saturation level and the group size on day , respectively); (3) we further assume for each group (and hence is constant), and set the saturation level of each group to () estimated from the curve fitting by Eq. 4. Note that the AA rule we used implicitly assumes is a constant for all the groups – otherwise the probability of joining a group should be proportional to – and hence the simulation focuses on studying the effect of the heterogeneity in the saturation level. We could have also set the of each group to be the value obtained from the fitting, which would make our simulation agree even better with the dataset, but our simulations show that the heterogeneity in is sufficient to reproduce the empirical observations. We did the simulation in a way that is similar to the PA simulations, as follows: On day , we first detect from the dataset which groups are alive (i.e. have at least one member) and denote them by a set, , and also obtain directly from the dataset the total number of new joining () events. Next we redistribute the new joining events to the alive groups by drawing a sample from the multinomial distribution, , where is a vector of probabilities whose values sum up to 1 and are proportional to . Finally for each group we consider the group leaving events by subtracting a value from its previous size; the value is drew from the Binomial distribution .
V Theoretical analysis of the group size distribution
We now study mathematically the group size distribution of such a system involving simultaneously the three mechanisms of group joining, leaving, and banning due to moderators. Since similar systems can exist in many other areas of human extremist activity, we cast our discussion in a more general sense and show the analytical results for both the case when the group joining follows the AA rule, and that when it follows the PA rule. Conventional rigorous treatment of the problem may involve solving the master equations, but become almost intractable – especially when taking into account the heterogeneity of the coefficients. Hence, we resort to an approximative method and study only two simple cases. Here are the details of the derivations.
v.1 The AA case
We start with the case of the group joining being an AA process. For simplicity, we consider the case that the groups are created with the same probability on each day during the periods of observation; hence the probability that a group created on day will remain alive (or observable) on day is , where is the probability that a group will get banned on a day, and . As mentioned earlier, in the case that the joining is an AA process then the size evolution of a group, (hereafter, we omit its label for convenience), is governed by . For simplicity, we consider the case when and are both timeindependent. In such a case, the solution of this differential equation when the initial condition is given by , from which we can inversely get , and hence . Consider the distribution of in Fig. 4, the detailed shape of which ( though not the main concern of this work) is probably irregular due to poor statistics or the irregular banning events of the moderators etc., but which very roughly distributes around a powerlaw whose exponent is . Inspired by this shape, we consider the simple case that the effective saturation level of all the groups follows a powerlaw distribution of , and that is the same for all groups. Then by ignoring the stochastic fluctuations in the group size evolution profiles (meaning that there is now a onetoone correspondence between and for a given ), the probability density that a group created at is observed at having size is given by . Expanding the integrand with respect to and keeping the terms up to , we obtain , i.e. it also follows a powerlaw distribution that scales approximately the same way with the distribution of . This similarity between the saturation level distribution and the group size distribution is also found in the empirical observation (Fig. 4).
v.2 The PA case
For comparison, we also calculate in a similar way the group size distribution for the case that the group joining is a PA process. Neglecting the finitepopulation effect and considering only the case when is a constant, we have , whose solution . Hence and . Therefore, the probability density of the group size is given by , where , which is greater than 1 for any and .
v.3 Numerical verification and discussion
We now check these analytical results since their derivations involve several approximations. We conducted stochastic simulations for both systems. The simulation settings comply to the basic assumptions we made when defining the systems for the analytic derivations. For the AA case, we initialize 1000 groups by assigning to each of them an initial group size of 1 and a saturation level () sampled from a Zipf’s distribution with a powerlaw exponent equaling . Next, on each day we generate the number of joining events for each group by sampling from the binomial distribution and add it to to obtain , where is the joining rate. Then we generate the number of leaving events for each group by sampling from , where is the leaving rate, and without loss of generality was set to be . Finally we ban each group with a probability (without loss of generality, we use a value of 0.02) by setting the number of members to 1, which means after banning, we immediately create another group so that the total number of groups is always a constant.
We run for a sufficiently large number (e.g. 2000) of steps (i.e. days) for each run. This means that the creation dates of the groups are effectively randomized automatically after a long period due to the randomness in the group banning, and the total number of followers will also converge to a stable level. In addition, to improve the precision, we repeat the simulation 10 times for every set of parameters and then obtain the average distributions by combining the sample from each run. For each run, without loss of generality, the sample is formed by the group sizes on each day during the 80100th day of the last 320 days in the simulation.
With respect to the simulation for the PA case, everything is the same as for the AA case, except that we don’t need anymore; on each day , the first step in the AA case (see above) is to generate the number of joining events for each group; the number is sampled from the binomial distribution and is added to to obtain .
Figure 6 shows that our analytical results agree well with the stochastic simulations. We can see that for the AA model, the distribution of the saturation level plays a pivotal role in determining the group size distribution. By contrast, for the PA model all three mechanisms matter. In addition for the PA case, when the group banning rate is small and the group leaving rate is significantly smaller than the joining rate, the powerlaw exponent is always around 1. We can also see that different microscopic usergroup interaction mechanisms (e.g. AA and PA) may result in the same powerlaw exponent in the macroscopic statistics, which means that how a macroscopic quantity scales could be insufficient to tell the microscopic mechanism. Hence user and group level analysis – as presented in this paper – becomes important.
Vi Conclusion
In conclusion, we have shown that the nonPA effect is crucial for explaining the rapid growth of groups observed during the early growth stage. The PA effect then becomes more appropriate when the group size is large, as indicated in Ref. Johnson et al. (2016)). We proposed a simple nonPA model that catches the effect of locality and microscopic heterogeneity in the dynamics of group formation, and which are overlooked by conventional PA models. While we concede that there could be alternative models, our model is arguably one of the simplest, and it explains well both the groupspecific and global statistical observations. Since the AA process is sizeindependent, it provides a novel and reasonable explanation for the cascading joining that results in a sharkfin shape observed in group size evolution. Such an observation can barely be described by a PA model without introducing an unusually high growth rate for the global total number of followers. In addition, since both PA and AA could produce a similar group size distribution, this work suggests the importance of a deeper understanding of the behaviors of individual users and groups.
There are still many open questions. For instance, the origin of the user and group level heterogeneity is still a mystery. More specifically, it is not clear if the broad distribution of the saturation levels originates from the topology of the network that influences the accessibility of the groups (i.e. the locality), or is due to some other more complicated mechanisms that result in the heterogeneity in users’ interests in the groups. Hopefully the present work draws researchers’ attention to these open questions, and serve as a steppingstone toward answering them.
Acknowledgments We are grateful to Yulia Vorobyeva, Andrew Gabriel and Anastasia Kuz for initial help with data collection and analysis. NFJ gratefully acknowledges funding under National Science Foundation (NSF) grant CNS1522693 and Air Force (AFOSR) grant FA95501610247. The views and conclusions contained herein are solely those of the authors and do not represent official policies or endorsements by any of the entities named in this paper.
References
 Johnson et al. (2016) N. Johnson, M. Zheng, Y. Vorobyeva, A. Gabriel, H. Qi, N. Velasquez, P. Manrique, D. Johnson, E. Restrepo, C. Song, et al., Science 352, 1459 (2016).
 Backstrom et al. (2006) L. Backstrom, D. Huttenlocher, J. Kleinberg, and X. Lan, in Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining (ACM, 2006) pp. 44–54.
 Mislove et al. (2007) A. Mislove, M. Marcon, K. P. Gummadi, P. Druschel, and B. Bhattacharjee, in Proceedings of the 7th ACM SIGCOMM conference on Internet measurement (ACM, 2007) pp. 29–42.
 Palla et al. (2007) G. Palla, A.L. Barabasi, and T. Vicsek, Nature 446, 664 (2007).
 Leskovec et al. (2009) J. Leskovec, K. J. Lang, A. Dasgupta, and M. W. Mahoney, Internet Mathematics 6, 29 (2009).
 Leskovec et al. (2008a) J. Leskovec, K. J. Lang, A. Dasgupta, and M. W. Mahoney, in Proceedings of the 17th international conference on World Wide Web (ACM, 2008) pp. 695–704.
 Traud et al. (2012) A. L. Traud, P. J. Mucha, and M. A. Porter, Physica A: Statistical Mechanics and its Applications 391, 4165 (2012).
 Backstrom et al. (2008) L. Backstrom, R. Kumar, C. Marlow, J. Novak, and A. Tomkins, in Proceedings of the 2008 International Conference on Web Search and Data Mining (ACM, 2008) pp. 117–128.
 Berger and Perez (2016) J. Berger and H. Perez, (2016), available at https://pdfs.semanticscholar.org/2eb1/e080058913c90a0ec1aa2ecabd1db9589c8f.pdf.
 Kairam et al. (2012) S. R. Kairam, D. J. Wang, and J. Leskovec, in Proceedings of the fifth ACM international conference on Web search and data mining (ACM, 2012) pp. 673–682.
 Albert and Barabási (2002) R. Albert and A.L. Barabási, Reviews of modern physics 74, 47 (2002).
 Barabâsi et al. (2002) A.L. Barabâsi, H. Jeong, Z. Néda, E. Ravasz, A. Schubert, and T. Vicsek, Physica A: Statistical mechanics and its applications 311, 590 (2002).
 Kumar et al. (2010) R. Kumar, J. Novak, and A. Tomkins, in Link mining: models, algorithms, and applications (Springer, 2010) pp. 337–357.
 Leskovec et al. (2008b) J. Leskovec, L. Backstrom, R. Kumar, and A. Tomkins, in Proceedings of the 14th ACM SIGKDD international conference on Knowledge discovery and data mining (ACM, 2008) pp. 462–470.
 Kelley et al. (2011) J. L. Kelley, L. J. Morrell, C. Inskip, J. Krause, and D. P. Croft, PloS one 6, e24280 (2011).
 Yang and Leskovec (2015) J. Yang and J. Leskovec, Knowledge and Information Systems 42, 181 (2015).
 Sorensen et al. (1987) C. Sorensen, H. Zhang, and T. Taylor, Physical review letters 59, 363 (1987).
 Gueron and Levin (1995) S. Gueron and S. A. Levin, Mathematical biosciences 128, 243 (1995).
 Alstott et al. (2014) J. Alstott, E. Bullmore, and D. Plenz, PloS one 9, e85777 (2014).
 Jones et al. (01 ) E. Jones, T. Oliphant, P. Peterson, et al., “SciPy: Open source scientific tools for Python,” (2001–), [Online; accessed on 20170925].