Interpretable Machine Learning as a Key to Understanding BBB Permeability

A model that says “crosses” or “does not cross” is unfinished until it also says which property of the molecule pushed the decision. On a blood-brain barrier table that question is answerable, because every row already has a label and a block of chemical descriptors.

The table is B3DB, the curated set of Meng, Xi, Huang, and Ayers [1]. The 2021 paper reports 7807 compounds with a categorical BBB+ or BBB− label, and a smaller subset of 1058 compounds with a numeric logBB. This note uses the categorical extended file, which carries 1613 Mordred descriptors. The code pins one commit of that file so the columns do not drift under you. Later releases of the repository have added compounds. The paper’s count and the file you download are allowed to differ. Say which one you used.

The shell line uses a notebook shortcut. In a terminal, drop the exclamation mark. The commit hash is the pin.

wget https://github.com/theochem/B3DB/raw/87240af2b4e585d56f9681a6426af6b7f2940e96/B3DB/B3DB_classification_extended.tsv.gz

import gzip
import shutil

import pandas as pd

with gzip.open("B3DB_classification_extended.tsv.gz", "rb") as src:
    with open("B3DB_classification_extended.tsv", "wb") as dst:
        shutil.copyfileobj(src, dst)

df = pd.read_csv("B3DB_classification_extended.tsv", sep="\t")
df.head()

The label, then the descriptor block

The label column is BBB+/BBB-. Map the two strings and drop anything else, so a blank cell cannot become a positive example by falling through an else.

The Mordred columns are the 1613 fields that begin at ABC and run to the end of the table (mZagreb2 in this file). Select them by that name. A positional slice that starts at column 6 starts on the label itself, takes in identifiers and comments, and stops near VE1_Dzpe, hundreds of columns before TopoPSA. A dependence plot for TopoPSA is only meaningful if TopoPSA was one of the columns the forest saw. In this file it is column 1568, inside the Mordred block and outside that positional slice.

dropna(axis=1) then removes any descriptor that is missing for even one molecule. That is a blunt filter. One failed calculation deletes the column for the whole set. It is applied after the descriptor block has been isolated, so it cannot delete the label or shift the remaining columns into the wrong place.

label_map = {"BBB+": 1, "BBB-": 0}
df = df[df["BBB+/BBB-"].isin(label_map)].copy()
df["labels"] = df["BBB+/BBB-"].map(label_map)

descriptors = df.loc[:, "ABC":].apply(pd.to_numeric, errors="coerce")
descriptors = descriptors.dropna(axis=1)
descriptors.head()

Fit on 70 percent

The split is stratified so the BBB+ and BBB− fractions match in the train and test sets. Thirty percent is held out. random_state is fixed on both the split and the forest, so a second run of this notebook hits the same rows and the same trees. The original forest call left the trees unseeded.

from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, roc_auc_score
from sklearn.model_selection import train_test_split

X = descriptors.reset_index(drop=True)
y = df["labels"].reset_index(drop=True)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, shuffle=True, stratify=y
)

rf = RandomForestClassifier(random_state=42)
rf.fit(X_train, y_train)

y_pred = rf.predict(X_test)
print(classification_report(y_test, y_pred))

roc_auc = roc_auc_score(y_test, rf.predict_proba(X_test)[:, 1])
print("ROC AUC:", roc_auc)

classification_report wants the true labels first and the predictions second. The report in Figure 4 was printed with those arguments reversed, so the precision and recall it shows are exchanged relative to scikit-learn’s definition. A rerun of the code above will not match that image cell for cell.

The ROC AUC is computed from the predicted probability of the positive class. An AUC computed from the hard 0/1 predictions is a different statistic. This page does not quote a number for either, because the value belongs to the fit you just ran. Read the printout.

Which descriptors the forest used

Impurity importance is the forest’s own score: how much each descriptor reduced impurity across the trees. It favors descriptors that can split in many places. It is a reason to look, not a mechanism.

ranking = (
    pd.DataFrame({
        "features": X.columns,
        "Imp_values": rf.feature_importances_,
    })
    .sort_values("Imp_values", ascending=False)
)
ranking.head(15)

SHAP

SHAP attributes a prediction to the descriptors. The calls below match the interface in which a tree explainer returns one array per class, and index 0 is the first class. If your installed SHAP returns a single three-dimensional array instead of that list, index the class axis explicitly. Do not assume [0] still means what it meant when these figures were drawn.

pip install shap

import shap

explain = shap.Explainer(rf)
shap_values = explain.shap_values(X_test)

shap.summary_plot(shap_values, X_test)

shap.summary_plot(shap_values[0], X_test)

shap.dependence_plot("TopoPSA", shap_values[0], X_test)

Bottom line

TopoPSA came out at the top of the importance table in the run these figures record, and the dependence plot is there so that ranking can be looked at rather than quoted. Polar surface area is already a familiar constraint on CNS exposure, which makes the result plausible. It does not make it automatic. A different split, a different forest, or a descriptor block that accidentally excluded TopoPSA would have told a different story.

The useful pair is the held-out score and the attribution. The score says whether the forest separated BBB+ from BBB− on molecules it did not train on. SHAP says which descriptors carried that separation. Impurity importance is only the shortlist.

Reference

[1] Meng, F., Xi, Y., Huang, J., and Ayers, P. W. A curated diverse molecular database of blood-brain barrier permeability with chemical descriptors. Scientific Data, 2021, 8, 289. https://doi.org/10.1038/s41597-021-01069-5