Action-based Modeling of Complex Networks

Complex networks can model a wide range of complex systems in nature and society, and many algorithms (network generators) capable of synthesizing networks with few and very specific structural characteristics (degree distribution, average path length, etc.) have been developed. However, there remains a significant lack of generators capable of synthesizing networks with strong resemblance to those observed in the real-world, which can subsequently be used as a null model, or to perform tasks such as extrapolation, compression and control. In this paper, a robust new approach we term Action-based Modeling is presented that creates a compact probabilistic model of a given target network, which can then be used to synthesize networks of arbitrary size. Statistical comparison to existing network generators is performed and results show that the performance of our approach is comparable to the current state-of-the-art methods on a variety of network measures, while also yielding easily interpretable generators. Additionally, the action-based approach described herein allows the user to consider an arbitrarily large set of structural characteristics during the generator design process.

• ABNG-D(·): This synthesis algorithm modifies ABNG-PA(·) by deleting edges from a starting network having more edges than the target. Again, the input parameter denotes the number of edge deletion queries in each time step during network synthesis. The computational complexity of this algorithm is proportional to the number of edges that need to be deleted from the starting network.
• ABNG-R: Another local operation that actions can perform is rewiring i.e. changing one end of an already existing edge. The rewiring technique is used in dk-random graphs [21] to sample from the entire ensemble of networks. In context of ABNG, ABNG-R starts with the original target network and rewires its edges using actions. Every edge of each node is queried for an action-based rewiring resulting in a total of 2m t action calls.
• ABNG-C: ABNG-C corresponds to a degree sequence preserving extension of ABNG-PA. It can be seen as a constrained version ABNG-PA where the degree sequence of the target network is preserved.
Each node i can add edge (v i , v j ) under the condition that the degree sequence is not violated. Joint degree distribution preserving synthesis algorithms is also a possibility. It can be seen that starting from a network with no edges, this synthesis algorithm needs to add m t edges leading to computational complexity directly proportional to m t . Table S1: Synthesis algorithms for ABNG: ρ is the probability an edge will be added by a given action matrix, and m s is the number of edges in the starting network.

Synthesis Algorithm Description
Action calls ABNG-PA(·) Action-based addition of edges to a starting network with m s m t It should be noted that the computational complexity of a synthesis algorithm can be upper-bounded by considering the complexity of most expensive action in the Action set. Also, as highlighted in Table S1, the computational complexity of a synthesis algorithm depends on the action matrix being used as input or ρ, which is the probability an edge will be added by a given action matrix.

S4 A Worked Example
In ABNG, nodes synchronously create edges in discrete time steps. For example, at t = 0 we are given a sparse starting network and actions for each node are evaluated based on the network given at t = 0, followed by creation of new networks at t = 1, which are used for evaluation of edges to be added at t = 2 and so on.
Also, the new edge is randomly chosen, independently of the others, according to a distribution depending only on a finite set of actions. To better understand the working of the ABNG framework, let us consider a simple example. Consider a starting network G 0 (shown in Figure S1a) at t = 0 with adjacency matrix A 0 , an action matrix M = [ 0.7 0.3 ], and ABNG with two hypothetical actions a 1 and a 2 . We compute matrices A 1 0 (for this example, it is assumed that a 1 is an action based on preferential attachment on node degree and is used for generating the matrix A 1 0 ) and A 2 0 , where the i th row of A 1 0 corresponds top i obtained by using action a 1 . The networks synthesized by ABNG (at t = 1) can be sampled using A 1 = 0.7A 1 0 + 0.3A 2 0 + A 0 , where each element of A 1 corresponds to the probability of existence of an edge. Three networks synthesized using A 1 are shown in S1b. The networks synthesized at t = 1 can now be used as starting networks for t = 2. We gain several key insights from this example: • The actions a 1 and a 2 belong to two different categories of actions, namely probabilistic and deterministic. In a deterministic action, a node v i selects another node v j to create an edge, whereas in a probabilistic action there exist probabilities of connecting to different nodes.
• For the rows of A 1 0 , it can be seen that the sum is < 1. This means it is feasible that a node might not create an edge even after choosing an action.
• It can be clearly seen that the condition for no multi-edges or self loops is enforced. All corresponding entries have zeros in A 1 0 and A 2 0 .
• It should be noted that the matrices must remain symmetric in case of undirected networks.    Figure S1: Network synthesis process of ABNG: S1a The initial starting network G 0 having adjacency matrix A 0 . S1b Three networks synthesized at t = 1 as explained in the example. It must be noted that the networks shown here have different edges.

S5 Algorithms
Algorithm S1 ABNG 1: Input: M, a q × (k + 1) matrix satisfying Equation 1, a starting network G = (V , E ), a target network T = (V, E) and a set A of k actions 2: while |E | < |E|, for every node v i in G do 3: choose a l for v i wp P * il l = 1, · · · , k 5: end for 6: end for is an estimation of the expected value.

S6 Proposed Implementation
In this section, a simple yet effective way of implementing the abstract concept of ABNG discussed in Section 2 is presented. The goal is to compute (estimate) an action matrix that optimizes ABNG based on its average case performance as defined by the optimization problem in Equation 1. To solve this multi-objective search problem, we implement Pareto Simulated Annealing (PSA) [11], as it is known to be a useful metaheuristic capable of global optimization in a large search space in a fixed amount of time (or iterations). It is also beneficial due to the fact that only one evaluation of the objective function is required at each iteration when compared with population-based GA algorithms, which require an evaluation for each member of the population.
The GP system developed in [4] was used in a meta-analysis in [14] to evaluate six network centrality measures. Results indicated that of the examined centrality measures, the degree distribution, betweenness centrality, and PageRank were the most effective for quantifying the (dis)similarity between the target and the synthesized networks. We will use these three measures; however, the framework allows for any user-desired measures.
The implementation of ABNG starts with the assumption that each node has the same probability distribution over actions. In other words, we assume that all nodes are homogeneous with respect to their preference over actions that in turn influences their choices for creating specific edges in the network. This implies that all rows in P are identical and the action matrix has dimensions 1 × (k + 1). Additional rows are dynamically added to M. This is discussed in detail with the optimization algorithm in this section.
The PSA approach explores the solution space by increasing (or decreasing) individual elements (at each iteration, an element is chosen uniformly at random) of the action matrix to find better solutions, while accepting worse solutions with a probability decreasing with the number of iterations.

S6.1 Assumptions
Before going into further details of the optimization algorithm, we list the simplifying assumptions in the current implementation of ABNG: 1. Input network types: It is assumed that all of the target and synthesized networks are simple graphs, i.e undirected with no self edges. The networks considered for experiments were also unweighted and unlabelled. This does not imply that ABNG is not applicable to such networks.

Network objectives:
No community-specific objectives were considered, although communities are likely in real-world networks. If needed, special objectives can be added to the optimization framework to synthesize networks with more specific community structure.
3. Starting network: The current implementation of ABNG needs a starting network with n nodes as input. Furthermore, this starting network cannot be empty because some actions can potentially become undefined because of lack of any network characteristics. For example, an action based on preferential attachment according to degree of a node would essentially be equivalent to randomly selecting a node in case of an empty network because each node will have a degree of 0. To tackle these issues, we create G 0 , a starting sparse network (edges nodes), as input.

Fitness:
As seen in Algorithm S1, even though the network is built over multiple iterations per node, ABNG does not evaluate characteristics for these interim networks. That is, only the target network characteristics are considered by the evaluation function.
5. Types of actions: As discussed in Section 5, the synthesis algorithm in the current implementation is restricted to adding edges only. This is a reasonable assumption because we only consider static target networks, which can be synthesized using this restricted subset of actions.
6. Stopping criteria: The network synthesis process of Algorithm S1 is terminated when the number of edges in the network being synthesized is equal to the number of edges in the target network.

S6.2 Optimizing the Action Matrix
Algorithm S3 Pareto Simulated Annealing  Owing to the multi-objective nature of the problem, we need to maintain a set S of potentially efficient and The optimization procedure could terminate at a local optima or some other non-optimal stationary point. To help overcome this issue, we use multiple starting points for the optimization process, i.e. the procedure starts with more than one M 0 . In the current implementation, we attempt to find the most simple generator (in terms of how nodes make decisions for connecting to other nodes) for the target network. In terms of the action matrix, the simplest generator can be obtained when all nodes have a common mechanism (same probability distribution over actions) for forming edges, and hence M is assumed to have dimensions 1 × (k + 1). This is followed by dynamically adding more 'node types' (or rows in the action matrix) to the current solutions. In the algorithm implementation, once a Pareto front (or a set of potentially efficient solutions) S is found for M : 1 × (k + 1), a solution is picked at random from S and a new row (generated at random from D) is added to M, making it a 2 × (k + 1) action matrix and so on. For initialization,P containing the probability of choosing rows in M is generated uniformly at random. For a q × (k + 1), q > 1 action matrix, we assume that only the newly added q th row andP need to be optimized because the remaining rows have been optimized for the target network in previous steps. In other words, this procedure of adding new rows to M implies that the previous rows contain information about how nodes selected actions when their choices were restricted and adding new rows will provide them with more choices.P allows nodes to choose among various rows, which may change when new rows are added and hence needs to be considered in the optimization framework. Though this procedure might restrict the search space for the action matrix, it suffices to provide some insights about the applicability and ability of the action-based approach to synthesize networks. As described in Section 5, each row of the action matrix reflects the mixed strategies used by nodes for making connections, and the aim is to find the minimum number of such rows that corresponds to the smallest parameter space capable of synthesizing the target network using ABNG.
More rows are added to the action matrix until at least one of the following criteria are met: • The newly added q th row has probability close to zero in the converged solutions, i.e.P q ≈ 0. This implies that the network structure can be explained equally well using a smaller action matrix. A threshold value of 0.05 is used in the present implementation.
• The newly added q th row is similar to one of the previous rows i.e., M(q, : Similarity between two rows can be defined using any vector similarity measure. This implies that the two strategies are practically equivalent. In case such an event occurs during optimization, theP values for both the rows are added (and assigned to the initial row) and the copy row is deleted i.e.
• Solutions from a smaller action matrix strictly dominate the new solutions, i.e. S(M b×(k+1) ) This means adding new rows does not improve the quality of the synthesized networks for the given set of actions.
The modified version of Pareto Simulated Annealing used here also incorporates an adaptable number of maximum iterations (Algorithm S3 lines [16][17][18][19][20]. After every 100 iterations, the previous 100 solutions are checked for improvement in any of the objectives. If S z = S z−100 , 100 more iterations are allowed, otherwise the process is terminated at a maximum of 1000 iterations. The strategy resulted from the peculiar behavior of the optimization process that can be observed in Figure S7, where convergence of solutions is observed before reaching the maximum number of iterations. This adaptable approach aids the algorithm in identifying potential local optima and consequently stopping the optimization process.

S6.3 The Action Set
An action for a node v i builds edges, all of which have v i as one endpoint. An action provides a well-defined strategy for selecting the other end v j of edge (v i , v j ). This implies that the actions in ABNG can only allow local topological changes in the network as a node is assumed to only build edges to other nodes but cannot create edges between other nodes. Every action for node v i returns another node v j with probabilityp j i to form edge (v i , v j ). So, for network G = {V, E}: where an edge This approach also enables an action (and hence ABNG) to insert edges that can be directed, weighted, self edges, etc. In this paper, it is assumed that actions are limited to adding a single undirected edge. Actions for removing, rewiring, or adding multiple edges can be easily added to the framework while following the guidelines given below: • Adding, removing, rewiring edges should be based on some criteria rather than being random. They can be inspired from some real-world phenomena for forming connections or based on some common criteria for application to particular domains.
• A node should only have the ability to change its local structure.
• The action should reveal some information about the network construction mechanism in the target network.
This raises the question: How many actions should be considered for synthesis of a network? There is no direct answer for such a question, however, having too many actions will increase the number of parameters (size of the action matrix) and hence make solving the optimization problem more difficult and time-consuming. This might even lead to degenerate actions and consequently complicate the generator.
Actions form an integral part of ABNG. They collectively serve as the mechanisms responsible for network synthesis. Hence, choosing a holistic or sufficient set of actions is crucial for implementing ABNG. The current implementation of ABNG uses eight actions belonging to four different categories: • Preferential attachment using network centrality measures -This is analogous to the Barabási-Albert Model [5], which uses node degree as an action to connect to different nodes. It connects to important nodes with higher probability, where importance is calculated using network centrality measures. We use degree (PAD), average neighbor degree (PAND), PageRank (PAPR) and betweenness (PAB) as centrality measures for four different actions. An important property of this set of actions is thatp is the same for all the nodes.
• Triadic closure -This action (TC) connects a node to another node that is a neighbor of its neighbor.
It captures the phenomenon: a friend of my friend is also my friend. Also, triangles form a commonly encountered structure in real-world networks [18]. In case a node has multiple second neighbors, each one has an equal probability of getting selected.
• Similarity-based actions -These provide a basis for nodes to connect to similar nodes, another phenomenon observed in most real-world networks [18]. Inverse log-weighted (SLW) and Jaccard similarity (SJ) measures are used in this implementation. Nodes connect to the most similar unconnected node for preventing multiple edges.
• There is another action (NA) that does not connect the current node to any other node. As we can see from Algorithm S1 , ABNG visits every node in the network to form new edges, and this action exempts a node v i from making a connection, i.e.p i = 0.

S7.2 ABNG for Real Networks
19 real-world networks were also considered and details are shown in the table. Real-world networks will likely not be simply described like the human-devised models in the previous experiments. Radar plots for the networks synthesized using action matrices obtained as solutions for these real-world networks can be seen in Figure S3. The estimated action matrices can potentially provide insights about the structure of these networks like, how many types of nodes exist in the network, how they weigh actions to form edges etc. Description of the action-based model for five of these networks can be seen in Section 3. Figure S3 also shows statistics when some other network generators, namely Chung-Lu [9,8] and ERGM [2], were used to fit the target network. From the comparison it can be concluded that: • Networks synthesized using ABNG show the most resemblance to the target networks when evaluated based on the network properties considered here.
• Unlike other models, the output parameters obtained using ABNG (action matrices) can provide a compact representation of structure of the target networks.

S7.3 Scaling with Network Size
An experiment was also performed to get empirical insights about scalability of the synthesis algorithm described in Algorithm S1. As described in Section S3, the complexity of any given synthesis algorithm depends on the input action matrix together with the size of the network. To understand the relation between network synthesis time and size of the network (number of nodes), an experiment was performed where each trial used a different action matrix and the size of the network was increased. Also, the mean degree of the networks was kept constant (d = 6) as the number of nodes in the network was varied. Mean degree of 6 was chosen based on observations of mean degree of real-world networks shown in Table S2.
Results shown in Figure S4

S7.4 Starting Network Variations
In the experiments conducted so far, it is assumed that the starting network is obtained by sampling 0.7 × n edges from the target network. In this section, different proportion of links are sampled from the target network to see how the generator performs when provided with different starting network. For this, we consider four different target networks and five different fraction of links, as shown in Figure S5. The In earlier experiments, a fixed starting network was used for synthesis throughout the optimization process (see Figure 1). While varying the fraction of links in the starting network, we also relaxed the assumption of using a fixed starting network and sampled a different starting network each time the algorithm was used to synthesize a network. Results indicate that ABNG can synthesize networks having similar Q when the assumption of a fixed starting network is relaxed, hence providing evidence about the robustness of the action-based approach.

S7.5 Spectral Goodness of Fit
A new statistic to evaluate how well a network generator explains the structure of the pattern of ties in the target network was proposed in [25]. The current version of the Spectral Goodness of Fit (SGOF) statistic is limited to only unlabeled and undirected networks, which matches with the type of networks synthesized using ABNG. Because of its simplicity, we use SGOF as a proxy measure for goodness of fit for the synthesized networks and do not consider it in the optimization process of ABNG.
The approach calculates the Euclidean Spectral Distance (ĒSD T ,G = ||λ T −λ G ||), whereλ T andλ G are the normalized spectra (of the Laplacian) of networks T and G. The spectral goodness of fit (SGOF) can be then obtained by: where N is the null model. For SGOF calculations, the Erdös-Rėnyi model is used as the null model. The SGOF measures the amount of observed structure (T ) explained by a fitted model ({G 1 , G 2 , · · · }), expressed as a percent improvement over a null model, where structure means deviation from randomness [25]. SGOF is bounded above by one, which means that the network generator (or fitted model) exactly describes the target network. Similarly, an SGOF of zero means that synthesized networks are only as good as the random networks, whereas a negative value signifies that the null model is a better approximation of the target network as compared to networks synthesized using the network generator. Figures S6a and S6b show the SGOF values obtained from 100 networks synthesized using different network generators for both human-devised and real-world networks. The networks for ABNG are synthesized using the action matrix corresponding to the point closest (based on 1-norm distance) to the origin in the Pareto front. Note that SGOF values close to zero are observed for the Erdös-Rėnyi network because it is itself used as the null model. To maintain consistency of using a dissimilarity metric ranging between 0 and 1, we transform the obtained SGOF value by using the function y = 1 − 2 x−1 , where y gives us the transformed dissimilarity and x is the SGOF value obtained using S4. This transformation is used in the heat maps of Table 2.

S7.6 Analyzing the Action Matrix
Here, the evolution of the action matrix and the objectives during the PSA iterations are examined to get a better understanding of the process of learning an action matrix in ABNG. Figure S7 illustrates two examples representative of typical evolution of solutions for 1-row and 2-row action matrices when using PSA for optimization in ABNG. The optimization process shows typical characteristics observed in evolutionary multi-objective optimization algorithms, where a lot of improvement is seen in the objective space in the first few iterations and the solutions seem to converge in the later iterations. Also, it can be seen that the best solutions are found before reaching the maximum number of iterations and can be observed when the curves become flat in Figure S7.  Figures S8 and S9 show different examples for the evolution of a 1 × 8 action matrix. Figure S9a shows a plot when the optimization process starts with a randomly generated starting point, while Figure   S9b shows the evolution of action matrix when the starting solution had equal probability for each action.
Similar to Figure    To visualize the diversity of the solutions based on the size of the action matrix, a 3D plot ( Figure S10) is shown for the three objectives used in the optimization process for the brain network with correlation threshold of 0.7. The example is representative of the distribution of Pareto optimal points obtained when ABNG is optimized for a target network. It is observed that the solutions are more spread out (or scattered) for smaller action matrices (1-row and 2-row) and the solutions seem to be more concentrated at a particular region when considering larger action matrices. Sensitivity of the action matrix was also tested for a 2-row action matrix obtained for the forest fire model network. In this scenario, the same approach was used to separately perturb the first row, second row andP associated with the action matrix. Again, the radar plots of Figure S12 provide evidence for the stability of the outcome even in the case of a bigger action matrix.
2. The one-at-a-time (OAT) approach does not fully explore the input space since it does not take into account the simultaneous variation of input variables. This means that the OAT approach cannot detect  properties that were used as objectives for optimization are underlined. Each line corresponds to average of 20 networks synthesized using a OAT perturbed action matrix where S12a is for the first row, S12b is for the second row and S12c is forP .
the presence of interactions between input variables. Instead, the second test varies the probabilities of all the actions simultaneously. Radar plots for this version of the sensitivity analysis can be seen in Figure S13 and S14 for the 1-row and 2-rows cases respectively. For each network, the test is performed five times with distinct variations in the action matrix. Clearly, even using this method there is a lot of overlap in the synthesized network properties for the different variations of the action matrix, especially for the properties considered as objectives for the optimization. The same is true for the case of a 2-row action matrix considered in Figure S14.
Results for both the cases reflected the robustness of the solutions obtained by showing that making small perturbations to the action matrix had little effect on the output of the synthesized networks. This provides preliminary evidence for the continuity in mapping the action matrix to the objective space and that the quality of the synthesized target networks is robust to small changes in parameters.  properties that were used as objectives for optimization are underlined. Each line corresponds to average of 20 networks synthesized using a simultaneously perturbed action matrix where S14a is for the first row, S14b is for the second row and S14c is forP .

S7.8 Datasets and Packages
The implementation of Action-based network generator has been done using the "igraph" package in R.
This involved synthesizing target networks from other human-devised network generators and performing statistical analysis on various networks. The "spectralGOF" package in R was used to compute SGOF values for different networks. Following real world networks were used: network of word adjacencies [19], US politics books sold on Amazon [1], a network of coappearances [16], network of American football games [12], collaboration network between Jazz musicians [13], a social network of dolphins [17], brain networks (with different correlation cut-offs) [6], two protein networks [7], a network of yeast protein interactome [15], three networks obtained from the Biogrid repository [26], US Airport network [10], Norwegian boards network [23], human protein interaction network [22], social network from an online community for students at University of California, Irvine [20], and a network representation of the topology of the Western States Power Grid of the United States [27].
The protein networks 1php and 1qop were obtained from the Protein Data Bank [7]. The .pdb files obtained from this database contains information about all atoms composing a given protein. This can be used to obtain contact maps containing key relations from the protein structure. The C-alpha atoms were chosen from the .pdb files and a contact map was obtained using a threshold of 8 Å, i.e. if the distance between two atoms i and j is less than 8 Å, then the undirected link (v i , v j ) exists. The networks Biogrid FRET, Far Western and Dosage Lethality were obtained from the biogrid database available at [26]. Brain fMRI data was obtained from the NITRC KKI dataset available at The Neuro Bureau ADHD-200 Preprocessed Repository [6]. The parcellation scheme described in [24] was used to define nodes for network analysis.
Correlation matrix for patient 1019436 was used at correlation cutoffs 0.7, 0.6 (brain 2), and 0.55 (brain 3) to generate the respective networks (largest component of the network was used as the target network).
The Norwegian boards network [23] shows interaction among board members of public companies in Norway as obtained from the data in August 2011. The human protein interaction network was obtained from [22] and the largest connected component was used as the target network. The US Airport network [10] shows connections between airports (nodes) if there is a direct flight between them.