Structural and Energetic Insights Into the Interaction of Niacin With the GPR109A Receptor

Abstract

The transmembrane G-protein coupled receptor GPR109A has been previously shown to function as a receptor for niacin in mediating antilipolytic effects. Although administration of high doses of niacin has shown beneficial effects on lipid metabolism, however, it is often accompanied by disturbing side effects such as flushing, liver damage, glucose intolerance, or gastrointestinal problems. Thus, it is important to understand niacin-GPR109A interactions, which can be beneficial for the development of alternate drugs having antilipolytic effects with less or no side effects. To get into the structural insights on niacin binding to GPR109A, we have performed 100 nanoseconds long all-atom MD simulations of five niacin-GPR109A complexes (automatically docked pose 0, and randomly placed niacin in poses 1 to 4 in the receptor crevice) and analyzed using binding free energy calculations and H-bond analysis. Steered MD simulations were used to get an average force for niacin translocation between the bulk and the external crevice of the wild type and mutant (N86Y, W91 S, S178I, and triple mutant of all three residues) GPR109A receptors, as well as GPR109B (as a control that does not bind niacin). The H-bond analysis revealed that TMH3 residue R111 interacts with niacin in a total of 4 (poses 0 to 3) complexes, while residues C177, S178, and S179 contact niacin in complex pose 4, and all these complexes were energetically stable. According to steered MD simulations, all the GPR109A mutants and GPR109B required greater force than that of wild-type GPR109A to translocate in the external crevice, suggesting increased sterical obstacles. Thus, the residues N86 (at the junction of TMH2/ECL2), W91 (ECL2), R111 (TMH3), and ECL3 residues (C177, S178, S179) play an important role for optimal routing of niacin entry and to bind GPR109A.

Graphical Abstract

Energy minimized structures of GPR109A-niacin complex.

Introduction

G-protein coupled receptors (GPRs) are the most important proteins involved in diverse biological processes including visual sense, sense of smell, immune system, and mood/behavior regulation. GPRs are traditionally transmembrane proteins, sharing a common three-dimensional topology, ¹ consisting of 7 transmembrane helices (7TMs), with N-terminus outside and C-terminus inside. The helices are connected by 3 extracellular loops (ECLs), 3 intracellular loops (ICLs), and 1 cytosolic helix. To date, the biggest challenge in the GPR research is the non-availability of the structural data at the atomic level, considering the abundance of sequence information on GPRs. The percentage identity between different families of GPRs is very low and the protein structure prediction based on the small amount of structural information on the superfamily of GPRs, namely, that of the X-ray crystallographic structure of bovine rhodopsin^2-5 and 2 human β₂-adrenoreceptor (β₂AR) structures, ⁶ leads to a large involvement of sequence/structure alignment correction and various molecular modeling and MD simulations approach to fully understand the action of several ligands to the GPRs.

The receptor designated GPR109A (HM74A in humans and PUMA-G in mice) belongs to the GPRs A family. In addition to GPR109A, 2 other receptors, namely, GPR109B and GPR81, are known to be activated by hydroxyl-carboxylic acids that are generated as intermediate ligands during energy metabolism. These receptors are primarily expressed in adipocytes and immune cells and couple to G proteins of the Gαi family which mediate anti-lipolytic effects. The inhibition of fat cell lipolysis via activation of Gαi-coupled receptors and subsequent inhibition of cAMP formation was observed through the GPR109A receptor, and activation by niacin results in reduced hydrolysis of triglycerides to fatty acids. Recent evidence also suggested the role of GPR109A in vascular inflammation by alternative signaling via β-arrestins ⁷ and also its role in inhibiting macrophage proinflammatory response. ⁸ Even though GPR109A and GPR109B are structurally similar with ~89% of sequence identity, selective binding of niacin to GPR109A receptor but not GPR109B^9-11 has been of much interest in pursuing GPR109A agonists.^12,13 Knockout mouse model of GPR109A has demonstrated an abolition of decrease in fatty acids and triglyceride plasma levels or the anti-lipolytic effects of niacin as well as the major side effects, such as cutaneous vasodilation called flushing. Similarly, high doses of niacin administration have also shown similar flushing side effects. ¹⁴ Thus, niacin interaction with the GPR109A channel is highly intriguing and its structural binding interaction studies might reveal new insights into designing novel therapeutic agents possessing cardiovascular benefits. Previous studies on class A GPRs including site-directed mutagenesis and molecular modeling techniques have shown that the binding site for most of the small molecule agonist (biogenic amines) are located on TMH 3, 4, 5, 6, and 7 and the ECL2 region. ³ ,^15-22 However, in the case of GRP109A, the binding site of niacin was found be located in the TMH 2, 3, 7 and ECL2 regions. ²³

Recently, human GPR109A was modeled using the X-ray crystal structure of rhodopsin (PDB ID: 1HZX) available in protein data bank (PDB). ²³ Previously, to gain insight into the GPR109A-niacin interactions, 1 ns MD simulations study for niacin-bound GPR109A in a water-vacuum-water box has been performed. ²³ However, underpinning all mutational results on the rhodopsin generated model of the GPR109A receptor for niacin binding had certain limitations. MD simulations of membrane proteins in the lipid bilayer are well developed and are readily available in many MD simulations packages: AMBER, ²⁴ CHARMM, ²⁵ NAMD, ²⁶ and GROMACS. ²⁷ Moreover, it is possible to run the MD simulations for a longer time scale, that is, up to a microsecond (μs) or more. ²⁸ Here, a key assumption is that, although the site-directed mutagenesis was carried out to recognize the binding site, the positions of the several identified functional residues were observed and interpreted only based on the vacuum-water-box of the GPR109A-niacin complex model. This setup does not exclude the possibility of misinterpreting the interaction of even the most conserved residues regardless of their identity or properties and their involvement in the entry mechanism of the ligand. Here, we review previously published experimental data by Tunaru et al ²³ and coupled with our new computational exploration of GPR109A using homology modeling and MD simulations of niacin binding.

Methods

Target sequence retrieval and template selection

To find the appropriate template, the primary amino acid sequence of GPR109A and GPR109B was submitted to different protein structure prediction servers: GPCR-I-TASSER, ²⁹ Phyre2, ³⁰ SWISS-MODEL, ³¹ and HHpred. ³² After careful inspection of all the models predicted from the webservers, based on the C-score, we chose the best structures predicted from GPCR-I-TASSER. A confidence score by GPCR-I-TASSER was used for estimating the quality of the predicted models. A higher value signifies a model with a high C-score and TM-score which specifies the quality of the model based on structural similarity. In parallel, we also generated GPR109A structural model from bovine rhodopsin (PDB ID: 1HZX) as mentioned in a previous study ²³ to observe the key structural differences between the GPCR-I-TASSER model and the rhodopsin-based model. The Modeller^33-35 program has been specifically used for generating rhodopsin-based 3D structure of GPR109A by the satisfaction of spatial restraints. For both GPCR-I-TASSER and Modeler-built protein structures, explicit hydrogens were generated automatically by the respective modeling software. The protonation state of every titratable residue was set based on pKa estimations by PROPKA. ³⁶ The refinement and validation of the models were performed by PROCHECK ³⁷ and the quality of the modeled structure based on nonbonded interactions was determined by ERRAT. ³⁸

Molecular docking

The binding potential of niacin against the GPR109A model built by GPCR-I-TASSER was computationally predicted by using the molecular docking method called AUTODOCK 4.2 to obtain the initial protein-ligand complexes. ³⁹ The initial three-dimensional structure of niacin required for docking studies with GPR109(A/B) receptors was retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/compound/938). During the docking process, the torsional angles for niacin were held flexible, whereas the torsional angles for the GPR109A were held rigid except for the residues Y86, S91, I178, R111, and R251 that have been previously implicated in ligand binding. ²³ By allowing the torsional degrees of freedom for the ligand, we intended to facilitate the exploration of the conformational space of the niacin and the binding region.^40,41 Polar hydrogen atoms and Kollman’s united atom partial charges were assigned for the protein using AutoDock tools.^35,39 The grid size was set to 70X70X70 points with a grid spacing of 0.375 Å. The grid box was centered by taking into consideration the coverage of the residues Y86, S91, I178, R111, and R251. For all the docking calculations, we applied the Lamarckian genetic algorithm (LGA) specified in AutoDock. A total of 25 000 000 steps of maximum energy evaluation were performed with a population size of 300, while the total number of independent runs was fixed to 150. To group the similar conformations or “clusters” based on their lowest energy conformations and their RMSD to one another, we used the default clustering algorithm described in the ADT/AutoDock tools. ³⁹ In the case where all the AutoDock clusters results were docked at 2 Å RMSD and the positions differed by less than 2 Å, these models were taken as identical and represented by the energetically top-ranked structures, as the energy differences within the docked structures placed in the same cluster were generally small under these assumptions.

MD simulations

The NAMD package ²⁶ was used to perform the MD simulations on the predicted protein-ligand complex. Assembly of the simulation cell, visualization, and analysis of the results were done using custom Tcl scripts in the VMD v1.9.3. ⁴² The lipid membrane used for the simulation run was taken from the pre-equilibrated (200 ns) “average composition” yeast membrane from Jeff Klauda lab (https://user.eng.umd.edu/~jbklauda/memb.html), which was very close to the “generic” plasma membrane of eukaryotes. The membrane consists of 270 lipid molecules, including 60 cholesterols, 100 DOPC, 20 DPPC, 20 POPA, 60 POPE, and 10 POPS residues. To preserve the proper lipid composition after embedding the protein into the membrane patch, we have kept all the lipids, so to make space for the protein, we have relocated all the lipids that were overlapping with the channel and moved those lipids laterally to other spots in the membrane. The total area of the membrane was adjusted (stretched) to accommodate the channel, resulting in a rectangular simulation cell of about 92.0 × 92.0 × 123.6 Å. Electro neutrality of the whole system was maintained by adding a total of randomly placed 87 K⁺ and 70 Cl⁻ ions to GPR109A-niacin complex and 87 K⁺ and 69 Cl⁻ ions to GPR109B-niacin complex up to the equivalent of 150 mM salt concentration. The systems were hydrated with ~20,000 water molecules, bringing the total system size to ~98,000 atoms. After insertion of the protein into the lipid bilayer, the system has been energy-minimized (5000 steps, conjugate gradient method) and simulated with harmonically restrained protein backbone for 10 ns for both the GPR109(A/B) protein and the niacin molecule. All the MD simulations were performed using the NPT ensemble using CHARMM36 force field and TIP3P water model. ⁴³ Constant pressure (1 atm) and 310 K temperature were supported by Langevin dynamics. A periodic boundary setting with a flexible cell was maintained with the cutoff distance applied for non-bonded interactions taken as 12 Å and particle mesh Ewald (PME) method was used to treat long-range electrostatic interactions with the switching distance of 10 Å. During the entire MD simulations, the coordinates of each system were saved for every 1 ps. Through the course of 10 ns simulation with the restrained backbone, lipid bilayer had established a stable contact with the protein for all the systems, and the levels of hydration of both the protein and the lipid had reached a steady state. The GPR109B system was the slowest to equilibrate, but even that had reached a reasonably equilibrated state for the medium (Figure S1A). After the release of protein backbone restraints, further equilibration was performed on all the systems for an additional 30 ns. The adjustments in both protein/lipid/water contacts were stabilized by around 15 ns (Figure S1B). Backbone RMSD reached a plateau on the same timescale as well, being in the region of 6 Å for the least stable receptor GPR109B (Figure S1C). Moreover, the transmembrane core of the receptors has been even more stable (RMSD deviation of the backbone stabilizing at the level of 4 Å over 30 ns), whereas the main dynamic regions were the extracellular domains connected by flexible linkers (Figure S1D).

Cluster analysis

To get the idea about the most dominant pose of niacin in the complex structure throughout the 100 ns MD trajectory, RMSD based cluster analysis was performed for niacin employing the CCPTRAJ module of AMBER16. ⁴⁴ Based on the RMSD of selected atoms, throughout the given trajectory, clusters of MD frames were constructed for similar conformations, and one frame of a cluster was considered as cluster representative structure of that particular cluster. For cluster analysis, DBScan algorithm which needs two parameters for clustering namely; distance cutoff to make a cluster ie, epsilon, and the minimum number of points to form a cluster ie, minpoint are considered. The obtained cluster representative structures were further compared with the respective energy minimized structure.

MM/PBSA analysis

In order to check the thermodynamic stability of all the complexes of GPR109A with niacin, MM/PBSA calculations employing a combination of programs CaFE, ⁴⁵ NAMD, ²⁶ VMD ⁴² and APBS ⁴⁶ was used. In MM/PBSA approach, molecular mechanics is combined with Poisson-Boltzmann and surface area continuum solvation to calculate the binding free energy of the protein-ligand complex. It is a post-processing method of MD simulations trajectory of receptor-ligand complex. It is also to be noted that MM/PBSA approach rapidly estimates the “calculated” binding free energy of protein-ligand complexes which can be well correlated with experiments but is not the absolute binding free energy. Since there are no major conformational changes in the receptors during MD simulations, a single trajectory method to calculate the binding free energy was preferred. For the calculation of binding free energy, a total of 5000 frames from the last 50 ns of 100 ns trajectory at an interval of 10 ns was selected. In the CaFE program, the gas phase energies (electrostatic and vdW) were calculated by NAMD, the polar solvation free energy was calculated by the PB equation implemented in APBS, and the non-polar solvation energy was estimated by a linear relation to the SASA (solvent accessible surface area).

Steered MD simulations

Steered Molecular dynamics (SMD) simulations method was also employed to explore the force profiles for moving the ligand into and out of the channel’s binding crevice. SMD simulations were performed using NAMD ²⁶ for the wild-type and the mutants of GPR109A, as well as for the GPR109B receptor. The starting conformation was set manually, with niacin positioned in random orientation at the outmost level of the extracellular vestibule—about 5 Å higher than the binding region predicted by automated docking. The starting position of the ligand was identical for all the systems. The steering harmonically restrained only z coordinate (normal to the membrane plane) of the center of mass of the ligand forcing it to change linearly with a speed of 1 Å / ns—first, toward the level of midplane of the membrane (10 Å over 10 ns) and then back outwards through the whole extracellular vestibule of the channel (20 Å, 20 ns). The restraining harmonic constant was set to 1 kcal/mol/Å. The ligand was free to change its conformation and move laterally in any direction in the pore crevice. The rest of the simulation parameters were the same as described above for unrestrained MD runs. The values for the applied force and ligand position were recorded every 1 ps and averaged with a running frame of 100 ps.

Results and Discussion

Model selection

Among the four models predicted from GPCR-I-TASSER (four each for GPR109A and GPR109B), the top most predicted structure based on the confidence score (C-score of −1.04 for GPR109A and −0.91 for GPR109B) was selected. The quality of GPCR-I-TASSER generated models was assessed by ERRAT (84.536 for GPR109A and 87.243 for GPR109B) and were scored high when compared to all the web server predicted models (HHPred, Phyre2, and SWISS-MODEL including the Modeler generated rhodopsin-based model). Similarly, the Ramachandran plot suggested that amino acid residues located in the favored regions were 94.8% and 95.5%, and other amino acid residues in the additional region were 5.2% and 4.5% for GPR109A and GPR109B, respectively, with no residues in the disallowed regions. All the other web server models and Modeler generated rhodopsin-based model have scored comparatively low percentage of amino acid residues located in the favored region and small percentage of amino acid residues located in disallowed regions. Thus, the current computational assessment suggested that GPCR-I-TASSER generated models demonstrated excellent quality over the other web server and Modeler generated model.

Molecular docking reveals distinct binding sites for niacin between GPR109A and GPR109B

The AutoDock program was implemented with the reference GPR109(A/B) structure with a grid box centered on the key amino acids involved in the niacin binding to GPR109A identified from Tunaru et al ²³ mutational studies. Notably, AutoDock predicted significant differences in the docking position of niacin with GPR109A and GPR109B along with the cluster results. The analysis for the GPR109A docking predicted 95% of the clusters in the docking results near the residues K166, S179, S178 and L176 with a predicted binding energy of −5.98 kcal/mol. Comparatively, in the GPR109B structure, 10 different clusters were formed with the top cluster predicted 37% near to the R111 with a predicted binding energy of −4.92 kcal/mol. Although the variation in the binding energies of niacin to GPR109A and GPR109B were favorable, the calculated values were not the absolute binding energies, but were only representative measures of how favorably niacin binds to the protein. Furthermore, the best docked poses from the highest scoring cluster for both GPR109A and GPR109B were selected based on the minimum binding energy (kcal/mol), number of stabilizing interactions such as hydrogen bonds and other weak interactions, docked scores, and cluster RMSD values. The docked poses were visualized using the VMD software package. ⁴² Finally, the most energetically favorable conformation of niacin was selected as an input for MD simulations.

Molecular dynamic simulations of GPR109(A/B) with niacin

In the current work, MD simulations were carried out to study GPR109A model with niacin positioned based on AutoDock predictions (for simplification we named it pose 0) in the water-lipid bilayer box and compared our results with previous studies of water-vacuum-water box system. ²³ Furthermore, to get the dynamic trends in niacin binding to GPR109A in the presence of an explicit medium, RMSD-based structural clustering of niacin/receptor complex considering 10000 frames taken at an interval of 10 ps from the whole 100 ns MD trajectory was performed. To notice the structural differences during the MD run, comparison between the cluster analysis-based representative structure to the starting energy minimized structure was considered. In the energy minimized structure of the GPR109A-niacin complex, the niacin carboxyl group makes a salt bridge with the amino acid K166 (Trans Membrane Helix 4, TMH4), a hydrogen bond with S179 (ECL2) backbone nitrogen and the pyridine nitrogen at a hydrogen bond (3.2 Å) distance with the W93 (ECL1) side chain (Figure 1A). In contrast, in the representative structure of the most populated binding cluster in MD simulations started with exactly same docked energy minimized structure, niacin changes position and makes 2 salt bridges—with R111 (TMH3) and R251 (TMH6) (Figure 1B). It seems that, it is the presence of strong positively charged guanidine group of these amino acid residues in the binding pocket that caused niacin to leave the starting location (at the entry route of the binding pocket in energy minimized structure) and move inside the pocket during MD simulations and form these 2 salt bridges (Figure 1A and B). For every simulation run, the H-bond interactions of niacin with GPR109A for 10 000 frames are taken at an interval of 10 ps from 100 ns MD trajectory. Among all the surrounding residues, niacin makes H-bond (Table 1) for maximum occupancy time with the following ones: R111 (TMH3) (24%) and R251 (TMH6) (38%). According to the H-bond analysis, niacin spends about one-third of the whole simulations time in position enabling simultaneous contact with both residues. In the representative structure of the highly dominant cluster (Figure S3A), the formation of these H-bonds further supports the presence of salt bridges between niacin and R111 (TMH3) and R251 (TMH6). Moreover, conformation of a membrane protein such as GPR109A, crucially depends on the surrounding medium. Even if the protein backbone would be restrained, a long-range electrostatic interaction with the medium might affect both the conformation of the protein side chains and the behavior of the polar ligand inside.

Figure 1.

LigPlot representation of residues interacting with niacin in GPR109A-niacin complex, pose-0, (A) energy minimized structure and (B) most populated binding cluster representative structure.

Table 1.

Hydrogen bond (H-bond) analysis of niacin with GPR109A during 100 ns MD simulation.

GPR109A amino acid residues	H-bond occupancy time (%) in GPR109A-niacin complexes
	pose-0	pose-1	pose-2	pose-3	pose-4
R111	24	48	34	21	–
K166	–	–	13	19	–
C177	–	–	–	–	45
S178	–	–	–	–	43
S179	–	–	–	38	15
R251	38	45	22	–	–

Abbreviation: MD, molecular dynamics.

According to previous MD studies on GPR109A-niacin complex ²³ (performed in water-vacuum-water box), the niacin pyridine ring embeds between sidechains of amino acids at the junction of TMH2/ECL1: W91, and residues of TMH7: F276 and Y284. In their system, the ECL2 residue: S178, makes H-bond with pyridine nitrogen of niacin. The TMH2 residue: N86, restrains the orientation of W91 (residue at the junction of TMH2/ECL1) as a results of sidechain atoms of both residues are involved in H-bonding. Another ECL2 residue: F180 restrains the orientation of F276 (TMH7) sidechain through aromatic interaction. In our MD simulations, niacin forms salt bridges with R251 (TMH6) in addition to R111 (TMH3) and the pyridine ring is directed on the opposite side of the crevice where the sidechain of amino acid residue W91 (at the junction of TMH2/ECL1), and the TMH7 residues: F276 and Y284 are oriented (Figure S2A). Moreover, among these 3, the nearest residue F276 (TMH7) has minimum distance between the F276 (TMH7) sidechain and niacin about 5.2 Å and there is no any direct interaction with niacin. Thus, in our simulated system, it is not possible to embed the pyridine ring at the junction of sidechains of residues W91 (at the junction of TMH2/ECL1), residues of TMH7: F276 and Y284. The orientation of F180 (ECL2) is outward of the binding pocket and is restricted by N-terminal residue; F10 stabilized the position of F180 through aromatic interactions, and because of that, the F180 (ECL2) sidechain is quite far from F276 (TMH7) and is not likely to restrict F276 dynamics (Figure S2A). As a result, S178 (ECL2) is unable to make interactions with the pyridine nitrogen atom of the ligand. There is an H-bond between the N86 (TMH2) and W91 (residue at the junction of TMH2/ECL1) but placed only between the backbone atoms instead of sidechain atoms (that has been observed in previous study ²³ ). The strong restriction of W91 (residue at the junction of TMH2/ECL1) sidechain orientation by N86 (TMH2) is also highly unlikely (Figure S2A). In addition, since the sidechain orientation of F180 (ECL2) is outward, it is unable to make aromatic interaction with F276 (TMH7), allowing niacin in occupying more space to interact with R111 (TMH3) and R251 (TMH6) and not allowing the pyridine ring of niacin to embed between the sidechains of the residues; at the junction of TMH2/ECL1: W91, residues of TMH7: F276 and Y284.

As described above, during MD simulations of the receptor with niacin that was initially placed at the docking-based-pose-0 on the entry route of GPRA109 binding pocket, niacin subsequently moves inside the receptor crevice and was found to bind to R111 (TMH3) and R251 (TMH6). Interestingly, it has also been shown by previous study ²³ that R111 (TMH3) is an important residue for niacin binding. As the next step, and whether R111 (TMH3) is the niacin binding residue to get a final docking pose, we ran MD simulations of 4 more GPR109A-niacin complexes (pose-1, pose-2, pose-3, and pose-4), in which niacin was placed randomly in the GPR109A binding crevice and the entry route at the start of different MD simulations run. Similar to pose-0, after performing 100 ns MD simulations of all the 4 complexes, RMSD-based was performed and structural clustering of niacin and H-bond analysis (Table 1) of 10000 frames taken at an interval of 10 ps from the whole 100 ns MD trajectory and compared with respective energy minimized structures and also compared the orientation of key residues in the MD simulated structure with previous study. ²³ In the first random complex, pose-1, niacin was placed at an H-bond distance of R111 (TMH3), but not R251 (TMH6). In the energy minimized starting structure of pose-1, niacin makes salt bridge only with R111 (TMH3) (Figure 2A), whereas in the cluster representative structure of the most dominant cluster, it makes salt bridges with both R111 (TMH3) and R251 (TMH6) (Figure 2B) similar to what was observed in complex pose-0. Moreover, according to H-bond analysis of the trajectory, residues that make H-bond with niacin for the maximum occupancy (48% and 45% throughout the MD trajectory) are R111 (TMH3) and R251 (TMH6) as shown in Table 1. This again supports the presence of salt bridge by these 2 residues with niacin in the cluster representative structure of the most dominant cluster. Since niacin was placed at a H-bond distance (2.9 Å) to R111 (TMH3), the guanidine positive charge group of R111 (TMH3) and R251 (TMH6) quickly attract the niacin carboxyl group to form a salt bridge during the MD simulations. Comparison of key residues and their sidechain orientation have shown that these side chain orientations are similar to the ones that are observed in simulations as the starting structure of complex pose-0 (Figure S2B). Another complex, pose-2, in which niacin was placed at start of the run at a distance of ~4.6 Å from both R111 (TMH3) and R251 (TMH6) (Figure 2C), niacin ends up moving closer and making H-bonds with these 2 residues with an occupancy time of 34% and 22% (Table 1), respectively. It also makes another H-bond with K166 (TMH4) with an occupancy time of 13%. In the cluster representative structure of the most dominant cluster (Figure 2D), there are 2 salt bridges of niacin with R111 (TMH3) and K166 (TMH4). This indicates that in this conformation, while the salt bridges with R111 (TMH3) and K166 (TMH4) coexist, the bridges of niacin to R111 (TMH3) and R251 (TMH6) rather alternate. Nevertheless, all 3 contribute to the stability of the ligand at this location. In the receptor-ligand complex, pose-2, due to interaction with K166 (TMH4) niacin pyridine ring is embedding between sidechains of R111 (TMH3), and residues of TMH7: Y284 and S276 (Figure S2C) and so is unable to reach anywhere near to W91 (residue at the junction of TMH2/ECL1). Furthermore, in complex, pose-3, niacin was placed at a distance of 3.1 Å, from W93 (ECL1) (Figure 3A). The H-bond analysis shows that there is no H-bond between niacin and residue at the junction of TMH2/ECL1; W91, during MD simulations run, but the ligand rather formed H-bonding with R111 (TMH3), K166 (TMH4), and S179 (ECL2) with the occupancy of 38%, 21%, and 19%, respectively (Table 1). In the cluster representative structure of the most dominant cluster, niacin is at an H-bond distance of S179 (ECL2) and K166 (TMH4) (Figure 3B) but there is no interaction with R111 (TMH3). Since niacin moves away from R111 (TMH3) after 45 ns and makes H-bond with S179 (ECL2) for most of MD simulations time, thus, in cluster representative structure there is no interaction with R111 (TMH3). In the complex of pose-3, sidechain orientation of all the key residues of binding pocket is similar to complex pose-0, and pyridine ring of niacin in the final frame is not embedded in-between the sidechain of key residues (Figure S2D). In the last complex, pose-4, niacin was placed near Y269 (Figure 3C), and it makes 2 H-bonds with ECL2 residues: C177 (backbone atom) and S178 (sidechain atom) with the occupancy time of 45% and 42% (Table 1), respectively. Also, in the cluster representative structure of the most dominant cluster, these 2 residues make H-bond with niacin (Figure 3D). Thus, the presence of these H-bonds in the cluster representative structure are complemented by H-bond analysis. Similar to other complexes, niacin does not embed in-between the sidechain of any key residues as suggested in the previous study ²³ (Figure S2E).

Figure 2.

Residues interacting with niacin in GPR109A-niacin complex. (A) and (C) are energy minimized structure and (B) and (D) are most populated binding cluster representative structure of pose-1 and pose-2, respectively.

Figure 3.

Residues interacting with niacin in GPR109A-niacin complex. (A) and (C) are energy minimized structure and (B) and (D) are most populated binding cluster representative structure of pose-3 and pose-4, respectively.

Overall, the niacin-GPR109A interaction analysis shows that niacin makes interaction with R111 (TMH3) in 4 out of 5 complexes; pose-0, pose-1, pose-2, and pose-3 with the different occupancy time of the whole 100 ns MD trajectory. In the complex, pose-0, when it docked at the entry route of the binding pocket during MD simulations it moved inside the pocket and make salt bridges with R111 (TMH3) and R251 (TMH6). When randomly placed at a H-bond distance and ~4.6 Å away from R111 (TMH3) in complexes, pose-1 and pose-2, respectively, niacin again makes salt bridges with R111 (TMH3) and R251 (TMH6). In complex, pose-3, where it was placed near to H-bond distance of W91 (residue at the junction of TMH2/ECL1), it makes H-bond with R111 (TMH3) in addition to other residues K166 (TMH4) and S179 (ECL2).

MM/PBSA analysis

After analyzing GPR109A-niacin interactions in all the complexes, MM/PBSA calculation was performed to calculate the binding free energy of niacin to GPR109A. The calculated binding free energies of all these complexes are mentioned in (Table 2). The calculated binding free energies of complexes, pose-0 and pose-1 are −6.5 and −9.5 kcal/mol, respectively. According to H-bond analysis, both the complexes have similar interactions (Table 1) during MD simulations which is easily related with similar calculated binding free energies. For another 2 complexes, pose-2 and pose-3, the calculated binding free energies are −16.9 and −18.3 kcal/mol, respectively. In addition to interaction with R111 (TMH3) and R251 (TMH6), niacin also interacts with K166 (TMH4) in the complex pose-2, which correlates well with the more energetic stability of pose-2 than pose-0 and pose-1. The complex pose-3, has only interactions with K166 (TMH4) and R251 (TMH6) but both the interactions present throughout the 50 ns trajectory consider for MM/PBSA calculation. Thus, pose-3 is energetically more stable like complex pose-2. The complex, pose-4, has binding free energy of −8.5 kcal/mol. In this complex, niacin makes 3 H-bonds with ECL2 residues: C177, S178, and S179. Though there are total of 3 interactions between niacin and GPR109A in pose-4, and all are H-bonds, it is not as much energetically stable complex as pose-2 and pose-3. According to the calculated binding free energies, all the complexes are energetically stable and niacin makes interactions with residues R111 (TMH3), K166 (TMH4), S178 (ECL2), and R251 (TMH6) highlighted in the previous study. ²³ Although the MM/PBSA method does not determine the absolute binding energy, but it only predicts calculated binding free energy based on the difference of free energies of complex, receptor, and ligand for given number of MD trajectory frames. The ligand binding affinities has been estimated by MM/PBSA and MM/GBSA are having correlation coefficient of r ²  = 0.0 to 0.9 compared to experiments. ⁴⁷ Thus, in GPR109A protein, residues R111 (TMH3), K166 (TMH4), S178 (ECL2), and R251 (TMH6) are involved in making interactions with niacin in different complexes and are responsible for these complexes to be energetically stable.

Table 2.

MM/PBSA calculation for GPR109A-niacin complexes.

MM-PBSA binding free energy (SD)
pose-0	pose-1	pose-2	pose-3	pose-4
−6.2 (5.1)	−9.4 (6.5)	−16.9 (5.7)	−18.3 (6.0)	−8.5 (4.8)

Abbreviations: MM, molecular mechanics; PBSA, Poisson-Boltzmann solvent accessible surface area.

SMD results

To estimate whether the mutations known to affect binding of niacin in GPR109A affect the ligand movement through the channel crevice comparing to WT (GPR109A) protein, SMD was employed. Similarly, this approach was also used to assess niacin translocation through the homolog GPR109B that does not bind niacin. With this approach, the ligand was placed at the level of the external rim of GPR109A vestibule and slowly steered it with constant velocity (1 Å / ns), first down to the bottom of the upper crevice and then back out of the channel. The position of the ligand’s center of mass was only restrained along the pore axis, while allowing it to change conformation and choose any favorable lateral position. The recorded force along the steering direction reflects the effort needed to propagate the ligand.

While the exact shape of the force profile varies among the mutants (Figure S3), the results consistently show that all the mutants require much higher force to move the ligand through the GPR109A crevice. For steering from the vestibule down to the bottom of the crevice, the mutants needed an extra force ranging from 20% to 199% more (Figure 4A). For the subsequent translocation out of the crevice, the mutants required 482% to 711% more force comparing to WT (Figure 4B). Since the force was required both on the way into and out of the protein crevice, it suggests that the nature of impediment for ligand motion was of sterical rather than electrostatic nature. Notably, GPR109B homolog is known not to bind niacin and in SMD simulations it also showed 84% higher force needed to translocate the ligand into the crevice (Figure 4A) and 421% higher force to get it out. SMD results for GPR109B, similar to GPR109A mutants, are consistent with sterical obstacles for the ligand.

Figure 4.

The force required in steered molecular dynamics to reposition niacin in the binding site of GPR109A/B and different mutants of GPR109A. (A) Steering the ligand down to the bottom of the crevice (B) Steering the ligand (translocation) out of the crevice.

Conclusion

The present study provides new structural and energetic insights into the interaction of niacin with GPR109A/B in comparison with the previously published work from Tunaru et al. ²³ Even though 1 ns MD simulations of GPR109A-niacin complex in water-vacuum-water box provides some crucial information, use of an explicit lipid in MD simulations of a membrane receptor affects the system both by making direct interactions with the protein and by the remote effect of the electrostatic field of the membrane on the position and conformation of the ligand and receptor side chains. In the current 100 ns MD simulations, TMH3 residue R111 is involved in making interactions with niacin in 4 out of 5 complexes (pose-0, pose-1, pose-2, and pose-3). While, in the last complex, pose-4, niacin makes interactions with ECL2 residues, C177, S178, and S179, present at the entry route of GPR109A which was not identified by Tunaru et al. ²³ Furthermore, our calculated average force for niacin translocation through the GPR109A vestibule in wild type and mutant (N86Y, W91 S, S178I, and triple mutant of these 3 residues) complexes through steered MD simulations shows an increase in the force in mutant(s) compared to wild type, which might result from sub-optimal location and/or conformation of the ligand due to the mutations. The residues at the junction of TMH2/ECL2: N86, W91 (ECL2), and ECL3 residues: C177, S178, and S179, and TMH3 residue R111 are important for niacin entry and binding to GPR109A, and might orchestrate the proper conformation dynamics of the ligand and enable a favorable binding.

Supplemental Material

sj-docx-1-bbi-10.1177_11779322211056122 – Supplemental material for Structural and Energetic Insights Into the Interaction of Niacin With the GPR109A Receptor

Supplemental material, sj-docx-1-bbi-10.1177_11779322211056122 for Structural and Energetic Insights Into the Interaction of Niacin With the GPR109A Receptor by Kiran Kumar Adepu, Sangita Kachhap, Andriy Anishkin and Sree V Chintapalli in Bioinformatics and Biology Insights