Network-principled deep generative models for designing drug combinations as graph sets.

MOTIVATION: Combination therapy has shown to improve therapeutic efficacy while reducing side effects. Importantly, it has become an indispensable strategy to overcome resistance in antibiotics, antimicrobials and anticancer drugs. Facing enormous chemical space and unclear design principles for small-molecule combinations, computational drug-combination design has not seen generative models to meet its potential to accelerate resistance-overcoming drug combination discovery.
RESULTS: We have developed the first deep generative model for drug combination design, by jointly embedding graph-structured domain knowledge and iteratively training a reinforcement learning-based chemical graph-set designer. First, we have developed hierarchical variational graph auto-encoders trained end-to-end to jointly embed gene-gene, gene-disease and disease-disease networks. Novel attentional pooling is introduced here for learning disease representations from associated genes' representations. Second, targeting diseases in learned representations, we have recast the drug-combination design problem as graph-set generation and developed a deep learning-based model with novel rewards. Specifically, besides chemical validity rewards, we have introduced novel generative adversarial award, being generalized sliced Wasserstein, for chemically diverse molecules with distributions similar to known drugs. We have also designed a network principle-based reward for disease-specific drug combinations. Numerical results indicate that, compared to state-of-the-art graph embedding methods, hierarchical variational graph auto-encoder learns more informative and generalizable disease representations. Results also show that the deep generative models generate drug combinations following the principle across diseases. Case studies on four diseases show that network-principled drug combinations tend to have low toxicity. The generated drug combinations collectively cover the disease module similar to FDA-approved drug combinations and could potentially suggest novel systems pharmacology strategies. Our method allows for examining and following network-based principle or hypothesis to efficiently generate disease-specific drug combinations in a vast chemical combinatorial space.
AVAILABILITY AND IMPLEMENTATION: https://github.com/Shen-Lab/Drug-Combo-Generator. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

1 Introduction

Drug resistance is a fundamental barrier to developing robust antimicrobial and anticancer therapies (Housman ; Taubes, 2008). Its first sign was observed in 1940s soon after the discovery of penicillin (Abraham and Chain, 1940), the first modern antibiotic. Since then, drug resistance has surfaced and progressed in infectious diseases such as HIV (Clavel and Hance, 2004), tuberculosis (TB) (Dooley ) and hepatitis (Ghany and Liang, 2007) as well as cancers (Holohan ). Mechanistically, it can emerge through drug efflux (Chang and Roth, 2001), activation of alternative pathways (Lovly and Shaw, 2014) and protein mutations (Balbas ; Toy ) while decreasing the efficacy of drugs.

Combination therapy is a resistance-overcoming strategy that has found success in combating HIV (Shafer and Vuitton, 1999), TB (Ramón-García ), cancers (Bozic ; Sharma and Allison, 2015) and so on. Considering that most diseases and their resistances are multifactorial (Kaplan and Junien, 2000; Keith ), multiple drugs targeting multiple components simultaneously could confer less resistance than individual drugs targeting components separately. Examples include targeting both mitogen-activated protein kinase kinase (MEK) and B-Raf proto-oncogene serine/threonine-protein kinase (BRAF) in patients with BRAF V600-mutant melanoma rather than targeting MEK or BRAF alone (Flaherty ; Madani Tonekaboni ). The effect of drug combination is usually categorized as synergistic, additive or antagonistic depending on whether it is greater than, equal to or less than the sum of individual drug effects (Chou, 2006). Synergistic combinations are effective at delaying the beginning of the resistance; however, antagonistic combinations are effective at suppressing expansion of resistance (Saputra ; Singh and Yeh, 2017), representing offensive and defensive strategies to overcome drug resistance. In particular, offensive strategies cause huge early causalities but defensive ones anticipate and develop protection against future threats (Saputra ).

Discovering a drug combination to overcome resistance is, however, extremely challenging, even more so than discovering a drug which is already a costly (∼billions of USD) (DiMasi ) and lengthy (∼12 years) (Van Norman, 2016) process with low success rates (3.4% phase-1 oncology compounds make it to approval and market) (Wong ). An apparent challenge, a combinatorial one, is in the scale of chemical space, which is estimated to be 1060 for single compounds (Bohacek ) and can ‘explode’ to for K-compound combinations. Even if the space is restricted to around 103 U.S. Food and Drug Administration (FDA)-approved human drugs, there are 105–106 pairwise combinations. Another challenge, a conceptual one, is in the complexity of systems biology. On top of on-target efficacy and off-target side effects or even toxicity that need to be considered for individual drugs, network-based design principles are much needed for drug combinations that effectively target multiple proteins in a disease module and have low toxicity or even resistance profiles (Billur Engin ; Martínez-Jiménez and Marti-Renom, 2016).

Current computational models in drug discovery, especially those for predicting pharmacokinetic and pharmacodynamic properties of individual drugs/compounds, can be categorized into discriminative and generative models. Discriminative models predict the distribution of a property for a given molecule, whereas generative models would learn the joint distribution on the property and molecules. For instance, discriminative models have been developed for predicting single compounds’ toxicities, based on support vector machines (Darnag ), random forest (Svetnik ) and deep learning (Mayr ). Whereas discriminative models are useful for evaluating given compounds or even searching compound libraries, generative models can effectively design compounds of desired properties in chemical space. Recent advance in inverse molecular design has seen deep generative models such as SMILES representation-based reinforcement learning (RL) (Popova ) or recurrent neural networks as well as graph representation-based generative adversarial networks (GANs), RL (You ) and generative tensorial RL (Zhavoronkov ).

Unlike single-drug design, current computational efforts for drug combinations are exclusively focused on discriminative models and lack generative models. The main focus for drug combination is to use discriminate models to identify synergistic or antagonistic drugs for a given specific disease. Examples include the Chou–Talalay method (Chou, 2010), integer linear programming (Pang ) and deep learning (Preuer ). However, it is daunting if not infeasible to enumerate all cases in the enormous chemical combinatorial space and evaluate their combination effects using a discriminative model. Not to mention that such methods often lack explainability.

Directly addressing aforementioned combinatorial and conceptual challenges and filling the void of generative models for drug combinations, in this study, we develop network-based representation learning for diseases and deep generative models for accelerated and principled drug combination design (the general case of K drugs). Recently, by analyzing the network-based relationships between disease proteins and drug targets in the human protein–protein interactome, Cheng et al. proposed an elegant principle for FDA-approved drug combinations that targets of two drugs both hit the disease module but cover different neighborhoods. Our methods allow for examining and following the proposed network-based principle (Cheng ) to efficiently generate disease-specific drug combinations in a vast chemical combinatorial space. They will also help meet a critical need of computational tools in a battle against quickly evolving bacterial, viral and tumor populations with accumulating resistance.

To tackle the problem, we have developed a network principle-based deep generative model for faster, broader and deeper exploration of drug combination space by following the principle underling FDA-approved drug combinations. First, we have developed hierarchical variational graph auto-encoders (HVGAE) for jointly embedding disease–disease network and gene–gene network. Through end-to-end training, we embed genes in a way that they can represent the human interactome. Then, we utilize their embeddings with novel attentional pooling to create features for each disease so that we can embed diseases more accurately. Second, we have also developed a reinforcement-learning-based graph-set generator for drug combination design by utilizing both gene/disease embedding and network principles. Besides those for chemical validity and properties, our rewards also include (i) a novel adversarial reward, generalized sliced Wasserstein distance (GSWD), that fosters generated molecules to be diverse yet similar in distribution to known compounds (ZINC database and FDA-approved drugs); and (ii) a network principle-based reward for drug combinations that are feasible for online calculations.

The overall schematics are shown in Figure 1, and they are detailed in Section 3.

Fig. 1.

Overall schematics of the proposed approach for generating disease-specific drug combinations

2 Data

2.1 Human interactome and its features

We used the human interactome data (a gene–gene network) from Menche that feature 13 460 proteins interconnected by 141 296 interactions.

We introduced edge features for the human interactome based on the biological nature of edges (interactions). The interactome was compiled by combining experimental support from various sources/databases including (i) regulatory interactions from TRANSFAC (Matys ); (ii) binary interactions from high-throughput (including Rolland ) and literature-curated datasets [including IntAct (Aranda ) and MINT (Ceol )] as well as literature-curated interactions from low-throughput experiments [IntAct, MINT, BioGRID (Stark ), and HPRD (Keshava Prasad )]; (iii) metabolic enzyme-coupled interactions from Lee ; (iv) protein complexes from CORUM (Ruepp ); (v) kinase-substrate pairs from PhosphositePlus (Hornbeck ); and (vi) signaling interactions. In summary, an edge could correspond to one or multiple physical interaction types. So we used a six-hot encoding for edge features, based on whether an edge corresponds to regulatory, binary, metabolic, complex, kinase and signaling interactions.

We also introduced features for nodes (genes) in the human interactome based on (i) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa ) (336 features) queried through Biopython (Cock ); (ii) Gene Ontology (GO) terms (Ashburner ) including biological process (30 769 features), molecular function (12 183 features) and cellular component (4451 features), mapped using the NCBI Gene2Go dataset; (iii) disease–gene associations from the database OMIM (Mendelian Inheritance in Man) (Hamosh et al., 2005) and the results from genome-wide association studies (Mottaz ; Ramos ) (299 features). The last 299 features correspond to 299 diseases represented by the Medical Subject Headings (MeSH) vocabulary (Rogers, 1963).

After removing those genes without KEGG pathway information, the human interactome used in this study has 13 119 genes and 352 464 physical interactions.

2.2 Disease–disease network

We used a disease–disease network from Menche with 299 nodes (diseases), created based on human interactome data (as detailed earlier), gene expression data (Su ), disease–gene associations (Hamosh ; Mottaz ; Ramos ), GO (Ashburner ), symptom similarity (Zhou ) and comorbidity (Hidalgo ). The original disease–disease network is a complete graph with real-valued edges. The edge value between two diseases shows how much they are topologically separated from each other. A positive/negative edge weight indicates that that two disease modules are topologically separated/overlapped. Therefore, we used zero-weight as the threshold and pruned positive-valued edges, which results in a disease–disease network of 299 nodes and 5986 edges (without weights).

2.3 Disease–gene associations

We used disease–gene associations from the database OMIM (Hamosh ). These associations bridge aforementioned gene–gene and disease–disease networks into a hierarchical graph of genes and diseases, based on which gene and disease representations will be learned.

2.4 Disease classification

For the purpose of assessment, we used the Comparative Toxicogenomics Database (CTD) (Davis ) to classify diseases into eight classes based on their Disease Ontology (DO) terms (Schriml ) where diseases are represented in the MeSH vocabulary (Rogers, 1963). In the CTD database only 201 of the 299 diseases have a corresponding DO term. Therefore, for the 98 diseases with missing DO terms we considered the majority of their parents’ DO terms, if applicable, as their DO terms. With this approach, we assigned DO terms to 66 such diseases and classified 267 of the 299 diseases. The 32 diseases with DO terms still missing are usually at the top layers of the MeSH tree.

2.5 FDA-approved drugs and drug combinations

To assess our deep generative model for drug combination design (to be detailed in Section 3.2), we consider a comprehensive list of US FDA-approved combination drugs (1940–2018.9) (Das ). The dataset contains 419 drug combinations consisting of 328 unique drugs, including 341 (81%), 67 (16%) and 11 (3%) of double, triple and quadruple drug combinations.

We also utilized the curated drug–disease association from CTD database (Davis ).

3 Materials and Methods

We have developed a network-based drug combination generator which can be utilized in overcoming drug resistance. Representing drugs through their molecular graphs, we recast the problem of drug combination generation into network-principled graph-set generation by incorporating prior knowledge such as human interactome (gene–gene), disease–gene, disease–disease, gene pathway and gene–GO relationships. Furthermore, we formulate the graph-set generation problem as learning an RL agent that iteratively adds substructures and edges to each molecular graph in a chemistry- and system-aware environment. To that end, the RL model is trained to maximize a desired property Q (e.g. therapeutic efficacy for drug combinations) while following the valency (chemical validity) rules and being similar in distribution to the prior set of graphs.

As shown in Figure 1, the proposed approach consists of (i) embedding prior knowledge (different network relationships) through HVGAE and (ii) generating drug combinations as graph sets through an RL algorithm, which will be detailed next.

Notations: As both gene–gene and disease–disease networks can be represented as graphs, notations are differentiated by superscripts ‘g’ and ‘d’ to indicate gene–gene and disease–disease networks, respectively. Drugs (compounds) are also represented as graphs and notations with ‘k’ in the superscript indicate the kth drug (graph) in the drug combination (graph set).

3.1 HVGAE for representation learning

Suppose that a gene–gene network is represented as a graph , where is the adjacency tensor of the gene–gene network with nodes and n edge types (k-hot encoding of six types of aforementioned physical interactions such as regulatory, binary, metabolic, complex, kinase and signaling interactions). We also define to be elementwise of . Furthermore, denotes the mth set of node features for gene–gene network where M (5 in the study) represents different types of node features such as pathways, three GO terms and gene–disease relationship. We also suppose that the disease–disease network is represented as graph , where is the adjacency matrix of the disease–disease network with nodes; and represents the set of node features for the disease–disease network.

We have developed a hierarchical embedding with two levels. In the first level, we embed the gene–gene network to get the features related to each disease and then we incorporate the disease features within the disease–disease network to embed their relationship. We infer the embedding for each gene and disease jointly through end-to-end training. The proposed HVGAE perform probabilistic auto-encoding to capture uncertainty of representations which is in the same spirit as the variational graph auto-encoder models introduced in (Hajiramezanali ; Hasanzadeh ; Kipf and Welling, 2016).

3.1.1 First level: gene–gene embedding

The inference model for variational embedding of the gene–gene network is formulated as follows. We first use M graph neural networks (GNNs) to transform individual nodes’ features in M types and then concatenate the M sets of results () into : where is an aggregation function combining output features of s for each node. We used a two layer fully connected neural network with ReLU activation functions followed by a single linear layer in our implementation. We then approximate the posterior distribution of stochastic latent variables [containing for where (32 in this study) is the latent space dimensionality for the ith gene], with a multivariate Gaussian distribution given the gene–gene network’s aggregated node features and adjacency tensor : where ; is the matrix of mean vectors ; and is the matrix of standard deviation vectors ().

The generative model for the gene–gene network is formulated as: and is the logistic sigmoid function. The loss for gene–gene variational embedding is represented as a variational lower bound (ELBO): where is the Kullback–Leibler divergence between and . We take the Gaussian prior for and make use of the reparameterization trick (Kipf and Welling, 2016) for training.

3.1.2 Second level: disease–disease embedding

The inference model for variational embedding of the disease–disease network is similar to that of the gene–gene network except that the disease–disease network’s aggregated node features, , are derived through parameterized attentional pooling of , latent variables of genes associated with the rth disease (a subset of ): where capture the importance of genes related to the rth disease for calculating its latent representations and is the latent space dimensionality of a disease.

Once , the disease–disease network’s aggregated node features for all diseases, are derived; we again define for the posterior distribution of stochastic latent variables similarly to what we did in Equation (2) except that AGG functions are removed since disease–disease network has one binary adjacency matrix; give the generative decoder for embedding the disease–disease network similarly to what we did in Equation (3); and calculate the variational lowerbound (ELBO) loss for the disease–disease network similarly to what we did in Equation (4). Details can be found in Supplementary Section S1.1.

Both levels of our proposed HVGAE, i.e. gene–gene and disease–disease variational graph representation learning, are jointly trained in an end-to-end fashion using the following overall loss:

3.2 RL-based graph-set generator for drug combinations

In this section, we introduce the RL-based drug combination generator. We will detail (i) the state space of graph sets (K compounds) and the action space of graph-set growth; (ii) multi-objective rewards including chemical validity and our generalized sliced Wasserstein reward for individual drugs as well as our newly designed network principle-based reward for drug combinations; and (iii) policy network that learns to take actions in the rewarding environment.

3.2.1 State and action space

We represent a graph set (drug combination) with K graphs as . Each graph where is the adjacency matrix, the node feature matrix, the edge-conditioned adjacency tensor and n the number of vertices for the kth graph; and is the number of features per nodes and ϵ the number of edge types.

The state space is the set of all K graphs with different numbers and types of nodes or edges. Specifically, the state of the environment s at iteration t is defined as the intermediate graph set generated so far which is fully observable by the RL agent.

The action space is the set of edges that can be added to the graph set. An action a at iteration t is analogous to link prediction in each graph in the set. More specifically, a link can either connect a new subgraph (a single node/atom or a subgraph/drug-substructure) to a node in or connect existing nodes within graph . The actions can be interpreted as connecting the current graph with a member of scaffold subgraphs set C. Mathematically, for , graph k at step t, the action is the quadruple of .

3.2.2 Multi-objective reward

We have defined a multi-objective reward R to satisfy certain requirements in drug combination therapy. First, a chemical validity reward maintains that individual compounds are chemically valid. Second, a novel adversarial reward, generalized sliced Wasserstein GAN (GS-WGAN), enforces that generated compounds are synthesizable and ‘drug-like’ by following the distribution of synthesizable compounds in the ZINC database (Irwin and Shoichet, 2005) or FDA-approved drugs. Third, a network principle-based award would encourage individual drugs to target the desired disease module but not to overlap in their target sets. Toxicity due to drug–drug interactions (DDIs) can also be included as a reward. It is intentionally left out in this study so that toxicity can be evaluated for drug combinations designed to follow the network principle.

When training the RL agent, we use different reward combinations in different stages. We first only use the weighted combination of chemical validity and GS-WGAN awards learning over drug combinations for all diseases; then we remove the penalized logP (Pen-logP) portion of chemical validity and add adversarial loss again while learning over drug combinations for all diseases; and finally use the combination of the three rewards as in the second stage but focusing on a target disease and possibly on restricted actions/scaffolds (in a spirit similar to transfer learning). The three types of rewards are detailed as follows.

Chemical validity reward for individual drugs. A small positive reward is assigned if the action does not violate valency rules. Otherwise a small negative reward is assigned. This is an intermediate reward added at each step. Another reward is on penalized logP (lipophilicity where P is the octanol–water partition coefficient) or Pen-logP values. The design and the parameters of this reward are adopted from You without optimization.

Adversarial reward using GSWD. To ensure that the generated molecules resemble a given set of molecules (such as those in ZINC or FDA-approved), we deploy GANs. GANs are very successful at modeling high-dimensional distributions from given samples. However, they are known to suffer from training instability and cannot generate diverse samples (a phenomenon known as mode collapse).

Wasserstein GANs (WGAN) have shown to improve stability and mode collapse by replacing the Jenson–Shannon divergence in original GAN formulation with the Wasserstein distance (WD) (Arjovsky ). More specifically, the objective function in WGAN with gradient penalty (Gulrajani ) is defined as follows: where p is the data distribution, λ is a hyper-parameter, R is the Lipschitz continuity regularization term, is the critic with parameters , and is the policy (generator) with parameters θ.

Despite theoretical advantages of WGANs, solving Equation (7) is computationally expensive and intractable for high-dimensional data. To overcome this problem, we propose and formulate a novel GS-WGAN which deploys GSWD (Kolouri ). GSWD, first, factorizes high-dimensional probabilities into multiple marginal 1D distributions with generalized Radon transform (GRT). Then, by taking advantage of closed form solution of WD in 1D, the distance between two distributions is approximated by the sum of WDs of marginal 1D distributions. More specifically, let represent GRT operator. The GRT of a probability distribution which is defined as follows: where is the 1D Dirac delta function, is a scalar, ψ is a unit vector in the unit hyper-sphere in a d-dimensional space (), and f is a projection function whose parameters will be learned in training. Injectivity of the GRT (Beylkin, 1984) is the requirement for the GSWD to be a valid distance. We use linear project here and can easily extend to two nonlinear cases that maintain the GRT-injectivity (circular nonlinear projections or homogeneous polynomials with an odd degree).

GSWD between two d-dimensional distributions and is therefore defined as:

The integral in the above equation can be approximated with a Riemann sum. Knowing the definition of GSWD, we define the objective function of GS-WGAN as follows: where the parameters and notations are the same as defined in Equation (7).

We note that x and y in Equation (10) are random variables in , which is not a reasonable assumption for graphs. To that end, we use an embedding function g that maps each graph to a vector in . We use graph convolutional (GC) layers followed by fully connected layers to implement g. We deploy the same type of neural network architecture for . We use as the adversarial reward used together with other rewards and optimize the total rewards with a policy gradient method (Section 3.2.3).

Network principle-based reward for drug combinations. Proteins or genes associated with a disease tend to form a localized neighborhood disease module rather than scattering randomly in the interactome (Cheng ). A network-based score has been introduced (Menche ) to efficiently capture the network proximity of a drug (X) and disease (Y) based on the shortest-path length d(x, y) between a drug target (x) and a disease protein (y): where is the shortest path distance; and σ are the mean and standard deviation of the reference distribution which is corresponding to the expected network topological distance between two randomly selected groups of proteins matched to size and degree (connectivity) distribution as the original disease proteins and drug targets in the human interactome. Z-score being negative (Z < 0) implies network proximity of disease module and drug targets which is desirable. From the drug combination perspective, it has been shown that the drug–drug relationship of complementary exposure has the least side effect and the most clinical efficacy (Cheng ). Complementary exposure for two drugs (X1 and X2) means that the X1 drug targets (x1) and X2 drug targets (x2) are not in the same neighborhood and have the least overlapping. Therefore, Cheng have proposed a network-separation score which is formulated as follows: where is the mean shortest path distance between drugs X1 and X2; and are the mean shortest path distance within drug targets X1 and X2, respectively. The separation score being positive (s > 0) implies that two networks are separated from each other, which is desirable. We have extended and combined these scores for general drug combination therapy where we have a set of k drugs and disease Y:

However, the exact online calculation of the reward is infeasible while training across all the diseases and the whole human interactome with more than 13K nodes and 352K edges. Therefore, we have developed a relaxed version of the reward which is feasible for online calculation and correlates with the actual reward. Specifically, we consider the normalized exclusive or (XOR) of intersections of disease modules with drug targets:

The relaxed network principle-based reward is penalizing a drug combination if the overlap between drug targets in the disease module is high; therefore, it will prevent the adverse DDIs. We scaled the network score by a constant (equals 10) such that the score would be in the same range as Pen-logP and can use the same weight in the total reward as Pen-logP did in You .

For a generated compound, we predict its protein targets by DeepAffinity (Karimi ), judging by whether the predicted IC50 is below 1 μM.

3.2.3 Policy network

Having explained the graph generation environment (various rewards), we outline the architecture of our proposed policy network. Our method takes the intermediate graph set and the collection of scaffold subgraphs C as inputs, and outputs the action a, which predicts a new link for each of the graphs in (You ).

Since the input to our policy network is a set of K compounds or graphs , we first deploy some layers of graph neural network to process each of the graphs. More specifically, where is a multilayer-graph neural network. The link prediction-based action at iteration t is a concatenation of four components for each of the K graphs: selection of two nodes, prediction of edge type, and prediction of termination. Each component is sampled according to a predicted distribution (You ). Details are included in the Supplementary Section S1.2. We note that the first node is always chosen from while the next node is chosen from . We also note that infeasible actions (i.e. actions that do not pass valency check) proposed by the policy network are rejected and the state remains unchanged. We adopt proximal policy optimization (PPO) (Schulman ), one of the state-of-the-art policy gradient methods, to train the model.

4 Results

To assess the performance of our proposed model, we have designed a series of experiments. In Section 4.1, we first compare HVGAE to state-of-the-art graph embedding methods in disease–disease network representation learning and further include several variants of HVGAE for ablation studies. We then assess the performance of the proposed RL method in two aspects. In a landscape assessment in Section 4.2, we examine designed pairwise compound combinations for 299 diseases in quantitative scores of following a network-based principle (Cheng ). In Section 4.3, we focus on four case studies involving multiple diseases of various systems pharmacology strategies. Our method is capable of generating higher-order combinations of K drugs. As FDA-approved drug combinations are often pairs, here we design compound pairs from the scaffolds of FDA-approved drug pairs. We further delve into designed compound pairs to understand the benefit of following network principles in lowering toxicity from DDIs. We also do so to understand their systems pharmacology strategies in comparison to the FDA-approved drug combinations.

4.1 HVGAE representation compares favorably to baselines

4.1.1 Experiment setup

To assess the performance of our proposed embedding method HVGAE, we compare its performance in (disease–disease) network reconstruction with Node2Vec (Grover and Leskovec, 2016), DeepWalk (Perozzi ) and VGAE (Kipf and Welling, 2016), as well as some variants of our own model for ablation study. Node2Vec and DeepWalk are random walk-based models that do not capture node attributes; hence, we only used the disease–disease graph structure. For VGAE, we used identity matrix as node attributes as suggested by the authors.

For our HVGAE described in Section 3.1, we also considered two variants for ablation study: HVGAE-disjoint does not jointly embed gene–gene and disease–disease networks and does not use attentional pooling for disease embedding, whereas HVGAE-noAtt just does not use attentional pooling. Specifically, in HVGAE-disjoint, we, first, learned an embedding for gene–gene network, then used the sum of mean of the node representations of genes affected by a disease as its node attributes. In HVGAE-noAtt, we jointly learned the representations while using sum of mean of the node representations of genes as node attributes for disease–disease network.

In node2vec and DeepWalk, the walk length was set to 80, the number of walks starting at each node was set to 10, and the nodes were embedded to a 16D space. The window size was 10 for node2vec while it is set to 10 in DeepWalk. All models were trained using Adam optimizer. In VGAE, a 32D GC layer followed by two 16D layers was used for mean and variance inference. The learning rate was set to 0.01.

For HVGAE and its variants (for ablation study), we embed gene networks in 32D space using a single GC layer with 32 filters for each of the five types of input followed by a 64D GC layer and two 32D GC layer to infer mean and variance of the representation. We used a single 32D fully connected layer for attention layer. For disease–disease network embedding, we deployed a single 32D GC layer followed by two 16D layer for mean and variance inference resulting in 16D embedding for disease–disease network. Learning rates were set to 0.001. The models were trained for 1000 epochs choosing the best representation based on their reconstruction performance at each epoch.

4.1.2 Numerical analysis and ablation study for network embedding

Table 1 summarizes the reconstruction performance of the aforementioned methods. Compared to all baselines, our HVGAE showed the best performance in all metrics considered. Node2Vec and DeepWalk showed the worst performance as they only use the graph structure. The performance of VGAE was very close to DeepWalk. This is due to the fact that no attributes have been provided to VGAE despite having the capability of capturing attributes.

Table 1.

Graph reconstruction performances (unit: %) in the disease–disease network using our proposed HVGAE and baselines

Method	AUC-ROC	AP	F1-Macro	F1-Micro
Node2Vec	79.01	72.82	35.73	51.10
DeepWalk	79.32	73.77	40.28	53.30
VGAE	88.12	85.71	60.19	64.98
HVGAE-disjoint	91.45	90.72	73.45	74.77
HVGAE-noAtt	92.83	92.34	73.81	75.14
HVGAE	96.11	95.89	79.77	80.45

Note: F1 scores are based on 50% threshold. Method names in bold are our HVGAE and its variants. Numbers in bold correspond to the best performances.

Compared to VGAE, HVGAE-disjoint without joint embedding or attentional pooling still saw better performance, which suggests that the attribute generated by the gene–gene network contains meaningful features about the disease–disease network. The slight performance gain from HVGAE-disjoint to HVGAE-noAtt shows that joint learning of both networks hierarchically helps to render more informative features for the disease–disease network. Finally, HVGAE had another performance boost compared to HVGAE-noAtt and outperformed all competing methods, which shows the benefit of attentional pooling. Specifically, the attention layer of HVGAE allows the model to produce features that are specifically informative for the disease–disease network representation learning.

4.2 Our model generates drug combinations following network principles across diseases

4.2.1 Experiment setup

We have trained the proposed reinforcement model in three stages using different rewards, disease sets and action spaces to increasingly focus on a target disease while exploiting all diseases whose representations already jointly embed gene–gene, disease–disease and gene–disease networks. In the first stage, we train the model to only generate drug-like small molecules which follow the chemistry valency reward, lipophilicity reward (Pen-logP where P is the octanol–water partition coefficient) (You ), and our novel adversarial reward for individual compounds. In this study, we trained the model for 3 days (4800 iterations) to learn to follow the valency conditions and promote high logP for generated compounds.

In the second stage, we start from the trained model at the end of the first stage (‘warm-start’ or ‘pre-training’). And we continue to train the model to generate good drug combinations across all diseases. We do so by adding the network principle-based reward for compound combinations and sequentially generating drug combinations for each disease one by one. Then, we calculate the network-based score for the generated drug combinations at the last epoch across disease ontologies and compare them with the FDA-approved melanoma drug combinations’ network-based score. In this study, we trained the model for 1500 iterations to generate drug combinations across all 299 diseases. In each iteration, we generated eight drug combinations for a given disease. We adopted PPO (Schulman ) with a learning rate of 0.001 to train the proposed RL for both stages.

The last stage is disease-specific and is detailed in Section 4.3.

4.2.2 Numerical analysis

Across disease ontologies we quantify the performance of the proposed RL (stage 2 model first) using quantitative scores of compound combinations following a network-based principle (Cheng ). We consider the generated combinations in the last epoch (the last 299 iterations) and calculate the network score based on disease ontologies. We asses our model based on two versions of disease classification, original disease ontology and its extension, explained in Section 2.4. Table 2 summarizes the network-based scores for our model. Specifically, suppose that the set of targets for drugs 1 and 2 are represented by A and B whereas the disease module is the universal set Ω, we report the portion exclusively covered by drug 1 (), exclusively covered by drug 2 (), overlapped by both () and collectively by both (). As a reference, we calculated the corresponding network scores for three FDA-approved drug combinations for melanoma.

Table 2.

Network-based score for the generated drug combinations based on disease ontology classifications

	Disease Ontology				Disease Ontology extended
	η_A−B	η_B−A	η_A∩B	η_A∪B	η_A−B	η_B−A	η_A∩B	η_A∪B
Infectious disease	0.25	0.10	0.06	0.41	0.20	0.07	0.05	0.33
Disease of anatomical entity	0.27	0.12	0.10	0.49	0.26	0.11	0.09	0.48
Disease of cellular proliferation	0.25	0.09	0.07	0.42	0.25	0.10	0.08	0.44
Disease of mental health	0.22	0.11	0.10	0.43	0.22	0.11	0.10	0.43
Disease of metabolism	0.22	0.13	0.10	0.46	0.23	0.14	0.11	0.48
Genetic disease	0.23	0.15	0.11	0.4	0.23	0.15	0.11	0.49
Syndrome	0.22	0.11	0.11	0.44	0.22	0.11	0.11	0.44

Based on the results shown in Table 2, we note that across all disease classes, the designed compound combinations learned in an environment, where the network principle(Cheng ) was rewarded, did achieve the desired performances. Specifically, their overlaps in disease modules were low as fractions were around 0.1; whereas their joint coverage in disease modules was high as fractions were in the range of 0.4–0.5 for all diseases.

Compared to a few FDA-approved drugs for melanoma in Table 3, we notice that the designed compound combinations had similar exclusive coverage ( and ) as the drug combinations. However, the overlapping and overall coverage ( and ) were both much higher in FDA-approved drug combinations than the designed. Improvements could be made by training the RL agent longer, as these scores had already been improving during the limited training process under computational restrictions. More improvement can be made by adjusting the network-based reward as well.

Table 3.

Network-based scores for FDA-approved melanoma drug combinations

	η_A−B	η_B−A	η_A∩B	η_A∪B
Dabrafenib + Trametinib	0.05	0.21	0.55	0.81
Encorafenib + Binimetinib	0.21	0.05	0.53	0.86
Vemurafenib + Cobimetinib	0.05	0.27	0.36	0.68

4.3 Case studies for specific diseases

4.3.1 Experiment setup

In the third and last stage of RL model training, we start from the stage 2 model and generate drug combinations for a fixed target disease and can choose scaffold libraries specific to the disease. In parallel, we trained the model for 500 iterations (roughly 1 day) to generate 4000 drug combinations specifically for each of four diseases featuring various drug-combination strategies: melanoma, lung cancer, ovarian cancer and breast cancer. In all cases, we started with the Murcko scaffolds of specific FDA-approved drug combinations to be detailed next.

Melanoma: Different targets in the same pathway. Resistance to BRAF kinase inhibitors is associated with reactivation of the mitogen-activated protein kinase (MAPK) pathway. There is, thus, a phase 1 and 2 trial of combined treatment with Dabrafenib, a selective BRAF inhibitor, and Trametinib, a selective MAPK kinase (MEK) inhibitor. As melanoma is not one of the 299 diseases, we chose broader neoplasm as an alternative. To compensate the loss of focus on target disease, we design compound pairs from Murcko scaffolds of Dabrafenib + Trametinib.

Lung and ovarian cancers: targeting parallel pathways. MAPK and phosphoinositide 3-kinase (PI3K) signaling pathways are parallels important for treating many cancers, including lung and ovarian cancers (Bedard ; Day and Siu, 2016). Clinical data suggest that dual blockade of these parallel pathways has synergistic effects. Buparlisib (BKM120) and Trametinib (GSK1120212; Mekinist) are used as a drug combination therapy for the purpose. Specifically, Buparlisib is a potent and highly specific PI3K inhibitor, whereas Trametinib is a highly selective, allosteric inhibitor of MEK1/MEK2 activation and kinase activity (Bedard ).

Breast cancer: reverse resistance. Endocrine therapies, including Fulvestrant, are the main treatment for hormone receptor-positive breast cancers (80% of breast cancers) (Turner ). However, they could confer resistance to patients during or after the treatment. A Phase 3 study is using Fulvestrant and Palbociclib as a combination therapy to reverse the resistance. Fulvestrant and Palbociclib are targeting different genes in different pathways. Specifically, Fulvestrant targets estrogen receptor α in estrogen signaling pathway and Palbociclib targets cyclin-dependent kinases 4 and 6 (CDK4 and CDK6) in cell cycle pathway (Turner ).

4.3.2 Baseline methods for drug pair combination

Since our proposed method is the first to generate drug combinations for specific diseases, we consider the following baseline methods to compare with: (i) random selection of 1000 pairs from 8724 small-molecule drugs in DrugBank (Wishart ); (ii) 628 FDA-approved drug combinations curated by Cheng for hypertension and cancers (our case studies are on four types of cancers); (iii) random selection of 1000 pairs of FDA-approved drugs for the given disease, based on drug–disease dataset ‘SCMFDD-L’ (Zhang ).

4.3.3 Designed pairs follow network principles and improve toxicity

We first compare the compound combinations designed by our model and those from the baselines using the network score that reflects the network-based principle. Figure 2(a–d) shows that our designed combinations in all four cases, with higher network scores in distribution, respected the network principle more than the baselines (including the FDA-approved pairs not necessarily specific for the target disease). The observation is statistically significant with P-values ranging from 6E-74 to 7E-7 (one-sided Kolmogorov–Smirnov test; see more details in the Supplementary Tables S2 and S3). Such a result is thanks to the network-principled reward we introduced.

Fig. 2.

Comparison of network score and toxicity of RL-generated pairs of compounds (our proposed method) with three baselines, i.e. random pairs of DrugBank compounds, FDA-approved drug pairs and random pairs of FDA-approved drugs for four case-study diseases. Panels (a-h) are combinations of row measure (network score or toxicity) and column diseases (melanoma, lung cancer, ovarian cancer, and breast cancer)

We also examine whether drug combinations designed to follow the network principle could reduce toxicity from DDIs. DDIs are crucial when using drug combinations since they may trigger unexpected pharmacological effects, including adverse drug events. We used a deep-learning model DeepDDI (Ryu ) with a mean accuracy of 92.4% to predict for each combination the probabilities of 86 types of DDIs (we manually split them into 16 positives and 70 negatives; see details in the Supplementary Section S1.3). To summarize over the DDIs, we considered both maximum and mean probabilities of positive or negative ones. And we compared those distributions between our designed pairs and baselines in each disease.

Figure 2(e–h), using the mean probability among negative DDIs, shows that our compound pairs designed for all four diseases were predicted to have less chances of toxicity compared to the baselines. One-sided Kolmogorov–Smirnov tests attested to the statistical significance of the observation as P-values ranged between 2E-166 and 2E-53. More analyses can be found in the Supplementary Section S3.

Taken together, Figure 2 suggested that following the network principle in designing drug combinations would help reduce toxicity due to DDIs.

4.3.4 Designed pairs reproduce approved polypharmacology strategies

We next examine the DeepAffinity-predicted target genes of our designed pairs and compare them to the polypharmacology strategies outlined in Section 4.3.1 for each disease. Since improved network scores have been shown to correlate with lower toxicity, we used the scores to filter the 4000 combinations designed for each disease. Specifically, we retained combinations with network scores above 0.5 and below 0.1. These designs are shared along with the codes.

For melanoma, out of 69 combination designs retained, 26% were predicted to jointly cover BRAF and MEK genes in a complementary way. In other words, one molecule only targets BRAF and the other only targets MEK, according to our DeepAffinity (Karimi )-predicted IC50, echoing the systems pharmacology strategy of the drug combination of Dabrafenib and Trametinib. There were also other designs which demand further examination and potentially contain novel strategies. All retained designs were predicted to target the MAPK pathway to which BRAF and MEK belong.

For lung and ovarian cancers, the same filtering criteria retained 204 (896) compound combinations designed for ovarian (lung) cancer. As disease modules can be limited, MEK1/2 does not exist in the used modules for ovarian or lung cancer so a gene-level analysis cannot be performed as the melanoma case. Instead, we performed the pathway-level analysis and found that 50.9% (45.2%) of combination designs for ovarian (lung) cancer were predicted to jointly and complementarily cover the MAPK and PI3K signaling pathways, which echoes the combination of Buparlisib and Trametinib. Moreover, 99.5% (100%) of these retained designs were predicted to jointly target both pathways for ovarian (lung) cancer.

For breast cancer, 77 designed compound-combinations passed the filters. As CDK4/6 does not belong to the breast-cancer module due to the limitation of disease modules used, we again only performed a pathway-level analysis. Nine percent of the combinations were predicted to jointly and complementarily cover estrogen receptor-signaling and cell cycle pathways as Fulvestrant and Palbociclib do. Also, 74% of the retained combinations jointly cover these pathways. These two portions suggest that many designed combinations were predicted to simultaneously target both pathways (with possible overlapping genes). If we consider PI3K signaling rather than cell cycle pathway for CDK4/6, 15.5% of retained drug combinations were predicted to jointly and complementarily cover estrogen and PI3K signaling pathways and all of them did jointly.

4.3.5 Ablation study for RL-based drug-combination generation

Besides HVGAE for network and disease embedding, two of our novel contributions in RL-based drug set generations were network-principled reward and adversarial reward through GS-WGAN. To assess the effects of these contributions to our model, we performed ablation study for Stage 3 using the case of melanoma. We ablated the originally proposed model in two ways: removing the network-principled reward or replacing the GS-WGAN adversarial reward with the previously used GAN reward based on Jenson–Shannon divergence. Results in Figure 3 suggested that both rewards led to faster initial growth and higher saturation values in network-based scores.

Fig. 3.

Ablation study for RL: best network scores achieved by three variants of the proposed method over training iterations

5 Conclusion

In response to the need of accelerated and principled drug-combination design, we have recast the problem as graph-set generation in a chemically and net-biologically valid environment and developed the first deep generative model with novel adversarial award and drug-combination award in RL for the purpose. We have also designed hierarchical variation graph auto-encoders (HGVAE) to jointly embed domain knowledge such as gene–gene, disease–disease, gene–disease networks and learn disease representations to be conditioned on in the generative model for disease-specific drug combination. Our results indicate that HGVAE learns integrative gene and disease representations that are much more generalizable and informative than state-of-the-art graph unsupervised-learning methods. The results also indicate that the RL model learns to generate drug combinations following a network-based principle thanks to our adversarial and drug-combination rewards. Case studies involving four diseases indicate that drug combinations designed to follow network principles tend to have low toxicity from DDIs. These designs also encode systems pharmacology strategies echoing FDA-approved drug combinations as well as other potentially promising strategies. As the first generative model for disease-specific drug-combination design, our study allows for assessing and following network-based mechanistic hypotheses in efficiently searching the chemical combinatorial space and effectively designing drug combinations.

Click here for additional data file.