Multi-Target Drug Design

Polypharmacology — the deliberate design of compounds that engage multiple targets simultaneously — offers therapeutic advantages over single-target drugs for complex diseases such as cancer, neurodegeneration, and infection. A compound might need to inhibit kinase A while avoiding kinase B (selectivity), or intentionally inhibit both kinase A and protease C to overcome resistance mechanisms. ML is well-suited to learning these multi-target profiles from large ChEMBL bioactivity datasets.

This guide covers the full practical pipeline: curating clean multi-target data from ChEMBL, training single-task baselines, multi-task GNNs, and selectivity fingerprints, and prioritising compound sets using Pareto optimisation of the activity profile vector.

Data Curation from ChEMBL

import pandas as pd
import numpy as np
from chembl_webresource_client.new_client import new_client

def fetch_multitarget_data(
    target_chembl_ids: list[str],
    standard_types: list[str] = ["IC50", "Ki", "Kd"],
    pchembl_min: float = 5.0,  # pIC50 ≥ 5 means IC50 ≤ 10 μM
    n_per_target: int = 5000,
) -> pd.DataFrame:
    """
    Fetch bioactivity data for multiple targets from ChEMBL.

    Args:
        target_chembl_ids: List of ChEMBL target IDs (e.g. ['CHEMBL279', 'CHEMBL340'])
        standard_types: Measurement types to include
        pchembl_min: Minimum pCHEMBL value (excludes very weak binders)
        n_per_target: Max compounds per target

    Returns:
        Combined DataFrame with columns: smiles, target_id, pchembl_value, label
    """
    activity = new_client.activity
    records = []
    for target_id in target_chembl_ids:
        results = activity.filter(
            target_chembl_id=target_id,
            standard_type__in=standard_types,
            pchembl_value__gte=pchembl_min,
        ).only(
            "molecule_chembl_id", "canonical_smiles",
            "standard_type", "pchembl_value", "target_chembl_id"
        )[:n_per_target]

        for r in results:
            if r.get("canonical_smiles") and r.get("pchembl_value"):
                records.append({
                    "smiles": r["canonical_smiles"],
                    "target_id": r["target_chembl_id"],
                    "standard_type": r["standard_type"],
                    "pchembl_value": float(r["pchembl_value"]),
                    "label": 1 if float(r["pchembl_value"]) >= pchembl_min else 0,
                })
    return pd.DataFrame(records)


def build_multitarget_matrix(
    df: pd.DataFrame,
    targets: list[str],
    active_threshold: float = 6.5,
) -> pd.DataFrame:
    """
    Pivot the long-format data into a compound × target activity matrix.

    Values are pCHEMBL mediated by median; NaN for missing data.
    Binary activity labels computed at active_threshold.

    Returns:
        DataFrame indexed by SMILES with columns for each target's pchembl + binary label.
    """
    # Aggregate: take median pCHEMBL per (smiles, target) pair
    agg = df.groupby(["smiles", "target_id"])["pchembl_value"].median().reset_index()
    matrix = agg.pivot(index="smiles", columns="target_id", values="pchembl_value")
    matrix.columns.name = None
    # Add binary columns
    for t in targets:
        if t in matrix.columns:
            matrix[f"{t}_active"] = (matrix[t] >= active_threshold).astype(int)
    return matrix.reset_index()

Multi-Task GNN Training

A multi-task GNN shares a common molecular representation encoder across all targets, using separate output heads per target. Shared representations improve generalisation for data-sparse targets (a form of transfer learning between pharmacologically related targets).

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch_geometric.data import Data, Batch
from torch_geometric.nn import GCNConv, global_mean_pool
from typing import Optional

class MultiTaskGNN(nn.Module):
    """
    Multi-task Graph Neural Network for polypharmacology.

    Arguments:
        in_node_features: Atom feature dimensionality
        hidden_dim: Hidden GNN layer dimension
        n_layers: Number of GCN layers
        targets: List of target names (one output head per target)
        dropout: Dropout rate
    """
    def __init__(
        self,
        in_node_features: int = 74,
        hidden_dim: int = 256,
        n_layers: int = 3,
        targets: list[str] = ["target_1", "target_2"],
        dropout: float = 0.2,
    ):
        super().__init__()
        self.targets = targets
        self.dropout = dropout

        # Shared encoder: GCN layers
        self.convs = nn.ModuleList()
        self.bns = nn.ModuleList()
        in_dim = in_node_features
        for _ in range(n_layers):
            self.convs.append(GCNConv(in_dim, hidden_dim))
            self.bns.append(nn.BatchNorm1d(hidden_dim))
            in_dim = hidden_dim

        # Per-target output heads (regression: predict pCHEMBL)
        self.heads = nn.ModuleDict({
            t: nn.Sequential(
                nn.Linear(hidden_dim, 128),
                nn.ReLU(),
                nn.Dropout(dropout),
                nn.Linear(128, 1),
            )
            for t in targets
        })

    def encode(self, x, edge_index, batch):
        for conv, bn in zip(self.convs, self.bns):
            x = conv(x, edge_index)
            x = bn(x)
            x = F.relu(x)
            x = F.dropout(x, p=self.dropout, training=self.training)
        return global_mean_pool(x, batch)  # (batch_size, hidden_dim)

    def forward(
        self,
        x: torch.Tensor,
        edge_index: torch.Tensor,
        batch: torch.Tensor,
    ) -> dict[str, torch.Tensor]:
        h = self.encode(x, edge_index, batch)
        return {t: self.heads[t](h).squeeze(-1) for t in self.targets}


def multitask_loss(
    predictions: dict[str, torch.Tensor],
    labels: dict[str, torch.Tensor],
    masks: dict[str, torch.Tensor],
    task_weights: Optional[dict[str, float]] = None,
) -> torch.Tensor:
    """
    Compute masked multi-task MSE loss. Missing labels (NaN) are masked out.

    Args:
        predictions: {target: tensor of shape (batch,)}
        labels: {target: tensor of shape (batch,)} — may contain NaN
        masks: {target: bool tensor, True where label is available}
        task_weights: Optional per-task loss weights (default uniform)

    Returns:
        Scalar total loss
    """
    total_loss = 0.0
    n_tasks_with_data = 0
    for target in predictions:
        mask = masks[target]
        if mask.sum() == 0:
            continue
        pred = predictions[target][mask]
        true = labels[target][mask]
        task_loss = F.mse_loss(pred, true)
        weight = task_weights[target] if task_weights else 1.0
        total_loss = total_loss + weight * task_loss
        n_tasks_with_data += 1
    return total_loss / max(n_tasks_with_data, 1)

Selectivity Fingerprints

A selectivity fingerprint encodes the predicted activity profile of a molecule across a panel of targets. It can be used as a feature for clustering compounds by their polypharmacology profile, or as a constraint in generative design ("design compounds that look like this selectivity pattern").

import numpy as np
from sklearn.metrics.pairwise import cosine_similarity

def compute_selectivity_fingerprint(
    activity_vector: np.ndarray,
    target_names: list[str],
    active_threshold: float = 6.5,
) -> dict:
    """
    Convert a raw pCHEMBL activity vector into a selectivity fingerprint.

    Args:
        activity_vector: pCHEMBL predicted values for each target (NaN = unknown)
        target_names: Target names corresponding to each position
        active_threshold: pCHEMBL threshold for binary active/inactive

    Returns:
        Dict with binary fingerprint, selectivity ratio, and top active targets
    """
    binary = (activity_vector >= active_threshold).astype(int)
    active_mask = ~np.isnan(activity_vector) & binary.astype(bool)

    n_active = int(binary[~np.isnan(activity_vector)].sum())
    n_measured = int(~np.isnan(activity_vector)).sum()

    # Selectivity ratio: 1 means perfectly selective (1 of N targets),
    # n_measured means pan-active across all measured targets
    active_targets = [target_names[i] for i in np.where(active_mask)[0]]

    # Gini-style selectivity coefficient (Bento et al.)
    valid = activity_vector[~np.isnan(activity_vector)]
    if len(valid) > 1:
        sorted_vals = np.sort(valid) / valid.sum() if valid.sum() > 0 else valid
        n = len(sorted_vals)
        gini = (2 * np.sum((np.arange(1, n+1)) * sorted_vals) - (n + 1)) / n
    else:
        gini = float("nan")

    return {
        "binary_fp": binary.tolist(),
        "n_active": n_active,
        "n_measured": n_measured,
        "selectivity_ratio": round(n_active / n_measured, 3) if n_measured > 0 else None,
        "gini_selectivity": round(float(gini), 3) if not np.isnan(gini) else None,
        "active_targets": active_targets,
    }


def selectivity_profile_distance(
    fp1: np.ndarray,
    fp2: np.ndarray,
) -> float:
    """Cosine distance between two selectivity fingerprints."""
    sim = cosine_similarity(fp1.reshape(1, -1), fp2.reshape(1, -1))[0, 0]
    return 1.0 - float(sim)

Pareto Optimisation for Target Profile Selection

import numpy as np
from typing import Optional

def is_pareto_dominated(point: np.ndarray, candidates: np.ndarray) -> bool:
    """
    Check if a point is Pareto dominated by any candidate.
    Assumes maximisation of all objectives.

    A point is dominated if there exists another candidate that is
    better in at least one objective and at least as good in all others.
    """
    for candidate in candidates:
        if np.all(candidate >= point) and np.any(candidate > point):
            return True
    return False


def pareto_front(
    objective_matrix: np.ndarray,
    target_names: Optional[list[str]] = None,
) -> np.ndarray:
    """
    Compute the Pareto front (non-dominated set) for multi-objective optimisation.

    Args:
        objective_matrix: Shape (n_compounds, n_objectives).
                          Each column is an objective to maximise.
        target_names: Optional names for objectives

    Returns:
        Boolean mask of shape (n_compounds,) marking Pareto-optimal compounds.
    """
    n = len(objective_matrix)
    is_pareto = np.ones(n, dtype=bool)
    for i in range(n):
        if not is_pareto[i]:
            continue
        others = objective_matrix[is_pareto]
        if is_pareto_dominated(objective_matrix[i], others):
            is_pareto[i] = False
    return is_pareto


# Example: select compounds that are active on targets A and B
# but NOT active on off-target C (selectivity constraint)
def design_selectivity_filter(
    predicted_pchembl: np.ndarray,    # shape (n_compounds, n_targets)
    target_names: list[str],
    desired_targets: list[str],
    avoid_targets: list[str],
    desired_threshold: float = 6.5,   # pCHEMBL >= 6.5 = active
    avoid_threshold: float = 6.0,     # pCHEMBL < 6.0 = selective/inactive
) -> np.ndarray:
    """
    Filter compounds by desired multi-target activity profile.

    Args:
        predicted_pchembl: (n_compounds, n_targets) pCHEMBL predictions
        target_names: Names of all targets
        desired_targets: Must be active on these targets
        avoid_targets: Must be inactive on these targets

    Returns:
        Boolean mask of passing compounds.
    """
    n_cmpds = predicted_pchembl.shape[0]
    mask = np.ones(n_cmpds, dtype=bool)

    for target in desired_targets:
        col = target_names.index(target)
        mask &= (predicted_pchembl[:, col] >= desired_threshold)

    for target in avoid_targets:
        col = target_names.index(target)
        mask &= (predicted_pchembl[:, col] < avoid_threshold)

    return mask

Key Takeaways

• Multi-task GNNs consistently outperform single-task models for sparse targets: when a target has <500 active compounds in ChEMBL, sharing a representation with related targets improves AUC by 3–8%.
• Selectivity fingerprints (binary active/inactive over a target panel) enable clustering compounds by polypharmacology profile and seeding generative models conditioned on desired target patterns.
• The Gini selectivity coefficient is more informative than simple selectivity ratio — it captures whether a compound is partially active on many targets vs. strongly active on one.
• Pareto optimisation is essential for multi-target design: a compound optimally balancing 3–5 activity constraints can only be found by mapping the Pareto front, not by rank-ordering a single score.
• Always curate data carefully: ChEMBL contains measurements from different assay formats and labs. Use pCHEMBL ≥ 5 filter and median-aggregate duplicates to avoid noisy pseudo-activity cliffs.

Try it live

Explore a related model hands-on: open Molecule Screening Engine — Similarity Search Live and change its parameters while it runs. Everything stays in your browser.