Structural transition in social networks: The role of homophily
Abstract
We introduce a model for the formation of social networks, which takes into account the homophily or the tendency of individuals to associate and bond with similar others, and the mechanisms of global and local attachment as well as tie reinforcement due to social interactions between people. We generalize the weighted social network model such that the nodes or individuals have features and each feature can have different values. Here the tendency for the tie formation between two individuals due to the overlap in their features represents homophily. We find a phase transition as a function of or , resulting in a phase diagram. For fixed and as a function of the system shows two phases separated at . For large, homogeneous, and well separated communities can be identified within which the features match almost perfectly (segregated phase). When becomes larger than , the nodes start to belong to several communities and within a community the features match only partially (overlapping phase). Several quantities reflect this transition, including the average degree, clustering coefficient, feature overlap, and the number of communities per node. We also make an attempt to interpret these results in terms of observations on social behavior of humans.
Introduction
In human societies homophily, the tendency of similar individuals getting associated and bonded with each other, is known to be a prime tie formation factor between a pair of individuals [1]. This association and bonding can be related to one or more features including gender, race, age, education level, economic and social status, and many more. Consequently homophily has a large impact on a number of fundamentally important social phenomena like segregation, inequality, perception biases, and the transmission of information between groups of individuals [2, 3, 4, 5]. The mechanism of focal closure, i.e., the process of forming links to others with shared characteristics but without common acquaintances can be considered as a typical manifestation of homophily [6, 7].
The system that emerges due to interpersonal tie formation, constitutes a social network, which is one of the paradigmatic examples of complex networks [8, 9]. Noticeable features of social networks are the significantly high occurrence of triangles (clustering) and the structure of communities [10], which are the densely “wired” parts of the network as compared to their connectivity to the rest of the network. Communities themselves have a complex structure, characterized by overlaps, hierarchy, and multiplexity [11, 12, 13, 14, 15, 16, 17]. Moreover, the topological features are correlated with the intensity of the ties represented by the weights of the links. The “strength of weak ties” hypothesis proposed by Granovetter [18], stating that a stronger tie between two persons leads to larger overlap between their neighbors, has been validated on empirical data [19, 20]. It suggests that the communities consist of strong links while intercommunity links are weak. Studies have shown that microscopically a key inducing factor of such community structure is the mechanisms of cyclic closure [6, 7], the tendency of nodes in social networks to make links with a topologically close nodes, i.e., people often form social ties with a person sharing common friends.
The two main mechanisms for social tie formation, the focal and cyclic closure, amplify the similarity of nodes in a community [21, 1]. As said above the focal closure is considered being based on homophily, while the triadic closure means the tie formation between two acquaintances of an ego, when these meet at an occasion usually related to similarity in values or features. Thus even a relatively weak preference for homophilic relationships would tend to be amplified over time, via a cumulative advantage. As Kossinets and Watts discussed [21], this casts the question: To what extent can observed patterns of homophily be attributed to individual preferences and structural constraints? They also discussed that “A thorough answer to this question would require the use of simulation models, in which choice homophily as well as focal and cyclic closure biases could be systematically varied.” The purpose of this paper is to contribute to answering this question by studying the macroscopic consequences of the joint effect of focal and cyclic closure mechanisms.
To this end, we develop a model based on the weighted social network (WSN) model [7]. In the WSN model two main mechanisms for the tieformation process were used, namely local attachment (LA) and global attachment (GA), corresponding to triadic and focal closure, respectively, together with link strength reinforcement. With this minimal model the complex Granovetterian weighttopology relation of social networks could be successfully reproduced [7]. More recently, this model was generalized to temporal and multiplex networks [22, 23, 17].
In this paper the WSN model is generalized to enable homophilic relationships. We introduce features of the individuals (represented by the nodes of the network) in the spirit of the classical Axelrod model [24] by modeling the different traits of a node with components of an dimensional vector. The Axelrod model was intensively investigated by San Miguel and coworkers [25, 26, 27, 28] focusing mainly on the problem of cultural drift and fragmentation by applying dynamic rules to describe the changes of the features of the nodes upon interactions. In most cases fixed network geometries were applied, except in Ref. [27], where the effect of network adaptation was studied. Here we take a different approach by considering the features fixed and model the tie formation based on homophilic interactions as suggested by social network theory. We show that homophily has a major impact on the emerging network structure. Depending on the number of features there is a transition from a few large and homogeneous communities to many smaller and heterogeneous ones. The possible number of different values a feature can take is also a relevant parameter and we show that a scaling relationship between the average degree of nodes, , and is approximately valid.
The paper is organized as follows: After we define the model in the next section, we demonstrate that the model exhibits two regimes depending on the parameter values. We also show that the similarity between nodes is very much amplified in one of the phases, and show that a global segregation is observed. Finally, we present a summary and discussion in the last section.
Model
First we briefly review the original weighted social network (WSN) model proposed in Refs. [7, 17], which is a dynamic, undirected, weighted network model leading to a stationary state. We start with unconnected nodes and links between them are created and updated with the following three mechanisms:

global attachment (GA): A node is selected at random, its degree is . With probability it is connected with a new link of weight to a randomly chosen node. In the unlikely event that the two nodes are already connected, a new target node is chosen.

local attachment (LA): A randomly chosen node (with ) chooses one of its neighbors with the probability proportional to which stands for the weight of the link between nodes and . Then node chooses one of its neighbors but , say , randomly with probability proportional to and if nodes and are not connected, they are connected with probability with a link of weight . In addition, all the involved links increase the weights by , whether a new link is created or not.

link deletion (LD): Each link is removed from the system with probability at each time step. (Time is measured in sweeps.)
Sequential update is applied, first GA, then LA to the nodes and then LD to the links. Over the wide range of the parameters this model produces Granovetterian structures with strongly wired communities connected by weak ties [7, 17].
In order to incorporate the effect of homophily into the model we extend the WSN model in the following way. Similarly to the Axelrod model and its variants [24, 25, 26, 29, 30, 31, 27], we assign a vector of components to each node , corresponding to different features of the node (related e.g. to gender, ethnicity, language, religion, etc.). Each feature can take one of possible values. Then the feature vector of node is a set of integers: with for each . The size of the vector, , represents the social complexity of the population, i.e., the larger implies the greater number of cultural characteristics that are attributable to each individual. The number of traits per feature, , represents the heterogeneity of the population: The larger means the greater number of cultural options in the society. Although may depend on individuals and may be different for each feature in reality, for simplicity we assume all nodes to have the same and keep constant. These will be parameters of our model.
The feature values or traits of the nodes are initially chosen to be uniformly random and then kept unchanged. This setup is different from those in studies of the Axelrod model [24, 25, 26], where traits change as a consequence of interaction with similar neighbors. People in reality have both changeable and rigid traits, and the degree of the changeability may also depend on the individual and cultural properties. In this paper we study the case with fixed traits as it is the simplest setting to provide a baseline for better understanding of the impact of homophily on the network structure. We will show that segregation may occur even without a trait changing mechanism, just as a consequence of tie formation influenced by homophily.
In our homophilic WSN model, the network is updated similarly to the simple WSN model but the features of nodes are taken into account, as shown in Fig. 1. In each GA and LA step, a feature of the focal node is randomly chosen from features and it can make links only to the nodes sharing the same trait for the feature , i.e., only to the nodes satisfying .
More specifically, in a GA step, the node chooses randomly a node, say , sharing the same trait for the feature in the network. If the node has any neighbor sharing the same trait for the feature , then the link between and is created with probability . Otherwise, the link is created with probability one. The weight of the created link is given as . In the LA step, among ’s neighbors sharing the same trait for the feature with , choose a node, say , with probability proportional to their link weights. Then, the node chooses a random node, say , among ’s neighbors sharing the same trait for the feature but with probability proportional to the link weights. It implies that . All the links included in this triangle increase their weights by . If nodes and are not connected, a link between them is created with probability and its weight is set to .
Note that this model reduces to the original WSN model when , because every node pair shares the same trait, i.e., for all and . When , the network is decomposed into disjoint subgraphs since the nodes with different traits have no chance to be connected, while the model dynamics in each subgraph is equivalent to the original WSN model but with a smaller system size of .
Starting with a network with no links, GA, LA, and LD are applied in each time step. After a sufficiently long equilibration period, ranging from to Monte Carlo steps in our work depending on the relaxation time, statistically stationary states are reached.
Results
First, we show the dependence of the average degree on and in Fig. 2(a). For we find two regimes: The average degree decreases as increases up to a certain value, , then it increases for . Such nonmonotonic dependence is also observed for the clustering coefficient in Fig. 2(b). The clustering coefficient shows increasing and decreasing behaviors for and , respectively, hence it has a peak at around , implying that the network structure qualitatively changes at around .
The transition point depends on the WSN parameters but, more importantly, also on ; as it is clear from Fig. 2 the function is decreasing. Consequently, for a given a similar transition takes place as a function of thus one can define .
To investigate what is happening at the transition point, we measure how many of the features are shared between the neighboring nodes. We define the feature overlap between nodes and as the fraction of shared features:
(1) 
where is the Kronecker delta. Figure 3(a) shows the feature overlap averaged over the connected node pairs, , as a function of . stays close to in the low phase indicating that most of the links are created between the nodes having completely identical features. Hereafter, we denote the nodes having as “matching nodes”. On the other hand, as becomes greater than , the feature overlap quickly drops from unity and converges to a value slightly greater than . Due to the above results we denote the two regimes of and as “segregated” and “overlapping” phases, respectively, since nodes are segregated into looselyconnected homogeneous communities in the former phase, while the majority of the communities overlap with each other in the latter phase. The corresponding phase diagram is shown in Fig. 4(a).
Let us shed light onto the origin of the transition. Since the feature values are drawn from a uniform discrete distribution from , the number of nodes sharing a specific set of traits is on average . The transition occurs when becomes less than the typical degree for the original WSN model. For instance, the transition point for leads to , which is of the order of the average degree for the original WSN model with the same parameters. Therefore, we expect that is a scaling variable for the average density, see Fig. 4(b).
By the model definition, each node has a strong tendency to connect to matching nodes as long as a matching node can be found, i.e., in the segregated phase. In the overlapping phase, on the other hand, the chance to find matching nodes becomes small thus the feature overlap cannot be one anymore. The nodes compromise with partially matching nodes as few matching nodes exist in the whole population.
When , the average degree decreases as a function of as shown in Fig. 2(a). This is because the number of matching nodes , which is equal to the number of potential neighbors, decreases with . The opposite behavior is found in the other phase, where the nodes are connected to others that share at least one feature. The number of nodes sharing at least one feature is calculated as , which is an increasing function of . Since there is no chance for building communities of matching nodes for the nodes have more potential neighbors with increasing thus they will have higher degrees. The spike found in the clustering coefficient [Fig. 2(b)] can also be explained similarly. When gets closer to in the segregated phase, the number of potential neighbors decreases and eventually becomes the same order as the degree. Since nodes are partitioned into groups of matching nodes and its size is comparable to the average degree, the nodes within a group are mostly connected to each other, while the number of links bridging the groups is small. Thus, the network is mostly composed of cliquelike components, yielding a high clustering coefficient as well as a high assortativity (not shown).
This transition behavior cannot be simply explained by the naive assumption that the probability of making a link is proportional to . Let us consider a null model in which links are randomly created between nodes and with a probability proportional to . The expected number of links between the node pairs having is given by the binomial form as
(2) 
where is a coefficient for the probability of creating links. Therefore, the average feature overlap in the network is calculated as
(3) 
This null model predicts a smooth monotonic decrease without a transition with increasing [see Fig. 3(a)], which is significantly different from the simulation results and with which it agrees only asymptotically.
Next we test a model with GA and LD but without LA by setting . The results in Fig. 3(a) do not show a transition and reproduce a curve similar to the null model. Thus, a key ingredient for the transition is the presence of LA rule. In other words, LA promotes a much stronger preference in link formation between matching pairs than predicted by the null model. Indeed, the link reinforcement in LA plays here a key role as initially minor differences between the probabilities of choosing neighbors with perfect matching and partial overlaps get strongly amplified due to the positive feedback produced by the weightdependence at the selection of the neighbors. As a consequence, the links created by triadic closures are almost always between matching nodes in the segregated phase. This argument has been numerically confirmed by disabling the link reinforcement by setting while keeping a positive , as shown in Fig. 3(a).
While the above amplification effect exists in the overlapping phase too, it cannot fully develop as the low number of available matching nodes results in a partial overlap. For instance, even when nodes and are both connected to and differ in their traits from those of ’s in one feature only, when and get connected by cyclic closure their traits may differ already in two features. Thus, the feature overlap between neighbors significantly decreases when nodes cannot connect to matching nodes.
The fact that in the segregated phase indicates that the nodes in the network form communities, each of which mostly consisting of matching nodes. Whether this segregation sets in or not depends on the value of : At the network structure changes from large homogeneous, segregated clusters to smaller more heterogeneous groups. This transition can be directly observed by applying the community detection to the generated networks. We use the multilevel Infomap algorithm [16] to detect a hierarchical structure of the communities. Since we are interested in the global segregation, we pick up the communities on the top of the hierarchy, and their overlap is measured by the number of communities per node as shown in Fig. 3(b) as a function of . We clearly see two distinct regimes. When , we have homogeneous communities consisting mainly of matching nodes thus the number of communities per node is close to one, indicating that there is little overlap of the membership between communities. For the feature overlap between linked nodes is suppressed. As a consequence, the membership overlap between communities increases, which works against segregation. This is why we named the two phases as segregated and overlapping phases.
These differences are illustrated in Fig. 5(a) where the egocentric network of a typical node is shown, and clearly all acquaintances have identical traits, i.e., matching nodes. There are, however, some () nodes which connect these large homogeneous communities. A typical example of these connecting nodes is shown in Fig. 5(b). Near the transition point, as shown in Fig. 5(c), a node typically has a smaller degree and shares less feature with its neighbors than those in the segregated phase. Figure 5(d) shows a sample from the overlapping regime. The ego is a part of more distinct communities that have small feature overlap because of a large , as shown in Fig. 3(a). As nodes are part of more and more communities the number of communities per node increases as shown in Fig. 3(b).
Another way to investigate the segregation is to see information spreading on the network. We simulate the susceptibleinfected (SI) dynamics [32] on the networks assuming that the transmission probability is simply proportional to link weights, irrespective of the features of nodes. The information starts to spread from a randomly chosen node and the fraction of infected nodes are calculated until all the nodes become infected. The SI model is able to detect the influence of the topology and the link weights onto the spreading speed [20, 33]. Figure 6 shows typical dynamics of spreading on the networks for (segregated phase) and (overlapping phase) when . These two networks have been chosen because they have similar average degrees ( and , respectively) and average link weights ( and , respectively). Although these two networks have comparable local network statistics, the spreading for the segregated phase is much slower than that for the overlapping phase. In the segregated phase, the spreading is trapped in a community of matching nodes: It takes a longer time to convey information to a node having different features than for the overlapping phase, in which spreading between different groups occurs more easily.
Summary and Discussion
In this paper we have generalized the weighted social network model to incorporate homophily and to study its effect on the emerging network structure. We introduced features for each node and every feature could have different values. The GA and LA steps were made then dependent on the overlap between the nodes’ features. The model shows two phases depending on the parameters and . For a given value of and the nodes in the whole network are segregated into nonoverlapping communities, consisting mostly of matching nodes. On the other hand, when , such segregation is not observed. Nodes partially share features with their neighbors yielding significantly larger amount of overlapping communities. Thus, if the selection mechanism of acquaintances is based on very few features ( is small) and/or if traits a feature can take are limited only to few options ( is small), then the society will be segregated.
What is the relevance of this finding, considering that in reality every person has a large number of features? We note that the features are not equally important in tie formation. There are cultures and situations, in which a few features get extreme importance. This is the case, in times of sharpening political situations, turmoils or wars. Another relevant example can be social groups in the online social networks, where people get into contact based on very few features. This (especially if amplified by the mentioned political conditions) leads then to the wellknown effect of “echo chambers” [34] and our results indicate that perceiving the complexity of the world or allowing for more subtle opinions can break up these bubbles.
The overlapping phase is thus a more realistic model of a normal society where one has acquaintances of different gender, age, location, political view, sport preference, etc. Each of these traits has its own communities in which there is no perfect feature overlap. For example, why should we do sports with people of the same political preference?
We demonstrated that the feature similarity between nodes can be much more amplified if the local attachment mechanism is taken into account. In other words, the homophily principle and the local attachment mechanism including reinforcement have a joint effect in increasing segregation in social networks. Our study implies strategies to avoid segregation and splitting of the society by maintaining cultural diversity.
The relationship between homophily and segregation has been recognized long ago. In the Schellingtype models [2, 35] people change their position in a fixed geometry (usually a lattice) to better satisfy their preferences. In another approach to social segregation, opinion dynamics is used so that similar people can influence each other [24, 25, 26, 36, 37]. In these models similarity is amplified by the changes in the opinions (or features). The model proposed in this paper has a different mechanism. Our aim was to investigate how homophily influences the social network if such elementary tie formation mechanisms like focal and cyclic closure (mimicked by our model as GA and LA, respectively) are taken into account. Our conclusion is that social segregation may be intensified by LA and link reinforcement even without changing the features in interactions. For more realistic modeling, we can of course think of a combined model incorporating both the feature updates and the triadic closure with reinforcement. Although we leave it for future studies, we expect that the segregation is observed in a wider range of parameter space as the feature updates promote segregation in general.
There are many further ways to extend our model and make it more realistic. A possible extension is to allow different values for each feature. This is a much more realistic setting and will be necessary, when comparing with a data. Also, individual differences in the number of features could be considered.
Acknowledgements
Y.M. acknowledges support from MEXT as “Exploratory Challenges on PostK computer (Studies of multilevel spatiotemporal simulation of socioeconomic phenomena)” and from Japan Society for the Promotion of Science (JSPS) (JSPS KAKENHI; grant no. 18H03621). H.H.J. acknowledges financial support by Basic Science Research Program through the National Research Foundation of Korea (NRF) grant funded by the Ministry of Education (2015R1D1A1A01058958). K.K. acknowledges the Rutherford Foundation Visiting Fellowship at The Alan Turing Institute, UK. The systematic simulations in this study were assisted by OACIS [38]. Y.M., H.H.J., J.T., and J.K. are thankful for the hospitality of Aalto University.
Author Contributions
All authors conceived and designed the experiments. Y.M. and J.T. conducted the experiments and analyzed the results. All authors reviewed the manuscript.
Data Availability
All data generated or analyzed during this study are included in this published article.
Competing Interests
The authors declare no competing interests.
References
 McPherson, M., SmithLovin, L. & Cook, J. M. Birds of a feather: Homophily in social networks. \JournalTitleAnn. Rev. Soc. 27, 415–444 (2001).
 Schelling, T. C. Models of segregation. \JournalTitleAm. Econ. Rev. 59, 488–493 (1969).
 Karimi, F., Génois, M., Wagner, C., Singer, P. & Strohmaier, M. Homophily influences ranking of minorities in social networks. \JournalTitleSci. Rep. 8, 11077 (2018).
 Lee, E., Karimi, F., Jo, H.H., Strohmaier, M. & Wagner, C. Homophily explains perception biases in social networks. \JournalTitlearXiv preprint arXiv:1710.08601 (2017).
 Halberstam, Y. & Knight, B. Homophily, group size, and the diffusion of political information in social networks: Evidence from twitter. \JournalTitleJournal of Public Economics 143, 73–88 (2016).
 Kossinets, G. & Watts, D. J. Empirical analysis of an evolving social network. \JournalTitleScience 311, 88–90 (2006).
 Kumpula, J. M., Onnela, J.P., Saramäki, J., Kaski, K. & Kertész, J. Emergence of communities in weighted networks. \JournalTitlePhys. Rev. Lett. 99, 228701 (2007).
 Newman, M. E. J. Networks: An Introduction (Oxford University Press, 2010).
 Barabási, A.L. Network Science (Cambridge University Press, 2016).
 Fortunato, S. & Hric, D. Community detection in networks: A user guide. \JournalTitlePhys. Rep. 659, 1–44 (2016).
 Palla, G., Derényi, I., Farkas, I. & Vicsek, T. Uncovering the overlapping community structure of complex networks in nature and society. \JournalTitleNature 435, 814 (2005).
 Girvan, M. & Newman, M. E. Community structure in social and biological networks. \JournalTitleProc. Nat. Acad. Sci. 99, 7821–7826 (2002).
 Lancichinetti, A., Fortunato, S. & Kertész, J. Detecting the overlapping and hierarchical community structure in complex networks. \JournalTitleNew J. Phys. 11, 033015 (2009).
 Ahn, Y.Y., Bagrow, J. P. & Lehmann, S. Link communities reveal multiscale complexity in networks. \JournalTitleNature 466, 761–764 (2010).
 Rosvall, M. & Bergstrom, C. T. Maps of random walks on complex networks reveal community structure. \JournalTitleProc. Nat. Acad. Sci. 105, 1118–1123 (2008).
 Rosvall, M. & Bergstrom, C. T. Multilevel compression of random walks on networks reveals hierarchical organization in large integrated systems. \JournalTitlePLoS ONE 6, e18209 (2011).
 Murase, Y., Jo, H.H., Török, J., Kertész, J. & Kaski, K. Modeling the role of relationship fading and breakup in social network formation. \JournalTitlePLoS ONE 10, e0133005 (2015).
 Granovetter, M. S. The strength of weak ties. \JournalTitleAm. J. Sociol. 78, 1360–1380 (1973).
 Onnela, J.P. et al. Analysis of a largescale weighted network of onetoone human communication. \JournalTitleNew J. Phys. 9, 179 (2007).
 Onnela, J.P. et al. Structure and tie strengths in mobile communication networks. \JournalTitleProc. Nat. Acad. Sci. 104, 7332–7336 (2007).
 Kossinets, G. & Watts, D. J. Origins of homophily in an evolving social network. \JournalTitleAm. J. Soc. 115, 405–450 (2009).
 Jo, H.H., Pan, R. K. & Kaski, K. Emergence of bursts and communities in evolving weighted networks. \JournalTitlePLoS ONE 6, e22687 (2011).
 Murase, Y., Török, J., Jo, H.H., Kaski, K. & Kertész, J. Multilayer weighted social network model. \JournalTitlePhys. Rev. E 90, 052810 (2014).
 Axelrod, R. The dissemination of culture: A model with local convergence and global polarization. \JournalTitleJ. Confl. Res. 41, 203–226 (1997).
 Centola, D., GonzalezAvella, J. C., Eguiluz, V. M. & San Miguel, M. Homophily, cultural drift, and the coevolution of cultural groups. \JournalTitleJ. Confl. Res. 51, 905–929 (2007).
 Vazquez, F., GonzálezAvella, J. C., Eguíluz, V. M. & San Miguel, M. Timescale competition leading to fragmentation and recombination transitions in the coevolution of network and states. \JournalTitlePhys. Rev. E 76, 046120 (2007).
 Min, B. & San Miguel, M. Fragmentation transitions in a coevolving nonlinear voter model. \JournalTitleSci. Rep. 7, 12864 (2017).
 Battiston, F., Nicosia, V., Latora, V. & San Miguel, M. Layered social influence promotes multiculturality in the axelrod model. \JournalTitleSci. Rep. 7, 1809 (2017).
 Vázquez, F. & Redner, S. Nonmonotonicity and divergent time scale in axelrod model dynamics. \JournalTitleEurophysics Letters 78, 18002 (2007).
 Gandica, Y., Charmell, A., VillegasFebres, J. & Bonalde, I. Clustersize entropy in the axelrod model of social influence: Smallworld networks and mass media. \JournalTitlePhys. Rev. E 84, 046109 (2011).
 Tilles, P. F. & Fontanari, J. F. Diffusion of innovations in axelrod’s model. \JournalTitleJ. Stat. Mech. 2015, P11026 (2015).
 PastorSatorras, R., Castellano, C., Van Mieghem, P. & Vespignani, A. Epidemic processes in complex networks. \JournalTitleReviews of Modern Physics 87, 925–979 (2015).
 Karsai, M. et al. Small but slow world: How network topology and burstiness slow down spreading. \JournalTitlePhys. Rev. E 83, 025102 (2011).
 Barberá, P., Jost, J. T., Nagler, J., Tuckler, J. A. & Bonneau, R. Tweeting from left to right: Is online political communication more than an echo chamber? \JournalTitlePsychological Science 26, 1531–1542 (2015).
 Henry, A. D., Prałat, P. & Zhang, C.Q. Emergence of segregation in evolving social networks. \JournalTitleProc. Nat. Acad. Sci. 108, 8605–8610 (2011).
 Deffuant, G., Neau, D., Amblard, F. & Weisbuch, G. Mixing beliefs among interacting agents. \JournalTitleAdv. Compl. Syst. 3, 87–98 (2000).
 Del Vicario, M., Scala, A., Caldarelli, G., Stanley, H. E. & Quattrociocchi, W. Modeling confirmation bias and polarization. \JournalTitleSci. Rep. 7, 40391 (2017).
 Murase, Y., Uchitane, T. & Ito, N. An opensource job management framework for parameterspace exploration: Oacis. In Journal of Physics: Conference Series, vol. 921, 012001 (IOP Publishing, 2017).