Deepening insights into cholinergic agents for intraocular pressure reduction: systems genetics, molecular modeling, and in vivo perspectives

Parasympathetic activation in the anterior eye segment regulates various physiological functions. This process, mediated by muscarinic acetylcholine receptors, also impacts intraocular pressure (IOP) through the trabecular meshwork. While FDA-approved M3 muscarinic receptor (M3R) agonists exist for IOP reduction, their systemic cholinergic adverse effects pose limitations in clinical use. Therefore, advancing our understanding of the cholinergic system in the anterior segment of the eye is crucial for developing additional IOP-reducing agents with improved safety profiles. Systems genetics analyses were utilized to explore correlations between IOP and the five major muscarinic receptor subtypes. Molecular docking and dynamics simulations were applied to human M3R homology model using a comprehensive set of human M3R ligands and 1,667 FDA-approved or investigational drugs. Lead compounds from the modeling studies were then tested for their IOP-lowering abilities in mice. Systems genetics analyses unveiled positive correlations in mRNA expressions among the five major muscarinic receptor subtypes, with a negative correlation observed only in M3R with IOP. Through modeling studies, rivastigmine and edrophonium emerged as the most optimally suited cholinergic drugs for reducing IOP via a potentially distinct mechanism from pilocarpine or physostigmine. Subsequent animal studies confirmed comparable IOP reductions among rivastigmine, edrophonium, and pilocarpine, with longer durations of action for rivastigmine and edrophonium. Mild cholinergic adverse effects were observed with pilocarpine and rivastigmine but absent with edrophonium. These findings advance ocular therapeutics, suggesting a more nuanced role of the parasympathetic system in the anterior eye segment for reducing IOP than previously thought.

1 Introduction

Cholinergic activation assumes a central role in orchestrating diverse physiological functions within the anterior segment of the eye (Pakala et al., 2023). Triggered by the neurotransmitter acetylcholine release by postganglionic parasympathetic neurons, its interaction with muscarinic acetylcholine receptors occurs with exquisite precision and speed, modulating ocular dynamics within a matter of milliseconds (Vidal et al., 2023). Most notably, the binding of acetylcholine to the M3 muscarinic receptor (M3R) in the iris sphincter muscle enables its contraction, inducing pupil constriction (miosis) (Figure 1). This physiological response adapts the eye to varying light conditions, thereby enhancing visual acuity. Moreover, the activation of M3R plays a pivotal role in the contraction of the ciliary muscle, thereby governing two vital functions within the anterior segment of the eye. First, this contraction relaxes the tension in the suspensory ligaments linking the ciliary muscle to the lens, subsequently alleviating stretching tension within the lens (Wu et al., 2022). This, in turn, facilitates the accommodation of the lens, allowing for focused vision at different distances. Second, the contraction of the ciliary muscle exerts a tractional force on the trabecular meshwork, facilitating an increased outflow of aqueous humor through the conventional pathway (Keller and Peters, 2022). This regulatory mechanism provides an avenue for adjusting the intraocular pressure (IOP) in the eye, offering a potential target for innovative therapeutics designed to reduce IOP in glaucoma patients (Casson, 2022; Chen et al., 2023).

Figure 1

Figure 1. Physiological roles of cholinergic activation in the anterior segment of the eye. M3 muscarinic receptor is the most predominant subtype present in the human iris sphincter and ciliary body. The molecular mechanism of this receptor involves the direct stimulation of phospholipase C–β (PLCβ) and indirect stimulation of Ras homolog family member A (RhoA) through Gq proteins that ultimately release the inhibitory effect of myosin-light-phosphatase (MLCP) on myosin-light-chain kinase (MLCK). This intricate signaling cascade impacts calcium levels and subsequently MLCK to induce smooth muscle contraction.

Muscarinic acetylcholine receptors comprise five distinct subtypes—M1, M2, M3, M4, and M5—each with specific roles and distributions in various tissues, including the eye (Thul and Lindskog, 2018; Ruan et al., 2021). These receptors operate through G protein-coupled receptor (GPCR) mechanisms, initiating intracellular signaling upon ligand binding (Kudlak and Tadi, 2023). M1 muscarinic receptors (M1R), primarily in the central nervous system, stimulate phospholipase C (PLC) via Gq proteins, affecting calcium levels and protein kinase C activation. M2 muscarinic receptors (M2R), prevalent in the heart, inhibit adenylate cyclase through Gi proteins, modulating cardiac activity and neurotransmitter release. M3Rs, which are widely distributed in smooth muscle and the eye, activate PLC through Gq proteins, influencing smooth muscle contraction and glandular secretion. M3Rs also indirectly activate Ras homolog family member A (RhoA) via Gq proteins, leading to the release of myosin-light-phosphatase (MLCP) inhibition on myosin-light-chain kinase (MLCK) (Balla et al., 2023). This complex signaling cascade influences calcium levels, ultimately affecting MLCK and inducing smooth muscle contraction. M4 muscarinic receptors (M4R), concentrated in the central nervous system, inhibit adenylate cyclase via Gi proteins, impacting neurotransmitter release. Lastly, M5 muscarinic receptors (M5R), found in various tissues, activate PLC through Gq proteins, influencing intracellular calcium levels and neuronal activity. Using radioligand binding, mRNA detection, and direct protein detection with monoclonal antibodies, research has unveiled the muscarinic receptor distribution in the human iris sphincter and ciliary body (Gil et al., 1997; Nietgen et al., 1999). This analysis highlights a prevalence of M3R, comprising approximately 60%–75% of the receptors. In contrast, lower proportions of M2R and M4R, ranging from 5% to 10%, are evident in these tissues. Notably, M1R, constituting 7%, is detected in the ciliary processes and iris sphincter. Curiously, M5R, which is typically localized to the central nervous system, was also detected in the iris sphincter, accounting for 5% of muscarinic receptors in this ocular region.

Primary open-angle glaucoma (POAG) is a leading global cause of blindness, characterized by progressive optic nerve damage and elevated IOP (Kwon et al., 2009). Steady-state IOP is generated by the balance of aqueous humor production by the ciliary body, drainage through the trabecular meshwork, and, to a lesser degree, the uveoscleral or nonconventional pathway (Llobet et al., 2003). The imbalance between aqueous humor production and drainage can lead to sustained pressure, compromising optic nerve health (Llobet et al., 2003; Tamm, 2009; Goel et al., 2010). Hence, therapeutic approaches targeting the predominant muscarinic receptor, M3R, have been developed to alleviate IOP in individuals with glaucoma. These strategies primarily involve two methods, both utilizing the endogenous ligand of muscarinic receptors, acetylcholine. The first method entails the direct activation of M3R using small molecules that mimic the chemical structure of acetylcholine, exemplified by pilocarpine—an FDA-approved M3R agonist—to reduce IOP (Chen et al., 2023). The second method involves inhibiting acetylcholine degradation at the synapse between the axon terminal of post-ganglionic parasympathetic neurons and cholinergic targets in the eye. This is achieved using small molecule acetylcholinesterase inhibitors, such as physostigmine, to transiently elevate acetylcholine levels and hyper-activate M3R, thereby lowering IOP (Andrade and Zafar Gondal, 2023). Despite their proven efficacy in lowering IOP, both pilocarpine and physostigmine exhibit substantial adverse effects that constrain their clinical application. These drugs are associated with potent systemic cholinergic side effects, including excessive sweating, lacrimation, urinary incontinence, diarrhea, muscle weakness, airway constriction, and cardiac arrhythmia (Arens and Kearney, 2019; Panarese and Moshirfar, 2023). Consequently, cholinergic drugs have largely been supplanted by alternative IOP-reducing medications with distinct mechanisms, such as beta-blockers, prostaglandin analogs, and rho-kinase inhibitors (Liu et al., 2022). These alternatives avoid the issues of accommodation and pupillary reflex loss associated with cholinergic drugs. Nonetheless, these medications still present drawbacks that curtail their clinical utility. Beta-blocker IOP drops, for instance, have been linked to adverse effects, like headaches, dizziness, impotence, and other discomforts in glaucoma patients (Haga et al., 2022). In addition, FDA-approved prostaglandin analogs, like Vyzulta (latanoprostene bunod, 0.024%), Travatan Z (travoprost,0.004%), and Lumigan (bimatoprost, 0.03%), may cause permanent pigmentation of the eyelids, lashes, and iris (Llobet et al., 2003; Beidoe and Mousa, 2012; Shah et al., 2013). Lastly, while Rho kinase inhibitors are generally well-tolerated, they can induce conjunctival hyperemia, instillation site pain, and conjunctival hemorrhages, underscoring the imperative for safer drug options (Moshirfar et al., 2018). Hence, in this investigation, we revisit the parasympathetic nervous system within the anterior segment of the eye with the objective of repurposing existing FDA-approved or investigational drugs that selectively target this system to reduce IOP while minimizing potential side effects. This led our team to first conduct systems genetics analyses and molecular modeling to screen for potential lead compounds capable of IOP reduction. Subsequently, we employed mouse models to assess their efficacy in lowering IOP and to screen for potential cholinergic adverse effects.

2 Materials and methods

2.1 Correlation analyses of muscarinic acetylcholine receptor genes and IOP

Systems genetics is a methodology that investigates the genetic basis of disease traits by simultaneously analyzing extensive sets of genotypes and phenotypes (Andreux et al., 2012; Houtkooper et al., 2013; Wang et al., 2016). The BXD family, a well-characterized mouse genetic reference population, combines C57BL/6J and DBA/2J mice genomes via spontaneous recombination over 20+ generations of inbreeding, segregating for approximately 5.5 million sequence variants (Taylor et al., 1999). This family of experimental mice has been genotyped and phenotyped thoroughly at many levels, from messenger RNA levels to IOP (Chesler et al., 2005; Ashbrook et al., 2015; Williams et al., 2016). This makes the BXD family a valuable resource for glaucoma-related systems genetics studies. First, transcript data for Chrm1, Chrm2, Chrm3, Chrm4, and Chrm5 (corresponding to M1R, M2R, M3R, M4R, and M5R, respectively) were acquired from Gene Network 2.0 (https://genenetwork.org/) (Mulligan et al., 2017). The comprehensive transcript data are openly accessible at UTHSC BXD All Ages Eye RNA-Seq (Nov20) TPM Log2 - eye mRNA database. Similarly, IOP measurements from 51 strains of BXD mice of all ages (1–30 months old), both sexes, average of left and right eyes were obtained, which are openly accessible at BXD Published Phenotypes–IOP database. Pearson correlation matrix analysis (p-value <0.05) was then conducted, following the instructions in the user’s manual at Gene Network 2.0, to identify statistically significant correlations between IOP and the five genes of interest. Lastly, another round of analysis was performed to identify statistically significant correlations among the mRNA expressions of the five genes.

2.2 Homology modeling of human M3R

In drug discovery, molecular modeling is commonly used to visualize the binding interaction of a ligand to a known protein structure (Kim et al., 2023). The ligand is docked against a specific site on the protein of interest, and the best binding and docking score are calculated based on the energy reduction of the two entities before and after the docking. A more negative docking score means a higher energy reduction and thus a tighter binding between the ligand and the protein (Guedes et al., 2014). Ideally, molecular modeling is most effectively performed using a protein structure established through X-ray crystallography, electron microscopy, or nuclear magnetic resonance imaging, and published in the Protein Data Bank for reference. Within the protein data bank, there were a total of ten published structures corresponding to the entry name “Muscarinic acetylcholine receptor M3,” composed of four human M3R and six rat M3R structures. All the human M3R structures, however, underwent significant modifications, like mutagenesis, truncation, separation, and chemogenetic alterations, that limited their utilities for molecular modelling (Zhang et al., 2022). Hence, homology modelling was pursued to predict the full structure of the human M3R before proceeding with the docking analysis. Homology modeling is considered the most accurate among the computational structure prediction methods in this situation (Muhammed and Aki-Yalcin, 2019). The full amino acid sequence of human M3R was first acquired from UNIPROT database, and its molecular weight as well as isoelectric point were predicted based on its amino acid sequence using SWISS-EXPASY MW/pI calculator (Audain et al., 2016) (Supplementary Figure S1). Human M3R structure was constructed via homology modeling using Sequence viewer function in Maestro (Schrödinger, Inc.) with a partially solved structure of human M3R in complex with iperexo, a muscarinic receptor super-agonist, (PDB: 8E9Z) as template structure with sequence identity of 100%. It is essential to note that our homology model incorporates a truncating modification in intracellular loop three spanning amino acid positions from 253 to 491, a common practice for membrane proteins with complex intracellular loops to enhance stability prior to MD simulation (Li et al., 2024; Stampelou et al., 2024). However, given that our study focuses on the transcellular binding pocket of M3R, we anticipate that the altered intracellular loops will not impact this interaction. Subsequently, assignment of bond orders and hydrogenation for the complex structures were conducted using Protein Preparation in Maestro. The ionization state of the human M3R suitable for pH 7.0 ± 2.0 was predicted using Epik (Shelley et al., 2007). H-bond optimization was conducted using PROPKA (Li et al., 2005). Energy minimization was conducted utilizing the OPLS3e force field to release the protein conformation from its initially crystallized state of M3R agonism associated with the template structure (Roos et al., 2019). Subsequently, loop structures were refined using Loop Refinement feature in Maestro. Lastly, a quantitative structural analysis of the homology model was performed on both the template structure and homology model using Alignment root-mean-square-deviation (RMSD) and Ramachandran plot analysis in Maestro. RMSD calculates the difference in atomic distance between template structure and homology model, where lower the RMSD the closer the homology model is to the template structure in 3D conformation. Each dot on the Ramachandran plot represents a peptide with corresponding phi-psi angles. Blue dots indicate “allowable” peptides with favorable phi-psi angles, while red dots indicate “disallowed” peptides due to unfavorable phi-psi angles. Peptides with phi/psi angles around −90/120 degrees correspond to beta sheets, −90/-30 degrees to right-handed alpha helices, and 90/60 degrees to left-handed alpha helices.

2.3 Molecular docking of compound libraries against human M3R

First, the binding site of the human M3R homology model was determined based on relevant literature sources (Martinez-Archundia et al., 2012; Kistemaker et al., 2019). Modified protocols from previous modelling studies and relevant literature were adopted for molecular docking (Saini et al., 2019; Li et al., 2021). The identified binding residues were subsequently utilized as input for the Receptor Grid Generation function in Maestro, creating a three-dimensional exploration space for molecular docking. The dimensions and center point coordinates of the grid box were set to 20 Å × 20 Å × 20 Å and 52.34, 95.03, 59.03 as coordinates for X, Y, Z-axis of the center point, respectively. Receptor grids were generated, and various classes of ligands were prepared for docking using the LigPrep function in Maestro (Wishart et al., 2018). First, a compound library was acquired from the DrugBank 6.0 database, consisting of six known agonists of M3R to serve as positive controls (pilocarpine, NGX267, bethanechol, cevimeline, xanomeline, methacholine) and 15 known antagonists of M3R to serve as negative controls (tiotropium, diphenidol, tridihexethyl, methanetheline, solifenacin, glycopyrronium, oxybutynin, oxyphencyclimine, promethazine, homatropine methylbromide, procyclidine, dicyclomine, olanzapine, ipratropium, tropicamide) for our docking study (Vacca et al., 1987; Greig et al., 2013; Knox et al., 2024). Following ligand preparation, the receptor grid files were employed in the Glide function, with extra-precision selected as the screening method and “amide-only” chosen for interaction settings. The torsional degrees of freedom of the ligand molecule were determined by its flexibility. Structural motifs associated with binding interactions were computed and generated using the “Ligand Interaction Diagram” software installed within Maestro interface (Basaran et al., 2024). This tool employs a sophisticated array of computational techniques to predict and visualize key interactions between proteins and ligands, including hydrogen bonds, van der Waals forces, covalent, and ionic bonds. The “Ligand Interaction Diagram” tool visualizes the results of Glide dock, which uses geometric criteria to identify hydrogen bonds, evaluating optimal distances and angles between donor and acceptor atoms, which are visualized as dashed lines in accompanying tools. Van der Waals interactions are modeled using Lennard-Jones potentials, highlighting favorable contacts through molecular surface displays or depicting close atomic contacts. While covalent bonding is generally predefined in the molecular structures, its visualization is managed through representations showing bonds as sticks or lines. Ionic bonds are detected through assessments of electrostatic potentials between charged groups, with visual representations often differentiated from hydrogen bonds by a mix of red and blue colors. This suite of computational predictions allows integrative assessment of various structural motifs predicted by Glide docking.

After completing the docking for positive and negative controls, a compound library consisting of 1,667 FDA-approved or investigational drugs was acquired from ZINC 5.0 database and prepared for docking as previously described. Due to the large size of the ligand sample, a series of molecular docking procedures was performed. Each successive series was conducted with higher precision but longer computational time. Glide docking operates at three precision levels: High Throughput Virtual Screening (HTVS), Standard Precision (SP), and Extra Precision (XP). HTVS is a less stringent method that discards ligands if their volume exceeds that of the binding site. SP is more computationally intensive, providing more accurate predictions by employing more thorough sampling and refinement processes. XP, on the other hand, performs even more extensive sampling than SP and uses a sophisticated scoring function that demands stricter ligand-receptor shape complementarity, thereby reducing false positives. XP also evaluates ligand fit to specific receptor conformations, which is useful for proteins with conformational flexibility. Within the context of our study, a compound library of 1,667 compounds underwent docking via HTVS mode with “amide-only” chosen for interaction settings, initially reducing the number of ligands to 82 compounds. Hit compounds generated underwent a secondary docking process using the SP mode in Glide for further refinement of the drug library, further reducing the library to 23 compounds. Finally, the docked ligands underwent a conclusive, third docking process using XP mode in Glide to complete the refinement to four compounds. These finalized hit compounds were then categorized as potential M3R agonists or antagonists based on the structural motifs identified from the initial docking process using the six positive and 15 negative controls. All four compounds housed a positively charged tertiary nitrogen, facilitating hydrogen and ionic bonding with M3R binding residues—a characteristic consistently observed across all 21 controls. Intriguingly, two of the compounds, sulpiride and profenamine, exhibited π-π stacking interactions with a tryptophan residue in M3R, a feature predominantly associated with 12 of the 15 M3R antagonists examined. Conversely, the other two compounds, rivastigmine and edrophonium, did not display this interaction, aligning instead with the structural themes of the agonists. To further explore the dynamic interactions and stability of these structural motifs, our team initiated MD simulations extending over 100 nanoseconds using pilocarpine and tiotropium, as known M3R agonist and antagonist controls, respectively.

2.4 Molecular dynamics (MD) simulation of human M3R in complex with pilocarpine and tiotropium

Utilizing Desmond (Schrödinger Inc., USA), MD simulations were executed to explore interactions and calculate binding free energy, aiming to elucidate the differences in binding modes for M3R agonism and antagonism. Modified protocols from previous studies and relevant literature were adopted (Saini et al., 2019; Kim et al., 2023). The binding sites of pilocarpine and tiotropium were verified according to the literature (Kruse et al., 2012). Configuration of all systems was performed using the “System Builder” tool in Maestro. The protein-ligand complex was positioned in an orthorhombic box with a 10 Å buffer distance, incorporating a hydration model with the SCP water model that was equilibrated with salt concentrations of 0.15 M Na⁺ and 0.15 M Cl⁻ at physiological pH. Studies with similar MD simulation workflow as our study reported that both SPC and TIP4P (Transferable Intermolecular Potential 4 Point) water models reproduce the best experimental values (Zielkiewicz, 2005). Choosing the SPC water model over the TIP4P model was based on considerations of computational efficiency and compatibility with other simulation components. SCP models, including variations like SPC/E, are less computationally demanding due to their simpler three-point charge system, which excludes the additional interaction site found in TIP4P. This makes SCP models faster, which is advantageous for large systems or extensive simulations (Steinberg et al., 2019). While TIP4P may offer superior accuracy for certain thermodynamic properties, SCP models still provide reasonable approximations for many applications where ultra-precise water behavior is not critical. Furthermore, the SCP model’s compatibility with widely used force fields and its straightforward implementation make it a practical choice, especially for those new to molecular dynamics or when working under computational constraints. This balance of efficiency, sufficient accuracy, and ease of use led our team to choose SCP over TIP4P.

Subsequently, key parameters, including van der Waals cut-off radius, time step, initial temperature, and system pressure, were set to 9 Å, 2.0 fs, 300 K, and 1.01325 bar, respectively. While the size of the receptor itself measured at approximately 11,000 atoms, the size of the system containing the protein-ligand complex with the buffer measured at approximately 36,000 atoms and total volume of 628.84 ų. MD simulations, lasting 100 ns under the NPT ensemble, were executed with a simulation sampling interval of 100 ps. A simulation time of 100 ns was selected for this study because pilocarpine and tiotropium are both known ligands of M3R. Our team determined that this duration would allow sufficient time for the protein and ligand to establish stable interactions. This choice aligns with other recent studies using molecular dynamics to simulate and analyze known interactions in real-time under physiological conditions (Gheidari et al., 2024; Roy et al., 2024). Subsequently, trajectories from each MD simulation were scrutinized using the Simulation Interactions Diagram tool in Maestro. Trajectory Frame Clustering in Maestro determined the most populated representative structure for each MD simulation before conducting molecular mechanics with generalized Born and surface area solvation (MM/GBSA) analysis for free energy calculation using Prime (Schrödinger Inc., USA). This entire procedure, from system configuration to MM/GBSA analysis, was performed three times for each protein-ligand complex. While the trajectory analysis yielded identical results across all replicates for each complex, we observed slight variations in the MM/GBSA analysis.

2.5 Molecular mechanics with generalized born and surface area solvation (MM/GBSA) analysis for free energy calculation

In MD simulation, free energy calculations provide quantitative production of protein–ligand binding energies. The binding energy (ΔG_bind) is calculated by Eq. 1:

$∆{\mathit{G}}_{\mathit{bind}}={\mathit{G}}_{\mathit{R}+\mathit{L}}-\left({\mathit{G}}_{\mathit{R}}+{\mathit{G}}_{\mathit{L}}\right)(1)$

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 (Froloff et al., 1997).

In the MM/GBSA approach (Genheden and Ryde, 2015), each free energy term in Eq. 1 is calculated using Eq. 2.

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

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 (Halim et al., 2021). The VSGB solvation model (Li et al., 2011) and OPLS3e force field were set for the calculation via MM/GBSA feature in Maestro.

2.6 Assessment of lead compounds for IOP reduction in mice

A total of six C57BL/6 J mice (n = 6; 3 males and 3 females, aged 2–6 months) were employed to assess the IOP-lowering effects of lead compounds identified through molecular docking and molecular dynamics simulations. As the drug molecules were soluble in aqueous media, they were easily prepared as aqueous solutions. However, to prevent rapid drainage from the eye surface before eliciting a pharmacological response, these molecules were incorporated into bioadhesive viscous Carbopol 981 eye drops. A 0.4% weight/volume (w/v) concentration of Carbopol 981 was soaked in Milli-Q water overnight for full swelling. Equal volumes of 0.4% w/v Carbopol 981 gel and drug aqueous solutions were mixed using a vortex mixer to create the final bioadhesive viscous eye drops, stored in closed, air-tight containers at 5 °C in darkness until use. The resulting products contained 0.6% w/v of each drug in 0.2% w/v Carbopol 981 gel. To ensure sterility, all tools and water used were sterile, and the procedures were carried out under aseptic conditions.

For each mouse, the tested ophthalmic formulation (10 μL) was applied to both right and left eyes. IOP of mice was subsequently measured every hour for the next 8 h using an induction-impact tonometer (Tonolab tonometer, Colonial Medical Supply, Franconia, NH) according to the manufacturer’s recommended procedures. Six consecutive IOP readings were averaged. All results were expressed as the mean percentage reduction in IOP from baseline (mean % reduction ± SEM). All procedures involving mice were approved by the Animal Care and Use review board of the University of Tennessee Health Science Center (UTHSC) and followed the Association of Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research in addition to the guidelines for laboratory animal experiments (Institute of Laboratory Animal Resources, Public Health Service Policy on Humane Care and Use of Laboratory Animals). Animals were housed under cyclic light (12 h on:12 h off) with 35% humidity in a specific pathogen-free facility at UTHSC and were allowed free access to water and food.

2.7 Quantitative structure-activity relationship (QSAR) and ADME (absorption, distribution, metabolism, excretion) analysis of cholinergic drugs

The 3D structures of pilocarpine, physostigmine, rivastigmine, edrophonium, sulpiride, and atropine were subjected to QSAR analysis using physics-based membrane permeability tool in Maestro (Leung et al., 2012; Leung et al., 2016). The physiochemical properties derived mathematically from QSAR were then used to predict the corneal permeability of the investigated cholinergic drugs (Karami et al., 2022). To calculate changes in free energy of partitioning (ΔG_o/w), we first obtained the distribution coefficient (cLogD at pH 7.4) and partition coefficient values (cLogP). The free energy of distribution value was then calculated based on mathematical equations adapted from the literature (Anderson et al., 1988; Leung et al., 2012; Karami et al., 2022). In the QSAR analysis, cLogD_pH7.4 represents the distribution coefficient for a drug at pH 7.4, which is the pH of the tear film. The higher the cLogD_pH7.4, the higher the hydrophobicity of the drug. ΔG_o/w represents the change in Gibbs free energy of distribution/partitioning of a drug, where higher the value the higher the hydrophilicity. Finally, TPSA represents the topological polar surface area of a drug, where the higher the TPSA the greater the physiochemical space occupied by the polar surface of the drug. The rule of thumb for optimal ocular permeability and absorption for topical eye drops proposed by Karami et al. (RO_x) suggests that drugs with cLogD_pH7.4 ≤ 4.0, ΔG_o/w ≤ 20 kcal/mol, and TPSA ≤ 250 Å are in favor of crossing the corneal membrane upon topical delivery (Karami et al., 2022). Lastly, the SDF files of pilocarpine, physostigmine, rivastigmine, edrophonium, sulpiride, and atropine were uploaded to SWISS-ADME, a freely accessible web tool for evaluating pharmacokinetics, drug-likeness, and medicinal chemistry properties of small molecules (Daina et al., 2017). Specifically, LogP represents the partition coefficient, indicating the compound’s lipophilicity and ability to cross the cell membrane, while LogS denotes its aqueous solubility, impacting absorption and bioavailability. GI absorption assesses the likelihood of a compound being absorbed through the gastrointestinal tract. CYP interactions highlight the compound’s interactions with cytochrome P450 enzymes, which are crucial for its degradation and metabolism. BBB (blood-brain barrier) permeability measures the compound’s ability to cross the blood-brain barrier, and P-gp substrate status identifies whether the compound is a substrate for P-glycoprotein, a transporter protein that can affect its absorption and excretion. These data collectively aimed to enhance our understanding of the pharmacokinetics of the cholinergic drugs in relation to their observed biological activities in mice.

2.8 Data and statistical analysis

The data and statistical analysis in this study comply with the recommendations of the British Journal of Pharmacology on experimental design and analysis in pharmacology. A change is considered statistically significant, if the statistical analysis yields a p-value <0.05. Pearson product-moment correlations were computed using the Correlation matrix tool within GeneNetwork 2.0. Denoted by “r,” the Pearson correlation assesses the strength of a linear association between two variables, specifically in our case, IOP (mmHg) and muscarinic acetylcholine receptor mRNA expressions. The calculated pharmacodynamic parameters for the lead compound encompass maximum percent IOP reduction, the time to reach the maximum decrease in percent IOP (T_max), the duration for IOP to return to baseline (i.e., end of drug effect; T_end), and the total area under the percent IOP reduction-versus-time curve (AUC). GraphPad Prism-10 software (GraphPad Software Inc., San Diego, CA, USA) was employed for all pharmacodynamic calculations, and the results were presented as the mean ± SEM.

2.9 Materials

All modelling operations in Maestro were performed using RedHat Enterprise Linux (RHEL) Version 8.4 Desktop computer containing two Linux X64 (AMD64/EM64T) Display GPU Drivers. Tested compounds, namely pilocarpine, rivastigmine, edrophonium, sulpiride, and atropine, were purchased from Selleckchem (New York, United States) with catalog numbers of S4231, S4713, S2087, S5900, S4655, and S4713, respectively. Meanwhile, our last tested compound, physostigmine, was purchased from MedChemExpress (New Jersey, United States) with a catalog number of HY-B1266.

3 Results

3.1 Correlation analyses of muscarinic acetylcholine receptor genes and IOP

Correlation matrix analyses revealed an overall negative correlation pattern between muscarinic receptor mRNA expression and IOP (Figure 2). However, the negative correlation attained statistical significance solely with Chrm3 mRNA expression, corresponding to M3R. This particular outcome aligns with the understanding that M3R stands out as the predominant muscarinic receptor subtype within the anterior segment of the eye. The specificity of the negative correlation with Chrm3 mRNA also reinforces the notion that the influence of muscarinic receptors on IOP is particularly associated with the M3 subtype, substantiating its pivotal role in modulating IOP dynamics.

Figure 2

Figure 2. Correlation analyses of muscarinic acetylcholine receptor genes and IOP. Correlation matrix analyses revealed an overall negative correlation between muscarinic acetylcholine receptor mRNA expression and IOP. Notably, the statistically significant negative correlation was observed only with Chrm3 mRNA expression, corresponding to M3R.

Meanwhile, a general positive correlation pattern emerged between Chrm3 and the mRNA expressions of other muscarinic acetylcholine receptor subtypes (Figure 3). Notably, all these positive correlations reached statistical significance, with the most robust correlation observed between Chrm2 and Chrm3 mRNA expressions. This finding suggests a potential regulatory role of M3R expression on other muscarinic subtypes, indicating a coordinated and interrelated regulation of their respective expressions. The intricate interplay among these receptors could play a role in modulating IOP, given all five muscarinic acetylcholine receptors are expressed in the anterior segment of the eye.

Figure 3

Figure 3. A positive correlation pattern manifested in the association between Chrm3 and the mRNA expressions of other muscarinic acetylcholine receptors. Remarkably, all these positive correlations attained statistical significance, and the most robust correlation was evident between Chrm2 and Chrm3 mRNA expressions.

3.2 Human M3R homology model and its binding site

The homology model of the human M3R housed its binding site within the hydrophobic core of its transmembrane domain (Figure 4). Molecular docking of pilocarpine, a known agonist, into this binding site revealed a docking score of roughly −10 kcal/mol. Particularly noteworthy in the binding mode were the hydrogen bond interactions from ASP148 and SER152, respectively, involving pilocarpine’s hydrogen atoms associated with its tertiary amine group. This is consistent with the binding mode of iperexo, an M3R super-agonist, in our template structure (PDB:8E9Z), in which the ligand’s tertiary amine group is seen associating with ASP148, SER152, and TYR530 via cation-π and hydrogen bond networks (Supplementary Figure S2 (Zhang et al., 2022). Both iperexo and pilocarpine displayed similar 3D conformations within the M3R binding site, with their positively charged tertiary amine oriented towards ASP148 and their aromatic rings towards VAL156. Furthermore, the binding modes of both compounds involved extensive π-cation and salt bridge interactions with the nitrogen ring via ASP148, TYR149, and TYR507, while there was no evidence of π-π stacking interaction with TRP504 in both compounds. Finally, pilocarpine exhibited a unique hydrogen bonding interaction between its carbonyl ring and TRP200.

Figure 4

Figure 4. The homology model of the human M3R and its binding site. The cyan structure on the left represents the ribbon diagram of human M3R, while the middle structure provides a van der Waals ribbon diagram, where red regions denote hydrophobic residues and blue regions hydrophilic residues. The right structure illustrates the contact surface diagram of human M3R, with red regions corresponding to hydrophobic residues and blue regions to hydrophilic residues. The molecular docking of pilocarpine, a recognized agonist, into this binding site resulted in a docking score of −9.839 kcal/mol. Particularly noteworthy in the interaction profile were the extensive hydrogen bonding interactions involving pilocarpine’s positively charged tertiary amine group and the polar structures of ASP148 and SER152. This structure subsequently underwent MD simulations to identify additional binding residues that may enhance the stability of this interaction.

Structural alignment of our homology model to the template structure revealed an all-atom RMSD of 4.283 Å and a central carbon (Cα) backbone RMSD of 0.005 Å. This discrepancy between the all-atom and Cα RMSD indicates that most of the variations in the peptide structures arise from the flexibility of the amino acid side chains. While an all-atom RMSD of less than 2.0 Å is ideal for homology models, an RMSD of less than 5.0 Å is also considered reasonably appropriate (Zhang and Skolnick, 2004; Kumar et al., 2019). Subsequent Ramachandran plot analysis compared the structural motif of our template structure (PDB: 8E9Z) to that of our homology model (Supplementary Figure S3. Overall, similar distributions of peptides were observed for all secondary structures, specifically beta-sheets and alpha-helices. However, our homology model displayed 30″disallowed” peptides, as opposed to the template model’s two “disallowed” peptides. In terms of structural quality of our homology model, an all-atom RMSD of 4.283 Å and a percentage “disallowed” peptides of 8.5% (30/352) suggest that there is still room for improvement in the quality of our homology model. This higher number of “disallowed” regions in the Ramachandran plot may be attributed to the methodology employed for constructing the model in Maestro. Specifically, the approach of attempting to model the entire protein in a single, unified step, rather than adopting a systematic, domain-wise modeling strategy followed by a careful integration of these domains using appropriate constraints, could have led to inaccuracies. Such a domain-wise approach often enables more precise control over the conformational space of each domain, potentially reducing the occurrence of structurally implausible regions (Xia et al., 2023). This underscores the need for refining our modeling strategies to enhance the structural accuracy of complex protein models. However, the residues within the canonical binding site of M3R, namely ASP148, TYR149, SER152, and TYR530, all had favorable phi-psi angles. Given that the focus of our study is within the transmembrane domain rather than looped domains with flexible amino acid side chains, our team proceeded with molecular docking and MD simulations.

3.3 Molecular dynamics simulation of human M3R in complex with pilocarpine and tiotropium

To delve into additional interactions associated with human M3R agonism and antagonism, we conducted triplicate molecular dynamics simulations for human M3R docked with pilocarpine and tiotropium, respectively. Readers are encouraged to view the video analysis of the simulation results for the hM3R-pilocarpine and hM3R-tiotropium complexes (Supplementary Videos S1, S2). The root-mean-square-deviation (RMSD) analysis was performed for both complexes, indicating stability at a low RMSD along the timeframes (Figure 5). Towards the end of the simulation, pilocarpine achieved an RMSD of approximately 2.5 Å, while tiotropium stabilized at an RMSD of around 3.0 Å, signifying stable binding for both compounds.

Figure 5

Figure 5. Trajectory analyses of molecular dynamics simulations for human M3R in complex with pilocarpine (A) and tiotropium (B). Root-mean-square-deviation (RMSD) analysis demonstrated stable interactions in both the pilocarpine-hM3R and tiotropium-hM3R complexes, with pilocarpine displaying a tighter binding, likely due to its lower molecular mass. Interaction fraction analyses were conducted to reveal all amino acid residue interactions observed during MD simulation. The bar diagrams for pilocarpine and tiotropium revealed TYR149 as a key contributor in both complexes, engaging in diverse interactions with pilocarpine and tiotropium. Notably, TRP504 exhibited a significantly higher interaction fraction with tiotropium compared to pilocarpine. The top binding residues differed in their locations, with pilocarpine-hM3R interactions upstream (148-152) and tiotropium-hM3R interactions downstream (504-508), as denoted by red boxes. Subsequent 2D Ligand-protein contact analysis was conducted to visualize these key binding interactions that occurred more than 30% of the simulation time. Red squares in the 2D diagram correspond to top binding residues identified from the interaction fraction analysis. Overall, these results underscore the significance of a positively charged tertiary amine in the binding process to human M3R, hinting at the potential of aromatic side groups to amplify antagonistic properties, especially through interactions with TRP504.

Interaction fraction analyses of the pilocarpine-hM3R and tiotropium-hM3R complexes were also conducted. An interaction fraction value of 1.0 implies that the ligand contacted a binding residue throughout 100% of the simulation time. Values exceeding 1.0 indicate multiple binding forces. For instance, in the interaction fraction analysis of pilocarpine-hM3R, TYR149 exhibited an interaction fraction value of 1.498, engaging with pilocarpine through 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 contrast, the top three contributors to the binding interaction for tiotropium were ASN508, TRP504, and TYR507 with scores of 1.755, 1.477, and 1.083, respectively. Notably, TRP504 exhibited an interaction fraction value with tiotropium nearly three times higher (1.477) than with pilocarpine (0.591). Furthermore, the top three binding residues with pilocarpine were upstream at positions 148-152, whereas they were downstream at positions 504-508 for tiotropium. Ligand-protein contact analysis of pilocarpine-hM3R revealed ASP148, TYR149, and SER152 intimately associated with the nitrogen ring of pilocarpine. Meanwhile, the ligand-protein contact analysis of tiotropium-hM3R showcased the preservation of π-π stacking exchange via TRP504 during the entirety of the simulation, aligning with the results of molecular docking. Together, these findings emphasize the importance of a positively charged tertiary amine for binding to human M3R and suggest that aromatic side groups may enhance antagonistic properties, particularly through π-π interaction with TRP504.

Finally, for a comparative assessment of the protein-ligand complex stability between pilocarpine and tiotropium, the MM/GBSA method was employed to calculate the binding free energy (ΔG_bind) (Supplementary Table S1). The MM/GBSA analysis indicated that the pilocarpine-hM3R and tiotropium-hM3R complexes exhibited predicted binding free energies of −46 kcal/mol and −96 kcal/mol, respectively. Remarkably, tiotropium demonstrated a more stable complex with hM3R than pilocarpine. Moreover, the nearly two-fold change in binding free energy suggests that once tiotropium binds to hM3R, its binding mode remains stable, and the transmembrane domain of human M3R containing the binding pocket likely maintains a stable, rigid conformation, aligning with the observations from the video analysis of the tiotropium-hM3R complex (Supplementary Video S2).

3.4 High-throughput virtual screening of compound libraries against human M3R homology model

The binding modes of the top three scoring known M3R agonists included the presence of cation-π exchange between a positively charged tertiary amine group and aromatic residues of human M3R, including TYR149, TRP504, and TYR507 (Figure 6). Additionally, hydrogen bonding between an oxygen group and ASN508 was another notable force. Similarly, the binding modes of the top three scoring known M3R antagonists exhibited cation-π exchange between a positively charged tertiary amine and residues TYR149, TRP504, and TYR507; ASN508 was also seen engaging through a hydrogen bond network. Notably, both M3R agonists and antagonists demonstrated a hydrogen bond network with ASN508, aligning with literature emphasizing the significance of this interaction across various muscarinic receptor subtypes (Kruse et al., 2013). However, in contrast to M3R agonists, the antagonists shared a consistent presence of π-π stacking exchange between their aromatic rings and TRP504 (Table 1). Notably, this interaction was absent in all tested hM3R agonists, as all of them lack aromatic side groups in their chemical structures except for pilocarpine. Meanwhile, 12 out of 15 M3R antagonists displayed π-π stacking exchange. Lastly, it is intriguing to observe that M3R antagonists overall achieved higher predicted affinities than agonists, likely attributed to the additional stability provided by π-π stacking exchange.

Figure 6

Figure 6. Top three scoring M3R agonists and antagonists upon molecular docking of six agonists and 15 antagonists of hM3R. (A) 3D binding modes of top three scoring M3R agonists and antagonists within human M3R binding site. Yellow dotted lines indicate hydrogen bond network; purple dotted lines salt bridges; green dotted lines cation-π interactions; blue dotted lines π-π stacking; and orange lines potential steric clash. Note that all three of them involve presence of extensive π-cation and salt bridge interactions with the nitrogen ring via ASP148, TYR149, TRP504, TYR507, and TYR530. (B) 2D visualization of relevant binding residues for the top three scoring M3R agonists and antagonists within human M3R binding site. In addition to the cation-π and salt bridge interactions found with the agonists, M3R antagonists all involved presence of extensive π-π stacking with TRP504, which may explain their improved docking scores compared to agonists.

Table 1

Table 1. High-throughput virtual screening of human M3R agonists (Casson, 2022) and antagonists (Goel et al., 2010).

Subsequently, high-throughput virtual screening was executed on human M3R, utilizing 1,667 FDA-approved or investigational drugs sourced from the ZINC database (Sterling and Irwin, 2015). After three rounds of molecular docking with Glide, the pool of drugs was narrowed down from 1,667 to four potential ligands for human M3R. The four lead compounds identified were edrophonium, rivastigmine, sulpiride, and profenamine (Figure 7). All four lead compounds exhibited extensive cation-π interactions involving their tertiary amine group and human M3R residues TYR149, TRP504, and TYR530. This aligns with our earlier findings, emphasizing the crucial role of the cation-π network in protein-ligand interactions, irrespective of agonism or antagonism. Notably, sulpiride and profenamine displayed additional π-π interactions between their aromatic rings and human M3R’s TRP504. Furthermore, sulpiride and profenamine demonstrated higher predicted affinities compared to edrophonium and rivastigmine. This is consistent with our previous observations that human M3R antagonists generally achieve higher predicted affinities than agonists, likely due to the added stability provided by π-π stacking exchange.

Figure 7

Figure 7. Potential human M3R agonists and antagonists upon high-throughput virtual screening of 1,667 FDA-approved or investigational drugs against human M3R downloaded from ZINC database. 3D analysis of binding interactions is accompanied by their associated 2D ligand interaction diagram for each hit compound. Recall that within the 3D analysis, yellow dotted lines indicate hydrogen bond network; purple dotted lines salt bridges; green dotted lines cation-π interactions; blue dotted lines π-π stacking; and orange lines potential steric clash. Note that both edrophonium and rivastigmine involve presence of cation-π with the tertiary nitrogen via TYR149, TRP504, and TYR530. Furthermore, their predicted affinities were comparable to the docking score of pilocarpine (−9.8 kcal/mol). Meanwhile, binding mode of potential antagonists of human M3R revealed the presence of π-π stacking with TRP504 in addition to hydrogen bonding and salt bridge interactions seen in M3R agonists.

3.5 Assessment of lead compounds for IOP reduction in mice

Rivastigmine and edrophonium were formulated to assess their impact on IOP in mice after topical dosing (Figure 8). Our positive controls for IOP reduction, pilocarpine (a known agonist of human M3R) and physostigmine (a recognized acetylcholinesterase inhibitor) achieved peak IOP reduction approximately 2 h after topical application. Remarkably, physostigmine demonstrated a substantial 30% reduction in IOP, compared to pilocarpine’s 22% reduction (Table 2). However, physostigmine-treated mice exhibited significant cholinergic adverse effects, such as excessive salivation, urination, defecation, muscle paralysis, and hypoventilation. In contrast, mice treated with pilocarpine displayed only mild cholinergic adverse effects, like profuse salivation and urination. For our lead compounds, both rivastigmine and edrophonium achieved maximum IOP reduction approximately 2.8 h after topical application, with effects lasting longer than pilocarpine or physostigmine. The T_end of Rivastigmine was at 5 h, and that of edrophonium was at 5.5 h. Both compounds achieved roughly 21% reduction in IOP, comparable to pilocarpine but less than physostigmine. In terms of cholinergic adverse effects, mice treated with rivastigmine displayed similar effects to pilocarpine, notably excessive salivation, and urination. Intriguingly, mice treated with edrophonium exhibited virtually no adverse cholinergic side effects despite achieving a similar level of IOP reduction as pilocarpine and rivastigmine. It is also noteworthy that AUC values for rivastigmine and edrophonium were higher than pilocarpine, likely due to their longer duration of action, despite similar levels of maximum IOP reduction. Meanwhile, our negative control, atropine, had no impact on IOP reduction. Similarly, our anticipated M3R antagonist, sulpiride, did not induce IOP reduction. To our knowledge, this is the first study to showcase not only the capacity of rivastigmine and edrophonium to lower IOP in mice but also to underscore the diverse potency and toxicity profiles among cholinergic drugs upon topical application. Together, these findings emphasize that the cholinergic activation of the anterior segment of the eye for IOP reduction is more nuanced than previously perceived.

Figure 8

Figure 8. Assessment of pilocarpine, physostigmine, rivastigmine, edrophonium, sulpiride, and atropine for IOP reduction in mice. Pilocarpine and physostigmine reached their maximum IOP reduction at 2 h post-topical application. Notably, physostigmine exhibited an impressive 30% reduction in IOP, surpassing pilocarpine’s 22% reduction. Meanwhile, rivastigmine and edrophonium achieved approximately 21% IOP reduction, which is comparable to pilocarpine but less than physostigmine. The T_end of Rivastigmine was at 5 h, while that of edrophonium was at 5.5 h, suggesting longer durations of action than pilocarpine’s T_end of 4.4 h. Lastly, AUC values for rivastigmine and edrophonium exceeded that of pilocarpine, further suggesting prolonged action despite similar levels of IOP reduction among these three drugs. Diverse potency and toxicity profiles were observed among the investigated cholinergic drugs, suggesting that the cholinergic activation of the anterior segment of the eye for IOP reduction is more nuanced than previously understood.

Table 2

Table 2. Pharmacodynamic parameters after application of drops containing 0.4% pilocarpine, atropine, physostigmine, rivastigmine, and edrophonium to mouse eyes.

3.6 Quantitative structure-activity relationships (QSAR) and ADME (absorption, distribution, metabolism, excretion) analysis of cholinergic drugs

Quantitative structure-activity relationships (QSAR) analysis was employed to predict the chemical properties of the investigated drugs at pH 7.4, corresponding to the pH of the tear film (Table 3). The calculation of cLogD_pH7.4, representing the distribution coefficient of a drug at pH 7.4, indicated that rivastigmine exhibited the highest distribution coefficient, suggesting its predominantly hydrophobic nature at tear film pH. Conversely, edrophonium chloride demonstrated the lowest distribution coefficient, indicating its predominance as a hydrophilic compound at tear film pH. Subsequent analysis of ΔG_o/w, representing the change in Gibbs free energy of drug distribution/partitioning, revealed that rivastigmine had the lowest free energy, aligning with its hydrophobic nature as suggested by its cLogD_pH7.4. Conversely, edrophonium chloride exhibited the highest Gibbs free energy, consistent with its hydrophilic nature at pH 7.4 as indicated by its cLogD_pH7.4. Meanwhile, sulpiride demonstrated the highest topological polar surface area (TPSA), while edrophonium chloride exhibited the lowest TPSA. Importantly, all tested drugs exhibited values of cLogD_pH7.4, ΔG_o/w, and TPSA that are favorable for crossing the corneal membrane, in accordance with Karami’s rule of thumb for optimal ocular permeability and absorption for topical eye drops (Karami et al., 2022). Notably, edrophonium chloride demonstrated the lowest distribution coefficient, highest Gibbs free energy of distribution/partitioning, and lowest TPSA among the tested drugs, all of which may be associated with its reduced cholinergic adverse effects observed in mice compared to rivastigmine, pilocarpine, and physostigmine.

Table 3

Table 3. Quantitative structure-activity relationships (QSAR) analysis of cholinergic drugs and their physiochemical properties related to rule of thumb for optimal ocular permeability and absorption for topical eye drops proposed by Karami et al.

SWISS-ADME was employed to better understand the pharmacokinetics of the cholinergic drugs in relation to their biological activities observed in mice (Table 4). Pilocarpine exhibited moderate lipophilicity (LogP 1.13) and good solubility (LogS −1.38), enabling high gastrointestinal absorption and effective blood-brain barrier penetration without CYP450 interactions or P-glycoprotein substrate activity. Physostigmine demonstrated moderate lipophilicity (LogP 1.65) and lower solubility (LogS −2.64), but maintained high GI absorption and BBB permeability, with potential CYP2D6 inhibition and P-gp substrate status suggesting potential drug interactions and complex elimination route. Rivastigmine, with higher lipophilicity (LogP 2.34) and lower solubility (LogS −2.82), achieved high GI absorption and BBB permeability without CYP interactions or P-gp substrate involvement, indicating favorable pharmacokinetic properties. On the other hand, edrophonium was the only ligand to be characterized by low lipophilicity (LogP −1.07) and low-to-moderate solubility (LogS 0.28). It also had high GI absorption and BBB penetration without CYP interactions, but its P-gp substrate status may influence absorption and excretion. Sulpiride displayed low lipophilicity (LogP 0.83) and solubility (LogS −2.66) yet maintained high GI absorption. It is predicted to have low BBB permeability with no CYP interactions, though its P-gp substrate status may affect its absorption and elimination. Lastly, atropine, with higher lipophilicity (LogP 2.09) and lower solubility (LogS −2.70), also achieved high GI absorption and BBB permeability, but its CYP2D6 inhibition and P-gp substrate status suggest potential drug interactions and complex distribution.

Table 4

Table 4. ADME (Absorption, Distribution, Metabolism, Excretion) Prediction and Analysis of Tested Cholinergic Drugs using SWISS-ADME (Daina et al., 2017). LogP represents the partition coefficient, indicating the compound’s lipophilicity and ability to cross cell membranes, while LogS denotes its aqueous solubility, impacting absorption and bioavailability. GI absorption assesses the probability of compound being absorbed through the gastrointestinal tract. CYP interactions highlight the compound’s interactions with cytochrome P450 enzymes, crucial for its metabolism. BBB permeability measures the compound’s ability to cross the blood-brain barrier, and P-gp substrate status identifies whether the compound is a substrate for P-glycoprotein, a transporter protein that can affect its absorption and excretion.

4 Discussion

The correlation matrix analyses presented in this study shed light on the intricate relationship between muscarinic receptor mRNA expression and IOP dynamics within the anterior segment of the eye. Our findings revealed a prevailing negative correlation pattern between muscarinic receptor mRNA expression and IOP. Notably, this negative correlation reached statistical significance exclusively with Chrm3 mRNA expression, corresponding to the M3 subtype of muscarinic receptor. The specificity of the negative correlation with Chrm3 mRNA underscores the pivotal role of the M3 subtype in modulating IOP dynamics, providing further support for the targeted investigation of M3R in the context of IOP regulation. Furthermore, a negative correlation between Chrm1 and IOP was observed, although it was not reported to be statistically significant (p-value = 0.0647). Note that the M1 receptor is the second most common muscarinic subtype expressed in the ciliary muscle and iris sphincter within the anterior segment of the eye (Gil et al., 1997; Nietgen et al., 1999). This suggests a compelling possibility that the M1 receptor may also play a role in regulating aqueous humor dynamics, alongside the M3 receptor. Another possibility is that Chrm1 is known for its expression in the CNS, including the retina. Recall that the systems genetics analysis was conducted using BXD mice aged 1-30 months, with some subtypes, such as BXD50 and BXD51, known to develop glaucoma spontaneously (Pilkinton et al., 2023). Consequently, optic neuropathy resulting from elevated IOP may lead to reduced CNS expression of Chrm1 in BXD mice and the subsequent negative correlation observed in systems genetics analysis of their whole eyes. In the meantime, an intriguing positive correlation pattern emerged when examining the relationship between Chrm3 and the mRNA expressions of other muscarinic acetylcholine receptors. Together, these findings signify a potential multifaceted regulatory network among muscarinic receptor subtypes that could contribute to the modulation of IOP.

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 lack aromatic side groups. Similar to pilocarpine, all the other agonists shared presence of the cation-π exchange. There were more variations in the presence of hydrogen bond network with ASN508. For instance, NGX267 formed only one hydrogen bond with ASN508, as opposed to pilocarpine forming two, likely because NGX267 has only one oxygen molecule from the carbonyl group. Cevimeline did not form any interaction likely because its oxygen ring was oriented away from ASN508. Curiously, we discovered a unique interaction in the binding mode of Xanomeline that has an aromatic sulfur ring that forms π-π stacking exchange with TRP200. 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 top 10 high-affinity antagonists with docking scores < −10 kcal/mol, oxyphencyclimine was the only molecule containing an aromatic side group that did not form π-π stacking exchange with TRP504. Finally, it is important to address the limitation of molecular docking. Docking scores are mere estimates of the binding free energy of a protein-ligand complex. Ligand binding does not necessarily mean activation or inhibition, as protein activity is governed by a myriad of factors, such as temperature, pH, solvation, cofactor binding, and chaperones (Deller et al., 2016). Hence, molecular docking is a powerful tool to visualize protein-ligand interaction but may not be used to accurately predict experimental Ki or binding affinity (Ramirez and Caballero, 2016; Pantsar and Poso, 2018).

To assess M3R selectivity, we compared our findings with published structures of other muscarinic receptor subtypes bound to ligands. In human M1R, where HTL9936 serves as a partial agonist, cation-π exchange involving a tertiary amine and aromatic residues (TYR106, TRP378, and TYR381) was identified (PDB: 6ZG4) (Brown et al., 2021). This is analogous to the cation-π interactions observed between a ligand’s tertiary amine group and human M3R residues TYR149, TRP504, and TYR530. Hydrogen bonding with ASN382 in M1R mirrored the hydrogen bonding seen with ASN508 in human M3R. However, a notable divergence was the replacement of GLN110 in M1R by ASN152 in human M3R, allowing for a more expansive conformational space in the M3R binding pocket. In human M2R, an agonist interaction involves PHE182 (PDB: 4MQS) (Kruse et al., 2013). The equivalent phenylalanine residue in human M3R (PHE202) had a different orientation, allowing for a broader conformational space exploited in the design of the hM3R-selective small molecule, BS46 (Liu et al., 2018). Human M4R featured an aromatic cage formed by TYR113, TYR278, TYR416, TYR439, and TYR443, enhancing the rigidity of its binding pocket (PDB: 7TRK) (Thal et al., 2016). The binding mode of iperoxo, a super-agonist for human M4R, revealed cation-π interaction between the ligand’s tertiary amine group and TYR278 in human M4R, which is similar to that seen between a ligand’s tertiary amine group and human M3R’s TYR149 (Vuckovic et al., 2023). Lastly, the crystal structure of human M5R exhibited a substitution of GLY408 in human M3R by a charged residue, LYS528. Nevertheless, the binding mode of tiotropium in complex with human M5R closely parallels that of human M3R, featuring cation-π interaction with TYR417 as well as a double hydrogen bond network involving ASN395 and tiotropium’s hydroxyl and ester groups (PDB: 6OL9) (Vuckovic etal., 2019).

Our in silico assays highlighted the significance of a tertiary amine with a positively charged nitrogen molecule and the ability to hydrogen bond with ASN508 for M3R binding, irrespective of agonism or antagonism. This aligns with acetylcholine’s structure, an endogenous M3R ligand featuring a positively charged tertiary amine and a carboxylic acid group capable of hydrogen bonding with ASN508. For agonism promotion, it is advisable that the ligand lacks aromatic side groups and cannot engage in π-π stacking exchange with TRP504. These structural insights led to identifying two cholinergic drugs from a pool of 1,667 FDA-approved or investigational drugs demonstrating the most favorable profile for reducing IOP: edrophonium and rivastigmine. Edrophonium, an FDA-approved acetylcholinesterase inhibitor for myasthenia gravis diagnosis, cardiac arrhythmia management, and neuromuscular blocking agent toxicity reversal, induces mild cholinergic effects, including bradycardia and bronchoconstriction (Pakala et al., 2023). In the past, edrophonium was tested in myasthenia gravis and normal patients to assess for its impact on IOP reduction. However, the study delivered the drug through systemic injection, involved only five patients, and produced inconclusive results, thus limiting the clinical utility of the data (Kornblueth et al., 1960). Otherwise, as of February of 2024, there are no studies published in PubMed, Scopus, and Embase databases testing its efficacy on lowering IOP in any other experimental models. Rivastigmine, another FDA-approved acetylcholinesterase inhibitor for neurological symptoms in Alzheimer’s and Parkinson’s disease, triggers cholinergic hyperactivation-related adverse effects, like nausea and vomiting (Varadharajan et al., 2023). Significantly, the efficacy of reducing IOP using rivastigmine was examined in rabbit eyes via topical formulation, resulting in a 15.2% maximal IOP reduction with the application of 1% rivastigmine (Goldblum et al., 2000). However, the study was published as a conference paper with incomplete information on specific topical formulation utilized for the study and did not investigate systemic cholinergic adverse effects in rabbits. Otherwise, as of February of 2024, there are no studies published in PubMed, Scopus, and Embase databases demonstrating its efficacy on lowering IOP in any other experimental models.

Most notably, as of June 2023, the FDA withdrew its approvals of majority of acetylcholinesterase inhibitors, including both edrophonium and rivastigmine, in the US due to reported concerns over their adverse systemic effects, specifically QT prolongation and Torsade de pointes (Verdin, 2023). It is essential to highlight that despite these withdrawals, edrophonium and rivastigmine still possess the potential to indirectly reduce IOP by inhibiting acetylcholine degradation through acetylcholinesterase inhibition. Moreover, our in silico assays demonstrated high-affinity interactions, revealing that both edrophonium and rivastigmine favor agonist binding modes within the M3R binding site. Consequently, in the pursuit of repurposing existing cholinergic drugs, edrophonium and rivastigmine emerged as the most promising candidates for lowering IOP through manipulation of the parasympathetic nervous system of the eye. This discovery prompted a comparative analysis of their IOP-lowering effects and cholinergic toxicities in mice.

The outcomes of our in vivo investigation into the impact of topically administered rivastigmine and edrophonium on IOP highlighted their potential as safer cholinergic drugs for IOP reduction, surpassing the established FDA-approved pilocarpine. Both rivastigmine and edrophonium demonstrated a comparable level of IOP reduction to pilocarpine, yet their effects exhibited a prolonged duration. This extended efficacy suggests the possibility of reduced dosing frequency in glaucoma treatment, enhancing patient adherence and therapeutic outcomes. Pilocarpine IOP drops are associated with cholinergic adverse effects, including increased sweating, muscle tremors, blurry vision, and nausea/vomiting (Pakala et al., 2023; Panarese and Moshirfar, 2023). Despite efforts to develop topical formulations to mitigate its cholinergic adverse effects, the washout of drugs into the tear duct remains a challenge. This drainage route can lead to drug entry into the nasopharynx and eventually the gastrointestinal system, causing systemic side effects over time (Akhter et al., 2022). While rivastigmine displayed adverse cholinergic effects similar to pilocarpine, a remarkable observation was made with edrophonium, which exhibited virtually no cholinergic side effects. Recall that IOP study utilized a single topical drop of each drug and monitored IOP changes in the subsequent 8 h. Unfortunately, we did not prepare enough solution or chemical for a topical hydrogel formulation towards multiple doses over an extended period to assess their long-term systemic side effects. While there are studies exploring the systemic adverse effects of oral cholinergic drugs, there were scarce literature reporting systemic adverse effects associated with repeat doses of topical cholinergic drops. Since adverse cholinergic effects significantly constrain the clinical use of cholinergic drugs, it is imperative to conduct future longitudinal studies investigating the effects of multiple doses of these drugs, before we further test their IOP-reducing effects in other experimental models of glaucoma.

One potential limitation of our findings lies in the fact that both edrophonium and rivastigmine are acknowledged inhibitors of acetylcholinesterase. While our modeling studies suggest that these compounds favor M3R agonism, determining whether the observed IOP reduction is directly attributed to M3R agonism, indirectly through diminished acetylcholine degradation, or a combination of both is challenging. Notably, pilocarpine, our positive control, primarily reduces IOP through M3R agonism, while physostigmine, another positive control, achieves IOP reduction by inhibiting acetylcholinesterase. However, in animal studies, both edrophonium and rivastigmine exhibited pharmacodynamic profiles distinct from our positive controls, indicating a more nuanced mechanism beyond M3R agonism, acetylcholinesterase inhibition, or a combination of both. Additionally, sulpiride, a dopamine D2 receptor antagonist used in schizophrenia management, also inhibits acetylcholinesterase (Dross, 1977). Modeling studies revealed π-π stacking with TRP504, favoring M3R antagonism over agonism. Animal studies confirmed that sulpiride, like atropine (a known muscarinic antagonist), failed to reduce IOP in mice, suggesting two possibilities. The first possibility stems from the results of our modelling that sulpiride’s predicted M3R antagonism may have overridden the increased acetylcholine levels resulting from its inhibition of acetylcholinesterase. Secondly, sulpiride, known for its diverse off-target effects on P-glycoproteins, carbonic anhydrases, and ATP-binding cassette sub-family G member 2 (ABCG2), may introduce complexities confounding its role in IOP modulation (Sulpiride, 2006). Cumulatively, all these findings suggest that the cholinergic modulation of the anterior eye segment for IOP reduction is more intricate than previously understood. Consequently, our team conducted additional in silico assays on the investigated cholinergic drugs through QSAR and ADME analysis.

The QSAR results highlighted that edrophonium chloride demonstrated the smallest TPSA among the tested drugs. This characteristic is favorable for optimal corneal permeability, suggesting that a significant portion of topically applied edrophonium likely traversed the corneal membrane rather than being drained through the nasolacrimal duct into systemic circulation. This phenomenon could potentially account for its diminished systemic cholinergic adverse effects observed in mice upon topical application. Meanwhile, sulpiride was identified as having the highest TPSA, potentially leading to decreased corneal permeability, providing a third explanation for its inability to lower IOP when applied topically. Notably, both QSAR and ADME analyses identified edrophonium as the most water-soluble compound among the tested drugs. This offers another potential explanation to why mice treated with edrophonium exhibited reduced systemic cholinergic adverse effects compared to other drugs: its aqueous solubility facilitated its entrapment in the aqueous core of our Carbopol 981 hydrogel formulation, allowing for controlled release with reduced washouts through drainage routes on the ocular surface.

Together, our findings from docking, affinity studies, QSAR analysis, and ADME studies collectively support the significance of specific molecular features associated with cholinergic activation as well as their optimal pharmacokinetic properties. Docking and affinity studies demonstrated that a tertiary amine with a positively charged nitrogen and the ability to hydrogen bond with ASN508 are critical for M3R binding, irrespective of agonism or antagonism. These structural insights led to the identification of edrophonium and rivastigmine as promising drugs for reducing IOP among 1,667 FDA-approved or investigational drugs. In vivo studies confirmed their comparable efficacy in reducing IOP to pilocarpine, with edrophonium showing prolonged effects and minimal systemic cholinergic adverse effects. The QSAR analysis further highlighted edrophonium’s small TPSA as critical descriptor influencing drug permeability across biological membranes. This characteristic is advantageous for ocular drug delivery, enhancing local bioavailability while minimizing systemic absorption. Additionally, both QSAR and ADME analyses identified edrophonium as the most water-soluble compound among the drugs tested. Once delivered directly to the corneal epithelium, edrophonium’s solubility may contribute to its entrapment in the aqueous core of our hydrogel vehicle, allowing for sustained release and reduced entry to systemic circulation through drainage routes on the ocular surface. These findings collectively support the potential of edrophonium to guide the development of derivatives or new cholinergic drugs that can safely lower IOP in glaucoma with reduced cholinergic adverse effects.

Overall, our discoveries indicate that the cholinergic activation of the anterior eye segment for IOP modulation involves more complexity than simple acetylcholinesterase inhibition or M3 receptor agonism. These insights prompt a re-evaluation of the parasympathetic system in the anterior eye segment for the development of safe and effective treatments for glaucoma. In terms of clinical implications, the combination of reduced systemic side effects and an extended duration of action positions edrophonium and its derivatives as attractive candidates for further investigation in experimental animal models of glaucoma.

5 Conclusion

In conclusion, our investigation into the molecular interactions, binding affinities, and in vivo effects of cholinergic drugs, specifically rivastigmine and edrophonium, has provided valuable insights into their potentials for IOP reduction with distinct pharmacokinetic profiles when compared to pilocarpine. While edrophonium has demonstrated potential as an IOP-reducing agent with fewer cholinergic adverse effects in mice, it is imperative to recognize its discontinuation in the United States and the reported adverse effects in patients that prompted this decision. Keeping this in mind, future research endeavors can utilize these insights to facilitate the clinical translation of edrophonium’s IOP-lowering effects, potentially leading to the development of its derivatives that may expand therapeutic options against glaucoma. This also underscores the potential need for re-opening the investigation of the cholinergic system in the eye, aiming to develop novel drugs capable of lowering IOP with a safer toxicity profile.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Ethics statement

All procedures involving mice were approved by the Animal Care and Use review board of the University of Tennessee Health Science Center (UTHSC) and followed the Association of Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research in addition to the guidelines for laboratory animal experiments (Institute of Laboratory Animal Resources, Public Health Service Policy on Humane Care and Use of Laboratory Animals. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

MK: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Validation, Visualization, Writing–original draft, Writing–review and editing. MM: Formal Analysis, Methodology, Resources, Visualization, Writing–review and editing. MJ: Conceptualization, Funding acquisition, Project administration, Resources, Software, Supervision, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. Funding was provided by a Challenge Grant from Research to Prevent Blindness to the Hamilton Eye Institute and R01EY021200 to MJ.

Acknowledgments

We extend our gratitude to Kirk Hevener and Wei Li for providing access to the molecular modeling station at the College of Pharmacy, University of Tennessee Health Science Center. We also extend our gratitude to Robert Williams and Lu Lu at the Department of Genetics, Genomics and Informatics, University of Tennessee Health Science Center. Funding was provided by a Challenge Grant from Research to Prevent Blindness to the Hamilton Eye Institute and R01EY021200 to MJ. Figures 1–4 were generated using BioRender.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2024.1423351/full#supplementary-material

SUPPLEMENTARY FIGURE S1The full amino acid sequence of human M3 muscarinic receptor (M3R) is presented with its major tertiary domains color-coded for clarity. Molecular weight (MW) and isoelectric point (pI) were calculated from this sequence using the SWISS-EXPASY pI/MW tool. It is important to note that the homology model used in our studies lacks intracellular loop 3 (indicated by a green box), which is absent from the structural representation. However, this omission is not expected to influence the interaction dynamics at the binding site, primarily located within the transmembrane domain.

SUPPLEMENTARY FIGURE S2Comparison of the binding mode of Iperexo (a muscarinic super-agonist) and that of pilocarpine. 3D binding modes of iperexo, pilocarpine, and their overlay within human M3R binding site were presented on the top row. Yellow dotted lines indicate hydrogen bond network; purple dotted lines salt bridges; green dotted lines cation-π interactions. Note that both iperexo and pilocarpine have similar conformations with their positively charged tertiary amine oriented towards ASP148 and their aromatic rings towards VAL156 at the other end. 2D protein-ligand interaction diagrams of iperexo and pilocarpine were presented on the bottom row. The binding modes of both compounds involved presence of extensive π-cation and salt bridge interactions with the nitrogen ring via ASP148, TYR149, and TYR507. Also note the presence of hydrogen bonding between pilocarpine’s carbonyl ring and human M3R’s TRP200. There was no sign of π-π stacking interaction with TRP504.

SUPPLEMENTARY FIGURE S3Ramachandran plot analysis of template M3R (PDB: 8E9Z) in comparison to human M3R homology model. Each dot on the plot represents a peptide with corresponding phi-psi angles. Blue dots indicate “allowable” peptides with favorable phi-psi angles, while red dots indicate “disallowed” peptides due to unfavorable phi-psi angles. Peptides with phi/psi angles around -90/120 degrees correspond to beta sheets, -90/-30 degrees to right-handed alpha helices, and 90/60 degrees to left-handed alpha helices. Overall, the analysis revealed a similar structural motif between the template and the homology model. The template had only 2 “disallowed” peptides, compared to 30 “disallowed” peptides in the homology model.

SUPPLEMENTARY VIDEO S1Video analysis of the 100 ns MD simulation result for hM3R-pilocarpine complex, where hM3R is depicted as green cartoon structure and pilocarpine as cyan stick structure.

SUPPLEMENTARY VIDEO S2Video analysis of the 100 ns MD simulation result for hM3R-tiotropium complex, where hM3R is depicted as orange cartoon structure and tiotropium as cyan stick structure.

SUPPLEMENTARY TABLE S1MM/GBSA binding free energy calculation of hM3R-pilocarpine and hM3R-tiotropium complexes (N = 3). All values are in kcal/mol. Standard deviations less than 0.01 kcal/mol are not depicted.

References

Akhter, M. H., Ahmad, I., Alshahrani, M. Y., Al-Harbi, A. I., Khalilullah, H., Afzal, O., et al. (2022). Drug delivery challenges and current progress in nanocarrier-based ocular therapeutic system. Gels. 8 (2), 82. doi:10.3390/gels8020082

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, B. D., Higuchi, W. I., and Raykar, P. V. (1988). Heterogeneity effects on permeability-partition coefficient relationships in human stratum corneum. Pharm. Res. 5, 566–573. doi:10.1023/a:1015989929342

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrade, O. A., and Zafar Gondal, A. (2023). Physostigmine. Treasure Island, FL: StatPearls.

Google Scholar

Andreux, P. A., Williams, E. G., Koutnikova, H., Houtkooper, R. H., Champy, M. F., Henry, H., et al. (2012). Systems genetics of metabolism: the use of the BXD murine reference panel for multiscalar integration of traits. Cell. 150 (6), 1287–1299. doi:10.1016/j.cell.2012.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Arens, A. M., and Kearney, T. (2019). Adverse effects of physostigmine. J. Med. Toxicol. 15 (3), 184–191. doi:10.1007/s13181-019-00697-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashbrook, D. G., Gini, B., and Hager, R. (2015). Genetic variation in offspring indirectly influences the quality of maternal behaviour in mice. Elife 4, e11814. doi:10.7554/eLife.11814

PubMed Abstract | CrossRef Full Text | Google Scholar

Audain, E., Ramos, Y., Hermjakob, H., Flower, D. R., and Perez-Riverol, Y. (2016). Accurate estimation of isoelectric point of protein and peptide based on amino acid sequences. Bioinformatics 32 (6), 821–827. doi:10.1093/bioinformatics/btv674

PubMed Abstract | CrossRef Full Text | Google Scholar

Balla, H., Borsodi, K., Orsy, P., Horvath, B., Molnar, P. J., Lenart, A., et al. (2023). Intracellular signaling pathways of muscarinic acetylcholine receptor-mediated detrusor muscle contractions. Am. J. Physiol. Ren. Physiol. 325 (5), F618–F628. doi:10.1152/ajprenal.00261.2022

CrossRef Full Text | Google Scholar

Basaran, E., Cakmak, R., Sahin, D., Kopru, S., Turkmenoglu, B., and Akkoc, S. (2024). Design, spectroscopic characterization, in silico and in vitro cytotoxic activity assessment of newly synthesized thymol Schiff base derivatives. J. Biomol. Struct. Dyn. 1–14. doi:10.1080/07391102.2024.2301747

CrossRef Full Text | Google Scholar

Beidoe, G., and Mousa, S. A. (2012). Current primary open-angle glaucoma treatments and future directions. Clin. Ophthalmol. 6, 1699–1707. doi:10.2147/OPTH.S32933

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, A. J. H., Bradley, S. J., Marshall, F. H., Brown, G. A., Bennett, K. A., Brown, J., et al. (2021). From structure to clinic: design of a muscarinic M1 receptor agonist with potential to treatment of Alzheimer's disease. Cell. 184 (24), 5886–5901.e22. doi:10.1016/j.cell.2021.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Casson, R. J. (2022). Medical therapy for glaucoma: a review. Clin. Exp. Ophthalmol. 50 (2), 198–212. doi:10.1111/ceo.13989

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Chen, Z., Deng, C., Chen, W., Zhang, H., and Wang, J. (2023). Changes to outflow structures after pilocarpine in primary open angle glaucoma compared with healthy individuals using optical coherence tomography. J. Glaucoma 32 (7), 593–599. doi:10.1097/IJG.0000000000002165

PubMed Abstract | CrossRef Full Text | Google Scholar

Chesler, E. J., Lu, L., Shou, S., Qu, Y., Gu, J., Wang, J., et al. (2005). Complex trait analysis of gene expression uncovers polygenic and pleiotropic networks that modulate nervous system function. Nat. Genet. 37 (3), 233–242. doi:10.1038/ng1518

PubMed Abstract | CrossRef Full Text | Google Scholar

Daina, A., Michielin, O., and Zoete, V. (2017). SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci. Rep. 7, 42717. doi:10.1038/srep42717

PubMed Abstract | CrossRef Full Text | Google Scholar

Deller, M. C., Kong, L., and Rupp, B. (2016). Protein stability: a crystallographer's perspective. Acta Crystallogr. F. Struct. Biol. Commun. 72 (Pt 2), 72–95. doi:10.1107/S2053230X15024619

PubMed Abstract | CrossRef Full Text | Google Scholar

Dross, K. (1977). The relevance of anticholinesterase properties to toxicity and neuromuscular effects of sulpiride (author's transl). Arzneimittelforschung 27 (2), 404–406.

PubMed Abstract | Google Scholar

Froloff, N., Windemuth, A., and Honig, B. (1997). On the calculation of binding free energies using continuum methods: application to MHC class I protein-peptide interactions. Protein Sci. 6 (6), 1293–1301. doi:10.1002/pro.5560060617

PubMed Abstract | CrossRef Full Text | Google Scholar

Genheden, S., and Ryde, U. (2015). The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin. Drug Discov. 10 (5), 449–461. doi:10.1517/17460441.2015.1032936

PubMed Abstract | CrossRef Full Text | Google Scholar

Gheidari, D., Mehrdad, M., and Bayat, M. (2024). Synthesis, docking, MD simulation, ADMET, drug likeness, and DFT studies of novel furo[2,3-b]indol-3a-ol as promising Cyclin-dependent kinase 2 inhibitors. Sci. Rep. 14 (1), 3084. doi:10.1038/s41598-024-53514-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Gil, D. W., Krauss, H. A., Bogardus, A. M., and Wolde Mussie, E. (1997). Muscarinic receptor subtypes in human iris-ciliary body measured by immunoprecipitation. Invest. Ophthalmol. Vis. Sci. 38 (7), 1434–1442.

PubMed Abstract | Google Scholar

Goel, M., Picciani, R. G., Lee, R. K., and Bhattacharya, S. K. (2010). Aqueous humor dynamics: a review. Open Ophthalmol. J. 4, 52–59. doi:10.2174/1874364101004010052

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldblum, D., Garweg, J. G., and Bohnke, M. (2000). Topical rivastigmine, a selective acetylcholinesterase inhibitor, lowers intraocular pressure in rabbits. J. Ocul. Pharmacol. Ther. 16 (1), 29–35. doi:10.1089/jop.2000.16.29

PubMed Abstract | CrossRef Full Text | Google Scholar

Greig, N. H., Reale, M., and Tata, A. M. (2013). New pharmacological approaches to the cholinergic system: an overview on muscarinic receptor ligands and cholinesterase inhibitors. Recent Pat. CNS Drug Discov. 8 (2), 123–141. doi:10.2174/1574889811308020003

PubMed Abstract | CrossRef Full Text | Google Scholar

Guedes, I. A., de Magalhaes, C. S., and Dardenne, L. E. (2014). Receptor-ligand molecular docking. Biophys. Rev. 6 (1), 75–87. doi:10.1007/s12551-013-0130-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Haga, C., Waller, T., Ahuja, A. S., Farford, B., Dawson, N., and Yin, M. (2022). There's danger in the drops: systemic effects of ophthalmic drops used to treat glaucoma. Cureus 14 (1), e20945. doi:10.7759/cureus.20945

PubMed Abstract | CrossRef Full Text | Google Scholar

Halim, S. A., Waqas, M., Khan, A., and Al-Harrasi, A. (2021). In Silico prediction of novel inhibitors of SARS-CoV-2 main protease through structure-based virtual screening and molecular dynamic simulation. Pharm. (Basel) 14 (9), 896. doi:10.3390/ph14090896

CrossRef Full Text | Google Scholar

Houtkooper, R. H., Mouchiroud, L., Ryu, D., Moullan, N., Katsyuba, E., Knott, G., et al. (2013). Mitonuclear protein imbalance as a conserved longevity mechanism. Nature 497 (7450), 451–457. doi:10.1038/nature12188

PubMed Abstract | CrossRef Full Text | Google Scholar

Karami, T. K., Hailu, S., Feng, S., Graham, R., and Gukasyan, H. J. (2022). Eyes on lipinski's rule of five: a new “rule of thumb” for physicochemical design space of ophthalmic drugs. J. Ocul. Pharmacol. Ther. 38 (1), 43–55. doi:10.1089/jop.2021.0069

CrossRef Full Text | Google Scholar

Keller, K. E., and Peters, D. M. (2022). Pathogenesis of glaucoma: extracellular matrix dysfunction in the trabecular meshwork-A review. Clin. Exp. Ophthalmol. 50 (2), 163–182. doi:10.1111/ceo.14027

CrossRef Full Text | Google Scholar

Kim, M. J., Martin, C. A., Kim, J., and Jablonski, M. M. (2023). Computational methods in glaucoma research: current status and future outlook. Mol. Asp. Med. 94, 101222. doi:10.1016/j.mam.2023.101222

CrossRef Full Text | Google Scholar

Kistemaker, L. E. M., Elzinga, C. R. S., Tautermann, C. S., Pieper, M. P., Seeliger, D., Alikhil, S., et al. (2019). Second M(3) muscarinic receptor binding site contributes to bronchoprotective by tiotropium. Br. J. Pharmacol. 176 (16), 2864–2876. doi:10.1111/bph.14707

PubMed Abstract | CrossRef Full Text | Google Scholar

Knox, C., Wilson, M., Klinger, C. M., Franklin, M., Oler, E., Wilson, A., et al. (2024). DrugBank 6.0: the DrugBank knowledgebase for 2024. Nucleic Acids Res. 52 (D1), D1265–D1275. doi:10.1093/nar/gkad976

PubMed Abstract | CrossRef Full Text | Google Scholar

Kornblueth, W., Jampolsky, A., Tamler, E., and Marg, E. (1960). Contraction of the oculorotary muscles and intraocular pressure. A tonographic and electromyographic study of the effect of edrophonium chloride (tensilon) and succinylcholine (anectine) on the intraocular pressure. Am. J. Ophthalmol. 49, 1381–1387. doi:10.1016/0002-9394(60)91359-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Kruse, A. C., Hu, J., Pan, A. C., Arlow, D. H., Rosenbaum, D. M., Rosemond, E., et al. (2012). Structure and dynamics of the M3 muscarinic acetylcholine receptor. Nature 482 (7386), 552–556. doi:10.1038/nature10867

PubMed Abstract | CrossRef Full Text | Google Scholar

Kruse, A. C., Ring, A. M., Manglik, A., Hu, J., Hu, K., Eitel, K., et al. (2013). Activation and allosteric modulation of a muscarinic acetylcholine receptor. Nature 504 (7478), 101–106. doi:10.1038/nature12735

PubMed Abstract | CrossRef Full Text | Google Scholar

Kudlak, M., and Tadi, P. (2023). Physiology, muscarinic receptor. Treasure Island, FL: StatPearls.

Google Scholar

Kumar, R., Jade, D., and Gupta, D. (2019). A novel identification approach for discovery of 5-HydroxyTriptamine 2A antagonists: combination of 2D/3D similarity screening, molecular docking and molecular dynamics. J. Biomol. Struct. Dyn. 37 (4), 931–943. doi:10.1080/07391102.2018.1444509

PubMed Abstract | CrossRef Full Text | Google Scholar

Kwon, Y. H., Fingert, J. H., Kuehn, M. H., and Alward, W. L. (2009). Primary open-angle glaucoma. N. Engl. J. Med. 360 (11), 1113–1124. doi:10.1056/NEJMra0804630

PubMed Abstract | CrossRef Full Text | Google Scholar

Leung, S. S., Mijalkovic, J., Borrelli, K., and Jacobson, M. P. (2012). Testing physical models of passive membrane permeation. J. Chem. Inf. Model. 52 (6), 1621–1636. doi:10.1021/ci200583t

PubMed Abstract | CrossRef Full Text | Google Scholar

Leung, S. S., Sindhikara, D., and Jacobson, M. P. (2016). Simple predictive models of passive membrane permeability incorporating size-dependent membrane-water partition. J. Chem. Inf. Model. 56 (5), 924–929. doi:10.1021/acs.jcim.6b00005

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Ibrahim, M. M., Chen, H., Li, W., and Jablonski, M. M. (2021). In Silico Screening and In Vivo evaluation of potential CACNA2D1 antagonists as intraocular pressure-reducing agents in glaucoma therapy. Pharm. (Basel) 14 (9), 887. doi:10.3390/ph14090887

CrossRef Full Text | Google Scholar

Li, H., Robertson, A. D., and Jensen, J. H. (2005). Very fast empirical prediction and rationalization of protein pKa values. Proteins 61 (4), 704–721. doi:10.1002/prot.20660

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Abel, R., Zhu, K., Cao, Y., Zhao, S., and Friesner, R. A. (2011). The VSGB 2.0 model: a next generation energy model for high resolution protein structure modeling. Proteins 79 (10), 2794–2812. doi:10.1002/prot.23106

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Zhang, X., Li, S., and Guo, J. (2024). Unraveling the interplay of extracellular domain conformational changes and parathyroid hormone type 1 receptor activation in class B1 G protein-coupled receptors: integrating enhanced sampling molecular dynamics simulations and markov state models. ACS Chem. Neurosci. 15, 844–853. doi:10.1021/acschemneuro.3c00747

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Hofmann, J., Fish, I., Schaake, B., Eitel, K., Bartuschat, A., et al. (2018). Structure-guided development of selective M3 muscarinic acetylcholine receptor antagonists. Proc. Natl. Acad. Sci. U. S. A. 115 (47), 12046–12050. doi:10.1073/pnas.1813988115

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, P., Wang, F., Song, Y., Wang, M., and Zhang, X. (2022). Current situation and progress of drugs for reducing intraocular pressure. Ther. Adv. Chronic Dis. 13, 20406223221140392. doi:10.1177/20406223221140392

PubMed Abstract | CrossRef Full Text | Google Scholar

Llobet, A., Gasull, X., and Gual, A. (2003). Understanding trabecular meshwork physiology: a key to the control of intraocular pressure? News Physiol. Sci. 18, 205–209. doi:10.1152/nips.01443.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez-Archundia, M., Cordomi, A., Garriga, P., and Perez, J. J. (2012). Molecular modeling of the M3 acetylcholine muscarinic receptor and its binding site. J. Biomed. Biotechnol. 2012, 789741. doi:10.1155/2012/789741

PubMed Abstract | CrossRef Full Text | Google Scholar

Moshirfar, M., Parker, L., Birdsong, O. C., Ronquillo, Y. C., Hofstedt, D., Shah, T. J., et al. (2018). Use of rho kinase inhibitors in Ophthalmology: a review of the literature. Med. Hypothesis Discov. Innov. Ophthalmol. 7 (3), 101–111.

PubMed Abstract | Google Scholar

Muhammed, M. T., and Aki-Yalcin, E. (2019). Homology modeling in drug discovery: overview, current applications, and future perspectives. Chem. Biol. Drug Des. 93 (1), 12–20. doi:10.1111/cbdd.13388

PubMed Abstract | CrossRef Full Text | Google Scholar

Mulligan, M. K., Mozhui, K., Prins, P., and Williams, R. W. (2017). GeneNetwork: a toolbox for systems genetics. Methods Mol. Biol. 1488, 75–120. doi:10.1007/978-1-4939-6427-7_4

PubMed Abstract | CrossRef Full Text | Google Scholar

Nietgen, G. W., Schmidt, J., Hesse, L., Honemann, C. W., and Durieux, M. E. (1999). Muscarinic receptor functioning and distribution in the eye: molecular basis and implications for clinical diagnosis and therapy. Eye (Lond). 13 (Pt 3a), 285–300. doi:10.1038/eye.1999.78

PubMed Abstract | CrossRef Full Text | Google Scholar

Pakala, R. S., Brown, K. N., and Preuss, C. V. (2023). Cholinergic medications. Treasure Island, FL: StatPearls.

Google Scholar

Panarese, V., and Moshirfar, M. (2023). Pilocarpine. Treasure Island, FL: StatPearls.

Google Scholar

Pantsar, T., and Poso, A. (2018). Binding affinity via docking: fact and fiction. Molecules 23 (8), 1899. doi:10.3390/molecules23081899

PubMed Abstract | CrossRef Full Text | Google Scholar

Pilkinton, S., Kim, Y., Wang, X., Hollingsworth, T., and Jablonski, M. M. (2023). BXD Murine recombinant models of glaucoma subtypes. Investigative Ophthalmol. Vis. Sci. 64 (8), 4712.

Google Scholar

Ramirez, D., and Caballero, J. (2016). Is it reliable to use common molecular docking methods for comparing the binding affinities of enantiomer pairs for their protein target? Int. J. Mol. Sci. 17 (4), 525. doi:10.3390/ijms17040525

PubMed Abstract | CrossRef Full Text | Google Scholar

Roos, K., Wu, C., Damm, W., Reboul, M., Stevenson, J. M., Lu, C., et al. (2019). OPLS3e: Extending Force Field Coverage for Drug-Like Small Molecules. J. Chem. Theory Comput. 15 (3), 1863–1874. doi:10.1021/acs.jctc.8b01026

PubMed Abstract | CrossRef Full Text | Google Scholar

Roy, P. K., Majumder, R., and Mandal, M. (2024). In-silico identification of novel DDI2 inhibitor in glioblastoma via repurposing FDA approved drugs using molecular docking and MD simulation study. J. Biomol. Struct. Dyn. 42 (5), 2270–2281. doi:10.1080/07391102.2023.2204371

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruan, Y., Patzak, A., Pfeiffer, N., and Gericke, A. (2021). Muscarinic acetylcholine receptors in the retina-therapeutic implications. Int. J. Mol. Sci. 22 (9), 4989. doi:10.3390/ijms22094989

PubMed Abstract | CrossRef Full Text | Google Scholar

Saini, G., Dalal, V., Savita, B. K., Sharma, N., Kumar, P., and Sharma, A. K. (2019). Molecular docking and dynamic approach to virtual screen inhibitors against Esbp of Candidatus Liberibacter asiaticus. J. Mol. Graph Model. 92, 329–340. doi:10.1016/j.jmgm.2019.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Shah, M., Lee, G., Lefebvre, D. R., Kronberg, B., Loomis, S., Brauner, S. C., et al. (2013). A cross-sectional survey of the association between bilateral topical prostaglandin analogue use and ocular adnexal features. PLoS One 8 (5), e61638. doi:10.1371/journal.pone.0061638

PubMed Abstract | CrossRef Full Text | Google Scholar

Shelley, J. C., Cholleti, A., Frye, L. L., Greenwood, J. R., Timlin, M. R., and Uchimaya, M. (2007). Epik: a software program for pK(a) prediction and protonation state generation for drug-like molecules. J. Comput. Aided Mol. Des. 21 (12), 681–691. doi:10.1007/s10822-007-9133-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Stampelou, M., Ladds, G., and Kolocouris, A. (2024). Computational workflow for refining AlphaFold models in drug design using kinetic and thermodynamic binding calculations: a case study for the unresolved inactive human adenosine A(3) receptor. J. Phys. Chem. B 128 (4), 914–936. doi:10.1021/acs.jpcb.3c05986

PubMed Abstract | CrossRef Full Text | Google Scholar

Steinberg, L., Russo, J., and Frey, J. (2019). A new topological descriptor for water network structure. J. Cheminform 11 (1), 48. doi:10.1186/s13321-019-0369-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Sterling, T., and Irwin, J. J. (2015). ZINC 15--ligand discovery for everyone. J. Chem. Inf. Model. 55 (11), 2324–2337. doi:10.1021/acs.jcim.5b00559

PubMed Abstract | CrossRef Full Text | Google Scholar

Sulpiride (2006). Drugs and lactation database (LactMed(R)). MD: Bethesda.

Google Scholar

Tamm, E. R. (2009). The trabecular meshwork outflow pathways: structural and functional aspects. Exp. Eye Res. 88 (4), 648–655. doi:10.1016/j.exer.2009.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, B. A., Wnek, C., Kotlus, B. S., Roemer, N., MacTaggart, T., and Phillips, S. J. (1999). Genotyping new BXD recombinant inbred mouse strains and comparison of BXD and consensus maps. Mamm. Genome 10 (4), 335–348. doi:10.1007/s003359900998

PubMed Abstract | CrossRef Full Text | Google Scholar

Thal, D. M., Sun, B., Feng, D., Nawaratne, V., Leach, K., Felder, C. C., et al. (2016). Crystal structures of the M1 and M4 muscarinic acetylcholine receptors. Nature 531 (7594), 335–340. doi:10.1038/nature17188

PubMed Abstract | CrossRef Full Text | Google Scholar

Thul, P. J., and Lindskog, C. (2018). The human protein atlas: a spatial map of the human proteome. Protein Sci. 27 (1), 233–244. doi:10.1002/pro.3307

PubMed Abstract | CrossRef Full Text | Google Scholar

Vacca, C., Maione, S., Tozzi, A., and Marmo, E. (1987). Benzamides and cholinesterases. Res. Commun. Chem. Pathol. Pharmacol. 55 (2), 193–201.

PubMed Abstract | Google Scholar

Varadharajan, A., Davis, A. D., Ghosh, A., Jagtap, T., Xavier, A., Menon, A. J., et al. (2023). Guidelines for pharmacotherapy in Alzheimer's disease - a primer on FDA-approved drugs. J. Neurosci. Rural. Pract. 14 (4), 566–573. doi:10.25259/JNRP_356_2023

PubMed Abstract | CrossRef Full Text | Google Scholar

Verdin, P. (2023). FDA new drug approvals in Q3 2023. Nat. Rev. Drug Discov. 22 (12), 946. doi:10.1038/d41573-023-00184-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidal, M., Onderdijk, K. E., Aguilera, A. M., Six, J., Maes, P. J., Fritz, T. H., et al. (2023). Cholinergic-related pupil activity reflects level of emotionality during motor performance. Eur. J. Neurosci. 59, 2193–2207. doi:10.1111/ejn.15998

PubMed Abstract | CrossRef Full Text | Google Scholar

Vuckovic, Z., Gentry, P. R., Berizzi, A. E., Hirata, K., Varghese, S., Thompson, G., et al. (2019). Crystal structure of the M5 muscarinic acetylcholine receptor. Proc. Natl. Acad. Sci. U. S. A. 116 (51), 26001–26007. doi:10.1073/pnas.1914446116

PubMed Abstract | CrossRef Full Text | Google Scholar

Vuckovic, Z., Wang, J., Pham, V., Mobbs, J. I., Belousoff, M. J., Bhattarai, A., et al. (2023). Pharmacological hallmarks of allostery at the M4 muscarinic receptor elucidated through structure and dynamics. Elife 12, e83477. doi:10.7554/eLife.83477

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Pandey, A. K., Mulligan, M. K., Williams, E. G., Mozhui, K., Li, Z., et al. (2016). Joint mouse-human phenome-wide association to test gene function and disease risk. Nat. Commun. 7, 10464. doi:10.1038/ncomms10464

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, E. G., Wu, Y., Jha, P., Dubuis, S., Blattmann, P., Argmann, C. A., et al. (2016). Systems proteomics of liver mitochondria function. Science. 352 (6291), aad0189. doi:10.1126/science.aad0189

PubMed Abstract | CrossRef Full Text | Google Scholar

Wishart, D. S., Feunang, Y. D., Guo, A. C., Lo, E. J., Marcu, A., Grant, J. R., et al. (2018). DrugBank 5.0: a major update to the DrugBank database for 2018. Nucleic Acids Res. 46 (D1), D1074–D1082. doi:10.1093/nar/gkx1037

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, F., Zhao, Y., and Zhang, H. (2022). Ocular autonomic nervous system: an update from anatomy to physiological functions. Vis. (Basel). 6 (1), 6. doi:10.3390/vision6010006

CrossRef Full Text | Google Scholar

Xia, Y., Zhao, K., Liu, D., Zhou, X., and Zhang, G. (2023). Multi-domain and complex protein structure prediction using inter-domain interactions from deep learning. Commun. Biol. 6 (1), 1221. doi:10.1038/s42003-023-05610-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, S., Gumpper, R. H., Huang, X. P., Liu, Y., Krumm, B. E., Cao, C., et al. (2022). Molecular basis for selective activation of DREADD-based chemogenetics. Nature 612 (7939), 354–362. doi:10.1038/s41586-022-05489-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., and Skolnick, J. (2004). Tertiary structure predictions on a comprehensive benchmark of medium to large size proteins. Biophys. J. 87 (4), 2647–2655. doi:10.1529/biophysj.104.045385

PubMed Abstract | CrossRef Full Text | Google Scholar

Zielkiewicz, J. (2005). Structural properties of water: comparison of the SPC, SPCE, TIP4P, and TIP5P models of water. J. Chem. Phys. 123 (10), 104501. doi:10.1063/1.2018637

PubMed Abstract | CrossRef Full Text | Google Scholar