Predicting Drug-Drug Interactions with Graph Neural Networks

Introduction

Drug-Drug Interactions (DDIs) can cause severe adverse effects when multiple medications are prescribed simultaneously. Predicting these interactions computationally is a massive challenge because drugs are complex molecular graphs, and their interactions form an even larger biomedical knowledge graph. In this tutorial, we will use PyTorch Geometric (PyG) to build a Graph Neural Network (GNN) that predicts potential DDIs by analyzing molecular structures and known interaction networks.

Architecture Overview

graph TD
    SMILES_A(Drug A SMILES) --> GraphA[Molecular Graph A]
    SMILES_B(Drug B SMILES) --> GraphB[Molecular Graph B]
    
    GraphA --> GCN(Graph Convolutional Network - Shared Weights)
    GraphB --> GCN
    
    GCN --> EmbA(Graph Embedding A - 64 dim)
    GCN --> EmbB(Graph Embedding B - 64 dim)
    
    EmbA --> Concat[Concatenate Embeddings]
    EmbB --> Concat
    
    Concat --> MLP(Multi-Layer Perceptron)
    MLP --> Pred(Probability of Interaction: 0.0 to 1.0)
        

Prerequisites

• Python 3.9+
• PyTorch, PyTorch Geometric
• RDKit (for molecular graph extraction)
• NetworkX (for graph visualization)

Step 1: Install Dependencies

pip install torch torchvision torchaudio
pip install torch_geometric
pip install rdkit networkx matplotlib

Step 2: Converting SMILES to Graphs

Drugs are typically represented as SMILES strings. We must parse these strings into mathematical graphs where nodes are atoms (with features like atomic number and hybridization) and edges are chemical bonds (with features like bond type).

from rdkit import Chem
import torch
from torch_geometric.data import Data

def smiles_to_graph(smiles):
    mol = Chem.MolFromSmiles(smiles)
    if not mol: return None
    
    # Node features: Atomic number
    node_features = []
    for atom in mol.GetAtoms():
        node_features.append([atom.GetAtomicNum()])
    x = torch.tensor(node_features, dtype=torch.float)
    
    # Edge indices: Bonds
    edges = []
    for bond in mol.GetBonds():
        i = bond.GetBeginAtomIdx()
        j = bond.GetEndAtomIdx()
        edges += [[i, j], [j, i]] # Undirected graph
    
    edge_index = torch.tensor(edges, dtype=torch.long).t().contiguous()
    
    return Data(x=x, edge_index=edge_index)

# Example: Aspirin
aspirin_graph = smiles_to_graph('CC(=O)OC1=CC=CC=C1C(=O)O')
print(f"Aspirin Graph: {aspirin_graph}")

Step 3: Building a Graph Convolutional Network (GCN)

We need an encoder that can learn a latent vector representation (embedding) for a given molecular graph. A Graph Convolutional Network (GCN) aggregates information from an atom's neighbors to build this representation.

import torch.nn.functional as F
from torch_geometric.nn import GCNConv, global_mean_pool

class MolecularGCN(torch.nn.Module):
    def __init__(self, num_node_features, hidden_channels):
        super(MolecularGCN, self).__init__()
        self.conv1 = GCNConv(num_node_features, hidden_channels)
        self.conv2 = GCNConv(hidden_channels, hidden_channels)

    def forward(self, x, edge_index, batch):
        # 1. Obtain node embeddings
        x = self.conv1(x, edge_index)
        x = x.relu()
        x = self.conv2(x, edge_index)
        x = x.relu()

        # 2. Readout layer: Average all node embeddings to get a single graph embedding
        x = global_mean_pool(x, batch)  # [batch_size, hidden_channels]
        
        return x

# Initialize encoder
num_features = 1 # We only used atomic number
hidden_dim = 64
encoder = MolecularGCN(num_node_features=num_features, hidden_channels=hidden_dim)

Step 4: Predicting Interactions and Training with Negative Sampling

To predict if Drug A interacts with Drug B, we pass both drugs through our GCN to get their embeddings. We then concatenate these embeddings and pass them through a Multi-Layer Perceptron (MLP) to output a probability. Because databases like DrugBank only list known interactions (positive samples), we must synthetically generate non-interacting pairs (negative sampling) to train the model properly using Binary Cross Entropy.

class DDIPredictor(torch.nn.Module):
    def __init__(self, encoder, hidden_dim):
        super(DDIPredictor, self).__init__()
        self.encoder = encoder
        self.mlp = torch.nn.Sequential(
            torch.nn.Linear(hidden_dim * 2, 128),
            torch.nn.ReLU(),
            torch.nn.Linear(128, 1),
            torch.nn.Sigmoid()
        )
        
    def forward(self, data_a, data_b):
        # Get embeddings for both drugs
        emb_a = self.encoder(data_a.x, data_a.edge_index, data_a.batch)
        emb_b = self.encoder(data_b.x, data_b.edge_index, data_b.batch)
        
        # Concatenate and predict
        combined = torch.cat([emb_a, emb_b], dim=1)
        return self.mlp(combined)

predictor = DDIPredictor(encoder, hidden_dim)
optimizer = torch.optim.Adam(predictor.parameters(), lr=0.001)
criterion = torch.nn.BCELoss()

# --- Training Loop Pseudocode with Negative Sampling ---
# for epoch in range(epochs):
#     predictor.train()
#     total_loss = 0
#     
#     for batch_a, batch_b, is_interacting in dataloader:
#         optimizer.zero_grad()
#         
#         # is_interacting is 1.0 for known pairs, 0.0 for randomly sampled negative pairs
#         predictions = predictor(batch_a, batch_b).squeeze()
#         loss = criterion(predictions, is_interacting)
#         
#         loss.backward()
#         optimizer.step()
#         total_loss += loss.item()
#     print(f"Epoch {epoch} Loss: {total_loss}")

Conclusion

Predicting DDIs is natively a graph problem. By utilizing PyTorch Geometric, we can map chemical structures directly into a mathematical space where interactions can be inferred. This approach forms the foundation of modern computational pharmacology, adverse event prediction systems, and computational drug repurposing frameworks.