Computational investigation of the impact of potential AT2R polymorphism on small molecule binding

Abstract For more than a century, the renin-angiotensin system (RAS) has been acknowledged for playing a crucial part in the physiological control of arterial pressure, as well as sodium and fluid balance. It is now generally acknowledged that one of the receptor of RAS system i.e. angiotensin type 2 receptor (AT2R) functions as a repair system during pathophysiologic circumstances and performs a significant protective role. Efforts have been made previously to design suitable agonist and antagonist molecules to potentially modulate AT2R. One of the agonists and antagonists, named C21 and EMA401, has been studied in a number of pathological conditions. Additionally, a wide panel of single nucleotide polymorphisms (SNPs) has been reported for AT2R, which might potentially affect the efficacy of these molecules. Therefore, computational investigations have been carried out to analyze all the SNPs (1151) reported in NCBI to find potential SNPs affecting the active site of AT2R, as this domain is still unexplored. Structures of these polymorphic forms were modeled, and in silico drug interaction studies with C21 and EMA401 were carried out. The two mutants (rs868939201 and rs1042852794) that significantly affect the binding affinity as that of the wild type were subjected to molecular dynamics simulations. Our analysis of native and mutant AT2R and their complexes with C21 and EMA401 indicated that the occurrence of these mutations affects the conformation of the protein and has affected the binding of these ligand molecules. The study’s findings will aid in the development of better, more versatile medications in the near future, and also in vitro and in vivo studies might be planned in accordance with recent findings. Communicated by Ramaswamy H. Sarma


Introduction
It is widely believed that hypertension is the most serious public health issue in the world due to its enormous prevalence (Singh et al., 2017).High blood pressure causes around 7.5 million deaths worldwide annually.It is predicted that up to 2025, there will be an increase in 1.56 billion adults with hypertension (Tabrizi et al., 2016).It is of two types that are primary and secondary.About 90 to 95% of cases are of primary hypertension, which means having high blood pressure for which no medical cause can be found.The remaining 5 to 10% of patients are from secondary hypertension, which is caused by other conditions which affect the kidneys, arteries, heart, or endocrine system (Carretero & Oparil, 2000).
The renin-angiotensin system (RAS) is considered the crucial regulatory system of the body that has been implicated in the pathogenesis of hypertension, renal disease, and cardiovascular disease.Since the discovery of renin as a pressor substance in 1898, the RAS system has been broadly studied because this system remains a primary contributor in the development and maintenance of hypertension (Lavoie & Sigmund, 2003).The RAS system genes mainly consist of renin, angiotensinogen (AGT), angiotensin-converting enzyme (ACE), angiotensin-converting enzyme 2 (ACE2), angiotensin type II receptor (AT2R) and angiotensin type I receptor (AT1R) (Barnas et al., 1997).
The AT2R receptors in fetal tissues are predominantly expressed, and their expression decreases in adult tissues, mainly restricted to a few organs (Levy, 2005).The gene that codes for the AT2R is present on the human Xq22-q2 chromosome, and this gene consists of two introns and three exons, and the region responsible for coding protein is present only on exon 3 (Zhou et al., 2007).In adult humans, the AT2R is primarily present in the uterus, brain, adrenal medulla, ovaries, kidney, and heart.The occurrence of this receptor in these organs indicates its involvement with the physiological implication in these organs (Matavelli & Siragy, 2015).The primary functions of the AT2R are to counteract the tasks performed by the AT1R by obstructing differentiation and cell proliferation, enhancing vasodilatation, and reducing oxidative stress and inflammation (Fatima et al., 2021).
In recent times, the RAS was projected to be implicated in many diseases as a probable drug target by blocking the system, either with AT1R blockers or by targeting angiotensinconverting enzymes with the angiotensin-converting enzyme inhibitors (ACEIs) (Gard et al., 1999).The other component of RAS, i.e., AT2R, is not exploited much as a drug target, but few potential agonist and antagonist molecules were reported for AT2R (Wan et al., 2004).Huge polymorphism in AT2R changes an amino acid in protein, which is expected to alter the function of AT2R and interactions with drugs to target the receptor for curing many diseases (Sharma et al., 2022).These SNPs exist in exons, introns, promoters, and UTR regions.In the case of AT1R, a previous study suggests that the presence of single nucleotide polymorphism (SNP) affects the efficacy of drugs (Arsenault et al., 2010) while no such reports were there for AT2R.
Therefore, computational methods in the current study were used to analyze all the SNPs of AT2R, initiated by finding and creating potential SNP mutants through homology modeling, and investigate their effects on the binding affinities of potential agonist and antagonist molecules for AT2R through molecular docking.The dynamic behavior of selected complexes was also analyzed using classical molecular dynamics simulations.The outcome of the study will help to study binding poses and design better drugs in the near future, which will also be effective in the case of SNP mutations.Furthermore, in vitro and in vivo studies can be designed according to current results to develop better therapeutics with broader coverage and subsequently validate the current approach.

Material and methods
The binding cavity's amino acid residues were examined for the presence of SNPs, and subsequently, SNPs were located using the stepwise method (Sharma et al., 2021).All the information was gathered from the National Center for Biotechnology Information to identify all SNPs in the human AT2R gene (https://www.ncbi.nlm.nih.gov/).For homology modeling, the crystal structure of AT2R (H.Zhang et al., 2017) was used as template.

AT2R SNPs analysis and mutant preparation
The human AT2R SNPs were downloaded from the NCBI SNP database.All SNP mutations were created manually on the mRNA sequence of human AT2R obtained from NCBI.Then all these manually mutated mRNAs were translated with the help of ExPASy TRANSLATION TOOL (https://web.expasy.org/translate/).The mutated self-transcribed proteins of AT2R were analyzed for the change of amino acids and their location in protein structure.Hence, the amino acid change at appropriate positions due to SNPs was identified and created on the wild-type PDB fasta file of AT2R (https://www.rcsb.org/structure/5UNF).Subsequently, all mutated sequences were saved in different text files.Furthermore, homology modeling was performed to generate the 3D structures for these mutants.

Homology modelling
The 3D structures of mutants of AT2R were created using homology modeling using the crystal structure of human AT2R (PDB ID: 5UNF) (H.Zhang et al., 2017).Template selection was the first step for predicting a structure using homology modeling.The well-defined structure for AT2R is available; hence the structure was used as a template for homology modeling.Modeller program was used to perform homology modeling.The template and query sequence were used as input in Modeller 9.14 (Kashyap et al., 2017).For alignment and modeling of mutant structures, python scripts were used.In the end, optimized 3D structures for all selected mutants were generated.

Ligand preparations
The newly proposed agonist molecule C21 and antagonist EMA401 were used as ligand molecules for AT2R.3D structure for both these ligands were downloaded from pubchem database (https://pubchem.ncbi.nlm.nih.gov/) in sdf format.These downloaded sdf files were optimized for the purpose of molecular docking and converted to pdb format using Open Babel and furthermore converted to pdbqt format by using Autodock tools 1.5.6.

Molecular docking analysis of AT2R
The agonist and antagonist, C21 and EMA401, respectively, were docked with wild type and mutated AT2R to check the effect of mutation on their binding affinities using autodock vina (Seeliger & De Groot, 2010;Sharma et al., 2021).The AT2R contains two chains, A and B, and both are identical; hence for simplification, the B chain was deleted before docking of AT2R.The active site information was obtained from the data provided in the protein's crystal structure (H.Zhang et al., 2017).The information from the crystal structure was used to set up a grid box to surround the area of interest (active site).The grid was positioned at 75.966, À 3.22, and 38.853 Å as the X, Y, and Z center grid values having a size of 40 � 40 � 40 Å along x, y, and z dimensions.All the parameters like binding energy, hydrogen bond interactions, and the orientation of the docked compounds within the active site were visualized.The criteria used for elucidating the inhibitory potential of ligands were the binding energies of ligands with the receptor.Lesser the binding energy of the ligand with the receptor more affinity is shown by the ligand for that particular receptor.After the molecular docking study had been done for all the mutants and wild-type AT2R, the Discovery studio visualizer was extensively used to study and visualize the resultant receptor-ligand complexes.

Molecular dynamics simulation
The stability of the protein-drug complex was evaluated by performing molecular dynamics (MD) simulation using Gromacs 2020.4 installed on an Intel Xenon workstation-E3-1245-8C, 3.50 GHz processor with 28 GB RAM.The workstation was powered by a NVIDIA Quadro P5000 GPU card (AlAjmi et al., 2018;Jabir et al., 2021).The protein topology was generated by the pdb2gmx command using CHARMM36-all atom forcefield, and TIP3P water molecules.The drug topology was generated by the CHARMM General Force Field (CGenFF) program (Vanommeslaeghe et al., 2010).The simulation was performed inside a dodecahedron box with the docked protein-drug complex positioned at the center and at least 1.0 nm from the edges.The box was solvated with TIP3P water molecules and neutralized by adding proper Na þ or Cl-ions.A salt (NaCl) of 150 mM concentration was added to mimic the physiological conditions.The energy of the system was minimized with a maximum of 50,000 steps of steepest descent minimization.The solvent and ions were equilibrated in different restrained phases.The isothermal-isochroic NVT ensemble was performed at 300 K temperature, while the isothermal-isobaric NPT ensemble was conducted at 1.0 bar pressure.The temperature and pressure of the system was maintained by the Berendsen thermostat and Parrinello-Rahman barostat, respectively (Parrinello & Rahman, 1981).Finally, a 100 ns production run on the equilibrated system.A time-step of 2 fs was fixed using Leap-frog integrator.LINCS algorithm was used as a constraint algorithm for NVT, NPT, and production runs.The short-range van der Waals' cutoff was set at 1.2 nm.The simulation results for molecular dynamics were examined for characteristics including root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), intermolecular hydrogen bonds, and analysis of interacting residues.All the experiments were performed independently in triplicates and the results are reported as mean ± standard errors.

Free energy calculation by prime/MM-GBSA
The free energy of interaction between the target proteins and ligands was evaluated using the Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) approach.Here, the trajectory in the last 10 ns of equilibrated MD simulations was employed to compute free energy.First, the docked complex was locally optimized using the Molecular Mechanics (MM) approach.Second, the OPLS3e force field along with the Generalized Born Surface Area (GBSA) continuum solvent was used to minimize their energies.The binding free energy was calculated using the following equations: where, DG Complex , DG Protein , and DG Ligand are the respective values of minimized energies of protein-ligand complex, protein and ligand.

Results
In the present study, all known SNPs of AT2R (1151) were analyzed to find out the potential SNPs affecting AT2R protein and subsequent interactions with agonist molecule C21 and antagonist EMA401 were analyzed.Additionally, to further strengthen our claims, molecular dynamics simulations were performed to assess the dynamics behavior of mutated AT2R upon ligand binding.

Screening of AT2R mutants for docking studies
In the NCBI SNP database, 1151 hits appeared for AT2R; out of 1151 hits, 187 SNPs created missense mutations.All these mutations were created on AT2R mRNA (NM_000686.5), and subsequently, all were self-translated.After this, self-transcribed proteins were compared with the fasta sequence of AT2R, and after that, all located mutations were made on the (5UNF) AT2R PDB fasta sequence to create SNP mutated proteins.In the end, only 151 were present in exons that were taken for further investigations.And lastly, 8 SNP variants were screened, which potentially affected the amino acid residues of the active site of AT2R, and three-dimensional structures for all these variants of AT2R were constructed using homology modeling (Table 1).

Homology modeling for structure prediction of AT1R and AT2R mutants
The three-dimensional structure of AT2R is available; hence it was used as a template for homology modeling.A Modeller program was used to perform homology modeling.The validation of all mutated models was done by Ramachandran plot and PROCHECK, which depicts that most amino acid residues were present in the favored region for all mutated modeled structures of AT2R (Table S1).

Molecular docking of AT2R and AT2R mutants with agonist and antagonist
Molecular docking of agonist molecule C21 and antagonist EMA401 into the active site of all eight mutated AT2R and normal AT2R was successful based on the formation of the complex of AT2R with these agonist and antagonist molecules.The binding energy, hydrogen bond interactions, the orientation of the docked compound, bond length, and active site residues were visualized.The binding affinity of different mutated variants of AT2R decreased compared to the normal AT2R for C21.Out of eight mutated SNP variants of AT2R, one of the mutated variants rs868939201 (A290G), having a change in amino acid from Methionine to Isoleucine, shows the least binding affinity of À 8.1 kcal/mol and the interacting residue forming the hydrogen bonding is Thr254 (Table 2) (Figure 1a).Other mutated SNP variant rs782087304 (T291A) having a change in amino acid lysine to asparagines also showed less binding affinity of À 8.2 kcal/mol and the interacting residue which forming the hydrogen bonding is Asn291 (Table 2).There were two other variants rs782235260, rs782685066(A) showing the binding affinity of À 8.3 kcal/mol which is also less as compare to normal AT2R with C21 which is an agonist of AT2R.Other variants also showed lesser binding affinity than that of the normal AT2R.The mutated variant rs1042852794 exhibited a binding affinity of À 8.9 kcal/mol, the highest among the selected AT2R mutants, and interacted strongly with the Tyr 127 (Figure 1b).
In case of antagonist molecule of AT2R i.e.EMA401 most of the mutated variants of AT2R shows increase in binding affinity as compare to the normal AT2R (Table 3).One of the SNP variant rs868939201 having a change in amino acid Methionine to isoleucine shows increased binding affinity of À 11.3 kcal/mol and the interacting residues forming the hydrogen bonding are Lys291, Arg258 (Figure 1c).Another variant rs781868698 having change in amino acid tryptophan to serine shows binding affinity of À 10.9 kcal/mol and strongly interacting residue is Lys291 (Table 3).Out of eight variants there was one variant which shows decrease in binding affinity as compared to normal AT2R.The SNP variant rs1042852794 having change in amino acid methionine to isoleucine shows binding affinity of À 9.5 which is less as that of normal AT2R (Figure 1d).

Analysis of molecular dynamics simulation
Molecular dynamics is extensively used to access the stability and dynamic behaviour of protein ligand complexes (Bhardwaj et al., 2020;Kumar et al., 2021).Here the two variants (rs868939201 and rs1042852794) were selected for classical molecular dynamics simulations with agonist C21 and antagonist EMA401.

Root mean square deviation (RMSD)
The RMSD of a protein-drug complex is generally determined to evaluate the stability of the complex and its dynamic nature.The RMSD is computed as a function of simulation time and the starting structure of a protein-drug complex (Bhardwaj et al., 2023;Muvva et al., 2014).Here, we have observed RMSD of wild type AT 2 R and the variants namely rs1042852794 and rs868939201 in the absence or presence of drugs namely C21 and EMA401 over a simulation period of 100 ns (Figure 2).The RMSD of AT 2 R, alone or in the presence of C21 and EMA401 attained equilibrium after 10 ns and remained consistent throughout the simulation (Figure 2A).The RMSD of AT 2 R alone varied in the range of 0.236-0.323nm range, with an average value of 0.272 ± 0.013 nm.Also, the RMSD of AT 2 R in the presence of C21 and EMA401 remained consistent throughout the simulation from the start, implying stable nature of the complex formation.During 10-100 ns, the RMSDs of AT 2 R-C21 and AT 2 R-EMA401 complexes varied in the range of 0.082-1.101nm, and 0.122-0.204nm, with an average value of 0.086 ± 0.007 nm and 0.153 ± 0.010 nm respectively (Figure 2A).The RMSD of C21 and EMA401 alone also remained steady throughout the simulation, with an average value of 0.034 ± 0.007 nm and 0.052 ± 0.005 nm respectively.

Analysis of interacting residues
During MD simulation, the amino acid residues involved in the interaction with a ligand was determined and the results are presented in Table 4. Evidently, AT 2 R interacted with C21 with LYS168:NZ through salt bridge 100% simulation time.

Analysis of free energy calculations
The strength of interaction between protein and ligand in the presence of the solvent is calculated by determining its free energy.We have computed the free energy of the target protein-ligand interaction using the MM-GBSA method (Table 5).The free energy of AT 2 R-C21 complex was the lowest (-40.10 kcal mol À 1 ), followed by AT 2 R-EMA401 complex (-37.17kcal mol À 1 ), rs1042852794-C21 complex (-34.85kcal mol À 1 ), rs1042852794-EMA401 complex (-34.71kcal mol À 1 ), rs868939201-EMA401 complex (-33.36 kcal mol À 1 ), and rs868939201-C21 complex (-33.29 kcal mol À 1 ).The major driving forces for the formation of a stable complex in all the cases were van der Waals' energy (DG vdW ), hydrogen bonding (DG H-bond ) and non-polar solvation or lipophilic energy (DG Sol_Lipo ).On the other hand, the main force opposing the protein-ligand interaction was Coulombic forces (DG Coulomb ), covalent interactions, (DG Covalent ) polar solvation energy (DG Solv or DG SolGB ).

Discussion
A huge number of polymorphism and SNPs were reported in AT2R.Some of the SNPs were reported to be associated with the implication of diseases like diabetic nephropathies, hypertension and renal diseases (Sharma et al., 2022).As we know, SNPs were associated with different diseases and changes in the amino acid results in reduction in efficacy of drugs/ligands and also alters the structure and functioning of AT2R, instead of it these all the SNPs were not studied for their pharmacological significance.Similar studies were reported for the other receptor i.e.AT1R while no such report of effect of SNPs on AT2R protein structure, function and ligand binding is reported till date (Arsenault et al., 2010;Sharma et al., 2021).Therefore, 1151 AT2R SNPs were analyzed and out of them only 8 were taken for further homology modeling and molecular docking studies which were resides inside the binding cavity of AT2R.
Through homology modelling, trustworthy structures were modelled for each of the 8 SNP mutations that were chosen.Since the AT2R crystal structure is accessible in the PDB database (PDB ID ¼ 5UNF) (H.Zhang et al., 2017), it was utilized as a model for all the altered structures and has an about 99% similarity in sequence.The residues were discovered in favorable locations for all the predicted structures, according to Ramachandran plot analyses (Table 2).
The effect of SNP mutations on the binding affinities of C21 and EMA401 were analyzed through molecular docking.C21 is potential drug molecule for AT2R and in case of AT2R the variant rs868939201 show decreased binding affinity than normal AT2R which may leads to failure of this molecule in the population group bearing this SNP variant.As far as for EMA401 antagonist molecule one of the SNP variant rs868939201 show increased binding affinity than the wild type AT2R and there was also a variant rs1042852794 that also shows decrease in binding affinity than the wild type.This increase and decrease in binding affinity in case of different variants may leads to different functioning of agonist and antagonist molecule which may also cause failure of these molecules in population group bearing SNPs which reduces its binding affinity.
The difficulty of the molecular docking approach to compute the conformational changes of both the ligand and receptor during the molecular recognition process is one of its limitations (Modestia et al., 2019).Therefore, these two variants (rs868939201, and rs1042852794) having significant impact on binding affinity along with wild type AT2R and also in complex with agonist (C21) and antagonist (EMA401) were subjected to classical molecular dynamics to access the impact of mutations on protein conformations and dynamics.
The steady behavior of RMSD of the target proteins in the presence of C21 and EMA401 suggested the formation of a stable interaction between protein and drug during molecular dynamics simulation.Moreover, it also suggested that there were no major conformational changes in the overall structure of protein-drug complex.But it is also evident from the RMSD pattern that the rs1042852794 mutant exhibit higher rmsd as compare to wild type .The binding of C21 to mutants rs868939201, and rs1042852794 also have much higher RMSD as that of C21 complexes with the wild type AT2R.
The RMSF is also consistent with the RMSD results for rs1042852794 mutant as much higher fluctuation were observed for C21 and EMA401 complexed with mutant as compare to the wild type.The higher fluctuation due to rs1042852794 might not allowed the C21 and EMA401 to attain proper conformation inside the binding cavity of the AT2R.The other mutant rs868939201 also showed significant fluctuations in RMSF FOR C21 and EMA401 complexes as that of the wild type.It was also noted that the fluctuations were much higher for C21 complex as that of EMA401 complex which suggested that the rs868939201 mutant had pronounced effect in case of C21.The radius of gyration tells about the compactness of protein structure (Hassan et al., 2018).The comparative results justified that the Rg for upon binding of C21 and EMA401 remains stable and also slight decrease in Rg was also observed in case of wild type AT2R upon binding of agonist and antagonist.In contrast, same is not the case for mutants as slight rise in Rg values was observed which depicted that the mutation might affect the proteins compactness upon binding of ligand molecules.The increasing trends of the SASA of the entire protein indicates the solvation of the hydrophobic core as unfolding proceeds and ultimately leads to the interruption of hydrophobic interactions among non-polar residues (D.Zhang & Lazim, 2017).The SASA for C21 complexed with rs1042852794 and rs868939201 mutants of AT2R also showed increased value as that of wild type complex.It also have impact on overall binding and ultimately the agonist molecule does not work properly.Even for EMA401 complexed with mutated AT2R (rs1042852794 and rs868939201) high SASA values were observed as compared to the wild type which impact overall binding of agonist molecule.The similar pattern of hydrogen bonds were computed for C21n and EMA401 in case of wild type and mutant (rs1042852794 and rs868939201) AT2R.Moreover to access the binding strength of C21 and EMA401 with wild type and mutant AT2R binding free energy was calculated which clearly suggested that the occurrence of mutations had increase the binding free energy values (decrease in affinity) as that of wild type.Overall the analysis suggested that the even a SNP mutation might impact the binding of agonist and antagonist molecules and here the two mutations were of atmost important which were considered for future designing of drug design and discovery experiments targeting AT2R.Furthermore the mutants must be tested in in vitro and in vivo conditions to have better idea about the efficacy of agonist and antagonist molecules.

Conclusion
Overall, the analysis suggested that even a SNP mutation might impact the binding of agonist and antagonist molecules and here the two mutations were of at most important which were considered for future designing of drug design and discovery experiments targeting AT 2 R. Furthermore, the mutants must be tested in in vitro and in vivo conditions to have better idea about the efficacy of agonist and antagonist molecules.

Figure 2 .
Figure 2. RMSD profile for wild type and mutated AT2R (rs868939201 and rs1042852794) alone and in complex with C21 and EMA401.

Figure 3 .
Figure 3. Rg profile for wild type and mutated AT2R (rs868939201 and rs1042852794) alone and in complex with C21 and EMA401.

Table 1 .
SNP variants affecting the active site of AT2R.

Table 2 .
Binding affinities and hydrogen bonding of AT2R and its SNP variants with C21.

Table 3 .
Binding affinities and hydrogen bonding of AT2R and its SNP variants with EMA401.

Table 4 .
Various amino acid residues interacting with ligands C21 and EMA401 during MD simulation (only amino acid residues interacting with ligands for more than 30% simulation time are represented).

Table 5 .
Free energy (MM-GBSA) analysis of the interaction between AT 2 R and mutant proteins in their complex form with C21 and EMA401.All the energies are in kcal mol À 1 .DG Coulomb , DG vdW , DG Covalent , DG Solv or DG SolGB , DG H-bond , DG SA or DG Sol-Lipo , and DG or DG Bind stands for minimized molecular mechanics energy, coulomb energy, van der Waals' energy, covalent binding energy, solvation energy, energy due to self-contact, energy due to H-bonds, lipophilic energy and binding energy, respectively.