Molecular docking is a computational procedure that attempts to predict noncovalent binding of a macromolecule (e.g., a membrane receptor) and a small molecule (e.g., a ligand) [6]. To study the interaction between OR2W1 and its ligands, we started with the receptor responses to 42 chemically diverse odorants reported previously by Adipietro et al. [4]. The amino acid sequence for OR2W1 (UniProt: Q9Y3N9) was obtained from UniProt KnowledgeBase (https://www.uniprot.org/) [7]. To pick out the amino acid sequences of 7TMRs with a higher BLAST score (Bits) than 45 or a lower E-value than 1e-07 [8,9], we analyzed the amino acid sequence of OR2W1 against the locally-built BLAST database of 7TMR amino acid sequences registered in the RCSB Protein Data Bank (PDB) [10,11]. Based on their similarity results for OR2W1, four PDB entries of three 7TMRs were selected as the experimental templates for the homology modeling of OR2W1: human serotonin receptor (HTR2A), bovine rhodopsin (Rho), and two turkey β1-adrenoceptor (ADRB1) PDB models.

In order to apply the activated structure of the receptor as a template, we adopted the PDB model of each receptor binding to an agonist or the corresponding G protein: HTR2A (PDB code: 6WHA; BLASTP score [E-value]: 50.4 bits [1e-08]), Rho (PDB code: 5TE3; BLASTP score [E-value]: 48.5 bits [4e-08]), ADRB1 (PDB codes: 6IBL, 6H7J; BLASTP score [E-value]: 48.5 bits [4e-08]). For these four PDB models, models were built from multiple templates using MODELLER (R. 9.25) [12]. Multiple sequence alignment (MSA) of the four PDB models was executed using Clustal Omega [13]. We applied the MSA result to MODELLER and the application automatically combined these four templates to build the model for OR2W1 using information from multiple templates to build the 3D structural model of OR2W1. After confirming the 3D model of OR2W1 with a Ramachandran plot, it was used for the subsequent docking experiment.

In addition to the confirmation of the 3D model of OR2W1, tools such as HMMTOP [14] and Phobius [15] were used to determine the transmembrane region of the OR. Hydrophobicity analysis was also conducted to verify the transmembrane region using the hydrophobicity scale from Kyte and Doolittle [16].

Automated molecular docking is widely used for the prediction of receptor-ligand complexes in interaction analysis and molecular design. There are several freely available programs for molecular docking analysis, including AutoDock4 [17], AutoDock Vina [6], idock [18], and smina. Smina was created as a fork of AutoDock Vina to provide enhanced support for minimization and scoring [19]. Based on a report that the above scoring methods fared worse, on average, than simply using the output from smina alone [20], subsequent docking and scoring experiments were performed using the smina program.

We utilized AutoDockTools4 (ADT4) [17], which accompanies AutoDock4, to adjust the experimental conditions of the 3D docking space inside the membrane receptor. The binding site grid box was visually defined for each receptor by employing the grid setting feature of ADT4. To compare whether the experiment using smina was performed properly, we calculated the agonist binding affinity and conformation of human adenosine receptor A2A (AdoRA2A) as well as OR2W1. We downloaded the 3D structural model of AdoRA2A (PDB code: 5G53) from RCSB PDB and removed its ligand to obtain a structural model of the receptor in an activated conformation. The AdoRA2A model was then subjected to molecular docking with N-ethyl-5'-carboxamidoadenosine (NECA) using smina.

Each run with smina was executed using the default parameters with the exception of the 3D coordinates of the search space, so that the program output nine docking poses for each run. As smina accepts a ligand in the SDF format, the 3D SDF files of 22 small molecules were downloaded from PubChem of NCBI and used for the docking experiment with the corresponding receptor.

We hypothesized that there would be a difference in binding affinity between OR2W1 agonists and inverse agonists, and tried to verify the significance of the difference. Based on the OR2W1 ligand binding data from the previously reported experimental results [4], we statistically compared the binding affinity scores of agonists and inverse agonists of OR2W1 using a protein modeling and docking experiment. We implemented the statistical evaluation using the R statistical package.

Among the ligands tested for binding to OR2W1, (+)-carvone produced the greatest response to the receptor [4]. The homology model of human OR2W1 generated by MODELLER is shown in Fig. 1. This structure was considered to be stable according to the Ramachandran plot of the preferred model structure (Supplementary Fig. 1). Additionally, the amino acid sequences of the seven transmembrane regions of the OR2W1 model were compared to those produced by sequence prediction tools such as HMMTOP [14]. The transmembrane domains from the OR2W1 model and the two predicted results of sequence prediction were aligned properly, supporting the validity of the 3D model of OR2W1. Even in the worst case of transmembrane region alignment, only four amino acid residues were short at the N-terminus of the second transmembrane helix compared to those of HMMTOP. Moreover, the amino acid residues of the extracellular side of the receptor transmembrane region almost matched the predicted results of the sequence prediction tool.

As a criterion for comparison, we carried out molecular docking of AdoRA2A and NECA using smina according to the method described above. The best binding affinity score obtained from the docking experiment of AdoRA2A and NECA using smina was ‒8.7 kcal/mol. This docked pose of NECA with AdoRA2A is shown in Fig. 2 (scaled ball and stick).

As shown in Fig. 2, the conformation of the docking model based on the results of smina showed almost the same structure as that previously obtained by X-ray crystallography [21]. The amino acid residue of Asn253 (6.55; Ballesteros-Weinstein nomenclature [22], in yellow) in AdoRA2A, which interacts with the amino group of NECA, was found to be critical for activation of the receptor in previous studies of the 3D structure of AdoRA2A [23-25].

The homology model of OR2W1 generated with MODELLER and its best docking conformation with (+)-carvone is presented in Fig. 3. The docking conformation of OR2W1 and (+)-carvone showed a minimum binding energy of ‒7.2 kcal/mol (Table 1).

As shown in Fig. 3, (+)-carvone is surrounded by the third, sixth, and seventh transmembrane domains of the receptor, and is located very close to the Tyr252 amino acid residue in the sixth transmembrane region. This Tyr252 residue is not only at a position similar to the Asn253 residue of AdoRA2A, but also appears to form pi-pi stacking with the agonist (+)-carvone.

Smina predicts the binding of a receptor to a ligand and outputs the corresponding binding affinity along with its positional conformation. Although the challenge of selecting the correct docked pose remains [19], the best model selected from each smina prediction not only had a zero distance of the lower and upper bound root-mean-square deviation from the best mode, but also the lowest binding affinity energy between receptor and ligand among the nine suggested modes. The example docking model of AdoRA2A-NECA shown in Fig. 2 was also adopted as described above. The binding affinity of the best model obtained from each docking experiment is shown in Table 1. The mode that scored best according to the smina scoring function was chosen as the representative mode and its affinity score was subjected to further statistical analysis.

The average binding affinities of the 12 agonists and 10 inverse agonists were ‒6.325 and -4.9, respectively. These two binding affinity groups were checked for normality using the Shapiro-Wilk test, and the resulting p-values of the agonist and inverse agonist groups were 0.5229 and 0.08436, respectively. Therefore, the binding affinity data from each group were considered to follow a normal distribution. The F-test was performed to determine whether the variances of the two groups were homogeneous, and as a result, the p-value was 0.1475. Accordingly, the two-sample t-test assuming equal variances was conducted to evaluate the difference between the average values of the two groups, and the difference was found to be significant (p = 0.01019, two-sided).

In addition, as a result of performing the Wilcoxon rank sum test under the assumption that the data did not follow a normal distribution, the p-value was found to be 0.01333. This suggests that the median values of the two distributions are not equal. From the above results, it seems reasonable that the modeled binding affinity value of the agonists to OR2W1 is lower than that of the OR2W1 inverse agonists.