PK/PD Modelling with ML

Approximately 40% of clinical failures are attributed to poor PK/PD properties — despite preclinical screening. ML PK prediction models can flag compounds for in vitro PK follow-up earlier in lead optimisation, but they are surrogates, not replacements, for measured data. The goal is triage: use predicted PK to guide which compounds are worth measuring, and which are chemically doomed.

Key PK Parameters and Their ML Predictability

Predicting Hepatic Clearance

import numpy as np
import pandas as pd
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline


PHYSICOCHEMICAL_PROPS = [
    "MolWt", "MolLogP", "TPSA", "NumHAcceptors", "NumHDonors",
    "NumRotatableBonds", "RingCount", "FractionCSP3",
    "NumAromaticRings", "HeavyAtomCount",
]


def compute_descriptors(smiles: str) -> dict | None:
    """
    Compute ECFP4 + physicochemical descriptors for PK prediction.

    Returns:
        Dict with fingerprint bits and RDKit descriptor values, or None if invalid.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    fp = list(AllChem.GetMorganFingerprintAsBitVect(mol, 2, 1024))
    desc = {}
    for prop in PHYSICOCHEMICAL_PROPS:
        try:
            desc[prop] = getattr(Descriptors, prop)(mol)
        except AttributeError:
            desc[prop] = 0.0

    return {"fp": fp, "descriptors": list(desc.values())}


def build_clearance_model(
    smiles_list: list[str],
    clearance_log: list[float],
) -> dict:
    """
    Train a hepatic clearance prediction model on intrinsic CLint data.

    Args:
        smiles_list: Canonical SMILES
        clearance_log: log10(CLint in µL/min/mg protein) values

    Returns:
        Dict with trained model and CV statistics.
    """
    features, targets = [], []
    for smi, cl in zip(smiles_list, clearance_log):
        desc = compute_descriptors(smi)
        if desc is not None:
            features.append(desc["fp"] + desc["descriptors"])
            targets.append(cl)

    X = np.array(features)
    y = np.array(targets)

    model = Pipeline([
        ("scaler", StandardScaler()),
        ("gbr", GradientBoostingRegressor(
            n_estimators=500, max_depth=5, learning_rate=0.05,
            subsample=0.8, min_samples_leaf=5, random_state=42,
        )),
    ])

    cv_scores = cross_val_score(model, X, y, cv=5, scoring="neg_mean_absolute_error", n_jobs=-1)
    model.fit(X, y)

    return {
        "model": model,
        "cv_mae_log_clearance": round(-cv_scores.mean(), 3),
        "n_compounds": len(y),
        "description": "log10(CLint µL/min/mg); > 1.5 = high clearance risk",
    }

Two-Compartment PK Simulation

from scipy.integrate import odeint
import numpy as np


def two_compartment_pk(
    dose_mg_kg: float,
    cl_l_h_kg: float,
    vd_central_l_kg: float,
    vd_peripheral_l_kg: float,
    k12: float,
    k21: float,
    body_weight_kg: float = 70.0,
    t_max_h: float = 24.0,
    n_timepoints: int = 200,
    route: str = "iv_bolus",
    ka: float = 1.0,
    f_oral: float = 1.0,
) -> dict:
    """
    Simulate two-compartment PK model for a drug candidate.

    Compartments: central (plasma) + peripheral (tissues)

    Args:
        dose_mg_kg: Dose in mg/kg
        cl_l_h_kg: Clearance in L/h/kg
        vd_central_l_kg: Central volume of distribution L/kg
        vd_peripheral_l_kg: Peripheral volume of distribution L/kg
        k12: Transfer rate constant central → peripheral (1/h)
        k21: Transfer rate constant peripheral → central (1/h)
        f_oral: Oral bioavailability (0–1)
        ka: Absorption rate constant (1/h) for oral route
        route: 'iv_bolus' or 'oral'

    Returns:
        Dict with time array and plasma concentration array.
    """
    BW = body_weight_kg
    CL = cl_l_h_kg * BW        # L/h total
    V1 = vd_central_l_kg * BW  # L central
    V2 = vd_peripheral_l_kg * BW  # L peripheral
    ke = CL / V1                # Elimination rate

    t = np.linspace(0, t_max_h, n_timepoints)
    dose_mg = dose_mg_kg * BW

    if route == "iv_bolus":
        def odes(y, t):
            c1, c2 = y
            dC1 = -(ke + k12) * c1 + k21 * c2 * V2 / V1
            dC2 = k12 * c1 * V1 / V2 - k21 * c2
            return [dC1, dC2]

        c0 = dose_mg / V1   # mg/L
        y0 = [c0, 0.0]

    else:  # oral
        def odes(y, t):
            a, c1, c2 = y
            absorbed = ka * a
            dA = -ka * a
            dC1 = (f_oral * absorbed / V1) - (ke + k12) * c1 + k21 * c2 * V2 / V1
            dC2 = k12 * c1 * V1 / V2 - k21 * c2
            return [dA, dC1, dC2]

        y0 = [dose_mg, 0.0, 0.0]

    sol = odeint(odes, y0, t)
    plasma_conc = sol[:, 0] if route == "iv_bolus" else sol[:, 1]  # mg/L

    auc = float(np.trapz(plasma_conc, t))
    cmax = float(plasma_conc.max())
    tmax_idx = int(plasma_conc.argmax())
    tmax = float(t[tmax_idx])

    return {
        "time_h": t.tolist(),
        "plasma_conc_mg_l": plasma_conc.tolist(),
        "auc_mg_h_l": round(auc, 3),
        "cmax_mg_l": round(cmax, 4),
        "tmax_h": round(tmax, 2),
        "t_half_h": round(0.693 / ke, 2),
    }

Key Takeaways

• ML PK models are most useful for rank-ordering series members (which analogue has lower clearance) rather than predicting absolute values. Absolute MAE of 0.3–0.5 log units corresponds to ~2–3× error in clearance — sufficient for triage decisions.
• Always correct clearance for protein binding before PK/PD integration: the pharmacologically active species is free (unbound) drug. A highly protein-bound compound with good total clearance may have very low free drug exposure.
• The two-compartment model is sufficient for most small-molecule lead optimisation PK analysis. More complex PBPK (physiologically-based PK) models are needed for tissue distribution predictions and drug-drug interaction assessment.
• Predicted half-life (t½ = 0.693 × Vd / CL) is a derived parameter — errors in both Vd and CL predictions amplify. Always verify predicted half-life with measured data before using for dose regimen design.

Try it live

Explore a related model hands-on: open Population PK/PD Bayesian Dosing and change its parameters while it runs. Everything stays in your browser.