Homology modeling of hM3R
Using crystal structure of rat M3R (PDB: 4DAJ), a homology model of hM3R was constructed, whose sequence identity was 62 %. The predicted structure of hM3R contains three different domains: extracellular (EC), transmembrane (TM), and cytoplasmic (CT) domains (Fig. 1A). Like most GPCRs, G protein interaction with hM3R is predicted to take place in the CT domain upon agonism31-33. Upon superimposition of rat M3R crystal structure to hM3R homology model, we observed conservation of TM and EC domains (Fig. 1B). Our homology model was further refined using Protein Preparation and Loop Refinement features in Maestro (Schrödinger Inc., USA). Binding site analysis of our homology model was done using Sitemap feature in Maestro (Fig. 1C), and the top-ranking binding site can be found in the EC domain (Fig. 1D). This is consistent with the crystal structure of rat M3R, where tiotropium is in complex within the EC domain. Upon superimposition, tiotropium from the crystal structure fits tightly inside the top-ranking binding site of hM3R identified by Sitemap (Fig. 2A). Upon molecular docking of tiotropium to hM3R homology model using Glide (Schrödinger Inc., USA), the predicted binding mode is nearly identical to the binding mode of tiotropium obtained from the crystal structure of rat M3R (Fig. 2B). Upon comparison of the binding forces of tiotropium in the rat crystal structure to our hM3R homology model, it was determined that they were nearly identical, including the presence of π-π stacking exchange between tiotropium’s aromatic ring and tryptophan residues (Fig. 2C and 2D). Together, the results indicate that molecular docking with our hM3R homology model can predict ligand binding with high fidelity.
High-throughput virtual screening (HTVS) of M3R muscarinic agonists and antagonists
The 3D chemical structures of M3R agonists and antagonists were obtained from Drug Banks database, and their structures were further refined using LigPrep in Maestro. The binding modes of top three scoring M3R agonists revealed that cation-π exchange can be found between a positively charged nitrogen ring and residues with aromatic sidechains, like TYR149, TRP504, and TYR507 (Fig. 3A). Another notable force is hydrogen bonding between an oxygen group and ASN508. The binding modes of top three scoring M3R antagonists also contained cation-π exchange between a positively charged nitrogen ring and residues with TYR149, TRP504, and TYR507 (Fig. 3B). ASN 508 can also be seen interacting with oxygen molecules. However, compared to M3R agonists, the antagonists all shared presence of π -π stacking exchange between their aromatic rings and TRP504. In fact, this interaction is absent in all of the hM3R agonists available in the database likely because they lack aromatic side groups in their chemical structures. Table 1 displays the results of HTVS that lists the docking scores and checks presence of π- π stacking exchange, which is found in 12 out of 15 M3R antagonists that we docked. It’s also interesting to note that M3R antagonists achieved higher docking scores than the agonists likely because of the added stability via π-π stacking exchange.
MD simulation for interaction analysis
To analyze additional interactions found in hM3R agonism and antagonism, we conducted MD simulations using docked complexes with the highest scoring agonist and antagonist, pilocarpine and tiotropium, respectively. The protein-ligand complexes with the highest population during 100 ns simulation were extracted for binding free energy calculation and structural comparisons using Trajectory Clustering tool in Maestro (Schrödinger Inc., USA). The video analysis of MD simulations for hM3R-pilocarpine and hM3R-tiotropium complexes can be found in Supplemental Fig. 1 and 2, respectively. The root mean square distance (RMSD) analysis of hM3R in complex with pilocarpine was performed, where lower the RMSD indicates tighter binding (Fig. 4A). Pilocarpine achieved an RMSD of roughly 2.5 Å towards the end of the simulation, which indicates stable binding34. The protein RMSD equilibrated to roughly 7.0 Å, which implies the protein assumed a stable conformation towards the end of the simulation. On the other hand, tiotropium achieved an RMSD of roughly 3.0 Å towards the end of the simulation (Fig. 4B). This value also indicates stable binding but not as tight as pilocarpine-hM3R complex likely because of the increased molecular mass of tiotropium. The protein RMSD equilibrated to roughly 10 Å, which indicates the protein assumed a stable conformation towards the end of the simulation
Interaction fraction analyses of pilocarpine-hM3R and tiotropium-hM3R complexes were also performed. Interaction fraction value of 1.0 means the ligand made contact with a binding residue during 100% time of the simulation. Values greater than 1.0 indicate more than one binding forces. For example, in the interaction fraction analysis of pilocarpine-hM3R (Fig. 5A), TYR149 has an interaction fraction value of 1.498 and makes contact with pilocarpine via hydrogen bonds, hydrophobic forces (cation-π exchange), and water bridges. The top three contributors to the binding interaction for pilocarpine were TYR149, ASP148, and SER152 with scores of 1.498, 1.217, and 0.914, respectively. In the interaction fraction analysis of tiotropium-hM3R (Fig. 5B), the top three contributors to the binding interaction were ASN508, TRP504, and TYR507 with scores of 1.755, 1.437, and 1.083, respectively. It’s interesting to note that TRP504 has almost 3-times the interaction fraction value with tiotropium (1.477) compared to pilocarpine (0.591). Furthermore, the top three binding residues with pilocarpine are upstream at positions 148-152, whereas they are downstream at positions 504-508 for tiotropium. Ligand-protein contact analysis of pilocarpine-hM3R revealed ASP148, TYR149, and SER152 intimately associating with the nitrogen ring of pilocarpine (Fig. 5C). Meanwhile, the ligand-protein contact analysis of tiotropium-hM3R displayed preservation of π-π stacking exchange via TRP504 during most of the simulation (Fig. 5D), which is consistent with the results of HTVS.
To analyze amino acid side chain stability, root mean square fluctuations (RMSF) analysis was performed, where higher the RMSF values indicates higher instability of amino acid side chain. In RMSF analysis of pilocarpine-hM3R complex (Fig. 6A), none of the amino acid side chains exceed RMSF value of 9.0 Å. On the other hand, in RMSF analysis of tiotropium-hM3R complex (Fig. 6B), there were three residues with RMSF values greater than 10.5 Å: LEU456 (11.03 Å), SER457 (12.16 Å), and LYS459 (11.89 Å). In fact, these three residues were much more stable in pilocarpine-hM3R complex: LEU456 (7.84 Å), SER457 (5.66 Å), and LYS459 (6.95 Å). It’s interesting to note that these residues span positions 456-459, which are near the top three binding residues (TRP504, TYR507, and ASN508) from the interaction fraction analysis of tiotropium-hM3R complex. Table 2 displays the difference in stability among shared binding residues of pilocarpine-hM3R and tiotropium-hM3R complexes. Overall, the binding residues of tiotropium have lower RMSF and hence higher stability. In fact, TRP504 had the highest difference in RMSF of 0.27 Å. This suggests TRP504 is in a more stable and rigid conformation when interacting with tiotropium likely due to added stability brought by π-π stacking exchange.
Through structural comparisons, we were able to locate TRP504 in helix 6 of hM3R (Fig. 7A and 7B). When we superimposed helix 6 of piloarpine-hM3R and tiotropium-hM3R (Fig. 7C), we identified a 2.1 Å shift in the α-carbon of TRP504 when interacting with tiotropium for π-π stacking. This led to an angle change of 18 ° upstream of TRP504, including the cytoplasmic domain, where the unstable residues LEU456, SER457, and LYS459 are directly located at. When we superimposed the cytoplasmic domain directly upstream of helix 6 (Fig. 7D), we could find LEU456, SER457, and LYS459 further away from the protein backbone in hM3R-tiotropium than in hM3R-pilocarpine, consistent with the high RMSF values we observed. It’s important to note that these three residues were more stable in hM3R-pilocarpine complex. Perhaps the increased fluctuations of these residues in the cytoplasmic domain are preventing G protein interactions and hence antagonism of hM3R.
Binding free energy analysis by MM/GBSA
To compare the stability of protein-ligand complex between pilocarpine and tiotropium, we used the MM/GBSA method to obtain the binding free energy (ΔGbind). Table 3 lists the ΔGbind of pilocarpine-hM3R and tiotropium-hM3R complexes, whose overall predicted binding free energies were −46.81 kcal/mol and -96.82 kcal/mol, respectively. This suggests that tiotropium forms a more stable complex with hM3R than pilocarpine does. Furthermore, near two-fold change in binding free energy indicates that once tiotropium binds to hM3R, its binding mode is stable, and the TM domain that contains the binding pocket likely remains in a stable, rigid conformation, which is consistent with the video analysis of tiotropium-hM3R complex in Supplemental Fig. 2.