Search tips
Search criteria 


Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
PLoS ONE. 2007; 2(6): e547.
Published online 2007 June 20. doi:  10.1371/journal.pone.0000547
PMCID: PMC1890306

Selective Control of the Apoptosis Signaling Network in Heterogeneous Cell Populations

Mark Cookson, Academic Editor



Selective control in a population is the ability to control a member of the population while leaving the other members relatively unaffected. The concept of selective control is developed using cell death or apoptosis in heterogeneous cell populations as an example. Control of apoptosis is essential in a variety of therapeutic environments, including cancer where cancer cell death is a desired outcome and Alzheimer's disease where neuron survival is the desired outcome. However, in both cases these responses must occur with minimal response in other cells exposed to treatment; that is, the response must be selective.

Methodology and Principal Findings

Apoptosis signaling in heterogeneous cells is described by an ensemble of gene networks with identical topology but different link strengths. Selective control depends on the statistics of signaling in the ensemble of networks, and we analyze the effects of superposition, non-linearity and feedback on these statistics. Parallel pathways promote normal statistics while series pathways promote skew distributions, which in the most extreme cases become log-normal. We also show that feedback and non-linearity can produce bimodal signaling statistics, as can discreteness and non-linearity. Two methods for optimizing selective control are presented. The first is an exhaustive search method and the second is a linear programming based approach. Though control of a single gene in the signaling network yields little selectivity, control of a few genes typically yields higher levels of selectivity. The statistics of gene combinations susceptible to selective control in heterogeneous apoptosis networks is studied and is used to identify general control strategies.

Conclusions and Significance

We have explored two methods for the study of selectivity in cell populations. The first is an exhaustive search method limited to three node perturbations. The second is an effective linear model, based on interpolation of single node sensitivity, in which the selective combinations can be found by linear programming optimization. We found that selectivity is promoted by acting on the least sensitive nodes in the case of weak populations, while selective control of robust populations is optimized through perturbations of more sensitive nodes. High throughput experiments with heterogeneous cell lines could be designed in an analogous manner, with the further possibility of incorporating the selectivity optimization process into a closed-loop control system.


Living cells carry out their functions, such as working, reproducing and dying, by appropriate response to extracellular and intracellular inputs to a complex network of signaling pathways. Genes which code for the proteins in these pathways are controlled by regulatory proteins which up-regulate or down-regulate these genes, depending on inputs to the signaling network. The enormous effort currently directed at understanding signaling networks may be subdivided into two areas, firstly extracting faithful wiring diagrams for the networks and, secondly developing methods to understand and control the messages which pass through them. In this contribution we develop the concept of selective control in diverse cell populations, and introduce computational methods which optimize selectivity for a particular signal for a designated member of a cell population.

We introduce the concept of selective control or selectivity in cell populations as the requirement of finding a set of inputs which induce one member of the population to produce a desired response while ensuring that the remaining members of the population have a minimal response. In the case of apoptosis or cell death, which we use as an illustrative example, we consider a population of cells and seek methods to kill a selected member of the population while ensuring the survival of the others in the population. Control of apoptosis is essential in a variety of therapeutic environments including cancer where cancer cell death is a desired outcome[1][4] and Alzheimer's disease where neuron survival is the desired outcome[5][7]. However in both cases these responses must occur with minimal response in other cells exposed to treatment that is, the response must be selective.

Though striking progress is occuring in the extraction of networks using a range of experimental data[8], [9], knowledge of signaling networks remains predominantly at the level of topology rather than detailed knowledge of the rate constants and non-linear message passing which occurs in the networks. Models to distinguish between members of a population of cells, for example different cancer cells and different normal tissue types, require differences in message passing parameters and/or expression levels of the genes in the network. Here, the computational procedures for selectivity in cell populations are elucidated using heterogeneous populations, where members of the population are distinguished by having message-passing efficiencies drawn from homogeneous random distributions.

Models of message passing in gene networks range from binary models with discrete message passing rules[10][12] to non-linear ordinary differential equations[13] and to stochastic spatio-temporal models[14] which are simulated using partial differential equations or Monte Carlo methods. Questions of interest also vary greatly, from generic questions about the number of attractors and their stability in random networks[15][19] to modeling the detailed dynamics of gene concentrations in particular pathways[20][22], and to the cellular response such as control of flagellar rotation in bacteria responding to chemotaxis. Some of the tools developed for the analysis of metabolic networks, both dynamically and using steady state flux balance approaches, can be profitably extended to signaling networks[23], [24]. In the flux balance approach (FBA), the ouput of a cell may be optimized with respect to an objective function and subject to the constraints of flux balance at each node in the network.

Though there are conceptual similarities between the FBA and our approach, there are also critical differences. Instead of flux balance, we require message passing rules which describe how a gene responds to the state of its neighbors. To illustrate the effects of different rules, we use an important example from systems biology, the apoptosis network. In particular we discuss the statistics of death signals produced by continuum and discrete message passing rules in this network. We also develop the concept of selective control. Rather than optimizing the objective for one cell or metabolic network, as occurs in FBA, we seek to optimize the response of one cell in a population while minimizing the response of other cells in the population. Selective control is demonstrated using exhaustive search over drug combinations in discrete models and using an approximate linear programming approach. Though drugs affecting specifically every node of the apoptosis network are not yet available, this is a very active field of pharmacological research and it is probably one of the biological networks where this ideal situation, from the control point of view, is closest to reality[25]. Moreover, the network we use is probably still an incomplete representation of the apoptosis network, both for the topology and for the kinetic parameters. Nevertheless, several authors have shown that useful results can be obtained from partially characterized models of biological networks[26], [27].


Statistics of signaling

We have modified the apoptosis network, hsa04210, of the Kegg database to that presented in Fig. 1, 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[20], [22], though the importance of this feedback is not universally accepted. There are 47 genes in Figure 1 and an additional node which we label the output node. The 47 genes may be catagorized (see Table 1) as: input genes (dashed circles in the Fig. 1), 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 Fig. 1) is added to represent the cumulative effect of many genes implicated in the onset of the cell death machinery.

Figure 1
The human apoptosis network modified, as described in the text, from that in the Kegg database (hsa04210).
Table 1
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 = 0,1…,M. Binary networks, where M = 1 have received the most attention, following the work of Kauffman[28]. 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 Figure 1. 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.

Continuous models

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,

equation image

where n(l) 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[29],

equation image

where c/d 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 = 1, c = 1, d = 0 is the majority rule signaling procedure, given by al = sl for sl positive and al = 0 for sl negative.

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 Fig. 1, 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[20], [22] induced by the link between CASP3 and CASP8 in Fig. 1. 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 Fig. 2 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 = si 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 Fig. 2a) is somewhat skew it is not too far from a normal distribution, however the statistics at NFkB (Fig. 2b) is highly skew and is almost lognormal. We now provide a simple explanation for the contrasting statistical behavior occurring for NFkB (34 in Fig. 1) and the output node (48 in Fig. 1).

Figure 2
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 Fig. 2 are presented in Figs. 3 and and4.4. Many paths enter the death node (48) and this is simplified to a set of independent parallel paths in Fig. 3. 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.

Figure 3
A parallel combination of signaling pathways with no series connections.
Figure 4
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 Fig. 1) and a simplified model of this connectivity is illustrated in Figure 4. In this case Eq. (1) yields,

equation image

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,

equation image

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 Fig. 2. One example, where we used the Michaelis-Menten limit of Eq. (2), is presented in Fig. 5. The statistics of the death node (Fig. 5a) 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.

Figure 5
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[20], [22]. This link leads to feedback as illustrated in the subgraph of Fig. 6. To elucidate the effect of feedback on signaling in the apoptosis network, we studied the response of the network in Fig. 6 with random weights on the edges and a range of signal strengths arriving at CASP8 (a1 in Fig. 6).

Figure 6
Feedback loop in the apoptosis network, following

The equations for the signals, si(t+1), arriving at time t+1 at the five genes in Fig. 6 are,

equation image
equation image
equation image
equation image
equation image

In this equation S is the input signal and ωij is the strength of the signaling between nodes i and j. 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(t+1), given by the Hill equation,

equation image

where b,c,d are model parameters. The typical activities of the five genes, ai(t), as a function of the input signal strength, S, are presented in Fig. 7 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[20], 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 Fig. 7c 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.

Figure 7
The steady state activity of genes in the feedback loop illustrated in Fig. 6.

We have also studied the effect of feedback on signaling statistics in heterogeneous cell populations using the full apoptosis network, Fig. 1. Signaling in this network is carried out by using the OSA procedure in combination with full iteration of the loopy sub-graph of Fig. 6. In the linear and Michaelis-Menten cases illustrated in the top two figures of Fig. 7, 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 Fig. 7c), the statistics of the death signal can become bimodal, as illustrated in Fig. 8. The signaling statistics can then be controlled by controlling the feedback and non-linearity in the network.

Figure 8
Death statistics in a heterogeneous population of 50,000 cells with co-operative non-linearity and feedback.

Discrete models

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 = 0,1,…M and |ωij| = 1,…M, 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,

equation image

where n(l) 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/n(l) where n(l) 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/n(l), 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,

equation image

where c1 is a constant which we usually take to be c1 = 1. An external file that holds a picture, illustration, etc.
Object name is pone.0000547.e023.jpg 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,

equation image

where the constant c2 is chosen so that the maximum signal corresponds to ml = M. In our case with M = 10, we take c2 = 2.17. The sigmoidal rule is given by,

equation image

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 Fig. 9a.

Figure 9
Discrete signal death statistics.

The death statistics resulting from these discrete models are presented in Figs. 9b–d. 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 Fig. 9a. In the next section we use these three discrete models in studies of selective control.

Statistics of selective control

Selective control aims at changing the life/death signal of one member of a population with a minimal change of the remaining members of the population. We address in this section the question of how network topology and general signal propagation properties can be used to design strategies for selective control. The control of the life/death signal is realized by acting with external perturbations (drugs) on the nodes and on the signaling flow. As we have seen in the previous sections, the characteristics of the signaling through the network can be strongly dependent on the OSA rule. We will see below how strategies for selectivity are affected by these rules. General strategies for selectivity can be inferred by analyzing the statistics of the nodes involved in selective perturbations. For a given network topology, we can identify the nodes that are more likely to appear in a selective perturbation, and analyze their correlations. Important insight, such as the role of balancing pro-apoptosis and anti-apoptosis perturbations in selectivity, can be obtained from this analysis. Moreover, the correlation analysis revealed that for robust populations selective nodes tend to be the ones that produce the stronger change in the output signal. The opposite is true in the case of weak populations, for which selectivity is improved by acting on nodes that produce weak signal changes in the output.

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 An external file that holds a picture, illustration, etc.
Object name is pone.0000547.e024.jpg and random strength of the links ωij[set membership][−M,−1] for inhibition, and ωij[set membership][1,M] for stimulation (M = 10 in the numerical calculations). The chosen population needs to represent living cells, i.e. it must have the property

equation image

where so 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<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 [lambda with overline] can be selectively controlled if we can find a perturbation on the gene expression levels δmi with the property

equation image
equation image

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−1)36, 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−1)2 perturbation combinations. For an arbitrary k-subset the number of combinations is An external file that holds a picture, illustration, etc.
Object name is pone.0000547.e025.jpg.

In Table 2 we present the selectivity results for ten different populations with 100 individuals each and three different OSA rules. The columns in the Table 2 represent different k-subsets for the three different rules. The value for the threshold so 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 so = 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 so = 1.

Table 2
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 Fig. 10 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.

Figure 10
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. Fig. 11 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 Fig. 11 was obtained using the same threshold (so = 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.

Figure 11
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 Fig. 12 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 Fig. 13 and Fig. 14 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, 11×106, and 26×106 for the linear, sigmoid and logarithmic OSA rules, respectively.

Figure 12
Correlation matrix for nodes appearing in the same selective perturbation in the linear OSA rule.
Figure 13
Correlation matrix for nodes appearing in the same selective perturbation in the sigmoidal OSA rule.
Figure 14
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[30]. 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

equation image

and depends in a nonlinear way on the perturbations of the internal nodes. A typical single node dependence is shown in Fig. 15 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

equation image

where the coefficient cλ,i can be estimated by approximating the single drug dependence such as the ones in Fig. 15 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[31], [32]. Once the coefficient cλ,i and the threshold value so have been fixed, the selectivity problem can be recast in the form of a linear programming optimization problem where we minimize the cost function

equation image

on the polytope defined by the constraint equations

equation image
equation image

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) [30].

Figure 15
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 Fig. 16 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.

Figure 16
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 An external file that holds a picture, illustration, etc.
Object name is pone.0000547.e026.jpg 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 Fig. 16.


Analysis of the statistical behavior of genes in the apoptosis network illustrates the importance of local connectivity on gene activity variations. Genes which receive signals from many parallel paths exhibit normal statistics while genes which lie at the end of a single pathway are prone to large statistical variations and highly skew statistics. Though non-linearity alone does not alter these broad conclusions, the combination of non-linearity and feedback may lead to bimodal statistics. This bimodality occurs due to combination of feedback and bistability[20], which in our signaling models is reflected in a sharp rise in activity with signal strength.

The concept of selective control in heterogeneous cell populations was developed, using the apoptosis network as an illustrative example. Selective control within a heterogeneous cell population is the ability to control one member of the population while leaving the other members relatively unaffected. General selective control strategies that are only dependent on the topology of the network and signaling rules can be inferred from analysis of networks with random link weights. For instance, linear and sigmoid rules identify the same set of nodes that are most efficient in selective control. These nodes can lead to a high degree of selectivity within a given population, especially by balancing pro- and anti- apoptosis perturbations. We have explored two methods for the study of selectivity. The first is an exhaustive search method limited to three node perturbations. The second is an effective linear model, based on interpolation of single node sensitivity, in which the selective combinations can be found by linear programming optimization. The two approaches identify the same strategies for selectivity. We have also identified a general rule that relates the life/death switching robustness of a population to the optimal selectivity strategy. Selectivity is promoted by acting on the least sensitive nodes in the case of weak populations, while selective control of robust populations is optimized through perturbations of more sensitive nodes. More generally selective control is a computational challenge in a broad range of systems biology problems where intervention needs to be directed at subsets of a diverse population.

At a more practical level, high throughput experiments with heterogeneous cell lines could be designed in such a way that the selectivity optimization process is part of a closed-loop control system[33]. Single drug measurements could be used to obtain directly the sensitivity in a given heterogeneous cell population. The linear programming optimization method can then use the sensitivity measurements to identify selective combinations in a model free way. The sensitivity/selectivity optimization can then be improved through iterated experiments in which various computational algorithms and experimental protocols are integrated.


Two programs were used to generate the signaling statistics and to analyse selectivity. One procedure was written in Mathematica and the other in C++. These independent programs were used to check the accuracy of the numerical analysis of the signaling procedures and the signaling statistics.

The software used to run the exhaustive search tests is written in C++ and compiled for both the Linux 32 bit and 64 bit operating systems. We used two clusters, available at the Burnham Institute for Medical Research, Falcon and Bsrc. Both of them are Portable Batch System (PBS) clusters, so that they have a native implementation of queues management enabling submission of the same program multiple times to achieve virtual parallelism. Falcon has 128 CPUs organized in 64 nodes and each of them is a x86 32 bit 2.4 GHz with 1 GB of private memory; Bsrc has 128 CPU organized in 32 nodes and each of them is a x86 64 bit 2.0 GHz with 2 GB of private memory. Exhaustive search requires an exponential computational time in the number of nodes and projecting from the time needed for combinations of 1,2,3 nodes, that is respectively 1 minute, 1.5 hours and 7.5 hours, we estimate that around 90 days would be needed to study all the combinations of 4 drugs. The calculations on the continuous models and the linear programming optimization were implemented using Mathematica on a personal computer. The NMinimize function in Mathematica finds the global minimum when the objective function and constraints are linear.


PMD thanks Leon Glass and Bard Ermentrout for useful discussions.


Competing Interests: The authors have declared that no competing interests exist.

Funding: The work was supported by the Burnham Institute and by Michigan State University. The funders had no role in the work described in the manuscript or in its preparation.


1. Fitzgerald J, Schoeber B, Nielsen UB, Sorger PK. Systems biology and combination therapy in the quest for clinical efficacy. Nature Chemical Biology. 2006;2(9):458–466.
2. Johnstone RW, Ruefli AA, Lowe SW. Apoptosis: A link between cancer genetics and chemotherapy. Cell. 2002;108:153–164. [PubMed]
3. Klein S, Levitzki A. Signal transduction therapy for cancer - whither now. Current Signal Transduction Therapy. 2006;1:1–12.
4. Petak I, Houghton JA, Kopper L. Molecular targeting of cell death signal transduction pathways in cancer. Current Signal Transduction Therapy. 2006;1:113–131.
5. Ankarcrona M, Winblad B. Biomarkers for apoptosis in Alzheimer's disease. International Journal of Geriatric Psychiatry. 2005;20:101–105. [PubMed]
6. Lu KP. Pinning down cell signaling, cancer and Alzheimer's disease. TRENDS in Biochemical Sciences. 2004;29:200–209. [PubMed]
7. Pompl PN, Yemul S, Xiang ZM, Ho L, Haroutunian V, et al. Caspase gene expression in the brain as a function of the clinical progression of Alzheimer's disease. Archives of Neurology. 2003;60:369–376. [PubMed]
8. Barrett CL, Palsson BO. Iterative reconstruction of transcriptional regulatory networks: An algorithmic approach. PLOS computational biology. 2006;2:429–438.
9. Perkins TJ, Jaeger J, Glass L. Reverse engineering the gap gene network of drosophila melanogaster. PLOS computational biology. 2006;2:417–428.
10. Dealy S, Kauffman S, Socolar J. A computational framework for the topological analysis and targeted disruption of signal transduction networks. Complexity. 2005;11:52–60.
11. Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cell-cycle network is robustly designed. Proceedings of the National Academy of Science U S A. 2004;101:4781–4786.
12. Shmulevich I, Kauffman SA, Aldana M. Eukaryotic cells are dynamically ordered or critical but not chaotic. Proceedings of the National Academy of Science U S A. 2005;102:13439–13444.
13. Aldridge BB, Burke JM, Lauffenburger DA, Sorger PK. Physiochemical modeling of cell signaling pathways. Nature Cell Biology. 2006;8:1195–1203.
14. Kholodenko BN. Cell-signalling dynamics in time and space. Nature Reviews Molecular Cell Biology. 2006;7:165–176.
15. Colizza V, Flammini A, Maritan A, Vespignani A. Characterization and modeling of protein-protein interaction networks. Physica A. 2005;352:1–27.
16. Drossel B. Number of attractors in random boolean networks. Physical Review E. 2005;72:16110.
17. Huang S, Eichler G, Bar-Yam Y, Ingber DE. Cell fates as high-dimensional attractor states of a complex gene regulatory network. Physical Review Letters. 2005;94:128701. [PubMed]
18. Klemm K, Bornholdt S. A computational framework for the topological analysis and targeted disruption of signal transduction networks. Physical Review E. 2005;72:55101.
19. Zhou H, Lipowsky R. Dynamic pattern evolution on scale-free networks. Proceedings of the National Academy of Science USA. 2005;102:10052–10057.
20. Bagci EZ, Vodovotz Y, Billiar TR, Ermentrout GB, Bahar I. Bistability in apoptosis: Roles of bax, bcl-2 and mitochondrial permeability transition pores. Biophysical Journal. 2006;91:1546–1549. [PubMed]
21. Buchler NE, Gerland U, Hwa T. Nonlinear protein degradation and the function of genetic circuits. Proceedings of the National Academy of Science USA. 2005;102:9559–9564.
22. Eißing T, Allgöwer F, Bullinger E. Robustness properties of apoptosis models with respect to parameter variations and intrinsic noise. IEEE Proceedings on Systems Biology. 2005;152:221–228.
23. Covert MW, Schilling CH, Palsson B. Regulation of gene expression in flux balance models of metabolism. Journal of Theoretical Biology. 2001;213:73–78. [PubMed]
24. Lee JM, Gianchandani EP, Papin JA. Flux balance analysis in the era of metabolomics. Briefings in Bioinformatics. 2006;7:140–150. [PubMed]
25. Reed JC, Pellecchia M. Apoptosis-based therapies for hematologic malignancies. Blood. 2005;106:408–418. [PubMed]
26. Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in drosophila melanogaster. Journal of Theoretical Biology. 2003;223:1–18. [PubMed]
27. Duarte NC, Becker SA, Jamshidi N, Thiele NI, Mo ML, et al. Global reconstruction of the human metabolic network based on genomic and bibliomic data. Proceedings of the National Academy of Science U S A. 2007;104:1777–1782.
28. Kauffman SA. The origins of order: Self-organization and selection in evolution. Oxford University Press; 1993.
29. Weiss JN. The hill equation revisited: uses and misuses. FASEB Journal. 1997;11:835–841. [PubMed]
30. Cormen TH, Leiserson CE, Rivest RL, Stein C. Introduction to Algorithms. MIT Press and McGraw-Hill; 2001.
31. Feng XJ, Rabitz H. Optimal identification of biochemical reaction networks. Biophysical Journal. 2004;86:1270–1281. [PubMed]
32. Hooshangi S, Feng XJ, Chen D, Li G, Weiss R, Rabitz H. Optimizing genetic circuits by global sensitivity analysis. Biophysical Journal. 2004;87:2195–2202. [PubMed]
33. Ku J, Feng X, Rabitz H. Closed-loop learning control of bio-networks. Journal of Computational Biology. 2004;11:642–659. [PubMed]

Articles from PLoS ONE are provided here courtesy of Public Library of Science