Selectivity Prediction for Kinase and GPCR Drug Discovery

Kinases make up ~530 members in the human kinome; GPCRs number over 800. A compound active against one family member will often show some activity against close homologs due to conserved ATP/orthosteric binding site geometry. Predicting selectivity computationally allows prioritising compounds that show both potency and specificity before expensive selectivity panel assays (KINOMEscan, PRISM, HTRF panels).

Proteochemometric (PCM) Modelling

PCM encodes both the compound (as ECFP4) and the target (as protein descriptors) as a combined feature vector. A single model learns to predict activity across all targets simultaneously, exploiting cross-target information.

import numpy as np
import pandas as pd
from rdkit import Chem
from rdkit.Chem import AllChem
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import GroupKFold, cross_val_score


def encode_target_sequence(sequence: str, window: int = 7) -> np.ndarray:
    """
    Encode a protein binding-site sub-sequence as a simple Z-scale descriptor vector.

    Uses 5 Z-scales per residue (Sandberg et al.): hydrophobicity, steric,
    polarity, hydrogen-bond, and electronic properties.

    Args:
        sequence: Single-letter amino acid sequence of binding site (e.g. kinase hinge)
        window: Number of residues to encode

    Returns:
        Flattened Z-scale vector of length window * 5.
    """
    # Z3 descriptors (simplified): [hydrophobicity, size, electronic]
    z_scales = {
        "A": [0.24, -2.32, 0.60], "C": [0.84, -1.67, 3.71], "D": [3.98, 0.93, 1.93],
        "E": [3.11, 0.26, -0.11], "F": [-4.22, 1.94, 1.06], "G": [2.05, -4.06, 0.36],
        "H": [2.47, 1.95, 0.26], "I": [-3.89, -1.73, -1.71], "K": [2.29, 0.89, -2.49],
        "L": [-4.28, -1.30, -1.49], "M": [-2.85, -0.22, 0.47], "N": [3.05, 1.62, 1.04],
        "P": [-1.66, 0.27, 1.84], "Q": [1.75, 0.50, -1.44], "R": [1.91, -0.59, -3.15],
        "S": [2.39, -1.07, 1.15], "T": [0.75, -2.18, -1.12], "V": [-2.59, -2.64, -1.54],
        "W": [-4.36, 3.94, 0.59], "Y": [-2.54, 2.44, 0.43],
    }
    default = [0.0, 0.0, 0.0]
    trimmed = (sequence[:window] if len(sequence) >= window else sequence.ljust(window, "G"))
    vec = []
    for aa in trimmed[:window]:
        vec.extend(z_scales.get(aa.upper(), default))
    return np.array(vec, dtype=float)


def build_pcm_features(
    smiles: str,
    target_sequence: str,
    fp_radius: int = 2,
    fp_nbits: int = 1024,
) -> np.ndarray:
    """
    Build a PCM feature vector by concatenating compound fingerprint and target descriptor.

    Args:
        smiles: Compound SMILES
        target_sequence: Target binding-site sequence (e.g. kinase hinge residues)
        fp_radius: ECFP radius
        fp_nbits: Fingerprint length

    Returns:
        Concatenated [compound_fp | target_zscales] feature vector.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return np.zeros(fp_nbits + 7 * 3)

    fp = list(AllChem.GetMorganFingerprintAsBitVect(mol, fp_radius, fp_nbits))
    target_feat = encode_target_sequence(target_sequence)
    return np.concatenate([fp, target_feat])


def train_pcm_model(
    smiles_list: list[str],
    target_sequences: list[str],
    activities: list[float],
    target_ids: list[str],
) -> dict:
    """
    Train a PCM model using compound + target features.

    Uses GroupKFold on target_ids: each target is a group, preventing
    the model from simply memorising per-target offsets.

    Args:
        smiles_list: Compound SMILES
        target_sequences: Binding-site sequences for each datapoint
        activities: pIC50/pKi values
        target_ids: Target identifiers (for group CV split)

    Returns:
        Dict with model, cross-target CV statistics.
    """
    X = np.array([
        build_pcm_features(smi, seq)
        for smi, seq in zip(smiles_list, target_sequences)
    ])
    y = np.array(activities)

    model = GradientBoostingRegressor(
        n_estimators=400, max_depth=5, learning_rate=0.05,
        subsample=0.8, random_state=42,
    )

    groups = np.array(target_ids)
    group_kfold = GroupKFold(n_splits=min(5, len(set(target_ids))))

    scores = cross_val_score(
        model, X, y,
        cv=group_kfold, groups=groups,
        scoring="neg_mean_absolute_error", n_jobs=-1,
    )

    model.fit(X, y)
    return {
        "pcm_model": model,
        "cross_target_cv_mae": round(-scores.mean(), 3),
        "cv_std": round(scores.std(), 3),
        "n_datapoints": len(y),
        "n_targets": len(set(target_ids)),
    }

Selectivity Entropy

import numpy as np


def selectivity_entropy(
    activities: dict[str, float],
    concentration: float = 10.0,
    kd_threshold: float = None,
) -> float:
    """
    Compute Selectivity Entropy (Sε) for a compound across a kinase panel.

    Sε = -Σ p_i * log(p_i) where p_i = Kd_i / Σ Kd_j

    Lower Sε = higher selectivity (activity concentrated on few targets).
    Sε = 0: perfectly selective (one target only).
    Sε = log(N): completely non-selective (equal affinity for all N targets).

    Reference: Uitdehaag & Zaman, PLOS ONE 2011.

    Args:
        activities: {target_name: Kd_value_nM} or {target_name: pKd}
        concentration: Screening concentration (nM) for binary conversion (unused here)
        kd_threshold: If set, only include targets with Kd <= threshold

    Returns:
        Selectivity entropy value.
    """
    kd_values = [v for v in activities.values() if v > 0]

    if kd_threshold is not None:
        kd_values = [v for v in kd_values if v <= kd_threshold]

    if not kd_values:
        return float("nan")

    kd_arr = np.array(kd_values, dtype=float)
    total = kd_arr.sum()
    if total == 0:
        return float("nan")

    p = kd_arr / total
    p = p[p > 0]  # Avoid log(0)
    entropy = float(-np.sum(p * np.log(p)))
    return round(entropy, 4)


def gini_selectivity(activities: dict[str, float]) -> float:
    """
    Alternative to entropy: Gini coefficient of activity distribution.
    Gini = 0 means perfectly equal (non-selective), 1 means perfectly selective.
    """
    values = np.sort(np.array(list(activities.values()), dtype=float))
    n = len(values)
    if n == 0 or values.sum() == 0:
        return 0.0
    cumsum = np.cumsum(values)
    gini = (n + 1 - 2 * np.sum(cumsum) / cumsum[-1]) / n
    return round(float(gini), 4)

KINOMEscan Data and Multi-Task Learning

Key Takeaways

• PCM models outperform single-target QSAR for selectivity prediction because they share information across related targets. Performance advantage is largest for targets with few datapoints that benefit from transfer from related targets.
• Selectivity entropy provides a single number summarising the breadth of activity across a panel. Pair it with the Gini coefficient: a compound can have low entropy while still hitting a few closely related targets.
• Multi-task neural networks (shared lower layers, target-specific upper layers) are state-of-the-art for kinase selectivity prediction. They consistently outperform independent single-target models at ~>50 targets.
• For GPCR functional selectivity (biased agonism), selectivity is not just potency but the ratio of G-protein vs. arrestin pathway activation. This requires separate functional assay endpoints beyond binding affinity alone.

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.