Synergistic regulation mechanism of iperoxo and LY2119620 for muscarinic acetylcholine M2 receptor

Muscarinic acetylcholine receptors are GPCRs that regulate the activity of a diverse array of central and peripheral functions in the human body, including the parasympathetic actions of acetylcholine. The M2 muscarinic receptor subtype plays a key role in modulating cardiac function and many important central processes. The orthosteric agonist and allosteric modulator can bind the pocket of M2. However, the detailed relationship between orthosteric agonist and allosteric modulator of M2 is still unclear. In this study, we intend to elucidate the residue-level regulation mechanism and pathway via a combined approach of dynamical correlation network and molecular dynamics simulation. Specifically computational residue-level fluctuation correlation data was analyzed to reveal detailed dynamics signatures in the regulation process. A hypothesis of “synergistic regulation” is proposed to reveal the cooperation affection between the orthosteric agonist and allosteric modulator, which is subsequently validated by perturbation and mutation analyses. Two possible synergistic regulation pathways of 2CU-I178-Y403-W400-F396-L114-Y440-Nb9 and IXO-V111-F396-L114-Y440-Nb9 were identified by the shortest path algorithm and were confirmed by the mutation of junction node. Furthermore, the efficiency of information transfer of bound M2 is significant higher than any single binding system. Our study shows that targeting the synergistic regulation pathways may better regulate the calcium channel of M2. The knowledge gained in this study may help develop drugs for diseases of the central nervous system and metabolic disorders.


Introduction
Muscarinic receptors are involved in a large number of physiological functions including regulation of heart rate and contractile forces, contraction of smooth muscles, and the release of neurotransmitters. 1 There are ve subtypes of muscarinic AChRs based on pharmacological activity: M1-M5 which regulate the activity of a diverse array of central and peripheral functions in the human body, including the parasympathetic actions of acetylcholine. 2 All ve are found in the CNS, while M1-M4 are also found in various tissues. 3 The activation of the M2 receptor in the heart is important for closing calcium channels in order to reduce the force and rate of heart contraction. 4 Hormonal, neurotransmitter, visual, and olfactory signaling is largely regulated by a versatile class of membrane receptors, referred to as G-protein coupled receptors (GPCRs). 5 The family of the muscarinic acetylcholine receptors has been widely studied as a model system for the interaction of allosteric modulators with GPCRs. 6 These GPCRs become important targets for agonist and antagonist drug design.
Muscarinic receptors have attracted particular interest due to their ability to bind small molecule allosteric modulators. 7 Many crystal structures were recently obtained for active states of the M2 and M3 muscarinic receptors, 8,9 including structures of a GPCR bound to an allosteric modulator. The crystal structure of M2 muscarinic receptor was released in 2013 (pdb code: 4MQS). 2 The structure ( Fig. 1) with POPC membrane includes an extracellular amino terminus, seven trans-membrane helices, a cytoplasmic carboxyl terminus, extracellular and intracellular loops, and two binding sites for orthosteric agonist iperoxo (IXO) and positive allosteric modulator of LY2119620 (2CU). The agonist IXO can enhance the effect of allosteric modulator 2CU. 10 However, the relationship between orthosteric agonist IXO and allosteric modulator 2CU is still unclear, as well as how the agonist and allostery information transfer the regulation information from binding site to G-protein interaction area. In this study we intend to answer the following questions. (1) Is there a synergistic regulation mechanism? (2) How does the information transfer from the regulation sites to the binding site of Nb9-8? And (3) what are the regulation pathways in M2?
To answer these questions, dynamics correlation networks were constructed based on all-atom molecular dynamics (MD) simulations of the wild type and mutant M2 receptor. From the comparisons of the networks between bound and free M2, a hypothesis of "synergistic regulation" and two possible synergistic regulation pathways were proposed to explain the agonist enhancing the allosteric effect, which are subsequently validated by both network weaken and mutation analyses.

Molecular dynamics simulation
The atomic coordinates of the active M2 in complex with Nb9-8 and IXO are extracted from the Protein Data Bank (pdb code: 4MQS), the coordinates of bound M2 from 4MQT, free M2 from 4MQS, and M2/2CU from 4MQT removing ligand. 2 D103E and Y403A are the mutants mentioned in the previous works to reveal the activation mechanism of IXO and 2CU. The mutant structures were constructed using SYBYL-X 2.1.1. 11 All structural visualizations were conducted in PyMOL 1.7.9. 12 All initial structures were rst minimized in SYBYL-X 2.1.1 to eliminate any possible overlaps or clashes. As the M2 receptor is a membrane protein, a protein/membrane complex was generated by CHARMM-GUI web-based Membrane Builder, M2 protein was loaded in POPC membrane before MD simulation. AMBER12 was used to perform MD simulations with periodic boundary conditions. 13 Hydrogen atoms were added using the LEaP module of AMBER12. Counterions were used to maintain system neutrality. All systems were solvated in a truncated octahedron box of TIP3P waters with a buffer of 10Å. The pairwise interactions (van der Waals and direct Coulomb) were computed with a cutoff distance of 8Å. Particle mesh Ewald (PME) was employed to treat long-range electrostatic interactions in AMBER12. 14 The improved ff99SBildn force eld was used for the intramolecular interactions. The Langevin thermostat was used in the preparation runs with a friction constant of 1 ps À1 and the Berendsen thermostat was used in the production runs. 12,13,15 All MD simulations were accelerated with the CUDA version of PMEMD in GPU cores of NVIDIA Tesla K80.
To relieve any further structural clash in the solvated systems, initial minimization with macromolecule frozen was performed using 500-step steepest descent minimization and 2000-step conjugate gradient minimization. Next the whole system was followed by 1000-step steepest descent minimization and 19 000-step conjugate gradient minimization. Aer minimization, a 400 ps heating up and a 200 ps equilibration in the NVT ensemble at 310 K were performed before MD simulation was conducted in the NPT ensemble at 310 K. To compare the difference among the chosen M2 systems, seven systems were simulated and the detailed simulation conditions are listed in Table 1.

Interaction analysis
Interaction analysis was handled with in-house soware. 16,17 Residues, substrates, or effectors are in hydrophobic interaction when mass centers of their side chains or ligands are closer than 6.5Å. A previous study has shown that residue chargecharge interactions up to 11Å were found to contribute to protein/protein binding free energies. 18,19 Thus, electrostatic (i.e., charge-charge) interactions are assigned when the distances between mass centers of their side chains are less than 11Å. Hydrogen bond is dened that the distance between the donor and acceptor is less than 3.5Å and the bond angle is larger than 120 . 17,[20][21][22] MM/PBSA free energy calculation Free-energy calculation together with MD simulations can provide qualitative predictions of M2-ligand binding energies. The free energy of binding (DG bind ) is estimated by eqn (1).
where G R+L , G R , and G L is the free energies of the complex, the isolated M2 and ligand, respectively. In the MM/PBSA approach, 20,23-25 each free energy term in eqn (1) is calculated as eqn (2). where E bond is the molecular mechanics bond energy, which is the sum of bond, angle and dihedral energies; E vdw is the molecular mechanics van der Waals energy contribution; E elec is the molecular mechanics electrostatic energy; G PB and G SA are polar and non-polar contributions to the solvation energy; T is the absolute temperature and S is the solute entropy. All the binding energies were calculated using the PBSA model in MMPBSA. 18,[23][24][25] Considering that the equilibration period may affect the result, only the last 50 ns of MD trajectories was used to calculate the binding free energy taking snapshots every 50 ps.

Correlation network analysis
Every amino acid or ligand was dened as one node for dynamics network. 26 The uctuation correlation between any pair of nodes was calculated with eqn (3).
wherer i ðtÞ ¼r i ðtÞ À hr i ðtÞ〉,r i ðtÞ is the position of node i at time t, and h$i represents time averaging, and r i (t) is the position of node i at time t. These elements were conveniently organized as a covariance matrix for a simulated system. In the current study, the covariance matrix for each system was constructed using snapshots (every 2 ps) in the last 50 ns of all simulated trajectories. Besides nodes, "edge" that transfers regulation information from one node to another is also a key concept in network construction. An edge is dened between any two nodes without covalent bond but with heavy atoms closer than 4.5Å over 75% sampling time. The strength of the edge is dened as the absolute value of the inter-node correlation (C ij ) between nodes i and j. The number of connected edges at each node is dened as the degree of the node. Correlation-weighted degree, which is the summation of strengths of all edges connected to a given node, indicates the importance of the node. Aer the network construction, network topological analyses were performed using Cytoscape3.1.1. 27 For community analyses, the Girvan-Newman algorithm was utilized. Most of the network analysis tools utilized here were developed by the Luthey-Schulten Group 27-29 and the strategy of process was successfully used in the allostery of riboswitch and protein. [30][31][32][33][34][35][36][37][38] Shortest-path length (SPL) analysis The Dijkstra algorithm was used to calculate the shortest path between any two nodes in the network. 39 The SPL is dened as the distance between all pairs of nodes in the network. The length of shortest path between two nodes is the shortest travelling distance from one node to another node along the network edges. This method has been used previously to predict residues important in allosteric signaling within protein families. 40

Stability of M2 in free and bound states
Ca RMSD relative to the initial structure shows that 160 ns simulations are sufficient for the equilibration of the wild types and mutant systems at 310 K (ESI Fig. S1  To quantitatively identify the conformational adjustment, the principal component analysis (PCA) was carried out on the four systems to study the motion. The rst three most dominant components (named PC1, PC2, and PC3) represent over 95% of the overall uctuations in each system (ESI Fig. S4 †), it showed the Ecl2, TM5-TM7 and Nb9-8 areas contribute most of the uctuations with shearing motion. To clearly display the motion, the results of each system are shown in Fig. 2. For the free M2, some conformational changes were observed. Upon binding allosteric effectors of IXO and 2CU, N terminal domain of bound M2 are in the opening motion mode. This motion was induced by the binding of allosteric effectors. We can also nd similar motion mode for N terminal domain of M2/2CU and M2/IXO.
To conrm the reliability and robustness of MD simulations for membrane protein M2, the binding free energies between  Table  2. Although there is limitations in MMPBSA approach, it is still possible to discover some interesting conclusions from the different systems. 41 The MMPBSA analysis shows that the binding energy between M2 and IXO is À41.87 kcal mol À1 , it is always much lower than that of mutant D103E with À37.73 kcal mol À1 . This is in qualitative agreement with experiment in that the D103E mutant signicantly decreases the binding affinity of IXO. 2 At the same time, the binding free energy between IXO and M2 for the Y403A mutant is about 3.82 kcal mol À1 higher than that of bound M2. This is also consistent with experiment in that Y403A mutant resulted in signicantly reduced affinity for the orthosteric agonists. 42 In summary, the MD simulation and binding free energy analysis are in good accord with previous experiment, indicating the good quality of the current model of M2.

Correlation networks of free and M2 complex are different
To reveal the allosteric mechanism, the dynamic correlation network analysis based on the covariance matrices calculated by in-house soware was performed to illustrate the residue uctuation correlation network and the network parameters for each system are listed in Table 3, the map of network is shown in ESI Fig. S5. † Table 3 indicates that the values of network parameter for bound M2 are the highest among these four systems for max average degree, max edge betweenness, the node betweenness, and the nodes number of degree > 10. This suggests that the network characters for bound M2 are signicant different from those of other systems Y403, W400, F396, L114, V129, and F274 with high degree were marked in the bound M2 network and most of them were found on the shortest-pathway to be discussed below. This indicates that the allosteric effectors indeed increase the correlation of nodes for bound M2.
To measure the distribution of weighted nodes with high degree and betweenness in each system, the distribution map is shown in ESI Fig. S6. † This gure indicates that the nodes of bound M2 with weighted degree larger than 10 are focused on Nb9-8 and transmembrane regions. The nodes of bound M2 with weighted betweenness are mostly located at the segments of TM3, TM6, TM7, and Nb9-8. These results are shown that the bound M2 has a signicant different information transfer network from other systems. The 2D-projection of correlation networks is shown in Fig. 3 for free M2, M2/IXO, M2/2CU, and bound M2 which reveal the whole network with nodes and edges.

Community networks are different
The analysis of correlation network shows that the networks are signicantly different among the four systems. To reveal the information transfer pathways, the Girvan-Newman was used to split these networks into communities. The community of four systems are shown in Fig. 4. This gure indicates that the community network of the free M2 has 14 communities including two isolated clusters, however the community network becomes more centralized upon binding to 2CU and/or IXO. There are only 12 communities with one isolated cluster for M2/2CU; 12 communities with one isolated cluster for M2/ IXO, and 10 communities without isolated cluster for the bound M2. This demonstrates that the effects of 2CU and IXO binding can optimize the integrity of the community in the complexes. 2 The information ow can freely transfer from 2CU and IXO regulation sites to the binding site of Nb9-8.
Furthermore, structural analysis indicates that there is one hydrogen bond, ve hydrophobic, and seven electrostatic interactions between 2CU and M2, and one hydrogen bond and six hydrophobic interactions between IXO and M2 with population higher than 40% (ESI Fig. S7 †) in the bound M2. Thus orthosteric agonist IXO and allosteric modulator 2CU introduce  strong interactions with M2. Furthermore, the binding free energy between 2CU and M2 for bound M2 is signicantly lower than that for M2/2CU (listed in Table 4). Therefore, the hypothesis of "synergistic regulation mechanism" was proposed to explain the regulation mode for M2. The previous work reports that 2CU shows strong positive cooperation with iperoxo and selectively enhances the affinity of the orthosteric agonist for the M2. 2

Modications to perturb the community networks
In order to validate the hypothesis of synergistic regulation, modications to weaken the interactions between IXO or 2CU and M2 for bound M2 were used to study the effects on the community networks. These modications were realized by just deleting the edges between 2CU or IXO and M2 nodes in the network. The community networks of weakened systems are shown in ESI Fig. S8. † Weakening operations will induce more fragmented community networks and decrease the efficiency of information transfer. For weakened IXO, the number of community increases from 10 to 15 including two isolated vertexes and the Nb9-8 community becomes fragmented. For weakened 2CU, the number of community number increases from 10 to 17 with two isolated nodes and the information could not transfer from 2CU to the effector site of Nb9-8. When weakening both IXO and 2CU, the community number increases from 10 to 17 with 4 isolated nodes. Furthermore, Nb9-8 domain are split into three communities and more than that of bound M2. In summary, these modications lead to signicant repartition of the community networks. This nding conrms that the interactions between 2CU/IXO and M2 indeed inuence the community network and further support our hypothesis of synergistic regulation mechanism.

Validation of regulation mechanism by mutation
In order to evaluate the receptor activation of M2 by IXO, the D103E mutation of M2 was analyzed in a previous experiment, 2 demonstrating that D103E abolished agonist-induced M2 receptor activation and reduced 380-fold affinities for iperoxo.
To further evaluate the hypothesis of "synergistic regulation mechanism", the network of D103E is constructed and shown in Fig. 5A, indicating that the network is different from that of wild type for bound M2. The community network was shown in Fig. 5B. This gure shows that the number of community increase from 12 to 15 with two isolated clusters in D103E compared with the wild type. This suggests that D103 plays a key role in information transfer of bound M2, this is consistent with the previous works that D103E resulted in signicantly reduced affinity for the orthosteric agonist. The interaction between IXO and M2 was shown in ESI Fig. S9 † which is signicant reduced. In summary, even if the information can transfer from the IXO orthosteric site to the effector site of Nb9-8, the efficiency is reduced due to mutant of D103E. These results further conrmed the hypothesis of synergistic regulation mechanism.

Synergistic regulation pathway for M2 complex
The network and community analyses of the wild type and mutants conrm that the IXO and 2CU binding induces  synergistic activation of M2. Next, it is natural to identify the synergistic regulation pathway based on the shortest path analysis which was used to identify the regulation pathway between the regulation sites and the binding site of Nb9-8. All he shortest-pathway length (SPL) was calculated from allosteric site to effector site of Nb9-8, the possible pathways with varies SPL were shown in ESI Fig. S10. † It indicates the information ow can transfer from the binding site to effector site via different pathways. The average shortest-pathway length (SPL) for each system is listed in Table 5. Two regulation pathways were identied in the bound M2: 2CU-I178-Y403-W400-F396-L114-Y440-Nb9 and IXO-V111-F396-L114-Y440-Nb9 (Fig. 6). This indicates that I178, Y403, W400, F396, L114, Y440, and V111 are key nodes for information transfer in bound M2. These strongly correlated residues may play critical roles in transferring the information of regulation. Note also that the SPL for the community network of bound M2 (7.80) is also shorter than the sum of SPLs of similar pathways in M2/IXO (4.40) and M2/2CU (6.60). This indicates that the efficiency of information transfer is 42.55% for M2/IXO and 27.07% for M2/2CU if the bound M2 efficiency is 100%. This further supports the hypothesis of synergistic regulation mechanism.
Mutation and network perturbation were used to validate the two regulation pathways. Y403 is one of the key residues on the shortest pathway from 2CU to Nb9. The dynamics correlation network of Y403A is shown in ESI Fig. S11. † This indicates that the network is signicant different from that of wild type. The SPL of mutant Y403A was also listed in Table 5, indicating that the SPL for Y403A was much longer than that for wild type. The efficiency of information transfer for Y403A is just 56.41% of wild type. This suggests that Y403 plays a key role in the information transfer and conrms the regulation pathway. This is consistent with the previous work that Y403A decreased the efficacy and functional affinity of agonist Ach. 42 At the same time, F396 is at the junction of the two proposed pathways and might be crucial for the information ow. In order to further conrm the key function of F396, mutant F396A was simulated for 160 ns. The binding free energy between IXO and M2 was about À28.40 kcal mol À1 and much higher than that of wild type. The community of F396A is shown in ESI Fig. S12 † with one isolated cluster. The number of community increases from 10 of bound M2 to 15. Therefore, the efficiency of information ow is signicantly reduced to 59.05% due to the fragmentation of the network. This suggests that the proposed regulation pathways, at least F396 before the junction, play key roles in the synergistic regulation. The previous works also conrm that the mutations of Y403, W400, and F396 in M2 decrease the binding affinity of agonist or efficacy.

Discussion
Previous work reports that 2CU can directly activate the M2 receptor and enhance the affinity of IXO. 2 Furthermore, 2CU has strong positive cooperative with IXO and enhances the potency of agonists. 43 D103E mutant experiment reduced the decrease of the binding affinities for IXO ($380-fold). 2 In simulation, the binding free energy increased about 1.67 kcal mol À1 for D103E mutant, consistent with the previous work. 2 The previous work reports that ECL2 mouth becomes smaller upon the activation of M2 by IXO or 2CU to prevent other ligand from moving into the closed, active extracellular vestibule. 8 TM5 and TM6 are closer to each other in M2-active state than those in the inactive state. 2 The conformer of ECL2 mouth from MD simulation is shown in Fig. 7, which indicates that the dimension of ECL2 mouth decreases from 7.0Å to 5.8Å upon the binding of IXO and 2CU. The close propensity could  M2/2CU  2CU-A191-W155-A194-V111-L115-F119-C124-K221-E223-Nb9  9 6.60 AE 1.00 M2/IXO IXO-S151-W400-N436-Y440-R121-Nb9 5 4.40 AE 1.10 F396A 2CU-W422-G425-L428-I431-T434-T37-G40-L43-F451-Nb9 9 5.63 AE 1.10 IXO-W400-Y196-V199-T203-W207-A212-R216-K218-E220-Nb9 9 7.58 AE 1.20 Y403A 2CU-Y104-A109-W148-I144-S64-Y60-N58-N444-Nb9 8 5.89 AE 1.10 IXO-N108-N113-S64-Y60-N58-N444-Nb9 prevent other ligand to enter into extracellular vestibule, this is consistent with the previous work. The alignment between inactive and active states for bound M2 is shown in ESI  Fig. S13, † it indicates that TM5 and TM6 are closer in the active state than in the inactive state. The previous work reports that 2CU forms two hydrogen bonds with TRP422, and ASN410 and induces conformation change in binding area of ext-TM4-TM5-TM6. 2 The hydrogen bond interaction between 2CU and M2 is shown in ESI Fig. S14. † This gure shows that 2CU has two hydrogen bonds with these residues to form a hydrogen bond network. These results are in agreement with the previous work.

Conclusion
Residue/residue uctuation correlation network and shortestpathway analysis were used to reveal the regulation mechanism of M2 upon binding to allosteric effectors of 2CU and agonist IXO. The results suggest that the dynamics correlation network of bound M2 has more hub nodes than that of free M2, M2/2CU, or M2/IXO. The community network of the bound M2 is clustered into a complete community without isolated cluster. However, there are two isolated clusters for free M2. The information ow can freely transfer from the IOX and 2CU regulator sites to the binding site of Nb9-8 for bound M2. The efficiency of information transfer for bound M2 is higher than that of M2/2CU or M2/IXO. The binding free energy between the antagonist IXO and M2 for the bound M2 is signicantly lower than that of M2/2CU or M2/IXO system. Therefore, we propose the hypothesis of synergistic regulation in bound M2 system and the hypothesis was validated through mutation and network perturbation. At the same time, two regulation pathways were found and conrmed by Y403A and F396A mutants.

Conflicts of interest
The authors declare that there is no conict of interest.