MATERIALS AND METHODS

Study area: This study was conducted at the Pharmacy School, Shiraz University of Medical Sciences, Iran. This project has been started in January, 2019 and ended in May, 2019.

Research procedure: A core i7 laptop on windows 8.1 was utilized for the preparation chemical structure of the compounds. Simulation processes were run by a 24 core computational server on Linux Ubuntu12. First, for the purpose of the homology modeling process, the UniProt database (www.uniprot.org) was used for drawing out the sequence of H₂ receptor (PDB code:P25021) as a FASTA template¹⁴. I-TASSER search engine (Iterative Threading ASSEmbly Refinement http://zhanglab.ccmb. med.umich.edu/ I-TASSER) determined the acquired sequence to recognize the templates from the Protein Data Bank.

The highest C-score obtained from the I-TASSER server was chosen as the best PDB structure to run the MD simulation process¹⁵.

In this study, 11 servers that mentioned below were utilized to predict transmembrane helices or alignment the model that possesses an accurate topology inside lipid bilayer:

In the next step, to get started molecular dynamics simulation, GROMOS96 53A6 force field was applied as performed in Gromacs 4.5.5 that describes the lipids derived by Berger lipid parameters²¹ and situates the receptor inside a lipid bilayer. A 128 DPPC (dipalmitoyl-phosphatidylcholine) and some commands are needed for the simulation. Also, VMD software was employed to align the intended procedure²².

A method called the Inflate GRO develops the system totally for the elimination of extra lipid residues and incorporation of the protein²³. To set the protein in the bilayer membrane, the mentioned method was employed in the role of an algorithm. In the next step, the most appropriate area per lipid for DPPC systems was prepared by repetitive runs of shrinking and minimization²⁴.

Area per lipid density during all steps of shrinking and minimization was calculated by employing GridMAT-MD_v2.0 Perl script²⁵. Following this, water and ions were added to the simulating system in order to deter permeation of water molecules inside the hydrophobic parts of the lipid membrane. In this step, the van der Waals radius of carbon atom was set to 0.375Å. A concentration of 0.15M NaCl was added on the system to simulate the physiological environment. After the Energy minimization step, the protein backbone was restrained. The simulation followed by 2 equilibration experiments including NVT and NPT ensemble. In these steps the system was subjected to a minimization step. While the Nose-Hover algorithm was used as an accurate thermostat in the NPT ensemble, a modified Berendsen was used in the NVT v-scale. Periodic Boundary Condition (PBC) was prepared during the simulation by using Particle Mesh Ewald (PME)long-range electrostatics. A 50 ns Molecular Dynamics Simulation was executed on the system to give enough time for conformational changes inside the receptor. In order to begin the docking process, the receptor should be sampled. For a sampling of several receptor conformations, TCL scripting was employed. Therefore, 200 PDB structures were adapted from the output file of MD trajectories covering all simulation time. Via assigning Gasteiger partial charges using MGLTOOLs 1.5.6 the obtained PDB files were converted to pdbqt files²⁶. In the next step, to prepare ligand, a collection of 28 molecule ligands were resumed from the CHEMBL database²⁷. In the next step, by using Open babel 2.3.2 the structures were converted to mol2 format. By the way, the Ki values were retrieved and saved into PKi (-logki). MGLTOOLS 1.5.6, for production 28 pdbqt files, was used to add the Gasteiger partial charges and merge non-polar hydrogen atoms.

Afterward, for the determination of the binding site, the sequence of H₂ receptor was subjected to RaptorX(http:// raptorx.uchicago.edu)²⁸.

Finally, in order to carry out the Molecular Docking process coordinate and size of the grid box must have been determined in advance. Based on the former literature, Coordinates of C_α for LYS121 were assumed as the grid box center. The size of the grid box computed according to the equation below:

Size x; y; z = 2×LDA = 30Å

LAD marks the spot that atomic distance is the largest in all data collections. Afterward, the AutoDock Vina 1.1.2 software²⁵ was utilized for 200*33 = 6600 cross-docking simulations²⁹. Exhaustiveness value was set to 100 to execute impressive docking simulations in AutoDock Vina for all outputs.

Statistical analysis: The Pearson correlation coefficient (R) as the statistical tool for analyzing the correlation of docking results, was computed between every single obtained docking energy value and the pKi values for each frame. Receptor-ligand interactions were shown based on some pattern gotten by protein-ligand interaction profiler (PLIP) (https://projects. biotec.tu-dresden.de/plip-web/plip/)³⁰.

RESULTS

Five templates were used during the modeling: 3MY9_A, 4U16_B, 32PQ_A, 2Y01_A and 2Y00_B.

The top 5 I-TASSER models based on their C-score is represented in Table 1. The model with the highest C-score was considered as the best one.

Table 1:	Top 5 I-TASSER models based on their C-score

Table 2:	Eleven methods were applied in order to predict TM regions of H2-receptor

Fig. 1:	Ramachandran plot for the protein was obtained from PROCHECK

The TM-score for the bestmodel is higher than 0.5 which means it has an acceptable topology for further modeling studies.

Ramachandran plot for the protein was obtained from PROCHECK Server (Fig. 1)³¹. The function of this plot is to verify the obtained model. More than 98% of the residues are placed in the favorable part of the plot, which means that conformational features of the modeled protein are tantamount to the native proteins.

11 methods were applied in order to predict TM regions of H₂-receptor. As presented in Table 2, the topology of membrane proteins has been predicted precisely and efficiently.

After obtaining an acceptable model for the 3D structure of the protein, the next step was to perform molecular dynamics (MD) simulations. The purpose of these simulations was to recognize the various conformations of the protein in the physiological conditions. A 50 ns MD simulation was conducted on the whole system to find out an equilibrate state for the receptor based on the RMSD variation of C_α and energy plot. The result of this step is depicted in Fig. 2a-b. The steady-state reaching time was 44 ns after the beginning of the process. Additionally, the energy diagram shows an equilibration of energy during the simulation.

Fig. 2(a-b):	A 50 ns MD simulation was conducted on the whole system to find out an equilibrate state for the receptor based on (a) RMSD variation of Cα and (b) Energy plot

To find out the most fluctuating parts of the receptor, a heat-map plot was delineated using the VMD software. This plot is shown in Fig. 3. Based on the heat-map, two extracellular domains of the protein showed a high level of fluctuation. These two domains are the residues 1-10 and 161-180. Additionally, a cytoplasmic domain revealed a high degree of fluctuation. This was the domain containing the residues 290-311.

For the docking simulations, it was necessary to find out the central residue of the grid box. Based on the RaptorX method and using of MOE, LYS121 was selected as the center. In order to opt the best frame for the docking simulation, the Pearson correlation of all docking frames with PKi values was attained. The most correlation was calculated for the frame 126 (|r| = 0.9) (Fig. 4).

To provide an outline for the binding mode of the ligands in the active site of H₂-receptor, the interaction maps of the receptor for all the compounds were obtained. The results of all ligand-receptor interactions have been represented in Table 3.

CHEMBL474991 has the best docking score (-10.1 kcal moL^–1). Figure 5 depicts the interactions between CHEMBL474991 and the active site of the target.

Mobility of binding site residues was also calculated by mapping the RMSD values of the relevant residues during the simulation. As illustrated in Fig. 6 variations of the binding site was fixed at the last steps of simulation with an RMSD average of 2Å.

DISCUSSION

Obtained results of this study shows a reliable template for human histamine-2 receptor. This 3D template has been verified by molecular dynamics simulation and cross-docking evaluation.

Fig. 3:	To find out the most fluctuating parts of the receptor, a heat-map plot was calculated using the VMD software

The original technique used in this study is homology modeling. Nowadays, homology modeling plays a major role in all levels of computational biology: from genes to macromolecules³². It can be called the most precise method for structural prediction⁶. It is also a beneficial tools for more comprehensive studies such as pathway detections³³. Homology modeling for proteins is not confined to finding the structure of receptors. it make us able to design more confident peptide drugs such as antibodies³⁴. Structure-based drug design requisite is a comprehensive library of receptors in human body. However, our current treasure of proteins is not as exhaustive as needed. Homology modeling can play a benevolent role in empowering this collection. As some previous studies have proven, homology modeling alongside molecular dynamics and docking simulations can be an inexpensive way to shed the light on the unknown structures^5,35. The output of this study can be seemed as an inchoate achievement.

Table 3:	All ligand-receptor interaction

Fig. 4:	In order to opt the best frame for the docking simulation, Pearson correlation of all docking frames with Pki values were attained

Fig. 5:	Interactions between CHEMBL474991 and the active site of the target

Fig. 6:	Mobility of binding site residues was also calculated by mapping the RMSD values of the relevant residues during simulation

This template can be recruited by further computational studies on H₂ receptor to become a definitely approved model.

CONCLUSION

The validity of obtained template for H₂ receptor was testified by means of several in silico techniques. A high level of similarity between the drawn model and the predicted structure was established by a Ramachandran plot. Cross-docking studies for recognized H₂ receptor antagonists were our second mechanism to justify the obtained model. These studies were performed on the frame 126 which had extracted by molecular dynamics simulations. According to the Pearson correlation test, one frame was accepted. The coefficient for this frame was 0.9 that represents a high degree of reliability for the results. Furthermore, our computational method could perfectly mimic the human H₂ receptor at the antagonist state.

SIGNIFICANCE STATEMENT

This research plays a crucial role in the prospective drug design studies in order to recognize new ligands. Due to the unavailability of the protein structure of H₂ receptor, the findings of this study pave the way for prospective computational studies on this receptor. The achieved structure based on this study enlarges the current library of human proteins, which is a strategic database for biological researches. Additionally, this study affirms the viability of homology modeling as an accessible, affordable and valuable technique in biological studies. Further unavailable protein structure can be designed by means of the methods utilized in the current study.

ACKNOWLEDGMENT

This research project was conducted from January-May, 2019. This Research Project was fully sponsored by Shiraz University of Medical Sciences with grant number 8915574058.