BACKGROUND
intracellular signaling pathways are important for the development and survival of living organisms, and many human diseases such as autoimmunity, cancer, diabetes and heart disease arise from aberrations in signaling  <cit> - <cit> .
in addition to the importance of intracellular signaling pathways, the last two decades of biomedical research have revealed their complex network-like structure  <cit> .
traditionally, signal transduction has been described as discrete linear pathways connecting cell surface receptors to cellular effectors, and through which information is transmitted and certain cellular responses are generated.
however, individual pathways are highly interconnected with each other, involving extensive cross-talks and feedback loops; thus propagation of cellular signals does not necessarily occur in a linear fashion  <cit> .
this complexity has been a major challenge in cell biology and an obstacle for the development of more effective therapeutic treatment for diseases such as cancer.
to tackle this challenge, systems biology approach has started attracting a great deal of attention through the development of mathematical models of biological systems by combining experimental, theoretical and computational techniques.
in this paper, we use a previously published boolean model of signal transduction in a generic cell to investigate the sensitivity of the network to molecule perturbations  <cit> .
in boolean networks the node activity can be described by two states,  <dig> and  <dig>  active and inactive, and each node is updated based on logical relationships with other nodes.
boolean networks can model a variety of real or artificial networks such as signal transduction networks  <cit> - <cit> , gene regulatory networks  <cit> - <cit> , or neural networks  <cit> .
there are several network/node attributes such as connectivity and types of boolean functions  that are responsible for the dynamics of boolean networks  <cit> , <cit> .
however, each of these individual attributes can be affected by any of the others.
thus, a true understanding of how these attributes combine to produce dynamical effects can only come from studying them in an integrated fashion.
the overall goal of this paper is to use a mathematical approach to fuse these attributes and determine which combination of attributes affects most the dynamics of the network.
so far, research has focused mostly on simplified views involving networks that obey either one single category of boolean functions, or are constructed using a fixed connectivity for all nodes or a widely used connectivity distribution  <cit> , <cit> .
little work has been done on assessing which of the attributes of the nodes has greater impact on the overall dynamics of the network.
we show that a network’s sensitivity to perturbations is the result of combinations of certain attributes of the nodes in various degrees using a mathematical method based on a nonadditive set functions and nonlinear integrals, in particular the so-called choquet integral  <cit> - <cit> .
it is intuitive that both topological and dynamical features have to be considered to understand sensitivity, and that has been illustrated in the literature under various network scenarios, e.g.
<cit>  or  <cit> .
our work adds to that conversation and uses the nonadditive set function approach to go one step further and identify the joint node attributes that are the most important in assessing the sensitivity of the network to molecule perturbations.
in particular we find that out of four topological and three dynamical node attributes under consideration, the out-degree of a node and the average sensitivity of the boolean function governing the dynamics of the node have the most significant joint impact on the overall dynamics of the network and its sensitivity to individual node mutations.
we provide an analysis of this result and how it compares to the roles of the other attributes under consideration.
the paper’s structure is as follows.
in section ‘network, node attributes, and target measure’ we present the important aspects of the network under consideration and introduce the main node attributes used in the nonadditive set function method.
we also introduce the target measure for the impact of node mutations on the dynamics of the network in order to compute the nonadditive set functions of the information fusion method.
more mathematical information on the attributes and the target measure are presented in section ‘methods’, together with the basic mathematical background for the information fusion method and procedure for generating the nonadditive set functions from a given data set.
we end section ‘methods’ with a description of the data sets and the approach used for our numerical investigations.
the method is validated in section ‘results and discussion’ and the best combination of attributes is identified.
an analysis of the results is also included.
we finalize with conclusions and ideas for future research in section ‘conclusions’.
<dig> network, node attributes, and target measure
we start this section by describing the network under consideration, followed by an overview of the node attributes used in information fusion, and the definition of the main target/quantity used to measure the impact of node mutations on the dynamics of the network in the mathematical procedure that is introduced in section ‘information fusion’.
the signal transduction network of a generic fibroblast cell considered in this paper consists of several main signaling pathways, including the receptor tyrosine kinase, the g-protein coupled receptor, and the integrin signaling pathway.
a boolean representation of this network has been provided in  <cit> , and has been studied further in  <cit> .
furthermore, the fully annotated signal transduction model is freely available for simulations and/or download via the cell collective software from www.thecellcollective.org <cit> , <cit> .
each node in the model represents a signaling molecule .
the model also contains nine external input nodes which represent extracellular stimuli, adding a stochastic component to the network, or a “background noise” that exists in all biological systems.
we consider the following node attributes: x1: in-degree .
x2: out-degree .
x3: minimum path length .
x4: average path length .
x5: bias .
x6: average sensitivity of the boolean function, which measures the likelihood that a single flip in an input vector generates a flip of the output of the boolean function .
x7: canalizing degree, which is a measure of the number of ways canalization occurs in the boolean functions .
the attributes considered in this paper represent only a few topological and dynamical node attributes that one could consider.
we have selected them due to being very common in analyses of boolean networks, keeping in mind also that our data-set has a fixed size that does not allow us to consider more than a few attributes at a time .
let us discuss them briefly pointing out relationships between them in relation to the sensitivity of the nodes to mutations and noise in general.
the in-degree is the typical parameter used in basically any boolean network model and indicates the level of connectivity of a node.
the in-degree distribution can distinguish between say random networks and scale-free networks.
a large connectivity may be related to an increased sensitivity to changes in the inputs.
a related attribute is the out-degree as a measure of the spread of influence of a given node.
the more outgoing links a node has, the wider the spread of its influence, so a change in its state may have more impact on the other nodes.
in the fibroblast network we have external stimuli that generate a “background noise” that exists in all biological systems as described in  <cit> .
one could expect that the behavior of nodes with a shorter path to the cell surface be more susceptible to the extracellular noise.
therefore we consider also the minimum and average path length to the cell surface.
a network with large connectivity and/or out-degree is expected to have smaller average path length.
in general the average path length in a network is an indicator of a small-world property.
however, one cannot separate the topological aspects from the dynamical attributes.
for example, if the boolean function of a node with a path length of one is the logical copy function , then the output of this node is the same as the input received from the extracellular stimulus.
therefore, if the extracellular stimulus is random, then the behavior of the node is also random.
however, if a node with a path length of one obeys a canalizing function, or a function with a bias close to  <dig> or  <dig>  then the external noise can be to some extent filtered out.
as a result, the behavior of the node can be rather predictable.
on the other hand, large path lengths may diminish the sensitivity of the network to perturbations.
thus it is important to consider also dynamical attributes such as bias and canalizing degree of the boolean functions, together with the topological attributes.
a large or small bias may induce more order or canalization in the network, thus leading to less sensitivity to disturbances; while a bias close to  <dig>  induces more randomness in the network behavior.
reduced sensitivity can result also from an increased degree of canalization.
the average sensitivity of a boolean function is a well studied complexity measure of how sensitive the output of the function is to changes in the inputs.
it is clearly influenced by the bias, the connectivity, or the degree of canalization.
this is a very natural choice of attribute when we are interested in assessing the impact of noise on the dynamics of the network.
on the other hand, a small connectivity may generate more canalization in the network.
thus the attributes, either topological or dynamical, may be closely related and one cannot completely separate them in any analysis.
further related attributes may be considered in the future, such as clustering coefficients or depth of canalization for fully or partially nested canalizing functions.
however, those are also related to the previously described attributes, so in this study, we take into account only seven basic attributes described in this paragraph and presented further in the next sections.
this way we cover all the most basic attributes, both topological and dynamical, without increasing the number of attributes that may be very closely related, and allowing us to stay within the limits of our data-set.
we create a data set of all these attributes for each node, and fuse this information using nonlinear integrals with respect to nonadditive set functions and establish a target measure that represents the impact of node mutations on the dynamics of the network.
we consider both activating and inactivating mutations, representing constitutively active nodes  and loss-of-function , respectively.
both of these scenarios are relevant to and common in many diseases, including cancer.
one of early examples of an activating mutation includes the protein ras which has been found to be constitutively active due to a mutation in the corresponding gene.
this mutation can result in a continuous stimulation of cell growth and, consequently, in cancer.
the target value we consider in this work is the average hamming distance, ahd, between a non-mutated/wild-type network and a mutated network .
the ahd will be computed under individual node mutations for the mathematical procedure described in what follows and depicted in the scheme of figure  <dig>
we end this section with a note on network updating schemes.
logical networks could be updated according to different update schemes that may generate disturbances or noise in the network.
real networks such as genetical or biological, physical, neural, chemical, or social, are subject to regular disturbances and have the ability to reach functional diversity and aim to maintain the same state under environmental noise .
intrinsic or environmental disturbances can be modeled using asynchrony by updating the value of certain nodes in the network at different time steps .
our goal is to focus on possible mutations within the network  which can be modeled by freezing the value of a node over time, and we choose a synchronous network in which all nodes update their values at each time step, thus controlling the possible environmental disturbances.
<dig> methods
<dig>  attributes
the average sensitivity x <dig> measures the likelihood that a single flip in an input vector generates a flip of the output of the boolean function.
more precisely, let g:{ <dig> }k→{ <dig> } be a boolean function with k inputs, and s the number of elements b∈{ <dig> }k such that h=∑j=1k|aj−bj|= <dig> and g≠g.
the quantity h is the hamming distance between the vectors a and b. the average sensitivity of g is defined as the average of all the quantities s with respect to a probability distribution of the input vectors a, so avs:=<s> <cit> , <cit> .
in  <cit>  it is shown that the average sensitivity is a suitable indicator of the dynamical behavior of a random boolean network, whose value can reflect the regime in which the network operates.
in this paper we consider a uniform probability distribution of the input vectors.
the canalizing degree x <dig> is a measure of canalizing inputs in the functions.
now, a boolean function g:{ <dig> }k→{ <dig> } is canalizing if at least for one value of one of the inputs the output is fixed, irrespective of the values of the other inputs.
let w be the number of ways canalization occurs for a given boolean function with k inputs, so that w∈{ <dig> ,…,k,2k}.
for example, in case of the logical and function all inputs are canalizing with canalizing value  <dig> and canalized output  <dig>  so that w=k.
on the other hand, a constant function has w=2k which is the largest possible value for w. the canalizing degree of a node, x <dig>  is the quotient w/∈ .
canalizing functions may induce stability in certain genetic networks  <cit> .
the nodes of the fibroblast network model  <cit>  obey various boolean rules including canalizing functions with one or more canalizing inputs, or nested/partially nested canalizing functions  <cit> .
<dig>  target value
the target value is the average hamming distance  between a non-muted/wild-type network and a mutated network.
more precisely, let a=, a <dig> …,an)∈{ <dig> }n be the state of the network at time t. then {a,a,…,a} is the trajectory of the network with initial state a up to time t. consider now the trajectory of the same network under the activation of node j, that is {b,b,…,b,… } where b=,…aj− <dig> ,aj+ <dig> …,an), and bj= <dig> for all t. for inactivation we replace  <dig> by  <dig>  then   ahd=1n·t∑i=1n∑t=1tai−bi
<dig>  information fusion
this section contains a brief overview of the mathematical method used for information fusion based on nonadditive set functions and nonlinear integrals.
the process of merging separate information sources using an aggregation tool into a one-dimensional datum on which we can make a prediction or classification decision, is called information fusion.
each information source in a data set is called an attribute and is denoted by xi,i= <dig> ,…,n .
the numerical information obtained from xi is an observation labelled f, and y denotes a numerical evaluation for the target .
to aggregate the numerical information obtained from various attributes, considering the interaction among the contribution rates of the attributes towards the aggregation target, a nonadditive set function defined on the power set of all aggregated attributes and a relevant nonlinear integral, such as the choquet integral, are needed to replace the classical additive measure and the lebesgue integral respectively  <cit> - <cit> .
the amount of combined contributions of various attributes toward the target may not equal the sum of individual contributions, since there may be a synergistic effect of the attributes toward the target.
the nonadditive set functions  describe such interactions among attributes and model their relative importance.
we provide the description of the choquet integral with respect to non-monotone nonadditive set functions used for this study.
let x={x <dig> x <dig> …,xn} be a nonempty finite set of attributes and p be the power set of x. the nonadditive set function μ is a mapping from p to dα+∫0∞μdα  where fα={x| f≥α} for α∈, if both terms on the right-hand side of the formula exist, and not both are infinite.
to incorporate the cases where a set of attributes has a negative contribution toward the target we use a signed nonadditive set function μ:p→.
when μ is a signed nonadditive set function, μ can be expressed as a difference of two nonadditive set functions μ+ and μ−:μ=μ+−μ−.
then the choquet integral with respect to a signed nonadditive set function is defined as   ∫fdμ=∫fdμ+−∫fdμ−.
this formula can be alternatively written as   ∫fdμ=∑j=12n−1zjμj
where μj=μ⋃i:ji=1{xi}=μ+⋃i:ji=1{xi}−μ−⋃i:ji=1{xi} if j is expressed in terms of binary digits jn,jn− <dig> …,j <dig> for every j and   zj=mini:frc∈[1/ <dig> )f−maxi:frc∈[ <dig> /2)f
if it is positive or j=2n− <dig>  and zj= <dig> otherwise.
here frc stands for fractional part.
if fk,k= <dig> ,…,l,j= <dig> ,…,n represent l observations of the n attributes, and yk,k= <dig> ,…,l are the l observations for a chosen target value, then let z=,m= and y= be the corresponding matrices.
then a linear model with respect to the unknown parameters can be created by setting zm=y.
now the optimal values of μ can be calculated by solving for m such that ||zm−y||=∑k=1l∫fkdμ−yk) <dig> is minimal.
the advantage of this alternative formula is that the value of the choquet integral is expressed as a linear function of the values of μ so that the optimal values of μ can be obtained by using an algebraic method.
the details of implementation and effectiveness of the algebraic method can be found in  <cit> .
in order to use this algebraic approach, we generate a data set with all the values of the attributes x <dig> …,x <dig>  this is basically a matrix with n= <dig> lines corresponding to the nodes of the network and n+ <dig> columns corresponding to the n= <dig> attributes and one target value y. the matrix basically contains the information of table  <dig>  where fi is the value of the j-th attribute of the i-th node.
the matrix is split in two parts: the first t lines represent the training set, while the remaining l=n−t lines represent the testing set.
the training data set is used to identify the nonadditive set functions μ. these are used to calculate the value of the choquet integral representing the estimated target values by using the numerical observations of attributes from the testing set.
the estimates are then compared to the actual target values from the testing data set to assess the accuracy of the prediction.
x
1
x
2
x
n
y
t
n
the first t lines correspond to the training set, while the remaining lines correspond to the testing set.
the attribute values are labeled with x, while the target values are labeled with y.
<dig>  data sets and the numerical procedure
in this section we describe the data sets and the numerical procedure, which combined with the mathematical formulae from the previous section, allow us to compute the nonadditive set functions.
to obtain the target values ahd, we generate a wild-type network and a mutated network with a single mutated node.
both networks are iterated  <dig> time steps starting with the same initial states.
after removing a transient part of length  <dig> time steps, we calculate ahd over the remaining  <dig> time steps according to .
this procedure is repeated  <dig> times to account for all individual mutations in the network.
furthermore, we also consider  <dig> different initial states of the network for each of the mutations above.
in our previous analysis of a boolean network model for the fibroblast network considered in this work  <cit> , we noticed that the different types of boolean rules governing the node dynamics lead to similar equilibria when modifying some of the parameters.
these rules include various canalizing functions, threshold functions, and bias functions.
the functions are classified in five large categories based on their impact on the activity level of the network ; this is done by plotting the activity level of the network, also known as the density of ones, at consecutive time points .
however, bifurcation diagrams along connectivity levels or proportions of these classes in the entire network indicate stability with a single fixed point, whose value depends on the modified parameters; the fixed points are reached in less than  <dig> time steps.
in our current simulations we use  <dig> time steps to reach the attractors.
thus, in the context of the signal transduction network under consideration, previous research indicates stability of the dynamics with a frozen state in the long run, regardless of the initial state .
this behavior has been previously observed in simulations of the actual network using cell collective  <cit> .
therefore, for the purpose of this work, we find it reasonable to choose a fairly small number of initial conditions from the entire state space to ease the computational burden.
as specified in  <cit> , nodes belonging to a given class of boolean functions may have a common behavior in the network.
for example, both proteins grb <dig> and nck belong to the same class in that paper, and they are both adaptor proteins .
by applying perturbations to various types of nodes or classes of nodes we may be able to associate biological and mathematical aspects of a particular class of nodes which may not be observed in laboratory.
for more information we refer the reader to  <cit>  and  <cit> .
each data set with  <dig> lines is split into the training set and the testing set.
the size of the training set must be at least 2n− <dig> where n is the number of attributes.
also according to our previous studies  <cit> , the prediction is more accurate when the size of a training set is much larger than the minimum size of a training set.
therefore, we have to combine less than seven attributes at a time.
we choose to use n= <dig> , <dig> attributes at a time to define nonadditive set functions, so that both the training and testing data sets are large enough.
since seven attributes of each node are considered, there are ∑i=357i= <dig> combinations of  <dig>   <dig> and  <dig> attributes.
note that combinations of fewer than three attributes are automatically considered in the mathematical approach even if n≥ <dig>  and repetitions of combinations occur as n is increased since for each n we consider all the possible combinations of  <dig> ,…n attributes.
by varying n we can make sure that the method works well even with less data for testing and training.
thus we test our method with three different training sets whose size increases with larger values of n. by using each of these  <dig> combinations of attributes, we create sub-data-sets  and apply the algebraic method to each of them as portrayed in figure  <dig>  thus we define nonadditive set functions and estimate the values of the choquet integral from  <dig> sub-data-sets using ahd as target.
figure  <dig> gives a simplified example of sub-data-sets created from an original sub-data set for more clarity.
we use the nonadditive set functions estimated from the training set to generate estimated target values for the testing set.
then we compare these estimates to the real target values computed for the testing set from the original data.
if they are close enough  then the method is validated and we can use it to identify the best combination of attributes.
the validation and analysis of the results based on the above procedure unfolds in several steps as described next.
<dig>  average error of target values
the average error  between the estimated target values and the original target values from the testing set is calculated as follows: consider the j-th combination of attributes, where j= <dig> ,…, <dig>  corresponding to the combination of n= <dig>  or  <dig> attributes as indicated in section ‘methods’.
let l be the length of the testing set; in case the training set has length  <dig>  we get l= <dig>  for each l= <dig> ,…,l , denote by yi,jl the l-th target value obtained from the testing set for the sub-data-sets , and ŷi,jl is the corresponding estimated target value obtained with the information fusion procedure.
in the simulations the values of the nonadditive functions are computed by applying a matlab numerical procedure for .
then the quantity   yi,j=1l∑l=1l|yi,jl−ŷi,jl|,i= <dig> ,…, <dig> j= <dig> ,…, <dig>
is the average error of the estimates over the testing set for the sub-data-set .
then we average these quantities over all the initial conditions to obtain the average error of the j-th combination of attributes, namely   ej=1100∑i=1100yi,j,
which we use to determine the accuracy of the method.
the smaller the average error, the better the estimates.
before we present our results, let us review briefly the main properties and hypotheses used for simulations.
the network under consideration of  <dig> nodes is governed by boolean functions of various types , and is updated synchronously.
the network is quenched, therefore the inputs and the boolean rules are fixed throughout the evolution of the network as in the signal transduction model  <cit> .
the information fusion procedure is based on the seven attributes described previously: in-degree, out-degree, minimum path length, average path length, bias, average sensitivity, and canalizing degree.
to assess the optimal non-linear combination of attributes we compute the nonadditive set functions μ, and using them, we determine the impact of individual mutations, both activating and inactivating, by computing the ahd as in .
we use these values in the results section to follow.
<dig>
RESULTS
it is important to validate our approach.
more precisely, once the nonadditive set functions are estimated from the training set , they are used to generate ahd values for the data in the testing set.
these estimated ahd values are then directly compared to the original ahd values of the testing set that are obtained with the actual data.
if the error of the estimation is small, then we accept the estimates as valid.
this validation is presented in section ‘validation of the method’.
once the method is validated, we can find the best combination of attributes by basically identifying the nonadditive set functions that minimize the error.
besides the simple error term used in validation we incorporate also an error term that insures consistency among the initial network states used in simulations.
the results of this approach are presented in section ‘best combination of attributes’.
at the same time we want to make sure the estimates are not dependent on the choice of training and testing sets by applying the method several times with modified training and testing sets.
we modify the length as well as the content of these sets to validate further the proposed method.
we also add an extra refinement to the choice of the best combination of attributes by choosing the  combination of attributes that leads to the highest values of the nonadditive set functions.
this is done in section ‘five-fold cross-validation’.
we finalize in section ‘analysis and discussion of combinations of attributes’ with an analysis of the relation between attributes based on all nonadditive set functions for all the combinations of attributes that have been used in this work, not only the best combination.
<dig>  validation of the method
we define nonadditive set functions by using the training set of size  <dig> out of the available  <dig> since with a bigger training set it is more likely to capture a limited effect of node perturbation.
as specified in section ‘methods’, the choice of n≤ <dig> insures that the algebraic procedure is meaningful even for a data set of  <dig> nodes.
on the other hand, we will use what we call a five-fold cross-validation procedure to show that the estimates for the nonadditive functions are not dependent on the choice of training and testing sets.
once the data are divided into training versus testing sets we find the average error  between the estimated target values and the original target values from the testing set.
the smaller the average error, the better the estimates .
as shown in figure  <dig> for activating mutations, the average error between the original and the estimated target values have small magnitudes.
thus, the estimated nonadditive set functions can capture the interaction among attributes and therefore can be used to model their relative importance.
similar results are obtained for inactivating mutations.
we have also checked the additivity of the estimated set functions.
if x={x <dig> x <dig> …,xn} is a set of n attributes and for a given set function μ we have that ∑i=1nμ≠μ, then the set function is not additive.
our computations indicate that for all the  <dig> set functions this condition is met, so they are indeed nonadditive.
this justifies the use of the choquet integral as an aggregation tool in the information fusion of the node attributes.
thus the fact that the joint contribution of the attributes has more impact on the targets than the stand alone attributes is taken into account.
more precisely, for each combination of attributes we compare the magnitude of the corresponding μ values.
for example, a scenario of only three attributes, say {x <dig> x <dig> x3}, yields μ, μ, μ, μ, μ, μ, μ. we note that in most cases, the maximum μ of all these values corresponds to a combination of attributes, rather than a single attribute.
for example in case of activating mutations, for  <dig>  %  of all combinations of attributes, the maximum μ value corresponds to combinations of at least two attributes.
thus combinations of various attributes of nodes may have an increased influence on the network sensitivity to node perturbation.
we can now apply the procedure to define the best combination of attributes and from that combination find which combination of attributes yields the most significant sensitivity to mutations.
<dig>  best combination of attributes
now we can find the best combination of attributes by basically identifying the nonadditive set functions that minimize the error.
however we need to make sure that our method is also consistent across the initial network states used in simulations before we decide what is the best combination.
we consider a numerical measure that can identify the combination of attributes that has the most impact on the sensitivity of the network to mutations.
aside from minimizing the error in , we want our estimates to be fairly consistent across the  <dig> random initial conditions for which we obtain  <dig> different values for average errors and nonadditive set functions.
therefore, we consider two other quantities besides the average error in : the standard deviation of the errors  and the standard deviation of the estimated nonadditive set functions , with the goal of finding the combination of attributes that minimizes the variation due to the choice of the initial states.
the consistency in target values is given by   ctj=stdyi,ji,j= <dig> ,… <dig>
where yi,j are given in , and stdi stands for standard deviation over the  <dig> different initial states.
to define the consistency in nonadditive set functions we first consider the quantities   si1i2…im=stdμi,m= <dig> ,…n  where the indices 1<=i1<=i2<=⋯<=im<=n identify the collection of attributes for which we compute the nonadditive set function, and n is the number of attributes in the jth combination of attributes.
for a given n there are 2n− <dig> such collections of indices.
then we average the quantities si1i2…im to obtain   cmj=12n−1∑i <dig> i <dig> …,imsi1i2…im
thus, for the consistency in nonadditive set functions, we average the individual standard deviations over all possible values of nonadditive set functions.
we define the best combination of attributes from among the given n attributes to be the one that minimizes the norm of the vector .
thus, the index j∈{ <dig> ,…,91} of the best combination of attributes is given by   j=argminjej2+ctj2+cmj <dig>
a fairly small norm insures that all three individual measures are small enough.
formula  represents our choice of measure for identifying the best combination of attributes in this paper.
alternative measures may be considered, however we have no basis for any speculation on the effect of the chosen measure.
our results reflect the measure defined in .
<dig>  five-fold cross-validation
the results obtained in the previous section are based on finding nonadditive set functions from a training set of size  <dig> out of the available  <dig>  so the testing set is of size  <dig>  in this section we show that the estimates are not dependent on the choice of training and testing sets.
for this purpose we modify the training and testing sets in length as well as in content.
this process is described next.
we also refine the choice of the best combination of attributes by choosing the one that leads to the highest μ values.
the choice is made from the minimal error combinations identified in the previous section.
we randomly partition the  <dig> data sets corresponding to the best combination of attributes into five subsets that are alternatively used as training versus testing sets.
one subset of length  <dig> is used for testing while the remaining subsets of length  <dig> for training.
this process is repeated five times with a different subset for testing.
this method has the advantage that every observation from the data set is used only once for testing and training.
going beyond five repetitions does not generate significant improvements.
we obtain a total of seven/fifteen values μ for activating/inactivating mutations, together with the average error for each of the five testing/training sets.
we calculate the overall mean of these average errors and the average of the μ values over the  <dig> initial conditions.
the overall mean of the average errors is  <dig> %.
although there is a slight increase in the average errors in comparison to those obtained with the training data of size of  <dig>  the errors are still quite small.
we also use different sizes of training and testing sets to investigate their impact on the accuracy of the average error changes and if there is a risk of under-testing when we use the five-fold cross-validation method.
the results are shown in table  <dig> for ahd under activating mutations.
in both cases, the average errors drop as the training data size increases.
however, the drop is very small once a certain threshold is reached.
for example, when the training set represents 50% of the data, the average error is almost as small as for 80% of the data, thus supporting the results obtained by the five-fold cross-validation method.
the impact of increased training set size on the average error for ahd under activating mutations is presented as a percentage.
in figure  <dig>  we graph the values of the nonadditive set functions obtained for the best combination {x <dig> x <dig> x6} for activating mutations, and {x <dig> x <dig> x <dig> x6} for inactivating mutations.
as mentioned earlier, the magnitude of μ indicates the amount of contributions made by an attribute or combination of attributes toward the target.
when the joint contribution of certain attributes is greater than the sum of an individual contribution of those attributes, there is a synergistic effect of those attributes toward the target.
on the other hand, if the joint contribution of attributes is smaller, contributions made by those attributes toward the target are being offset  <cit> .
as shown in figure  <dig>  the differences in the magnitude of μ’s except the greatest μ are relatively small.
note that the joint contribution of x <dig>  the out-degree, and x <dig>  the average sensitivity of the boolean function of a node, is the greatest and significantly larger than all the others.
this holds for either activating or inactivating mutations.
in other words, the out-degree and the average sensitivity of boolean function jointly serve as a good predictor for the impact on the network dynamics of a node mutation.
we find that this result is not only limited to the best combination of attributes {x <dig> x <dig> x6}.
there are  <dig> out of  <dig> combinations of attributes which include both x <dig> and x <dig>  actually,  <dig> % of these  <dig> combinations, include x <dig> and x <dig> in the combination with the greatest μ.
for a description of the whole method see the additional file  <dig> where we provide a detailed algorithm, given step by step.
<dig>  analysis and discussion of combinations of attributes
based on the results obtained from our mathematical method we can now discuss the relation between attributes with the help of the values of nonadditive set functions for all the combinations of attributes that have been used in this work.
we compare these results with some other results in the literature.
we go beyond the attributes that were find as the best combination to explore further possible relationships.
to make this analysis concise, we focus on activating mutations .
in the literature, the in-degree of a node, x <dig>  has been used as the main parameter for analyzing the network evolution.
for example, by using the nk boolean model proposed by kauffman, derrida and pomeau obtained theoretical results showing that the phase transition is controlled by the in-degree and the bias of a node  <cit>  in the context of the entire ensemble of nk networks.
we emphasize again that both topological and dynamical parameters have to be considered for a suitable analysis of boolean networks, even if there may be correlations between them.
our focus is on a numerical assessment  of a particular type of network with fixed links and boolean functions throughout the evolution of the network that represents a sample of a signal transduction network.
thus our statistical assessment focuses not on an entire ensemble of network, but on a very specific one, so our results pertain to the particular network under consideration.
in the future we plan on analyzing other related networks to identify potential commonalities and explore attributes that may be generally important for networks with certain characteristics.
in our simulations, 64% of all combinations of attributes which give the greatest μ values include x <dig>  while only 53% of combinations of attributes include x <dig> .
a mutated node with a large number of downstream nodes has potential for amplifying the effect of the mutation as it spreads to the rest of the network.
however, the out-degree alone may not be the best predictor for the impact of mutations, since the actual boolean function governing the dynamics of the node needs to be taken into account.
this fact is confirmed by our result that x <dig> alone does not yield the greatest μ values.
observe in figure  <dig>  that the joint nonadditive set functions have significantly larger values than the individual attributes.
also notice the negative impact of certain attributes.
going beyond the attributes used in figure  <dig> , there is also a possibility of the bias of a node, x <dig>  being a good indicator for the impact of mutating a node on the network.
in fact, x <dig> alone gives the greatest μ value for 12% of combinations of attributes which include x <dig>  furthermore, x <dig>  and x <dig> , which are bias-related, give the greatest μ value for 8% of combinations of attributes which include x <dig> and x <dig> respectively.
none of the other attributes by themselves gives the greatest μ value.
now let us briefly consider x <dig> .
if a mutated node has multiple canalizing variables that yield a canalized value  <dig>  then the effect of mutation on is greater because this node is supposed to have a lower probability of being active without mutation.
thus, the canalizing degree can indicate the effect of mutating a node to some extent, but it does not distinguish between canalized values  <dig> or  <dig>  which matters if we are interested in activating versus inactivating mutations.
we have shown that the bias and bias related attributes have a potential of being good predictors of the impact of mutating a node.
however, for 88% of combinations of attributes including x <dig> and 92% of those including x <dig> or x <dig>  the attributes alone do not give the greatest μ value.
when assessing the impact of mutating a node, we consider the gap between what is supposed to happen with and without mutation, together with how the effect of mutation spreads.
that can happen either through a few downstream nodes or through a large number of downstream nodes.
that is to say, the combination of attributes x <dig> with x <dig> x <dig> and/or x <dig> should be a good predictor for the impact of mutating a node on the network.
this argument is supported by the result that 57% of combination of attributes which include x <dig> give the greatest μ value when x <dig> is combined with x <dig> x <dig> and/or x <dig>  from this perspective, it is reasonable that the greatest μ value is obtained from the joint combination of out-degree and the average sensitivity of a boolean function.
this work supplements previous research showing that it is important for network analyses in the area of drug target discovery to consider not only the static properties  of individual nodes  of the network , but also properties that give rise to the underlying dynamics .
logical models provide new opportunities to improve the predictive capabilities of computational models as they are easy to construct, and analyses via computer simulations are efficient and capable of covering a relatively large number conditions, and have also been applied to drug target prediction  <cit> , <cit> , as well as prediction of potential drug side-effects and sensitivity  <cit> .
<dig>
CONCLUSIONS
in this paper, we apply a mathematical method of information fusion based on nonadditive set functions and the choquet integral to a boolean model of biochemical signal transduction.
nonadditive set functions are defined for all possible combinations of attributes and can be used to model their relative importance and interactions.
we use an algebraic method to identify the nonadditive set functions and investigate which attribute or combination of attributes is important in predicting the sensitivity of the network to molecule perturbations.
our results support the hypothesis that combinations of different attributes of a node in various degrees play an important role in determining a network’s sensitivity to perturbations, more than the individual attributes.
specifically, we find that the out-degree of a node and the sensitivity of a boolean function to perturbations have the most significant joint impact on the overall dynamics of the network and its sensitivity to perturbations.
the method presented in this work is applicable to any network whose characteristics can be used to generate databases similar to those considered in this paper.
of course, the associated computational burden can restrict the network size that one could consider; however, large networks offer expanded training and testing sets which automatically improve the accuracy of the results.
our existing simulation tools are relatively efficient; for example using cell collective  <cit>  we are currently able to easily simulate our largest model with approximately  <dig> nodes and more than  <dig>  interactions.
that model will be subject for future analysis using the algorithm of this paper.
also, matlab, which was used for our information fusion procedure, is also capable of handling very large matrices.
to our knowledge the work in this paper is the first attempt to assess the importance of nonlinear combinations of attributes in a boolean network.
previous studies on boolean networks mostly assume several parameters as important from the very beginning  and use those to identify dynamics or robustness properties.
on the other hand, the choquet integral has already been used in the area of information fusion for other biological or non-biological systems, so it is a tool that has already shown its usefulness in different contexts.
therefore we believe that the information fusion procedure proposed in this paper has the ability to open the discussion on various aspects of the dynamics and robustness of boolean networks using new non-linear tools from the area of information fusion.
the immediate continuation of this work would be to categorize the magnitude of the most significant attributes that lead to a certain magnitude of the sensitivity measures.
thus one can identify the types of molecules that could be targeted in therapies in order to generate a faster response to a disease treatment, and indicate in which way the attributes should be modified for this purpose.
delving further into the mathematical aspect of the information fusion on boolean networks and the usage of more advanced fusion techniques is one of our future goals.
for instance, a convincing way for the existence of interaction of attributes is to take the mobius transformation for the set function and to see whether its values at sets consisting of more than one attribute are significantly larger than zero.
similarly, one can improve the method by considering a linear or quadratic core of the choquet integral.
on the other hand, we are also considering alternative target measures for the information fusion procedure.
it might be important to look at attributes that capture biologically meaningful, possibly higher-level properties.
for example, in  <cit>  the authors use the mean first passage time from one state of the network to another  state, to identify the best candidate genes for intervention in a gene regulatory network governed by probabilistic boolean functions which are randomly chosen for each node at each time step.
the signal transduction network considered in this paper is governed by boolean rules fixed by the actual biological processes they represent.
however, it would be of interest to consider the alternative target measure of mean passage time to identify the similarities and differences with the ahd approach.
the similarities would point out some possible universal best attributes that are independent of the target measures, while the differences would allow one to adapt to the actual goal of the sensitivity analysis.
this could be identifying the effects of individual mutations, versus say identifying the mutations that lead to a certain dynamical behavior in the shortest amount of time.
at the same time, instead of comparing trajectories as in ahd it might also be useful to focus on reachability of biologically important attractors.
future research could also explore how outcomes will be different if we include attributes not considered in this study, such as feedback loop or clustering information.
we are also interested in trying different methods to preselect attributes to be used for the information fusion procedure instead of trying all combinations of attributes of a node in the network.
moreover, applying the method to an asynchronous network might reveal further aspects of the impact of intrinsic noise in the system in combination with node mutations.
our method is applicable to any type of update procedure.
at the same time, larger data sets of networks may provide more data for information fusion and results could become more accurate, since the training and testing sets would be expanded.
however, large networks have huge state spaces.
therefore the sample of initial conditions needs to be chosen carefully.
in general, it is reasonable to assume that the networks are non-ergodic, so not all states are equally likely at a given time point.
thus, by taking into account the long run activity of the network, it may be useful to consider sampling only the equilibrium state space or distribution, or a biologically meaningful wild-type behavior, thus severely restricting the state space.
a similar approach has been recently used for finding the critical condition for the average sensitivity of a boolean network model governed by nested and partially nested canalizing functions  <cit> , <cit> .
furthermore, exploring other biological networks could reveal some general attribute combinations that may be independent of the data set, thus identifying some underlying characteristics as universal predictors for the sensitivity to perturbations.
we believe that the outcomes of these studies have the potential for providing a better understanding of the underlying mechanism of complex signaling networks as well as identifying possible target molecules and/or attributes for disease treatment.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
nk, th, and jar conceived the biological background and approach, and provided the boolean network data.
nk and zw conceived the information fusion method.
nk and mtm conceived and performed the selection of attributes, the target measure, the validation, and the best attribute combination approach.
nk conducted the analysis and performed most of the simulations.
th and la performed some additional simulations and data analyses.
nk, mtm, and th wrote the manuscript.
mtm, th, jar, zw critically revised the manuscript.
all authors read and approved the final manuscript.
additional file
supplementary material
additional file 1
supplementary material.
click here for file
acknowledgements
this study is partially supported by a 2012- <dig> fund for investing in the research enterprise  grant at the university of nebraska at omaha.
furthermore, we appreciate the financial support from the national institutes of health through grant #5r01da <dig> to jar and th.

 intracellular signaling pathways are important for the development and survival of living organisms, and many human diseases such as autoimmunity, cancer, diabetes and heart disease arise from aberrations in signaling -. in addition to the importance of intracellular signaling pathways, the last two decades of biomedical research have revealed their complex network-like structure. traditionally, signal transduction has been described as discrete linear pathways connecting cell surface receptors to cellular effectors, and through which information is transmitted and certain cellular responses are generated. however, individual pathways are highly interconnected with each other, involving extensive cross-talks and feedback loops; thus propagation of cellular signals does not necessarily occur in a linear fashion. this complexity has been a major challenge in cell biology and an obstacle for the development of more effective therapeutic treatment for diseases such as cancer. to tackle this challenge, systems biology approach has started attracting a great deal of attention through the development of mathematical models of biological systems by combining experimental, theoretical and computational techniques. in this paper, we use a previously published boolean model of signal transduction in a generic cell to investigate the sensitivity of the network to molecule perturbations. in boolean networks the node activity can be described by two states, and active and inactive, and each node is updated based on logical relationships with other nodes. boolean networks can model a variety of real or artificial networks such as signal transduction networks -, gene regulatory networks -, or neural networks. there are several network/node attributes such as connectivity and types of boolean functions that are responsible for the dynamics of boolean networks,. however, each of these individual attributes can be affected by any of the others. thus, a true understanding of how these attributes combine to produce dynamical effects can only come from studying them in an integrated fashion. the overall goal of this paper is to use a mathematical approach to fuse these attributes and determine which combination of attributes affects most the dynamics of the network. so far, research has focused mostly on simplified views involving networks that obey either one single category of boolean functions, or are constructed using a fixed connectivity for all nodes or a widely used connectivity distribution,. little work has been done on assessing which of the attributes of the nodes has greater impact on the overall dynamics of the network. we show that a network’s sensitivity to perturbations is the result of combinations of certain attributes of the nodes in various degrees using a mathematical method based on a nonadditive set functions and nonlinear integrals, in particular the so-called choquet integral -. it is intuitive that both topological and dynamical features have to be considered to understand sensitivity, and that has been illustrated in the literature under various network scenarios, e.g. or. our work adds to that conversation and uses the nonadditive set function approach to go one step further and identify the joint node attributes that are the most important in assessing the sensitivity of the network to molecule perturbations. in particular we find that out of four topological and three dynamical node attributes under consideration, the out-degree of a node and the average sensitivity of the boolean function governing the dynamics of the node have the most significant joint impact on the overall dynamics of the network and its sensitivity to individual node mutations. we provide an analysis of this result and how it compares to the roles of the other attributes under consideration. the paper’s structure is as follows. in section ‘network, node attributes, and target measure’ we present the important aspects of the network under consideration and introduce the main node attributes used in the nonadditive set function method. we also introduce the target measure for the impact of node mutations on the dynamics of the network in order to compute the nonadditive set functions of the information fusion method. more mathematical information on the attributes and the target measure are presented in section ‘methods’, together with the basic mathematical background for the information fusion method and procedure for generating the nonadditive set functions from a given data set. we end section ‘methods’ with a description of the data sets and the approach used for our numerical investigations. the method is validated in section ‘results and discussion’ and the best combination of attributes is identified. an analysis of the results is also included. we finalize with conclusions and ideas for future research in section ‘conclusions’. network, node attributes, and target measure we start this section by describing the network under consideration, followed by an overview of the node attributes used in information fusion, and the definition of the main target/quantity used to measure the impact of node mutations on the dynamics of the network in the mathematical procedure that is introduced in section ‘information fusion’. the signal transduction network of a generic fibroblast cell considered in this paper consists of several main signaling pathways, including the receptor tyrosine kinase, the g-protein coupled receptor, and the integrin signaling pathway. a boolean representation of this network has been provided in, and has been studied further in. furthermore, the fully annotated signal transduction model is freely available for simulations and/or download via the cell collective software from www.thecellcollective.org,. each node in the model represents a signaling molecule. the model also contains nine external input nodes which represent extracellular stimuli, adding a stochastic component to the network, or a “background noise” that exists in all biological systems. we consider the following node attributes: x1: in-degree. x2: out-degree. x3: minimum path length. x4: average path length. x5: bias. x6: average sensitivity of the boolean function, which measures the likelihood that a single flip in an input vector generates a flip of the output of the boolean function. x7: canalizing degree, which is a measure of the number of ways canalization occurs in the boolean functions. the attributes considered in this paper represent only a few topological and dynamical node attributes that one could consider. we have selected them due to being very common in analyses of boolean networks, keeping in mind also that our data-set has a fixed size that does not allow us to consider more than a few attributes at a time. let us discuss them briefly pointing out relationships between them in relation to the sensitivity of the nodes to mutations and noise in general. the in-degree is the typical parameter used in basically any boolean network model and indicates the level of connectivity of a node. the in-degree distribution can distinguish between say random networks and scale-free networks. a large connectivity may be related to an increased sensitivity to changes in the inputs. a related attribute is the out-degree as a measure of the spread of influence of a given node. the more outgoing links a node has, the wider the spread of its influence, so a change in its state may have more impact on the other nodes. in the fibroblast network we have external stimuli that generate a “background noise” that exists in all biological systems as described in. one could expect that the behavior of nodes with a shorter path to the cell surface be more susceptible to the extracellular noise. therefore we consider also the minimum and average path length to the cell surface. a network with large connectivity and/or out-degree is expected to have smaller average path length. in general the average path length in a network is an indicator of a small-world property. however, one cannot separate the topological aspects from the dynamical attributes. for example, if the boolean function of a node with a path length of one is the logical copy function, then the output of this node is the same as the input received from the extracellular stimulus. therefore, if the extracellular stimulus is random, then the behavior of the node is also random. however, if a node with a path length of one obeys a canalizing function, or a function with a bias close to or then the external noise can be to some extent filtered out. as a result, the behavior of the node can be rather predictable. on the other hand, large path lengths may diminish the sensitivity of the network to perturbations. thus it is important to consider also dynamical attributes such as bias and canalizing degree of the boolean functions, together with the topological attributes. a large or small bias may induce more order or canalization in the network, thus leading to less sensitivity to disturbances; while a bias close to induces more randomness in the network behavior. reduced sensitivity can result also from an increased degree of canalization. the average sensitivity of a boolean function is a well studied complexity measure of how sensitive the output of the function is to changes in the inputs. it is clearly influenced by the bias, the connectivity, or the degree of canalization. this is a very natural choice of attribute when we are interested in assessing the impact of noise on the dynamics of the network. on the other hand, a small connectivity may generate more canalization in the network. thus the attributes, either topological or dynamical, may be closely related and one cannot completely separate them in any analysis. further related attributes may be considered in the future, such as clustering coefficients or depth of canalization for fully or partially nested canalizing functions. however, those are also related to the previously described attributes, so in this study, we take into account only seven basic attributes described in this paragraph and presented further in the next sections. this way we cover all the most basic attributes, both topological and dynamical, without increasing the number of attributes that may be very closely related, and allowing us to stay within the limits of our data-set. we create a data set of all these attributes for each node, and fuse this information using nonlinear integrals with respect to nonadditive set functions and establish a target measure that represents the impact of node mutations on the dynamics of the network. we consider both activating and inactivating mutations, representing constitutively active nodes and loss-of-function, respectively. both of these scenarios are relevant to and common in many diseases, including cancer. one of early examples of an activating mutation includes the protein ras which has been found to be constitutively active due to a mutation in the corresponding gene. this mutation can result in a continuous stimulation of cell growth and, consequently, in cancer. the target value we consider in this work is the average hamming distance, ahd, between a non-mutated/wild-type network and a mutated network. the ahd will be computed under individual node mutations for the mathematical procedure described in what follows and depicted in the scheme of figure we end this section with a note on network updating schemes. logical networks could be updated according to different update schemes that may generate disturbances or noise in the network. real networks such as genetical or biological, physical, neural, chemical, or social, are subject to regular disturbances and have the ability to reach functional diversity and aim to maintain the same state under environmental noise. intrinsic or environmental disturbances can be modeled using asynchrony by updating the value of certain nodes in the network at different time steps. our goal is to focus on possible mutations within the network which can be modeled by freezing the value of a node over time, and we choose a synchronous network in which all nodes update their values at each time step, thus controlling the possible environmental disturbances. methods attributes the average sensitivity x measures the likelihood that a single flip in an input vector generates a flip of the output of the boolean function. more precisely, let g:{ }k→{ } be a boolean function with k inputs, and s the number of elements b∈{ }k such that h=∑j=1k|aj−bj|= and g≠g. the quantity h is the hamming distance between the vectors a and b. the average sensitivity of g is defined as the average of all the quantities s with respect to a probability distribution of the input vectors a, so avs:=<s>,. in it is shown that the average sensitivity is a suitable indicator of the dynamical behavior of a random boolean network, whose value can reflect the regime in which the network operates. in this paper we consider a uniform probability distribution of the input vectors. the canalizing degree x is a measure of canalizing inputs in the functions. now, a boolean function g:{ }k→{ } is canalizing if at least for one value of one of the inputs the output is fixed, irrespective of the values of the other inputs. let w be the number of ways canalization occurs for a given boolean function with k inputs, so that w∈{,…,k,2k}. for example, in case of the logical and function all inputs are canalizing with canalizing value and canalized output so that w=k. on the other hand, a constant function has w=2k which is the largest possible value for w. the canalizing degree of a node, x is the quotient w/∈. canalizing functions may induce stability in certain genetic networks. the nodes of the fibroblast network model obey various boolean rules including canalizing functions with one or more canalizing inputs, or nested/partially nested canalizing functions. target value the target value is the average hamming distance between a non-muted/wild-type network and a mutated network. more precisely, let a=, a …,an)∈{ }n be the state of the network at time t. then {a,a,…,a} is the trajectory of the network with initial state a up to time t. consider now the trajectory of the same network under the activation of node j, that is {b,b,…,b,… } where b=,…aj−,aj+ …,an), and bj= for all t. for inactivation we replace by then ahd=1n·t∑i=1n∑t=1tai−bi information fusion this section contains a brief overview of the mathematical method used for information fusion based on nonadditive set functions and nonlinear integrals. the process of merging separate information sources using an aggregation tool into a one-dimensional datum on which we can make a prediction or classification decision, is called information fusion. each information source in a data set is called an attribute and is denoted by xi,i=,…,n. the numerical information obtained from xi is an observation labelled f, and y denotes a numerical evaluation for the target. to aggregate the numerical information obtained from various attributes, considering the interaction among the contribution rates of the attributes towards the aggregation target, a nonadditive set function defined on the power set of all aggregated attributes and a relevant nonlinear integral, such as the choquet integral, are needed to replace the classical additive measure and the lebesgue integral respectively -. the amount of combined contributions of various attributes toward the target may not equal the sum of individual contributions, since there may be a synergistic effect of the attributes toward the target. the nonadditive set functions describe such interactions among attributes and model their relative importance. we provide the description of the choquet integral with respect to non-monotone nonadditive set functions used for this study. let x={x x …,xn} be a nonempty finite set of attributes and p be the power set of x. the nonadditive set function μ is a mapping from p to dα+∫0∞μdα where fα={x| f≥α} for α∈, if both terms on the right-hand side of the formula exist, and not both are infinite. to incorporate the cases where a set of attributes has a negative contribution toward the target we use a signed nonadditive set function μ:p→. when μ is a signed nonadditive set function, μ can be expressed as a difference of two nonadditive set functions μ+ and μ−:μ=μ+−μ−. then the choquet integral with respect to a signed nonadditive set function is defined as ∫fdμ=∫fdμ+−∫fdμ−. this formula can be alternatively written as ∫fdμ=∑j=12n−1zjμj where μj=μ⋃i:ji=1{xi}=μ+⋃i:ji=1{xi}−μ−⋃i:ji=1{xi} if j is expressed in terms of binary digits jn,jn− …,j for every j and zj=mini:frc∈[1/ )f−maxi:frc∈[ /2)f if it is positive or j=2n− and zj= otherwise. here frc stands for fractional part. if fk,k=,…,l,j=,…,n represent l observations of the n attributes, and yk,k=,…,l are the l observations for a chosen target value, then let z=,m= and y= be the corresponding matrices. then a linear model with respect to the unknown parameters can be created by setting zm=y. now the optimal values of μ can be calculated by solving for m such that ||zm−y||=∑k=1l∫fkdμ−yk) is minimal. the advantage of this alternative formula is that the value of the choquet integral is expressed as a linear function of the values of μ so that the optimal values of μ can be obtained by using an algebraic method. the details of implementation and effectiveness of the algebraic method can be found in. in order to use this algebraic approach, we generate a data set with all the values of the attributes x …,x this is basically a matrix with n= lines corresponding to the nodes of the network and n+ columns corresponding to the n= attributes and one target value y. the matrix basically contains the information of table where fi is the value of the j-th attribute of the i-th node. the matrix is split in two parts: the first t lines represent the training set, while the remaining l=n−t lines represent the testing set. the training data set is used to identify the nonadditive set functions μ. these are used to calculate the value of the choquet integral representing the estimated target values by using the numerical observations of attributes from the testing set. the estimates are then compared to the actual target values from the testing data set to assess the accuracy of the prediction. x 1 x 2 x n y t n the first t lines correspond to the training set, while the remaining lines correspond to the testing set. the attribute values are labeled with x, while the target values are labeled with y. data sets and the numerical procedure in this section we describe the data sets and the numerical procedure, which combined with the mathematical formulae from the previous section, allow us to compute the nonadditive set functions. to obtain the target values ahd, we generate a wild-type network and a mutated network with a single mutated node. both networks are iterated time steps starting with the same initial states. after removing a transient part of length time steps, we calculate ahd over the remaining time steps according to. this procedure is repeated times to account for all individual mutations in the network. furthermore, we also consider different initial states of the network for each of the mutations above. in our previous analysis of a boolean network model for the fibroblast network considered in this work, we noticed that the different types of boolean rules governing the node dynamics lead to similar equilibria when modifying some of the parameters. these rules include various canalizing functions, threshold functions, and bias functions. the functions are classified in five large categories based on their impact on the activity level of the network; this is done by plotting the activity level of the network, also known as the density of ones, at consecutive time points. however, bifurcation diagrams along connectivity levels or proportions of these classes in the entire network indicate stability with a single fixed point, whose value depends on the modified parameters; the fixed points are reached in less than time steps. in our current simulations we use time steps to reach the attractors. thus, in the context of the signal transduction network under consideration, previous research indicates stability of the dynamics with a frozen state in the long run, regardless of the initial state. this behavior has been previously observed in simulations of the actual network using cell collective. therefore, for the purpose of this work, we find it reasonable to choose a fairly small number of initial conditions from the entire state space to ease the computational burden. as specified in, nodes belonging to a given class of boolean functions may have a common behavior in the network. for example, both proteins grb and nck belong to the same class in that paper, and they are both adaptor proteins. by applying perturbations to various types of nodes or classes of nodes we may be able to associate biological and mathematical aspects of a particular class of nodes which may not be observed in laboratory. for more information we refer the reader to and. each data set with lines is split into the training set and the testing set. the size of the training set must be at least 2n− where n is the number of attributes. also according to our previous studies, the prediction is more accurate when the size of a training set is much larger than the minimum size of a training set. therefore, we have to combine less than seven attributes at a time. we choose to use n=, attributes at a time to define nonadditive set functions, so that both the training and testing data sets are large enough. since seven attributes of each node are considered, there are ∑i=357i= combinations of and attributes. note that combinations of fewer than three attributes are automatically considered in the mathematical approach even if n≥ and repetitions of combinations occur as n is increased since for each n we consider all the possible combinations of,…n attributes. by varying n we can make sure that the method works well even with less data for testing and training. thus we test our method with three different training sets whose size increases with larger values of n. by using each of these combinations of attributes, we create sub-data-sets and apply the algebraic method to each of them as portrayed in figure thus we define nonadditive set functions and estimate the values of the choquet integral from sub-data-sets using ahd as target. figure gives a simplified example of sub-data-sets created from an original sub-data set for more clarity. we use the nonadditive set functions estimated from the training set to generate estimated target values for the testing set. then we compare these estimates to the real target values computed for the testing set from the original data. if they are close enough then the method is validated and we can use it to identify the best combination of attributes. the validation and analysis of the results based on the above procedure unfolds in several steps as described next. average error of target values the average error between the estimated target values and the original target values from the testing set is calculated as follows: consider the j-th combination of attributes, where j=,…, corresponding to the combination of n= or attributes as indicated in section ‘methods’. let l be the length of the testing set; in case the training set has length we get l= for each l=,…,l, denote by yi,jl the l-th target value obtained from the testing set for the sub-data-sets, and ŷi,jl is the corresponding estimated target value obtained with the information fusion procedure. in the simulations the values of the nonadditive functions are computed by applying a matlab numerical procedure for. then the quantity yi,j=1l∑l=1l|yi,jl−ŷi,jl|,i=,…, j=,…, is the average error of the estimates over the testing set for the sub-data-set. then we average these quantities over all the initial conditions to obtain the average error of the j-th combination of attributes, namely ej=1100∑i=1100yi,j, which we use to determine the accuracy of the method. the smaller the average error, the better the estimates. before we present our results, let us review briefly the main properties and hypotheses used for simulations. the network under consideration of nodes is governed by boolean functions of various types, and is updated synchronously. the network is quenched, therefore the inputs and the boolean rules are fixed throughout the evolution of the network as in the signal transduction model. the information fusion procedure is based on the seven attributes described previously: in-degree, out-degree, minimum path length, average path length, bias, average sensitivity, and canalizing degree. to assess the optimal non-linear combination of attributes we compute the nonadditive set functions μ, and using them, we determine the impact of individual mutations, both activating and inactivating, by computing the ahd as in. we use these values in the results section to follow.  it is important to validate our approach. more precisely, once the nonadditive set functions are estimated from the training set, they are used to generate ahd values for the data in the testing set. these estimated ahd values are then directly compared to the original ahd values of the testing set that are obtained with the actual data. if the error of the estimation is small, then we accept the estimates as valid. this validation is presented in section ‘validation of the method’. once the method is validated, we can find the best combination of attributes by basically identifying the nonadditive set functions that minimize the error. besides the simple error term used in validation we incorporate also an error term that insures consistency among the initial network states used in simulations. the results of this approach are presented in section ‘best combination of attributes’. at the same time we want to make sure the estimates are not dependent on the choice of training and testing sets by applying the method several times with modified training and testing sets. we modify the length as well as the content of these sets to validate further the proposed method. we also add an extra refinement to the choice of the best combination of attributes by choosing the combination of attributes that leads to the highest values of the nonadditive set functions. this is done in section ‘five-fold cross-validation’. we finalize in section ‘analysis and discussion of combinations of attributes’ with an analysis of the relation between attributes based on all nonadditive set functions for all the combinations of attributes that have been used in this work, not only the best combination. validation of the method we define nonadditive set functions by using the training set of size out of the available since with a bigger training set it is more likely to capture a limited effect of node perturbation. as specified in section ‘methods’, the choice of n≤ insures that the algebraic procedure is meaningful even for a data set of nodes. on the other hand, we will use what we call a five-fold cross-validation procedure to show that the estimates for the nonadditive functions are not dependent on the choice of training and testing sets. once the data are divided into training versus testing sets we find the average error between the estimated target values and the original target values from the testing set. the smaller the average error, the better the estimates. as shown in figure for activating mutations, the average error between the original and the estimated target values have small magnitudes. thus, the estimated nonadditive set functions can capture the interaction among attributes and therefore can be used to model their relative importance. similar results are obtained for inactivating mutations. we have also checked the additivity of the estimated set functions. if x={x x …,xn} is a set of n attributes and for a given set function μ we have that ∑i=1nμ≠μ, then the set function is not additive. our computations indicate that for all the set functions this condition is met, so they are indeed nonadditive. this justifies the use of the choquet integral as an aggregation tool in the information fusion of the node attributes. thus the fact that the joint contribution of the attributes has more impact on the targets than the stand alone attributes is taken into account. more precisely, for each combination of attributes we compare the magnitude of the corresponding μ values. for example, a scenario of only three attributes, say {x x x3}, yields μ, μ, μ, μ, μ, μ, μ. we note that in most cases, the maximum μ of all these values corresponds to a combination of attributes, rather than a single attribute. for example in case of activating mutations, for % of all combinations of attributes, the maximum μ value corresponds to combinations of at least two attributes. thus combinations of various attributes of nodes may have an increased influence on the network sensitivity to node perturbation. we can now apply the procedure to define the best combination of attributes and from that combination find which combination of attributes yields the most significant sensitivity to mutations. best combination of attributes now we can find the best combination of attributes by basically identifying the nonadditive set functions that minimize the error. however we need to make sure that our method is also consistent across the initial network states used in simulations before we decide what is the best combination. we consider a numerical measure that can identify the combination of attributes that has the most impact on the sensitivity of the network to mutations. aside from minimizing the error in, we want our estimates to be fairly consistent across the random initial conditions for which we obtain different values for average errors and nonadditive set functions. therefore, we consider two other quantities besides the average error in: the standard deviation of the errors and the standard deviation of the estimated nonadditive set functions, with the goal of finding the combination of attributes that minimizes the variation due to the choice of the initial states. the consistency in target values is given by ctj=stdyi,ji,j=,… where yi,j are given in, and stdi stands for standard deviation over the different initial states. to define the consistency in nonadditive set functions we first consider the quantities si1i2…im=stdμi,m=,…n where the indices 1<=i1<=i2<=⋯<=im<=n identify the collection of attributes for which we compute the nonadditive set function, and n is the number of attributes in the jth combination of attributes. for a given n there are 2n− such collections of indices. then we average the quantities si1i2…im to obtain cmj=12n−1∑i i …,imsi1i2…im thus, for the consistency in nonadditive set functions, we average the individual standard deviations over all possible values of nonadditive set functions. we define the best combination of attributes from among the given n attributes to be the one that minimizes the norm of the vector. thus, the index j∈{,…,91} of the best combination of attributes is given by j=argminjej2+ctj2+cmj a fairly small norm insures that all three individual measures are small enough. formula represents our choice of measure for identifying the best combination of attributes in this paper. alternative measures may be considered, however we have no basis for any speculation on the effect of the chosen measure. our results reflect the measure defined in. five-fold cross-validation the results obtained in the previous section are based on finding nonadditive set functions from a training set of size out of the available so the testing set is of size in this section we show that the estimates are not dependent on the choice of training and testing sets. for this purpose we modify the training and testing sets in length as well as in content. this process is described next. we also refine the choice of the best combination of attributes by choosing the one that leads to the highest μ values. the choice is made from the minimal error combinations identified in the previous section. we randomly partition the data sets corresponding to the best combination of attributes into five subsets that are alternatively used as training versus testing sets. one subset of length is used for testing while the remaining subsets of length for training. this process is repeated five times with a different subset for testing. this method has the advantage that every observation from the data set is used only once for testing and training. going beyond five repetitions does not generate significant improvements. we obtain a total of seven/fifteen values μ for activating/inactivating mutations, together with the average error for each of the five testing/training sets. we calculate the overall mean of these average errors and the average of the μ values over the initial conditions. the overall mean of the average errors is %. although there is a slight increase in the average errors in comparison to those obtained with the training data of size of the errors are still quite small. we also use different sizes of training and testing sets to investigate their impact on the accuracy of the average error changes and if there is a risk of under-testing when we use the five-fold cross-validation method. the results are shown in table for ahd under activating mutations. in both cases, the average errors drop as the training data size increases. however, the drop is very small once a certain threshold is reached. for example, when the training set represents 50% of the data, the average error is almost as small as for 80% of the data, thus supporting the results obtained by the five-fold cross-validation method. the impact of increased training set size on the average error for ahd under activating mutations is presented as a percentage. in figure we graph the values of the nonadditive set functions obtained for the best combination {x x x6} for activating mutations, and {x x x x6} for inactivating mutations. as mentioned earlier, the magnitude of μ indicates the amount of contributions made by an attribute or combination of attributes toward the target. when the joint contribution of certain attributes is greater than the sum of an individual contribution of those attributes, there is a synergistic effect of those attributes toward the target. on the other hand, if the joint contribution of attributes is smaller, contributions made by those attributes toward the target are being offset. as shown in figure the differences in the magnitude of μ’s except the greatest μ are relatively small. note that the joint contribution of x the out-degree, and x the average sensitivity of the boolean function of a node, is the greatest and significantly larger than all the others. this holds for either activating or inactivating mutations. in other words, the out-degree and the average sensitivity of boolean function jointly serve as a good predictor for the impact on the network dynamics of a node mutation. we find that this result is not only limited to the best combination of attributes {x x x6}. there are out of combinations of attributes which include both x and x actually, % of these combinations, include x and x in the combination with the greatest μ. for a description of the whole method see the additional file where we provide a detailed algorithm, given step by step. analysis and discussion of combinations of attributes based on the results obtained from our mathematical method we can now discuss the relation between attributes with the help of the values of nonadditive set functions for all the combinations of attributes that have been used in this work. we compare these results with some other results in the literature. we go beyond the attributes that were find as the best combination to explore further possible relationships. to make this analysis concise, we focus on activating mutations. in the literature, the in-degree of a node, x has been used as the main parameter for analyzing the network evolution. for example, by using the nk boolean model proposed by kauffman, derrida and pomeau obtained theoretical results showing that the phase transition is controlled by the in-degree and the bias of a node in the context of the entire ensemble of nk networks. we emphasize again that both topological and dynamical parameters have to be considered for a suitable analysis of boolean networks, even if there may be correlations between them. our focus is on a numerical assessment of a particular type of network with fixed links and boolean functions throughout the evolution of the network that represents a sample of a signal transduction network. thus our statistical assessment focuses not on an entire ensemble of network, but on a very specific one, so our results pertain to the particular network under consideration. in the future we plan on analyzing other related networks to identify potential commonalities and explore attributes that may be generally important for networks with certain characteristics. in our simulations, 64% of all combinations of attributes which give the greatest μ values include x while only 53% of combinations of attributes include x. a mutated node with a large number of downstream nodes has potential for amplifying the effect of the mutation as it spreads to the rest of the network. however, the out-degree alone may not be the best predictor for the impact of mutations, since the actual boolean function governing the dynamics of the node needs to be taken into account. this fact is confirmed by our result that x alone does not yield the greatest μ values. observe in figure that the joint nonadditive set functions have significantly larger values than the individual attributes. also notice the negative impact of certain attributes. going beyond the attributes used in figure, there is also a possibility of the bias of a node, x being a good indicator for the impact of mutating a node on the network. in fact, x alone gives the greatest μ value for 12% of combinations of attributes which include x furthermore, x and x, which are bias-related, give the greatest μ value for 8% of combinations of attributes which include x and x respectively. none of the other attributes by themselves gives the greatest μ value. now let us briefly consider x. if a mutated node has multiple canalizing variables that yield a canalized value then the effect of mutation on is greater because this node is supposed to have a lower probability of being active without mutation. thus, the canalizing degree can indicate the effect of mutating a node to some extent, but it does not distinguish between canalized values or which matters if we are interested in activating versus inactivating mutations. we have shown that the bias and bias related attributes have a potential of being good predictors of the impact of mutating a node. however, for 88% of combinations of attributes including x and 92% of those including x or x the attributes alone do not give the greatest μ value. when assessing the impact of mutating a node, we consider the gap between what is supposed to happen with and without mutation, together with how the effect of mutation spreads. that can happen either through a few downstream nodes or through a large number of downstream nodes. that is to say, the combination of attributes x with x x and/or x should be a good predictor for the impact of mutating a node on the network. this argument is supported by the result that 57% of combination of attributes which include x give the greatest μ value when x is combined with x x and/or x from this perspective, it is reasonable that the greatest μ value is obtained from the joint combination of out-degree and the average sensitivity of a boolean function. this work supplements previous research showing that it is important for network analyses in the area of drug target discovery to consider not only the static properties of individual nodes of the network, but also properties that give rise to the underlying dynamics. logical models provide new opportunities to improve the predictive capabilities of computational models as they are easy to construct, and analyses via computer simulations are efficient and capable of covering a relatively large number conditions, and have also been applied to drug target prediction,, as well as prediction of potential drug side-effects and sensitivity.  in this paper, we apply a mathematical method of information fusion based on nonadditive set functions and the choquet integral to a boolean model of biochemical signal transduction. nonadditive set functions are defined for all possible combinations of attributes and can be used to model their relative importance and interactions. we use an algebraic method to identify the nonadditive set functions and investigate which attribute or combination of attributes is important in predicting the sensitivity of the network to molecule perturbations. our results support the hypothesis that combinations of different attributes of a node in various degrees play an important role in determining a network’s sensitivity to perturbations, more than the individual attributes. specifically, we find that the out-degree of a node and the sensitivity of a boolean function to perturbations have the most significant joint impact on the overall dynamics of the network and its sensitivity to perturbations. the method presented in this work is applicable to any network whose characteristics can be used to generate databases similar to those considered in this paper. of course, the associated computational burden can restrict the network size that one could consider; however, large networks offer expanded training and testing sets which automatically improve the accuracy of the results. our existing simulation tools are relatively efficient; for example using cell collective we are currently able to easily simulate our largest model with approximately nodes and more than interactions. that model will be subject for future analysis using the algorithm of this paper. also, matlab, which was used for our information fusion procedure, is also capable of handling very large matrices. to our knowledge the work in this paper is the first attempt to assess the importance of nonlinear combinations of attributes in a boolean network. previous studies on boolean networks mostly assume several parameters as important from the very beginning and use those to identify dynamics or robustness properties. on the other hand, the choquet integral has already been used in the area of information fusion for other biological or non-biological systems, so it is a tool that has already shown its usefulness in different contexts. therefore we believe that the information fusion procedure proposed in this paper has the ability to open the discussion on various aspects of the dynamics and robustness of boolean networks using new non-linear tools from the area of information fusion. the immediate continuation of this work would be to categorize the magnitude of the most significant attributes that lead to a certain magnitude of the sensitivity measures. thus one can identify the types of molecules that could be targeted in therapies in order to generate a faster response to a disease treatment, and indicate in which way the attributes should be modified for this purpose. delving further into the mathematical aspect of the information fusion on boolean networks and the usage of more advanced fusion techniques is one of our future goals. for instance, a convincing way for the existence of interaction of attributes is to take the mobius transformation for the set function and to see whether its values at sets consisting of more than one attribute are significantly larger than zero. similarly, one can improve the method by considering a linear or quadratic core of the choquet integral. on the other hand, we are also considering alternative target measures for the information fusion procedure. it might be important to look at attributes that capture biologically meaningful, possibly higher-level properties. for example, in the authors use the mean first passage time from one state of the network to another state, to identify the best candidate genes for intervention in a gene regulatory network governed by probabilistic boolean functions which are randomly chosen for each node at each time step. the signal transduction network considered in this paper is governed by boolean rules fixed by the actual biological processes they represent. however, it would be of interest to consider the alternative target measure of mean passage time to identify the similarities and differences with the ahd approach. the similarities would point out some possible universal best attributes that are independent of the target measures, while the differences would allow one to adapt to the actual goal of the sensitivity analysis. this could be identifying the effects of individual mutations, versus say identifying the mutations that lead to a certain dynamical behavior in the shortest amount of time. at the same time, instead of comparing trajectories as in ahd it might also be useful to focus on reachability of biologically important attractors. future research could also explore how outcomes will be different if we include attributes not considered in this study, such as feedback loop or clustering information. we are also interested in trying different methods to preselect attributes to be used for the information fusion procedure instead of trying all combinations of attributes of a node in the network. moreover, applying the method to an asynchronous network might reveal further aspects of the impact of intrinsic noise in the system in combination with node mutations. our method is applicable to any type of update procedure. at the same time, larger data sets of networks may provide more data for information fusion and results could become more accurate, since the training and testing sets would be expanded. however, large networks have huge state spaces. therefore the sample of initial conditions needs to be chosen carefully. in general, it is reasonable to assume that the networks are non-ergodic, so not all states are equally likely at a given time point. thus, by taking into account the long run activity of the network, it may be useful to consider sampling only the equilibrium state space or distribution, or a biologically meaningful wild-type behavior, thus severely restricting the state space. a similar approach has been recently used for finding the critical condition for the average sensitivity of a boolean network model governed by nested and partially nested canalizing functions,. furthermore, exploring other biological networks could reveal some general attribute combinations that may be independent of the data set, thus identifying some underlying characteristics as universal predictors for the sensitivity to perturbations. we believe that the outcomes of these studies have the potential for providing a better understanding of the underlying mechanism of complex signaling networks as well as identifying possible target molecules and/or attributes for disease treatment. competing interests the authors declare that they have no competing interests. authors’ contributions nk, th, and jar conceived the biological background and approach, and provided the boolean network data. nk and zw conceived the information fusion method. nk and mtm conceived and performed the selection of attributes, the target measure, the validation, and the best attribute combination approach. nk conducted the analysis and performed most of the simulations. th and la performed some additional simulations and data analyses. nk, mtm, and th wrote the manuscript. mtm, th, jar, zw critically revised the manuscript. all authors read and approved the final manuscript. additional file supplementary material additional file 1 supplementary material. click here for file acknowledgements this study is partially supported by a 2012- fund for investing in the research enterprise grant at the university of nebraska at omaha. furthermore, we appreciate the financial support from the national institutes of health through grant #5r01da to jar and th.
an algebraic method for information fusion based on nonadditive set functions is used to assess the joint contribution of boolean network attributes to the sensitivity of the network to individual node mutations. we find that for a biochemical signal transduction network consisting of several main signaling pathways whose nodes represent signaling molecules, the algebraic method provides a robust classification of attribute contributions. the algebraic method is robust against the choice of initial conditions and partition of data sets in training and testing sets for estimation of the nonadditive set functions of the information fusion procedure. the impact of node mutations is assessed using as target measure the average hamming distance between a non-mutated/wild-type network and a mutated network. the node attributes or characteristics under consideration are: in-degree, out-degree, minimum and average path lengths, bias, average sensitivity of boolean functions, and canalizing degrees. this method indicates that for the biochemical network, the most significant impact is generated mainly by the combined effects of two attributes: out-degree, and average sensitivity of nodes. the results support the idea that both topological and dynamical properties of the nodes need to be under consideration. information fusionnode attributessignal transductionnonadditive set functionschoquet integralsensitivity
 an algebraic method for information fusion based on nonadditive set functions is used to assess the joint contribution of boolean network attributes to the sensitivity of the network to individual node mutations. the node attributes or characteristics under consideration are: in-degree, out-degree, minimum and average path lengths, bias, average sensitivity of boolean functions, and canalizing degrees. the impact of node mutations is assessed using as target measure the average hamming distance between a non-mutated/wild-type network and a mutated network.  we find that for a biochemical signal transduction network consisting of several main signaling pathways whose nodes represent signaling molecules, the algebraic method provides a robust classification of attribute contributions. this method indicates that for the biochemical network, the most significant impact is generated mainly by the combined effects of two attributes: out-degree, and average sensitivity of nodes.  the results support the idea that both topological and dynamical properties of the nodes need to be under consideration. the algebraic method is robust against the choice of initial conditions and partition of data sets in training and testing sets for estimation of the nonadditive set functions of the information fusion procedure. information fusionnode attributessignal transductionnonadditive set functionschoquet integralsensitivity 