Chemical Space Exploration

Drug-like chemical space is often cited as containing 10²³–10⁶⁰ possible stable organic molecules with MW < 500 Da. Even the largest commercial screening libraries contain only ~10⁸ compounds, leaving the vast majority of space unexplored. Intelligently exploring and sampling this space — to bias towards drug-like, synthesisable, and diverse compounds — is a core skill in computational drug discovery.

This guide covers the practical toolkit for chemical space analysis: physicochemical filters to define "drug-like" regions, diversity-based subset selection, dimensionality reduction for visualisation, and metrics to quantify coverage relative to known bioactive compounds.

Drug-Likeness Filters

Before exploring, it is useful to filter molecular databases to drug-like (or lead-like) regions. Several rule sets define these boundaries:

from rdkit import Chem
from rdkit.Chem import Descriptors, Lipinski, rdMolDescriptors, FilterCatalog

def drug_likeness_profile(smiles: str) -> dict:
    """
    Compute drug-likeness properties and filter results for a molecule.
    Returns descriptors, filter pass/fail, and PAINS classification.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return {"error": "Invalid SMILES", "smiles": smiles}

    mw = Descriptors.ExactMolWt(mol)
    clogp = Descriptors.MolLogP(mol)
    hbd = Lipinski.NumHDonors(mol)
    hba = Lipinski.NumHAcceptors(mol)
    rot = rdMolDescriptors.CalcNumRotatableBonds(mol)
    tpsa = rdMolDescriptors.CalcTPSA(mol)
    n_rings = rdMolDescriptors.CalcNumRings(mol)
    n_stereo = len(Chem.FindMolChiralCenters(mol, includeUnassigned=True))

    # Lipinski Ro5
    ro5_viols = sum([mw > 500, clogp > 5, hbd > 5, hba > 10])

    # Veber
    veber_pass = rot <= 10 and tpsa <= 140

    # Lead-like (Ro4)
    leadlike = mw <= 400 and clogp <= 4 and hbd <= 4 and hba <= 8

    # PAINS check using RDKit FilterCatalog
    params = FilterCatalog.FilterCatalogParams()
    params.AddCatalog(FilterCatalog.FilterCatalogParams.FilterCatalogs.PAINS)
    catalog = FilterCatalog.FilterCatalog(params)
    pains_match = catalog.GetFirstMatch(mol)
    is_pains = pains_match is not None
    pains_alert = pains_match.GetDescription() if is_pains else None

    return {
        "smiles": smiles,
        "mw": round(mw, 2),
        "clogp": round(clogp, 2),
        "hbd": hbd,
        "hba": hba,
        "rotatable_bonds": rot,
        "tpsa": round(tpsa, 2),
        "n_rings": n_rings,
        "n_stereocenters": n_stereo,
        "lipinski_violations": ro5_viols,
        "lipinski_pass": ro5_viols <= 1,
        "veber_pass": veber_pass,
        "lead_like": leadlike,
        "is_pains": is_pains,
        "pains_alert": pains_alert,
    }


def filter_library(
    smiles_list: list[str],
    max_lipinski_violations: int = 1,
    require_veber: bool = True,
    remove_pains: bool = True,
) -> list[str]:
    """
    Filter a SMILES library to drug-like compounds.
    Returns filtered list of SMILES.
    """
    results = []
    for smi in smiles_list:
        profile = drug_likeness_profile(smi)
        if "error" in profile:
            continue
        if profile["lipinski_violations"] > max_lipinski_violations:
            continue
        if require_veber and not profile["veber_pass"]:
            continue
        if remove_pains and profile["is_pains"]:
            continue
        results.append(smi)
    return results

Diversity Selection: MaxMin Algorithm

When you have a filtered library of N compounds and want to select a representative diverse subset of k molecules for screening, MaxMin (also called MaxMin dissimilarity selection) is the standard algorithm. It picks each new compound to maximise its minimum Tanimoto distance to the already-selected set.

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs
import numpy as np
import random

def maxmin_diversity_selection(
    smiles_list: list[str],
    n_select: int,
    fp_radius: int = 2,
    fp_bits: int = 2048,
    seed: int = 42,
) -> list[int]:
    """
    MaxMin diversity selection algorithm (Shenkin & McDonald, 1994).

    Selects indices of n_select compounds that maximally cover chemical space.

    Args:
        smiles_list: Input SMILES list
        n_select: Number of compounds to select
        fp_radius: Morgan radius for ECFP fingerprint
        fp_bits: Fingerprint bit length
        seed: Random seed for initial compound selection

    Returns:
        List of selected indices into smiles_list
    """
    # Generate fingerprints
    fps = []
    valid_idx = []
    for i, smi in enumerate(smiles_list):
        mol = Chem.MolFromSmiles(smi)
        if mol:
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, fp_radius, nBits=fp_bits)
            fps.append(fp)
            valid_idx.append(i)

    if len(fps) < n_select:
        return valid_idx

    # Pick first compound randomly
    random.seed(seed)
    selected = [random.randint(0, len(fps) - 1)]

    # Precompute distances from each fp to all selected
    # Maintain min_distances array for efficiency
    min_distances = np.ones(len(fps))

    def update_min_distances(new_selected_idx: int) -> None:
        new_fp = fps[new_selected_idx]
        for j, fp in enumerate(fps):
            sim = DataStructs.TanimotoSimilarity(new_fp, fp)
            dist = 1.0 - sim
            if dist < min_distances[j]:
                min_distances[j] = dist

    # Initialise with first selection
    update_min_distances(selected[0])
    # Cannot re-select already chosen
    min_distances[selected[0]] = -1.0

    while len(selected) < n_select:
        # Pick compound with max min-distance to selected set
        next_idx = int(np.argmax(min_distances))
        selected.append(next_idx)
        update_min_distances(next_idx)
        min_distances[next_idx] = -1.0  # Mark as selected

    return [valid_idx[i] for i in selected]

Butina Clustering for Diverse Subset Selection

Butina clustering is a fast, distance-based clustering algorithm that does not require specifying the number of clusters upfront. It uses a Tanimoto similarity threshold to assign molecules to clusters, and selecting one compound per cluster gives a diverse representative subset.

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs
from rdkit.ML.Cluster import Butina

def butina_cluster(
    smiles_list: list[str],
    cutoff: float = 0.4,
    fp_radius: int = 2,
    fp_bits: int = 2048,
) -> dict:
    """
    Cluster molecules using the Butina (Taylor-Butina) algorithm.

    Args:
        smiles_list: Input SMILES list
        cutoff: Distance cutoff (1 - Tanimoto). 0.4 = ~60% similarity threshold
        fp_radius: Morgan fingerprint radius
        fp_bits: Fingerprint bits

    Returns:
        Dict with cluster assignments and cluster centroids (representative SMILES)
    """
    fps = []
    valid_smiles = []
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol:
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, fp_radius, nBits=fp_bits)
            fps.append(fp)
            valid_smiles.append(smi)

    # Compute lower-triangle distance matrix
    n = len(fps)
    dists = []
    for i in range(1, n):
        sims = DataStructs.BulkTanimotoSimilarity(fps[i], fps[:i])
        dists.extend([1.0 - s for s in sims])

    # Run Butina clustering
    clusters = Butina.ClusterData(dists, n, cutoff, isDistData=True)

    # clusters is a tuple of tuples (centroid first in each)
    cluster_assignments = {smi: -1 for smi in valid_smiles}
    centroids = []
    for cluster_id, cluster in enumerate(clusters):
        centroid_idx = cluster[0]
        centroids.append(valid_smiles[centroid_idx])
        for mol_idx in cluster:
            cluster_assignments[valid_smiles[mol_idx]] = cluster_id

    return {
        "n_clusters": len(clusters),
        "cluster_sizes": sorted([len(c) for c in clusters], reverse=True),
        "centroids": centroids,
        "assignments": cluster_assignments,
    }

PCA and UMAP Visualisation of Chemical Space

import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import pandas as pd

def compute_physicochemical_descriptors(smiles_list: list[str]) -> pd.DataFrame:
    """
    Compute a standard set of 9 physicochemical descriptors for visualisation.
    """
    records = []
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            continue
        from rdkit.Chem import Lipinski, rdMolDescriptors
        records.append({
            "smiles": smi,
            "mw":     Descriptors.ExactMolWt(mol),
            "clogp":  Descriptors.MolLogP(mol),
            "hbd":    Lipinski.NumHDonors(mol),
            "hba":    Lipinski.NumHAcceptors(mol),
            "tpsa":   rdMolDescriptors.CalcTPSA(mol),
            "rot":    rdMolDescriptors.CalcNumRotatableBonds(mol),
            "rings":  rdMolDescriptors.CalcNumRings(mol),
            "arom":   rdMolDescriptors.CalcNumAromaticRings(mol),
            "csp3":   rdMolDescriptors.CalcFractionCSP3(mol),
        })
    return pd.DataFrame(records)


def pca_chemical_space(
    smiles_list: list[str],
    n_components: int = 2,
) -> pd.DataFrame:
    """
    Project molecules into 2D PCA space using physicochemical descriptors.

    Returns DataFrame with smiles, PC1, PC2, and explained variance info.
    """
    desc_df = compute_physicochemical_descriptors(smiles_list)
    feature_cols = [c for c in desc_df.columns if c != "smiles"]
    X = desc_df[feature_cols].values

    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    pca = PCA(n_components=n_components)
    coords = pca.fit_transform(X_scaled)

    result = desc_df[["smiles"]].copy()
    for i in range(n_components):
        result[f"PC{i+1}"] = coords[:, i]
    result["pc1_var"] = pca.explained_variance_ratio_[0]
    result["pc2_var"] = pca.explained_variance_ratio_[1]
    return result


def umap_chemical_space(
    smiles_list: list[str],
    fp_bits: int = 2048,
    n_neighbors: int = 15,
    min_dist: float = 0.1,
    metric: str = "jaccard",   # Jaccard = 1 - Tanimoto for binary fps
) -> pd.DataFrame:
    """
    Project molecules into 2D UMAP space using Morgan fingerprints.

    Args:
        smiles_list: Input SMILES
        fp_bits: Morgan fingerprint size
        n_neighbors: UMAP n_neighbors (controls local vs global structure)
        min_dist: UMAP min_dist (controls point packing)
        metric: Distance metric — 'jaccard' is equivalent to Tanimoto distance

    Returns:
        DataFrame with smiles, umap1, umap2
    """
    import umap  # pip install umap-learn
    fps = []
    valid_smiles = []
    for smi in smiles_list:
        mol = Chem.MolFromSmiles(smi)
        if mol:
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=fp_bits)
            arr = np.zeros(fp_bits, dtype=np.uint8)
            from rdkit.DataStructs import ConvertToNumpyArray
            ConvertToNumpyArray(fp, arr)
            fps.append(arr)
            valid_smiles.append(smi)

    X = np.vstack(fps)
    reducer = umap.UMAP(
        n_components=2,
        n_neighbors=n_neighbors,
        min_dist=min_dist,
        metric=metric,
        random_state=42,
    )
    coords = reducer.fit_transform(X)
    return pd.DataFrame({
        "smiles": valid_smiles,
        "umap1": coords[:, 0],
        "umap2": coords[:, 1],
    })

Chemical Space Coverage Metrics

Once you have a 2D or high-dimensional projected library, you can quantify how well it covers a reference space (e.g., approved drugs or known bioactives from ChEMBL).

import numpy as np
from scipy.spatial import ConvexHull
from sklearn.neighbors import NearestNeighbors

def nearest_neighbour_coverage(
    library_fps: np.ndarray,
    reference_fps: np.ndarray,
    tanimoto_threshold: float = 0.4,
) -> dict:
    """
    Compute what fraction of reference compounds are 'covered' by the library
    (i.e., have at least one library neighbour within the Tanimoto threshold).

    Uses Jaccard (= Tanimoto for binary) via sklearn NearestNeighbors.

    Args:
        library_fps: Binary fingerprint matrix (n_library x fp_bits)
        reference_fps: Binary fingerprint matrix (n_reference x fp_bits)
        tanimoto_threshold: Coverage threshold (default 0.4 = 60% similarity)

    Returns:
        Dict with coverage_fraction and n_covered
    """
    distance_threshold = 1.0 - tanimoto_threshold
    nbrs = NearestNeighbors(
        n_neighbors=1, metric="jaccard", algorithm="ball_tree"
    ).fit(library_fps)
    distances, _ = nbrs.kneighbors(reference_fps)
    min_distances = distances[:, 0]
    covered = (min_distances <= distance_threshold).sum()
    return {
        "n_reference": len(reference_fps),
        "n_covered": int(covered),
        "coverage_fraction": round(covered / len(reference_fps), 4),
        "mean_min_distance": round(float(min_distances.mean()), 4),
    }


def internal_diversity(fps_matrix: np.ndarray, sample_size: int = 1000) -> float:
    """
    Estimate mean pairwise Tanimoto distance within a library (internal diversity).
    Samples random pairs for efficiency on large libraries.

    Returns mean Tanimoto distance (0 = identical, 1 = fully diverse).
    """
    n = len(fps_matrix)
    rng = np.random.default_rng(42)
    pairs = rng.integers(0, n, size=(sample_size, 2))
    distances = []
    for i, j in pairs:
        if i == j:
            continue
        a, b = fps_matrix[i].astype(float), fps_matrix[j].astype(float)
        intersection = np.sum(a * b)
        union = np.sum(np.clip(a + b, 0, 1))
        jaccard_sim = intersection / union if union > 0 else 0.0
        distances.append(1.0 - jaccard_sim)
    return float(np.mean(distances)) if distances else 0.0

Comparing Library Coverage Against Known Drugs

Key Takeaways

• Drug-like space is defined by multiple overlapping rule sets (Ro5, Veber, REOS); using Ro5 + PAINS removal is a standard starting point for any HTS or virtual screen.
• MaxMin diversity selection efficiently selects a maximally diverse subset in O(N × k) time and is suitable for selecting <10,000 compounds from million-scale libraries.
• Butina clustering at Tanimoto threshold 0.6 (cutoff 0.4) is the standard method for reducing redundancy in screening collections before modelling.
• UMAP with Jaccard distance on Morgan fingerprints gives more chemically meaningful embeddings than PCA on physicochemical descriptors, but PCA is faster and more interpretable.
• Nearest-neighbour coverage at Tanimoto ≥ 0.4 is the standard metric for quantifying how well a library represents a reference set like DrugBank or ChEMBL.

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.