Abstract
Background
Drug combination therapy is commonly used in clinical practice. Many methods including Bliss independence method have been proposed for drug combination design based on simulations models or experiments. Although Bliss independence method can help to solve the drug combination design problem when there are only a small number of combinations, as the number of combinations increases, it may not be scalable. Exploration of system structure becomes important to reduce the complexity of the design problem.
Results
In this paper, we deduced a mathematical model which can simplify the serial structure and parallel structure of biological pathway for synergy evaluation of drug combinations. We demonstrated in steady state the sign of the synergism assessment factor derivative of the original system can be predicted by the sign of its simplified system. In addition, we analyzed the influence of feedback structure on survival ratio of the serial structure. We provided a sufficient condition under which the combination effect could be maintained. Furthermore, we applied our method to find three synergistic drug combinations on tumor necrosis factor αinduced NFκB pathway and subsequently verified by the cell experiment.
Conclusions
We identified several structural properties underlying the Bliss independence criterion, and developed a systematic simplification framework for drug combiation desgin by combining simulation and system reaction network topology analysis. We hope that this work can provide insights to tackle the challenging problem of assessment of combinational drug therapy effect in a large scale signaling pathway. And hopefully in the future our method could be expanded to more general criteria.
Background
Drug combination therapy is commonly used in clinical practice [13]. For example, herbal remedies in traditional Chinese medicine are believed to have synergism effect [4,5]. How to define the drug synergism has been a longstanding controversy amongst pharmacologists, toxicologists and biologists [6,7]. Among existing methods, under the assumption that two drugs acting by independent mechanisms, Bliss independence model is used to define combined effect of two drugs [8]; given that two similar drugs competitively acting on a target, Loewe additivity model is used to predict the combined effect of two drugs [9]. ChouTalalay further proposed the Combination Index (CI) theorem, serving as a general expression and quantification of drug interaction based on the massaction law in biophysics and biochemistry [1012]. These models are widely used in invitro and insilico experiments for drug targets design and doseresponse relationship analysis to instruct the selection of drugs and design of combination scheme [1315]. Although improvements in the application scope and sensitivity of synergy evaluation techniques allow a greater exploitation of drug combination studies, it is unlikely that experimental techniques will be sufficient to completely screen the vast space of drug combinations in a costeffective and timely manner. Hence, finding a way to delimit this space and obtain a manageable set of synergistic combinations is still an ongoing challenge.
To meet this challenge, we present a new method developed from the original Bliss independence criterion to analyze the relationship between structures and effects for combinational drug targets design from a mathematical aspect. Since Bliss model is relatively simple and still widely used by some researchers recently [1618], we start our studies from this simple model and operationally provide a combination effect assessment index inspired by the Combination Index theorem [1012]. With a foundational property of this index, the structure information could be used to help the analysis of drug combinational effects and design of combiantional drug targets. Under this frame, we study two classic structures (serial and parallel structures) in biological signaling systems and propose simplification rules which are helpful for analyzing drug combination effects on the original system. Furthermore, analysis of the feedback structures, which is also very common in signaling pathways, is conducted as an expansion to an original structures without feedback. The usefulness of all the results is demonstrated by numerical experiments.
As a concrete example, we applied our method to an inflammatory angiogenesisrelated pathway, the tumor necrosis factor α (TNFα)induced NFκB pathway. The comprehensive research of this pathway has accumulated abundant exprimental data. This allows us to construct a TNFαinduced NFκB pathway model. Here, we further extended previous model to endothelial cells to construct a more accurate model for drug efficacy prediction. With this new NFκB model in hand, we simualted the combined effects of three important inhibitors, namely Aldehyde, Geldanamycin and PS1145, in NFκB pathway. The simulation results suggested that three inhibitor combinations yeilded significant synergism and were validated the simulated results by cell experiments.
Methods
Original Bliss Independence Criterion
Bliss independence [8] or fractional product method [19] is the index for calculating the expected doseresponse relationship for drug combination therapy as compared to monotherapy. It assumes that the two inhibitors act via independent mechanisms. Then drug combination can be represented as the union of two probabilistically independent events. And this criterion is identical to the mutually nonexclusive case [20]. The combined effect of two inhibitors (F_{UA}) is computed as the product of individual effects of the two inhibitors, F_{UA1 }and F_{UA2}.
where F_{UA }is the remaining enzyme activity (fractional unaffected).
Based on the definition, F_{UA }is the expected combined effect. If the actual combined effect of the two inhibitors is equal to F_{UA}, it is additive effect case and there is no interaction between the two inhibitors. If the actual combined effect is lower than F_{UA}, it is called antagonism. If the actual combined effect is higher than F_{UA}, it is called synergism which leads many possible favorable outcomes like increasing or maintaining drug efficacy as decreasing dosage and provides fundations to the combination therapy [20].
Survival ratio
We use survival ratio as representation of the effect (fractional unaffected) and define it as the ratio of component concentrations before and after intervention:
where a, b are parameters that could be affected. Often they have relationship with the inhibitor doses. a_{0 }and b_{0 }represent the normal values, which are the values before being inhibited; a and b represent the values after being inhibited. The output of a system is usually defined as the concentration of some components.
Inspired by the Combination Index theorem offered by ChouTalalay [10,11], here we introduce an operational concept, "Synergism Assessment Factor", for addressing the interaction of drug combination. Then the Bliss independence criterion could be rewritten as:
Where r(a, b) is the actual combined effect and r(a_{0}, b)·r(a, b_{0}) is the expect combined effect calculated by the product of individual effects. S denotes Synergism Assessment Factor. Eq. 3 is identical to the fraction product equation of Webb [19] and the mutually nonexclusive case in [21]. Compared with the critical point (CI = 1) of ChouTalalay's Combination Index, we used S = 0 as the critical point to determine whether there is sysnergistic effect. Under the first order mutually nonexclusive case, using Eq.3 will get the same conclusion on combination effect as using Combination Index.
Survival ratios of individual invention and combined invention can be measured through invitro or insilico experiments, so it is convenient to verify whether synergism is generated under specific drug combination with this criterion [22]. Therefore, it is widely used in combination therapy design [1618]. However, it is hard to predict the combined effect of two inhibitors without experiments according to this model itself. In order to predict the proper dose range to generate synergism, we have to gain the doseresponse relationship. The doseresponse relationship could be get through series of invitro experiments costly under different doses. Sometimes the doseresponse relationship is assumed to have some special form like Hill equation to reduce the experiment costs [22]. It is feasible when possible targets number is small. As the targets number increases, it will face the combinatorial explosion to choose targets and proper doses, and the experiments cost also increases. New method to narrow down the possibilities in searching targets and doses generating synergism by experiments needs to be developed.
Extended Bliss Independence Criterion
Here, we extend Bliss independence criterion with the sensitivity information of Synergism Assessment Factor. The system with some special structures could be simplified for synergismgenerating targets with doses based on this criterion. We define:
where, r(a, b) is the survival ratio of the system; a and b are parameters affected by inhibitor (ususally the reaction velocity constants [13]), a_{0 }and b_{0 }are separately the value of unaffected. Here DS denotes Synergism Assessment Factor Derivative. Actually DS is the secondorder patial derivative of S. Our criterion is stated as
This criterion is based on the following observations. The connection between DS and S is:
It is easy to see that if DS < 0 (for the interested drug dose ranges), then S is guaranteed to be smaller than 0. It means that synergism is generated over the parameter ranges (a_{0}, a), (b_{0}, b). Similarly, the condition to generate antagonism is also intuitive.
This extended Bliss independence criterion could be seen as the derivative form of the original Bliss independence criterion.
Fundamental property of synergism assessment factor derivative
The extended Bliss independence criterion introduced in Eq. 5 enables us to find special structures of systems that simplification is possible for synergismgenerating targets with doses. Here we provide a basic property of the synergism assessment factor derivative DS. It is the foundation for our main results.
If the inhibitors individually affect some intermediate processes, then the inhibition on these processes could be taken as directly on the products of these processes. That is to say, x = Φ (a) is the product of the process where parameter a is affected; y = Ψ (b) is the product of the process where parameter b is affected. Then the inhibition effects to the output of system on a and b could be taken as on x and y. The analysis of combination effects on the original system could then be limited on a simplified system with x and y as parameters.
Lemma 1
a and b are system parameters that will be affected by inhibitors. a is of processes that will produce product x (x = Φ (a)), and b is of processes that will produce product y (y = Ψ (b)). Then the synergism assessment factor derivative DS of the original system satisfies
Where , x = Φ (a), y = Ψ (b), x_{0 }= Φ (a_{0}), y_{0 }= Ψ (b_{0}). Details of proof are given in Additional file 1.
Additional file 1. Proof of system simplification rules. Proof of Lemma 1 (fundamental property of synergism assessment factor derivative), Corollary 1 (simplification rule for serial structure) and 2 (simplification rule for parallel structure).
Format: DOC Size: 293KB Download file
This file can be viewed with: Microsoft Word Viewer
DS is the synergism assessment factor derivative of the original system while DS' could be seen as that of a simplified system. Usually the signs of derivatives and are easy to know, then the sign of DS could be decided on the sign of DS'; and we only need to analyze the sign of DS', which is the combination effect of the simplified system. Meanwhile, with the sensitivity information of intermediate process to inhibitors, we can compare the DS values of different inhibitors that have the same structure properties and select proper drug combinations with synergism property.
Methods for the case study of TNFαinduced NFκB pathway
Model construction and drug selection (see details in Results and Discussion)
We constructed the TNFαinduced NFκB pathway in Human umbilical vein endothelial cell (HUVEC). The model was developed based on literatures [2326]. We tweaked the parameter values in terms of the experimental data derived from HUVEC [27]. The new pathway model yielded a better simulation of NFκB activation in HUVEC. The details of Ordinary Differential Equations (ODEs) model of the NFκB pathway could be found in the Additional file 2.
Additional file 2. Summary of the model of TNFαinduced NFκB pathway. Tables S1  S4 list details of the mathematical model of TNFαinduced NFκB pathway including reactions and rate equations, initial components concentrations, kinetic parameters and ordinary differential equations.
Format: DOC Size: 117KB Download file
This file can be viewed with: Microsoft Word Viewer
According to simulation results, we made a short list of NFκB activity inhibitors covering three key nodes in the pathway including Proteasome Inhibitor II Aldehyde, HSP90 inhibitor Geldanamycin and IKKβ inhibitor PS1145 (Additional file 3, Table S5). To determine the doses of these inhibitors in our experiment, we refer to the relative IC50 values of these inhibitors taken from the published experimental and clinical data [2830]. In our model, Intercellular cell adhesion molecule1 (ICAM1) is very sensitive to TNFα stimulation. It is directly regulated by activated NFκB and become output index of downstream of this pathway.
Additional file 3. Summary of algorithm for TNFαinduced NFκB pathway. Implementation of the algorithm for TNFαinduced NFκB pathway.
Format: DOC Size: 139KB Download file
This file can be viewed with: Microsoft Word Viewer
In the simulation, the changes on relative reaction velocity constants were taken as the inhibition influence on the targets [31,32]. According to Lemma 1, we simplified the system and considered the synergism assessment factor on this simplified system. Through simulations with changing the reaction velocity constants in a wide range that responded to a wide dose range, the drug combination effects of the simplified system were gained. Additional file 3 provides a detailed description of the simulation algorithm.
Cell culture and reagents
To evaluate the computational results of our method, we conducted a cell experiment as follows. HUVECs were isolated from freshly obtained human umbilical cords by established methods [33]. The cells were grown onto gelatincoated 10 cm^{2 }culture dishes in a standard endothelial cell medium (ECM) (ScienCell Research Laboratories). ECM consists of 500 ml of basal medium, 25 ml of fetal bovine serum (FBS), 5 ml of endothelial cell growth supplement (ECGS) and 5 ml of penicillin/streptomycin solution (P/S). Cells used for this study were from passages 4 to 8 in ECM at 37°C in a 5% CO_{2}/humidified air incubator and starved for 6 hours in 0.1% FBS medium before each assay. All experiments were carried out when the cells were 80% confluent. Proteasome Inhibitor II Aldehyde and HSP90 inhibitor Geldanamycin were purchased from Alexis Biochemicals (San Diego, CA). IκB kinase (IKKβ) inhibitor PS1145 was from SigmaAldrich (St. Louis, MO). These inhibitors were dissolved in DMSO and stored at 20°C until use. The recombinant human TNFα (rhTNFα) protein was purchased from Cell Signaling Technology Inc. (CST, Beverly, MA). Antibodies to ICAM1 and βactin were obtained from CST.
Western blotting
To detect the effect of three combinations on output index of this pathway, ICAM1 directly regulated by NFκB was investigated by western blotting as follows. HUVECs were treated with 100 nM Aldehyde, 100 nM Geldanamycin, 100 nM PS1145 and various combinations of Aldehyde, Geldanamycin and PS1145 at the dose of 100 nM for 2 hours followed by 10 ng/ml TNFα. After 6 hours of treatment, wholecell extracts from treated cells and immunoblotting were prepared as previously described [34]. Whole cell lysates were subjected to SDSPAGE 10% gels. Proteins were transferred to nitrocellulose blotting membranes (BioRad Laboratories, Hercules, CA), and immunoblotted 4°C overnight with antiICAM1 and antiβactin Abs (typically 1:1,000 dilution) followed by secondary antibody conjugated with horseradish peroxidase (1:10,000 dilution). The SuperSignal^{® }West Dura (Thermo Fisher Scientific, Rockford, Ill, USA) was used for detection according to the manufacturer's instruction.
Results and Discussion
Most signaling pathways are constructed from three types of structures  serial, parallel and feedback, so our study will be focusing on these three structures. Based on Lemma 1, the simplifications of serial structures and parallel structures are studied. We also analyze the influence of feedback structure on serial structures. Numerical examples are also given; the simulation results of both original and simplified systems showed potency of the method on analyzing the combination effects. All these structures and figures are from Fitzgerald's work [1].
Simplification rule for serial structure
In Figure 1, A activates B, then B activates C. The concentration of activated C is the output of this system. The sequential activation processes construct a typical serial structure. To illustrate the drug combination effects, we designate that I_{1 }and I_{2 }are inhibitors affecting the two activation processes separately.
Figure 1. Serial structure. The left subfigure illustrate the general serial structure. The right subfigure is an example from [1]. The shaded areas are the simplified systems. Based on the simplification rules the virtual inhibitors on the simplified systems in this example target on the reaction BP+C→CP (or it also could be seen as target on component BP) and component C.
Corollary 1
The sign of the synergism assessment factor derivative DS of original serial structure in Figure 1 is opposite to the sign of the synergism assessment factor derivative DS' of the simplified structure (shaded area in Figure 1). That is to say,
Simplification rule for parallel structure
In Figure 2, there are parallel activation processes: A_{1 }activates B_{1}, A_{2 }activates B_{2}. As in serial structure, I_{1 }and I_{2 }are inhibitors affecting those two activation processes separately. Both B_{1 }and B_{2 }can activate C. The output of this system is the concentration of activated C. The relation between B_{1 }and B_{2 }could be demonstrated as logic OR.
Figure 2. Parallel structure. The left subfigure illustrate the general parallel structure. The right subfigure is an example from [1]. The shaded areas are the simplified systems. Based on the simplification rules the virtual inhibitors on the simplified systems in this example target on the reactions B_{1}P+C→CP (or it also could be seen as target on component B_{1}P) and B_{2}P+C→CP (or target on component B_{2}P).
Corollary 2
The sign of the synergism assessment factor derivative DS of original parallel structure in Figure 2 is the same as the sign of the synergism assessment factor derivative DS' of the simplified structure (shaded area in Figure 2). That is to say,
By Corollary 1 and 2, whether the simplified systems could generate synergism under equivalent inhibition can help to determine whether the original systems can generate synergism effect under actual inhibition (see details of proofs in Additional file 1). In Figure 1, the original serial structure under I_{1 }and I_{2 }inhibition can be simplified as the system (shaded) under equivalent I_{1' }and I_{2 }inhibition from the view of combination effects; It is the same to simplify the original parallel structure under I_{1 }and I_{2 }inhibition as a system under equivalent I_{1' }and I_{2' }inhibition. This could simplify the system structures, meanwhile simplify the drug combination analysis. Besides, it becomes easier to find drug combinations that could generate synergism on the systems based on these conclusions.
Combination effect preservation rule for systems with feedback
Feedback structures are common regulatory structures in biological systems, especially in signaling pathways. Positive and negative feedback structures could shape the signaling responses in time and space [35], like performancing as oscillators or bistable switches. Feedback structures increase the complexity of system structures, and make it more difficult to analyze the drug combination effects on the systems. From the viewpoint of drug combination therapy, sometimes we just need to know how the feedback structures influence the drug combination effects.
Considering an original system without feedback loop that can generate synergism effect under some drug combinations (left subfigure in Figure 3), if after adding a feedback loop (right subfigure in Figure 3), the output of the new system decreases compared to that of original system, then the feedback loop can be believed that it strengthens the drug combination effects. Upon that, it could only consider combination effect analysis on the original system.
Figure 3. Illustration of feedback structures. R is the input of x. y is the output of the system. f is the intermediate process from x to y. g is the feedback function. a and b are affected parameters.
In general feedback structures (as shown in right subfigure in Figure 3) and the system without the feedback loop can be modeled by Ordinary Differential Equations (ODEs) separately as
Lemma 2
Adding feedback loop to a system will not decrease the drug combination effects of the original system if . That is,
1) Negative feedback (g(·) < 0), and > 0
or
2) Positive feedback (g(·) < 0), and < 0 is satisfied. Details of proof are given in Additional file 4.
Additional file 4. Proof of influence of feedback. Proof of Lemma 2 for influence of feedback structure to original serial system based on comparison principle that is used to compute bounds on solutions of differential equations.
Format: DOC Size: 36KB Download file
This file can be viewed with: Microsoft Word Viewer
Examples of serial structure and parallel structure
We adapted numerical examples in [1] as shown in Figure 1 and Figure 2 respectively to verify the Corollary 1 and 2 above. All parameters are referred to [1]. In addition, all the reactions are modeled as MichaelisMenten equations. Direct evaluation of the synergism asscessment factor S by simulation leads to the results shown in Figure 4.
Figure 4. S curves of Figure 1 and 2. In the simulation, we adopted the same doses of both inhibitors like [1]. The left subfigure is the simulation results of the serial structure in Figure 1. The right subfigure is the simulation results of the parallel structure in Figure 2. The S curves of original systems are in red. The blue curves are the S curves recovered from the results of simplified systems based on Lemma 1 and the two corollaries. In fact the S curve of the simplified serial system has opposite sign to that of the original serial system. Obviously the S curves recovered agreed with the original S curves. The signs of S curves are negative showing the synergistic effect of the serial and parallel structures. In left subfigure the absolute values of S are much small (10^{3 }order) to those of the S in right subfigure (10^{1 }order). Combination therapy of the parallel structure could be more effective. These results are conincident with the results in [1].
In Figure 4, left figure shows the results of serial structure and right figure shows the results of parallel structure. The blue curves referred to the S curves recovered from results of the simplified systems, while the red lines referred to the results of the original systems. It is obvious that the recovered S curves agreed with the original S curves. Actually, the S curve of the simplified serial system has opposite sign to that of the original serial system.
S curves with negative signs indicate that both the original serial and parallel structures generate synergistic effect. In left subfigure the absolute values of S are smaller (10^{3 }order) than those of the S in right subfigure (10^{1 }order). Combination therapy of the parallel structure could be more effective. These results are conincident with the results in [1]. As a comparison from these results, the new method enables us to evaluate the combination effects of original systems by analyzing the effects of the corresponding simplified systems.
An example of negative feedback structure
Figure 5 from [1] shows the example of negative feedback structure. This negative feedback system is constructed based on the serial system in Figure 1. From the results in [1] the influence of feedback on combination effects is verified to have nonnegative impact under the condition 1) in Lemma 2.
Figure 5. Negative feedback structure. This system with negative feedback loop (from [1]) is constructed based on the serial system in Figure 1 by adding a feedback loop (orange color) which is from the output CP to the activation of B to the serial structure. The output of this system is still the concentration of activation C.
Based on the ODEs of the serial structure, the ODEs of this feedback structure are as follows:
In the model, x_{1 }and x_{2 }are concentrations of BP and CP (activated C) respectively. x_{2 }is the output of the system. x_{1t }and x_{2t }are initial concentrations of B and C respectively. R is the concentration of A. K_{m1 }and K_{m2 }are MichaelisMenten constants of the B→BP and C→CP activation and are affected by I_{1 }and I_{2}, respectively. K_{m5 }are MichaelisMenten constants of BP→B activation mediated by CP. From the simulation results in [1], it is obvious that the output of feedback structure example is smaller than the output of serial structure. This comparison verified Lemma 2 that feedback loop can strengthen the drug combination effects.
It should be pionted out, although our method leads to simplified systems, that does not mean the analysis of the simplified systems is easy. Usually simulations or experiments are still needed for analyzing these systems due to their complex dynamics [14,3638]. Our method also relies on the faithful modeling of the systems which may not be trivial since identifying system structure in general is still a challenging work [3941].
Results and discussion for the case study of TNFαinduced NFκB pathway
TNFα is both a proinflammation cytokine and a proangiogenic factor [42]. It is responsible for inflammatory angiogenesis and tumorigenesis. The induction of NFκB signaling pathway by TNFα can regulate the transcriptional expression of several genes in vascular endothelial cells that lead to angiogenesis [42]. Several standard drug targets such as HSP90 and IKKβ are among the key molecules involved in the pathway responsible for generating angiogenic factors. However, essentially all singletarget inhibitors have low therapeutic effects in inflammatory and angiogenic diseases. To find more efficient drug combination solutions, we constructed the model of an inflammatory angiogenesisrelated pathway, the TNFαinduced NFκB pathway (Figure 6A) and employed our method in this paper to search for synergistic drug combination solutions within this pathway model.
Figure 6. Case study of TNFαinduced NFκB pathway. A. The kinetic model of TNFαinduced NFκB pathway. The dashed box area indicates the simplified system. B. Biological schematic of TNFαinduced NFκB pathway. Not all arrows represent direct physical interactions. Drugs used in this study and their key targets are highlighted. C. HUVECs were pretreated with three combinations at the same concentration 100 nM for 2 hours before being stimulated with 10 ng/ml TNFα for 6 hours. Cells were lysed and output signaling protein was probed by western blotting. The level of βactin was used as an internal control for normalization. The gray values are recorded at the bottom. D. Ratio of ICAM1 to βactin was determined with densitometry. C, control, T, 10 ng/ml TNFα, G, Geldanamycin, A, Aldehyde, P, PS1145.
As a result of our method, IκB degradation proteasome, HSP90 and IKKβ as important targets in the TNFαinduced NFκB pathway are selected. In the simulation, changes on the relative reaction velocity constants were taken as the inhibition influence on targets. It should be pointed out that Geldanamycin as a HSP90 inhibitor dose not target RIP1 directly but could indirectly regulate the binding rate of RIP1, so we changed the reaction velocity constants of RIP1 binding with TNFR1 complex. The change ratio of the reaction velocity constants ranged from 0.9 to 0.0001 fold to cover a wide dose range.
As shown in Figure 6B, the TNFαinduced NFκB pathway has a serial structure with feedback. The three targets locate on the serial path (RIP1 and IKKβ) and feedback path (IκB degradation), especially IKKβ is on the joint of serial path and feedback path. Since TRAF2RIP1TRADDTNFR1TNFα activates IKKβ, the inhibition effect on RIP1 could be seen as finally acted on IKKβ. According to Lemma 1, the system could be simplified as the dashed box area in Figure 6A. The corresponding reactions inhibited by the three inhibitors in original and simplified systems were shown in Table 1. We then employed our method to calculate the synergism assessment factor on the simplified system. The sign of synergism assessment factor decides whether there is synergistic effect or not. The value of synergism assessment factor reflects the efficacy extent of drug combination from Aldehyde, Geldanamycin and PS1145 (shown in Figure 6B). We screened out the synergistic drug combinations according to the sign of synergism assessment factor in simplified system. Then we did the simulation on the original system for the synergistic drug combinations. The simulation results are shown in Table 2. It shows that all the three drug combinations could generate synergistic effect. And the signs of synergism assessment factors in simplified system were consistent with those of the synergism assessment factors in original system.
Table 1. Comparison of reactions inhibited by drugs in original and simplified systems.
Table 2. Simulation and experimental results for the case study of TNFαinduced NFκB pathway.
ICAM1, an intercellular adhesion molecule expressed in endothelial cells, is a common cellular readout of TNFα induced signaling pathway [43]. To verify the simulation results, we observed the effect of Aldehyde, Geldanamycin, PS1145 and related combinations on ICAM1 expression level by western blotting analysis. In the experiment, we used Bliss independent model to assess whether three combinations induced synergistic inhibition effect on ICAM1 expression. As shown in Figure 6C and Figure 6D, all the three drug combinations generate synergistic effect. Taken together, this result suggests that the synergistic combinations predicted by our method are qualitatively consistent with the experimental observations. Currently, the method can only make qualitative prediction on drug combinational effects. This still could provide some clues for drug combination design based on mechanisms. Besides, as the model we provided here is specific for TNFαinduced NFκB pathway in HUVEC, the general applicability of our method still needs further investigation.
Conclusions
In this paper we presented a new method based on an extended Bliss independence criterion to analyze the relationship between structures and effects for combination drug targets design from a mathematical aspect. We analyzed two classic structures, serial structure and parallel structure, and showed in steady state the sign of the synergism assessment factor derivative of the original system can be predicted by the sign of its simplified system. In addition, we analyzed the influence of feedback structure on survival ratio of the system, and showed that the feedback structure could not destroy the drug combination effect of the system without feedback under some conditions. We demonstrated by numerical examples that these results are useful for reducing the amount of computatioal load if system reaction network topology knowledge is available. In the case study, the effects of inhibitor combinations predicted by our method were experimentally validated by measuring the output (ICAM1 expression) of TNFαinduced NFκB pathway. Hopefully, this work can provide some insights to tackle the challenging problem of assessment of combination drug therapy effct in a large scale signaling pathways. As we point out, Bliss model is relatively simple, so in this paper we focused on this simple model. With deep understanding of doseeffect curves, we hope in the future our method could be expanded to more general criteria such as the law of mass action [1012,21].
Authors' contributions
HY: study conception, research design, mathematical analysis, manuscript writing; BZ: model constructing, experimental study, manuscript writing; SL: directed the design of the model constructing, experimental study, study conception, helped with the bioinformatics analyses, participated the revision of the manuscript; QZ: formulated and directed the design of the study. All authors read and approved the final manuscript.
Note
In the following figures (Figure 1, 2, 3, 4, 5), A, B, C, A_{1}, A_{2}, B_{1}, B_{2 }are components of the systems, I_{1 }and I_{2 }are inhibitors. The shaded areas represent simplified systems for the original systems. I_{1' }and I_{2' }are virtual inhibitors on the simplified systems equivalent to I_{1 }and I_{2 }.
Acknowledgements
We thank Prof. Tong Zhou and Mr. Tao Ma for helpful discussions. This work was supported by NSFC Grant Nos. (60574067, 60736027, 60721003, 60934004) and the Program of Introducing Talents of Discipline to Universities (111 International Collaboration Project of China (B06002).
References

Fitzgerald JB, Schoeberl B, Nielsen UB, Sorger PK: Systems biology and combination therapy in the quest for clinical efficacy.
Nat Chem Biol 2006, 2:458466. PubMed Abstract  Publisher Full Text

Barrera NP, Morales B, Torres S, Villalon M: Principles: mechanisms and modeling of synergism in cellular responses.
Trends Pharmacol Sci 2005, 26:526532. PubMed Abstract  Publisher Full Text

Araujo RP, Liotta LA, Petricoin EF: Proteins, drug targets and the mechanisms they control: the simple truth about complex networks.
Nat Rev Drug Discov 2007, 6:871880. PubMed Abstract  Publisher Full Text

Li XJ, Zhang HY: Synergy in natural medicines: implications for drug discovery.
Trends Pharmacol Sci 2008, 29:331332. PubMed Abstract  Publisher Full Text

Li S: Network systems underlying traditional Chinese medicine syndrome and herb formula.
Current Bioinformatics 2009, 4:188196. Publisher Full Text

Tallarida RJ: Doseresponse analysis. In Drug Synergism and Dose Effect DataAnalysis. Edited by Tallarida RJ. New York: Chapman & Hall/CRC Press; 2000:2139. Publisher Full Text

Greco WR, Faessel H, Levasseur L: The search for cytotoxic synergy between anticancer agents: a case of Dorothy and the ruby slippers?
J Natl Cancer Inst 1996, 88:699700. PubMed Abstract  Publisher Full Text

Bliss CI: The toxicity of poisons applied jointly.
Ann Appl Biol 1939, 26:585615. Publisher Full Text

Loewe S: The problem of synergism and antagonism of combined drugs.

Chou TC, Talalay P: Analysis of combined drug effects: a new look at a very old problem.
Trends Pharmacol Sci 1983, 4:450454. Publisher Full Text

Chou TC, Talalay P: Quantitative analysis of doseeffect relationships: the combined effects of multiple drugs or enzyme inhibitors.
Adv Enzyme Regul 1984, 22:2755. PubMed Abstract  Publisher Full Text

Chou TC: Drug Combination Studies and Their Synergy Quantification Using the ChouTalalay Method.
Cancer res 2010, 70:440446. PubMed Abstract  Publisher Full Text

Araujo RP, Petricoin EF, Liotta LA: A mathematical model of combination therapy using the EGFR signaling network.
Biosystems 2005, 80:5769. PubMed Abstract  Publisher Full Text

Boik JC, Newman RA, Boik RJ: Quantifying synergism/antagonism using nonlinear mixedeffects modeling: a simulation study.
Stat Med 2008, 27:10401061. PubMed Abstract  Publisher Full Text

Lee JJ, Kong M, Ayers GD, Lotan R: Interaction index and different methods for determining drug interaction in combination therapy.
J Biopharm Stat 2007, 17:461480. PubMed Abstract  Publisher Full Text

Morris JR, Boutell C, Keppler M, Densham R, Weekes D: The SUMO modification pathway is involved in the BRCA1 response to genotoxic stress.
Nature 2009, 462:886890. PubMed Abstract  Publisher Full Text

Petraitis V, Petraitiene R, Hope WW, Meletiadis J, Mickiene D, Hughes JE, Cotton MP, Stergiopoulou T, Kasai M, Francesconi A, Schaufele RL, Sein T, Avila NA, Bacher J, Walsh TJ: Combination therapy in treatment of experimental pulmonary aspergillosis: in vitro and in vivo correlations of the concentration and dose dependent interactions between anidulafungin and voriconazole by Bliss independence drug interaction analysis.
Antimicrob Agents Chemother 2009, 53:23822391. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hegreness M, Shoresh N, Damian D, Hartl D, Kishony R: Accelerated evolution of resistance in multidrug environments.
Proc Natl Acad Sci USA 2008, 105:1397713981. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Webb JL: Enzyme and Metabolic Inhibitors. New York: Academic Press; 1963:5579.

Chou TC: Theoretical basis, experimental design, and computerized simulation of synergism and antagonism in drug combination studies.
Pharmacol Rev 2006, 58:621681. PubMed Abstract  Publisher Full Text

Chou TC, Talalay P: Generalized equations for the analysis of inhibitions of MichaelisMenten and higherorder kinetic systems with two or more mutually exclusive and nonexclusive inhibitors.
Eur J Biochem 1981, 115:207216. PubMed Abstract  Publisher Full Text

Goldoni M, Johansson CA: Mathematical approach to study combined effects of toxicants in vitro: evaluation of the Bliss independence criterion and the Loewe additivity model.
Toxicol In Vitro 2007, 21:759769. PubMed Abstract  Publisher Full Text

Hoffmann A, Levchenko A, Scott ML, Baltimore D: The IκBNFκB Signaling Module: Temporal Control and Selective Gene Activation.
Science 2002, 298:12411245. PubMed Abstract  Publisher Full Text

Cho KH, Shin SY, Lee HW, Wolkenhauer O: Investigations into the analysis and modeling of the TNFαmediated NFkappaB signaling pathway.
Genome Res 2003, 13:24132422. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Cho KH, Shin SY, Kolch W, Wolkenhauer O: Experimental design in systems biology, based on parameter sensitivity analysis using a Monte Carlo method: A case study for the TNFαmediated NFkappa B signal transduction pathway.
SimulationTransactions of the Society for Modeling and Simulation International 2003, 79:726739. Publisher Full Text

Park SG, Lee T, Kang HY, Park K, Cho KH, Jung G: The influence of the signal dynamics of activated form of IKK on NFkappaB and antiapoptotic gene expressions: a systems biology approach.
FEBS Lett 2006, 580:822830. PubMed Abstract  Publisher Full Text

Zhou Z, Connell MC, MacEwan DJ: TNFR1induced NFkappaB, but not ERK, p38MAPK or JNK activation, mediates TNFinduced ICAM1 and VCAM1 expression on endothelial cells.
Cell Signal 2007, 19:12381248. PubMed Abstract  Publisher Full Text

Yoshida S, Ono M, Shono T, Izumi H, Ishibashi T, Suzuki H, Kuwano M: Involvement of interleukin8, vascular endothelial growth factor, and basic fibroblast growth factor in tumor necrosis factor alphadependent angiogenesis.
Mol Cell Biol 1997, 17:40154023. PubMed Abstract  PubMed Central Full Text

Schauer SL, Bellas RE, Sonenshein GE: Dominant Signals Leading to Inhibitor kappaB Protein Degradation Mediate CD40 Ligand Rescue of WEHI 231 Immature B Cells from ReceptorMediated Apoptosis.
J Immunol 1998, 160:43984405. PubMed Abstract  Publisher Full Text

Karin M, Yamamoto Y, Wang QM: The IKK NFkappaB system: a treasure trove for drug development.
Nat Rev Drug Discov 2004, 3:1726. PubMed Abstract  Publisher Full Text

Suresh Babu CV, Yoon S, Nam HS, Yoo YS: Simulation and sensitivity analysis of phosphorylation of EGFR signal transduction pathway in PC12 cell model.
Syst Biol (Stevenage) 2004, 1:213221. PubMed Abstract  Publisher Full Text

Hatakeyama M, Kimura S, Naka T, Kawasaki T, Yumoto N, Ichikawa M, Kim JH, Saito K, Saeki M, Shirouzu M, Yokoyama S, Konagaya A: A computational model on the modulation of mitogenactivated protein kinase (MAPK) and Akt pathways in heregulininduced ErbB signaling.
Biochem J 2003, 15:451463. Publisher Full Text

Jaffe EA, Nachman RL, Becker CG, Minick CR: Culture of human endothelial cells derived from umbilical veins. Identification by morphologic and immunologic criteria.
J Clin Invest 1973, 52:27452756. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Guttridge DC, Albanese C, Reuther JY, Pestell RG, Baldwin AS Jr: NFkappaB controls cell growth and differentiation through transcriptional regulation of cyclin D1.
Mol Cell Biol 1999, 19:57855799. PubMed Abstract  PubMed Central Full Text

Brandman O, Meyer T: Feedback loops shape cellular signals in space and time.
Science 2008, 322:390395. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wong PK, Yu F, Shahangian A, Cheng G, Sun R, Ho CM: Closedloop control of cellular functions using combinatory drugs guided by a stochastic search algorithm.
Proc Natl Acad Sci USA 2008, 105:51055110. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yang K, Bai H, Ouyang Q, Lai L, Tang C: Findingmultiple target optimal intervention in diseaserelated molecular network.
Mol Syst Biol 2008, 4:228. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sung MH, Simon R: In Silico Simulation of Inhibitor Drug Effects on Nuclear FactorκB Pathway Dynamics.
Mol Pharmacol 2004, 66:7075. PubMed Abstract  Publisher Full Text

Riel NAW, Sontag ED: Parameter estimation in models combining signal transduction and metabolic pathways: the dependent input approach.
IEE ProcSyst Biol 2006, 153(4):263274. Publisher Full Text

Quach M, Brunel N, d'AlchéBuc F: Estimating parameters and hidden variables in nonlinear statespace models based on ODEs for biological networks inference.
Bioinformatics 2007, 23:32093216. PubMed Abstract  Publisher Full Text

Jong H, Ropers D: Strategies for dealing with incomplete information in the modeling of molecular interaction networks.
Brief Bioinform 2006, 7:354363. PubMed Abstract  Publisher Full Text

Vasilevskaya IA, O'Dwyer PJ: Effects of Geldanamycin on Signaling through ActivatorProtein 1 in Hypoxic HT29 Human Colon Adenocarcinoma Cells.
Cancer Res 1999, 59:39353940. PubMed Abstract  Publisher Full Text

Gururaja TL, Yung S, Ding RX, Huang JN, Zhou XL, McLaughlin J, DanielIssakani S, Singh R, Cooper RD, Payan DG, Masuda ES, Kinoshita ET: A Class of Small Molecules that Inhibit TNFαInduced Survival and Death Pathways via Prevention of Interactions between TNFαRI, TRADD, and RIP1.
Chem Biol 2007, 14:11051118. PubMed Abstract  Publisher Full Text