Structure-Kinetics Relationships of Opioids from Metadynamics and Machine Learning

The nation’s opioid overdose deaths reached an all-time high in 2021. The majority of deaths are due to synthetic opioids represented by fentanyl. Naloxone, which is a FDA-approved reversal agent, antagonizes opioids through competitive binding at the μ-opioid receptor (mOR). Thus, knowledge of opioid’s residence time is important for assessing the effectiveness of naloxone. Here we estimated the residence times of 15 fentanyl and 4 morphine analogs using metadynamics, and compared them with the most recent measurement of the opioid kinetic, dissociation, and naloxone inhibitory constants (Mann, Li et al, Clin. Pharmacol. Therapeut. 2022). Importantly, the microscopic simulations offered a glimpse at the common binding mechanism and molecular determinants of dissociation kinetics for fentanyl analogs. The insights inspired us to develop a machine learning (ML) approach to analyze the kinetic impact of fentanyl’s substituents based on the interactions with mOR residues. This proof-of-concept approach is general; for example, it may be used to tune ligand residence times in computer-aided drug discovery.


List of Figures
S9 Interaction energies between selected mOR residues and the R1 group . . S23 S10 Interaction energies between selected mOR residues and the R2 group . . S24 S11 Interaction energies between selected mOR residues and the R4 group . . S25 S12 Model qualities of the tree-based regression models (all 24 features) . . . . S26 S13 Model qualities of the tree-based regression models (15 features) . . . . . . S27 S14 Correlation between the interaction energies and the calculated residence

Methods and protocols
Ligand similarity matrix calculations. The structural similarities of the opioids were evaluated using Tanimoto coefficients (T c ) S1 calculated by the RDKit package S2 We first converted the SMILES (Simplified molecular-input line-entry system) strings of molecules to the 1,024-bit Extended-Connectivity fingerprints with the radius of 2 bonds. S3 Once the fingerprints were obtained, T c between two molecules was calculated as S4  and MOP) were studied in this work. As in our previous metadynamics study of FENbound mOR, S5 the starting structures of mOR in complex with the fentanyl analogs were modified from a representative snapshot taken from the global free energy minimum sampled by the WE simulations of FEN-bound mOR. S6 The latter simulations were initiated from the MD relaxed docked structure generated using the X-ray structure of the agonist BU72-bound mOR (PDB entry 5C1M) S7 as a template and root-mean-square deviation (RMSD) from it as a progress variable. The same X-ray structure was used as a template to dock MOP and BUP into mOR. For docking NLX and NTX into mOR, the antagonist β-FNA bound X-ray crystal structure of mOR (PDB entry 4DKL) S8 was used as a template. Although the crystal structures are of mouse's mOR, the amino acid sequences of the human and mouse receptors are nearly identical. The protonation states of mOR were taken from our previous work, which were determined by the membrane-enabled hybrid-solvent continuous constant pH molecular dynamics (CpHMD) method with the pH replica-exchange protocol. S9-S11 It was found that D114 (D 2.50 ) is charged and H297 (H297 6.52 ) samples the neutral Hid and Hie tautomer states at physiological pH. Since the previous metadynamics simulations of FEN-mOR dissociation S5 found that the Hid state S4 gave one order of magnitude slower kinetics than the Hie state, resulting in a better agreement with experiment, we adopted it in the current simulations. The rest of the titratable sites were found to be in the standard protonation states (Asp/Glu is charged, His is in the neutral Hie state, Lys is charged, and Cys is neutral) and were adopted as such. All fentanyls and morphinans carry a charge of +1. Following the previous protocols, S5,S6 each of the aforementioned ligand-mOR complex structures was embedded in the POPC (1-palmitoyl-2-oleoyl-glycero-3-phosphocholine) lipid bilayer and solvated with water and Na + /Cl − ions. The resulting system underwent stepwise MD equilibration and finally unrestrained MD to relax the conformation for at least 200 ns prior to the metadynamics simulations (details see our previous work S5,S6 ).

Molecular dynamics protocol.
All simulations were performed with NAMD2. S12 The protein and lipids were represented by the CHARMM36 protein and lipid force fields, respectively. S13,S14 The CHARMM modified TIP3P model was used to represent water. S15 The fentanyls and morphinans were represented by the CGenFF force field (version 3.0.1) obtained through the Paramchem server. S16,S17 All bond and angles involving hydrogen atoms were constrained using the SHAKE algorithm S18 to allow an integration timestep of 2 fs. The simulations were performed under periodic boundary and constant NPT conditions in a flexible cell with a constant xy ratio. The temperature was maintained at 310 K by Langevin dynamics with a damping coefficient γ of 1 ps −1 , and the pressure was controlled at 1 atm by the Nosé-Hoover Langevin piston method. S19,S20 The van der Waals interactions were smoothly switch off from 10 to 12 Å using a switching function.
The particle mesh Ewald (PME) method S21 was used to calculate long-range electrostatic energies with a sixth-order interpolation and a grid spacing of 1 Å. Most of the analyses were performed using VMD S22 and PLUMED. S23 Metadynamics simulations of ligand-mOR dissociation. Metadynamics simulations were performed using the Colvars module S24 in NAMD2, S12 which implements the well-S5 tempered metadynamics method. S25,S26 The ligand-mOR dissociation process was accelerated through the deposition of a time-dependent Gaussian potential along two collective variables CV1 and CV2. Following our previous work, S5 the CV1 was the ligand-mOR contact number, and CV2 was the ligand center-of-mass (COM) z position relative to that of the C α atoms in the putative binding pocket according to the X-ray structure of BU72bound mOR. S7 These CVs do not presume a particular unbinding pathway. The residues in the binding pocket include Y75, Q124, N127, W133, L144, D147, Y148, M151, F152, L232, K233, V236, A240, W293, I296, H297, V300, W318, H319, I322, and Y326. The ligand-mOR contact number (CN) is defined using a differentiable switching function where d ij is the distance between the heavy atom i in ligand and j in mOR and R cut is the effect cutoff distance which was set to 4.5 Å. Following the previous ligand-protein dissociation simulations, S5 15 independent metadynamics simulations with different random velocity seeds were carried out for each ligand. A cutoff of 15 Å for CV2 (z position) was used to define the unbound state. The robustness of this cutoff in range of 10-20 Å was examined previously. S5 There were a total of 285 trajectories, with the aggregate sampling time of ∼15 µs. The Gaussian weight parameter W (hill height) was set to 0.5 kcal/mol and the width σ was set to 5 for CV1 and 0.5 Å for CV2. The tempering parameter ∆T was set to 14T, where T is the simulation temperature 310 K. The PLUMED program S23 was used to calculate the unbiased free energy profiles and the unbiased dissociation times (see below).
Estimation of ligand-mOR residence times. Following Casanovas S27 and others, S28 we first converted the individual metadynamics exit time t exit,meta to the real time t real by S6 summing up the time steps rescaled at each step : where i is the step number for bias deposition, N is the total number of bias deposition steps, δt is the time step for bias deposition and V meta (t i ) is the bias potential at time t i .
Following Salvalaglio et al., S29 the dissociation time (τ ) was estimated by fitting the empirical cumulative distribution (CDF) function which represents the probability of observing at least one dissociation event by time t, to the theoretical CDF for a Poisson process, where τ is the estimated residence time. The Kolmogorov-Smirnov (KS) test was used to test the null hypothesis that the sample of dissociation times extracted from metadynamics and a large sample of times randomly generated according to the theoretical probability density reflect the same underlying distribution. S29 The null hypothesis is rejected if pvalue < α, where α is typically chosen as 0.05. To estimate the errors of τ , we performed bootstrapping analysis with 10,000 samples. Following Palacio-Rodriguez et al, S30 we used the samples that passed the KS test, defined by p-value ≥0.05, for calculation of the mean and standard error (SE) of τ .
Calculations of protein-ligand or residue-substituent interaction energies. The interaction energy calculations were performed using the in-house tcl scripts for VMD. S22 The interaction energy was defined as the sum of electrostatic and van der Waals energies between the protein (or an individual residue) and ligand (or a substituent group).
All atoms and a dielectric constant of 4 were used in the calculations. For each ligand, the 15 trajectories (saved every 50 ps), with the initial 5 ns and the dissociated portion (CV2 greater than 15 Å, see above) discarded, were combined, resulting in 7,500-16,000 S7 frames.
Feature engineering. For ligand, the aforementioned trajectory frames were randomly sampled 100 times and each gave 2,000 snapshots. For these snapshots, the interac- V300-R4, W318-R4, and I322-R4) as features for a supervised ML analysis.

Machine learning analysis of kinetics modulators. For each ligand, 2000 frames
were randomly sampled from the 37,500-80,000 combined trajectory frames. From these frames the interaction energies of the 24 residue-substituent pairs were calculated and the average energies obtained from the reweighted distributions using the reweighting protocol in PLUMED (version 2.5) S23 were saved. This resulted in 24 average energies (features) for each ligand. Since we only have 15 residence times (for 15 ligands), a data augmentation "trick" was applied in which we repeat the above protocol 100 times, resulting in 2,400 (24×100) energies for 1500 ligands. At this stage, there are more parameters but we will address this issue later.
The above energies were used as features to train tree-based regression models for predicting the log 10 transformed calculated residence times log 10 (τ cal ) of fentanyls. The ML was performed using PyCaret S31 which implements six popular tree-based regression methods, Extra Trees, Random Forest, Decision Tree, Adaptive Boosting, Gradient Boosting, and Extreme Gradient Boosting (Table S3). The training and test sets were split in a 80:20 ratio and a 10-fold cross validation was used for model validation and hyperparameter tuning. After tuning the hyperparameters using random grid search for 1000 steps, feature importance scores were calculated. To reduce possible dataset bias, we ran the ML procedure (training, testing, and feature importance calculation) for a total of 100 trials; in each trial, the splitting for training and test sets was random and independent. Among those tested models (Table S3 and Fig. S12), the Extra Trees, Random Forest, Gradient Boosting, and Extreme Gradient Boosting gave similarly good performances based on multiple evaluation metrics, including the mean average error (MAE), root mean squared error (RMSE), coefficient of determination (R 2 ), and mean absolute percentage error (MAPE); it gave the second best score for the root mean squared logarithmic error (RMSLE). Thus, we calculated the importance scores and SHAP values of S9 the features for these four types of models.

S10
Supplemental tables S11  S30 The mean and standard error of the mean (SEM) values of τ cal were estimated from the success samples. The p-value refers to the average p-value calculated from the success samples. τ exp , K d and K i,NLX values were obtained by Mann et al. S32 While the SE values of K d and K i,NLX for each of the compounds was reported, the one of its τ exp was estimated from the reported 95% confidence intervals, S32 and is defined as (upper limit -lower limit)/3.92.   109 The model quality metrics listed are mean average error (MAE), root mean squared error (RMSE), coefficient of determination (R 2 ), root mean squared logarithmic error (RMSLE), and mean absolute percentage error (MAPE). These values were obtained based on 100 trials through PyCaret. S31 The first four methods are similar in performance although the extra trees method is slightly better.   Figure S4: Interactions between fentanyls and the mOR residues. a) Contact profiles for fentanyls with R4 and/or R1 modification with the relevant mOR residues. b) Contact profiles for fentanyls with R2 modification with the relevant mOR residues. P is the fraction of the analyzed trajectories that a compound interacted with a residue defined by a heavyatom distance cutoff of 5 Å. P max is the maximum value of contact fractions (P ) among all compounds.  Figure S5: Approximate free energy surfaces illustrating the D147 salt bridge and the H297 hydrogen bond with the rest of fentanyls with R4 and/or R1 modification. a-f) R piper-D147 , which is the distance between the piperidine nitrogen and carboxylate carbon atom of D147, illustrates the formation of the D147 salt bridge. R piper-H297 , which is the distance between the piperidine nitrogen and the nearest imidazole nitrogen of H297, illustrates the formation of the H297 hydrogen bond. The approximate free energies were calculated using the reweighting protocol S35 Figure S9: Violin plots of the reweighted interaction energies between selected mOR residues and the R1 group of different fentanyls. The maximum, median, and minimum are indicated by upper, middle and lower bars, respectively.

S18
S23 Figure S10: Violin plots of the reweighted interaction energies between the selected mOR residues and the R2 group of different fentanyls. The maximum, median, and minimum are indicated by upper, middle and lower bars, respectively.
S24 Figure S11: Violin plots of the reweighted interaction energies between the selected mOR residues and the R4 group of different fentanyls. The maximum, median, and minimum are indicated by upper, middle and lower bars, respectively.

Figure S12: Model qualities of the tree-based regression models (all 24 features).
Histograms of the root-mean squared error (RMSE), mean absolute error (MAE), and R 2 values were calculated over 100 independent ML trials.