Toxicity Prediction with ML

Toxicity prediction is one of the most critical and challenging applications of ML in drug discovery. A compound must not only be potent against its target, it must also be safe to administer. Toxicity assays are expensive and time-consuming, and late-stage attrition due to unexpected toxicity remains a major problem for pharmaceutical pipelines. Computational toxicity prediction allows early triaging of large compound libraries to flag likely toxic compounds before in vitro or in vivo testing.

This guide covers the toxicity landscape in drug discovery, key benchmark datasets, featurization and model architectures, strategies for handling class imbalance, a worked code example with RDKit and scikit-learn, and the regulatory context for computational toxicity prediction.

Types of Toxicity in Drug Discovery

Toxicity is not a single property but a family of endpoints that arise from different mechanisms at different biological scales:

Each toxicity type requires different modeling approaches. For example, Ames mutagenicity is well-studied with classical rule-based and ML models achieving AUROC above 0.85. DILI (drug-induced liver injury) is much harder because it involves complex metabolic and idiosyncratic mechanisms that are not well captured by chemical structure alone.

Key Datasets

Tox21 (12 Assays)

The Tox21 challenge dataset contains around 8,000 compounds tested across 12 in vitro assays covering nuclear receptors (NR-AR, NR-AR-LBD, NR-AhR, NR-ER, NR-ER-LBD, NR-Aromatase, NR-PPAR-gamma) and stress response pathways (SR-ARE, SR-ATAD5, SR-HSE, SR-MMP, SR-p53). The Tox21 challenge in 2014 established deep multitask networks as competitive with Random Forests for this benchmark.

Key statistics:

• ~8,000 compounds across 12 assay tasks
• Positive rate per assay ranges from 3% to 19%
• Significant class imbalance in most tasks
• Official scaffold split provided for reproducible benchmarking

ClinTox

ClinTox contains approximately 1,478 drugs with two binary labels: (1) clinical trial toxicity failure and (2) FDA approval status. It is a challenging dataset because the positive class is small and the relationship between chemical structure and idiosyncratic clinical toxicity is indirect. ClinTox is used to evaluate whether models can distinguish between clinically failed and approved drugs based purely on structure.

SIDER (Side Effect Resource)

SIDER contains information about 1,427 marketed drugs and their side effects organized by organ system MedDRA terms. It is often used for multi-label toxicity classification, where the goal is to predict which side effect categories a new drug is likely to exhibit. Because the drugs in SIDER are all approved, the dataset represents a biased subset of chemical space.

Ames Mutagenicity

The Hansen mutagenicity dataset contains around 6,500 compounds with Ames test labels (mutagenic/non-mutagenic). The Ames test is a regulatory gold standard for genotoxicity screening and is required by ICH S2(R1) for all new chemical entities. QSAR models for Ames mutagenicity have a long history and are included in commercial tools like DEREK Nexus.

Featurization: Morgan Fingerprints and Physicochemical Descriptors

For Tox21-sized datasets, Morgan fingerprints (ECFP4, 2048 bits) are the standard featurization. They capture local structural features that are relevant to reactivity and binding, which underlies most toxicity mechanisms. Physicochemical descriptors add global context:

• MolLogP: lipophilicity affects membrane penetration and accumulation
• TPSA: polar surface area relates to membrane permeability and distribution
• NumHeavyAtoms: molecular size correlates with non-specific toxicity risks
• FractionCSP3: complexity index correlates with metabolic liability

Fragment-based alerts (structural alerts for specific toxicophores) can also be added as binary features. RDKit provides access to the FilterCatalog module with curated structural alert collections (PAINS, BMS, Dundee, GSK).

Models: Random Forest, XGBoost, Multitask Neural Networks

Random Forest on ECFP4

For individual Tox21 tasks, a Random Forest with class weight balancing on ECFP4 fingerprints achieves AUROC in the range of 0.72-0.82. This is the practical industry baseline for toxicity classification.

Multitask Neural Networks

The 2014 Tox21 challenge winner used a deep multitask neural network that shared lower layers across all 12 tasks. The multitask setting is especially beneficial for toxicity because:

• Many Tox21 assays measure related pathways (e.g., multiple nuclear receptor assays)
• Shared representations learn general reactivity patterns useful across tasks
• Effective data augmentation through tasks that have more positive examples

A typical multitask architecture: ECFP4 input (2048 bits) - Dense(2000, ReLU) - Dense(1000, ReLU) - Dense(500, ReLU) - 12 sigmoid output heads with task-specific losses.

Handling Class Imbalance in Toxicity Datasets

Toxicity datasets are highly imbalanced. In Tox21, most assays have 5-15% positive rate. This means a classifier predicting all-negative achieves 85-95% accuracy while being completely useless. Several strategies address this:

Always evaluate with AUROC and AUPRC (average precision). AUROC summarizes performance across all thresholds but can be misleading with severe imbalance. AUPRC is more informative for rare positive classes.

Code Example: Tox21 with RDKit and scikit-learn

import numpy as np
import pandas as pd
from rdkit import Chem
from rdkit.Chem import AllChem
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score, average_precision_score
from sklearn.model_selection import cross_val_score
import warnings
warnings.filterwarnings('ignore')

# -------------------------------------------------------
# 1. Load Tox21 dataset from DeepChem MoleculeNet
# -------------------------------------------------------
import deepchem as dc

tasks, datasets, transformers = dc.molnet.load_tox21(
    featurizer='ECFP',
    splitter='scaffold'
)

train_dataset, valid_dataset, test_dataset = datasets

# Extract SMILES and labels
# DeepChem ECFP featurizer returns pre-computed bit vectors
X_train = train_dataset.X  # shape (n_train, 1024)
y_train = train_dataset.y  # shape (n_train, 12)
w_train = train_dataset.w  # sample weights (mask for missing labels)

X_test  = test_dataset.X
y_test  = test_dataset.y
w_test  = test_dataset.w

print(f"Train: {X_train.shape}, Test: {X_test.shape}")
print(f"Tasks: {tasks}")

# -------------------------------------------------------
# 2. Train per-task Random Forest models
# -------------------------------------------------------
aurocs = []
auprcs = []

for task_idx, task_name in enumerate(tasks):
    # Mask out compounds with missing labels for this task
    train_mask = w_train[:, task_idx] > 0
    test_mask  = w_test[:, task_idx] > 0

    X_tr = X_train[train_mask]
    y_tr = y_train[train_mask, task_idx].astype(int)
    X_te = X_test[test_mask]
    y_te = y_test[test_mask, task_idx].astype(int)

    # Skip tasks with too few positives
    if y_tr.sum() < 10:
        continue

    pos_ratio = y_tr.mean()
    clf = RandomForestClassifier(
        n_estimators=200,
        max_depth=None,
        class_weight='balanced',
        n_jobs=-1,
        random_state=42
    )
    clf.fit(X_tr, y_tr)

    y_pred_proba = clf.predict_proba(X_te)[:, 1]
    auroc = roc_auc_score(y_te, y_pred_proba)
    auprc = average_precision_score(y_te, y_pred_proba)

    aurocs.append(auroc)
    auprcs.append(auprc)
    print(f"  {task_name:20s}: AUROC={auroc:.3f}, AUPRC={auprc:.3f}, "
          f"pos_rate={pos_ratio:.1%}")

print(f"\nMean AUROC: {np.mean(aurocs):.3f}")
print(f"Mean AUPRC: {np.mean(auprcs):.3f}")

Multitask Neural Network with PyTorch

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
import numpy as np

class MultitaskToxModel(nn.Module):
    """
    Multitask neural network for toxicity prediction.
    Shared layers learn a general molecular representation.
    Task-specific output heads predict each toxicity endpoint.
    """
    def __init__(self, input_dim: int, n_tasks: int, dropout: float = 0.25):
        super().__init__()
        self.shared = nn.Sequential(
            nn.Linear(input_dim, 2000),
            nn.ReLU(),
            nn.BatchNorm1d(2000),
            nn.Dropout(dropout),
            nn.Linear(2000, 1000),
            nn.ReLU(),
            nn.BatchNorm1d(1000),
            nn.Dropout(dropout),
            nn.Linear(1000, 500),
            nn.ReLU(),
            nn.BatchNorm1d(500),
            nn.Dropout(dropout),
        )
        # Separate output head per toxicity task
        self.heads = nn.ModuleList([nn.Linear(500, 1) for _ in range(n_tasks)])

    def forward(self, x):
        h = self.shared(x)
        return torch.cat([head(h) for head in self.heads], dim=1)


def train_multitask_tox(X_train, y_train, w_train, n_tasks, epochs=50, lr=1e-3):
    """Train multitask model with masked binary cross-entropy loss."""
    model = MultitaskToxModel(input_dim=X_train.shape[1], n_tasks=n_tasks)
    optimizer = optim.Adam(model.parameters(), lr=lr)

    X_t = torch.tensor(X_train, dtype=torch.float32)
    y_t = torch.tensor(y_train, dtype=torch.float32)
    w_t = torch.tensor(w_train, dtype=torch.float32)

    dataset = TensorDataset(X_t, y_t, w_t)
    loader  = DataLoader(dataset, batch_size=128, shuffle=True)

    for epoch in range(epochs):
        model.train()
        total_loss = 0.0
        for xb, yb, wb in loader:
            logits = model(xb)                        # (batch, n_tasks)
            bce    = nn.functional.binary_cross_entropy_with_logits(
                logits, yb, reduction='none'
            )
            # Mask missing labels (w=0 means no measurement)
            loss   = (bce * wb).sum() / (wb.sum() + 1e-8)
            optimizer.zero_grad()
            loss.backward()
            optimizer.step()
            total_loss += loss.item()

        if (epoch + 1) % 10 == 0:
            print(f"Epoch {epoch+1}/{epochs}: loss={total_loss/len(loader):.4f}")

    return model

# Usage (assuming X_train, y_train, w_train from DeepChem)
model = train_multitask_tox(
    X_train.astype(np.float32),
    y_train.astype(np.float32),
    w_train.astype(np.float32),
    n_tasks=len(tasks)
)

Evaluation: AUROC, AUPRC, and Calibration

For toxicity prediction, the evaluation must go beyond accuracy:

• AUROC: Measures ranking ability. The probability that a randomly chosen toxic compound is ranked higher than a randomly chosen non-toxic compound. Useful for comparing models regardless of threshold.
• AUPRC: Average precision summarizes the precision-recall curve. More sensitive to model quality on minority classes than AUROC. Use this as the primary metric when positive rate is below 10%.
• Calibration: Predicted probabilities should match observed toxic frequencies. Use sklearn.calibration.calibration_curve and Platt scaling or isotonic regression to calibrate. Well-calibrated probabilities enable compound prioritization by risk level rather than just binary classification.

When reporting results, always specify the splitting strategy (scaffold vs. random) and the missing-label handling policy, as different implementations yield very different numbers on the same Tox21 dataset.

Regulatory Context

Computational toxicity prediction is not just an academic exercise. Regulatory agencies explicitly recognize in silico toxicology tools:

• FDA (US): FDA Guidance on Genotoxicity Testing (ICH M7) accepts computational approaches for impurity qualification. The FDA has published guidance on the use of QSAR models for Ames mutagenicity prediction.
• EMA (EU): The European Medicines Agency follows ICH S2(R1) and M7, both of which allow QSAR tools as part of genotoxicity assessment.
• REACH (EU Chemicals): The Registration, Evaluation, Authorisation and Restriction of Chemicals regulation requires toxicity assessment for industrial chemicals and explicitly permits QSAR predictions as alternatives to animal testing.
• ECHA QSAR Toolbox: The European Chemicals Agency provides a free toolbox for regulatory QSAR-based toxicity prediction, based on structural alerts and analog-based prediction.

For regulatory submissions, the model must be validated according to OECD principles for QSAR model validation: defined endpoint, unambiguous algorithm, defined applicability domain, appropriate goodness-of-fit and predictivity, and mechanistic interpretation where possible.

Further Reading

• Mayr et al. (2016): DeepTox and Tox21 Challenge
• NIH Tox21 Program
• DeepChem MoleculeNet Tox21
• ECHA QSAR Toolbox
• ICH M7 Guideline (EMA)

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.