Predicting and characterizing selective multiple drug treatments for metabolicdiseases and cancer

Abstract

Background

In the field of drug discovery, assessing the potential of multidrug therapies isa difficult task because of the combinatorial complexity (both theoretical andexperimental) and because of the requirements on the selectivity of the therapy.To cope with this problem, we have developed a novel method for the systematic insilico investigation of synergistic effects of currently available drugs ongenome-scale metabolic networks.

Results

The algorithm finds the optimal combination of drugs which guarantees theinhibition of an objective function, while minimizing the side effect on the othercellular processes. Two different applications are considered: finding drugsynergisms for human metabolic diseases (like diabetes, obesity and hypertension)and finding antitumoral drug combinations with minimal side effect on the normalhuman cell. The results we obtain are consistent with some of the availabletherapeutic indications and predict new multiple drug treatments. A clusteranalysis on all possible interactions among the currently available drugsindicates a limited variety on the metabolic targets for the approved drugs.

Conclusion

The in silico prediction of drug synergisms can represent an important tool forthe repurposing of drugs in a realistic perspective which considers also theselectivity of the therapy. Moreover, for a more profitable exploitation ofdrug-drug interactions, we have shown that also experimental drugs which have adifferent mechanism of action can be reconsider as potential ingredients of newmulticompound therapeutic indications. Needless to say the clues provided by acomputational study like ours need in any case to be thoroughly evaluatedexperimentally.

Background

In spite of the advances in molecular and computational biology, the discovery of newdrugs still remains a very challenging task which requires a very long period ofresearch and development before any new compound can be commercialized. A possiblealternative to the search of new active compounds is to make use of the unexploitedproperties of already available drugs, since a wide knowledge about both theirtherapeutic and toxicity effects has already been gathered during the study for theirapproval. In this perspective, a natural approach to broaden the range of applicationsof the existing drugs is to try to combine them in multiple drug therapies [1–3]. However, even though both the financial burden of conducting trials as wellas the risk of adverse events in trial populations is expected to be sensibly lower foralready approved drugs, so far the experimental investigations of multicomponenttherapies have been quite limited [4, 5]. Major obstacles to this approach are the high number of possiblecombinations but also our limited understanding of the complex mode of action of amultidrug treatment. Indeed, multiple perturbations can show three types of interaction,which have been classified as synergistic, antagonistic and additive [6] (alternatively called, aggravating, buffering and non-epistatic [7]). We have focused our attention on the first type, where the use of drugcombinations represents an enhancement with respect to the superposition of the singleperturbations.

In order to identify synergistic effects, Ref. [4] investigated all pairs of a set of known drugs at different doses, obtaininga functional classification of the interactions by looking at their inhibitory effect.Various computational approaches, based on reconstructed genome-scale metabolicnetworks, have been also developed in Refs. [7–10], in order to identify the synergistic effects triggered by multiple drugs ormultiple genetic perturbations (for example the so-called “syntheticlethality”). Unfortunately, a systematic evaluation of the effects of all possiblecombinations of drugs is unfeasible, because their number scales exponentially with thenumber of chemicals taken into account (for an exhaustive search over 40 drugs, morethan a trillion possible combinations should be tested). Moreover, a drug profile isgiven by both its therapeutic effect and by its side effect, the latter being related toits selectivity. For example, drugs such as anticancer agents have to selectively actonly on tumor cells [11]. Similarly, metabolic diseases are induced by the imbalance of key metabolicpathways which, if modulated without affecting other vital functions, can rescue fromthe pathology [12, 13]. Referring to these specific needs, we have developed an algorithm based onthe metabolic network of humans and on the comparison between the metabolic networks ofhuman and cancer cells (as reconstructed recently in [14]) aimed at expanding the spectrum of applications of the existing drugs to newselective treatments against metabolic diseases and tumors.

The algorithm here proposed is based on Flux Balance Analysis (FBA) [15, 16]: in spite of its simple formalism, FBA has already proven to be reliable inproviding quantitative understanding of cell metabolism. The computational methodpresented in this paper is based on a bilevel optimization which, after reformulationthrough duality theory, allows the algorithm to efficiently search the interactionsbetween drugs. With respect to the available literature [7–9, 17–19], the procedure we are proposing presents at least three importantdifferences: (i) the synergisms are efficiently explored over all drug combinationswithout limiting only to pairwise combinations but without doing an exhaustive search,thanks to the application of duality theory; (ii) the multiple drug treatments suggestedby the method guarantee both the inhibition of the chosen target (efficacy) and aminimal side effect on the other cellular functions (selectivity); (iii) in ourprocedure, any metabolic process of the network can be chosen as possible disease andphenotype readout, not only cell growth as common in the FBA literature (a more detailedcomparison with the current literature is reported in the Additional file 1). Inspired by works such as Refs. [20, 21], we treat the inhibition of a metabolic reaction by a drug as the silencingof the gene which codes for the catalyzing enzyme. In this simplified framework, thesynergistic effect resulting from multiple perturbations of the metabolic network isstill well captured [7, 22] (see Figure 1 for a simple example). The selectivity ofany drug treatment is correlated to its side effect, which is estimated as number ofstopped reactions (see Figure 1 (C) _obj vs (D)) plus a correction term for the known non-metabolic targets. Furtherdetails are described in the Methods Section. It is worth noting that, because of thesteady state assumption of FBA, this formalism does not identify synergisms betweendrugs that manifest themselves as alterations of kinetic parameters and consequently ofconcentration of metabolites, like the case of drugs which act on the same linear path(as e.g. Trimethoprim and Sulfamethoxazole do on folate synthesis [23]). Similarly, interactions where for instance one drug inhibits thebiodegradation of the others cannot be found by FBA-based methods [24].

Figure 1

Example of drug synergism in FBA. For the toy network depicted in(A) the aim is to stop the objective reaction v₁₀ (inred) by choosing a combination of drugs (the three valves “⋈”)while blocking the minimum number of reactions other than v₁₀.In the drawing, blue arrows indicate active fluxes while gray arrows refer tostopped reactions; the valve is red if the drug is used, gray otherwise. Panel(B) shows how the use of a single drug (v₄ orv₇) does not stop the objective reaction, while the drug atv₂ blocks the objective reaction v₁₀but it also blocks all fluxes of the network (panel (C)). Therefore, theoptimal drug combination blocking the objective function v₁₀with minimal side effect is given by the synergism of the two drugs acting atv₄ and v₇ (panel (D)). Thecomparison of panels (B) and (D) shows how a synergism is a behaviorwhich cannot be simply inferred by the superposition of the effects of the singledrugs, but that structurally depends from the topology of the network.

In spite of the limitations of restricting our analysis only to metabolic networks, theresults presented in the following provide a promising example of how in silicomodels can be used as practical tools for exploring the functional interactions betweendrugs and of the (little explored) potential offered by synergistic drugcombinations.

Results and discussion

In order to explore its potentialities and limitations, the algorithm (described in theMethods Section) has been applied to two different case studies:

• 1. finding drug synergisms for metabolic diseases (like diabetes, obesity and hypertension) on the human network [25];

• 2. finding antitumoral drug combinations with minimal side effect on the normal human cell (using the cancer network of [14] to model the metabolism of a human tumor).

Some features of these two metabolic networks are listed in Table 1, together with the number of drugs currently approved (from [26]). In particular, in order to generate realistic solutions, the availableinformation about these existing drugs has been carefully filtered, selecting onlyinhibitions of metabolic human targets that have been experimentally proven (moredetails in Methods, Additional file 1: Table S1 and Additionalfile 2).

Table 1 Features of the metabolic networks considered in the paper

Human metabolism

In this Section, we consider the inhibition of specific functions of the humanmetabolism obtained without impairing other vital processes. First, the inhibitoryeffect of each single drug on the whole network has been calculated. Then, thefollowing screening is performed: we systematically consider each reaction of thenetwork as a potential objective function and we apply the algorithm, searching forthe most selective synergism capable of blocking this objective reaction. Comparingthis solution with the single-drug effects evaluated in advance, we distinguish threecases:

• a) the drug combination leads to a new inhibition since no singledrug can stop the objective reaction;

• b) the objective reaction can be stopped also by a single drug but the drugcombination is more selective (has a minor side effect);

• c) the objective reaction can be stopped also by a single drug and themultiple drug solution is less selective (this solution is not interestingbecause it triggers a larger side effect).

In all cases where a single or a multidrug solution is found, also all suboptimalsolutions are hierarchically identified, iterating the procedure while excluding thecurrent optimum, until the problem becomes unfeasible (i.e. no more solutions exist,capable of blocking that objective reaction). At the end of the screening, weobtained a set of 32 multicomponent solutions, ranging from combinations of two up tofour compounds (see Table 2). The following characterization ofthe synergistic effects is performed. For each combination we identify the set Y of metabolic reactions which cannot be stopped by any single drug of thecombination, but which are stopped when all these drugs are used together. Then, thesynergism is described by the vector $\mathbf{\text{s}}\in {\{0,1\}}^{N_r}$, where s _j  = 1 if reaction j belongs to Y (N _r is the number of reactions in the metabolic network). From the vectorss of the 32 multiple drug solutions a matrix of distances can beconstructed and a cluster analysis performed on these distances; the resultingdistance-based tree (similar to a phylogenetic tree) is drawn in Figure 2 (upper panel). The synergisms are clustered in six classes (withclearly identifiable subclasses in some of them) labeled from “A” to“F”. This classification can be used to build also a proximity networkfor the drugs, linking those that belong to the same synergistic interaction. Theoutcome is drawn in Figure 2 (bottom panel) and shows that thesame clustering applies to the drugs involved in the synergisms. The resulthighlights how drugs can often be used in alternative one to the others: for example,the synergistic pairs of class C contain one drug among those labeled with the number7 or 55 (Rosiglitazone or Acetylsalicylic acid) in combination with one drug amongnumber 4 or 43 or 60 (Pravastatin or Naftifine or Tioconazole; see Additional file1: Table S1 for all correspondences between names andnumbers). Being the cardinality of class C equal to 6, we can deduce that thesesolutions are generated only by the combination of the pair and the triplet justmentioned. Three exceptions to the sharp clusterization of Figure 2 are represented by drugs labeled with the numbers 7, 22 and 50,respectively Rosiglitazone, Indomethacin and Pentoxifylline. Indeed Rosiglitazonetargets many metabolic reactions (60, all in the fatty acid metabolism) which allowtwo types of interaction: class A for fatty acid activation and class C forcholesterol metabolism. On the other hand, Indomethacin causes only 7 inhibitions,some belonging to glycerolphospholipids metabolism and some others to pyruvatepathways: the first interact synergistically with drugs which target fatty acidreactions (class A), whereas the second can be combined with drugs acting on pyruvatemetabolism (like Fomepizole, drug number 75, in class B). Finally, Pentoxifyllineinhibits reactions both in the salvage pathway for nucleotides (which give synergismsin class E) and in pyrimidine catabolism (class F).

Figure 2

Classification of the synergisms for the human metabolic network. Top panel: Each leaf of the tree represents a multidrugsolution that we have found. The layout of the graph is obtained through thesame method used for phylogenetic trees (a distance tree, see text) andmanifestly shows the clustering of these synergisms; the six clearly visibleclasses have been labeled with letters (from “A” to“F”). Names of the pathways mainly affected by each classare reported near the clusters. Bottom panel: This network of drugsrepresents a detailed characterization of the classes of synergisms. Each drugis indicated by a circle (whose radius is proportional to the number of itsdirect targets; drugs are labeled with numbers according to Additional file1: Table S1). Each synergism is drawn as a coloredline which connects the drugs involved (each synergism has its own color andthe line thickness is proportional to the number of stopped reactions). Even inthis more detailed representation, the six classes are still visible. Somesubclasses can also be identified: drug pairs (7, 62) and (7, 65) in class Aand drug pairs (42, 58) and (51, 58) in class E (indicated with broken lines)exploit part of the synergism of the entire class; indeed these 4 pairs are theisolated leaves in the corresponding clusters in the top panel. Note the roleof drugs 7, 22 and 50 in bridging classes A-C, A-Band E-F. Details of the metabolic functions to which theseclasses of synergisms correspond are given in Figure 3.

Table 2 List of all drug synergisms

The analysis of the complete results obtained from the screening over all metabolicreactions is shown in Figure 3, where the stoppable reactionsare grouped on the basis of the metabolic pathway to which they belong. Figure 3 reports also the class and the targets of the drug combinationswhich induces the inhibition. As one can see, some synergisms occur between reactionswhich belong to different pathways: in particular, sphingolipids subsystems, CoA andpyrimidine biosynthesis contain reactions whose inhibitions are caused byinteractions situated in other pathways, since none of the combined drugs havetargets on them. Moreover, among the multiple drug solutions of Figure 3, there are several new inhibitions and a few moreselective cases (a comparison of the side effects induced by single andmultiple drug solutions is reported in Additional file 1:Figure S1).

Figure 3

Drug synergisms for the human metabolic network. Left panel:For each affected pathway, the histogram reports the number of objectivereactions which can be stopped; gray-scale bars represent reactions stoppedonly by a Single drug or multidrug solution, classified as Newinhibition (meaning that no single drug is capable of triggering theinhibition), More selective and Less selective inhibitions(referring to the case where both single and multiple drug treatments arepossible and the multiple one has respectively a lower and a higher sideeffect). Right panels: The two plots refer to multiple drug solutionsonly. For the same pathways as in the left panel, we report here the fractionof the direct drug targets and the fraction of the synergistic inhibitionswhich are induced by the six classes of synergisms (shown in Figure 2): the comparison between the two stacks shows thatsynergistic interactions can occur on pathways that are not direct targets ofthe drugs.

Among the cases of new inhibitions, the case of Guanylate kinase, although of notherapeutic interest, represents an easily visualizable example of the nonlinearityin the superposition of the effects as anticipated in the hypothetical situationpresented in Figure 1. We consider the phosphorylation of GMPinto GDP catalyzed by guanylate kinase as objective reaction. Since the blockage ofGMP production will cause also the arrest of any transcription process, thisinhibition constitutes only a toy example of synergism devoid of any practical value.For this problem, the algorithm proposes the combination of Mercaptopurine,Dipyridamole and Mycophenolic acid: the synergism takes place through thesimultaneous inhibition of guanine phosphoribosyltransferase, 3^′,5^′-cyclic-nucleotide phosphodiesterase and IMP dehydrogenase (see Figure 4 for a representation of the corresponding subnetwork). Indeed,these reactions are alternative ways of GMP biosynthesis. When and only when they areall blocked, Guanylate kinase lacks its substrate and stops as well.

Figure 4

Nonlinearity in the synergism: the example of Guanylate kinase. The partof the human network here represented shows the nonlinear interaction when thethree drug targets (the three valves, with the name of the drugs) aresimultaneously inhibited: when this is the case, the objective reaction ofGuanylate kinase (in red) is stopped. Gray arrows and gray circles indicaterespectively stopped reactions and metabolites which become unavailable.

The complete list of the objective reactions, with relative pathways and synergisticinhibitions, is reported in Additional file 1: Table S2:this list contains many inhibitions in the fatty acid, cholesterol and carnitinetransport pathways, which may represent solutions for obesity. In particular,concerning the case of hyperlipidemia diseases, the algorithm finds the combinationof Rosiglitazone and Cerulenin (class A) as inhibitor of many reactions in thecarnitine transferase and fatty acid desaturase pathways. This synergism has beenreported in the literature for being active versus the biosynthesis of fatty acids inprostate tumors [27]: the mean IC₅₀ are 45μ M and 32μ Mfor Rosiglitazone and Cerulenin alone, whereas it reduces to 5μ M whenthey are combined: the authors claimed that this effect comes from the reducedproduction of fatty acids preventing the growth and differentiation of prostatecells. It is worth noting that we predict this combination also in three anticancersolutions (together with an additional target on palmitate conversion, see nextsubsection). Moreover, this pair is part of other six synergisms on the humanmetabolism, all concerning the same pathways (i.e. a very similar mechanism ofinteraction).

Another significant example is represented by the inhibition of Dihydroceramidedesaturase. Indeed, ceramide is the hydrophobic membrane anchor of sphingolipids andis involved as a bioactive molecule in cell growth regulation, apoptosis, senescence,and diverse cell responses, particularly those linked to stress situations [28, 29]; moreover, recent studies have shown the role for ceramide biosynthesis inbody weight regulation, energy expenditure, hence in the metabolic obesity syndrome [30]. For these reasons Dihydroceramide desaturase has been proposed as apromising potential target for metabolic diseases. Currently some specific inhibitorsof Dihydroceramide desaturase are under investigation and development (for exampleGT11, XM462 and analogous [31, 32]) although there is no approved drug yet. Indeed in our model no singledrug can stop this reaction. Our algorithm finds some possible multidrug treatmentsthat block this reaction: among them, there are the synergistic pairs ofRosiglitazone plus Simvastatin (Pravastatin), and Acetylsalicyclic acid plusAtorvastatin (Pravastatin). Concerning the first synergism, clinical experiments haveshown that combining these two drugs a significant reduction (about 30% less) of theintracellular accumulation of lipid is achieved [33]. Moreover, the authors of [34] investigate the adverse effect of single and combined therapies(hypoglycemia, body weight increase) and claim that adverse events are generallysimilar (the safety profile of Rosiglitazone was not adversely affected by theaddition of Atorvastatin). Also our results predict a limited worsening of theadverse effect: indeed, after the combination with Atorvastatin the side effect ofRosiglitazone passes from 252.7 (Additional file 1: TableS1) to 288.8 (Table 2), i.e. it increases of about 14% only.For the same therapeutic purpose, the pair of Acetylsalicyclic acid and Atorvastatinhas been also studied. Clinical trials are currently ongoing [35] and some of them have already shown promising results [36]. The rationale for this approach is based on the restoration of plateletsensitivity by reduction of the cholesterol levels. Moreover, as mentioned above, itis known that ceramide is involved in apoptosis. Indeed, this combination has beentested for the treatment of prostate cancer: the results have shows a linearsynergism between these two drugs [37].

Human vs Cancer

A similar approach can be used to improve the selectivity andspecificity of the treatment when dealing simultaneously with more thanone type of cells. With a small adjustment, our procedure can force the solution topreserve the metabolism of one cell while inhibiting an objective reaction of another(see Additional file 1 for a detailed formulation). Indeed,a drug interaction can explore the differences in the topologies of the metabolicnetworks and, in this way, bypass the restrictions caused for instance by targetshomology. This is crucial in case of anticancer therapy since tumoral and normalcells share the same genes.

For this purpose, we use the metabolic network of a generic human cancer assembled in [14] and the human metabolic network. We apply the modified version of thealgorithm to the biomass reaction of the cancer network (which must be stopped),while minimizing the side effect on the regular human metabolism. As in the previoussection, the procedure is iterated until the problem becomes unfeasible. The resultsare shown in Figure 5 and Table 3. Thesolutions are mainly single drugs which differ one from the other in terms of sideeffect on the human network. Many are known chemotherapeutic agents such asFloxuridine, Mycophenolic acid, Methotrexate, Pemetrexed, Ribavirin, Myo-Inositol,Simvastatin, Leflunomide, Indomethacin, Hydroxyurea, Arsenic trioxide, Gemcitabine [38–49]. Since only one synergism is present, between Fomepizole and Auranofin,these results suggest that approved drugs do not seem to induce significantinteractions at the level of the metabolism.

Figure 5

Results on cancer vs human selectivity problem. Left panel: The iterative application of the algorithm to thecancer vs human networks finds 21 solutions before becoming unfeasible(magenta line). Including the possibility of inhibiting an additional target,other 31 solutions are found (blue line). Right panel: The bars countthe number of solutions which stop the biomass metabolite in the cancermetabolism (same color code). Solutions have been classified according to thenecessity or less of the inhibition of an additional target (see Tables 3, Additional file 3: Table S3 formore details).

Table 3 Solutions for the cancer vs human problem

In order to increase the range of putative synergisms, we have tried to use thealgorithm to set up a search for targets potentially interacting with the currentlyavailable drugs. Instead of searching for the reactions which are syntheticallylethal (as partially done, for instance, in [14]) and proposing them as potential new targets, we look for the reactionswhose single inhibition may give a lethal synergism with any combination of theapproved drugs. The problem is formally equivalent to the one we have alreadydescribed: we search again for the optimal drug combination after having deleted areaction in the cancer network (the same reaction is removed also from the humannetwork; this reaction is called “additional target” since it will be thetarget of an additional new drug). We systematically consider each reaction of thecancer network as an additional target; the results of this screening are reported inAdditional file 3: Table S3. After a search in literatureof possible inhibitors of the additional targets of the results, we identify someinteresting solutions. For instance, the inhibition of methenyltetrahydrofolatecyclohydrolase (combined with the use of Mimosine) can be induced by the experimentaldrug 5,6,7,8- tetrahydro-N⁵,N¹⁰- carbonylfolic acid [50]; therefore, its combination with Mimosine could represent a potentialantitumor therapy. Also for other four additional targets there exist experimentalinhibitors whose activity is reported in the literature [51–53] (details are shown in Additional file 3: TableS3). Concerning the drug solutions identified by the algorithm, the pairRosiglitazone plus Cerulenin is proposed in combination with three possibleadditional targets. One is the palmitate fatty acid conversion: it is worth notingthat the validated antitumoral activity of this pairs versus prostate cancer cells(as mentioned above) is due to the reduction of the synthesis of fatty acids [27]. In our prediction, indeed, among the metabolites which are no longeravailable because of this inhibition, there are cholesterol ester, mono-, di- andtri-acylglycerol. The other two additional targets are related tophosphatidylcholine. Part of the phosphatidyilcholine pathway has been alreadyidentified as synthetic lethal [14], but without mentioning any possible exploitation. Our results suggest theuse of a combined drug therapy as possible way to take advantage of this syntheticlethality.

Conclusions

The field of drug combinatorics is largely unexplored experimentally and the potentialof combined drug therapies is difficult to assess, mostly for lack of suitablesystematic methodologies. To try to fill this gap, we have developed an algorithm whichis capable of exploring efficiently the optimal synergisms among all possible drugcombinations and of characterizing them in terms of side effect and selectivity. Indeed,the success of a drug discovery process depends on multiple aspects, not least on thefulfillment of requirements regarding selectivity and toxicity: for metabolic diseasesthe modulation of the key pathways without affecting the other vital functions can beinstrumental for rescuing from the pathology. Similarly, anticancer compounds shouldonly kill cancer cells without affecting normal ones. These requirements are rarelytaken into account in standard computational approaches. The results we obtained byapplying our algorithm to the human and tumor vs human metabolic networks showthe possibility to take advantage of drug synergisms in proposing new therapies: thepotentialities lay in the possibility to intervene with a different mechanism of actionwith respect to those that are currently available. In this enlarged repertoire ofpossibilities, we have identified examples of drug repurposing (some of them werepreviously demonstrated experimentally), a procedure which is becoming more and moreattractive thanks to the reduced costs on the preclinical and clinical steps.

One of the main features of FBA-based knockout studies is that metabolic networks appearto be robust [1, 17], meaning that there seem to be an high degree of redundancy of the pathwaysinside a network (property alternatively reported as “nonessentiality” ofthe gene in Refs. [8, 19]). In the context of drug synergism, this property reflects into the presenceof optimal solutions consisting of many drugs (in our case up to four, or even more ifwe consider the number of inhibition targets of each drug). For the same reason, theresults show the necessity to extend the search to all possible drug combinationswithout limiting to those of low cardinality. This fact becomes significant especiallywhen the drugs to combine present a high similarity in terms of inhibited targets;indeed, the characterization of the synergisms we have found shows a limited variety ofpossible interactions between the available drugs (only six classes were identified, seeFigure 2). However, beside increasing the cardinality of thesolutions, a very strong robustness may also reduce the total number of solutionsbecause it makes more difficult to induce the simultaneous inhibition of all theredundant pathways; indeed, if the screening we have run on the human metabolism isapplied also to the less robust cancer network (seen as stand-alone network), the numberof possible inhibitions is much higher (see Additional file 1:Figure S2 in comparison with Figure 3). Moreover, when we want tostop the biomass reaction of the cancer, mainly single drug solutions are found: this isagain an index of the low redundancy of the cancer network (see Table 3) and of the limited variety in the metabolic targets for the availabledrugs.

Nevertheless, adding the possibility to inhibit an extra target, we could identifiedsome experimental compounds (other than the drugs from DrugBank) which may be used asanticancer in a combination with the approved drugs (see Additional file 3: Table S3). These examples show that predictive tools like themethod we are proposing become more important if one considers also the possibility ofcombining active compounds which are not yet approved but for which a minimalcharacterization of the mechanism of action is available. In this perspective,experimental compounds which inhibit additional targets that are different from thoseaffected by the approved drugs may represent a good chance for improving, throughsynergism, the spectrum of the whole set of the currently available drugs. Moreover, theapplication of our method can be extended to situations where multiple networks arecompared and contrasted. It is expected that problems like this will become important assoon as tissue-specific networks of human metabolism and cancer-specific networks willbecome available in the near future.

In a broader perspective, if instead of confining our study only to the human network weconsider also the metabolism of microorganisms, the exploitation of drug synergismsobtained with our algorithm can be useful in investigating a wide range of situations:(i) when specific enzyme inhibitors are not currently available, multiple drug solutionscould represent an example of the reprofiling of existing drugs for newtherapeutic indications; (ii) when the target enzyme has undergone a mutation renderingineffective the original therapy, a synergistic solution may bypass the resistanceacting on other enzymes and therefore help in fighting resistance; (iii) whenthe optimal synergism has no lethal impact while the single drug solution has. Thischange in lethality can be important for instance in cases of human-hostedbacteria producing toxic by-products: in order to save the useful symbiosis with thesecommensal bacteria, a selective (but not lethal) inhibition of the toxic processes mustbe pursued [54, 55].

Apart from the specific cases we have studied, the main objective of this paper is topropose an efficient method to single out drug combinations with potentially interestingtherapeutic effect. Given the exponential character of the combinatorics involved,unguided “fishing expeditions” are intrinsically ill posed. Hence we expectthat methods like the one proposed in this paper should be of help in dealing with sucha complex problem. Needless to say, as the predictions are based on in silicomodels of metabolic pathways, the validity of the results must be assessedexperimentally, but this is beyond the scope of the work.

Methods

The drug-synergism algorithm is developed in the framework of FBA (see SupportingInformation for a brief comparison with the available literature on similar methods).The problem deals with the following sets (and numbers):

$$\begin{array}{ll}\mathcal{R}& =\{1,\dots ,N_r\}=\text{set of reactions;}\\ \mathcal{M}& =\{1,\dots ,N_m\}=\text{set of metabolites;}\\ \mathcal{D}& =\{1,\dots ,N_d\}=\text{set of drugs;}\\ {\mathcal{T}}_j& =\{1,\dots ,N_{t,j}\}=\text{set of drugs having the reaction}\\ j\text{as a target.}\end{array}$$

Then, the optimal synergism problem can be stated as follows:

Problem: Given:

• a metabolic network, which means a stoichiometric matrix $\mathit{S}\in {\mathbb{R}}^{N_m\times N_r}$ and a vector of fluxes v with upper-bounds U , both laying in ${\mathbb{R}}^{N_r}$ ;

• an objective reaction flux v _obj which has to bestopped;

• the set $\mathcal{D}$ of drugs together with their inhibition targets;

we want to find the subset of drugs $D\subseteq \mathcal{D}$ such that D blocks v _obj causing the minimal side effect, i.e. a minimum perturbation on the overallreaction fluxes.

Of course, we are not interested in procedures which perform an exhaustive search in thespace of all drugs combinations.

By decomposing any reversible reaction in a couple of irreversible reactions, we canalways assume that fluxes have non-negative values. Then, the space H of allpossible steady state fluxes v as defined by FBA is

$$H:=\{\mathbf{v}:\mathbf{Sv}=0,0\le v_j\le U_j\forall j\in \mathcal{R}\},$$

where vj and U _j are the j-th components of the vectors v and Urespectively. Given the inhibition targets of each drug, we assume that a drug inhibitscompletely the enzymes responsible for the targeted reactions, hence stopping therelative fluxes. Therefore, adding a drug combination D to the problem meansforcing (directly or indirectly) to zero some of the fluxes. We introduce the binaryvariables ${\left\{d_k\right\}}_{k\in \mathcal{D}}$ such that

$$d_k=\left\{\begin{array}{c}0\text{if drug}k\text{is used (i.e.}k\in D);\hfill \\ 1\text{if drug}k\text{is not used (i.e.}k\notin D).\hfill \end{array}\right.$$

Then, if the k-th drug targets the j-th reaction we write:

$$v_j\le U_j{d}_k.$$

This reduces the space of feasible steady state fluxes to a subsetH(D) ⊂ H (notice thatH(∅) ≡ H). We can now introduce adefinition which quantifies the side effect, σ(D), of adrug combination D on the reaction fluxes.

Clearly, there are many possible definitions. If we were dealing with microorganisms wecould adapt the analogous MOMA (Minimization of Metabolic Adjustment) or ROOM(Regulatory On/Off Minimization) [20, 21] criteria modeling the perturbation effect (induced by knockout effect, in theliterature) on the fluxes with respect to the “wild type” fluxes. Thisrequires however to know the unperturbed fluxes (the wild type reference) which for amicroorganism corresponds to using the biomass production as cost function of the FBAproblem [56]. Unfortunately, for human metabolic network such a commonly accepted FBAcriterion is unavailable, hence it is not possible to determine the metabolic fluxes forthe unperturbed network (“wild type”). Consequently, also the calculation ofthe perturbed fluxes and the quantification of the side effect are more ambiguous. Ourchoice in this paper is to quantify the side effect as number of stopped reaction (moreprecisely the number of reactions that cannot take place because of the perturbation).This bypasses the lack of the reference fluxes and allows to quantify the number ofcellular functions which are no longer available, regardless to the type oftissue to which the human cell belongs.

Since the available drugs we selected do not have only metabolic targets, a measure ofthe side effect based exclusively on metabolic reactions disregards the perturbationinduced by the drugs on other cellular functions (for example signaling cascades,protein synthesis, etc). This additional information can be incorporated in the modelweighting each drug variable d _k in the objective function according to its non-metabolic effects (for which onlya pure superposition is considered because of the lack of more quantitative models).Following this approach, we define:

$$\sigma \left(D\right):=\sum _{j=1}^{N_r}(1-y_j)+\sum _{k=1}^{N_d}{\beta }_k(1-d_k);$$

(1)

where the parameter β _k is an estimation of the non-metabolic perturbation induced by drug k andy _j is a binary variable such that y _j  = 0 when the flux is lower than a threshold ϵ(ϵ = 0. 1 in this paper); this can be expressed bymeans of the following linear constraints:

$$\begin{array}{ll}ϵy_j\le v_j& \forall j\in \mathcal{R};\\ U_j{y}_j\ge v_j& \forall j\in \mathcal{R}.\end{array}$$

In order to avoid a double count, for reversible reactions we impose the additionalconstraint

$$\begin{array}{ll}{y}_j+y_l\le 1& \forall j\text{and}l\in \mathcal{R}\text{s.t.}{v}_j\\ \text{and}{v}_l\text{are opposite fluxes}.\end{array}$$

Although the choice of the weights β _k is quite arbitrary, we have tried to evaluate both terms of (1) using ahomogeneous criteria, relating the values of β _k to the number of non-metabolic targets of each drug. In particular we set:

$${\beta }_k:=\bar{\beta }·\text{[\# non-metabolic targets of drug}k],$$

where the numerical coefficient $\bar{\beta }$ captures the spread of the perturbation across the non-metabolicsystems. In analogy with the metabolic network, $\bar{\beta }$ is equal to the mean number of metabolic reactions that are stoppedwhen a single metabolic target is inhibited. Referring to the human metabolic networkand averaging over all drugs we selected, we obtain $\bar{\beta }=7.7$.

The optimal solution can be described as follows:

Solution: For any subset of drugs $D\subseteq \mathcal{D}$ , we can find the set H(D) and the minimal perturbation σ(D). By restricting to drug combinations which inhibit the objective reaction(i.e. such that ${\text{max}}_{\mathbf{v}\in H\left(D\right)}\left(v_{\mathit{\text{obj}}}\right)=0$ ) the optimal solution Y is

$$Y=arg\underset{\left[\begin{array}{c}\hfill D\subseteq \mathcal{D}\hfill \\ {\text{max}}_{\mathit{v}\in H\left(D\right)}\left(v_{\mathit{\text{obj}}}\right)=0\hfill \end{array}\right]}{\text{min}}\sigma \left(D\right).$$

(2)

The resulting bilevel optimization is a min-max integer linear programming problem [57]. The inner problem adjusts the fluxes in order to achieve themaximum flow for the objective reaction when all fluxes are subjected to theinhibitions (drugs) imposed by the outer problem and to the stoichiometric constraints.The outer problem selects the combination of drugs which minimizes the sideeffect, restricting to those solutions of the inner problem which guarantee no flow forthe objective reaction.

The bilevel optimization problem is the following:

$$\begin{array}{cc}\text{Minimize}\hfill & {\sum }_{j=1}^{N_r}{\alpha }_j(1-y_j)+{\sum }_{k=1}^{N_d}{\beta }_k(1-d_k)-b{\sum }_{k=1}^{N_d}{d}_k\mathit{\text{“outer problem”}}\hfill \\ \text{such that}\hfill \\ \left[\begin{array}{ccc}\text{Maximize}\hfill & v_{\mathit{\text{obj}}}\hfill & \mathit{\text{“innerproblem”}}\hfill \\ \text{such that}\hfill \\ {\sum }_{j=1}^{N_r}{S}_{i,j}{v}_j=0\hfill & \forall i\in \mathcal{M}\hfill \\ v_j\le U_j\hfill & \forall j\in \mathcal{R}\hfill \\ v_j\le U_j{d}_k\hfill & \forall j\in \mathcal{R},k\in {\mathcal{T}}_j;\hfill \end{array}\right]\hfill \\ v_{\mathit{\text{obj}}}=0;\hfill \\ ϵy_j\le v_j\forall j\in \mathcal{R};\hfill \\ U_j{y}_j\ge v_j\forall j\in \mathcal{R};\hfill \\ y_j+y_l\le 1\forall j\in \mathcal{R}\text{such that}{v}_l\text{is the opposite flux of}{v}_j,\hfill \end{array}$$

where the parameter b ≪ 1(b = 0. 001 in this paper) is introduced in theobjective function of the outer problem in order to exclude the combinations containingredundant inhibitions and therefore avoiding an “over-selection” ofdrugs.

To solve this bilevel optimization we apply the strong duality theorem which consists inappending a list of constraints corresponding to the dual of the inner problem andsetting the primal objective function equal to the dual [58, 59]. This leads to a single minimization problem. By calling the dual variablesas follows

$$\begin{array}{ll}{\mu }_1,\dots ,{\mu }_{N_m}\in \mathbb{R}:& \text{associated to the first}{N}_m\\ \text{constraints in the inner problem};\\ {\lambda }_1,\dots ,{\lambda }_{N_r}\in {\mathbb{R}}_{+}:& \text{associated to the second set of the}\\ \text{constraints of the inner problem};\\ {\delta }_1,\dots ,{\delta }_{N_t}\in {\mathbb{R}}_{+}:& \text{associated to drugs targets (third}\\ \text{set of constraints;}{N}_t:=\sum _{j=1}^{N_r}{N}_{t,j}),\end{array}$$

the final optimization problem becomes:

Minimize

$$\sum _{j=1}^{N_r}{\alpha }_j(1-y_j)+\sum _{k=1}^{N_d}{\beta }_k(1-d_k)-b\sum _{k=1}^{N_d}{d}_k$$

such that

$$\begin{array}{lc}\hfill \sum _{j=1}^{N_r}{S}_{i,j}{v}_j=0& \forall i\in \mathcal{M};\hfill \\ \hfill v_j\le U_j& \forall j\in \mathcal{R};\hfill \\ \hfill v_j\le U_j{d}_k& \forall j\in \mathcal{R},k\in {\mathcal{T}}_j;\hfill \\ \sum _{i=1}^{N_m}{S}_{i,j}{\mu }_i+{\lambda }_j+\sum _{i=1}^{N_{t,j}}{\delta }_i\ge 0& \forall j\in \mathcal{R}\setminus \left\{\text{obj}\right\};\hfill \\ \sum _{i=1}^{N_m}{S}_{i,j}{\mu }_i+{\lambda }_j+\sum _{i=1}^{N_{t,j}}{\delta }_i\ge 1;& \text{for}j=\text{obj};\hfill \end{array}$$

$$\begin{array}{l}{v}_{\mathit{\text{obj}}}=\sum _{i=1}^{N_r}{U}_i\left({\lambda }_i+{\delta }_i\sum _{k\in {\mathcal{T}}_i}{d}_k\right);\\ \hfill v_{\mathit{\text{obj}}}=0;\\ \hfill ϵy_j\le v_j& \forall j\in \mathcal{R};\hfill \\ \hfill U_j{y}_j\ge v_j& \forall j\in \mathcal{R};\hfill \\ \hfill y_j+y_l\le 1& \forall j\in \mathcal{R}\text{such that}{v}_l\text{is}\hfill \\ \text{the opposite flux of}{v}_j.\hfill \end{array}$$

(3)

The key simplification is that the nonlinear terms δ _i d _k  =: z _ik in (3) (the strong duality theorem equality) are exactly linearizable as follows:

$$0\le z_{\mathrm{ik}}\le {\delta }_i^{\max }{d}_k$$

(4)

$${\delta }_i-{\delta }_i^{\max }(1-d_k)\le z_{\mathrm{ik}}\le {\delta }_i$$

(5)

where ${\delta }_i^{\max }$ is the upper bound for the dual variable δ _i (chosen arbitrarily big in the implementation).

The codes for the algorithm have been developed in MATLAB (MathWorks R2010b) and areavailable in Additional files 4, 5,6, 7, 8and 9. All Mixed Integer Linear Optimizations have beenperformed using the ILOG-IBM CPLEX 12.1, under free academic license.

Selection procedure for the drugs

In order to generate realistic solutions, drugs have been accurately selected fromhttp://www.drugbank.ca[26] according to the following procedure (the output of the query issummarized in Additional file 1: Table S1):

• 1. the whole database contains 6708 drugs;

• 2. only approved drugs have been picked out, restricting the search to 1570 drugs;

• 3. we filter for drugs which act on human enzymes (identified by the EC number): the set reduces to 473 drugs;

• 4. we select only drugs for which an inhibitory effect on at least one enzyme of the human metabolism has been experimentally proven. The set reduces to 267 drugs (with the EC numbers of the inhibited enzymes). For these drugs, we count the number of non-metabolic targets, including also cases of agonism, antagonism and activation, since they all represent a perturbation to the regular functioning of the target (activation is not considered for metabolic targets because it does not affect the number of stopped reaction in our FBA formulation).

• 5. the reactions of the metabolic network directly inhibited by each drug are identified through the available correspondence between EC numbers of the inhibited enzymes and the gene codes first, and then through the correspondence between gene codes and metabolic reactions. During this step, it may happen that many genes, and hence many reactions, are associated to the same EC number. Therefore, although in the original database a drug inhibits only a single or a few targets, the number of metabolic reactions affected by the drug can be high. For instance, this is the case of Rosiglitazone which inhibits a single target, the long-chain-fatty-acid-CoA ligase an enzyme responsible for the binding of the acyl-CoA group to a long fatty acid chain. Since the substrate of this enzyme can be any carbon chain, regardless of the unsaturation (presence or not of double bonds C=C) and of the exact length (it just requires a chain longer than 12 carbons atoms), we end up with a drug which inhibits up to 60 metabolic targets.

• 6. drugs which have exactly the same metabolic targets are grouped; the final 85 groups are listed in Additional file 1: Table S1. For each group only a representative is reported, namely the drug which has the minimal number of targets outside the metabolism. The full list of drugs in each group is reported in the Additional file 2.