A Fully In Silico Protocol to Understand Olfactory Receptor-Odorant Interactions

Abstract

Understanding olfactory receptor (OR)-odorant interaction is crucial for unraveling the molecular intricacies of smell, a sense that is essential for health and survival and has potential therapeutic applications. Nevertheless, the absence of comprehensive experimental data concerning ORs has significantly impeded the understanding of the structural dimensions of olfaction, thereby necessitating innovative approaches to elucidate the structural intricacies of ORs. In this study, we developed an in silico protocol to predict OR structures and study relevant odorant interactions using the OR51E2-propionate complex as a reference. We also developed a hybrid homology modeling strategy leveraging homologous Alphafold structures. This approach resulted in structures with better stability than Alphafold predicted models, as validated through molecular dynamics simulations. Our pipeline accurately replicated experimental findings for OR51E2 and was extended to three homologous ORs: OR51E1, OR51D1, and OR51G2. We used a total of 217 molecules from the M2OR database and key food odorants and applied K-nearest neighbor clustering, selecting a total of 78 representative molecules based on their proximity to cluster centroids for molecular docking studies. Our computational pipeline successfully verified over 25 previously established odorant-OR relationships, including the identification of potential interactions between OR51G2 and molecules such as trans-2-nonenal and acetyl glutamic acid. This framework provides an efficient method for predicting and characterizing potential OR-odorant pairs, streamlining the discovery process prior to experimental confirmation and advancing our understanding of OR binding mechanisms.

1

Flow-chart representing the methodology pipeline, divided into 5 phases - Structure Prediction, Binding Site Prediction, Ligand Selection, Molecular Docking, and Molecular Dynamics Simulations.

2

(a) RMSD values for Alpfafoldv2 and generic homology modeling performed using SwissModel and the cryo-EM structure of an activated Cholecystokinin A receptor (CCKAR)-Gi complex. On the right, showing the hybrid homology modeling structure for OR51E2, modeled using the Olfr78 structure from AFv2, hybrid structures from mouse and human ORs, and AFv3 structures with an AFv3-hybrid comparison. Pruned pairs are the subset of atom pairs retained after removing mismatched/unsuitable atoms, ensuring more accurate RMSD calculation (e.g., RMSD between 219 pruned atom pairs is 1.019 Å and across all 302 pairs: 2.238). (b) Ramachandran plots for the two protein models in focus, OR51G2. AFv2 models show outliers­(in red), while the hybrid model has its residue data within the constraints.

3

(a) Predicted binding sites in OR51E2 via PUResNet; (b) predicted binding sites of OR51E2 via various binding site prediction tools. (c) Principal component analysis (PCA) plot for the molecules shows the 5 clusters and their centroids for each data set obtained i.e. KFOs and M2OR.

4

(a) Binding energies of all odorants (key food odorants and molecules from M2OR db on each hybrid receptor model. (b) Representing the predicted binding interaction of propionate with OR51E2, CYS:178, and HIS:180, forming an H-bond is apparent in the 2-D diagram on the left, and the 3-D docking pose shows the formation of these interactions. (c) The left figure compares GNINA docking energies across four olfactory receptors, where lower scores indicate stronger binding. OR51E2 shows significantly weaker energies. Statistical analysis (Kruskal–Wallis, p = 1.72 × 10^–12) confirms significant receptor differences, with Dunn’s test showing OR51E2 binds significantly weaker than OR51D1 and OR51G2 (p < 0.0001). The Plot on the right presents the ROC curves for OR51E1 and OR51E2, the predictive performance of GNINA docking scores for receptor responsiveness. OR51E1 achieves a high AUC of 0.78, showing a strong correlation between docking scores and experimental responsiveness, whereas OR51E2 shows a significantly lower AUC of 0.35.

5

(a) Representation of the binding interaction of OR51D1 and pentyl ropionate, the highest interaction energy was obtained for this interaction and is solidified with the presence of two hydrogen-bond interactions, i.e., between the carbonyl-ARG:276, and the alkoxy-GLN:195. (b) Binding interaction between OR51G2 and 3,6-nonadienal. Regardless of the highest interaction energy with OR51G2, the OR-Odorant pair forms an unfavorable acceptor-acceptor bond; contrary to that, the same carbonyl oxygen forms a hydrogen bond with HIS:110, suggesting an avenue for stability. (c) Binding interaction between OR51E2 and butyric acid, carbonyl forms a hydrogen bond with HIS:180 and a carbon–hydrogen bond with VAL:179.

6

(a) Representation of the binding interaction of OR51E1 and citral, the highest interaction energy was obtained for this interaction and is solidified with the presence of two hydrogen-bond interactions, i.e., between the carbonyl-GLN:185, and the alkoxy-HIS:184. (b) Binding interaction between OR51E1 and geranylacetate. Regardless of the highest interaction energy with OR51G2, the OR-odorant pair forms an unfavorable acceptor-acceptor bond; contrary to that, the same carbonyl oxygen forms a carbon–hydrogen bond with SER:111, suggesting an avenue for stability. (c) Binding interaction between OR51G2 and acetyl glutamic acid, carbonyl forms a hydrogen bond in all its interactions and thus shows strong binding energy compared to other molecules.

7

(a) Each graph shows the RMSD for three runs each for both models in comparison to each other across a length of 100 ns simulations. The system for AFv2 is observed to be more flexible initially in the first 30,000 ps compared to an overall stable hybrid model. (b) Radius of gyration across the axes of the models shows that the hybrid model tends to have more movement at the beginning and end of the simulation. (c) Root mean squared fluctuations across the receptor sequence within the lipid bilayer show increased flexibility in the ECL regions of the hybrid model, suggesting increased dynamics of the loop structures that may facilitate interactions with other proteins or the extracellular environment. (d) SASA calculations across the length of the simulation signify that the hybrid model is slightly more exposed to the solvent in the environment i.e., water.

8

(a) Ramachandran plots for post-MD structures of the hybrid and Alphafold models. (b) Binding sites for hybrid and Alphafold models before and after MD simulations. The binding site for the hybrid model has increased relatively more than the Alphafold model by both volume and area.