Shared binding modes of M3 muscarinic agonists and antagonists may incentivize development of novel IOP-reducing drugs: findings from in-silico assays

doi:10.21203/rs.3.rs-1662403/v1

PDF

Article

Shared binding modes of M3 muscarinic agonists and antagonists may incentivize development of novel IOP-reducing drugs: findings from in-silico assays

https://doi.org/10.21203/rs.3.rs-1662403/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

M3 muscarinic acetylcholine receptors (M3R) are G-protein-coupled receptors expressed in the ciliary muscle. Pilocarpine, an M3R agonist, reduces intraocular pressure (IOP) by causing contraction of the ciliary muscle, which leads to widening of the spaces in the trabecular meshwork and facilitates aqueous humor flow out of the eye. Despite their IOP-reducing potentials, there are only four M3R agonists that have been FDA-approved, and the chemical features that differentiate agonists from antagonists are poorly understood. In this study, we created a homology model of human M3R (hM3R) to analyze the binding modes that differentiate M3R agonists from antagonists, using such in-silico assays as high-throughput virtual screening and molecular dynamics. Our in-silico results suggest that a nitrogen ring with a positive charge is crucial for forming cation-π exchange with hM3R, regardless of agonism or antagonism. However, we discovered π-π stacking exchange that is exclusively found among M3R antagonists. Structural analysis suggests that this π-π stacking exchange leads to structural changes in one of the α-helices and causes instability in the cytoplasmic domain, where G protein interactions occur. Together, these finding may prove useful in developing drugs that not only target but also selectively agonize hM3R to reduce IOP in glaucoma patients.

Primary open angle glaucoma (POAG) is a leading cause of blindness worldwide¹. POAG is characterized by progressive optic nerve damage from retinal ganglion cell death, which is typically associated with elevated intraocular pressure (IOP). Steady-state IOP is generated by the balance of aqueous humor (AH) production by the ciliary body (CB) and its drainage through the trabecular meshwork and to a lesser degree the uveoscleral or nonconventional pathway. An imbalance between the inflow and outflow of AH leads to change in IOP, which could lead to persistent pressure in the eye and compromise optic nerve health^2–4. Hence, decreasing IOP is often considered as the first-line treatment for POAG.

Muscarinic acetylcholine receptors are G-protein-coupled receptors (GPCR) with five different subtypes, all of which are expressed in the anterior chamber of the eye^5,6. In particular, M3 muscarinic acetylcholine receptors (M3R) accounts for nearly 60% of all muscarinic receptors that are expressed in the ciliary muscle⁶. It has been experimentally validated that M3R activation causes contraction of the ciliary muscle, which leads to widening of the spaces in the trabecular meshwork and facilitates aqueous humor flow out of the eye^7–9. M3R activation can be achieved by small molecules that resemble its endogenous ligand, acetylcholine. For example, pilocarpine is an FDA-approved drug that has been used to reduce IOP by agonizing M3R⁹. On the other hand, M3R antagonists, such as atropine, have little effect on IOP and are contraindicated in glaucoma patients prescribed with pilocarpine^10–13. There is experimental evidence suggesting that muscarinic agonists and antagonists share the same binding pocket^14,15, yet the binding modes that differentiate each other are poorly understood.

In drug discovery and lead optimization, the computational method is a highly effective technique and has been widely applied^16–20. Specifically, molecular dynamics (MD) simulation considers the flexibility of protein structures based on Newtonian principles and is often used to analyze the protein-ligand interaction or the binding free energy collected during MD simulation^21–24. These results can then be utilized for compound optimization. Molecular mechanics and generalized born surface area (MM/GBSA) method is also used to predict the binding free energy of a compound binding to macromolecules^25–27. MM/GBSA based on MD simulation results can be used to study stability of a protein-ligand complex and analyze the variation in the sensitivity of a ligand caused by mutations in amino acid residues of a protein²⁸. When the structure of the protein is unknown, which is the case in our study, it needs to first be predicted. Homology modeling is considered the most accurate among the computational structure prediction methods in this situation²⁹. In homology modeling, a template of the target protein, which is the peptide sequence of another protein with high sequence identity, must first be found. A protein sequence with over 30% identity to a known structure can often be predicted with an accuracy equivalent to a low-resolution X-ray structure³⁰. The two sequences are then aligned, and a homology model is constructed. Post-modification of the model is often needed for more reliable results.

In this study, we developed a homology model of human M3R (hM3R) based on the crystal structure of rat M3R in complex with tiotropium (a known M3R antagonist). We then docked tiotropium in our homology model and compared its binding mode with tiotropium in rat M3R crystal structure to test the fidelity of our homology model. Subsequently, we performed high-throughput virtual screening (HTVS) of 21 compounds, composed of 6 M3R agonists and 15 M3R antagonists, to analyze differences in their binding modes. Finally, we performed MD simulations with the top scoring agonist and antagonist to analyze variations in the mode of interaction with hM3R. Our in-silico results suggest that a tertiary amine with a positive charge is crucial for formation of cation-π exchange with hM3R regardless of agonism or antagonism. However, the M3R antagonists formed π-π stacking exchange with hM3R, which is absent in M3R agonists because they lack aromatic structures. Structural comparison of hM3R suggests that this π-π stacking exchange leads to structural changes in one of the α-helices and causes instability in the cytoplasmic domain of hM3R, where G protein interactions occur.

Homology modeling of hM3R

Using crystal structure of rat M3R (PDB: 4DAJ), a homology model of hM3R was constructed, whose sequence identity was 62 %. The predicted structure of hM3R contains three different domains: extracellular (EC), transmembrane (TM), and cytoplasmic (CT) domains (Fig. 1A). Like most GPCRs, G protein interaction with hM3R is predicted to take place in the CT domain upon agonism^31-33. Upon superimposition of rat M3R crystal structure to hM3R homology model, we observed conservation of TM and EC domains (Fig. 1B). Our homology model was further refined using Protein Preparation and Loop Refinement features in Maestro (Schrödinger Inc., USA). Binding site analysis of our homology model was done using Sitemap feature in Maestro (Fig. 1C), and the top-ranking binding site can be found in the EC domain (Fig. 1D). This is consistent with the crystal structure of rat M3R, where tiotropium is in complex within the EC domain. Upon superimposition, tiotropium from the crystal structure fits tightly inside the top-ranking binding site of hM3R identified by Sitemap (Fig. 2A). Upon molecular docking of tiotropium to hM3R homology model using Glide (Schrödinger Inc., USA), the predicted binding mode is nearly identical to the binding mode of tiotropium obtained from the crystal structure of rat M3R (Fig. 2B). Upon comparison of the binding forces of tiotropium in the rat crystal structure to our hM3R homology model, it was determined that they were nearly identical, including the presence of π-π stacking exchange between tiotropium’s aromatic ring and tryptophan residues (Fig. 2C and 2D). Together, the results indicate that molecular docking with our hM3R homology model can predict ligand binding with high fidelity.

High-throughput virtual screening (HTVS) of M3R muscarinic agonists and antagonists

The 3D chemical structures of M3R agonists and antagonists were obtained from Drug Banks database, and their structures were further refined using LigPrep in Maestro. The binding modes of top three scoring M3R agonists revealed that cation-π exchange can be found between a positively charged nitrogen ring and residues with aromatic sidechains, like TYR149, TRP504, and TYR507 (Fig. 3A). Another notable force is hydrogen bonding between an oxygen group and ASN508. The binding modes of top three scoring M3R antagonists also contained cation-π exchange between a positively charged nitrogen ring and residues with TYR149, TRP504, and TYR507 (Fig. 3B). ASN 508 can also be seen interacting with oxygen molecules. However, compared to M3R agonists, the antagonists all shared presence of π -π stacking exchange between their aromatic rings and TRP504. In fact, this interaction is absent in all of the hM3R agonists available in the database likely because they lack aromatic side groups in their chemical structures. Table 1 displays the results of HTVS that lists the docking scores and checks presence of π- π stacking exchange, which is found in 12 out of 15 M3R antagonists that we docked. It’s also interesting to note that M3R antagonists achieved higher docking scores than the agonists likely because of the added stability via π-π stacking exchange.

MD simulation for interaction analysis

To analyze additional interactions found in hM3R agonism and antagonism, we conducted MD simulations using docked complexes with the highest scoring agonist and antagonist, pilocarpine and tiotropium, respectively. The protein-ligand complexes with the highest population during 100 ns simulation were extracted for binding free energy calculation and structural comparisons using Trajectory Clustering tool in Maestro (Schrödinger Inc., USA). The video analysis of MD simulations for hM3R-pilocarpine and hM3R-tiotropium complexes can be found in Supplemental Fig. 1 and 2, respectively. The root mean square distance (RMSD) analysis of hM3R in complex with pilocarpine was performed, where lower the RMSD indicates tighter binding (Fig. 4A). Pilocarpine achieved an RMSD of roughly 2.5 Å towards the end of the simulation, which indicates stable binding³⁴. The protein RMSD equilibrated to roughly 7.0 Å, which implies the protein assumed a stable conformation towards the end of the simulation. On the other hand, tiotropium achieved an RMSD of roughly 3.0 Å towards the end of the simulation (Fig. 4B). This value also indicates stable binding but not as tight as pilocarpine-hM3R complex likely because of the increased molecular mass of tiotropium. The protein RMSD equilibrated to roughly 10 Å, which indicates the protein assumed a stable conformation towards the end of the simulation

Interaction fraction analyses of pilocarpine-hM3R and tiotropium-hM3R complexes were also performed. Interaction fraction value of 1.0 means the ligand made contact with a binding residue during 100% time of the simulation. Values greater than 1.0 indicate more than one binding forces. For example, in the interaction fraction analysis of pilocarpine-hM3R (Fig. 5A), TYR149 has an interaction fraction value of 1.498 and makes contact with pilocarpine via hydrogen bonds, hydrophobic forces (cation-π exchange), and water bridges. The top three contributors to the binding interaction for pilocarpine were TYR149, ASP148, and SER152 with scores of 1.498, 1.217, and 0.914, respectively. In the interaction fraction analysis of tiotropium-hM3R (Fig. 5B), the top three contributors to the binding interaction were ASN508, TRP504, and TYR507 with scores of 1.755, 1.437, and 1.083, respectively. It’s interesting to note that TRP504 has almost 3-times the interaction fraction value with tiotropium (1.477) compared to pilocarpine (0.591). Furthermore, the top three binding residues with pilocarpine are upstream at positions 148-152, whereas they are downstream at positions 504-508 for tiotropium. Ligand-protein contact analysis of pilocarpine-hM3R revealed ASP148, TYR149, and SER152 intimately associating with the nitrogen ring of pilocarpine (Fig. 5C). Meanwhile, the ligand-protein contact analysis of tiotropium-hM3R displayed preservation of π-π stacking exchange via TRP504 during most of the simulation (Fig. 5D), which is consistent with the results of HTVS.

To analyze amino acid side chain stability, root mean square fluctuations (RMSF) analysis was performed, where higher the RMSF values indicates higher instability of amino acid side chain. In RMSF analysis of pilocarpine-hM3R complex (Fig. 6A), none of the amino acid side chains exceed RMSF value of 9.0 Å. On the other hand, in RMSF analysis of tiotropium-hM3R complex (Fig. 6B), there were three residues with RMSF values greater than 10.5 Å: LEU456 (11.03 Å), SER457 (12.16 Å), and LYS459 (11.89 Å). In fact, these three residues were much more stable in pilocarpine-hM3R complex: LEU456 (7.84 Å), SER457 (5.66 Å), and LYS459 (6.95 Å). It’s interesting to note that these residues span positions 456-459, which are near the top three binding residues (TRP504, TYR507, and ASN508) from the interaction fraction analysis of tiotropium-hM3R complex. Table 2 displays the difference in stability among shared binding residues of pilocarpine-hM3R and tiotropium-hM3R complexes. Overall, the binding residues of tiotropium have lower RMSF and hence higher stability. In fact, TRP504 had the highest difference in RMSF of 0.27 Å. This suggests TRP504 is in a more stable and rigid conformation when interacting with tiotropium likely due to added stability brought by π-π stacking exchange.

Through structural comparisons, we were able to locate TRP504 in helix 6 of hM3R (Fig. 7A and 7B). When we superimposed helix 6 of piloarpine-hM3R and tiotropium-hM3R (Fig. 7C), we identified a 2.1 Å shift in the α-carbon of TRP504 when interacting with tiotropium for π-π stacking. This led to an angle change of 18 ° upstream of TRP504, including the cytoplasmic domain, where the unstable residues LEU456, SER457, and LYS459 are directly located at. When we superimposed the cytoplasmic domain directly upstream of helix 6 (Fig. 7D), we could find LEU456, SER457, and LYS459 further away from the protein backbone in hM3R-tiotropium than in hM3R-pilocarpine, consistent with the high RMSF values we observed. It’s important to note that these three residues were more stable in hM3R-pilocarpine complex. Perhaps the increased fluctuations of these residues in the cytoplasmic domain are preventing G protein interactions and hence antagonism of hM3R.

Binding free energy analysis by MM/GBSA

To compare the stability of protein-ligand complex between pilocarpine and tiotropium, we used the MM/GBSA method to obtain the binding free energy (ΔG_bind). Table 3 lists the ΔG_bind of pilocarpine-hM3R and tiotropium-hM3R complexes, whose overall predicted binding free energies were −46.81 kcal/mol and -96.82 kcal/mol, respectively. This suggests that tiotropium forms a more stable complex with hM3R than pilocarpine does. Furthermore, near two-fold change in binding free energy indicates that once tiotropium binds to hM3R, its binding mode is stable, and the TM domain that contains the binding pocket likely remains in a stable, rigid conformation, which is consistent with the video analysis of tiotropium-hM3R complex in Supplemental Fig. 2.

In this study, we created a homology model of hM3R based on published structures of rat M3R since there were no crystal structure of hM3R published yet. There were minor differences in the cytoplasmic domains between hM3R and rat M3R, and this was expected because the cytoplasmic domain consists of several loop structures that are not as rigid as α-helices or β-sheets. Meanwhile, the TM and EC domains were well-conserved. Since the tiotropium binding pocket was found in TM and EC domains, we expected the binding site of hM3R to be conserved, which was confirmed when we ran Sitemap analysis and successfully superimposed tiotropium from rat M3R crystal structure to our homology model’s top-ranking predicted binding site. Molecular docking of tiotropium to hM3R revealed that the binding mode of tiotropium was nearly identical to that of rat M3R crystal structure, confirming the fidelity of our homology model for molecular docking.

When we ran HTVS of muscarinic agonists and antagonists of M3R, we discovered none of the agonists displayed π-π stacking exchange and lacked the capacity to form one since they do not have any aromatic side groups. On the other hand, of the 15 antagonists we screened, 12 of them shared presence of π-π stacking exchange with TRP504. Among the three antagonists that did not share any π-π stacking exchange, dicyclomine was the only molecule that did not have an aromatic side group. Furthermore, among the high-affinity antagonists (docking score < -10 kcal/mol), oxyphencyclimine was the only molecule, whose aromatic side group did not form π-π stacking exchange with TRP504. MD simulation of the top-scoring agonist and antagonist revealed that the π-π stacking exchange causes TRP504 to shift and lead to angle changes in helix 6 that is directly associated with the cytoplasmic domain, where G protein interactions occur. RMSF analysis revealed increased fluctuations of three residues found in the cytoplasmic domain that may interfere with G protein interactions.

Together, our in-silico assays suggest that a tertiary amine with a positively charged nitrogen molecule is crucial for binding to M3R regardless of agonism or antagonism. This is consistent with the structure of acetylcholine, an endogenous ligand of M3R that also houses a positively charged tertiary amine. However, to promote agonism, it’s crucial the ligand does not contain aromatic side groups and has no capacity to form π-π stacking exchange with TRP504. Because there are only 6 muscarinic agonists that can be found in Drug Banks database (compared to 61 antagonists), our in-silico data may prove useful in medicinal chemistry effort to develop drugs that not only target but also selectively agonize hM3R to reduce IOP in glaucoma patients.

Preparation for human M3R homology model

Human M3R structure was constructed via homology modeling using Sequence viewer function in Maestro and were based on 4DAJ structure. Assignment of bond orders and hydrogenation for the complex structures were performed using Protein Preparation in Maestro. The ionization state of the hM3R suitable for pH 7.0 ± 2.0 was predicted using Epik³⁵. H-bond optimization was conducted using PROPKA³⁶. Energy minimization was performed using the OPLS3e force field³⁷. Finally, loop structures were refined using Loop Refinement feature in Maestro. HTVS was performed using Glide in Maestro.

Md Simulation

MD simulations for interaction analysis and binding free energy calculation were performed using Desmond (Schrödinger Inc., USA). All systems were set up using “System Builder” software in Maestro. The pilocarpine-hM3R and tiotropium-hM3R complexes were placed in the orthorhombic box with a buffer distance of 10 Å to create a hydration model using the SCP water model. The following parameters were set: cut-of radius for van der Waals, time step, initial temperature, and pressure of the system set to 9 Å, 2.0 fs, 300 K, and 1.01325 bar, respectively. The sampling interval during the simulation was set to 100 ps. Finally, we performed MD simulations under the NPT ensemble for 100 ns.

Interaction Analysis And Trajectory Clustering For Mm/gbsa

Simulation Interactions Diagram tool in Maestro was used to perform an interaction analysis of our MD simulation, where our ligand was given 100 ns to equilibrate to stable conformation in relation to hM3R. Trajectory Frame Clustering function in Maestro was used to estimate the most populated representative structure for each MD simulation.

Mm/gbsa Free Energy Calculation

In MD simulation, free energy calculations give quantitative production of protein–ligand binding energies. The binding energy (ΔG_bind) was calculated by Equation [1]:

[1]\(\varDelta {G}_{bind}= {G}_{R+L}-({G}_{R}+{G}_{L})\)

where \({G}_{R+L}\) represents the hM3R in complex with ligands, while \({G}_{R}\) and \({G}_{L}\) represent the optimized free hM3R and optimized free ligand Gibbs energies, respectively³⁸.

In the generalized born surface area (MM/GBSA) approach²⁵, each free energy term in Equation [1] was calculated using Equation [2].

[2]:\(G={G}_{coulomb}+{G}_{vDW}+{G}_{covalent}+{G}_{solv}+{G}_{self-contact}+ {G}_{H-bond}+ {G}_{lipo}+ {G}_{packing}\)

where \({G}_{coulomb}\) represents coulombic energy, \({G}_{vDW}\) Van der Waals energy, \({G}_{covalent}\) covalent binding energy, \({G}_{solv}\) Generalized Born electrostatic solvation energy, \({G}_{self-contact}\) self-contact energy, \({G}_{H-bond}\) hydrogen-bonding energy, \({G}_{lipo}\)lipophilic energy, and \({G}_{packing}\) pi-pi packing energy.

The performance of the MM/GBSA algorithm is based on the specificity of the forcefield and ligand partial charges, the specificity of protein–inhibitor complex, MD simulation, inner dielectric constant, and the docking pose number based on top scoring³⁹. The VSGB solvation model⁴⁰ and OPLS3e force field were set for the calculation via MM/GBSA feature in Maestro.

Consent for publication

All authors consent for publication.

Data Availability

Data and materials will be available upon request.

Acknowledgement

This work was supported by NIH grants EY021200 (MMJ) and EY029950 (MMJ) and an unrestricted grant from Research to Prevent Blindness.

Author information

Affiliation: Department of Ophthalmology, Hamilton Eye Institute, University of Tennessee Health Science Center, Memphis, TN 38163, USA

Names: Minjae J. Kim & Monica M. Jablonski

Contributions: Work conceived by M.J.K. and M.M.J. Experiments performed by M.J.K. Data graphed and analyzed by M.J.K. M.J.K. and M.M.J. wrote the manuscript. Project financed by M.M.J. All authors edited and approved final manuscript.

Corresponding author: Correspondence to Monica M. Jablonski

Ethics Declaration

Competing interests: The authors declare no competing interests.

Kwon, Y. H., Fingert, J. H., Kuehn, M. H. & Alward, W. L. Primary open-angle glaucoma. N Engl J Med 360, 1113–1124, doi:10.1056/NEJMra0804630 (2009).
Goel, M., Picciani, R. G., Lee, R. K. & Bhattacharya, S. K. Aqueous humor dynamics: a review. Open Ophthalmol J 4, 52–59, doi:10.2174/1874364101004010052 (2010).
Llobet, A., Gasull, X. & Gual, A. Understanding trabecular meshwork physiology: a key to the control of intraocular pressure? News Physiol Sci 18, 205–209, doi:10.1152/nips.01443.2003 (2003).
Tamm, E. R. The trabecular meshwork outflow pathways: structural and functional aspects. Exp Eye Res 88, 648–655, doi:10.1016/j.exer.2009.02.007 (2009).
Gil, D. W., Krauss, H. A., Bogardus, A. M. & WoldeMussie, E. Muscarinic receptor subtypes in human iris-ciliary body measured by immunoprecipitation. Invest Ophthalmol Vis Sci 38, 1434–1442 (1997).
Zhang, X., Hernandez, M. R., Yang, H. & Erickson, K. Expression of muscarinic receptor subtype mRNA in the human ciliary muscle. Invest Ophthalmol Vis Sci 36, 1645–1657 (1995).
Lutjen-Drecoll, E. Structural factors influencing outflow facility and its changeability under drugs. A study in Macaca arctoides. Invest Ophthalmol 12, 280–294 (1973).
Lutjen-Drecoll, E., Wiendl, H. & Kaufman, P. L. Acute and chronic structural effects of pilocarpine on monkey outflow tissues. Trans Am Ophthalmol Soc 96, 171–191; discussion 192 – 175 (1998).
Skaat, A. et al. Effect of Pilocarpine Hydrochloride on the Schlemm Canal in Healthy Eyes and Eyes With Open-Angle Glaucoma. JAMA Ophthalmol 134, 976–981, doi:10.1001/jamaophthalmol.2016.1881 (2016).
Ho, M. C., Hsieh, Y. T., Shen, E. P., Hsu, W. C. & Cheng, H. C. Short-term refractive and ocular parameter changes after topical atropine. Taiwan J Ophthalmol 10, 111–115, doi:10.4103/tjo.tjo_110_18 (2020).
Rengstorff, R. H. & Doughty, C. B. Mydriatic and cycloplegic drugs: a review of ocular and systemic complications. Am J Optom Physiol Opt 59, 162–177 (1982).
Sanchez-Gonzalez, J. M., De-Hita-Cantalejo, C., Baustita-Llamas, M. J., Sanchez-Gonzalez, M. C. & Capote-Puente, R. The Combined Effect of Low-dose Atropine with Orthokeratology in Pediatric Myopia Control: Review of the Current Treatment Status for Myopia. J Clin Med 9, doi:10.3390/jcm9082371 (2020).
Wu, T. E. et al. Evaluating the Effect of Topical Atropine Use for Myopia Control on Intraocular Pressure by Using Machine Learning. J Clin Med 10, doi:10.3390/jcm10010111 (2020).
Birdsall, N. J., Burgen, A. S., Hiley, C. R. & Hulme, E. C. Binding of agonists and antagonists to muscarinic receptors. J Supramol Struct 4, 367–371, doi:10.1002/jss.400040307 (1976).
Birdsall, N. J. et al. Muscarinic receptor subtypes and the selectivity of agonists and antagonists. Pharmacology 37 Suppl 1, 22–31, doi:10.1159/000138503 (1988).
Chiba, S. et al. Identification of potential inhibitors based on compound proposal contest: Tyrosine-protein kinase Yes as a target. Sci Rep 5, 17209, doi:10.1038/srep17209 (2015).
Wakui, N., Yoshino, R., Yasuo, N., Ohue, M. & Sekijima, M. Exploring the selectivity of inhibitor complexes with Bcl-2 and Bcl-XL: A molecular dynamics simulation approach. J Mol Graph Model 79, 166–174, doi:10.1016/j.jmgm.2017.11.011 (2018).
Yasuo, N. & Sekijima, M. Improved Method of Structure-Based Virtual Screening via Interaction-Energy-Based Learning. J Chem Inf Model 59, 1050–1061, doi:10.1021/acs.jcim.8b00673 (2019).
Yoshino, R. et al. Pharmacophore modeling for anti-Chagas drug design using the fragment molecular orbital method. PLoS One 10, e0125829, doi:10.1371/journal.pone.0125829 (2015).
Li, H., Ibrahim, M. M., Chen, H., Li, W. & Jablonski, M. M. In Silico Screening and In Vivo Evaluation of Potential CACNA2D1 Antagonists as Intraocular Pressure-Reducing Agents in Glaucoma Therapy. Pharmaceuticals (Basel) 14, doi:10.3390/ph14090887 (2021).
Hayes, R. L. et al. Magnesium fluctuations modulate RNA dynamics in the SAM-I riboswitch. J Am Chem Soc 134, 12043–12053, doi:10.1021/ja301454u (2012).
Ingolfsson, H. I. et al. Lipid organization of the plasma membrane. J Am Chem Soc 136, 14554–14559, doi:10.1021/ja507832e (2014).
Sekijima, M., Motono, C., Yamasaki, S., Kaneko, K. & Akiyama, Y. Molecular dynamics simulation of dimeric and monomeric forms of human prion protein: insight into dynamics and properties. Biophys J 85, 1176–1185, doi:10.1016/S0006-3495(03)74553-6 (2003).
Yildirim, A., Sharma, M., Varner, B. M., Fang, L. & Feig, M. Conformational preferences of DNA in reduced dielectric environments. J Phys Chem B 118, 10874–10881, doi:10.1021/jp505727w (2014).
Genheden, S. & Ryde, U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin Drug Discov 10, 449–461, doi:10.1517/17460441.2015.1032936 (2015).
Kollman, P. A. et al. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc Chem Res 33, 889–897, doi:10.1021/ar000033j (2000).
Srinivasan, J., Miller, J., Kollman, P. A. & Case, D. A. Continuum solvent studies of the stability of RNA hairpin loops and helices. J Biomol Struct Dyn 16, 671–682, doi:10.1080/07391102.1998.10508279 (1998).
Kammarabutr, J., Mahalapbutr, P., Nutho, B., Kungwan, N. & Rungrotmongkol, T. Low susceptibility of asunaprevir towards R155K and D168A point mutations in HCV NS3/4A protease: A molecular dynamics simulation. J Mol Graph Model 89, 122–130, doi:10.1016/j.jmgm.2019.03.006 (2019).
Muhammed, M. T. & Aki-Yalcin, E. Homology modeling in drug discovery: Overview, current applications, and future perspectives. Chem Biol Drug Des 93, 12–20, doi:10.1111/cbdd.13388 (2019).
Xiang, Z. Advances in homology protein structure modeling. Curr Protein Pept Sci 7, 217–227, doi:10.2174/138920306777452312 (2006).
De Lean, A., Stadel, J. M. & Lefkowitz, R. J. A ternary complex model explains the agonist-specific binding properties of the adenylate cyclase-coupled beta-adrenergic receptor. J Biol Chem 255, 7108–7117 (1980).
Maguire, M. E., Van Arsdale, P. M. & Gilman, A. G. An agonist-specific effect of guanine nucleotides on binding to the beta adrenergic receptor. Mol Pharmacol 12, 335–339 (1976).
Weis, W. I. & Kobilka, B. K. The Molecular Basis of G Protein-Coupled Receptor Activation. Annu Rev Biochem 87, 897–919, doi:10.1146/annurev-biochem-060614-033910 (2018).
Ramirez, D. & Caballero, J. Is It Reliable to Take the Molecular Docking Top Scoring Position as the Best Solution without Considering Available Structural Data? Molecules 23, doi:10.3390/molecules23051038 (2018).
Shelley, J. C. et al. Epik: a software program for pK(a) prediction and protonation state generation for drug-like molecules. J Comput Aided Mol Des 21, 681–691, doi:10.1007/s10822-007-9133-z (2007).
Li, H., Robertson, A. D. & Jensen, J. H. Very fast empirical prediction and rationalization of protein pKa values. Proteins 61, 704–721, doi:10.1002/prot.20660 (2005).
Roos, K. et al. OPLS3e: Extending Force Field Coverage for Drug-Like Small Molecules. J Chem Theory Comput 15, 1863–1874, doi:10.1021/acs.jctc.8b01026 (2019).
Froloff, N., Windemuth, A. & Honig, B. On the calculation of binding free energies using continuum methods: application to MHC class I protein-peptide interactions. Protein Sci 6, 1293–1301, doi:10.1002/pro.5560060617 (1997).
Halim, S. A., Waqas, M., Khan, A. & Al-Harrasi, A. In Silico Prediction of Novel Inhibitors of SARS-CoV-2 Main Protease through Structure-Based Virtual Screening and Molecular Dynamic Simulation. Pharmaceuticals (Basel) 14, doi:10.3390/ph14090896 (2021).
Li, J. et al. The VSGB 2.0 model: a next generation energy model for high resolution protein structure modeling. Proteins 79, 2794–2812, doi:10.1002/prot.23106 (2011).

Table 1. HTVS of hM3R agonists and antagonists.

HM3R Agonists	Docking Scores (kcal/mol)	π-π stacking with TRP504?
Pilocarpine	-8.993	No
NGX267	-8.692	No
Cevimeline	-8.108	No
Bethanechol	-6.539	No
Xanomeline	-6.019	No
Methacholine	-3.746	No
HM3R Antagonists	Docking Scores (kcal/mol)	π-π stacking with TRP504?
Tiotropium	-11.036	Yes
Diphenidol	-10.626	Yes
Tridihexethyl	-10.606	Yes
Methanetheline	-10.400	Yes
Solifenacin	-10.378	Yes
Glycopyrronium	-10.264	Yes
Oxybutin	-10.159	Yes
Oxyphencycline	-10.034	No
Promethazine	-9.850	Yes
Homatropine methylbromide	-9.79	No
Procyclidine	-9.748	Yes
Dicycloamine	-9.45	No
Olanzapine	-9.443	Yes
Ipratropium	-9.391	Yes
Tropicamide	-9.064	Yes

Table 2. RMSF values (Å) of shared binding residues of hM3R for pilocarpine and tiotropium, from highest differences to lowest.

Binding Residues	RMSF (Å) with pilocarpine	RMSF (Å) with tiotropium	Difference (Å)
TRP504	1.18	0.91	0.27
Ser152	1.07	0.81	0.26
ALA236	1.95	1.77	0.18
CYS533	1.30	1.15	0.15
ASN508	1.60	1.58	0.02
ASP148	1.14	1.12	0.02
TYR534	0.73	0.73	0
TYR149	1.10	1.15	-0.05
TYR530	1.70	1.84	-0.14

Table 3. MM/GBSA binding free energy calculation of hM3R-pilocarpine and hM3R-tiotropium complexes. All values are in kcal/mol.

Protein-Ligand Complex

ΔG_coulomb

ΔG_vdW

ΔG_covalent

ΔG_solv

ΔG_self-contact

ΔG_H-bond

ΔG_Lipo

ΔG_Packing

ΔG_Bind

hM3R-

Pilocarpine

-0.99

-37.43

0.89

12.94

-1.82

-20.39

-46.81

hM3R-

Tiotropium

4.73

-60.93

2.72

-2.74

-1.22

-36.04

-3.34

-96.82

No competing interests reported.

SupplementalVideo1.mpeg
Supplemental Video 1
SupplementalVideo2.mpeg
Supplemental Video 2

PDF

Reviews received at journal
27 Jun, 2022
Reviewers agreed at journal
22 Jun, 2022
Reviewers invited by journal
21 Jun, 2022
Editor assigned by journal
21 Jun, 2022
Editor invited by journal
23 May, 2022
Submission checks completed at journal
23 May, 2022
First submitted to journal
16 May, 2022

You are reading this latest preprint version

Shared binding modes of M3 muscarinic agonists and antagonists may incentivize development of novel IOP-reducing drugs: findings from in-silico assays

Status:

Version 1

Abstract

Figures

Introduction

Results

Discussion

Methods

Preparation for human M3R homology model

Md Simulation

Interaction Analysis And Trajectory Clustering For Mm/gbsa

Mm/gbsa Free Energy Calculation

[1]\(\varDelta {G}_{bind}= {G}_{R+L}-({G}_{R}+{G}_{L})\)

[2]:\(G={G}_{coulomb}+{G}_{vDW}+{G}_{covalent}+{G}_{solv}+{G}_{self-contact}+ {G}_{H-bond}+ {G}_{lipo}+ {G}_{packing}\)

Declarations

References

Tables

Competing Interests

Supplementary Files

Status:

Version 1