Aqueous Solubility Prediction

Aqueous solubility — the ability of a compound to dissolve in water — is one of the most important physicochemical properties in drug discovery. A compound that precipitates in aqueous buffer will give unreliable assay results, and a drug that is insoluble in the GI tract will have poor oral bioavailability regardless of its potency. Predicting solubility in silico from molecular structure allows projects to filter out insoluble compounds before synthesis and to guide structural modification toward more soluble chemical space.

The standard measure is logS: the base-10 logarithm of molar aqueous solubility at pH 7.4 and 25 °C. A logS of −1 corresponds to 0.1 mol/L (highly soluble), while −6 corresponds to 1 μmol/L (practically insoluble for oral drugs). Most lead-like drugs fall between logS −1 and −5.

Key Datasets

Solubility prediction models are only as good as their training data. Several public datasets are commonly used as benchmarks:

The Delaney dataset is small but clean and universally used for benchmarking. For production models, AqSolDB or a curated ChEMBL extract is preferred. Always check for measurement protocol heterogeneity — shake-flask vs. nephelometry vs. kinetic solubility can differ by 1–2 log units for the same compound.

The ESOL Model: A Physicochemical Baseline

The ESOL (Estimated Solubility) model from Delaney (2004) is the classic linear regression baseline. It predicts logS using six interpretable descriptors:

logS = 0.16 − 0.63 × clogP − 0.0062 × MW + 0.066 × Arom − 0.74 × RB + 1.17 × HBD + 0.44 × PSA_fraction

Where: clogP is Crippen logP, MW is molecular weight (Da), Arom is aromatic proportion, RB is rotatable bonds, HBD is H-bond donors, PSA_fraction is polar surface area fraction.

ESOL achieves RMSE ≈ 0.9–1.0 log units, which is close to the inter-laboratory reproducibility limit for solubility measurements (~0.6–0.8 log units).

Implementing ESOL in Python

from rdkit import Chem
from rdkit.Chem import Descriptors, rdMolDescriptors
import numpy as np

def esol_features(smiles: str) -> dict:
    """
    Compute the six physicochemical descriptors used in the ESOL model.
    Returns a dict of feature name -> value, or None if SMILES is invalid.
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    mw = Descriptors.MolWt(mol)
    clogp = Descriptors.MolLogP(mol)
    n_hbd = Descriptors.NumHDonors(mol)
    n_rotatable = Descriptors.NumRotatableBonds(mol)
    psa = Descriptors.TPSA(mol)

    # Aromatic proportion: fraction of heavy atoms in aromatic rings
    n_arom = sum(1 for atom in mol.GetAtoms() if atom.GetIsAromatic())
    n_heavy = mol.GetNumHeavyAtoms()
    arom_prop = n_arom / n_heavy if n_heavy > 0 else 0.0

    # PSA fraction
    psa_fraction = psa / (psa + 26.0)  # Delaney's normalization

    return {
        "MolWt": mw,
        "clogP": clogp,
        "NumHBD": n_hbd,
        "NumRotBonds": n_rotatable,
        "AromaticProportion": round(arom_prop, 4),
        "PSA_fraction": round(psa_fraction, 4),
    }


def esol_predict(smiles: str) -> float:
    """
    Predict logS using the ESOL linear model (Delaney 2004 coefficients).
    Returns predicted logS or np.nan for invalid SMILES.
    """
    feats = esol_features(smiles)
    if feats is None:
        return np.nan
    logS = (
        0.16
        - 0.63 * feats["clogP"]
        - 0.0062 * feats["MolWt"]
        + 0.066 * feats["AromaticProportion"]
        - 0.74 * feats["NumRotBonds"]
        + 1.17 * feats["NumHBD"]
        + 0.44 * feats["PSA_fraction"]
    )
    return round(logS, 2)


# Test
for smi, name in [
    ("c1ccccc1", "benzene"),
    ("CC(=O)Oc1ccccc1C(=O)O", "aspirin"),
    ("CC12CCC3C(C1CCC2O)CCC4=CC(=O)CCC34C", "testosterone"),
]:
    print(f"{name}: logS = {esol_predict(smi)}")

Feature Engineering for ML Models

Beyond ESOL's six features, a richer descriptor set significantly improves ML model performance. Common feature groups for solubility prediction:

Full Feature Matrix Construction

from rdkit import Chem
from rdkit.Chem import Descriptors, AllChem
import numpy as np
import pandas as pd

DESCRIPTOR_NAMES = [
    "MolWt", "MolLogP", "MolMR", "TPSA",
    "NumHDonors", "NumHAcceptors", "NumRotatableBonds",
    "RingCount", "NumAromaticRings", "FractionCSP3",
    "HeavyAtomCount", "NumValenceElectrons",
    "MaxPartialCharge", "MinPartialCharge",
]

def compute_descriptors(smiles: str) -> np.ndarray | None:
    """Return a 1D array of physicochemical descriptors, or None if invalid."""
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None
    values = []
    for name in DESCRIPTOR_NAMES:
        try:
            values.append(getattr(Descriptors, name)(mol))
        except Exception:
            values.append(np.nan)
    return np.array(values, dtype=np.float32)


def smiles_to_feature_matrix(
    smiles_list: list,
    include_fingerprint: bool = True,
    fp_bits: int = 1024
) -> np.ndarray:
    """
    Build feature matrix combining physicochemical descriptors
    and optionally Morgan fingerprints.
    """
    desc_list, fp_list = [], []
    valid_mask = []

    for smi in smiles_list:
        desc = compute_descriptors(smi)
        mol = Chem.MolFromSmiles(smi)
        if desc is None or mol is None:
            valid_mask.append(False)
            continue
        desc_list.append(desc)
        valid_mask.append(True)
        if include_fingerprint:
            fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=fp_bits)
            arr = np.zeros(fp_bits, dtype=np.uint8)
            fp.ToNumpyArray(arr)
            fp_list.append(arr)

    if not desc_list:
        return np.empty((0, len(DESCRIPTOR_NAMES) + fp_bits))

    X_desc = np.vstack(desc_list)
    if include_fingerprint and fp_list:
        X_fp = np.vstack(fp_list)
        return np.hstack([X_desc, X_fp])
    return X_desc

Training a Random Forest Solubility Model

import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score
from sklearn.metrics import mean_squared_error, r2_score

# Load the Delaney/AqSolDB dataset
# Available from: https://github.com/deepchem/deepchem/tree/master/datasets
df = pd.read_csv("delaney.csv")       # columns: smiles, measured_log_solubility_in_mols_per_litre
print(f"Dataset: {len(df)} compounds")
print(df["measured_log_solubility_in_mols_per_litre"].describe())

# Build feature matrix
X = smiles_to_feature_matrix(df["smiles"].tolist(), include_fingerprint=True)
y = df["measured_log_solubility_in_mols_per_litre"].values

print(f"Feature matrix: {X.shape}")

# 5-fold cross-validation
rf = RandomForestRegressor(
    n_estimators=300,
    max_features=0.4,
    min_samples_leaf=2,
    random_state=42,
    n_jobs=-1
)

rmse_scores = np.sqrt(-cross_val_score(rf, X, y, cv=5, scoring="neg_mean_squared_error"))
r2_scores = cross_val_score(rf, X, y, cv=5, scoring="r2")
print(f"5-fold CV RMSE: {rmse_scores.mean():.3f} ± {rmse_scores.std():.3f}")
print(f"5-fold CV R²:   {r2_scores.mean():.3f} ± {r2_scores.std():.3f}")

# Train final model
rf.fit(X, y)

# Feature importance (top 10)
desc_names = DESCRIPTOR_NAMES + [f"fp_{i}" for i in range(1024)]
importances = pd.Series(rf.feature_importances_, index=desc_names)
print("\nTop 10 important features:")
print(importances.nlargest(10).to_string())

Benchmark Results: What to Expect

GNN Approach with DeepChem

import deepchem as dc
import numpy as np

# Load the Delaney dataset with scaffold split (DeepChem built-in)
tasks, datasets, transformers = dc.molnet.load_delaney(
    featurizer="GraphConv",  # Graph-based featurization
    splitter="scaffold"      # Scaffold-based split for realistic evaluation
)
train_dataset, valid_dataset, test_dataset = datasets

print(f"Train: {len(train_dataset)} | Valid: {len(valid_dataset)} | Test: {len(test_dataset)}")

# Build an AttentiveFP model
model = dc.models.AttentiveFPModel(
    n_tasks=1,
    mode="regression",
    num_layers=4,
    num_timesteps=3,
    graph_feat_size=200,
    dropout=0.2,
    learning_rate=1e-3,
    batch_size=64,
    model_dir="solubility_afp"
)

# Train
model.fit(train_dataset, nb_epoch=60, checkpoint_interval=20)

# Evaluate
metric = dc.metrics.Metric(dc.metrics.pearson_r2_score)
rmse_metric = dc.metrics.Metric(dc.metrics.rms_score)

print("Test R²:  ", model.evaluate(test_dataset, [metric], transformers))
print("Test RMSE:", model.evaluate(test_dataset, [rmse_metric], transformers))

# Predict new compounds
new_smiles = ["CC(=O)Oc1ccccc1C(=O)O", "c1ccc(F)cc1", "CC(N)Cc1ccccc1"]
new_dataset = dc.data.NumpyDataset(
    X=np.array(new_smiles),
    ids=new_smiles
)
featurizer = dc.feat.MolGraphConvFeaturizer()
X_new = featurizer.featurize(new_smiles)
new_dataset = dc.data.NumpyDataset(X=X_new, ids=new_smiles)
predictions = model.predict(new_dataset)
for smi, pred in zip(new_smiles, predictions.flatten()):
    print(f"{smi}: predicted logS = {pred:.2f}")

Dataset Curation Pitfalls

Solubility data is notoriously noisy. Common issues you must address before training:

• Protocol heterogeneity: nephelometric (kinetic) solubility can be 1–2 orders of magnitude higher than thermodynamic (shake-flask) solubility. Never mix them without a correction term or separate models.
• pH dependence: ionisable compounds have different solubility at pH 2 (stomach) vs. pH 7.4 (blood). Always check the pH of measurement.
• Duplicate SMILES with inconsistent values: the same SMILES from different sources may have different measurement values. Cluster duplicates and take the median, or flag high-variance entries.
• Salts and counterions: many databases list Na+, Cl− etc. in SMILES. Strip to the largest organic fragment before modeling.
• Tautomers and standardisation: standardise all SMILES with a tool like MolVS or RDKit's rdMolStandardize before computing features.

from rdkit import Chem
from rdkit.Chem.MolStandardize import rdMolStandardize
import pandas as pd

def standardize_smiles(smiles: str) -> str | None:
    """
    Standardize a SMILES string:
    1. Parse to molecule
    2. Remove salts — keep largest fragment
    3. Neutralize charges
    4. Canonicalize
    """
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        return None

    # Remove salts
    remover = rdMolStandardize.FragmentRemover()
    mol = remover.remove(mol)

    # Neutralize
    uncharger = rdMolStandardize.Uncharger()
    mol = uncharger.uncharge(mol)

    # Tautomer canonicalization
    te = rdMolStandardize.TautomerEnumerator()
    mol = te.Canonicalize(mol)

    return Chem.MolToSmiles(mol)


# Apply to raw dataset
df = pd.read_csv("raw_solubility.csv")
df["std_smiles"] = df["smiles"].apply(standardize_smiles)
df = df.dropna(subset=["std_smiles"])

# Remove duplicates — keep median value per canonical SMILES
df = df.groupby("std_smiles")["logS"].median().reset_index()
print(f"After deduplication: {len(df)} unique compounds")

Deployment Considerations

• Applicability domain: solubility models are not reliable for molecules outside their training domain. Use leverage-based (hat matrix) applicability domain, or simply check Tanimoto similarity to the nearest training set neighbour. Flag predictions for compounds with T < 0.3 as unreliable.
• logS vs. kinetic solubility: if your assay gives nephelometric (kinetic) data, train a separate kinetic solubility model. Do not mix with thermodynamic logS models.
• Temperature and co-solvents: DMSO from solvent stocks affects solubility. Meta-data on measurement conditions should be stored alongside predictions.
• Ensemble predictions: serving the mean and standard deviation of an ensemble of 10 RF models provides uncertainty estimates cheaply.

Key Takeaways

• The Delaney ESOL six-descriptor linear model (RMSE ~1.0) is the classic baseline — implement it as a sanity check before adding ML complexity.
• Combining Crippen logP and Morgan fingerprints in a Random Forest closes most of the gap with deep learning at minimal complexity cost.
• GNNs (AttentiveFP) achieve RMSE ~0.6 on Delaney scaffold splits; the realistic improvement over RF is 0.1–0.15 log units.
• Data curation is the highest-leverage activity: mixing kinetic and thermodynamic solubility measurements is the most common source of poor model performance.
• Always report scaffold-split performance, not random-split performance — random splits inflate R² by ~0.05–0.10 due to scaffold leakage.

Try it live

Explore a related model hands-on: open Solubility Product Ksp — Precipitation & the Common-Ion Effect and change its parameters while it runs. Everything stays in your browser.