Dynamic characteristics rather than static hubs are important in biological networks

Biological processes are rarely a consequence of single protein interactions but rather of complex regulatory networks. However, interaction graphs cannot adequately capture temporal changes. Among models that investigate dynamics, Boolean network models can approximate simple features of interaction graphs integrating also dynamics. Nevertheless, dynamic analyses are time-consuming and with growing number of nodes may become infeasible. Therefore, we set up a method to identify minimal sets of nodes able to determine network dynamics. This approach is able to depict dynamics without calculating exhaustively the complete network dynamics. Applying it to a variety of biological networks, we identified small sets of nodes sufficient to determine the dynamic behavior of the whole system. Further characterization of these sets showed that the majority of dynamic decision-makers were not static hubs. Our work suggests a paradigm shift unraveling a new class of nodes different from static hubs and able to determine network dynamics.


Introduction 28
Recently, there has been a shift in biological research from studying single compounds 29 towards complex regulatory networks and their dynamic behavior (Kitano, 2002). Hence, 30 network analysis has emerged as a powerful tool to understand complex biological 31 processes (Barabási et al., 2011). In accordance, diseases or the development of cancers 32 are rarely a consequence of a mutation of a single component within a network, but rather its 33 perturbation (Barabási et al., 2011; Barabási and Oltvai, 2004). Thus, to understand the 34 development of diseases, we have to follow the network dynamics (Sverchkov and Craven,35 2017). However, a significant limitation is that interaction graphs cannot adequately 36 represent temporal changes of a pathway or network. Among models that investigate 37 dynamics, Boolean network models (Kauffman, 1969) can approximate simple features of 38 interaction graphs integrating also dynamics. Here, known regulatory interactions can be 39 directly formulated as Boolean functions (4) (Naldi et al., 2015) which can then be used to 40 obtain the dynamics of the network. The dynamic simulation of a Boolean network leads to 41 periodic sequences of states, called attractors, describing the long-term behavior of the 42 model (Xiao, 2009). Attractors depict the activity of each component within the system at a 43 specific point in time and can be associated with phenotypes (Kauffman, 1993). In this 44 sense, the sequence of attractor states and the compounds involved can be seen as a 45 model of the activities specific for a biological state.

47
Biological networks exhibit a scale-free topology ( Barabási and Oltvai, 2004). That is, a small 48 number of components called hubs are highly connected, whereas the majority of nodes 49 have few connections (Guimerà and Nunes Amaral, 2005). Often hubs are seen as the 50 master regulators of biological processes (Borneman et al., 2006; He and Zhang, 2006) 51 whose removal correlates with lethal phenotypes (Jeong et al., 2001). However, studies 52 state that it is almost impossible to eliminate the activity of hub nodes by targeted inhibition 53 with small molecules (Lu et al., 2007;Song et al., 2017). The most striking difference 54 between both observations is the fact that one side considers the complete loss of proteins 55 while the other side thinks about inhibition. Again, these studies are mainly based on static 56 interaction graphs. Since master regulators are supposed to affect the behavior in time of 57 biological processes, it is appealing to investigate them from a perspective involving network 58 dynamics. Nevertheless, investigating and detecting master regulators by their impact on 59 network dynamics comes with the limit given by the complexity of the system. This might 60 lead to unfeasibility of dynamic investigation or to reduction procedures, that might in turn 61 alter the dynamic behavior. Hence, strategies to bridge the gap needed to unravel the 62 dynamics of biological networks are of crucial interest. Which nodes are then crucial for the 63 dynamics of a biological network? How can they be detected? Are these dynamic regulators 64 actually the already widely studied static hubs? Can these nodes be relevant as intervention 65 targets? To uncover these questions, we studied the dynamics in biological systems ( Figure  66 1).

68
Here, we present techniques on how to identify nodes sufficient to determine the dynamic 69 behavior of the network and discuss the implications of targeting these nodes for intervention 70 experiments 71 72 73

74
Small sets of compounds determine the dynamics 75 The state of a Boolean network is given by a vector of values (0/1) assigned to their nodes 76 (n). A state change is a result of the application of a Boolean function for each node 77 (Kauffman, 1969 for the exhaustive search to calculate an implicant set of size k of a network with n nodes. 97 Furthermore, our analyses revealed that the typical size of minimal implicant sets is small 98 (Figure 2A). The cardinality of implicants ranges from 1 to 9 for an average size of 19 nodes 99 per network. Here, the exhaustive approach identified a mean of 4.4 implicant nodes and the 100 heuristic a mean number of 5.1 implicant nodes. Only a low correlation between the network 101 size and the size of the implicant set could be found (Pearson correlation at 0.32 using the 102 exhaustive strategy and 0.39 using the heuristic strategy). Furthermore, the relation between 103 the set size of the implicant set and the network size is not linear ( 2 = 0.104 using the 104 exhaustive strategy and 2 =0.153 using heuristic strategy).

106
Majority of implicants are lowly connected but behave like hubs 107 Next, we studied the connectedness of the dynamic influencing nodes of the implicant sets 108 to uncover if they are hub nodes. To ensure that the links of nodes within the considered 109 Boolean network models are comparable to their biological representatives, we initially 110 compared them to the interactions of their representatives found on BioGRID (Stark et al., 111 2006). Here, we did not find any differences (p = 1.0), concluding that the static interactions 112 in the Boolean network models are comparable to biological interactions ( Figure 2B).

113
Across all considered Boolean network models, only 3.2% of compounds could be identified 114 as static hubs (z-score >2.5 (Guimerà and Nunes Amaral, 2005)). Comparing these nodes 115 with our implicant sets showed that they differ significantly (p = 2.2•10 -56 ) ( Figure 2C). 116 Based on these results, we analyzed the effect of targeted interventions, i.  To sum up, literature comparison of our simulations of interventions with implicants could 166 confirm our results independent of the cell context. Based on these results, it can be 167 reasoned that even nodes with only a few connections in a network structure can change the 168 phenotype of biological processes. Our new method helps to identify exactly these 169 dynamically important nodes.

171
Implicant screening identifies promising candidates for treating colorectal cancer 172 patients 173 Based on an extensive literature search, we could already show that our method is capable 174 of identifying promising intervention targets. In the next step, we now guide laboratory 175 experiments based on a modelling approach and the analysis of its dynamics. For this 176 purpose, we developed a new Boolean network model (Appendix 2- Table 1)  Our method identified seven implicants in this newly established Boolean network model. 180 Two of them, ERK and TCF/LEF are static network hubs, while the other five (APC, CIP2A, 181 AKT1, RAC1 and GSK3 _deg) are lowly connected components ( Figure 3A). 182 Based on the availability of small molecules for human treatment, the requirement that the 183 target has not been tested in clinical trials for colorectal cancer patients and that no 184 resistances are known for the inhibitor, we selected the most promising intervention 185 candidates out of the list of implicants and studied their effect. Loss of functions of the 186 implicant APC is one of the most frequently occurring alterations in colorectal cancers 187 (Kinzler and Vogelstein, 1996 (Kawasaki et al., 2003). 199 Therefore, we studied the influence of our interventions on these effects.

200
Treatment with either BVD-523 or with TD-52 reduced the proliferative potential of SW480 201 by 2-fold (mean values, Figure 3B) within 24 hours in comparison to untreated or DMSO 202 treated controls without increasing apoptosis ( Figure 3C). By combining both approaches of 203 inhibition, an even stronger mean inhibitory effect of 3.4-fold was achieved ( Figure 3B). In 204 addition, the migratory potential of ERK or CIP2A treated SW480 cells was significantly 205 reduced ( Figures 3D and 3E) and E-cadherin was restored at the cell membrane ( Figure 3F  206 and Appendix 4- Figure 1) thereby supporting that cell adhesion further reduces the 207 migratory potential. Moreover, these effects are further enhanced when the treatment is 208 combined. 209 These initial in vitro experiments, which were performed based on implicants identified in the 210 network structure of a Boolean network model, already show significant results in the 211 inhibition of proliferation and migration of colorectal cancer cellsleading to the conclusion 212 that the consideration of dynamic characteristics has an added value for the identification of 213 promising therapeutic targets. Thus, we presume that dynamic characteristics should be 214 given more significant consideration and should be taken into account when planning 215 interventions. In this study we used Boolean network models. In this paper, we introduce the notion of implicants, defining a set of components that can 319 determine the dynamics of a whole network. Intervention studies with these components 320 support the idea that they are promising candidates for therapeutic approaches. Thus, we 321 believe that dynamic characteristics rather than static hubs are essential factors to be 322 considered when developing drugs or planning laboratory intervention studies. With our 323 approach, we approximate personalized medicine by creating models that are individually 324 adapted to the patient and identify the best possible intervention candidates through 325 simulation. 326 327 328

329
Boolean networks. For the definition of Boolean networks (Kauffman, 1993(Kauffman, , 1969 and Boolean functions = { 1 , … , }, : → (Kauffman, 1993(Kauffman, , 1969. If the value of Fi 335 depends on 1 , 2 , … . , let denote the function defined on these inputs, i.e. ( ) = 336 ( 1 , 2 , … , ). is also called the Boolean function for the i-th position (e.g. gene) of F. 337 To analyze the dynamics of Boolean networks over time, the state of the networks ( ) = 338 ( 1 ( ), … . , ( )) is defined by the states of all variables at point in time t. In synchronous 339 Boolean network models, all Boolean functions are updated at the same time to proceed 340 from state at time t to its successor state at time t + 1 which is defined by ( + 1) = 341 ( 1 ( ( )), … , ( ( ))) of x(t). The dynamics of Boolean network models can be viewed in a 342 state transition graph linking each state of the state space to its successor state. The state-343 space of Boolean networks with n nodes is finite with 2 possible states. Thus, the model will 344 eventually enter a recurrent sequence of states called attractors (cycles) depicting the long-345 term behavior of the network. In a biological context, attractors are associated with 346 phenotypes. 347 All Boolean network simulations were performed with R v3 networks whose dynamics could not be investigated in feasible time, networks with too many 355 attractors and networks in which nearly all nodes are input nodes. Additionally, we excluded 356 non-scale-free networks and those which could be reduced up to their input nodes. In total, 357 we analyzed 35 networks with a size ranging from 5 to 40 nodes (Table 1).

359
Test for scale-freeness. If a Boolean network has a scale-free network architecture, it can 360 be described by the power law distribution ( ) ∝ − where is the power law scaling 361 parameter. To identify scale freeness, we tested if the degree distribution of the network can 362 plausibly be described by the power law distribution by using the R-package poweRlaw 363 v. 70.2 (Gillespie, 2015).

365
Reducing network size. The size and the complexity of Boolean networks increases rapidly 366 with each additional node and exhaustive search methods cannot be applied to larger 367 networks. Therefore, we reduced the search space of large Boolean networks to accelerate 368 the analysis. This was achieved by removing nodes which do not regulate other nodes. The 369 procedure was repeated until all superfluous nodes were removed (Appendix 1-Algorithm 370 3).

372
Identify minimal set of nodes determining the dynamics. Two strategies were applied to 373 determine a minimal set of compounds determining the dynamics of the complete network. 374 The heuristic is defined in terms of significance. Here, the significance of a node is maximal 375 if its transition function depends on its own value. Otherwise the significance of node g is 376 equal to the number of nodes whose transition function depends on g. Therefore, in each 377 iteration the heuristic selects a compound g with the highest significance until a set of 378 implicant nodes is found (Appendix 1-Algorithm 1). For reference we use an exhaustive 379 search algorithm to find minimal implicant sets G of size k for increasing values of k. The 380 algorithm terminates for the smallest value of k such that some subset : = 1 , … , is 381 an implicant set, i.e. for every assignment a network state is observable (Appendix 1-382 Algorithm 2). 383 384 Partial assignments. A partial assignment defines the value of some nodes of the Boolean 385 network. The value of the other nodes is undefined. The transition function 386 = { 1 , … , }, : → can be applied to a partial assignment as follows: If the value of 387 the transition function of a node i is uniquely determined from nodes which are defined 388 under the partial assignment then the node is assigned this value. The value of all other 389 nodes remains unchanged. Repeated application of the transition function to a partial 390 assignment will define the value of not necessarily all nodes. If the assignment can be 391 extended to all nodes and hence defines a state of the network, we say that the network 392 state is observable from the partial assignment.