Chemical combination effects predict connectivity in biological systems

Chemical synergies can be novel probes of biological systems. Simulated response shapes depend on target connectivity in a pathway. Experiments with yeast and cancer cells confirm simulated effects. Profiles across many combinations yield target location information.


Introduction
Living organisms are built of interacting components, whose function and dysfunction can be described through dynamic network models (Davidson et al, 2002). Systems Biology involves the iterative construction of such models (Ideker et al, 2001), and may eventually improve the understanding of diseases using in silico simulations. Such simulations may eventually permit drugs to be prioritized for clinical trials, reducing potential risks and increasing the likelihood of successful outcomes.
Owing to the staggering complexity of biological systems, efforts to model them require large and diverse sets of data on connections between components and responses to system perturbations. Among the most advanced models are those developed for baker's yeast, Saccharomyces cerevisiae (Zhang et al, 2005), which are built upon protein-DNA  and protein-protein (Ho et al, 2002) associations, supplemented by correlated changes in gene expression (Hughes et al, 2000) or protein abundances (Gygi et al, 1999) under differing conditions. Information gleaned from targeted synergies, such as paired mutations (Tong et al, 2004) and gene-drug interactions , have proven to be especially useful for revealing functional connections between components. Chemical combinations also show promise, and a proliferation experiment with yeast mutants in the presence of probe mixtures (Haggarty et al, 2003) has found that chemical profiles correlate with genetic similarity. This potential is confirmed by recent experiments using antibacterial combinations (Yeh et al, 2006) that show a relationship between synergy and chemical target relatedness.
Combination responses to varying concentrations of compounds provide a more detailed look at synergistic perturbations. Combination therapies have been used increasingly over the past century, and comprehensive reviews (Berenbaum, 1989;Greco et al, 1995) describe the experimental designs and combination analyses employed. Combinations of two or more agents can be tested using either exhaustive or efficient designs (Carter and Wampler, 1986), and the most widely used is the factorial design (also 'checkerboard' or 'dose matrix') where combinations are tested in all possible permutations of serially diluted single agent doses (Figure 1). A dose-matrix experiment comprehensively samples the underlying response surface with few assumptions about its shape. We have previously reported an approach for high-throughput dose-matrix screening of chemical combinations (Borisy et al, 2003;Keith et al, 2005;Zimmermann et al, 2007) in cell-based assays that preserve disease-relevant biological connections. Such screens yield a variety of response surfaces, with distinct shapes for combinations that work through different known mechanisms, suggesting that combination effects may contain information on the nature of functional connections between drug targets.
The past study of drug combinations has focused mainly on the question of whether a combination is more potent than equally effective doses of its constituents (Greco et al, 1995).
Synergy over this level is especially important when justifying clinical uses, as it defines the point at which the combination can provide additional benefit over simply increasing the dose of either agent. This most widely used dose additivity model (Loewe, 1928) represents the expected response if both agents are actually the same compound. In that case, a slice through the response surface at any chosen iso-effect level (or 'isobole') should show a linear relationship between the doses of the two agents. For example, if 50% inhibition (I¼T/U for treated and untreated samples) is achieved separately by 1 mM of drug A or 2 mM of drug B, a combination of 0.5 mM of A and Figure 1 The morphology of cell-based dose-matrix responses to chemical combinations differs between mechanisms. Data for two synergistic antibacterial combinations using a methicillin-resistant Staphylococcus aureus strain are shown in three-dimensional projections. In Augmentin s (A), clavulanate disables a bacterial defence against penicillin drugs, whereas both agents in Bactrim s (B) inhibit enzymes in folate metabolism. The mechanisms of action differ between the two combinations, as do the shapes of their response surfaces. Dose-matrix response surfaces show the inhibition of growth (1Àtreated/untreated) for each pairwise permutation of the serially diluted single agent doses. The 2D maps (C, D) show a top-down view of the surface, using the same colours for inhibition levels, to better display mathematical descriptions of the shape without obscuring any data, and permit many matrices to be shown together with clarity.
1 mM of B should also inhibit by 50%. Formally, the response at combined concentrations X,Y is the inhibition I Loewe that satisfies (X/X I ) þ (Y/Y I )¼1, where X I and Y I are the single agent effective concentrations that produce I Loewe . Deviations from Loewe additivity are usually quantified using the combination index (Chou and Talalay, 1983) CI¼(X/X I ) þ (Y/Y I ), which is essentially a ratio of total effective drug dose (combination versus single agents) required to achieve a given effect level. Loewe additivity is experimentally demanding because determining the effective concentrations requires well-sampled single agent curves. Moreover, the dose-additive constraint is not in closed form; so Loewe additive response surfaces generally must be computed using iterative root finding (Berenbaum, 1985). The other frequently used reference model is Bliss independence (Bliss, 1939), which is the expectation (multiplicative probabilities) for independent yet competing elimination agents, like bullets aimed at a limited set of targets. For inhibitions, I Mult ¼I X þ I Y -I X I Y , where I X and I Y are the single agent inhibition levels at concentrations X and Y. Although this model is of questionable relevance to drugs in biological settings, Bliss independence has often been favoured because it can be directly calculated from minimally sampled experiments without single agent response curve interpolation or iterative root finding. A third useful reference is the highest single agent (HSA) model, or Gaddum's non-interaction (Berenbaum, 1989), where the expected combination effect is simply the maximum of the single agent responses at corresponding concentrations.
The relationship between combination effects and drug mechanism has been intuitively recognized for some time, but attempts to theoretically model this association (Ashford and Cobby, 1974;Harrap and Jackson, 1975;Jackson, 1993) have proven problematic. These studies simulated a branched metabolic pathway as a system of linked ordinary differential equations, and the final reaction velocity was calculated in response to varying concentrations of inhibitors targeting different modelled enzymes. Synergy over Loewe additivity was assessed for each pair of enzyme inhibitors by computing simulated dose matrices and reporting the smallest combination index at 50% inhibition. However, these efforts could not unambiguously predict synergy or antagonism, because it was possible to generate any combination index for most inhibitor pairs by varying one or more of the many kinetic parameters in the pathway model.
We approached the problem by focusing not just on the presence or absence of synergy relative to some reference level, but by using a set of shape models to characterize the morphology of full response surfaces. Here we show that there is a clear relationship between dose-matrix response shape for a combination and the connectivity of its targets (the type of connection between them), using simulated pathways and combination experiments. The simulations were used to define a reference set of shape models that represent connectivity-related responses, and the experiments demonstrated that the shapes predicted by the simulations do indeed occur in biological settings. We also use the data to show that the synergy profile produced by a drug across many combinations in a screen depends on the drug's mechanism, and can also be used to infer connectivity relationships between drug targets.

Results
We considered four models for combination effect morphology ( Figure 2 and Materials and methods) that reflect historical combination analyses and that represent many of the responses observed in our therapeutic combination screens. The HSA model describes a simple superposition of the single agent curves. Loewe additivity (Loewe, 1928) is the drug-withitself reference to represent dose-additive pairings. The 'Bliss boosting' model extends Bliss Independence (Bliss, 1939) to allow variable boosts in effect at high concentrations. Finally, the 'Potentiation' model describes combinations like Bactrim s (Figure 1), where one agent's curve is shifted with a power-law slope and superposed on the enhancer's own activity. These models are not an optimally designed basis set of shape descriptors because they are not orthogonal (they produce

Bliss boosting Potentiation
Combination response shape models that describe many of the observed response morphologies. Each model (shown using the same colour scale as Figures 1 and 3) is used to calculate an expected combination effect I model at any concentration X,Y, based on the single agent response curves. (A) HSA is a superposition of the X and Y single agent responses, calculated from the inhibitions I X at X and I Y at Y. (B) Loewe additivity (Loewe, 1928) is the drug-with-itself reference for synergy, where I Loewe at X,Y yields additive doses relative to the components' effective concentrations X I ,Y I at I Loewe . (C) Bliss boosting describes combinations with a variable boost b above E max (the greater of the single agent limiting efficacies E X ,E Y ), at high combined concentrations. Useful reference levels for Bliss boosting are 'cancelling', 'suppressing', 'masking', 'multiplicative' corresponding to Bliss independence (Bliss, 1939), and 'saturating' (see Materials and methods). Finally, (D) potentiation can characterize responses similar to those of Bactrim s (Figure 1), where one single agent's curve I X (C) is shifted with a power-law slope p above an enhancer concentration Y pot , and superposed on the enhancer's own activity. These models can be extended to higher-order combinations, and used in the same form (with adjustments to Bliss boosting) for any type of measurement, provided that the activities of both agents vary monotonically with concentration. Of these models, only Loewe additivity has an a priori mechanistic basis.
degenerate, or indistinguishable, shapes for some parameter settings), and also are not complete, in that some observed surfaces do not fit any of them. Nevertheless, this set does cover many of the observed effects in our therapeutic screens, and permits a quantitative exploration of different synergy types. Each observed dose matrix can be analysed by fitting all the shape models to the data by means of least-squares minimization across any free parameters, and determining which shape models produce consistent fits. It is also useful to calculate a volume V HSA between the data and the HSA surface, normalized to the concentration sampling, to characterize the overall strength of combination effects.

Simulations of metabolic pathways
We used numerical simulations of metabolic pathways to show that pairs of targets with differing connectivity produce distinct combination effects ( Figure 3). We constructed a branched, unregulated pathway of linked Michaelis-Menten reactions (Supplementary 1), and computed the inhibition of final reaction velocity after applying varying concentrations of paired competitive enzymatic inhibitors (affecting K m ). Each pair of inhibitors was applied repeatedly over a dosing matrix, and we determined the response morphology by fitting our shape models to each simulated dose matrix, and finding which shape provided the best match over the whole surface ( Figure 3). We found that when both inhibitors targeted the same enzyme, Loewe additive surfaces were produced as expected. When the inhibitors had separate targets within this pathway, the resulting surfaces were best fit by Bliss boosting, with differing levels depending on the target connectivity. These differences in combination effect were robust to plausible variations of the simulated pathway that did not affect the network topology. Random perturbations of up to two orders of magnitude in all the reaction parameters recovered the original best-fit shapes (Supplementary Figure  S2), suggesting that the response morphology is affected by only such extreme changes in kinetic parameters that alter the network topology. Similar Monte Carlo simulations perturbing specific parameter classes or individual enzymes showed no strong sensitivities to particular kinetic parameters or targets for this pathway (Supplementary 3). Converting the inhibitors to non-competitive (affecting V m ) kinetics likewise recovered the original best-fit shapes, altering mainly the single agent curves. An exception is that unbypassed serial combinations (e.g., R1xR6; Figure 3) yielded HSA-like responses for noncompetitive inhibitors (Supplementary 2). Adjusting the end point to something more like proliferation by means of an exponential transform (see Materials and methods) distorted the resulting boost levels, but still did not change which shape model was preferred. All of these factors could strongly affect the single agent curves (e.g., potency, steepness, and limiting efficacy) but the best-fit shape model for each simulated combination effect was very stable to perturbation.
Introducing regulation by negative feedback did have strong effects on the combination responses. Pathway B was added to investigate synergy across pathways and to predict effects for sterol biosynthesis. Whereas inhibitors with serially placed targets in an unregulated linear pathway produced multiplicative boosts, enclosing such targets in a negative feedback loop yielded strong potentiation (although like unregulated serial inhibitors, these yield weaker effects for non-competitive inhibitors). The synergies under negative feedback are not immediately intuitive (Supplementary Figure S3) and an analytical derivation of the simulated power-law potentiation for this model would be challenging, underscoring that simulation is the most practical means for deriving even the simplest emergent properties of complex systems.
Another change affecting topology is to alter the type of junctions in a pathway. For the branches in Pathway A, the reaction could proceed down either side to the end point; so parallel pairs of inhibitors across these junctions (e.g., R3 with R5) are functionally analogous to logical OR operations. Inhibiting across such OR junctions produced saturating Bliss responses in our simulations. Replacing such junctions with functional AND operations (here implemented by pairing targets across Pathway A and Pathway B, both required for R AB ) resulted in masking surfaces (close to HSA or Bliss masking) for combinations that straddled the junction. Such masking effects were produced even for cross-pathway combinations involving R11 and R12, which are bypassed in Pathway B to ensure clear distinctions between masking, multiplicative, and saturating boosts.
The best-fit shapes for all the combinations in both simulated pathways are summarized in Figure 4, using symbols whose size reflects the HSA volume and where the colour and symbol denote the best-fit model. The ambiguity between the shape models for each combination is illustrated in Supplementary Figure S4. In these simulations, the fitted shapes clearly correspond to distinct mechanistic configurations. It is notable that all the cross-pathway combinations produced either HSA surfaces or other masking surfaces with weaker combination effects (with |V HSA |B1, between the measured effect and the HSA model surface) than are seen within each of the pathways (|V HSA |B3-10).
These simulations suggest specific hypotheses in terms of our shape models for experimental combinations of chemical inhibitors of metabolic pathways. Agents sharing a primary target should, of course, produce Loewe additive responses. Branched unregulated pathways should produce various levels of Bliss boosting, depending on the target connectivity, and linear pathways regulated by negative feedback should show power-law potentiation. Combinations with targets that cross between such pathways are likely to produce masking effects (similar to the HSA model or Bliss masking).

Yeast combinations experiment
To test some of the hypotheses from our simulations, we undertook a chemical combination screen focused on the sterol metabolism in yeast. This pathway is very well characterized with known inhibitors at a number of enzymes (Wills et al, 2000), and is regulated by negative feedback (Gardner et al, 2001). We selected 10 antifungal drugs (Table I), six with targets along a linear section of the pathway and four targeting other cellular functions. Because sterol biosynthesis is much the same across yeast species (Wills et al, 2000), we selected Candida glabrata, a frequent focus of antifungal treatment that we had used in past therapeutic screens, and measured proliferation using a metabolic assay (Materials and methods). After determining single agent activities for all the drugs, we designed dose-matrix experiments with concentration samples centred on each drug's effective concentration at 50% effect (EC 50 ) when possible.
The dose-matrix responses for each combination are presented in Figure 5 and Supplementary 4. The combination effects are visible towards the top right-hand corner of each response surface, especially between the sterol pathway inhibitors. Some of the single agents produced variable responses, and this is most obvious for those drug-with-itself dose matrices along the diagonal of the grid which show only one side appearing to be active. The self-combination of Figure 3 Simulations of a multiply inhibited network yield distinct response shapes that depend on target connectivity. Substrates (black nodes) are metabolized through a series of Michaelis-Menten reactions (grey circles) from sources with constant reaction velocities to a limitless sink at the end-point. We calculated dose-matrix inhibitions of the end point velocity at R E to paired competitive inhibitors at various enzymes. The response surfaces for some combinations are shown, with the inhibitors indicated by joined markers. For Pathway A, same-target pairs produced Loewe additivity, and separated targets led to various Bliss boosts depending on where the targets were placed. Pathway B was added to investigate combination effects across pathways, with a bypass to ensure clear distinctions between cross-pathway boosting levels. Because both inputs are required for R AB , the junction is equivalent to a logical AND function. The negative feedback in Pathway B was introduced to represent processes like sterol biosynthesis in yeast. Inhibiting across pathways produced HSA-like masking effects, and inhibiting within the negative feedback loop yielded potentiation effects. rapamycin illustrates how a spurious synergy can arise from instability in the measured potency of a drug, and the apparent synergy of miconazole and itraconazole shows evidence of single agent instability-its inhibition errors are exceptionally large, indicating inconsistency between the replicate blocks.
Much more striking are the synergies between sterol inhibitors, which span a greater range of concentrations and which recapitulate the known antifungal synergy between statins and azoles (Lorenz and Parks, 1990).
The shape and magnitude of each combination effect is summarized by symbols whose size represents the total synergy (as measured by V HSA ) and whose colour denotes the best-fit shape model. Degeneracies between the models are measured by the comparative chi-squared (w 2 ) goodness-of-fit estimates (Supplementary 4 and Supplementary Figure S5). Although the data quality did not permit a clear distinction between the HSA model and Loewe additivity, the two most similar models in the set, most of the drug pairs sharing the same targets and agent-with-self combinations produced responses that were consistent with Loewe additivity. The exceptions all had relatively weak synergy, with only one, Itraconazole with Miconazole, exceeding the estimated systematic error level (V HSA 41; see Materials and methods). By contrast, all but one of the cross-target sterol combinations exceeded that level, and all produced surfaces that were best fit by potentiation. Across pathways, the combination effects were more variable, and strong synergies (V HSA 41; see Materials and methods) were relatively rare (5/28 versus 10/22 for same pathway), with an B8% Poisson probability (Press et al, 1997, y14.3.3).
The observed effects confirm predictions from our pathway simulations. We expected same-target combinations to produce Loewe additivity and cross-target sterol pairings to give rise to potentiation. Among the 11 same-target and 11 crosstarget sterol combinations tested, the prediction accuracy is only 54% if we insist on unambiguously correct shape classifications (the predicted model has w 2 min , and Dw 2 4w 2 min for the next best fit). However, almost all the failures are due to the degeneracy of the models for HSA-like combinations. If we count as a success any surface whose best fit was indistinguishable from the predicted shape, with Dw 2 ow 2 min , the prediction accuracy rises to 72% with an uncertainty of B22% due to the sample size, assuming Poisson statistics within each confusion matrix class. Across pathways, the relative diversity of combination effects makes sense, because their target connectivities are probably more varied, and the relative rarity of strong synergies accords with the Figure 4 The best-fit shape models for our simulated pathway ( Figure 3) consistently depend on target connectivity. Combination effects between competitive inhibitors are shown as symbols whose area is proportional to the HSA volume and whose colour indicates the best-fit model shape: black for HSA, blue for Loewe additivity, green for Bliss boosting, red for potentiation, and grey when two or more models could not be distinguished (see Supplementary Figure  S4). The Bliss boosting model is subdivided using shape to indicate the boosting level (square for saturating, diamond for multiplicative, circle for masking, and open circle for suppressing). The bold lines separate Pathway A from Pathway B combinations, and the thin lines partition Pathway A into its component sections. All inhibitor-with-self pairings fit Loewe additivity, as expected. Other combinations within Pathway A consistently fit the Bliss model with distinct boosting levels. Cross-target combinations within the negative feedback in Pathway B consistently show potentiation, and cross-pathway combinations always have masking responses (close to Bliss masking, additivity, or HSA).  Terrell and Hughes (1992). f Douglas et al (1997).
Chemical combinations and target connectivity J Lehár et al predominantly masking effects in our simulated cross-pathway combinations.

Combination screen in human tumour cells
There is further support for our simulations from a prior combination screen aimed at discovering potential combination therapeutics. Human tumour cells (HCT116) were tested using a proliferation assay for responses to all pairwise combinations of 90 drugs and probes, most with known mechanism (Supplementary 5). This screen provides a broader sampling of possible target connectivities than had our C. glabrata experiment, and permits us to investigate the applicability of our shape models to cellular signalling networks. Because this screen had been designed for discovering potent synergies, a four-fold dilution factor had been used (compared to two-fold for our yeast experiment), and many compounds were either inactive alone or had not been sampled far above their EC 50 concentrations, often obscuring distinctions between fitted shape models. Nevertheless, we expect strong potency shifts to be nearly always detectable, and other shape classes can be distinguished when the single agents have been appropriately sampled.
For each combination, we obtained an HSA volume V HSA and best-fit shapes as described above. We also examined 'synergy profiles' constructed from each agent's V HSA values across all the other probes, and recorded the correlation coefficient between the synergy profiles for each pair of drugs in the screen. Each pair of drugs was classified by target relatedness, using the known target annotations (Supplementary 5), as 'Identical' for drug-with-itself, 'Same' for drugs sharing a target, 'Related' for pairs with distinct targets in the same pathway or function, 'Different' for pairs targeting distinct functions, and 'Unknown' for pairs involving at least one uncharacterized drug. Inter-kinase combinations were classified separately as 'Kinases'. Supplementary 6 presents the resulting scores and correlations for all the tested combinations, and Supplementary 7 shows their distributions grouped by target similarity.
The shape models in Figure 2 provide reasonable coverage of the observed combination effects. Only B30% of 4092 drug pairs in the screen produced a formally acceptable fit (w 2 min o2) to any of the models (Supplementary 7a). However, given the systematic errors produced by the experimental process (Materials and methods), a more fair assessment of the uncertainties can be obtained by examining the w 2 distribution of the Loewe additive model for the 87 drug-with-self combinations, as a negative control set. These 'Identical' combinations have a median additive w 2 B9, suggesting that the uncertainties are roughly three times larger than the formal estimates from replicate dose matrices, due to systematic errors that were not accounted for by the scatter between replicates. Adopting this w 2 level as more representative of the errors extends the coverage (w 2 min o18) of our models to most of the  (Table I), six with known targets on the sterol pathway (inhibitor markers). The results are shown (centre) with concentrations increasing from the bottom left of each drug pair's dose matrix. The combination effect symbols (right) summarize the observed response for each pair in terms of the models presented in Figure 2. Size shows synergy as measured by V HSA , the volume between the data and HSA surface. The best-fit shape model is indicated by colour: black for HSA, blue for Loewe additivity, green for Bliss boosting, red for potentiation, and grey when two or more models could not be distinguished (see Supplementary Figure S5). Drug pairs sharing a target mostly produced weak responses (V HSA o1) consistent with Loewe additivity, and cross-target sterol combinations produced potentiation, as expected for a linear pathway under negative feedback (B70% prediction accuracy). Combinations involving non-sterol drugs showed fewer strong combination effects and their synergy profiles across all the drugs were more variable.
Chemical combinations and target connectivity J Lehár et al screen, leaving only B10%, which are dominated by those combinations that are equally inconsistent with all four models. For an analysis of the observed shape distributions, we selected a set of 'optimally sampled' combinations, where at least two concentration samples had been tested above the EC 50 for both single agents, enabling distinctions between strong combination effect shapes. Among these 253 combinations, B60% fit masking models (HSA, additive, or Bliss masking), B30% were best fit by potentiation, with the remaining B10% accounted for by multiplicative or saturating Bliss boosts. These proportions hold even for the 91 optimal combinations with sub-saturated single agent activities (E max o90%), for which masking and multiplicative effects would be clearly distinguishable. For all target connectivity groupings, the optimally sampled combinations were dominated by masking and potentiation effects.
Across the whole screen, combination effects and synergy profile correlations ( Figure 6A and B) support the trends observed in our pathway simulations and the C. glabrata experiment. When we examined synergy score distributions (Supplementary 7), we found that drug pairs targeting different pathways or involving unknown targets had significantly more variable scores (Materials and methods) with a slight trend to lower synergy levels. The synergy profile correlations show a more striking trend. Pairs with the same or related targets ought to have similar connectivities with the rest of the set, because they are likely to have similar mechanistic relationships to the rest of the drug targets in the screen. Indeed there are significant differences between the profile correlation distributions with a stronger trend towards increasing correlation with target similarity.
The responses to pairings of kinase inhibitors also provide support for an association between combination effects and connectivity. Whereas combinations between kinases in general had similar distributions to those for other related target pairs, combinations between drugs targeting the phosphatidylinositol-3 kinase (PI3K), protein kinase C (PKC), or mitogen-activated protein kinase (MAPK) signalling pathways (Table II) showed exceptionally strong synergies and correlations. HCT116 cells maintain an activating KRAS mutation that drives cellular growth and survival (Shirasawa et al, 1993), and are especially reliant upon epidermal growth factor receptor (EGFR)-mediated signalling for proliferation (Awwad et al, 2003). Because these pathways are central to EGFR signalling and interconnected (Schlessinger, 2004), we would expect HCT116 cells to be very sensitive to these combinations, and that inhibitors of those pathways should have similar synergy profiles. Moreover, the combination effects between the pathways show differences ( Figure 6C) that are consistent across mechanistic replicates, despite limited concentration sampling and evidence of secondary target activity. Such differences may eventually provide novel insights into how these pathways interact in the context of HCT116 proliferation (e.g., that the MAPK/PKC connection may be less direct than the other target pairings).

Discussion
The simulations and experiments presented here add a new dimension to earlier studies investigating the effects of multiple perturbations to biological systems. Early studies of statistical epistatic associations between random mutagenesis Figure 6 Combination effect statistics from the HCT116 screen show expected trends, and the exceptionally strong synergies between inhibitors of EGFR signalling (Table II) reflect known connections. The distributions of HSA volumes (A) and synergy profile correlations (B) are shown for drug pairs in the HCT116 tumour cell screen (Supplementary 6), grouped by target similarity ('Identical' for drug-with-self combinations; 'Same' for distinct drugs sharing a target; 'Related' for differing targets sharing a common functional group; 'Different' for pairs with different functions; and a separate group 'Kinases' for combinations between kinase inhibitors). There is a significant trend towards more variable synergy and lower profile correlations for more distantly related targets (Supplementary 7), and the most synergistic and highly correlated pairs within the 'Kinases' set are dominated by those targeting the EGFR signalling pathways, which are major drivers of HCT116 proliferation. (C) Combinations between MAPK/PI3K, MAPK/PKC, and PI3K/PKC targets show different levels of response that are consistent across mechanistic replicates, despite limited concentration sampling (only TUBC and ROTT were optimal) and evidence of off-target effects (synergy between PI3K inhibitors and discrepant ROTT effects compared to the other PKC inhibitors). Shape model degeneracies are shown in Supplementary Figure S6. and chemical stresses (Elena and Lenski, 1997;Kishony and Leibler 2003) were followed by systematic studies of mutant libraries in combination with chemical probes (Giaever et al, 2004;Lum et al, 2004;Parsons et al, 2004) and synergies between targeted chemical probes (Haggarty et al, 2003;Yeh et al, 2006). This study explores the additional information available from variable dose response matrices, and relates quantitative combination effects, as measured response shape classification or by HSA volumes, to target connectivity. We use the term 'connectivity' to capture both the distance and topology of the connection between the targets, rather than the more common 'interaction' , denoting the presence of a functional connection. Our pathway simulations show that biological systems are expected to produce a useful variety of combination effects that are correlated with connectivity, and our experimental results demonstrate that the simulations are biologically relevant and predictive.
There was only partial consistency between specific chemical synergies in our yeast experiment and corresponding interactions observed in earlier genetic screens (Supplementary 8). The observed synergies between sterol inhibitors are reflected in homozygous synthetic lethalities (Stark et al, 2006) that are internal to the pathway, and rapamycin's synergy with terbinafine agrees with the reported interaction between TOR1 and ERG1. However, the mutant screen of caspofungin's target, FKS1, showed no correlation with our combination effects. Moreover, there is little correspondence between our results and the cross-activities identified in drug sensitivity screens of heterozygous yeast mutants (Giaever et al, 2004;Lum et al, 2004;Parsons et al, 2004). The mixed results of these comparisons are not very surprising. Aside from experimental errors (see Materials and methods), some of the disagreement can be ascribed to the quantitative nature of dose-matrix response data, because the responses to varying doses (Hartman and Tippery, 2004), or even timeresolved experiments (Warringer et al, 2003), can detect the effects of non-essential genes that would have been missed in lethality screens. Another contribution may come from likely differences between genetic and protein networks (Ozier et al, 2003). Chemical combinations can interfere directly with protein interactions without mediation through gene expres-sion, and although the protein and gene networks are intertwined, connectivity differences between the two types of targets should be expected.
There was more agreement between our simulation results and a prior flux balance analysis (FBA) of the yeast metabolic network (Segrè et al, 2005) in response to paired mutations. The simulated mutations should correspond roughly to the high-concentration limit of our dose matrices. Both approaches predict strong synergies (synthetic lethality) when the system is inhibited across parallel alternative pathways, and find masking effects (buffering) when unrelated pathways are impaired, provided that the two pathways are not competing for a common supply of precursors (Supplementary Figure S2 in Segrè et al). The main difference is that the FBA simulations produced more multiplicative than masking pairings, especially for cross-pathway combinations. This may be due to the partial view of combination effects provided by the epistasis scores (Supplementary Figure S1), or to an implicit assumption, in the FBA simulations, that competition for the same metabolic ingredients is widespread. The masking responses found in our cross-pathway simulations are supported experimentally, by the weaker synergy levels seen in our C. glabrata experiment, and by the large fraction of masking effects (B60%) in both the HCT116 screen and prior antibacterial combination experiments (Yeh et al, 2006).
Thus the dose-matrix response shapes and synergy profiles from systematic combination screens can provide valuable new insights into network connectivity. A possible application in chemical genetics would be to extend phenotypic profiling for mechanism inference (Perlman et al, 2004;Macdonald et al, 2006), using combinations effectively as extra assays. Query drugs with unknown targets could have their synergy profiles compared to a previously assembled profile library for a diverse set of chemical probes. The known targets for those probes whose profiles correlate most closely with a query drug's are likely candidates for the drug's mechanism. For constructing biological network models, observed combination effects can be compared to expectations from target connectivity through existing protein interaction networks, and inconsistencies can be used to improve the models by means of a prediction-validation procedure (Ideker et al, 2001). Combining chemicals is much simpler than the process of constructing double mutant strains, and chemical combinations can yield detailed information about connections between proteins that are not accessible to single agents or mutagenesis, enabling the study of disease networks like those involved in cancer (Schoeberl et al, 2002) and inflammation (Bouwmeester et al, 2004) signalling.
Of course, there are limitations to this approach for probing biological networks. The collection of full dose matrices is expensive, and will necessitate a trade-off between the added information obtained from response shapes and the loss of assay fidelity from the use of high-throughput measurements (e.g., single time points for growth rather than time-resolved growth curves). The mechanistic information that can be extracted will depend strongly on the generality and specificity of the set of shape models being used. The four models presented here represent only a subset of the observed responses, and the pathway simulations sampled only a small minority of possible target connectivities. Even with more  Gschwendt et al (1994).
Chemical combinations and target connectivity J Lehár et al extensive simulations, there will certainly be mechanistic degeneracies associated with each distinct response shape. For example, it will probably be possible to produce potentiation like that seen in Bactrim using several distinct target configurations. Thus, response shapes cannot be used to uniquely define combination mechanisms, but merely to exclude inconsistent connections. Finally, chemical probes often have poorly characterized mechanisms or multiple targets, and combination screens will have to account for this by covering biological targets with several chemically distinct probes. Despite all these limitations, the present study has shown that even with a limited set of shape models, dosematrix data can provide useful connectivity information.
The results of this initial exploration suggest several directions for expanding the response surface approach to combination analysis. First, it is essential to extend the current shape model set by using more comprehensive simulations of larger networks. The most obvious way to do this would be to define a basis set that represents simple observed network connectivities, including serial and parallel arrangements, bifurcations, both AND and OR junctions, and some regulation mechanisms (Milo et al, 2002). Combination responses from pathway simulations covering the connectivity basis could then be used to define shape models for comparison with simulated responses for larger networks, as well as with experimental data. Another approach would be to systematically explore perturbations of an existing biological network, such as S. cerevisiae metabolism (Famili et al, 2003), and to model drug inhibition by partially restricting the reaction rates at metabolic enzymes, essentially the Segrè et al's FBA for partial inhibitors. The challenge with this approach would be to identify meaningful connectivity measures that correlate with simulated synergy types. A third very promising area for future exploration is higher order combinations. The emergence of resistance to treatment for afflictions like cancer (Fojo and Bates, 2003) can be effectively addressed by combinations of three or more drugs (Komarova and Wodarz, 2005). Combination approaches can also offer improved control by perturbing multiple components of the disease network (Csermely et al, 2004), and higher order combinations should provide yet more details on mechanistic connections between their targets. The response surface approach presented here can be readily adapted to three or more drugs, as can the specific shape models shown in Figure 2, and a metabolic FBA of higher order combinations should provide important insights into the effectiveness of multi-target treatments. Finally, because the simulated combination responses were insensitive to the reaction details, this methodology may even be extendable to non-biological network problems where responses to targeted perturbations can be measured (Milo et al, 2002).

Combination response shape models
Experimental dose-matrix responses include both potency shifts, where a combination alters the apparent potency of the single agents, and efficacy boosts, where the joint effect exceeds levels possible for either of the constituents. Many of these responses can be described by simple models that use the single agent curves to predict the combination effect (Figure 2), and whose parameters provide quantitative measures. Strong potency shifts relative to a simple HSA surface can be described using a power-law potentiation model, and efficacy boosts are quantifiable using a model that extends Bliss Independence (Bliss, 1939) to variable combination effect levels at high concentration. These models are expressed in terms of inhibitions, but can be adapted to other types of measurement, or for combinations involving more than two perturbing agents, without altering their qualitative properties. The shapes presented are not a complete orthogonal set covering all expected combination responses, but were chosen to represent some of the variety seen in our therapeutic screens, as well as to reflect historical reference models used in the literature.
The HSA model, or Gaddum's non-interaction (Berenbaum, 1989), is based simply on the intuition that if a combination's effect exceeds those of its constituents, there must be some interaction. Mathematically, the HSA model is a superposition of the single agent curves where, at any combined concentration (X,Y), the inhibition I HSA ¼max(I X ,I Y ), where I X and I Y are effects produced by the single agents at (X,0) and (0,Y). The model values are calculated at each dosematrix point, where I X and I Y were determined using sigmoidal fits to the single agent response data. The volume V HSA ¼S X,Y ln f X ln f Y (I data ÀI HSA ) between the data and the HSA surface, adjusted for variable dilution factors f X , f Y , can be readily calculated and provides a practical overall measure of synergy (in units of inhibition), summed over dimensionless log-concentration space.
We also considered Loewe additivity (Loewe, 1928), the drug-withitself standard. At each combined concentration (X,Y), we used an iterative approach (Berenbaum, 1985) to find the inhibition I Loewe that satisfies (X/X I ) þ (Y/Y I )¼1, where X I and Y I are the single agent effective concentrations. Starting with a guess that I¼I HSA , we interpolated the single agent curves to find X I ,Y I that produce I, calculated the corresponding combination index, and used bisection (Press et al, 1997) to converge on a value of I with combination index CI¼1. Loewe additivity reduces to HSA at concentrations that are very different from X I ,Y I , as well as for I above an agent's limiting activity (effectively X I ,Y I -N).
To model boosts in efficacy at high concentrations different from what the single agents can achieve, we used a 'Bliss boosting' model, adapting the Bliss independence model (Bliss, 1939) for boost levels other than multiplicative. Mathematically, I Bliss ¼I X þ I Y þ (bÀE min )(I X I Y /E X E Y ), where E min is the lesser of E X and E Y , the limiting single agent efficacies. The one free parameter b, in units of effect, determines the amount of boosting above E max , the greater of the single agent efficacies. For inhibition effects, there are a number of useful reference levels: b¼ÀE max produces 'cancelling', with zero effect at high concentration; b¼E min ÀE max is 'suppressing', where the less effective agent prevails; b¼0 yields 'masking', where the more effective agent prevails; b¼E min (1ÀE max ) is multiplicative (Bliss independence); and b¼(1ÀE max ) produces 'saturating' at 100% inhibition. At the high concentration limit, the suppressing, masking, multiplicative, and saturating levels correspond to 'suppression', 'buffering', 'no epistasis', and 'synthetic lethal' in recent classifications for epistasis (Elena and Lenski, 1997;Segrè et al, 2005;Yeh et al, 2006). Bliss boosting can be used in this form for any type of effect measure provided that the single agent responses increase monotonically with concentration. Only the saturating and multiplicative reference levels need to be adjusted according to the measurement scale. For example, if we consider fitness ratios (treated over untreated growth ratios) with f¼ln(T/U) in place of inhibitions, there is no saturating level and b¼E min indicates multiplicative. Note that Bliss masking is similar to both HSA and Loewe additivity for very high or very low concentrations where the single agent effects are almost flat.
Finally, we used a power-law 'potentiation' model to describe strong shifts in potency similar to that seen for Bactrim s (Figure 1). For such combinations, an active single agent's response curve shows an increase in potency as the enhancer is titrated in, which is seen as linear iso-effect contours in serially diluted dose matrices (logarithmic concentration space). The inhibition I Potent ¼max(I X (C), I Y ), where I X (C) is the single agent response curve of the potentiated compound, at a shifted concentration C¼X[1 þ (Y/Y pot ) |p| ] sign(p) , where sign(p) is a unit sign function evaluated at (À1, 0, þ 1) corresponding to the sign of its argument. Although the functional form presented here is awkward, it was the simplest way we could find to achieve power-law potentiation and depotentiation above a threshold enhancer concentration, with the enhancer's own activity superimposed. The two free parameters for this model are the threshold concentration Y pot above which potentiation takes effect, and the potentiation slope p (synergy for positive and antagonism for negative p). There is no potentiation for p¼0, where this model reduces to an HSA surface. Just as for HSA and Loewe additivity, the form of this model is identical for any type of effect measure.
These models are not orthogonal; so we defined a fitting order to ensure an unambiguous choice. The shape models have degeneracies where two or more models can describe similar response shapes. Loewe additivity and HSA, for example, are largely degenerate, in that only when the single agents increase slowly over the sampled concentrations can a distinction between the two models be made. Nevertheless, both can be readily distinguished from Bliss boosting and potentiation, which themselves are distinct over a large range of their parameters, and many observed surfaces fit only one model effectively (Supplementary 2). To ensure an unambiguous choice that applies the principle of Occam's razor when fitting several models to an observed surface, we ordered them by increasing complexity (HSA, Loewe, Bliss, potentiation) and declared as best fit the first model that was consistent with the observed surface. Within Bliss boosting, we classified combinations as cancelling, suppressing, masking, multiplicative, or saturating, depending on which fell closest to the best-fit b for the observed response shape. Note that for some pairings, the boost levels cannot be separated. For example, if one single agent reaches 100% inhibition at high concentration, it will be impossible to distinguish saturating from multiplicative or masking levels. If both single agents reach 100%, only cancelling effects will be distinguishable from the others. For such boost level degeneracies, the first consistent level in the sequence above was chosen.
We note that these models were used primarily to characterize response surface shapes. Indeed, of the four shape models discussed here, only Loewe additivity has an a priori physical foundation as a predictive combination model. Although the original Bliss independence model (corresponding to the multiplicative level of Bliss boosting) is based on statistical independence, that foundation has no clear relevance to biological systems, where widely separated targets are neither statistically independent nor cleanly competing for the same end point. The Bliss boosting and potentiation models were included in our analysis to describe distinct shapes observed in our screens, and for well-sampled dose matrices neither is very good at fitting the effects that the other was designed to describe.

Pathway simulations
Simulations of metabolic pathways were performed using the XPPAUT ordinary differential equation simulation software (Ermentrout, 2002), on Boston University's Biowulf computing cluster. The network shown in Figure 3 was modelled as a system of Michaelis-Menten enzymatic reactions, each with three kinetic parameters: the limiting reaction rate V m , the rate constant K m , and an exponential degradation rate D g . The two pathways were driven by supplying substrates S1 and S11 from limitless sources with constant reaction velocities, and converting them through the enzymatic reactions before passing through the endpoint reaction R E into a limitless sink. For symmetry and simplicity, we set all V m ¼1 and K m ¼1, and the degradations D g to a lower level of 0.01 in most cases. We used D g ¼1 for substrates at bifurcations and for the final product of Pathway A to ensure illustrative differences in the single agent inhibitions (reducing response shape degeneracies) and to balance the flux through both pathways (for stable simulations). The system was initially started with all substrate levels set to 0.1 and run to steady state using an RK4 integration algorithm calculated at time intervals of 0.001 over 1000 000 iterations, and the resulting substrate levels and parameter settings are shown in Supplementary 1. The combination response for each pair of inhibitors was determined using 100 simulations over a 10 Â10 dose matrix, where for each dosing point the targeted reactions were restricted by competitive inhibitors (that affect K m ) and the simulation was run for a further 600 000 iterations using a step size of 0.005, before the inhibition I¼1-T/U was calculated on the treated end-point reaction velocity T, relative to the untreated velocity U. Optimal 9-point inhibitor concentration samples were chosen for each single agent by using a bisection search to find the single agent's EC 10 and EC 90 concentrations, defining six samples with a fixed dilution factor f¼(EC 90 /EC 10 ) 0.2 to cover that range, and extending the total range with two additional concentration steps above the EC 90 and one more below the EC 10 at the same dilution factor f.
The resulting combination response data are tabulated in Supplementary 2. Each surface was compared to the models presented in Figure 2 (assuming 1% inhibition errors), and the best fit was determined to be HSA, Loewe additivity, Bliss boosting, or potentiation, in that order, depending on which was first consistent with the best-fitting model's w 2 (where consistency meant that its w 2 was less than twice the best fit's). Bliss boosting models were further classified as cancelling, suppressing, masking, multiplicative, or saturating, where it could be uniquely specified given the single agent limiting efficacies. The shape classifications, chi-squared, and best-fit parameter values are listed in Supplementary 2, and the shapes are also summarized in Figure 3. Supplementary Figure S1 compares the shape classifications and target connectivities from this simulation to epistasis calculations (Segrè et al, 2005;Yeh et al, 2006) based on single concentration points taken from the matrices. We varied the model parameters to explore the sensitivity of our results to the kinetic parameters. First, we performed Monte Carlo tests with 100 simulations, each time allowing all of the parameters to randomly vary with uniform probability density within a multiplicative range 10 7D , and counted how often the best-fit model surface agreed with the original fit before perturbation (Supplementary Figure  S2). The Monte Carlo tests were repeated with steadily increasing D, to determine the level at which the unperturbed simulations failed to predict the outcome, both with all V m , K m , and D g varied simultaneously, and with each parameter class perturbed separately. Finally, we explored sensitivity to individual nodes by perturbing each enzyme in Pathway A separately.
We also explored some simulation assumptions to test the robustness of our results. First, we replaced the competitive inhibitor kinetics with non-competitive (affecting V m ) inhibition, performed the simulations, and repeated the model fitting to determine the best-fit response shapes. Second, we addressed the concern that a metabolic reaction velocity might provide a poor comparison to proliferation, the most accessible experimental end point. As proliferating cells undergo exponential doubling, a linear reduction of velocity by a factor T/U for a growth-limiting metabolic reaction should slow the doubling time by a proportional factor. Over a time interval corresponding to n uninhibited doublings, the resulting inhibition of cell population becomes I growth B1-2 nIv , in terms of the velocity-based inhibition I V . We applied this transformation to our simulated dose-matrix inhibitions for Pathway A, assuming nB10 doublings, and repeated the surface fitting procedure on the transformed dose matrices to determine the best-fit response shapes.

Combination screening experiments
The chemical library was archived in robotically accessible vials, to which diluent (dimethyl sulphoxide or water) was added in preparation for addition to 384-well master plates by a Tecan Freedom liquiddispensing robot. Liquid transfers to dilution and assay plates were handled using a Perkin-Elmer MiniTrak station adapted for the combination high-throughput procedure. Each 384-well assay plate contained six 6 Â 6 dose-matrix blocks, with four serial dilutions (twofold for yeast and four-fold for HCT116) of the top concentration for each agent. Additional wells were reserved for transfer and untreated control wells. The compound mixtures were then added to the biological component of the assay as appropriate.
For the yeast experiment, C. glabrata (ATC #90030) cultures were cultured overnight in flasks with RPMI media in a 321C shaker at 250 r.p.m. The next day, the cultures were diluted into media until they had an absorbance equivalent to 50% of McFarland standard, corresponding to a density of B8000 cells/ml, and dispensed with automated dispensers into 384-well assay plates (35 ml wells containing media and test compounds) to yield B160 cells/well. The plates were incubated in the presence of drugs and media at 321C for 18 h (B10 doublings). Upon completion, cell populations were measured at a single time point, with a metabolic Alamar Blue fluorescence reagent (at 590 nm emission) using a Perkin-Elmer Fusion plate reader after excitation at 535 nm.
For the tumour cell line experiment, the HCT116 (ATC #CCL-247) cells were cultured in flasks with RPMI-1640 media for 2-20 passages, drawing test populations along the way. Extracts were trypsinized, counted, and seeded into 384-well plates at a density of 1500 cells/well in 35 ml media using automated dispensers. To ensure adherence, the plated cells were incubated at 371C, with 5% CO 2 overnight, after which the compounds were added and the plates were incubated again at 371C with 5% CO 2 for 72 h (or about 2-3 doublings). After incubation, 10.5% Alamar Blue fluorescence dye (in 40 ml media) was added and plates were incubated for 6 h. Cell population measurements were made at a single time point upon completion, using a Perkin-Elmer Victor II plate reader (excitation at 535/590 nm emission).
This experimental procedure has a variety of sources for error in final proliferation measurements. Each compound is dissolved once in a single stock tube, and delivered into X-and Y-diluted master plates. Aliquots are then transferred from the master plates to each assay plate to produce the combinations. While this approach permits orthogonal serial dilutions and aids uniformity across instances of each master plate, there can be systematic errors in compound plating between the X and Y masters for less stable compounds. Cells are seeded in media at high densities, to reduce the well-to-well variability, but occasional clogs in a pipette can lead to non-uniform distribution. During incubation, slight variations in temperature or humidity can lead to varying growth rates across wells on a plate and between plates. Because the single cell population measurement is made during the exponential growth phase, variations can be amplified in the final readout, and this variation can be exaggerated when cells are sensitized by test compounds. The result of these various factors is that cell population measurements vary by B1-3% on each plate, across untreated wells, and the variations do not fit a normal distribution. Moreover, some 1-2% of wells on each plate show occasional spikes that are very different from their neighbours.
The contents of each plate were tracked in an automated laboratory information management system, using integrated barcode scanners in the liquid handling equipment, and stored in an Oracle database. Plates lacking compound transfers or with insufficient dynamic range (signal-to-noise ratio o5 between untreated controls and a cell-free background) were rejected and repeated. The remaining plates were inspected using custom quality control software. Individual wells that fell outside the expected range for normal assay readouts, or which were discontinuous with their neighbours, were marked for exclusion. Finally, the single agent wells in each combination block were visually inspected for consistency across the experiment, and combination blocks containing the most discrepant single agents were marked for exclusion.

Data analysis
Dose matrices were assembled from replicate combination blocks on experimental plates. Some of the data showed systematic variations across the plate, most likely owing to temperature or humidity gradients during incubation. In order to correct these effects, our assay plates had 40 untreated wells arranged around and between the blocks. The systematic variations were removed by fitting a smooth function of plate location to the untreated well data values, and dividing out the modelled spatial variation. After plate effect correction, fluorescence counts T from each treated well were converted to inhibitions I data ¼(UÀT)/U relative to the median U of 20 untreated wells. Replicate blocks for each dose matrix were merged, using the median inhibition at each dosing point. Standard error estimates s data for each median inhibition were also calculated, based on the quadrature sum of a minimum acceptable 3% error, the median absolute deviation (MAD) of the corrected untreated data on each plate normalized to their median, and the MAD between replicate inhibition data between plates, using an empirical conversion from MAD to standard deviations (Filliben, 2005).
Each dose matrix was scored for synergy, and combination effects were compared to our response morphology models (Figure 2). Sigmoidal dose responses I¼EC a /(S a þ C a ) as a function of concentration C were fitted to the single agent data, where E is the limiting efficacy at high concentration, S is the effective concentration, and a is the sigmoidicity or steepness of transition. These fitted curves were used for combination effect models. The HSA, Bliss boosting, and potentiation models were calculated, and the Loewe additivity model was solved at each dosing point by iterative root finding (Berenbaum, 1985). As a quantitative estimate of the combination activity, we computed HSA volume excess scores V HSA ¼ P X,Y ln f X ln f Y (I data -I HSA ), summed over all C X ,C Y , along with an associated error estimate s HSA ¼sqrt[ P X,Y ln 2 f X ln 2 f Y s 2 data ], both expressions using adjustments to account for differing dilution factors f X , f Y . This volume measures the total excess inhibition, summed over dimensionless log-concentration space. Data surfaces were compared to the shapes presented in Figure 2 by varying the model parameters, using a Nelder-Mead simplex algorithm (Press et al, 1997, y10.4) to minimize w 2 ¼S data [(I data -I model )/s data ] 2 /N data , the reduced chi-squared. All four models were explored, and the overall best fit with the least chi-square w 2 min was chosen to represent the shape of the combination. A best fit was chosen from HSA, additivity, boosting, or potentiation, in that order, depending on which was first consistent with the w 2 min (Dw 2 ow 2 min ). Finally, systematic errors that were not captured by our error estimates from replicates can be estimated by examining the drugwith-self combinations as a negative control for combination effect. For both screens, the median Loewe additivity w 2 was B9, suggesting that the true reproducibility of drug response was roughly three times the estimated levels from the scatter between replicates. This suggests systematic errors of B10% at each dose-matrix point, leading to integrated volume error estimates of s HSA B0.4. Thus, only combination effects with V HSA 41 are strong enough to be reliably distinguishable from drug-with-itself combinations.
For the HCT116 combination screen, synergy profiles were constructed for each drug by collecting into a vector all the V HSA scores involving that chemical. These profiles were then compared for all pairs of drugs in the screen by means of the Pearson correlation coefficient. The synergy profile correlations and synergy scores were grouped by the target similarity (Supplementary 5) of the two drugs in each combination, and classified as follows: 'Identical' for drug-withitself; 'Same' for distinct drugs inhibiting the same target; 'Related' for differing targets sharing a common functional group; 'Kinases' for combinations of kinase inhibitors; 'Different' for unrelated pairs of drugs; and 'Unknown' for pairs involving one or more unknown mechanism. The target similarity group distributions were compared for significant differences using a chi-square test that assumes only Poisson counting statistics within each bin in the histograms (Press et al, 1997, y14.3.3).

Supplementary information
Supplementary information is available at the Molecular Systems Biology website (www.nature.com/msb).