Quantitative modeling of dose–response and drug combination based on pathway network

Background Quantitative description of dose–response of a drug for complex systems is essential for treatment of diseases and drug discovery. Given the growth of large-scale biological data obtained by multi-level assays, computational modeling has become an important approach to understand the mechanism of drug action. However, due to complicated interactions between drugs and cellular targets, the prediction of drug efficacy is a challenge, especially for complex systems. And the biological systems can be regarded as networks, where nodes represent molecular entities (DNA, RNA, protein and small compound) and processes, edges represent the relationships between nodes. Thus we combine biological pathway-based network modeling and molecular docking to evaluate drug efficacy. Results Network efficiency (NE) and network flux (NF) are both global measures of the network connectivity. In this work, we used NE and NF to quantitatively evaluate the inhibitory effects of compounds against the lipopolysaccharide-induced production of prostaglandin E2. The edge values of the pathway network of this biological process were reset according to the Michaelis-Menten equation, which used the binding constant and drug concentration to determine the degree of inhibition of the target protein in the pathway. The combination of NE and NF was adopted to evaluate the inhibitory effects. The dose–response curve was sigmoid and the EC50 values of 5 compounds were in good agreement with experimental results (R2 = 0.93). Moreover, we found that 2 drugs produced maximal synergism when they were combined according to the ratio between each EC50. Conclusions This quantitative model has the ability to predict the dose–response relationships of single drug and drug combination in the context of the pathway network of biological process. These findings are valuable for the evaluation of drug efficacy and thus provide an effective approach for pathway network-based drug discovery.


Background
The dose-response relation is a key topic in pharmacology. How to predict the efficacy of a compound for a system (protein, biological process, cell, tissue, organ and the body) is critical for drug discovery. The drugs (magic bullets) developed in the past decades were designed to selectively target a specific protein. However, when a single drug is administered and enters the body, interaction with 1 or more cellular targets is possible [1,2]. A drug may produce multiple effects in the system through interacting with multiple cellular targets, which is called "polypharmacology" [2,3]. The pathogenesis of complex diseases such as cardiovascular disorders and diabetes is related to a lot of genetic and environmental factors [4,5]. The human body is a complicated, integrated and networked biological system. And drugs which selectively target 1 protein cannot treat complex diseases effectively due to the robustness and redundancy of the biological system [6][7][8][9][10][11]. Meanwhile, multi-target drug therapies may be more effective than individual high-affinity drugs for complex diseases [12]. Nevertheless, the more drugs administrated or more targets with which drugs can interact, the more complicated the mechanism would be. Therefore, the prediction of drug efficacy is a challenge, especially for complex systems.
When 2 or more drugs are administrated in combination, the interactions among drugs would add a further complication to the prediction of the dose-response relation of drug combination. Generally, the drug interaction would generally produce 1 of 3 different effects: synergism, antagonism and additive effect [13,14]. Synergism means that drug combination could produce exaggerated effect, and antagonism could reduce the total effect. Synergism is especially important in clinical applications since it allows the use of smaller amounts of drugs and thus reduces the adverse effect or toxicity [14][15][16][17].
We have developed a pathway network-based approach to evaluate the efficacy of a compound against biological processes, such as blood clotting [18] and platelet aggregation [5]. Recently, we used this method for virtual screening of active compounds for the inhibition of lipopolysaccharide (LPS)-induced prostaglandin E2 (PGE2) production [19]. In this work, we demonstrate an advance in the quantitative modeling of dose-response and drug combination based on the pathway network of LPS-induced PGE2 production. PGE2 is the principal inflammation mediator, which could participate in many pathological processes [20][21][22]. The production of PGE2 can be regarded as a biomarker of inflammation. Generally, the pathway of LPSinduced PGE2 production was modeled as a network. And the binding affinity of a compound to a protein in the pathway network was assessed by molecular docking. The docking results had influence on the edge weights by relating them to enzyme efficiency via Michaelis-Menten kinetics. The effect of a compound on the entire network, and thus, in this particular case, the production of PGE2, was assessed by using the network connectivity measures (network efficiency and network flux). By integrating molecular docking and network analysis, network efficiency, network flux and their combination can quantitatively describe the inhibitory effects of compounds. Moreover, the efficacy and synergism or antagonism of the combination between 2 compounds were also evaluated.

Pathway network of LPS-induced PGE2 production
The pathway network of LPS-induced PGE2 production ( Fig. 1) comprised 30 nodes and 38 edges (arrows), where nodes represented proteins and small molecules involved in the process of LPS-induced PGE2 production, and edges meant that the node in front of the arrow was downstream in the pathway. This network was a scale-free and smallworld network, which were 2 typical characteristics of biological networks [23]. It indicated that the network can have strong stability and can resist random attacks [23,24]. However, it would be vulnerable for targeted attacks, such as selective drugs. Therefore, it offered an opportunity for us to develop drugs to treat inflammation-related diseases, especially multi-target drugs and drug combination to simultaneously block multiple targets with varying degrees.

Network efficiency and network flux
Network efficiency was first proposed by Latora V. and Nagurney A. to measure the importance of a node in a network when the node was removed [25,26]. Then it was adopted by our lab to evaluate the efficacy of a drug against blood clotting cascade [18] and platelet aggregation [5]. NE Fig. 1 The pathway network of LPS-induced PGE2 production. Circle and diamond represent small molecule and protein, respectively. The target proteins for molecular docking are marked as green diamonds was a global measure of the connectivity of a network and can reflect the integrity of the network. NE was defined as the sum of the shortest path lengths between each node in the network. Thus in the calculation of NE, all shortest paths between a node and other N-1 nodes counted (N was the number of nodes in the network). However, it cannot reflect the different weightiness of the node because all nodes were calculated N-1 times [5]. Actually the further downstream a node located in the pathway network, the more important it would be. Thus network flux was proposed to calculate the shortest paths between all other nodes and the exit node of the pathway network in our previous work [5]. The combination of NE and NF took into account the different importance of the node, thus it can be used to predict the potencies of compounds against platelet aggregation [5].
The degree of decrease of NE (NEd) and NF (NFd), and the geometric mean of NEd and NFd (NEF) were all indicators of the network connectivity. In this work, we studied 5 active compounds ( Fig. 2) extracted from Reduning Injection which was a widely used Chinese medicine prescription [19]. The activities of 5 compounds against LPS induced PGE2 production at different concentrations were predicted. However, NEd, NFd and NEF had different accuracies for the predictions, as shown in Fig. 3. Caffeic acid and Scopoletin were the 2 of the most potent compounds, so the 2 compounds were picked out as examples. In the case of Caffeic acid (Caa), the predictions of NEd were lower than the in vitro experimental results, while the predicted inhibition rates by NFd were higher than experimental results (Fig. 3a). It was more complicated for Scopoletin (Sco): the model had higher predictions at low concentrations and lower predictions at high concentrations by NEd, while it was the direct opposite of predictions by NFd (Fig. 3b). However, the model gave good agreements between predictions by NEF and experiment results in all cases. Therefore, NEF was a better evaluation indicator for this system and was used in further evaluations.

Dose-response curve
The shape of the dose-response curve was important to evaluate the efficacy of a compound. All predicted and experimental dose-response curves were sigmoid. Table 1 listed the parameters of fitted dose-response curves of 5 active compounds according to the predictions by NEF. The E max and E min were close to 100 % and 0, respectively. And the correlation coefficients were higher than 0.999, especially for Caa and Sco (Figs. 4a, b). Moreover, the predicted EC50 values of 5 active compounds well matched with the experiment values (R 2 = 0.93, Table 1), which indicated that the prediction model by NEF would be reliable.

Drug combination and dose-response surface
This approach can also evaluate the combination of 2 or more compounds. There would be 3 effects for drug combination: synergism, additive effect and antagonism. The combination index (CI) proposed by Chou T. C. was adopted to quantify the synergic degree of drug combination [13]. When 2 drugs both existed in the system, it would produce a dose-response surface (DRS). Fig. 5a showed the DRS of the combinations between Caa and Sco in different doses. A series of drug combination can have the same effects on the system, which can be described by the isobologram (Fig. 5b). Each dose pairs on the isobologram represented possible combinations that produced the equivalent effect. In the case of combination of Caa and Sco, it was a typical synergistic effect. However, each drug pair in the isobologram differed in degree of synergism. The stars in Fig. 5b pointed out the optimal combination for each degree of inhibition. And we found that the dose ratios of 2 compounds for maximal synergism were nearly the same with the ratio between each EC50. Therefore, the experimental inhibition rates of 6 combinations of Caa and Sco with the ratio between each EC50 value were determined. And the results agreed well with the predictions (R 2 = 0.84, Fig. 5c).

Discussion
The pathway of a biological process is a minimal biosystem with a specific function and can be abstracted as a network. The network properties can relate with the state of the biosystem to a certain extent, especially for biomarkers related to diseases [5]. Thus the influence on the pathway network of a compound can be used to evaluate the efficacy [5,18,23,[27][28][29]. These results demonstrated that the degree of the decreases of NE, NF and NEF were measures of inhibition of a drug against LPS-induced PGE2 production.
A drug could target multiple proteins in the biological pathway, and a drug combination could produce synergistic or antagonistic effect in different extent through multitarget interactions [14,19,30]. Synergism is useful in illuminating mechanisms of drug action and exploring computational models to predict new drug leads in drug discovery. The dose-response surface and isobologram are 2 practical tools. Moreover, the combination index is a convenience for researchers to determine whether synergism, additive effect, or antagonism exists in a drug combination. The exact mechanisms of the inhibitions of drugs against LPS-induced PGE2 production are unclear now since the biological system of the cells or organism is complicated. Our recent works indicated that most active compounds would have polypharmacology according to drug-target network [30]. However, the computational approach in this work generally does not need to know about the exact mechanism, which could broaden the scope of application, especially for complex systems.
Although in the above text the predictions of this model agreed well with the experimental results, it would be necessary to note that this approach should require several conditions to obtain reasonable predictions. First, the pathway of a biological process should be as fully accurately as possible. The calculations of   NE and NF are heavily dependent on the completeness of the pathway network. The information of the pathway of LPS induced PGE2 production is abundant and this system has been reviewed in literatures and databases (described in the section of Construction of the pathway network of LPS-induced PGE2 production). However, when others apply this approach in this work to other systems, they need to be very careful. When a pathway has low completeness or lacks sufficient annotations, the predictions would have large deviations. Second, the structures (determined by X-ray or NMR) of most target proteins in the pathway should be known. It's best to obtain ligand-protein complex structure and thus the binding site can be defined as the space which was occupied by the ligand. Beyond that, the binding energy calculated by molecular docking or molecular dynamics simulation should be accurate. The binding predictions from molecular docking are subject to a margin of error due to the principle and method of the calculation. However, we can try to reduced the error. For example, we used the most commonly used software AutoDock4 and adopted the validated protocols whose predictions had been validated by experimental results in our previous works [5,18]. Big errors of binding predictions may affect the calculation of edge values and then reduce the accuracies of the calculation of NE and NF. Furthermore the predictions of inhibition rates of drugs would have big errors. Finally, the prediction model should be validated by experiments. The known active compounds can be used as training set to adjust the parameters of the predicting model.

Construction of the pathway network of LPS-induced PGE2 production
The pathway network of LPS-induced PGE2 production was constructed in our recent work [19] according to the information extracted from KEGG pathway database [31], Reactome [32], and literatures [33][34][35][36][37][38][39][40]. Generally, LPS-induced PGE2 production was involved in 2 pathways: Toll-like receptor signaling pathway (ID: map04620 in KEGG pathway database) and NF-kappa B signaling pathway (map04064). First, LPS can interact with CD14 and the complex facilitates the recognition of LPS stimulation by TLR4. Then the signaling is divided into MyD88-dependent and TRIF-dependent pathways. MyD88 and TRIF can activate transcription factors such as IRF-5, NF-kappa B and AP-1 in the downstream pathway. In particular, Eliopoulos A. G. and colleagues contributed the pathway of CREB, a key regulator of COX-2 transcription [39]. Finally, the pathway comprised 30 nodes and 38 edges (Fig. 1). Cytoscape 2.8 was used to visualize the network and calculate the network properties by Network Analysis plugin [41].

Molecular docking
There were 14 important proteins (  [42,43] according to the protocols described in previous works [5,19,44]. The energy grid was a 20 × 20 × 20 Å cube centered on the occupied space of the original ligand with a spacing of 0.375 Å between the grid points. The maximum number of energy evaluations was set to 2.5 × 10 7 . The AD4score function was used to evaluate the affinity between compound and protein, and the docking score was pK i .

Calculation of network efficiency and network flux
According to Michaelis-Menten equation and the law of mass action, the rate equation in presence of 1 noncompetitive inhibitor I was: where K m , K I , [S] and [I] are Michaelis constant, inhibition constant of I, concentration of substrate S and inhibitor I; v 0 is the activity of the enzyme without inhibitor [45]. Thus we defined the fraction of affection (f a ) to quantify what the percentage the enzyme was inhibited: When 2 mutually exclusive inhibitors (X and Y) both existed in the system, f a would be [45]: In the pathway network, the value of an edge (EV) represented the resistance in signal transduction. That is, when a target protein was inhibited, the value of the edge which came out from the target protein would enlarge to accommodate it. We arbitrarily set the initial (default) EV for each edge as 1. We arbitrarily assumed that the most potent inhibitor can block the target 99.5 %, thus we defined the highest EV as 200 (1/(1-99.5 %)). Accordingly, the EV at different concentration of inhibitor I would be: when 2 inhibitors both existed. Network efficiency and network flux were both measures of the network connectivity [5,18,26]. NE was defined as the sum of the reciprocals of the shortest path lengths between all pairs of nodes in the pathway network: NF was defined as the sum of the reciprocals of the shortest path lengths between other nodes and the exit of the pathway network: where d ij is the length of the shortest path between nodes i and j. The calculation programs of NE and NF were written in C++ language using the Dijkstra algorithm.
In order to evaluate the influence of a compound on the pathway network, we defined the NEd as the degree of decrease of NE as following: where NE max and NE min are the maximal and minimal NE when all EVs are set as 1 and 200, respectively. Similarly, the NFd was defined accordingly: Finally, we defined the NEF as the geometric mean of NEd and NFd to evaluate the impact of a compound on the pathway network comprehensively: Fitting of dose-response curve Typically, the dose-response relation can be simulated by the following equation [46]: where E max and E min are the top and bottom asymptotes of the response, y is the inhibition rate when the concentration of the inhibitor is [I], EC50 is the concentration of inhibitor at half-maximal effect, and n is the slope parameter like the Hill coefficient [47]. The fitting of computational efficacy or experimental results versus the concentration of the inhibitor was performed and the correlation coefficient was used to evaluate the reliability of the model.

Combination index
The effect of the combination of 2 drugs may be simple additive, exaggerated (synergistic) or attenuated (antagonistic). In order to quantify the synergism or antagonism for 2 drugs (D 1 and D 2 ), Chou T. C. introduced the term combination index [13]: where CI <1, =1 and >1 indicated synergism, additive effect and antagonism, respectively. (D x ) 1 and (Dx) 2 represented the concentrations when D 1 and D 2 alone can inhibit the system x %. (D) 1 and (D) 2 were the concentrations when D 1 and D 2 in combination can inhibit the system x %. . The experiments were repeated 3 times at each concentration of each compound. The inhibition rate of a compound against LPSinduced PGE 2 production was calculated by: where C(control) was the background concentrations of PGE2. C(drug) and C (model) represented the concentrations of PGE2 when the RAW246.7 cells were incubated with drug or DMSO and then stimulated by LPS, respectively.

Conclusions
In this study, we developed a quantitative model to predict the dose-response curves of single drug and drug combination based on the pathway network of LPS-induced PGE2 production. The network efficiency and network flux are both measures of the connectivity of the pathway network. And thus the degrees of the decrease of NE, NF and NEF could evaluate the efficacy of a drug to the biological systems. By integrating molecular docking and network analysis, the dose-response relationships of 5 compounds against LPS-induced PGE2 production were evaluated and the predictions agreed well with experimental results. Furthermore we explored the dose-response relationships of drug combinations to study the synergism. Moreover, identifying novel effective drug combinations or multi-target agents is a new trend in drug discovery. Actually, the herb medicines are the natural combinations of active compounds. This computational method can be used to quantitatively evaluate the efficacy of a mixture of 2 or more drugs, even herb medicines. The increasing of the complexity of multiple omics data sets requires more effective approaches for drug screening. And computational modeling is an indispensable tool for understanding dose-response relationship and mechanisms of a drug or drug combination. This work would provide a new computational approach to evaluate drug efficacy before clinic trials and screen optimal combination for drug discovery when a biological system/process has a well-defined pathway.