Predicting Multi-Drug Inhibition Interactions based on Signaling Networks and Single Drug Dose-Response Information

Inhibiting multiple pathways/targets with combination drug therapies is widely used to maximize therapeutic benefit in many diseases. Therapeutic success can often be improved by concurrently and/or sequentially inhibiting multiple pathways/targets by combining drugs with different mechanisms of action. It is critical to understand potential drug interactions in developing and optimizing new multidrug regimens. We used data from experiments of single drugs (and few combinations) and existing signaling network knowledge from sources such as KEGG, to develop a statistical rescaling model to describe the effects of drugs on network topology. The model comprises a Hill equation for signals arriving at each receptor, a generic enzymatic rate equation to transmit signals among connecting genes, and a regression model to represent the cumulative effect of genes implicated in activation of the cell death machinery. We performed in silico experiments to derive a global sensitivity index of each term in the functional ANOVA of the dose-response model by generating doses of the drugs with the Quasi Monte-Carlo method. Only terms with large (principal) global sensitivity indices remain in the functional dose-response model. Thus, we reduce a high dimensional problem into one that can be managed experimentally. The model generates an index that indicates potential interactions of the drugs and the principal combinations, those associated with large global sensitivity indices, can be considered as candidates for further experimental validation.

Keywords: Drug combinations; Drug interactions; Signaling network; Statistical models; Function ANOVA; In silico experiments; Quasi Monte-Carlo

Combination therapy is the hallmark of therapies for cancer, viral or microbial infections, hypertension, and other diseases involving complex biological networks. For example, the transformation of a normal cell to a cancer cell involves multiple genetic mutations and up to 12 pathways [1]. Therefore, therapeutic success can often be improved by concurrently and/or sequentially inhibiting multiple pathways/targets by combining drugs with different mechanisms of action. Synergistic drug combinations, which are more effective than predicted from summing the effects of individual drugs, often achieve increased efficacy with lower doses and have reduced toxicity [2]. Since many molecularly-targeted agents exhibit synergy when used in combination with ≥2 cytotoxic drugs [2,3], significant interest remains in developing quantitative methods to detect drug synergy [4-9]. However, combination drug therapy targeting just a few gene products may be ineffective [4,10,11]. Increasing the number of agents in a combination may provide better outcomes. In the past decade, the approach to cancer therapy has been revolutionized by the identification of a variety of novel signal transduction targets amenable to therapeutic intervention. These targets were identified based on improved understanding of the molecular mechanisms of action of second messengers, other components of signal transduction pathways, and systems biology. These advances have also made available large number of potential agents and call for new quantitative approaches for combination therapy [6,12,13]. Despite the changing paradigm to target multiple pathways, methodological advances in accurately identifying drug interactions have fallen behind, as shown by a paucity of literature on the design and analysis of multi-drug combinations.

The challenge presented by multi-drug combinations is exceptional. For example, with 10 drugs, each with only 3 doses, the number of potential combinations reaches 1,048,545, a high dimensional statistical problem. Even with 5 drugs, the number of combinations is 1008. The exponentially growing number of combinations with even a few selected doses per drug quickly precludes laboratory testing. Despite the biological advances mentioned above and the importance of multi-agent combinations, current methods are mostly topological as opposed to quantitative, and do not account for high dimensionality and proper model assumptions [14,15]. Consequently, many multidrug combination studies use suboptimal experimental design and synergy analysis, where only pairwise combinations are studied and/or the dose of one or more drugs is fixed.

Currently no published methods for screening large number of combinations use both network information and single drug experimental data. One recent attempt to use network information is a search algorithm by Calzolari et al. [16] based on a deterministic model. The method can only accommodate 2 pre-fixed dose levels of each drug with no more than 6 constituent drugs. Furthermore, the method provides neither dose selection nor sufficient information on dose-response, partly because the method ignores the single drug dose-response data that is a prerequisite for drug combination studies. A method is urgently needed to enable optimal experimental design that provides the concentrations or doses of the combinations, and allow an experimentally feasible exploration of the dose-effect surface with the smallest possible sample size. We have developed a novel method to screen the large number of combinations and identify an experimentally manageable experimental design by using the dose-response data for single drugs and pathway/network knowledge to obtain an estimate of the functional structure of the dose-response relationship.

For the purposes of this study, biological networks controlled/regulated by a receptor(s) are viewed as comprising connecting genes and output nodes that are implicated in determining activation of the cell death machinery. Figure 1 presents a typical example network - apoptosis related signaling from the KEGG database (hsa04210). Different nodes have various signal propagation rules. For a given dose-level x = (x₁,x₂,…x_s)^T of drugs A₁,A₂,…A_S, denote a_0i (x) as the signal of receptor i obtained (i = 1, 2, …, r)and a_i (x) as the signal connecting gene i obtained(i = 1, 2, …, r). Gene activity levels often exhibit a non-linear relationship to their upstream regulatory signals. Typically, a Hill equation [17] can be used to model the activity a_0i (x)at receptor i,

where α_i and β_i = (β_i1,β_i2,….β_is)^T are the parameters to be estimated. To characterize the transmission of signals among connecting genes, the generic enzymatic rate equations can be used to adjust for possible feedback loops. Such equations have been motivated by various computational and biological considerations, a result of the close interaction between experimental and computational efforts [18,19]. Let a_j (x) be the activity at gene j and a_(i,j)(x) the signal sending from gene j to gene i. The activity a_i(x) at gene i is defined to be the summation of all signals a_(i,j)(x) for gene j linked up gene i, and the generic enzymatic rate equation then suggests that

where n(i) is the set of genes that signal to gene i, and ω is the expected steady state parameter. V_Fj and V_Bj are the forward and backward parameters, respectively. From the biological point of view, V_Fj can be interpreted as the steady state biomass flux sent out from gene i when the reaction is forward, while V_Bj is interpreted as the steady state biomass flux sent out from gene i when the reaction is backward. When the reaction from gene j to gene i is irreversible in the backward direction, V_Fj = 0. The number of parameters V_Fj and V_Bj may become large if many connecting genes exist in the network. The forward and backward parameters V_Fj and V_Bj of connecting gene i may differ with those of connecting gene j (i ≠ j). Statistical variations typically occur when signals pass though the network because of link instability, stochastic noise inherent in the signal propagation rules, and/or chaos phenomena from the presence of loops. To model the network efficiently, therefore it is reasonable to assume that and (i = 1,2,…,) are random effects that are independently and identically distributed (i.i.d.) normal random variables with mean µ₁ and variance σ₁².

A linear model is used to represent the cumulative effect of genes implicated at activation of the cell death machinery. For a given dose-level x = (x₁,x₂,…x_s)^T of drugs A₁,A₂,…A_S, let Y(x) be the observed viability and a(x)=(a_il(x),…a_ih(x))^T be the vector of the activities at genes i₁,….,i_h which activate the output, then we have

Where the subscript k is the k-th replication at dose-level x = (x₁,x₂,…x_s)^T, ∈_k (x)~N(0,(σ(x))²) is the measurement error of the observed data, and the standard deviation σ(x)of the measurement error may depend on the dose-level x = (x₁,x₂,…x_s)^T. u₀ is the intercept parameter, and u = (u₁,..,u_h)^T is the vector of regression parameters to be estimated. The positive parameter u_i indicates promotion by gene i; the negative parameter u_j indicates inhibition by gene j.

A training dataset is needed to model the biological network. In combination studies, the data from single drug experiments are usually available a priori. To use the network to discover multi-drug interactions, data from single drug responses cannot ensure the identification of parameters in models (1)-(3). Some limited data on the drug combinations is needed to find a model. The IC₅₀ of each drug can be estimated from the single drug dose-response curves. We recommend a relatively simple experiment of the drug combinations with each drug at its individual IC₅₀. To explore the dose-response in high dimensions, which is not well estimated by the single drug experimental data, experimentally 2s higher order combinations are considered for s drugs. Table 1 provides an example using 10 combinations of 5 drugs. The training dataset then consists of the data from single drug experiments and from limited higher order combination experiments.

Based on the training data, we can estimate the parameters in equations (1)-(3) with the maximum likelihood approach. Let β = (β₁^T,…, β_r^T )^T, α= (α₁,…, α_r)^T, and θ= (β^T, α^T, ω, u^T, σ₀², μ₁, σ₁²)^T and be the vector of all parameters to be estimated. Since a(x) in model (3) equals to zero when x = 0, the intercept μ₀ should be 100% cell viability if there is no drug intervention on the network.

Suppose that there are n distinct inputs x₁,.., x_n, and k_i replications at each input x_i, the corresponding output is Y_ij for j = 1,2,.., k_i; and i = 1,2,…, n. For given μ₁, σ₁² and a sample V_Fi and V_Bi (i =1, 2,…) from the normal distribution N(μ₁, σ₁² ), the ECM algorithm [20] can be applied to obtain the maximum likelihood estimation of β, α, ω, u, σ₀². Furthermore, for given β, α, ω, u, σ₀², we can obtain n samples of V_Fi and V_Bi using equation (3) with n distinct inputs x₁,…, x_n. The estimation of μ₁ and σ₁² can then be obtained. Details for parameter estimation are provided in the Appendix.

After we obtain the estimated parameters in models (1)-(3), the drug dose-response surface y = g(x) can be estimated in silico. A detailed and rigorous definition of a drug dose-response surface can be found in the paper by Kong and Lee [7]. The drug dose-response surface is an (s+1)- dimensional response surface in s drug combinations which is used to describe the dose-effect relationship. In this work, we would not expect the precise drug dose-response surface to be adequately estimated with limited data and in silico experiments. However, we can get sufficient information of drug interactions using the functional ANOVA [21-23], which is similar to functional principal component analysis.

Recall that x = (x₁, x₂,…x_s)^T, is the dose-level of s drugs A₁, A₂,…, A_S, and y = g(x) is the corresponding dose-response. Let g₀ = ∫_[0,1]^s g(x)dx be the overall mean of g(x). Then there is a unique decomposition

which satisfies

Moreover, (4) has orthogonal components, i.e., if (i₁,…,i_u) ≠ (j_i,…, j_v),

The functions on the right-hand side of the equation (4) are determined by

and so on. From equations (4)-(6), we have

The total and partial variances can be defined by

respectively. Note that The ratios

are called global sensitivity indices [21-23]. The integer k is called the order of the index. All R_i1,…,ikʹs are non-negative and their sum

The equality R_i1,…,ik = 0 implies that g_i1,…,ik = 0 and so the interaction of drugs A_i1…A_ik is not significant. Significance of the interaction of drugs A_i1…A_ikdecreases with decreasing R_i1,…,ik. Hence the dose-response model can be reduced if we only retain the principal terms with the largest global sensitivity indices, an approach similar to principal component analysis. It is also expected that the number of terms in the dose-response functional ANOVA representation will be reduced significantly because the cumulative global sensitivity indices of the first few terms usually contribute a dominant portion (say, 80%) of the total variation [24]. To obtain the numerical values of the global sensitivity indices, the Quasi-Monte Carlo methods for approximating the integrals can be adopted. For more details, please refer to Fang et al. [24].

To illustrate the proposed methods, we consider the apoptosis network (hsa04210) as shown in Figure 1. Two simulation experiments are conducted to investigate the effectiveness of the optimal network simulator for the discovery of multidrug interactions using the apoptosis signaling network. The first example involves a combination study of 5 drugs; the second example considers as many as 10 drugs.

Consider a combination study of 5 drugs A₁,…,A₅. The dose-levels of A₁,…,A₅ are denoted by x₁,…,x₅, respectively. Without loss of generality, assume that x₁,…,x₅ Є [0,1] . Based on the KEGG network (hsa04210), the “true” dose-response is assumed to be

Where x = (x₁,…., x₅)^T. Decompose the true dose-response (9) using the functional ANOVA representation (7) and (8), the global sensitivity indices are 38.06% for variable term x₁x₂x₃x_4, 29.61% for variable term x₁, 15.78% for variable term x₂ and 8.40% for variable term x₃, respectively. The variable terms are ignored if their global sensitivity indices are less than 1%. The response in (9) can then be approximated by

With the relative mean squared prediction error (RMSPE) given by RMSPE = {∫_(0,1)⁵ [(y(x) − y*(x))/y(x)]² dx} × 100% = 8.25%, which indicates that y*(x) is a good approximation of the true dose-response. The “true” response surface in (9) is mainly contributed by the interaction of A₁A₂A₃A₄ and drugs A₁, A₂ (the total of global sensitivity indices is 83.45%), 236 and the resulting global sensitivity indices are shown in the left pie of Figure 2.

We conducted the simulation experiments to investigate whether our method can identify correctly these true significant single drugs and drug interactions. From equation (9), the single drug dose-responses can be obtained by setting other variables to zero. For each drug, 8 concentration levels {0.0100; 0.0193; 0.0372; 0.0717; 0.1382; 0.2664; 0.5135; 0.9900} are chosen and 8 replications at each concentration with the random error ε ~ N(0,1). The generated data of single drug experiments are listed in Table S1 (Supplement).

From the single drug dose-responses, the respective IC₅₀ of each drug is estimated to be 0.07, 0.10, 0.14, 0.19, and 0.26 for drugs A₁, A₂, A₃, A₄ and A₅. To generate the data of drug combinations, 10 mixtures as shown in Table 1 are chosen and three replications at each mixture with the random error ε ~ N(0,1). Data generated from combination experiments are listed in Table S2 of Supplement.

Based on the training data, we obtained the parameter estimation for models (1)-(3), and then the dose-response function y = g(x) can be simulated for any dose-level x = (x₁,…,x₅)^T Є [0,1]⁵ of drugs A₁,…..,A₅. To estimate the global sensitivity index R_i1,…,ik, the Quasi Monte-Carlo method is used because of its improved efficiency, fast convergent rates in approximating the integrals [25,26]. With the Quasi Monte-Carlo method, we generate a set of 5-dimensional points {x^(j) : j = 1,2,…,N}, which is space filling over the experimental domain [0,1]⁵. With properly chosen x^(j) s, e.g., the uniformly scattered points, the quantity

converges to D_i1,…,ik much faster (O(N^-1(logN)^5-1)) than randomly generated x^(j) s(O(N^-1/2)). Hence, we can obtain the estimates of all global sensitivity indices R_i1,…,ik. In this simulation, 10,000 Monte-Carlo samples are generated and the estimation of the global sensitivity indices is repeated 100 times. The average of the estimated global sensitivity indices (and their corresponding standard deviations) are shown in the right pie of Figure 2. The small standard deviations indicate that the integrations are at convergence. On average, 34.04% global sensitivity index for variable term x₁x₂x₃x₄, 37.35% for variable term x₁, and 21.60% for variable term x₂. The results show that for the dose-response of the 5 drug combinations, the most significant contributions are the interaction of A₁A₂A₃A₄, drugs A₁ and A₂ (the total of global sensitivity indices is about 90%), which is consistent with the global sensitivity indices from the true dose-response.

To further illustrate the proposed methods, we considered a combination study using more drugs for the simulation. The dose-levels of 10 drugs A_1,….., A₁₀ are denoted by x_1,……x₅, respectively. Based on the KEGG network (hsa04210), the “true” dose-response is assumed to be

where x = (x₁,…,x₁₀)^T Є [0,1]¹⁰. Decompose the true dose-response (12) using the functional ANOVA representation (7) and (8), the global sensitivity indices are 41.69% for variable term x₁x₂x₃x₄x₅ 21.23% for variable term x₁x₂x₃x₄x₅x₆x₇, 11.90% for variable term x₁x₂x₃, 11.01% for variable term x₁, and 7.21% for variable term x₂, respectively. The variable terms are ignored if their global sensitivity indices are less than 1%. The response in (12) can then be approximated by

with RMSPE = 5.21%, which indicates that y*(x) provides a good approximation of the true dose-response. Thus, the “true” response surface in (12) is mainly contributed by the interaction of A₁A₂A₃A₄A₅A₆A₇, A₁A₂A₃A₄A₅, A₁A₂A₃ and drug A₁ (the total of global sensitivity indices is 85.83%); the resulting global sensitivity indices are shown in the left pie of Figure 3.

Simulation experiments are conducted to investigate whether our method can correctly identify these true significant single drugs and drug interactions in this setting. From equation (12), the single drug dose-responses can be obtained by setting other variables to zero. The corresponding IC₅₀ of each drug is estimated to be 0.07, 0.08, 0.10, 0.11, 0.14, 0.16, 0.19, 0.22, 0.26, and 0.31 for drugs A₁, …, A₁₀, respectively. For each drug, 8 concentration levels {0.01; 0.02; 0.04; 0.07; 0.14; 0.27; 0.51; 0.99} are chosen with 8 replications at each concentration with the random error ε ~ N(0,1). The data generated for the single drug experiments are listed in Table S3 (Supplement). To generate the data for drug combinations, 20 mixtures are chosen with three replications at each mixture with the random error ε ~ N(0,1) , and the corresponding data of combination experiments are listed in Table S4 of the Supplement.

Similar to Simulation 1, the estimation of parameters in models (1)-(3) is obtained using the training data. In this simulation, 10,000 Quasi Monte-Carlo samples are generated and the estimation of the global sensitivity indices is repeated 100 times. The average of the estimated global sensitivity indices (and their corresponding standard deviations) are shown in the right pie of Figure 3. The small standard deviations indicate that the integrations are at convergence. On average, 45.39% global sensitivity index for variable term x₁x₂x₃x₄x₅, 25.91% for variable term x₁x₂x₃x₄x₅x₆x₇, 9.27% for variable term x₁x₂x₃, and 8.97% for variable term x₁, respectively. The results show that for the dose-response of the 10 drug combinations, the most significant contributions are the interactions of of A₁A₂A₃A₄A₅, A₁A₂A₃A₄A₅A₆A₇, A₁A₂A₃ and drug of A₁ (the total of global sensitivity indices is about 90%), which is also consistent with the global sensitivity indices from the true dose-response.

To get an idea of the computing efficiency of the proposed method, Table 2 below provides the CUP times for Simulations 1 and 2. The CUP time for each simulation can be divided into two parts: the first part is for the ECM algorithm (parameter estimation), the second part is for the functional ANOVA (integral calculation). It can be seen from Table 2 that the 10-drug simulation is much more time-consuming than the 5-drug simulation. However, the CPU time for the 10-drug simulation is still acceptable, say that about 1 hour for the ECM algorithm and about 5 minutes for the functional ANOVA on a regular PC.

Cancer cells carry out their functions following appropriate responses to the extracellular and intracellular inputs to their complex network of signaling pathways. Many genes that code for proteins in these pathways are controlled by regulatory proteins that up-regulate or down regulate these genes depending on the inputs to the signaling network. Though striking progress is being made in extracting networks using a range of experimental data [27,28], knowledge of signaling networks remains predominantly at the level of topology rather than details on the rate constants and nonlinear message passing that occurs within the networks. Models to distinguish between members of a population of cells, for example, different cancer cells from different normal tissue types, require differences in message passing parameters and/or expression levels of the genes in the network.

Using data from experiments with single drugs (and some drug combinations) and the existing network information, we proposed statistical models to describe the drug effects on the network. Through these statistical models, we conducted computer experiments (in silico) to derive a global sensitivity index of each term in the functional ANOVA of dose-response model by generating doses of the drugs with the Quasi Monte-Carlo method. Then, we can predict the main effects that occur with combinations of multiple drugs. Two simulation studies illustrate the superior performance of our methods. The principal global sensitivity indices generally select 3 to 4 terms of multidrug combinations in the functional ANOVA model if the true dose-response function is smooth. Therefore, the model generates an index that indicates potential interactions of the drugs and the principal combinations, those associated with large global sensitivity indices, can be considered as candidates for further experimental validation.

In some signaling networks there may be loops due to feedback links among genes. Several procedures can simulate signal propagation passing networks with loops [29-31]. For any specific loop, an effective algorithm for optimal signaling with feedback loops should be developed. Furthermore, to investigate drug interactions it is important to understand the mechanisms of each drug in drug combination studies. If the relationships among drugs and receptors/or genes are known, the number of parameters in models (1)-(3) can be decreased, and some constraints can be applied to achieve high prediction accuracy. It is worth noting that the number of parameters in our network model is usually large even there are some constraints on them. The large number of parameters may cause over-fitting. One possible solution to address the potential over-fitting problem is to use the penalized likelihood approach or regularized regression machine learning, i.e., a penalty function such as Lasso [32] or SCAD [33] can be added to (A3) to achieve sparse estimation.

In this paper we have developed a novel method to screen the large number of combinations and identify an experimentally manageable experimental design by using the dose-response data for single drugs and pathway/network knowledge to obtain an estimate of the functional structure of the dose-response relationship. It is highly beneficial in bringing forth a framework for selecting drug interactions, and developing experimental designs and statistical procedures to estimate the high dimensional dose-response surface.

The authors are grateful to the reviewer for his/her constructive comments. The research of Drs. Fang, Huang and Tan is partially supported by the National Cancer Institute (NCI) grant R01CA164717. Dr. Clarke’s research is supported in part by NCI system biology grant 5U54CA149147.

The authors declare that they have no competing interests.

HBF has contributed to the problem formulation, developing the statistical models, and writing the manuscript. HH has contributed to data analysis, implementing and simulating the models. RC has contributed to the problem formulation, and writing the manuscript. MT has contributed to the problem formulation, developing the statistical methods and writing the manuscript. All authors have read and approved the final version of the manuscript.

>