A formal model for analyzing drug combination effects and its application in TNF-α-induced NFκB pathway
© Yan et al. 2010
Received: 19 September 2009
Accepted: 25 April 2010
Published: 25 April 2010
Skip to main content
© Yan et al. 2010
Received: 19 September 2009
Accepted: 25 April 2010
Published: 25 April 2010
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.
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.
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.
Drug combination therapy is commonly used in clinical practice [1–3]. 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 long-standing 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 ; given that two similar drugs competitively acting on a target, Loewe additivity model is used to predict the combined effect of two drugs . Chou-Talalay further proposed the Combination Index (CI) theorem, serving as a general expression and quantification of drug interaction based on the mass-action law in biophysics and biochemistry [10–12]. These models are widely used in in-vitro and in-silico experiments for drug targets design and dose-response relationship analysis to instruct the selection of drugs and design of combination scheme [13–15]. 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 cost-effective 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 [16–18], we start our studies from this simple model and operationally provide a combination effect assessment index inspired by the Combination Index theorem [10–12]. 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 angiogenesis-related 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 PS-1145, in NFκB pathway. The simulation results suggested that three inhibitor combinations yeilded significant synergism and were validated the simulated results by cell experiments.
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 .
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.
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  and the mutually non-exclusive case in . Compared with the critical point (CI = 1) of Chou-Talalay's Combination Index, we used S = 0 as the critical point to determine whether there is sysnergistic effect. Under the first order mutually non-exclusive 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 in-vitro or in-silico experiments, so it is convenient to verify whether synergism is generated under specific drug combination with this criterion . Therefore, it is widely used in combination therapy design [16–18]. 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 dose-response relationship. The dose-response relationship could be get through series of in-vitro experiments costly under different doses. Sometimes the dose-response relationship is assumed to have some special form like Hill equation to reduce the experiment costs . 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.
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.
The extended Bliss independence criterion introduced in Eq. 5 enables us to find special structures of systems that simplification is possible for synergism-generating 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.
Where , x = Φ (a), y = Ψ (b), x 0 = Φ (a 0), y 0 = Ψ (b 0). Details of proof are given in Additional file 1.
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.
We constructed the TNF-α-induced NFκB pathway in Human umbilical vein endothelial cell (HUVEC). The model was developed based on literatures [23–26]. We tweaked the parameter values in terms of the experimental data derived from HUVEC . 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.
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 PS-1145 (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 [28–30]. In our model, Intercellular cell adhesion molecule-1 (ICAM-1) is very sensitive to TNF-α stimulation. It is directly regulated by activated NFκB and become output index of downstream of this pathway.
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.
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 . The cells were grown onto gelatin-coated 10 cm2 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% CO2/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 PS-1145 was from Sigma-Aldrich (St. Louis, MO). These inhibitors were dissolved in DMSO and stored at -20°C until use. The recombinant human TNF-α (rh-TNF-α) protein was purchased from Cell Signaling Technology Inc. (CST, Beverly, MA). Antibodies to ICAM-1 and β-actin were obtained from CST.
To detect the effect of three combinations on output index of this pathway, ICAM-1 directly regulated by NFκB was investigated by western blotting as follows. HUVECs were treated with 100 nM Aldehyde, 100 nM Geldanamycin, 100 nM PS-1145 and various combinations of Aldehyde, Geldanamycin and PS-1145 at the dose of 100 nM for 2 hours followed by 10 ng/ml TNF-α. After 6 hours of treatment, whole-cell extracts from treated cells and immunoblotting were prepared as previously described . Whole cell lysates were subjected to SDS-PAGE 10% gels. Proteins were transferred to nitrocellulose blotting membranes (Bio-Rad Laboratories, Hercules, CA), and immunoblotted 4°C overnight with anti-ICAM-1 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.
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 .
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 I1 and I2 inhibition can be simplified as the system (shaded) under equivalent I1' and I2 inhibition from the view of combination effects; It is the same to simplify the original parallel structure under I1 and I2 inhibition as a system under equivalent I1' and I2' 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.
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 , 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.
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
2) Positive feedback ( g(·) <0), and <0 is satisfied. Details of proof are given in Additional file 4.
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 . 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.
In the model, x 1 and x 2 are concentrations of B-P and C-P (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 Michaelis-Menten constants of the B→B-P and C→C-P activation and are affected by I1 and I2, respectively. K m5 are Michaelis-Menten constants of B-P→B activation mediated by C-P. From the simulation results in , 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, 36–38]. 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 [39–41].
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.
Comparison of reactions inhibited by drugs in original and simplified systems.
in original system
in simplified system
RIP1 binding with TNFR1 complex
IKKK-P activates IKK-β
IKK-β-P phosphorylates IκB
IKK-β-P phosphorylates IκB
Simulation and experimental results for the case study of TNF-α-induced NFκB pathway.
Synergism Assessment Factor (Simulation)
Synergism Assessment Factor (Experiment)
ICAM-1, an intercellular adhesion molecule expressed in endothelial cells, is a common cellular readout of TNF-α induced signaling pathway . To verify the simulation results, we observed the effect of Aldehyde, Geldanamycin, PS-1145 and related combinations on ICAM-1 expression level by western blotting analysis. In the experiment, we used Bliss independent model to assess whether three combinations induced synergistic inhibition effect on ICAM-1 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.
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 (ICAM-1 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 dose-effect curves, we hope in the future our method could be expanded to more general criteria such as the law of mass action [10–12, 21].
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).
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.