Method for identification of sensitive nodes in Boolean models of biological networks.

Biological systems are often represented as Boolean networks and analysed to identify sensitive nodes which on perturbation disproportionately change a predefined output. There exist different kinds of perturbation methods: perturbation of function, perturbation of state and perturbation in update scheme. Nodes may have defects in interpretation of the inputs from other nodes and calculation of the node output. To simulate these defects and systematically assess their effect on the system output, two new function perturbations, referred to as 'not of function' and 'function of not', are introduced. In the former, the inputs are assumed to be correctly interpreted but the output of the update rule is perturbed; and in the latter, each input is perturbed but the correct update rule is applied. These and previously used perturbation methods were applied to two existing Boolean models, namely the human melanogenesis signalling network and the fly segment polarity network. Through mathematical simulations, it was found that these methods successfully identified nodes earlier found to be sensitive using other methods, and were also able to identify sensitive nodes which were previously unreported.

1 Introduction

Biological pathways are often represented as networks, the nodes being the biomolecules and edges being the connections. Dynamic models can explain how abundances of biomolecules change over time due to their interactions. Dynamic modelling approaches can be continuous or discrete. In continuous dynamic modelling, the number of nodes and reactions is limited by sparse data leading to limited identifiability of kinetic parameters [1-3]. Boolean modelling is the simplest type of discrete dynamic modelling with abundances represented by 0 (absent/low) and 1 (present/high). It does not require knowledge about the kinetic details of the interactions. The only information needed is the logic of regulatory interactions such as the activating or inhibitory nature of genetic regulations. In a Boolean network, which is a rule‐based binary network, the interaction between nodes is represented using logic rules. Synchronous or asynchronous updating is used to update node states and hence simulate the system dynamics from a given set of initial node states [4]. Boolean networks with a varying number of nodes from <10 [5] to approximately 100 [6, 7] have been used to investigate biological systems.

It is often of interest to identify sensitive nodes in a regulatory network that when perturbed lead to a significant change in the network output. For instance, in models for signalling or metabolic pathways in pathogens, nodes that disproportionately affect survival are potential drug targets. The same motivation exists for identifying sensitive nodes in cancer cell pathways. The robustness (or otherwise) of a signalling network can be assessed from identification and analysis of all sensitive nodes. To estimate sensitivity, a perturbation is applied to every node or edge, and the effect on a set of node states predefined as the system output is calculated.

Previous studies on robustness of Boolean networks have used perturbation methods that can be classified in three broad classes: state perturbations, function perturbations and update rule perturbations. A vast majority of studies use state perturbation to explain system properties including node sensitivity. Shmulevich et al. [8] explored the effect of random gene state perturbation on entire network, i.e. any gene can flip its value for only one‐time point from 0 to 1 or vice versa with probability p. Lee et al. [9] performed node control analysis (constitutive state perturbation) to identify an effective target to reduce skin pigmentation. In this method, the state value of each internal regulatory node is fixed at either ‘0’ for inhibition or ‘1’ for constitutive activation and then the steady‐state activity of output nodes is measured. Fauré et al. [10] simulated the effect of loss of function and gain of function mutation in mammalian cell cycle by constraining selected node within specific value intervals. Subramanian and Gadgil [11] showed that transient state perturbation in Drosophila Melanogaster segment polarity network leads to an ectopic expression pattern. Saadatpour et al. [12] introduced dynamic perturbation that entails setting the node's status opposite to the existing state (diseased state) and normally updating other nodes.

Function perturbations change the normal truth table for a node or set of nodes. Function perturbations have also been used to estimate sensitivity. Garg et al. [13] assume that one gene (or one function) can have a fault at a given time. At a different time in the same trajectory, another gene (or function) can be faulty. The node faults [stochasticity in nodes (SINs)] are interpreted as a change of the current state at that time; moreover, the function faults [stochasticity in function (SIF)] are interpreted as using a different truth table at that time point. They find that the SIN approach predicts biologically implausible behaviour, whereas the SIF approach predicts more biologically relevant robustness. Qian and Dougherty [14] take into account 1 bit function perturbation which entails flipping the value of a single row in the truth table of a probabilistic Boolean model. Another study by the same authors used a similar approach along with the change in probabilistic parameter, i.e. change in the probability of selecting each constitutive Boolean network in the probabilistic Boolean model and changing the perturbation probabilities [15].

Change in updating scheme as a means of assessing robustness has been used by a few researchers. Chaves et al. [16] considered the effect of a perturbation in synchronous update scheme on the dynamics of the model for the D. melanogaster segment polarity genes. Perturbation in the time scales or using different kinds of updating schemes in combination with knockout strategies or state perturbation is also an effective way to identify sensitive nodes [10, 16–18]. Other studies demonstrate different kinds of perturbations not easily classifiable into these three categories. Structural perturbation strategies have been developed [12, 19] to identify essential nodes in a static network whose disruption can reverse the abnormal state of the signalling network. Here, topological intervention involves ranking of the nodes by the effects of their loss (knockout) on the connectivity between the network's inputs and outputs. There are also many reports studying the effect of function perturbation on an ensemble of Boolean networks but not on a specific Boolean network [20, 21].

Here, we introduce two new function perturbation methods and use them to identify sensitive nodes in two specific networks. Our methods are general and applicable to any individual Boolean model or probabilistic Boolean model. These perturbations were applied to the existing melanogenesis signalling network [9] and D. melanogaster segment polarity network [22]. Mathematical simulations revealed that for melanogenesis network, nodes identified as sensitive by the new function perturbation methods are in agreement with state perturbation. Similarly, for D. melanogaster segment polarity network, results of gene mutation performed by Albert and Othmer [22] and transient state perturbation [11] coincide with the results of function perturbation. The nodes identified by each method individually as sensitive nodes are elements of the union of the set of sensitive nodes identified through constitutive activation and constitutive inhibition perturbations. In addition, the new methods identify new nodes as sensitive. We discuss the experimental support for the sensitivity of the newly identified nodes.

2 Materials and methods

2.1 Simulation of existing Boolean models

2.1.1 Melanogenesis signalling network

The melanogenesis network constructed by Lee et al. [9] contains two main modules – the keratinocyte and the melanocyte. There are a total of 62 nodes and 113 links (80 activating and 33 inhibiting links). Of the 62 nodes, there is one external input node [ultraviolet B (UVB) radiation]. The objective was to identify safe and effective targets in the network for reduction of pigmentation as measured by the state of the output nodes. To this end, constitutive activation and constitutive inhibition of each node was simulated by setting the node state to 1 and 0, respectively. The UV input was varied from 0% (always off) to 100% (always on). A ‘wild‐type (WT)’ profile of the average state of each node at each UV level was obtained in the absence of any perturbation. Sensitivity of each node was estimated by calculating the post‐perturbation change in the profile of the three output nodes: B‐cell chronic lymphocytic leukemia/lymphoma 2 (Bcl2) in keratinocyte (Bcl2K), Bcl2 in melanocyte (Bcl2M) and melanin. Nodes whose constitutive activation or inhibition results in a significant reduction in the melanin node activity without significantly affecting Bcl2 activity are reported as potential targets in this paper.

We look at node‐wise sensitivity for each of the three output nodes identified in this paper, and identify nodes whose constitutive activation/repression has a significant effect on each of the three output node profiles, as quantified by the magnitude of the (negative) correlation coefficient between the WT and perturbed profile or the Euclidean distance between the two profiles. The set of nodes thus identified includes the nodes identified by Lee et al. [9] as depigmentation targets. Next, we apply each of the new function perturbation methods to the network, keeping other simulation parameters constant; and assess the overlap between the nodes identified as sensitive by the new methods, and those identified as sensitive using constitutive activation/inhibition perturbations. Simulations were carried out using the network and rules reported in [9]. The intensities of input node UVBs were set to 0, 25, 50, 75 and 100% through a random (non‐cyclic) input with the corresponding probability of being ON. Average of node state values was calculated for each UV level as an average of the last 100 of 1000 time steps for each of 100 random initial conditions. We verified that the results are robust to change in the simulation parameters.

2.1.2 Segment polarity network

The D. melanogaster segment polarity gene expression is defined and maintained through spatiotemporal interactions between gene products including secreted proteins, receptors and transcription factors expressed by cells in a parasegment. A continuous state model was developed by von Dassow who concluded that the patterning was robust to the choice of reaction rate constants [23]. This idea was taken to its logical limit by Albert and Othmer [22] who developed a Boolean model of the regulatory network, thereby obviating the need for any rate parameter. Their model [22] represents 14 cells spread across four parasegments. The first and last parasegments consists of three cells, whereas the second and third parasegments consist of four cells. Each cell has 15 nodes, of which one (SLP) is treated as an input. A parasegment thus has 56 nodes.

Simulations were carried out using Boolean updating rules, initial conditions and parameters specified by Albert and Othmer [22]. Simulations were carried out till attractor state was attained. The simulations were repeated for perturbations, where individual nodes were subjected to constitutive activation/inhibition. Simulations were also carried out after applying each of the two new perturbations introduced here. The node was identified as sensitive if the system steady state on perturbation was either a qualitatively different attractor (i.e. not a point attractor) or the node states differed by 20% from the WT.

All Boolean model simulations were carried out using MATLAB 2015b.

3 Results

First, we describe the two new function perturbation methods to identify sensitive nodes, followed by application of these methods on two existing Boolean networks.

3.1 Two new function perturbation methods

Each node j in the network is associated with variable x(t) which describes its expression level at time t. In Boolean models with synchronous updating, the future state of node j, denoted by x(t + 1), is defined by a logic rule involving the current states of its regulators (inputs), i.e. x(t + 1) = F[ (t)], where F is a Boolean rule and represents the vector of all node states.

Biological processes occur in an inherently noisy environment. Here, we simulate the effect of two permanent defects in the regulatory network. Defects in nodes due to misinterpretation of one or more input signals or miscalculation of the output even when the inputs are received correctly are captured by the function of not (FoN) and not of function (NoF) perturbations (Fig. 1). Biological regulatory networks have been compared with electrical circuits. Nodes are components of a digital circuit that read inputs and emit an output depending on the input. The NoF perturbation simulates a defective node that reads inputs correctly but gives the incorrect output. For example, consider a node where binding of two components results in activation (Fig. 1). In a malfunctioning node with an NoF perturbation, deactivation of the otherwise active node would result from the presence of the two inputs. The FoN perturbation simulates a situation, where the node logic is functioning properly (activation when there are two non‐zero inputs) but there is an error in reading the inputs such that the presence of either component is misinterpreted as absence. This results in a node that is active only when both inputs are absent. These perturbations are incorporated by flipping the output of the function (NoF) or by flipping all the inputs to the function (FoN). Fig. 1 illustrates the two new function perturbation methods of flipping the output (FoN) and misreading the inputs (NoF).

Fig. 1

Methods of function perturbation with an example

Equations for FoN and NoF, Example model, Both a and b are required to activate c but ‘NoF’ perturbation (c 1) will result in inhibition of c, i.e. not of output. Similarly, with ‘FoN’ perturbation (c 2) only inputs are flipped

The effect of perturbing the nodes on the output is compared with the output when there is no perturbation (WT) by using similarity measures such as correlation coefficient and Euclidean distance. Correlation and distance between the WT and the perturbed network are calculated for the steady‐state values/pattern of output, before and after applying the function perturbation using equations where nature of perturbation k∈{0, 1, 2}, i∈{output nodes} and j∈{all nodes}. R is the correlation coefficient, D is Euclidean distance and AvgD is average distance across all the n levels of UV or in general over the length of the vector . is a vector representing steady‐state activity of output node i for all levels of UVB input when node j is subjected to a perturbation of type k, where k = 0, 1, 2 represents no perturbation (WT), NoF perturbation and FoN perturbation, respectively. denotes steady‐state activity of output node i for WT condition (k = 0) for all levels of inputs. We apply these perturbations to two existing Boolean models. The identified sensitive nodes were then compared with those given in previous studies for constitutive state perturbation.

3.2 Effect of function perturbations in melanogenesis network

Effect of function perturbations in melanogenesis network is calculated by measuring a change in the activity profile of the output nodes melanin, Bcl2M and Bcl2K relative to the respective unperturbed output profiles, for the inputs of 0, 25, 50, 75 and 100% UVB. Change in the activity profile of the output is calculated using the measures described previously. The nodes having either the top five Euclidean distance score or correlation coefficient value <−0.8 were selected as sensitive (Table 1). Sensitive nodes with correlation coefficient <−0.8 are depicted in Fig. 2. We also performed constitutive activation and inhibition perturbations of each node in the network, as described in [9], and checked its effect on the outputs using similar measures. When function perturbation is applied to the sensitive nodes, it resulted in a Euclidean distance score of 126–180 for melanin. This is equivalent to AvgD > = 25. Similarly, for Bcl2M and Bcl2K the AvgD score was in the range of 21–24 and 20–30, respectively. AvgD values for constitutive activation and inhibition falls in the range of 30–36 and 7–9 for melanin, 30–32 and 12–20 for Bcl2M and 22–33 and 14–22 for Bcl2K. The comparison showed that most of the nodes identified as sensitive by function perturbations are in agreement with the results of state perturbation analysis but with a few exceptions (non‐underlined nodes in Table 1). There are literature reports suggesting the importance of these nodes. For instance, an experimental study by Jost et al. [24] showed that inhibition of MAPK/ERK kinase (MEK) enzymatic activity in keratinocyte is associated with down‐regulation of Bcl‐2 expression and increased susceptibility to cell death induction. There is also literature evidence for nodes identified as sensitive by both the new methods as well as constitutive activation/inhibition perturbations. For instance, activation of cAMP response element‐binding protein (CREB) is known to activate the microphthalmia‐associated transcription factor (MITF) promoter that promotes melanogenesis [25].

Table 1

Sensitive nodes in the melanogenesis signalling network

	Constitutive activation	Inhibition	NoF	FoN
melanin	MITFproteinM, PKCM, RasM, ET1K, ETRM, SGM, IL1K	ERKM, AktM, PI3KM, RafM, PDK1M, MEKM, MITFproteinM, bcateninM, CREBM, IL1K, ASK1M, MITFmRNAM, p38M, MKK6M	MITFproteinM, SGM, IL1K, RasM, ET1K, ETRM, PKCM	RasM, ET1K, MITFproteinM, ETRM, PKCM
Bcl2M	AktM, PI3KM, PDK1M, PKCM, ET1K, ETRM, RasM, SGM, IL1K	ASK1M, p38M, MKK6M, AktM, PI3KM, PDK1M, CREBM	SGM, IL1K, AktM, RasM, ET1K, ETRM, PDK1M, PI3KM	RasM, ET1K, ETRM, PKCM, PDK1M, AktM
Bcl2K	ASK1K, MKK6K, p38K, MKK4K, JNKK, p53K, RasK ERKK, RafK, SGK	PDK1K, PI3KK, EGFRK, AktK, ASK1K	MKK6K, p38K, ASK1K, JNKK, MKK4K, RasK, SGK	MKK6K, p38K, JNKK, MEKK, RasK

The table lists nodes which highly influence melanin, Bcl2M and Bcl2K activities when perturbed with constitutive activation, inhibition, NoF and FoN. Shown are the nodes with the top five Euclidean distance score and correlation coefficient value <−0.8 for each method of perturbation. Underlined nodes are identified as sensitive both by at least one method of function perturbation and one of constitutive activation/inhibition.

Fig. 2

Effect of function perturbation (left column – NoF and right column – FoN) of selected nodes on UVB‐induced skin pigmentation and Bcl‐2 expressions. Each data point represents average steady‐state activity of output node at 0, 25, 50, 75 and 100% of UVB. Error bar represents standard deviation. Steady‐state value of WT and Lee et al. WT (values of WT outputs mentioned in [9]) fall under each other's standard deviation

Positive relationship between UVB and melanin synthesis (WT) is seen for k = 0. A negative relationship is observed when function perturbations are applied to certain nodes, implying a large average distance and negative correlation. Effect of perturbing PKCM, ETRM and ET1K on melanin are quantitatively identical, ( Nodes that greatly affect WT activity of Bcl2M on ‘NoF’ perturbations. Lee et al. WT values are not available for Bcl2M, A negative relation between UVB and Bcl2K activations. Perturbation in the growth factor receptor‐bound protein 2 (Grb2) and Son of Sevenless (SOS) complex (SG) and GTPase (Ras) in keratinocytes results in a quantitatively overlapping positive relationship

3.3 Effect of function perturbations in D. melanogaster segment polarity network

To identify sensitive nodes in segment polarity network, Boolean function was perturbed for each node in all the cells in all parasegments. It is observed that, for a few nodes, when the logic function was perturbed, the system tends to approach a cyclic attractor. This condition highly differs from WT pattern, where system approaches a point attractor. Therefore, such nodes were classified as sensitive (Table 2). For others (node perturbations resulting in point attractors), the Euclidean distance was calculated and if the distance was >3.35, then the node was assigned as sensitive. This is equivalent to the condition that expression value for at least 20% nodes (11 out of 56) should be changed when a particular node is perturbed. Constitutive inhibition analysis performed previously has shown that null mutation of selected genes in segment polarity network result in alternate steady‐state patterns such as ‘no segmentation pattern’ and ‘broad stripes pattern’ [22]. Similarly, critical and benign nodes identified by transient perturbation are those in which a perturbation leads to the ‘broad stripes pattern’ or alternate steady state [11]. We were successfully able to identify most of the nodes previously identified as sensitive using the new function perturbation methods. There are also a few exceptions where nodes identified as sensitive previously are not identified as such by FoN or NoF perturbations; and new nodes are identified as sensitive (Table 2). Our analysis indicates that patterning is sensitive to perturbation of Cubitus interruptus. This is in contrast to the results of Albert and Othmer. Interestingly, the experimental literature also seems to be divided, with one report suggesting that there is no requirement for cubitus interruptus (CI) before embryonic stage 11 [26] as well as another suggesting that there is an ‘absolute requirement’ for CI in hedgehog signalling [27]. It seems likely that the node sensitivity changes under different conditions, thus supporting the use of multiple perturbation methods to assess sensitivity.

Table 2

Sensitive nodes in the D. melanogaster segment polarity network

	NoF	FoN	Transient perturbation	Gene mutation (knockout)
nodes with cyclic attractor	en(1), en(2), EN(1), EN(2), wg(2), wg(3), wg(4), WG(3), WG(4), PTC(1), hh(1), hh(2), HH(1), SMO(4)	EN(1), EN(2), wg(3), wg(4), WG(3), WG(4), hh(1), HH(1)	wg(1), wg(3), wg(2), WG(1),WG(3), en(2), en(4),EN(2), EN(4), hh(2), hh(4), HH(2), HH(4), ptc(1), PTC(1), PTC(3), SMO(3)	wg(1), wg(2), wg(3),wg(4), en(1), en(2), en(3), en(4), hh(1), hh(2), hh(3), hh(4),ptc(1), ptc(2), ptc(3), ptc(4)
nodes with point attractor	en(4), EN(4), wg(1), WG(1), PTC(3), CI(1), hh(4), HH(2), HH(4), CIR(1), CIR(4), SMO(3)	en(1), en(4), EN(4), WG(1), PTC(3), CI(1), HH(2), HH(4), SMO(3)

Comparison of nodes identified as sensitive using ‘NoF’, ‘FoN’ with those identified as critical nodes by performing transient perturbation and gene mutation. Numbers in parenthesis indicate the cell to which they belong (1–4) with respect to the parasegment. Underlined nodes are identified as sensitive both by one of the methods of function perturbation and one of transient perturbation/gene mutation.

4 Conclusion

We developed two new methods of dynamic function perturbation, namely FoN and NoF and applied it to two existing Boolean models – melanogenesis signalling network and segment polarity network. To our knowledge the perturbation methods closest to FoN and NoF are the SIF and SIN derived by Garg et al. However, there are critical differences. As implemented, NoF results in a change of state of the output node; moreover, it is equivalent to SIN for a given updating time when that node has a defect. In the NoF approach, there is a ‘permanent’ defect that persists through the simulation, whereas in SIN, different nodes may be defective at different times during a single simulation instance. We also find that predictions using NoF are consistent with other perturbation studies in contrast to the ‘implausible’ results obtained using SIN [13]. SIF perturbations assume that defects only arise in active nodes, and use as an example the unlikeliness of transcription without activation. However, leaky transcription is known to occur, and switch – on defects are possible. We include both switch‐on and switch‐off defects in our FoN perturbation. Conceptually, we regard such a perturbation as a defect in interpreting all the input signals. An interesting follow‐up study would be to consider a defect in interpreting one input of a multi‐input node. However, this would complicate the comparison between nodes with differing number of inputs. We have applied each of the two perturbations on two models of differing size and complexity (the segment polarity network with 13 nodes each in 4 cells and the melanogenesis network with 64 nodes in two cell types). We found that as both methods do not require addition of additional inputs to nodes but just involve changing the existing truth tables, implementation was not difficult. We have also verified that it is possible to use these methods to perturb a larger network (∼100 nodes, results not shown). We have analysed the perturbed networks solely from the perspective of identification of sensitive nodes. The effect of small perturbations on attractor states and trajectories can lead to better insights about the stability of the phenotype [28]. Further analysis of the effects of such perturbations on the attractors and their basins of attraction would be desirable but very challenging for large networks even of the order of magnitude of the melanogenesis network. Larger Boolean networks [6, 7] have been used to analyse signalling pathways and cancer pathways. However, analysis of the attractor states of the large unperturbed networks is itself challenging, especially for networks with asynchronous updating [29]. Although methods [30, 31] for identifying attractors for large networks such as the cancer pathways network have been presented, we have focused on identification of node sensitivity and not carried out an analysis of the state space of the perturbed networks in this paper. Although we have used synchronous updating, the methods are equally applicable to asynchronous updating since they both involve a time‐invariant change to the truth table. The updating order and frequency for a particular asynchronous updating scheme can be applied to the modified truth table corresponding to NoF and FoN perturbations.

To examine the ability of FoN and NoF perturbations to identify sensitive nodes in the network, we compared the results obtained by our methods with those of existing ones. Sensitivity is expected to be a function of the nature of the perturbation applied to the network. However, in the case of melanogenesis, a few nodes were robustly identified as sensitive irrespective of perturbation method applied to them, e.g. melanin activity was found to be sensitive to MITFproteinM; moreover, PDK1M and AktM were found to be critical for maintaining Bcl2M level. FoN perturbation results in identification of MEKK as a sensitive node, important for Bcl2K activity, consistent with a reported experimental result. In the case of segment polarity network, we assigned nodes as ‘sensitive’ if perturbation in them results in variation from WT pattern. We were able to identify few new nodes which are important in maintaining WT steady‐state pattern, e.g. imbalance between cubitus interruptus transcriptional activator (CIA) and cubitus interruptus transcriptional repressor (CIR) in posterior cells of parasegment leads to the mutant state [16]. We were able to identify CI as a sensitive node, which is in agreement with some (but not other) experimental reports. As contradictory reports are likely to indicate that the sensitivity differs depending on the experimental condition tested, such data suggests that multiple perturbation methods may capture differing biological situations, and hence a comprehensive determination of node sensitivity may require different perturbation methods to be applied to the network.

Most of the nodes identified as sensitive by our methods are also identified as such by constitutive activation/inhibition. This suggests that our method is more stringent than constitutive activation/repression; and relaxing the Euclidean distance cut‐off criteria would result in identification of more sensitive nodes (inclusion of non‐underlined nodes in Tables 1 and 2). There is no theoretical result to our knowledge suggesting an optimal perturbation method. Indeed, a variety of stochastic perturbations resulting from intrinsic and external sources are encountered by individual cells and developing organisms. Hence, depending on the question sought to be answered, either a specific perturbation corresponding to a specific experiment (for instance gene knockout) is applied to assess the effect on the network output or a suite of perturbation methods is applied to study node and network robustness. In this paper, we have presented two methods that we believe would be useful additions to this suite of perturbation methods for Boolean networks. These methods result in further support for the nodes previously identified as sensitive by other perturbation methods. However, more interestingly, they also lead to the identification of sensitive nodes not identified as such by existing perturbation methods assessed here. This suggests that these new methods query system dynamics and response in a way differing from existing methods. Hence, these methods are expected to be a useful addition to the set of perturbations used to assess node and network sensitivities.