Abstract
Background
Cells constantly sense many internal and environmental signals and respond through their complex signaling network, leading to particular biological outcomes. However, a systematic characterization and optimization of multisignal responses remains a pressing challenge to traditional experimental approaches due to the arising complexity associated with the increasing number of signals and their intensities.
Results
We established and validated a datadriven mathematical approach to systematically characterize signalresponse relationships. Our results demonstrate how mathematical learning algorithms can enable systematic characterization of multisignal induced biological activities. The proposed approach enables identification of input combinations that can result in desired biological responses. In retrospect, the results show that, unlike a single drug, a properly chosen combination of drugs can lead to a significant difference in the responses of different cell types, increasing the differential targeting of certain combinations. The successful validation of identified combinations demonstrates the power of this approach. Moreover, the approach enables examining the efficacy of all lower order mixtures of the tested signals. The approach also enables identification of systemlevel signaling interactions between the applied signals. Many of the signaling interactions identified were consistent with the literature, and other unknown interactions emerged.
Conclusions
This approach can facilitate development of systems biology and optimal drug combination therapies for cancer and other diseases and for understanding key interactions within the cellular network upon treatment with multiple signals.
Background
Understanding how multiple signals affect cellular functions is necessary in order to be able to understand and control these functions. Extensive studies have been done to address how the activation/inhibition of a particular cellular signaling pathway leads to a specific response. Several challenges limit the ability to study the simultaneous effects of multiple signaling. The complexity and lack of detailed knowledge of cellular systems prevent, in many cases, accounting for the effects of some unknown interactions among pathways or among nonprimary signal targets. In addition, genetic or epigenetic alterations between otherwise similar cells can cause a significant difference in their responses. This places additional constraints on the experimental outcomes obtained by analyzing individual components. Furthermore, a critical challenge in the investigation of the effects of multiple signals is the arising complexity associated with the increasing number of signals and their various intensities. Without a systematic approach to replace a large number of time and resource consuming experimental tests, it is difficult to characterize the effects of these signals, to identify appropriate signal combinations.
There has been an increasing interest in examining how various biological activities are regulated by multiple interacting signals [14]. Berenbaum introduced a direct search method to optimize cancer chemotherapy regimens [5]. Recently, a method based on stepwise direct search for identifying optimal combination of drugs for pain treatment has been introduced [6]. The method can also be applied in clinical research. More recently, a biased random walk approach called the "modified Gur game" approach was introduced to identify potent drug combinations [7,8]. It is applied towards an objective with a "small" number of experimental trials. While the goal of these studies is to achieve optimization with minimal number of tests, the approach in these studies has several limitations including sensitivity to the design of the automatons driving the random walk and sensitivity to initial conditions. Its capacity to compare the performance of multiple systems will be limited due to the limited amount of obtained information. Moreover, the approach does not guarantee convergence to local or global maxima. In [8], the modified Gur game approach was used to identify a wide range of drug concentrations for which a stochastic search algorithm, differential evolution, was used to maximize an objective function. Although this approach converges to better combinations, the determination of the range of drugs to be used in the combination is sensitive to initial conditions. Another recent and very similar approach uses a deterministic search algorithm for optimization of drug combinations [9]. Determining optimal combinations for systems where a mechanistic model based on massaction kinetics was recently presented [10]. The use of search algorithms as well as other systems approaches that include the mechanistic massaction models were reviewed in [11]. Another limitation of these approaches is that they require repetition of the experiment in case the optimization parameters are to be modified or there is a change in the objective function. This limitation becomes significant when considering multiobjective optimization functions in which the objective function is dependent on subjective parameters, resulting in the need to carry over several experiments to determine a suitable set of parameters. Additionally, the convergence of the experimentally applied search algorithms depends on the appropriate selection of the algorithm's parameters. In the absence of a model to enable reasonable selection of the algorithm's parameters, convergence of these algorithms might be compromised. Furthermore, the experimentally applied algorithms can converge to an optima (local or global) that is not very robust to "small" variations in the input signals.
Other recent work on identifying drug combination focuses on identifying mixtures of drugs where the the search space is reduced to use only a single concentration while the search space is increased by searching through a larger number of potential drugs [12,13]. The latter example describes an approach and applies it towards finding promising mixtures for lung cancer.
In this work, we achieve the desired goals through the integration of datadriven mathematical tools with biological measurements to generate quantitative models of cellular functions (Figure 1). Instead of mapping physical interactions, the resulting model is a quantitative model mapping particular signals to their cellular process responses. The responses represent the net change in certain cellular activities caused by signal interactions within a large and complex network. The model is generated using a suitable mathematical approximation method, which relies on testing a relatively small subset of all possible signal combinations and is capable of predicting the response to the complete set of signal combinations. Through running in silico experiments, the model enables analyzing the response of the system to various combinations and determining or selecting subsets of signal combinations that can yield desired cellular responses. The determination of these subsets can be achieved using tools such as stochastic search algorithms and cluster analysis. The proposed approach will facilitate the understanding of fundamental cellular responses, which are system responses reflecting the activity of a complex signaling network controlled by multiple internal and external signals. This approach can promote efficient understanding of cellular functions without intermediates. In addition, the approach allows multiple cell types or other biological systems to be quantitatively characterized, modeled, and compared in parallel. The maximal difference or similarity can be identified using a computational search. It can facilitate the development of drug combination therapies for various types of cancers [14,15].
Figure 1. Summary of approach. (A) The proposed approach uses inputoutput data to generate mathematical models capable of predicting the cellular responses to input combinations. The models enable using other mathematical tools for analyzing the cellular responses and for selecting the appropriate combinations of the input signal to drive the system to respond favorably. (B) The desired drugs are combined in certain concentrations and a few combinations are chosen and evaluated experimentally. A predictive mathematical model is generated that can predict the response to all possible combinations. The model can be used to analyze and predict drug interactions and their effects on the observed cellular response and can also be used to determine effective combinations.
In comparison to the approaches previously mentioned, the approach introduced in this work overcomes their limitations by identification of the complete response function and carrying out the optimization in silico. The cost of carrying such in silico experiments is significantly less and is generally faster. The identification of the system response function also provides additional information regarding the potential of using a smaller number of drugs and on key systemlevel signaling interactions.
Results
We utilize mathematical tools to characterize the effects of three and four agents on the differential response of cancer and normal cells. The cellular ATP levels of a nonsmall cell lung cancer cells A549 and primary lung fibroblast culture of AG02603 cells in response to combinations of three chemical agents are measured. Mathematical tools are used to construct predictive models of the cellular ATP levels in response to the combinations. We examine the ability to utilize relatively small numbers of combinations for model generation. The results are extended to study systems of higher complexity with the addition of a fourth chemical agent. The resulting models are used to compare the responses of normal and cancer cell to the same set of combinations. We show that a properly selected combination can result in a significant difference in the respective performance of normal and cancer cells. We also experimentally validate the selected combinations. Furthermore, we show that a combination of chemical agents, if properly chosen, can be more effective than a single agent in inducing a differential response between normal and cancer cells. Using the models, we will examine all the possible lower order mixtures of the four drugs. Moreover, we extract key systemlevel signaling interactions and compare these interactions between different cell types. We also compare these interactions to known interactions between the drugs.
SignalCellular Response Modeling with a Complete Data Set
Inhibition of cell survival and proliferation has been a widelyused approach in cancer treatment [16]. We investigated the combined effect of several drugs that target critical cellular signaling pathways for cell survival and proliferation. Three drugs AG490, U0126, and indirubin3'monoxime (I3M), which target three distinct while connected signaling pathways critical to most cancer and noncancer cells, were chosen in our study (Figure 2).
Figure 2. Simplified pathway and drug interactions. Shown are simplified pathways targeted in the three and fourdrug combination treatment of nonsmall cell lung cancer cells A549 and primary fibroblast AG02603 cells that are already known or reported. The dashed arrows indicate indirect connections.
One of the goals of this work is to identify differences in the responses between cancer and noncancer cells upon the drug treatment. The interactions in Figure 2 are an oversimplified set of interactions of the drugs used. The simplified diagram serves to illustrate some of the known interactions within the cell upon treatment with various drugs. AG490 is a tyrosine kinase inhibitor; U0126 is a MEK inhibitor; and indirubin3'monoxime is a cyclindependent kinase inhibitor. The drugs are also known to inhibit other targets in addition to the intended target enzymes and as such can lead to unknown interactions. Moreover, each pathway has various interactions with more pathways that are not depicted. The drugs may activate some pathways such as some of the stress responses, directly or indirectly. Thus, the interactions of these drugs and their combinatory effects are difficult to predict.
Cellular ATP represents one of the most common and essential markers for all live cells. Measuring ATP is a generally accepted quantitative and sensitive assay for assessing the inhibition of cellular growth, proliferation, and induction of cell killing by drugs. The cellular ATP level, which is regulated by multiple cellular pathways, was experimentally quantitated [17]. Total cellular ATPlevels of A549, a nonsmall cell lung cancer cell line, and AG02603, a normal fibroblast cell culture, were measured 72 hours after drug treatment and normalized to untreated cellular ATPlevels. The drug response dose curves were measured for each of the three drugs (Figure 3). The three drugs (each with eight different concentrations including zero) comprised 512 possible combinations in total. The doses of individual drugs were chosen based on the individual doseresponse curves and covered the concentration ranges that resulted in minimal to maximal cell inhibitory effect (Table 1).
Figure 3. Singledrug dose response curves. Shown are the experimental singledrug dose response curves for the four drugs used in the study. The data was used to identify the drug concentrations to be used in combination studies.
Table 1. Concentrations of the drugs used in the three and four drug treatments.
The ATP level in response to all of the 512 drug combinations was experimentally measured in lung cancer A549 cells and in primary lung fibroblast AG02603 cells. The fibroblast cells were derived from normal healthy tissue and are not cancer cells. There are several mathematical methods that can be used to generate models of inputoutput data. Here we provide a comparison of some of these methods in view of the function approximation problem considered. The methods include two neural network structures and two linear regression models. The neural network structures are a single layer multilayer perceptron (MLP) and a cascaded neural network [18,19]. We have examined different numbers of neurons per layer for each of these neural network structures. The results below show that a fourneuron singlelayer MLP is sufficient for the purposes of this work. For the cascaded network, two layers with a single neuron per layer were sufficient. Networks with more neurons per layer also produced satisfactory results. The two linear regression models involve different nonlinear regressors. The first one uses interaction terms that are pairwise and kwise products of all the concentrations of the drugs. The second is a quadratic response surface that uses only pairwise products and quadratic terms of the concentrations (See the Methods Sections).
The different models were trained against 80 out of 512 points with the goal of minimizing the mean square error of prediction. The outputs of the models are processed through a saturation function to limit outputs to the interval [0,1]. The trained models can predict the responses to all 512 combinations with high fidelity. The correlation coefficients between the predicted normalized ATP levels and their corresponding experimentally measured values are higher than R = 0.91 (Figure 4A). Looking at only the points that were not used for training, i.e., the 432 points, the correlation coefficients between the predicted normalized ATP levels and their corresponding experimentally measured values were also high (Table 2). We have also examined how the different models compare to each other by measuring the correlation between their predicted values (Figure 4B). Examining the correlation between the different models, we see that there is also strong correlation with Rvalues higher than 0.94.
Figure 4. Threedrug model fitting. The figure shows an evaluation of modeling the responses of A549 and AG02603 to combinations of three drugs. (A) The panels show a plot of the model predicted cellular ATP levels versus the experimentally measured values. Cellular ATPlevel predictive models for A549 and AG02603 cells were developed using a number of different methods. a. A linear regression model that uses pairwise products of concentrations and quadratic terms (QRF). b. A linear regression model with nwise products of concentrations (LR). c. A cascaded neural network with two singleneuron layers (Cascaded NNet). d. A fourneuron single layer multilayer perceptron artificial neural network (MLP). The models are based on fitting 80 out of 512 combinations and the figure shows the predicted versus experimental values for all 512 combinations. The correlation between the experimentally tested cellular ATPlevel (xaxis) and the predicted cellular ATPlevel (yaxis) is shown. The circles in the graphs represent individual data points. The diagonal line represents a perfect fit between the experimental and predicted data. (B) Comparison between the predicted normalized ATP levels of different models. The predicted ATP levels for different models are plotted against each other. a. QRF versus LR. b. Cascaded NNet versus MLP. c. QRF versus MLP. d. LR versus MLP. The correlation coefficients between the different methods are shown.
Table 2. Correlation coefficients between the 432 points not used for training and the corresponding experimentally measured values.
Experimental testing of all possible combinations can be a costly process. Whenever the response of the biological system is smooth enough, we can utilize a smaller number of combinations to map the entire response surface. In this regard, we have examined the effects of using varying numbers of points to fit the models on the accuracy of prediction. The four methods discussed above were considered and models using 10, 20, 40, 80, 160, and 320 points were fitted. To mimic an actual experimental setup, the points were randomly selected out of the 512 possible combinations using a uniform distribution. The mean square error of prediction for each of the methods and fitting data for both cell types indicate that increasing the number of points reduces the mean square error (Figure 5).
Figure 5. Model accuracy as a function of modeling method and data size. The figure shows the effects of using different sizes of data sets for fitting the models and the effects on the model accuracy as measured by the mean square error of prediction. The mean square error of predicted data using four different methods and using 10, 20, 40, 80, 160, and 320 data points, for both A549 and AG02603 cells are displayed. The four methods are two linear regression methods and two neural network methods. The linear regression methods include one that uses pairwise products of concentrations and quadratic terms of the concentrations (QRF), and the other uses the nwise products of concentrations of drugs as interaction terms (LR). The two neural network methods are a cascaded neural network with two singleneuron layers, and a fourneuron singlelayer multilayer perceptron (MLP). As the number of points used to generate the model increases, the mean square error decreases.
However, no significant improvement in the errors are observed for models with more than 80 points. Using a small number of points results in poor prediction with the linear regression models, the interaction model and the quadratic models. In the absence of postprocessing of the data by passing through a saturation function, the mean square error of prediction of the linear regression models become significantly worse. Between the two linear regression models, the quadratic model performs better than the interaction model, potentially due to the added quadratic terms, suggesting the nonlinearity of the response of these cells to the drug combinations used. These models also have higher mean square errors than the neural network models. The differences between the mean square error of these models tend to diminish as the number of points used increases. The two neural network models were comparable overall.
The number of data points required to generate a valid model relies on factors including intrinsic signalresponse relationships for individual cell cultures and the experimental measurement error. In addition, the smoothness of many signalresponse relationships enables the modeling to rely on less dense mapping over small ranges of signal concentrations. Our results suggest that with a proper mathematical modeling method, the effect of signal combinations can be systematically described through randomly testing a relatively small percentage of signal combinations within specified concentration ranges.
Characterization of SignalCellular Response Relationships in Systems of Higher Complexity
The simulated models capable of systematically describing the signalcellular response relationships for various cells enable the comparison of the celltype specific differences in cellular responses to multidrug treatments. An interesting question is whether our approach can simulate more complex systems with a relatively small set of experimental data, and whether efficient multidrug combinations that lead to a high level of A549 cell inhibition while preserving AG02603 cells can be identified among more drugs. In this regard, we added GF 109203X, a PKC inhibitor, to the three drugs tested in the above section. The addition of a fourth drug increases the search space to 2401 combinations (four drugs with seven concentrations each). Models representing the relationship between the four drugs and the ATPlevels of A549 and AG02603 cells were generated based on 148 experimentally tested combinations. We fitted a singlelayer neural network with four neurons using the experimental data. The correlation coefficients between the predicted data from this model and the experimental data were 0.98 for A549 and 0.97 for AG02603 (Figure 6). Overall, there is a very high agreement between the predicted and experimental results (Figure 7).
Figure 6. Fourdrug model fitting. The figure shows the modelpredicted ATPlevels versus the experimentally obtained values for four drugs combinations. The model used is a four neuron single layer neural network model that was fitting using a data set of size 148 out of 2401 combinations. The correlation between the experimentally tested cellular ATPlevel (xaxis) and the predicted cellular ATPlevel (yaxis) is shown. The circles in the graphs represent individual data points. The diagonal line represents a perfect fit between the experimental and predicted data.
Figure 7. Analysis of drug combination differential responses. The figure presents the analysis results to identify combinations of drugs that can maximize a performance function reflecting the desired objective of minimizing the killing of AG02603 cells and maximizing the killing of A549 cels (see methods section for the details on the performance functions. (A) Results of kmeans clustering analysis of the performance of fourdrug combinations. Points are grouped together to minimize the distance to the mean of that group. Each group is given a identification number (ID). The xaxis shows cluster IDs; the yaxis shows the performance score of each drug combination. The heat map is a function of the performance for each point. (B) Cell ATPlevel (xaxis: AG02603; yaxis: A549) upon fourdrug combination treatments. The red stars represent the experimentally tested 148 drug combinations. The squares represent the whole set of 2401 possible drug combinations with the heat map colors reflecting the performance of each combination.
We have also investigated the ability of these models to predict the experimental data obtained in the three drug combinations experiment. These combinations correspond to fourdrug combinations where the concentration of the fourth drug is set to zero.
Correlation coefficients between the 512 data points and their predicted values were 0.86 for A549 cells and 0.88 for AG02603. These correlation coefficients are quite reasonable given that the two experiments were conducted several months apart and the model was trained on a separate and independent data set. This provides evidence of the ability of our models to predict cellular responses despite the relatively small number of data points used to generate these models.
Effects of Single Signal Versus Multiple Signals
Although all four drugs inhibit cellular functions common in both cell types, combinations of these drugs may result in a significant difference in cellular ATPlevels between A549 and AG02603. Inhibition of A549 cells and preservation of AG02603 cells using the same drug combination constitute two conflicting objectives. Our goal of identifying the drug combinations that effectively satisfy the above objectives can be realized by utilizing a multiobjective search or optimization technique. A performance function combining the relative importance of each of the two conflicting objectives is introduced (Materials and Methods). The drug combinations resulting in the highest performance are considered the best drug combinations of a given set of drugs with corresponding concentrations.
Identification of the best performing subset is achievable through various methods. Enumeration and sorting of all possible combinations and their corresponding performances is one method. Alternatively, we can use a clustering algorithm such as a kmeans clustering algorithm to group combinations with similar performances [20,21].
Clustering the points into 20 different groups (Figure 7A), we find that the points with the highest performance are associated with low A549 ATPlevel and moderate to high AG02603 cellular ATPlevel (Figure 7B). The heat map on both panels is a function of performance and the best performing combinations are highlighted in dark red (Figure 7). The effect of drug combinations on the difference in the ATPlevel of A549 and AG02603 cells can be clearly seen in Figure 8 which shows the predicted individual doseresponse curves of both cell types in the presence of different levels of the other three drugs. The levels chosen are zero (low) concentrations, medium concentrations, maximum concentrations, and a selected combination from the subset of combinations identified to maximize the introduced performance function (Table 3). A significant difference in the response was observed at the selected combinations. This result illustrates how a properly chosen combination can result in a response unachievable individually by any of the drugs. As evident in Figure 8, the response to any individual drug on its own was small. However, a properly selected combination of the same four drugs yields a significant increase in the response.
Figure 8. Singledrug response curves. The figure shows an evaluation of the modelpredicted ATPlevels of A549 and AG02603 when the concentrations of three drugs are held constant while the concentration of the fourth drug is varied. The values of the concentration at which the three drugs are held constant are shown in Table 3. Shown are the singledrug effects of increasing concentrations of (A) AG490, (B) U0126, (C) indirubin3'monoxime and (D) GF 109203X on cellular ATPlevel with respect to different levels of the other three drugs in the treatment. Here the second row fourth column subfigure corresponds to the cellular ATPlevels of A549 and AG02603 to varying concentrations of U0126 while the remaining three other drugs are held at the selected concentration in Table 3.
Table 3. Shown are the concentrations used for for generating the drug response curves.
Additionally, we investigated the effects of pairwise combinations of drugs on the responses of both A549 and AG02603 cells. The remaining two drugs were fixed at one of two cases, zero or a selected optimal concentration (Table 3) from the set of combinations identified to maximize the performance functions. Similar to single drug responses, there is a significant difference between combinations of two drugs and combinations of 4 drugs (Figure 9). The data illustrates that there is a significant difference between normal and cancer cell ATP levels when two drugs are used with the other two fixed at a selected combination concentration versus zero concentration. In addition, the data shows that using a four drug combination increases the effective range or therapeutic windows of the two drugs when compared to two drug combinations.
Figure 9. Twodrug interactions. The figure shows an evaluation of the modelpredicted ATPlevels of A549 and AG02603 when the concentrations of two drugs are held constant while the concentration of the other two drugs are varied. The first and third columns (A) and (C) correspond to ATP levels of AG02603 in response to variations in different pairs of drug. The concentrations of the other two drugs are fixed at zero for column (A) and a selected combination (Table 3) for column (C). The second and fourth columns (B) and (D) correspond to ATP levels of A549 in response to variations in different pairs of drug (rows). The concentrations of the other two drugs are fixed at zero for column (B) and a selected combination (Table 3) for column (D). The figure shows that at the selected combination, the difference between the response surface of AG02603 and A549 is significantly higher than when the difference at zero concentrations.
Examining the Efficacy of Mixtures of Two, Three, and Four drugs
The availability of the model enables examining all lower order mixtures of drugs and their potential performances. Testing four or more drugs at time can enable efficient identification of effective lower order mixtures of the drugs. Based on the model, we evaluated the performances of all possible mixtures of drugs at varying concentrations. We identified the best achievable performance for each mixture (Figure 10). The results show that there is an improvement in the achievable performance as more drugs are used.
Figure 10. Evaluation of the best combinations of all lower order mixtures of the drugs. (A) The performance of the best combination for single drugs as well as combinations of two, three, four drugs are shown. Also shown are the performances of the worst possible combinations. (B) The normalized ATP level for the best combination for all single, two, three, and fourdrug mixtures are shown.
However, this improvement becomes less as four drugs are used instead of three. The results also show that the improvement is more significant for certain mixtures of drugs, e.g., the performance of the mixture of AG490 and I3M is better than the performance of AG490 and U0126. A similar observation can be made about threedrug mixtures where the mixture of AG490, I3M and GF109203X performs better than the three other threedrug mixtures. This information is not known a priori and without this approach, determining the best mixtures requires testing of each mixture independently. This can be a lengthy and costly process. Moreover, if one elects to randomly select combinations for the mixtures then there is no guarantee that the combination or mixture can have a good performance. In fact the random combination and mixture can perform quite poorly (Figure 10A). Hence, an effective approach for identifying mixtures of drugs and effective combinations can be very useful in this regard.
The potency of our approach can be illustrated by examining the performance of this approach versus more simple approaches. An example of such approaches is combining the drugs at the best singledrug concentrations. The best concentrations for each single drug were 30uM, 30uM, 30uM, and 100uM for AG490, U0126, I3M, and GF109203X respectively. Combining the drugs at the best singledrug concentrations results in a poor performance of 44.48 corresponding to zero normalized ATP levels of A549 and AG02603. One can argue that this combination can be very toxic because of the combined high concentration of each of the drugs. However, there is no simpler way to reduce the concentration of some or all drugs to achieve a better performance. Additionally, an argument can be made for the selection of the mixture of two drugs by picking the two best single drugs. In our case, this would correspond to mixing U0126 and I3M. This mixture does not perform as well as the mixture of AG490 and I3M (Figure 10). Again, there is no simpler approach for the selection of such mixtures. The approach we introduced attempts to answer to this need and provides a systematic approach that can be used to identify the best combinations for mixtures of two, three, and four drugs. In a more general setting, the approach enables identification of efficient combinations and mixtures for any number of ndrugs.
Validation of selected top performing combinations (30 uM of AG490, 0.3 uM of U0126, 10 uM of I3M, and 0.3 uM of GF109203X), (30 uM of AG490, 0.3 uM of U0126, 30 uM of I3M, and 1 uM of GF109203X), and (30 uM of AG490, 1 uM of U0126, 10 uM of I3M, and 0.3 uM of GF109203X) showed that the experimental ATP levels of A549 and AG02603 were 0.2 and 0.65, 0.13 and 0.51, and 0.13 and 0.5 respectively. This is in agreement with the predicted values for A549 and AG02603 of 0.21 and 0.59, 0.08 and 0.51, and 0.2 and 0.58 respectively.
DrugDrug Interactions Affect the Observed System Responses
Modeling of multisignal induced cellular responses can also be used to study key systemlevel signaling interactions between the applied signals and their effects on the system outputs. To that end, consider the linear regression model with interaction terms. This model uses a total of 15 regressors that are the concentrations of the drugs, six pairwise interaction terms, four threedrug interaction terms, and one fourdrug interaction term. This is a large number of regressors and, potentially, only a portion of these regressors are necessary to describe the variation in the data. Reduction of the regressors matrix X can be achieved using principal components analysis. Partial least squares can also be used and in this case the reduction in the regressors matrix takes into account variation in the output. However, the components obtained with principal component analysis or partial least squares models would be hard to interpret given the large number of regressors. Instead, we pursue a subset selection algorithm based on all the possible subset regressions [22]. The algorithm provides the best models of 1, 2, 3, . . ., 15 regressors. In total, the algorithm provides the best 15 models out of 2^{15 } 1 possible models.
The residual sum of squares of the best models shows that there is no significant reduction in the residual sum of squares for models with more than 10 regressors for the A549 data, and for models with more than 7 regressors for the AG02603 data (Figure 11A,B). The regressors used in these models correspond to singledrug concentrations, pairwise interactions, and three drug interactions for the the A549 data. Whereas the regressors for the AG02603 data included only singledrug concentrations and pairwise interactions (Figure 11C,D). The effects of these interactions vary between positive and negative. The individual concentrations and threedrug interactions have a negative influence, and the pairwise interactions have a positive influence on ATP levels. Examining the interaction terms, these interaction terms can be grouped into three categories, interactions that only occur in A549 cells, interactions that only occur in AG02603 cells, and interactions that are common to both cell types (Figure 12).
Figure 11. Analysis of models with varying numbers of regressors. Shown are the results of best subset regression algorithm which evaluates the best models of 1, 2, . . ., n regressors. (A) The residual sum of squares of the best models of A549 data as a function of the number of components (regressors). There is no significant decrease in the residual sum of squares for models with more than 10 regressors. (B) The residual sum of squares of the best models of AG02603 data as a function of the number of components (regressors). There is no significant decrease in the residual sum of squares for models with more than 7 regressors. (C) The regression coefficients for the best 10 component model of A549 data. The components involve single drug concentrations, pairwise and threedrug interactions. (D) The regression coefficients for the best 7 component model of AG02603 data. The components involve single drug concentrations and pairwise interactions.
Figure 12. Systemlevel signaling interactions in A549 and AG02603 cells. Shown are some modelpredicted systemlevel signaling interactions upon treatment with multiple drugs. The figure shows the interactions that are specific to A549, the interactions that are specific to AG02603, and the interactions common to both A549 and AG02603 cells. Direct arrows between the drugs and the ATP levels exist and are omitted to simplify the diagram.
Discussion and Conclusions
We demonstrated the integration of an efficient mathematical approach for a systematic quantitative characterization of the effect of multisignal combinations in two different cell types. Our method enables the establishment of accurate models to directly connect multisignal combinations and their effects through a learning process. Only a small percentage of total data points are required to be experimentally tested to establish a predictive model that is capable of simulating the effect of all possible signal combinations.
The resulting predictive model is also able to systematically reveal the interdrug interactions which are often nonlinear relationships. Such a model is necessary for multidrug combination optimization as the optimal combination, upon the addition of a new drug, cannot be achieved simply by testing various amounts of the new drug added to the previously optimized combination. Moreover, the approach allows for examining all lower order mixtures of the drugs and for evaluating the effectiveness without added experimental overhead. This enables efficient selection of drug combinations of different drugs particularly as the interactions between the different drugs vary. The approach enables systematic selection of input signals (drug combinations) that can achieve desired therapeutic goals. Experimental validation of selected top performing combinations presents additional evidence of the validity and utility of the approach introduced.
Our approach requires a relatively small number of experimental measurements. With the development of large scalehigh throughput measurement systems, our approach will be necessary and more efficient, e.g., when the number of inputs increases. Additionally, it will allow for integration into automated machines for testing and analysis of various biological systems. The addition of other cellular outputs can be another factor in favor of using this approach.
Mathematical Modeling Methods
We demonstrated the utility of various mathematical modeling approaches for the purposes of modeling and predicting multisignal induced cellular responses. Linear regression with models including interaction and quadratic terms were capable of producing powerful predictive models. The usefulness of linear regression methods and subset selection algorithms was also demonstrated as they enabled determination of key systemlevel signal interactions that result in the observed cellular responses. Alternately, neural network models performed generally better for fitting the data particularly when fewer data points were available to fit the models. Additionally, different structures of neural networks were almost equally powerful in producing the models. This present ample opportunity to exploit different structure of networks to mimic more realistic interactions within the cell given such a priori data. The choice of the best method for fitting the model is dependent on many factors. While in some cases simple linear regression can be used, in other cases nonlinear regression methods are necessary to get satisfactory models. Neural networks represent powerful computational models with great flexibility. Different types and configuration renders them useful for a wide range of problems including dynamic prediction using, e.g., recurrent neural networks. This presents an additional use for neural network to build on a wide range of applications in biological and medical research [2334].
It is important to note here that the novelty in this work is not this particular use of artificial neural networks. Rather, it is the use of a systemlevel modelbased approach to study the effects of a large number of signals on various cellular processes and to use that as a basis for selection of the retrospective optimal input signals. Our method resembles a general systems biology approach that can be utilized to address a broad range of biological questions. The method enables comparison of multiple system performances through modeling of their responses. The advantage of our method becomes critical when input signal combinations are characterized for the development of effective in vivo therapies, in which case the limited experimental scale imposes restrictions on the number of possible drug combinations to test. By largely reducing the number of tests required, our approach can greatly facilitate the development of clinically applicable treatments.
For the drugs and cell types chosen, our results showed that four drugs did not provide a significant improvement over three drugs. However, this is not a general result and the causes of this observation are not well known and pose an interesting problem to be examined in a separate study. A potential cause is that the signaling pathways affected by some of the drugs are saturated potentially due to sharing of a common target between two of the drugs. In other studies we are working on with different drugs and cells, four drugs result in an improvement over three drugs and five drugs also improve on four drugs.
Mechanistic Reasoning
Our modeling approach enabled identification of key systemlevel signaling interactions that contribute to the observed responses. These interactions highlighted potential differences in the signaling networks of different cell types as demonstrated. While some signaling interactions were common to the two cell types studied, there were other interactions that were specific to each cell type. Many of these interactions were known a priori. For example, the interactions between U0126 and GF109203X, I3M and GF109203X, AG490 and I3M, as well as AG490 and U0126 and GF109203X. There are other interactions that were not known a priori such as the interaction between U0126 and I3M. Thus the prediction of such interactions can guide experimental identification of the molecules mediating the interactions.
In the work presented, we used mathematical tools to construct a predictive model of cellular outcomes. The method was developed in experimental systems that only involved single output measurements. However, the method is a general method that can be used to construct models of the effects of multiple signals on various cellular outcomes, including signaling intermediates. Molecules from various cellular pathways, "intermediates", can serve as candidates for measurement and quantification [1,35,36]. Moreover, measuring these cellular outcomes and intermediates at various time points enables the construction of predictive dynamic models. The incorporation of time provides the model with predictive capability similar to those of ordinary differential equation models, though without a priori assumptions or knowledge about the molecular interactions. The introduction of more cellular outcomes presents the opportunity to utilize additional tools that can infer sets of rules which can provide, at least in part, descriptive reasoning of some of the internal interactions within the cell [3739].
Implications on the Design of Combination Therapies
Combinational therapies have attracted significant research efforts. Combination therapy for cancer is one example. The challenges associated with cancer treatment range from a lack of understanding of fundamental pathogenic mechanisms to practical experimentation limitations. Cancer is usually caused by multiple mutations and alterations of multiple signaling pathways which pose an extra challenge when defining the mechanisms underlining cancer development [40]. In addition, there is extensive heterogeneity of tumors among individual patients. Thus there are multiple potential targets for cancer therapy, resulting in an increasing need for distinct therapeutic agents and their combinations. Furthermore, drug resistance, arising from a subpopulation of original tumor cells or from subsequent mutations under the selection pressure of drug treatment, has been a major hurdle in cancer treatment. Drug combination therapy is emerging as a potentially effective approach to prevent drug resistance, as well as achieve higher efficacy and lower toxicity [4143].
However, the design of combination therapies based on studies of the response to individual drugs might not lead to the desired outcomes as interactions between the various drugs can lead to unknown outcomes. The emerging toxicities can be a major hurdle to the development of combination therapies. Therefore, it is important to consider system level signals or outputs that are representative of the state of the cell in addition to cellular intermediates targeted by the drugs. As such, it is important to investigate whether and how a combination of drugs can lead to a better system response compared to any of the individual singledrug responses.
Development of drug combination therapies for cancer can lead to more effective therapies to overcome drug resistance and achieve maximal drug efficacy. The demanding task is to find a combination that is maximally effective on cancer cells with minimal effect on noncancer cells. The method introduced in this work provides the ability to study the systemlevel effects of combinations of drugs and use the data to select combinations that can lead to desired outcomes by comparing the mathematical models of multiple cellular types. Furthermore, the method is capable of studying and characterizing experimental problems of higher complexity such as the order and timing of administration of multiple inputs into the system, which may help further reduce cytotoxicity and enhance efficacy with the same drugs utilized in clinical treatment of diseases. Moreover, employing our method with emerging technologies such as microfluidic devices or "labonchip" will enable highthroughput investigation of multidrug combinations, and become a promising platform for developing personalized medicine.
Methods
Cell Culture and Experimental Measurements
A549 cells were cultured in RPMI1640 medium (CellGro) supplemented with 10% heatinactivated FBS, 100U/ml penicillin and 100g/ml streptomycin (CellGro). AG02603 cells were in MEM (Gibco) supplemented with nonessential amino acid (Gibco), 15%heatinactivated FBS, 100U/ml penicillin, and 100g/ml streptomycin. Cells were seeded one day before drug treatment at a density of 1000cells/well for A549, and 8000cells/well for AG02603, in 96 well plates. ATP assays were conducted with the ATPlite 1step assay system (PerkinElmer) following the manufacturer's instructions and the luminescence signal was measured by Lmax microplate luminometer (Molecular Devices).
Neural Networks Models
We used two types of neural networks to fit the data. The first is a single layer fourneuron multilayer perceptron and the other is a twolayer cascaded neural network with a single neuron per layer. The networks were constructed and trained using the neural network toolbox of MATLAB [44]. The training was done using Bayesian Regularization.
Linear Regression Models
We used two linear regression models to fit the data. The two models are
Performance Function
In our work, we used a single performance function reflecting the relative importance of each criteria to our setup. The criteria are: maximize AG02603 cellular ATPlevel and minimize AG02603 cellular ATPlevel. The performance function used is defined as follows. Let x_{nc }be the ATPlevel of AG02603 and x_{cc }be the ATPlevel of A549 cells. Then, the performance function used is given by Perf(x_{nc}, x_{cc})
where
μ is the column vector . The performance function range is adjusted to be within the interval [0,100] using the relation Perf = (Perf  0.25)/0.14 * 100. This function can be easily generalized to accommodate more decision criteria.
Authors' contributions
IA, FY, and RS wrote the paper. IA, FY, and RS analyzed the data. IA, CMH, and JSS devised the mathematical solutions. FY, JF, and KY conducted the wet experiments. IA carried out the computational experiments. IA, FY, SD, CMH, JSS, and RS conceived the idea. All authors read and approved the final manuscript.
Acknowledgements
We thank Shiva Portonovo and Drs. Christopher Denny, Desmond Smith, Helen Brown, Elliot Landaw, and Henry Huang for critical reading of the manuscript; all Sun lab members for helpful discussions. This study is supported by NIH (CA091791 and CA127042), the Stop Cancer Foundation, Burroughs Wellcome Fund, NSF grant (ECS0501394) and NIH Roadmap for Medical Research (PN2 EY018228).
References

Janes KA, Albeck JG, Gaudet S, Sorger PK, Lauffenburger DA, Yaffe MB: A systems model of signaling identifies a molecular basis set for cytokineinduced apoptosis.
Science 2005, 310(5754):16461653. PubMed Abstract  Publisher Full Text

Janes KA, Yaffe MB: Datadriven modelling of signaltransduction networks.
Nat Rev Mol Cell Biol 2006, 7(11):820828. PubMed Abstract  Publisher Full Text

Irish JM, Kotecha N, Nolan GP: Mapping normal and cancer cell signalling networks: towards singlecell proteomics.
Nat Rev Cancer 2006, 6(2):146155. PubMed Abstract  Publisher Full Text

Natarajan M, Lin KM, Hsueh RC, Sternweis PC, Ranganathan R: A global analysis of crosstalk in a mammalian cellular signalling network.
Nat Cell Biol 2006, 8(6):571580. PubMed Abstract  Publisher Full Text

Berenbaum M: Direct search methods in the optimisation of cancer chemotherapy regimens.
Br J Cancer 1990, 61:101109. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Curatolo M, Sveticic G: Drug combinations in pain treatment: a review of the published evidence and a method for finding the optimal combination.
Best Practice & Research Clinical Anaesthesiology 2002, 16(4):507519. PubMed Abstract  Publisher 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(13):510510. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sun CP, Usui T, Yu F, AlShyoukh I, Shamma J, Sun R, Ho CM: Integrative systems control approach for reactivating Kaposi's sarcomaassociated herpesvirus (KSHV) with combinatory drugs.
Integrative Biology 2009, 1:123130. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Calzolari D, Bruschi S, Coquin L, Schofield J, Feala J, Reed J, McCulloch A, Paternostro G: Search Algorithms as a Framework for the Optimization of Drug Combinations.
PLOS Computational Biology 2008, 4(12):e1000249. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Iadevaia S, Lu Y, Morales FC, Mills GB, Ram PT: Identification of Optimal Drug Combinations Targeting Cellular Networks: Integrating PhosphoProteomics and Computational Network Analysis.
Cancer Research 2010, 70(17):67046714. PubMed Abstract  Publisher Full Text

Feala J, Cortes J, Duxbury P, Piermarocchi C, McCulloch A, Paternostro G: Systems approaches and algorithms for discovery of combinatorial therapies.

Mellor D, Prieto E, Mathieson L, Moscato P: A Kernelisation Approach for Multiple dHitting Set and Its Application in Optimal MultiDrug Therapeutic Combinations.
PLoS ONE 2010, 5(10):e13055. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zinner RG, Barrett BL, Popova E, Damien P, Volgin AY, Gelovani JG, Lotan R, Tran HT, Pisano C, Mills GB, Mao L, Hong WK, Lippman SM, Miller JH: Algorithmic guided screening of drug combinations of arbitrary size for activity against cancer cells.
Molecular Cancer Therapeutics 2009, 8(3):521532. PubMed Abstract  Publisher Full Text

Mukherjee JS, Rich ML, Socci AR, Joseph JK, Virú FA, Shin SS, Furin JJ, Becerra MC, Barry DJ, Kim JY, Bayona J, Farmer P, Fawzi MCS, Seung KJ: Programmes and principles in treatment of multidrugresistant tuberculosis.
Lancet 2004, 363(9407):474481. PubMed Abstract  Publisher Full Text

Nature 2001, 410(6831):9951001. PubMed Abstract  Publisher Full Text

Evan GI, Vousden KH: Proliferation, cell cycle and apoptosis in cancer.
Nature 2001, 411(6835):342348. PubMed Abstract  Publisher Full Text

Bellamy WT: Prediction of response to drug therapy of cancer. A review of in vitro assays.
Drugs 1992, 44(5):690708. PubMed Abstract  Publisher Full Text

Jain AK, Mao J, Mohiuddin KM: Artificial Neural Networks: A Tutorial.

Gupta MM, Homma N, Jin L: Static and Dynamic Neural Networks: From Fundamentals to Advanced Theory. New York, NY, USA: John Wiley & Sons, Inc; 2003.

Afifi A, Clark VA, May S: ComputerAided Multivariate Analysis. fourth edition. Chapman & Hall/CRC; 2003.

Hofmann M, Gatu C, Kontoghiorghes E: Efficient algorithms for computing the best subset regression models for largescale... [http://linkinghub.elsevier.com/retrieve/pii/S0167947307001168] webcite

AgatonovicKustrin S, Beresford R: Basic concepts of artificial neural network (ANN) modeling and its application in pharmaceutical research.
J Pharm Biomed Anal 2000, 22(5):71727. PubMed Abstract  Publisher Full Text

Bourquin J, Schmidli H, van Hoogevest P, Leuenberger H: Application of artificial neural networks (ANN) in the development of solid dosage forms.
Pharm Dev Technol 1997, 2(2):11121. PubMed Abstract  Publisher Full Text

Bourquin J, Schmidli H, van Hoogevest P, Leuenberger H: Basic concepts of artificial neural networks (ANN) modeling in the application to pharmaceutical development.
Pharm Dev Technol 1997, 2(2):95109. PubMed Abstract  Publisher Full Text

Chen HY, Chen TC, Min DI, Fischer GW, Wu YM: Prediction of tacrolimus blood levels by using the neural network with genetic algorithm in liver transplantation patients.
Ther Drug Monit 1999, 21:506. PubMed Abstract  Publisher Full Text

de Matas M, Shao Q, Richardson CH, Chrystyn H: Evaluation of in vitro in vivo correlations for dry powder inhaler delivery using artificial neural networks.
Eur J Pharm Sci 2008, 33:8090. PubMed Abstract  Publisher Full Text

Lin CC, Wang YC, Chen JY, Liou YJ, Bai YM, Lai IC, Chen TT, Chiu HW, Li YC: Artificial neural network prediction of clozapine response with combined pharmacogenetic and clinical data.
Comput Methods Programs Biomed 2008, 91(2):919. PubMed Abstract  Publisher Full Text

Papadourakis GM, Gaga E, Vareltzis G, Bebis G: Use of artificial neural networks for clinical decisionmaking (Maldescensus testis).
Neural Networks, 1992. IJCNN. International Joint Conference on 1992, 3:159164.

Payne SJ, Arrol HP, Hunt SV, Young SP: Automated classification and analysis of the calcium response of single T lymphocytes using a neural network approach.
IEEE Trans Neural Netw 2005, 16(4):94958. PubMed Abstract  Publisher Full Text

Penny W, Frost D: Neural networks in clinical medicine.
Med Decis Making 1996, 16(4):38698. PubMed Abstract  Publisher Full Text

Sun Y, Peng Y, Chen Y, Shukla AJ: Application of artificial neural networks in the design of controlled release drug delivery systems.
Adv Drug Deliv Rev 2003, 55(9):120115. PubMed Abstract  Publisher Full Text

Xie H, Gan Y, Ma S, Gan L, Chen Q: Optimization and evaluation of timedependent tablets comprising an immediate and sustained release profile using artificial neural network.
Drug Dev Ind Pharm 2008, 34(4):36372. PubMed Abstract  Publisher Full Text

Yamamura S, Kawada K, Takehira R, Nishizawa K, Katayama S, Hirano M, Momose Y: Prediction of aminoglycoside response against methicillinresistant Staphylococcus aureus infection in burn patients by artificial neural network modeling.
Biomed Pharmacother 2008, 62:538. PubMed Abstract  Publisher Full Text

MillerJensen K, Janes KA, Brugge JS, Lauffenburger DA: Common effector processing mediates cellspecific responses to stimuli.
Nature 2007, 448(7153):604608. PubMed Abstract  Publisher Full Text

Irish JM, Hovland R, Krutzik PO, Perez OD, Bruserud Φystein, Gjertsen BT, Nolan" GP: Single cell profiling of potentiated phosphoprotein networks in cancer cells.
Cell 2004, 118(2):217228. PubMed Abstract  Publisher Full Text

Andrews R, Diederich J, Tickle AB: Survey and critique of techniques for extracting rules from trained neural networks.

Jacobsson H: Rule extraction from recurrent neural networks: A taxonomy and review.
Neural Computation 2005, 17:12231263. Publisher Full Text

Kahramanli H, Allahverdi N: Rule extraction from trained adaptive neural networks using artificial immune systems.
Expert Systems with Applications 2009, 36:15131522. Publisher Full Text

Sjoblom T, Jones S, Wood LD, Parsons DW, Lin J, Barber TD, Mandelker D, Leary RJ, Ptak J, Silliman N, Szabo S, Buckhaults P, Farrell C, Meeh P, Markowitz SD, Willis J, Dawson D, Willson JKV, Gazdar AF, Hartigan J, Wu L, Liu C, Parmigiani G, Park BH, Bachman KE, Papadopoulos N, Vogelstein B, Kinzler KW, Velculescu VE: The Consensus Coding Sequences of Human Breast and Colorectal Cancers.
Science 2006, 314(5797):268274. PubMed Abstract  Publisher Full Text

Vogelstein B, Kinzler KW: Cancer genes and the pathways they control.
Nat Med 2004, 10(8):789799. PubMed Abstract  Publisher Full Text

Sawyers CL: Cancer: Mixing cocktails.
Nature 2007, 449(7165):993996. PubMed Abstract  Publisher Full Text

Ramaswamy S: Rational design of cancerdrug combinations.
N Engl J Med 2007, 357(3):299300. PubMed Abstract  Publisher Full Text

Demuth H, Beale M, Hagan M: Neural Network Toolbox 6 User's Guide.