Statistics of signaling
We have modified the apoptosis network, hsa04210, of the Kegg database to that presented in , where network complexes consisting of several genes are split into individual nodes. Recent work emphasizes the importance of positive feedback between CASP3 and CASP8 denoted by the long dashed arrow in the figure
, though the importance of this feedback is not universally accepted. There are 47 genes in and an additional node which we label the output node. The 47 genes may be catagorized (see ) as: input genes (dashed circles in the ), which transmit signals to the network from the other parts of the cellular network; membrane genes which code for membrane proteins and complexes bound to membrane proteins; death genes which signal onset and execution of apoptosis; life genes which reduce apoptotic signals and are upregulated in many cancer cells; and finally the remaining genes in the network. The output, or “death” node (48 in ) is added to represent the cumulative effect of many genes implicated in the onset of the cell death machinery.
The human apoptosis network modified, as described in the text, from that in the Kegg database (hsa04210).
Typical roles of genes in the signaling network.
We consider two types of models, those where the gene activities, ai
, are continuous and those with discrete gene activities, mi
. We use models with continuous activities to illustrate the generic statistics of signal propagation in the apoptosis network, while discrete models are more convenient for the exhaustive search methods used in the selectivity studies. In the discrete models the gene activity has discrete values up to a maximum value M
, so that ml
. Binary networks, where M
1 have received the most attention, following the work of Kauffman
. We considered three discrete cases, M
1,2,10, though here we focus upon M
10 which is closer to the non-linear continuous behavior observed in experiment.
Each gene receives signals from the genes that it is connected to in the signaling network of . The signal arriving at a gene depends on the strength of the connections to its neighbors in the network. We define the strength of these connections to be ωij between the ith and jth genes. Since the network is directed, ωij≠ωji. The values of ωij are poorly characterized even in metabolic networks where they correspond to reaction rates. In the absence of detailed knowledge about these connections we take them to be random variables and in this way develop a generic understanding of signal propagation in heterogeneous cell populations. The link weights, ωij, have positive random values for promotion links and negative random values for inhibition links.
In the continuous models, the edge weights ωij
are uniform continuous random variables and each gene has activity ai
which is a continuous variable. The signal arriving at a gene is given by the sum,
) is the set of genes which send signals to gene l
. The signal sl
arriving at gene l
may be positive or negative, where a negative signal implies inhibition. However, the activity of a gene must be positive or zero so that a negative signal at a gene implies complete inhibition and the gene is switched off, so that its activity is set to zero. This is a basic non-linearity in signaling networks.
In addition, gene activity levels are often observed to depend in a non-linear way on the signals arriving at the gene. A common approximation to the nonlinear response in gene activity is the Hill equation
is the saturation value of the gene activity, d
determines the onset of saturation and the exponent b
is the cooperativity index. The case, b
1, is Michaelis-Menten behavior characteristic of a chemical reaction in the presence of a substrate. The simplest case b
0 is the majority rule signaling procedure, given by al
positive and al
0 for sl
There are several procedures for simulating signal propagation through networks. In binary networks there has been considerable study of synchronous as opposed to asynchronous updates, where in the former case the gene activity levels at time t
are used to update all of the activity levels at time t
+1. In contrast asynchronous methods update gene activity randomly, for example one randomly chosen gene at a time, to model stochastic behavior. The number of attractors found in binary networks appears to depend on the update procedure 16]. In the absence of the feedback link between nodes 42 and 28 in , the network presented there has no loops. In this case signal propagation through the network is deterministic and non-chaotic. Furthermore signal propagation through the loopless network can be carried out in one sweep of the network by ordering the nodes according to their distance from the inputs. The nodes are then updated in order of their distance from the input, a procedure which we denote the optimal signaling algorithm(OSA). As described later, this procedure can be modified to take into account the feedback
induced by the link between CASP3 and CASP8 in . In the absence of this link, the longest path from the inputs to the output has 15 links and hence the OSA algorithm is completed in 15 time steps.
We directly tested the effect of asynchronous, synchronous and OSA procedures on signaling in the loopless apoptosis network and found that they produce essentially the same results. This result is at first counterintuitive as undirected random networks show more complex behavior, such as chaos, and larger sets of attractors. However directed networks are deterministic so that for a given set of inputs in a network with fixed edge weights, and using deterministic message passing rules, there is a unique output. Statistical variations do occur however when the link weights are varied or stochastic noise is added to the message passing rules. Since we are interested in heterogeneous populations, which are analogous to quenched disorder, we consider variations in signaling due to variations in the edge weights.
Typical statistical behavior of signals passing through the loopless apoptosis network are presented in for an important gene, NFkB, in the interior of the network and also for the cumulative death signal at the output. These distributions are found by simulating 50000 different networks, where each network has links (values of ωij
) having weights drawn from a uniform distribution on the interval [0,1]. These simulations were carried out for a model with continuous activities ai
, where Eq. (1) is used to find the total signal arriving at a gene and we use the linear relation ai
for positive signals and ai
0 for negative signals. The input genes in the network were assigned random values on the interval [0,1]. Positive values of the signal arriving at the output node indicate cell death, while negative values denote life. Although the output node statistics (see ) is somewhat skew it is not too far from a normal distribution, however the statistics at NFkB () is highly skew and is almost lognormal. We now provide a simple explanation for the contrasting statistical behavior occurring for NFkB (34 in ) and the output node (48 in ).
The distribution of gene activities in heterogeneous cell populations with a population size of 50,000.
Simplified models elucidating the origin of the signaling statistics observed in are presented in and . Many paths enter the death node (48) and this is simplified to a set of independent parallel paths in . According to Eq. (1), the death signal is then a sum of random variables and it is well known in that case that the statistics of the signal should be a normal distribution, in the asymptotic limit. The observed near normal behavior observed for the death node is then due to the fact that many parallel paths enter the death node. Deviations from the ideal normal distribution are expected for several reasons, including the fact that the activities cannot be negative, the presence of correlations in the signals entering the death node, and due to the fact that we are far from the asymptotic limit.
A parallel combination of signaling pathways with no series connections.
A signaling pathway with four steps in series and with no parallel connections.
In contrast, NFkB is at the end of a chain of connections (see ) and a simplified model of this connectivity is illustrated in . In this case Eq. (1) yields,
This is a random multiplicative process, so that if the link variables ω have random noise, then the output signal, aout
asymptotically obeys log-normal statistics. The log normal distribution, in the variable x
, is given by,
which is typically highly skew and exhibits large fluctuations. Here σ, μ are parameters in the distribution. The statistical variations of signals arriving at genes in complex networks clearly varies a great deal depending on the local connectivity of the genes. In cases where there are mostly linear pathways, as is believed to occur in some cancer cells, there is a greater potential for strong fluctuations in signaling statistics.
Non-linearity is a hallmark of genetic response and we have studied the effect of a variety of non-linear activity-signal behaviors on signaling in the apoptosis network. We tested the effect of various Hill equation parameters on the activity statistics of the apoptosis network. In these calculations, each gene has the same non-linear behavior given by Eq. (2). We found that the generic behavior was similar to that presented in . One example, where we used the Michaelis-Menten limit of Eq. (2), is presented in . The statistics of the death node () remain close to a normal distribution, while the statistics of NFkB remains close to log-normal. The geometry of the network thus controls the signaling statistics even in the presence of non-linearity.
Majority rule signaling with non-linearity.
An important feature absent from the Kegg apoptosis network is feedback. The heavy dashed connection between CASP3 and CASP8 (nodes 42 and 28) produces feedback which has recently been found to be important in the apoptosome
. This link leads to feedback as illustrated in the subgraph of . To elucidate the effect of feedback on signaling in the apoptosis network, we studied the response of the network in with random weights on the edges and a range of signal strengths arriving at CASP8 (a1
Feedback loop in the apoptosis network, following
The equations for the signals, si
+1), arriving at time t
+1 at the five genes in are,
In this equation S
is the input signal and ωij
is the strength of the signaling between nodes i
. In heterogeneous population studies these links are taken to be random. The activity of each gene at time t
+1 is found using a non-linear relation to the signal si
+1), given by the Hill equation,
are model parameters. The typical activities of the five genes, ai
), as a function of the input signal strength, S
, are presented in for three types of Hill equation parameters, linear (top figure), Michaelis-Menten (middle figure) and co-operative (bottom figure). As observed in modeling using ODE's
, co-operative signaling leads to new behavior and a sharp onset of a transition between a low activity state and a high activity state. The behavior of is typical of the co-operative case, and the location of the jump discontinuity and the values of the gene activities depend on the parameter values used in the simulations. One example is presented in this figure. We found that for each parameter set there is a steady state response at long times and this is the value that is plotted in the figures.
The steady state activity of genes in the feedback loop illustrated in .
We have also studied the effect of feedback on signaling statistics in heterogeneous cell populations using the full apoptosis network, . Signaling in this network is carried out by using the OSA procedure in combination with full iteration of the loopy sub-graph of . In the linear and Michaelis-Menten cases illustrated in the top two figures of , feedback amplifies the signal, but does not qualitatively change the statistics of death signals. However in the case of co-operative non-linearity where bi-stability and strong sensitivity to the input signal occurs (see ), the statistics of the death signal can become bimodal, as illustrated in . The signaling statistics can then be controlled by controlling the feedback and non-linearity in the network.
Death statistics in a heterogeneous population of 50,000 cells with co-operative non-linearity and feedback.
Discrete models enable exhaustive search over the activity states of the genes and, as elucidated below, identification of the optimal combinations for selective control. In these models we discretize the weights and the activities into the same number of discrete levels, so that mi
, for a model where the activity level of a gene has a maximum integer value of M
. We find that the generic behaviors for large values of M
are similar to that of the continuous models of the last subsection. On the other hand, for binary networks where each gene has activity zero or one, the behavior is quite different, with a key novelty the fact that many genes have zero incoming signal and hence a decision must be made about their activity in this case. Using the momentum rule, where the state is maintained unless changed by an incoming signal, leads to a strong dependence on initial conditions as the initial state is unchanged unless a signal is received to change it. As in Eq. (1), the signal, sl
arriving at gene l
is a linear superposition given by,
) is the set of genes which signal directly to gene l
. In the discrete models, the signal, sl
, produced by the superposition rule above is then normalised to sl
) where n
) is the number of neighbors of gene l
. We use several different linear and non-linear relations to find the discrete activity of a gene, ml
, from the normalized signal sl
), as described below. In all cases, if the signal arriving at a gene is negative, the gene is completely inhibited and ml
0, which is a basic non-linearity in both continuous and discrete signaling models.
In the linear rule, the discrete activity is found from the normalized signal using,
is a constant which we usually take to be c1
is the ceiling function, which raises a floating point number, x
, to its next largest integer value. For the logarithmic rule the discrete activity is given by,
where the constant c2
is chosen so that the maximum signal corresponds to ml
. In our case with M
10, we take c2
2.17. The sigmoidal rule is given by,
with the parameters α
20 and β
10. This has a form similar to the Fermi function in physics and to dose-response curves in radiation therapy. For M
10 the maximum normalized signal which can arrive at a gene is 100. The functions (12–14) are constructed so that for small signals the gene activity is one, while signals of maximum value yield gene activity 10, ensuring that all possible relations between signal and activity are represented. These behaviors are presented in .
Discrete signal death statistics.
The death statistics resulting from these discrete models are presented in . The discrete linear model leads to unimodal statistics, however the sigmoidal and logarithmic functions lead to bimodal statistics. This is due to the fact that the latter functions amplify small signals, as is evident in . In the next section we use these three discrete models in studies of selective control.
Exhaustive search in discrete models
In this subsection we carry out an exhaustive search of selective perturbations by discretizing the control parameters and signaling variables. In the next section, we will show how some of the key features of selectivity statistics can be captured with a simplified method based on linear programming optimization, which is less demanding from a computational point of view.
We start by generating a population of different apoptosis networks with the same topology, but random values for the initial gene expressions
and random strength of the links ωij
,−1] for inhibition, and ωij
] for stimulation (M
10 in the numerical calculations). The chosen population needs to represent living cells, i.e.
it must have the property
is the signal life/death threshold value, λ is the index labeling individuals in the population, and so,λ
is the output node signal. We take homeostasis into account by adding the constraint that each individual λ remains alive under fluctuations on the input nodes. After so,λ
has been calculated for a given input, we let the input nodes fluctuate and we recalculate its value. If the output of a network is less than a threshold value, that is so,λ
, for ten random fluctuations of the input nodes, then we keep the individual λ in the living population. In the analysis of selectivity, described below, a population of 100 was chosen as the number of different cell types in the human body is of this order. Preliminary studies indicate that, as expected, selective control is easier for smaller populations, provided the number of control nodes is fixed.
Once we have created a living population, we can start to study the effect of external perturbations on the nodes. These perturbations represent the effect of drugs that stimulate or inhibit one or more gene. We will therefore represent the effect of the drug by changing the gene expression levels in the nodes by δmi
. We will say that an individual
can be selectively controlled
if we can find a perturbation on the gene expression levels δmi
with the property
Though there are 48 nodes, 11 are input nodes and one is the output node, so there are 36 control nodes in the network. There are M
−1 possible perturbations on each control node so the total number of possible perturbations on all the nodes, (M
, which is too large to be explored exhaustively. Therefore, we considered perturbations that act only on k
-subsets of all possible perturbations. For k
1 this means that we are considering only single node perturbations, which requires a search over 36(M
−1) possibilities. For k
2 we are considering pairs of nodes, with 630(M
perturbation combinations. For an arbitrary k
-subset the number of combinations is
In we present the selectivity results for ten different populations with 100 individuals each and three different OSA rules. The columns in the represent different k
-subsets for the three different rules. The value for the threshold o
was set to 1. In the linear case, the average percentage of individuals that can be selectively killed within the k
3 subset is about 63%. We have found that this average depends strongly on the value of the threshold, and it decreases for higher threshold values. For instance, by setting the threshold at o
7 the average selectivity (within the k
3 manifold) is reduced to 8.3%. The logarithmic and sigmoid rules give a overall higher selectivity, with an average selectivity of 63.5% and 73.6%, respectively, for death threshold o
Number of individuals that can be selectively controlled in 10 populations with 100 individuals each.
We have studied the statistics of nodes entering in selective combinations. In the linear case, single-node selectivity is obtained by acting on BAD (30), IkBa (33), NFkB (34), BAX (40), CASP3 (42), CASP7 (43) or TP53 (44). These nodes can be divided into those that are pro- and anti-apoptosis, based on their average effect on the output signal. Single drug selective apoptosis is generally induced by stimulating a pro-apoptosis node, or by inhibiting an anti-apoptosis node. For instance, nodes IkBa (33) and NFkB (34) have on average an anti-apoptosis behavior, so they are mainly associated with negative δm. The dotted line in shows the distribution of nodes entering in all the selective combinations found using the linear OSA rule. Notice the peaks corresponding to the same nodes that can induce single-drug selectivity discussed above. In the figure, we also plot the distribution of nodes for the sigmoid (dashed line) and the logarithmic (solid line) rules. The sigmoid and the linear rule identify a very similar set of nodes that are more likely to enter in selective perturbations. Both BAD (30) and CASP3 (42) are strongly present in the linear and sigmoid model. The sigmoid rule slightly enhances the number of combinations in which the inhibition of anti-apoptosis genes represents the mechanism for selectivity, see e.g. CHUCK (32), IkBa (33), NFkB (34), BIRC2 (35) and BCL2 (37). Both linear and sigmoid rules suggest that selectivity might be better achieved by a direct up-regulation of caspases, or by acting on pro- or anti- apoptosis proteins of the BCL2 family. The logarithmic rule results in a different distribution, suggesting a different strategy for selectivity. Most of the selective combinations in the logarithmic case involve nodes upstream in the signaling network, indicating that the best strategy for selectivity is an action on cell-membrane FAS or TNF pathways.
Distribution of nodes entering in selective combinations for the Linear (dotted), Sigmoid (dashed), and Logarithmic (solid) OSA rules.
Interesting correlations can be observed between the statistics of selective and mortal perturbations. Any perturbation that kills at least one member of the population is defined as mortal
. Selective perturbations are mortal perturbations, but there are many more nonselective mortal combinations which kill more than one individual in the population. shows the distribution of nodes entering in selective and mortal combinations in the case of linear (top panel) and logarithmic OSA rules (lower panel). The linear rule shows a strong correlation between nodes that appear in selective combinations and nodes that appear in mortal combinations. The correlation between the mortality and and selectivity distribution is 0.88. The logarithmic case is completely different (lower panel) and exhibits a strong anti-correlation between selectivity and mortality (−0.80). The anti-correlation in the logarithmic case can be reduced by increasing the threshold, and switches to correlation for values of the threshold larger than the average value of the output signal. This behavior can be understood in the following way. For a very high life/death threshold it is very difficult to find individuals that can be killed. In that case we can say that the population is very robust with respect to external perturbations. Therefore, if one individual can be killed, the perturbation that kills that individual is very likely to be selective. We have observed that by pushing the threshold to very high values the correlation between selectivity and mortality approaches one. The nodes that are involved in this case are the ones that are able to produce the strongest change in the output, and are the same nodes usually involved in many mortal combinations. On the other hand, if the population is weak, nodes that are highly mortal are likely to kill more than one individual at the same time, therefore selectivity is associated with nodes that are less mortal, i.e. those that produce small changes in the output signal. The robustness or weakness of a population is determined not only by the value of the life/death threshold, but also by the OSA rules giving different signaling statistics as discussed in the previous section. In fact, the behavior in was obtained using the same threshold (o
1). There, the correlation/anticorrelation is a consequence of the higher sensitivity of the logarithmic rule to external perturbations, which implies that the population is much weaker compared to the linear case.
Distribution of nodes appearing in selective and mortal perturbations in the linear (top panel) and logarithmic (lower panel) OSA rules.
We have analyzed correlations between two nodes appearing in selective combinations. For the linear rule, this is shown in using a correlation matrix plot. Darker dots indicate a higher number of selective combinations containing the two nodes given by the row and column of the matrix. Notice the strong presence of the mitochondrial BAD (30) and BAX (40), and the caspases CASP3 (42) and CASP7(43). However, these nodes often appear in combination with other nodes, many of which have an anti-apoptosis character. The balance of pro and anti-apoptosis perturbations is the key element which increases the selectivity from 1.5% in the k
1 case to the 63% in the three-node perturbation. Notice also that in this linear case the nodes involved in selective combinations are often downstream (i.e. close to the output node) in the signaling network. We show in and the same correlation matrices in the case of the sigmoid and logarithmic OSA rules. In these models more nodes are involved in the selective combinations. The correlation pattern for the sigmoid rule shows strong similarities to the linear case. However, the logarithmic rule results in a qualitatively different pattern in the correlation matrix. Notice for instance that the role of BAD (30), which was dominant in the linear and sigmoid rule, is strongly reduced in this case. In contrast to the linear and sigmoid case, nodes in long pathways dominate in the selective combinations. Overall, the presence of nonlinearity in the OSA rules seems to enhance the possibility of selective control. This is also confirmed by the trend in the total number of k
3 selective combinations, being 9×106
, and 26×106
for the linear, sigmoid and logarithmic OSA rules, respectively.
Correlation matrix for nodes appearing in the same selective perturbation in the linear OSA rule.
Correlation matrix for nodes appearing in the same selective perturbation in the sigmoidal OSA rule.
Correlation matrix for nodes appearing in the same selective perturbation in the logarithmic OSA rule.
Linear programming methods
The exhaustive search method, discussed in the previous subsection, for selective combinations becomes very demanding when combinations with a large number of drugs are involved. A different approach consists in defining an effective model for the dependence of the output signal on the perturbations. The effective model is derived from the original OSA approach using a sensitivity analysis, and some of the selective perturbations found with the effective model are also selective perturbations for the original OSA problem. The advantage of the effective model is that selective combinations can be efficiently obtained by linear programming methods
. We will analyze below the statistics of the selective combinations for the linear OSA method obtained through the effective model. The form of the solutions in many cases involves a number of nodes larger than three, in contrast to the exhaustive method discussed in the previous section. The effective method therefore provides a different sampling of the full space of selective combinations. However, we will see that this different sampling leads to very a similar statistics for the selective nodes.
The output signal from a OSA can be written as
and depends in a nonlinear way on the perturbations of the internal nodes. A typical single node dependence is shown in for a pro-apoptosis BAD (30) and an anti-apoptosis IkBa (33) node for two different individuals (solid and dashed lines) as a function of the strength of the single node perturbation. Notice that the single drug dependence follows a step-like dependence as a consequence of the discreteness of the gene expression levels. Moreover, due to the constraint that the activity must remain positive, there are regions in which the dependence on the external perturbation saturates. The effective linear model can be defined as
where the coefficient cλ,i
can be estimated by approximating the single drug dependence such as the ones in with a linear dependence using interpolation methods. The linear interpolation works well only for some nodes and individuals, since in many cases the dependence is non-monotonic and highly nonlinear. Moreover, notice that we are neglecting higher order many-node effects that are included in the OSA approach. However, we will see below that this does not affect considerably the statistics of the genes that are more likely to appear in selective combinations. The interpolation method also allows us to identify nodes that lead to the highest variations of the output, and the selective combinations can be restricted to nodes within that set. The smallest coefficients cλ,i
can therefore be neglected. The reduction of the control parameter phase space by sensitivity analysis is often a key element in global optimization problems. This was shown explicitly in the case of parameter identification in biochemical reaction networks
. Once the coefficient cλ,i
and the threshold value o
have been fixed, the selectivity problem can be recast in the form of a linear programming optimization problem where we minimize the cost function
on the polytope defined by the constraint equations
The solution provided by the linear programming method is optimal in the sense that it gives the global minimum of the function in Eq. (20) 
Output signal as a function of the perturbation strength on a pro-apoptosis (BAD, red increasing line) and anti-apoptosis node (IkBa, green decreasing line).
We show in the statistics of the nodes found in the linear programming optimization (solid line) compared to the same statistics obtained with the exhaustive search described in the previous section (dotted line for linear, dashed line for sigmoid). The two approaches identify basically the same nodes as the ones that are more likely to be present in selective combinations. Notice the strong presence of BAD (30), BAX (40) and Caspases (42–44) in the three approaches. Also, the relatively strong peaks at IL-3R (17), PI3K (20), and AKT3 (26) are captured by the effective approach. These peaks suggest the possibility of selective control by acting on the AKT signaling pathway.
Distribution of nodes entering in selective combinations obtained with the Linear OSA (dotted), sigmoid (dashed) and with the effective linear programming method (solid line).
We also have used the optimal control parameter
obtained with the linear programming on the original nonlinear OSA and checked the system for selectivity. Typically, we found that the linear programming solution is only partially selective for the original OSA rule. The drug combinations found from the linear programming method can selectively kill only 20% of the 100 individuals in the survivor population, which is considerably smaller than that found using exhaustive search. However, many selective combinations found by the linear programming approach are quasi-selective
, in the sense that they kill two/three individuals in the population rather than one which is requirement for selectivity. This quasi-selectivity captured by the linear programming method is at the origin of the strong similarity in the distributions of .