The screening workflow

The general steps of NIB docking and VS are (1) the selection and preparation of the protein target 3D structure, (2) preparation of ligand 3D structures, (3) selection of the protein binding site centroid, (4) generation of the negative image of the binding site, (5) shape and/or charge comparison of the negative image and ligand 3D structures, and finally (6) either evaluation of the model and the workflow with benchmarking sets or investigation and selection of screened molecules for experimental testing (V Fig. 1).

The ligand 3D structure preparation step (2) is of high importance as the Panther/ShaEP NIB screening uses rigid structures of the ligands in the shape/charge comparison step. The ligand 3D structures may be prepared with a number of software and workflows and the ligands can be in either single or multiconformer sets. Four ligand preparation workflows with commercial, academic free or free software were found to produce vast differences in both the benchmarking results (V Tables 2 and 3) and the number of ligands that survive through both the ligand preparation (2) and screening (5) steps (V Table 1). The use of single conformations of the ligands can greatly increase the speed of the ligand preparation and screening phase (V Fig. 3) depending on the size of the multiconformer set. Surprisingly, in this COX-2 benchmarking example, the single conformer sets also fared better than multiconformer sets, although the effect was clearly case-specific (V Tables 2 and 3). However, the observation may be specific for the particular COX-2 benchmarking example and the settings used in the ligand preparation workflows. Even so, it might prove beneficial in NIB screening benchmarking to compare different ligand preparation workflows and the use of single versus multiconformer ligand 3D structures especially in the case of suboptimal results in the first rounds of benchmarking.

The selection and preparation of the protein target 3D structure in step (1) defines the frame for the NIB model. The quality of the protein structure and its binding cavity affect the shape and charge properties of the NIB model and how well it represents the target. The binding site geometry (V Fig. 5C), H-bonding

50

groups (V Fig. 1, 4 A, 5 B) and possible cofactors, ligands (V Fig. 1 and 4D–E) or water molecules can all be considered in the model generation. Protonation can be considered rigidly if the protons are added to the protein prior to the Panther NIB model generation, or when applicable, alternative rotations of protons are considered in the absence of protonation (V Fig. 5 B). The benchmarking demonstrated that the screening yields vary between NIB models created with different settings and NIB models of identical settings from the two structures (V Tables 2 and 3), the latter being a phenomenon that is common in molecular docking and structure-based VS.

Creating the initial NIB model with Panther is fairly simple in steps (3) and (4), and there are many options to modify and optimize the model. It was demonstrated with a few settings that varying the settings can affect the screening yield case-specifically (V Fig. 4). For example, the NIB model can be restricted to the volume of the bound ligand in the used protein structure, which can shape the NIB model more alike the ligand (V Fig. D–E). Changing the packing method (V Fig. 4G–H) or cavity detection radius (V Fig. 5 E) can also result in a different shape of the model as the model may reach different nooks in the binding cavity. One of the advantages of NIB Panther models is that the definition of the binding cavity can be very specific if the user wants to target a precisely defined area of the protein binding cavity or surface pocket (V Fig. 5A).

As the pocket representation is atomic, the model can also be manually augmented with fragments of small molecules in order to define more precisely the desired properties of the screened compounds (Jokinen et al. 2019).

The shape and charge comparison in step (5), in other words docking and screening or rescoring, can be performed using ligand conformers created ab initio with a separate ligand preparation workflow or from a prior molecular docking run. Here, the rigid NIB docking and screening using ab initio ligand conformers produced consistently better AUC values than prior flexible docking (Kurkinen et al. 2018) using PLANTS (Korb et al. 2007, 2009) (V Tables 3 and 4).

However, a combination of the two approaches using NIB rescoring (R-NIB) or consensus scoring of ligand conformers from PLANTS docking is even more successful than either NIB or PLANTS screening alone (V Tables 1–3) as shown also by prior studies (Kurkinen et al. 2018).

The utilized workflow from the ligand and protein structure preparation to NIB screening and rescoring was provided as command line inputs and help scripts in the supplementary material of the original publication. The example workflow is based on the DUD benchmarking set for COX-2. The provided workflow utilizes freely available Open Babel (O’Boyle et al. 2011), Panther (Niinivehmas et al. 2015), ShaEP (Vainio et al. 2009) and Rocker (Lätti et al. 2016) for ligand preparation, NIB model generation, docking/screening and VS benchmarking calculation/visualization steps, respectively. The Open Babel ligand preparation did not produce the best results from the benchmarking (V Tables 2 and 3) and, accordingly, it may be necessary to use another ligand preparation workflow for actual VS studies. However, the workflow is excellent for demonstrative or tutorial purposes as it is performed with open and relatively easy-to-use tools. Ligand conformers of previous (Kurkinen et al. 2018) PLANTS

(Korb et al. 2007, 2009) docking are also included for R-NIB rescoring and consensus scoring. Together with the thorough discussion of the NIB methodology, the provided workflow inputs can provide a starting point for NIB VS for anyone new to the method.

Computational protein structure-based methods were applied to evaluate and predict CYP mediated metabolism. The application of docking and MD simulations on small data sets of small molecules provided insight on the structural determinants of substrate and SOM selectivity on various CYP enzymes, and even offered a glimpse of the tentative structural basis of the differing inhibition mode of one ligand on two similar CYP enzymes. Metrics of 1) binding free energy using MMGBSA, 2) ligand site accessibility to the heme iron using atomic distances, and 3) ligand stability using the STDs of the atomic distances are suggested for future studies of MD-based CYP ligand binding mode and substrate SOM prediction. A detailed MD-based protocol is not provided as the predictions here were partly based on manual analysis of the simulations due to large changes of the binding poses during the simulations.

Novel profluorescent tool molecules for CYP enzymes were discovered, with two of them being isoform-selective in certain tissues. Molecular modelling aided in the identification of the most potent target CYPs for the tool molecules. Docking offered a view of the likely interactions that would facilitate the reaction leading to the fluorescent product. MD simulations explained differences in the catalytic activity of the fluorescence-producing reaction of the compounds with the aid of the metrics proposed for binding mode and SOM prediction. Finally, the MD simulations suggested that water interactions and access channel composition have an essential role in ligand binding in the CYP1 family.

MMGBSA offered an adequate frame for the prediction of ligand binding poses and substrate SOMs in CYP enzymes. The model systems of CYP2A6 and 2A5 in complex with small coumarin derivatives offered insights into how the MD/MMGBSA method fares on a data set where the enzymes and ligands are highly similar. However, it was demonstrated that binding energy estimation alone does not necessarily predict a metabolically active binding mode and thus a correct SOM prediction. The binding mode with the best interaction energetics might not be metabolically active as also the chemical nature, 3D location and the general orientation of the molecule also affect whether a reaction occurs in a certain binding mode (Cruciani et al. 2005, 2013, Huang et al. 2013, Tyzack et al.