How to build and validate 3D-QSAR models - Insights from Xu et al., 2020
A 3D-QSAR model is a map from aligned molecular fields to a potency. The fit tells you whether that map predicts compounds you held out. The contour map tells you where changing the molecule might help. Neither result survives an alignment that put the ligands on top of one another in the wrong frame.
The workflow here is the standard CoMFA and CoMSIA stack: dock or otherwise pose the ligands, calculate fields on a shared grid, fit them with partial least squares, and refuse the model if a scrambled label looks just as good. Xu and colleagues used that stack, together with docking and molecular dynamics, on thieno[3,2-b]pyrrole-5-carboxamide inhibitors of LSD1 [1]. The field settings below are the usual ones this kind of study is built from.
Put the potency on one scale
Collect IC50 values from assays you are willing to treat as comparable. Then convert to pIC50 before any fitting:
$$ \text{pIC}_{50} = -\log_{10}(\text{IC}_{50} \text{ in mol/L}) $$The unit is the whole conversion. An IC50 of 1000 nM is \(10^{-6}\) mol/L, so the pIC50 is 6. The same number left in nanomolar and passed to a base-10 logarithm is a different, wrong, scale. Do this conversion once, and do not mix Ki or percent inhibition into the same column without saying so.
Hold out about 20 to 25 percent of the compounds as the external test set. The training set is the remaining 75 to 80 percent. The test compounds have to sit in the same chemical space as the training compounds. A test set that is a different series is an extrapolation, and a pretty \(r^2\) on it does not mean what it says.
Alignment is the model
CoMFA and CoMSIA descriptors exist on a grid in space. Two ligands contribute to the same grid point only if you have already decided how their atoms correspond. That decision is the alignment.
A common route is to take the protein structure from the Protein Data Bank, dock the series with a program such as SYBYL or AutoDock, and align every ligand to a reference pose. Redock a ligand whose experimental pose you know. An RMSD of 2.0 Å or less between the docked pose and that experimental pose is the usual acceptance check for this workflow. It is a check on the docking setup, not a guarantee that every other ligand found the right pose.
If the reference pose is wrong, every field shifts with it. The PLS model will still quote an \(r^2\). The contours will describe the error.
Fields
CoMFA
Comparative molecular field analysis uses two fields, steric and electrostatic.
- Grid spacing: 2.0 Å.
- Probe: an sp³ carbon with a charge of +1.
- Energy cutoff: 30 kcal/mol, so a clash at one grid point cannot dominate the column.
CoMSIA
Comparative molecular similarity indices analysis adds hydrophobic, hydrogen-bond donor, and hydrogen-bond acceptor fields to the steric and electrostatic ones. The attenuation factor is 0.3. It sets how fast a field falls off with distance from the atoms. Changing it changes which grid points still see the ligand.
PLS, and the numbers it produces
Partial least squares is used because the grid columns are many, collinear, and noisy. Leave-one-out cross-validation on the training set chooses the optimal number of components (ONC) and produces \(q^2\), the cross-validated \(r^2\). That is an internal number. It is not the test set.
The model is then refit on the full training set with that component count. Report:
- \(r^2\), the fit to the training compounds.
- SEE, the standard error of that fit.
- The F statistic for the regression.
A training \(r^2\) far above \(q^2\) is a model that has used components the left-out compounds do not support.
Tests the model is allowed to fail
Predict pIC50 for the compounds that were held out, and compute the external predictive \(r^2\), written \(r^2_\text{pred}\).
Tropsha’s diagnostics look at that external prediction more closely than a single \(r^2\). The set includes \(r^2_0\), the slopes \(k\) and \(k'\), \(r_m^2\), and \(\Delta r_m^2\). In this workflow a model is treated as predictive when \(k\) is near 1 and \(r_m^2\) is above 0.5. Those two cuts do not replace looking at the plot of predicted against measured pIC50.
Then scramble the biological activities, rebuild, and repeat. A real relationship collapses: the randomized trials should show low \(q^2\) and low \(r^2\). If they do not, the original fit was willing to explain noise.
Reading the contours
Only after those checks do the colors become a design suggestion.
Hydrophobic and hydrogen-bond contours, when you calculated CoMSIA fields, are read the same way: a color marks a region where that property, in this alignment, moved with potency.
The numbers, in one place
| Metric | What it is answering |
|---|---|
| \(q^2\) | Did leave-one-out on the training set still predict? |
| \(r^2\) | How well does the final model fit the training compounds? |
| \(r^2_\text{pred}\) | How well does it predict the compounds you held out? |
| SEE | How large are the training residuals? |
| F | Is the regression distinguishable from no relationship? |
| \(r_m^2\), \(\Delta r_m^2\) | Tropsha’s closer look at the external predictions |
| RMSD | Did the docking setup reproduce a known pose? |
| Y-randomization | Does the fit collapse when the labels are scrambled? |
Bottom line
A 3D-QSAR study is not a hunt for a large \(r^2\). It is an aligned field model that has to predict a held-out set, fail when the labels are scrambled, and only then offer a contour as a place to change the molecule. The LSD1 paper is one worked example of that order [1]. The order is the part worth copying.
Reference
[1] Xu, Y., He, Z., Liu, H., Chen, Y., Gao, Y., Zhang, S., … and Yang, M. 3D-QSAR, molecular docking, and molecular dynamics simulation study of thieno[3,2-b]pyrrole-5-carboxamide derivatives as LSD1 inhibitors. RSC Advances, 2020, 10, 6927–6943.
