Introduction

Skeletal muscle nicotinic acetylcholine receptors (AChRs) mediate synaptic transmission between motor neurons and muscle fibers. These pentameric ligand-gated ion channels (pLGICs) are composed of four different subunits (2 α1, β1, δ, and ε) arranged symmetrically around a central ion-conducting pore (Rahman et al., 2020; Unwin, 1993). Binding of the neurotransmitter acetylcholine (ACh) or other agonists to two target sites in the extracellular domain (ECD) increases the probability of a global conformational change that opens a distant (∼60Å) channel ‘gate’ in the transmembrane domain (TMD), to allow water and cations to cross the membrane.

At the start of the channel-opening isomerization the agonist binding sites switch from low to high affinity (LA→HA), and at the end of the isomerization the gate switches from closed to open (C_LA→O_HA) (Purohit et al., 2013). Without a bound agonist, the probability of the O shape (P_O) is small and baseline activity is negligible (Jackson, 1986; Nayak et al., 2012; Purohit & Auerbach, 2009). Agonists (A) increase P_O because they provide more favorable binding energy to the target in its active vs resting conformation, making O_HA more stable than C_LA. Agonists also bind strongly to the transition state that separates O from C and so also increase the channel-opening rate constant (Grosman et al., 2000).

The binding free energy difference between high and low affinity (ΔG_HA-ΔG_LA) drives receptor activation and determines agonist efficacy, the high-concentration limit of the dose-response curve. In contrast, the binding free energy ratio (ΔG_LA/ΔG_LA) generates a correlation between the curve’s maximum and midpoint and determine agonist efficiency, the fraction binding energy converted into movements that trigger receptor conformational change.

Adult-type muscle AChRs have an unfavorable unliganded isomerization free energy (+8.3 kcal/mol at −100 mV) (Nayak et al., 2012) and 2 approximately equal and independent binding sites a (Nayak & Auerbach, 2017). For ACh, ΔG_HA=-10.2 kcal/mol and ΔG_HA =-5.1 kcal/mol per site, so adding their differences at 2 sites to that of the intrinsic isomerization, gives −1.9 kcal/mol. A free energy difference is proportional to the logarithm of an equilibrium constant. Hence, after binding 2 neurotransmitters the AChR gating equilibrium constant increases from 7.4×10⁻⁷ to 25, and P_O increases from ∼10⁻⁶ to 0.96 (Fig. 1A). Binding energy from 2 neurotransmitters turns the receptor approximately from off to on.

Fig. 1:

Our goal is to associate ΔG_LA→ΔG_HA changes with underlying structural rearrangements at the binding site. Here, we use MD simulations to explore the conformational dynamics of the LA→HA transition at the α-δ neurotransmitter binding site of the Torpedo AChR (6UVW; (Rahman et al., 2020)). Binding free energies calculated in silico aligned with those estimated experimentally by using electrophysiology, so it was possible to interpret the agonist energy changes with restructuring of the binding pocket. It is difficult to draw sweeping conclusions because we examined only 4 agonists and only one wild-type site. However, in the LA→HA transition for all constructs the agonist somersaults, loop C droops to reposition αY190 and rearrange an aromatic pocket that contracts and becomes more hydrophobic. Also, in all systems the initial LA and final HA conformations of the binding site were connected by a brief intermediate state. For some but not all agonists, loop F moves in and a conserved salt bridge shortens. The results suggest that that agonist mobility is important determinant of the increase in binding free energy, and that rather than acting as a lid that closes, loop C is a positioning device that inserts an affinity-increasing side chain into the pocket core.

Results

Hold

Traditionally, receptor activation is divided into binding (formation of a ligand-protein complex) and gating (receptor isomerization). Here, we are concerned only with movements at binding sites that connect these two steps. In binding, the agonist (A) diffuses to interact with the target (A+C=^AC), after which a local rearrangement establishes the LA complex (^AC=AC_LA) (Jadey & Auerbach, 2012). In gating, the site again rearranges to produce the HA complex (AC_LA=AC_HA), after which the rest of the protein domains (ECD, TMD and gate) restructure to establish AO_HA. We call the two local rearrangements of the binding site ‘catch’ (free energy change ΔG_LA) and ‘hold’ (ΔG_HA-ΔG_LA) (Fig. 1B). The end of binding is AC_LA and the start go gating is AC_HA, so hold is the local rearrangement that connects these two activation steps.

The free energy changes associated with the other events in binding and gating, namely diffusion to the site (the chemical potential) and restructuring of receptor domains, are the same for the agonists we have examined and so are not relevant to this study. Furthermore, energy coupling between the binding pocket and distant residues (>15 Å) typically is weak (Gupta et al., 2017) (Fig. S1A). Only catch and hold determine differences between agonists with regard to affinity (ΔG_LA, ΔG_HA), efficacy (ΔG_HA-ΔG_LA) and efficiency (1-ΔG_LA/ΔG_HA) (Nayak et al., 2019). In AChRs, ΔG_LA and ΔG_HA have been measured experimentally for 23 agonists and 53 binding site mutations (Indurthi & Auerbach, 2023) (Fig. S1B). The structures generated in the MD simulations were identified as AC_LA or AC_HA by aligning ΔG_LA and ΔG_HA values calculated in silico with those from experiments.

The cryo-EM structures of AChRs include apo, resting-C with a bound toxin and high-affinity desensitized (Rahman et al., 2022; Rahman et al., 2020; Zarkadas et al., 2022). Binding sites of desensitized AChRs have a similar (perhaps identical) high affinity as do O sites (Auerbach, 2020). The end states of the hold rearrangement are short-lived and have not been captured as structures.

Because ΔG_LA and ΔG_HA are independent of energy changes outside the binding pocket, in the simulations we removed the TMD and studied only α-Δ subunit ECD dimers. The starting structure for the hold transition was (PDB: 6UWZ) after toxin removal that is essentially the same as the apo structure (Zarkadas et al., 2022)) (Fig. 1C). The 4 agonists we examined were acetylcholine (ACh), carbamylcholine (CCh), epibatidine (Ebt) and epiboxidine (Ebx).

Docking

Binding energy changes in hold, AC_LA=AC_HA (Fig 1B). To understand the structural underpinnings, we began by removing α-bungarotoxin from the high-resolution (2.69 Å) structure 6UWZ (Rahman et al., 2020) and docking agonists into the empty pocket at the α-δ subunit interface (Fig. 1C). The docking simulations had strong binding scores, ranging from −5.82 to −9.18 kcal/mol, suggesting stable interactions. However, unlike the agonist orientation apparent in cryo-EM structures of desensitized AChRs (Gharpure et al., 2019; Noviello et al., 2021; Rahman et al., 2022; Walsh et al., 2018; Zarkadas et al., 2022) (Fig. 1D), the agonist’s tail pointed away from (cis) rather than towards (trans) the complimentary (non-α) subunit. This was the case for the top 3 docking scores for all 4 ligands (Fig. 2A). Because the cis orientation was not observed in any cryo-EM (HA) structure, we examined the top 200 docking scores for CCh. All poses were cis, with the tail pointing towards the α subunit. The cis orientation was also apparent in docking of agonist into a putatively LA binding site of a homology model (Tripathy et al., 2019). To start, we associated the agonist orientation in the hold end states as cis in AC_LA versus trans in AC_HA.

Figure 2.

LA→HA transition

We performed 200 ns MD simulations for each of the docked complexes, in triplicate (with different seed values) (Fig. 2B top). The MD simulations for all three replicates demonstrated similar root-mean-square-deviation (RMSD) patterns, with differences consistently below 1Å (Table S1) suggesting that the trajectories were robust and reproducible. Although there were initial deviations, all systems became stable by ∼120 ns. The root-mean-square-fluctuation (RMSF) of the C_α atoms, employed to identify regions of flexibility, showed instability mainly in loops C and F around the pocket (Fig. S2). The fluctuations observed at the base of the ECD were anticipated because the TMD that offers stability here was absent.

In addition to the fluctuations in these loops, during the 200 ns simulations the ligand orientation changed dramatically, cis-to-trans. At the start of the simulations the agonist’s tail points towards α, and at the end it points towards δ (Fig. 2B bottom). This rearrangement represents a substantial rotation (a ‘flip’) about the agonist’s main cationic center (N⁺) that remains close to αW149 throughout the somersault. The final MD configuration (with CCh) was well-aligned with the CCh-bound cryo-EM desensitized structure (7QL6.PDB) (RMSD < 0.5 Å) (Fig. S3), further demonstrating that the simulation had converged.

To test whether rearrangements at the binding pocket in the simulations themselves enforce the trans configuration, we removed the bound CCh from the final state MD structure and re-docked all four agonists. The preferred poses for all 4 agonists were all trans (Fig. 2C). In addition, binding energies from docking scores were higher for the converged MD structure compared to the initial binding site structure. As a further control we carried out MD simulations of a pentamer docked with ACh and found similar structural changes at the binding pocket compared to the dimer (Fig. S4). These results indicate that a dimer is appropriate to study local rearrangements of an AChR neurotransmitter binding site and suggest that in the MD simulations the binding site restructures from AC_LA to AC_HA conformation (Fig. 2D).

Principal component analysis

To identify the most prominent conformational states of the receptor during the MD simulations, principal component analysis (PCA) was carried out based on the most significant fluctuations (principal components). Considering the first two principal components (PC-1 and PC-2) that capture the most pronounced C_α displacements, the protein backbone fluctuations revealed three energy minima (m1, m2, and m3) that mainly reflect the stabilities of loops C and F (Fig. 3). For all four agonists, the trajectories quickly reach m1 (≤20 ns), then transition through m2, and culminate in m3 (180-200 ns). In the simulations progress in time was sequential, m1→m2→m3.

Figure 3.

Table 1.

Accordingly, we hypothesized that m1 corresponds to AC_LA and m3 corresponds to AC_HA. For all agonists, the m2 population has greater variance (more instability) compared to m1 and m3. This, along with its position in the sequence, suggests that m2 represents an unstable intermediate configuration within hold, between AC_LA and AC_HA. For all agonists and trajectories, m3 had the least variance (was most stable), again supporting convergence by 200 ns.

Transitions between the 3 energy minima varied with the agonist. Always, m1 was apparent first (at 5-14 ns for ACh; 15-24 ns for CCh; 21-41 ns for Ebt; 20-56 ns for Ebx), followed by m2 (20-50 ns; 50-60 ns; 60-100 ns; 100-120 ns). All systems converged to m3 (180-200 ns) (Fig. 3). The bottom of the m1, m2 and m3 energy wells represent the most stable configurations of each population and serve as a reference point for the rest of energy well.

Although the bottom of the well for 3 energy minima from PCA represent the most stable overall conformation of the protein, they do not convey direct information regarding agonist stability or orientation. To incorporate this information, we performed a cluster analysis of the ligand using frames from the bottom of each PCA well as inputs. The top three clusters, each having RMSD of ≤ 1Å, shared a similar ligand orientation (Fig. S5) and were selected to compute binding free energy using MM-PBSA. The remaining frames were disregarded. The fraction of frames accepted for free energy calculations ranged from a minimum of 20% (ACh, m1) to a maximum of 71% (Ebt, m3). On average, the fraction of frames selected from each energy minimum was ∼50%, showing the dynamic nature of the ligand in the pocket (Table 1).

Computed vs experimental binding energies

The above procedure defined populations of structures that we hypothesize represent the end states of hold, AC_LA (m1) and AC_HA (m3), plus an intermediate (m2). We used MM-PBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) computations of the selected structures to estimate binding energy of each agonist to each population. MM-PBSA is particularly useful to compare binding affinities between different states and among agonists (Kollman et al., 2000; Rastelli et al., 2010; Sun et al., 2014). In our use, MM-PBSA computations included both enthalpic (ΔH) contributions from molecular interactions (van der Waals, electrostatics, solvation) and entropic (TΔS) contributions calculated using normal mode analysis. By capturing the dynamic interplay between the ligand and receptor as well as the influence of the surrounding solvent environment this approach enables computation of a Gibbs free energy (ΔG). However, there are limitations with regard to MM-PBSA accurately predicting absolute binding free energies (Genheden & Ryde, 2015; Hou et al., 2011) that depends on parameterization of the ligand (Oostenbrink et al., 2004).

The calculated and experimental ΔG values are shown in Table 1. For ACh and CCh, there is good agreement between ΔG_m1 and ΔG_LA and between ΔG_m3 and ΔG_HA. For these agonists, in silico values overestimated experimental ones only by ∼8% and ∼25%. The agreement was not as good for the other 2 agonists, as calculated values overestimated experimental ones by ∼45% (Ebt) and ∼130% (Ebt). However, the fractional overestimation was approximately the same for ΔG_LA and ΔG_HA.

Despite these discrepancies in absolute values, the trend of an increase in binding free energy, m1 versus m3, is consistent with the assignment, AC_LA versus AC_HA. For all agonists, the stability of the ligand-receptor complex increases significantly in the m1→m3 transition (Fig. 4A). Further support for the assignment is from comparing calculated vs experimental efficiencies (η) that depend on ΔG_LA/ΔG_HA (Nayak et al., 2019). Table 2 shows that η values calculated in silico agree closely (≤5%) with those measured experimentally, for all agonists. Apparently, MD simulations of a single pdb file successfully model the hold transition. The match between calculated and experimental efficiencies is evidence that that m1 represents low-affinity AC_LA and m3 represents high-affinity AC_HA.

Figure 4.

In the AC_LA→AC_HA conformational change, all agonists show an increase in van der Waals (vdW) interaction energy. The increase was most pronounced with Ebt (−9.77 kcal/mol) and least for Ebx (−3.54 kcal/mol), with CCh (−8.81 kcal/mol) and ACh (−6.71 kcal/mol) falling in between (Fig. 4, Table S2). Also, in AC_LA→AC_HA, restructuring of the binding pocket further contribute to generating a smaller, more-tightly packed pocket. This is shown in Fig. 4C as a reduction in binding pocket volume and in Fig. 4D as a decrease in the number of water molecules. The agonist’s flip into the pocket results in de-wetting and a more hydrophobic environment.

Flip-Flop

The ligand somersault (flip) was prominent in all simulations of hold. All agonists undergo a large rotation about the N⁺-αW149 fulcrum (Fig. 5; Video 1). The extent of rotation was 131° for CCh, 156° for ACh, 148° for Ebt and 176° for Ebx. Another consistent and prominent restructuring event in AC_LA→AC_HA was the downward displacement (flop) of loop C toward the core of the pocket. This movement is conserved in related receptors and thought to play a role in receptor activation by agonists (Basak et al., 2020; Hansen et al., 2005; Hibbs et al., 2009) (see below). In the simulations, we quantified loop C flop by measuring the displacement of its tip (C_α of αC192) relative to its starting position (at 0 ns). The distance traveled at 200 ns varied with the agonist and was 9.2 Å for CCh, 5.5 Å for Ebt, 5.2 Å for ACh and 3.3 Å for Ebx (Fig. 5). We do not discern a pattern with regard to agonist affinity, efficacy or efficiency.

Figure 5.

In hold, residues that form the pocket cavity move, apparently to accommodate the re-orientation of the ligand. In AC_LA αY190, αY93, and δW57 are spaced apart, providing an open slot that accommodates the agonist’s tail in the cis orientation. However, in AC_HA αY190 is positioned closer to αY93, effectively creating a unified surface that disallows the cis position of the tail (Fig. S6). However, in the simulations the event sequence is not clear so perhaps flip precedes flop, and agonist re-orientation first creates the space for αY190 to occupy.

αW149 (in loop B) makes a cation-π bond with N⁺ that is important for ligand stabilization (Xiu et al., 2009; Zhong et al., 1998). In the simulations, the position and orientation of αW149 relative to N⁺ remained almost unchanged for all 4 agonists, in both AC_LA (cis) and AC_HA (trans) configurations (Fig. 6). Apparently, the cation-π bond acts as a slippery anchor that maintains stability of the N⁺ position even as the tail rotates about it.

Fig. 6:

Loop C flop results in re-positioning of αY190 into the aromatic pocket. αY190 is the main ligand binding energy source, nearly twice as strong as αW149 (Purohit et al., 2014), perhaps because deprotonation of its hydroxyl group generates a cation-anion bond with N⁺ (Bruhova & Auerbach, 2017). An interaction between αY190 and the αK145-αD200 salt bridge, suggested to be the trigger for the gating isomerization (Bruhova & Auerbach, 2017; Mukhtasimova et al., 2005). In AC_LA→AC_HA, loop C flop restructures the pocket to put N⁺ between the αY190 and αW149 rings, and αY190 closer to the salt bridge (Fig. 6).

Other hold rearrangements

In AC_LA (m1), the cis orientation is stabilized by the hydroxyl of αY93 that is proximal to a functional group on the agonist tail. Also, the tail nitrogen in CCh, Ebt, and Ebx (amide or amine) forms a hydrogen bond with the hydroxyl group of αY198, but only in the intermediate m2 state. This interaction appears to guide the flip and provide transient stability to the intermediate state. The absence of a secondary nitrogen in ACh possibly explains its alternate re-orientation route, as only this ligand was observed to flip counterclockwise (Fig. 5). It may also explain why the m2 basin in the PCA plots for ACh is especially broad and shallow (Fig. 3).

In AC_HA (m3), the αY190/αY198/δW57 rings and the αW149 indole interact with N⁺, and the hydroxyl group of αY93 forms a hydrogen bond with the ligand head group (Celie et al., 2004; Hansen et al., 2005). However, Ala substitutions of these 5 aromatic residues have quite different consequences and result ACh binding free energy losses of +4.4, +1.9, +0.1, +2.3 and +0.9 kcal/mol, respectively (Purohit et al., 2014). In AC_HA, the tail functional moiety (carbonyl or secondary ring amine) in the trans orientation is sandwiched between αC192, αC193 and αY198 (loop C), αT150 (loop B) and δC108, δN109, δL111, and δL121 (loop E) (Fig. S7).

αK145 is complex and undergoes agonist-dependent hold rearrangements. In AC_HA this residue i) makes a salt bridge with αD200, ii) has a cation-π interaction with αY93, iii) contacts αY190 (possibly through the hydroxyl group) and iv) approaches loop F of the complementary subunit (Zarkadas et al., 2022). Further, the mutation αK145A lowers the efficiencies of ACh and CCh (from 0.51 to 0.41) but has no effect on those of Ebt and Ebx (0.41) (Indurthi & Auerbach, 2023).

This agonist dependence with regard to function has a parallel in structure. In the hold simulations, the αK145-αD200 distance (Å) increases from 4.4 in apo to 6.1 or 5.1 in AC_LA (CCh or ACh). The spread may relate to the agonist tail being inserted into the αY93 slot in the cis orientation. In AC_LA →AC_HA, for some agonist the αK145-αD200 distance (Å) shortens substantially to 2.6 or 2.8, making the bridge more rigid (ACh or CCh). However, with Ebt and Ebx the αK145-αD200 distance remains unstable, fluctuating between 3-6 Å throughout the course of the simulation. Likewise, αK145-αY190 distance stabilizes at 5.3 Å or 6.1 Å in AC_HA with ACh and CCh but continues to fluctuate with Ebt and Ebx. In contrast, the αY190-αD200 distance remains stable in AC_HA for all ligands. This suggests that the instability in the salt αK145-αD200 and αK145-αY190 distances can primarily be attributed to αK145 (Fig. S8).

In AC_HA, the functional group found in the tail of each agonist -- the carbonyl group in CCh and ACh and the secondary ring amine in Ebt and Ebx - forms a hydrogen bond with the carbonyl of δN109 and the amine of δL121 in the complementary δ subunit via a structural water (Fig. 6), as reported elsewhere (Blum et al., 2010; Gharpure et al., 2019; Zarkadas et al., 2022). This water-mediated bonding appears to be crucial for stabilizing the agonist in the trans orientation. The start of the agonist’s flip, m1→m2, begins the establishment of this water-mediated bond that is completed in m2→m3, with all 4 agonists (Fig. 6).

Discussion

In diffusion-catch-hold activation of AChRs (; Fig. 1B), the binding free energy increase in the last step generates the force that initiates the receptor isomerization. We identified in MD simulations (Fig. 3) the hold end states, AC_LA (m1) and AC_HA, (m3), as well as an intermediate configuration (m2). Our approach to understanding the conformational dynamics of hold is to compare calculated with experimental binding free energies, ΔG_LA and ΔG_HA. There are many experimental measurements for many different AChR agonists and several binding sites (Indurthi & Auerbach, 2023; Jadey & Auerbach, 2012; Purohit et al., 2014), but here we describe the comparison for only 4 agonists and only the wt α-δ site.

The alignment between experimental and calculated free energies was closer for ACh and CCh than for Ebt and Ebx (Table 1). The match improved when free energy ratios rather than absolute values were compared. Individual in silico values overestimated experimental ones by approximately the same factor (different for each ligand) that, however, cancelled in the ratio. We do not know the origin of this factor but speculate that it could be caused by errors in ligand parameterization. Regardless, the results show that agonist efficiency (η) can be estimated from only agonist and receptor structure (one PDB file), and that η is a useful metric for associating structure with function.

In the AC_LA→AC_HA simulations, two conspicuous binding site rearrangements were a cis→trans somersault of the agonist (a flip) and an up→down droop of loop C (a flop). The loop C movement is well known, but its functional role has been puzzling. The extent of the flop is correlated with agonist potency in ACh binding proteins (Basak et al., 2020; Du et al., 2015; Polovinkin et al., 2018), but voltage-clamp fluorometry of other receptors shows it to be similar with agonists and antagonists (Chang & Weiss, 2002; Munro et al., 2019; Pless & Lynch, 2009). In AChRs, truncation of loop C prevents activation by agonists but has no effect on unliganded activation, but agonist sensitivity is maintained if αY190 and αY198 are intact (Purohit & Auerbach, 2013). These results suggest that the downward displacement of loop C in AChRs is mainly associated with the generation of strong binding energy.

Perhaps the most puzzling result regarding loop C is that the agonist association rate constant k_on is slower to AC_LA (up position) than to AC_HA (down position) (Grosman & Auerbach, 2001; Nayak & Auerbach, 2017). Moreover, whereas k_on to AC_LA is correlated with agonist potency (the dissociation rate constant k_off is approximately the same for all agonists), k_on and k_off to AC_HA both are agonist dependent. The flop of loop C in hold removes a barrier to make k_on approximately diffusion-limited, and also reduces discrimination between ligands (molecular recognition.

Loop C droop resembles the ‘clamshell closure’ apparent at other receptors binding site (Chen et al., 2017; Hansen et al., 2005; Hibbs & Gouaux, 2011). However, in the light of the above experimental results, this label is misleading because ‘closing’ (or ‘capping’) a binding pocket is expected to slow rather than speed k_on. For this reason, we use less functionally loaded words to describe the downward displacement of loop C, namely a ‘flop’ or ‘droop’. In AC_LA→AC_HA loop C indeed moves closer to αC of αW149, to cover the pocket and create higher affinity. However, it does not do so by imposing a steric barrier. We did not study the catch rearrangement and do not speculate about the nature of the barrier that limits the agonist association rate constant to C_LA that has loop C in an up position.

Importantly, loop C flop repositions αY190, the main determinant of both LA and HA binding energies (Bruhova & Auerbach, 2017; Purohit et al., 2014). In AC_LA→AC_HA, this side chain fills the slot adjacent to αY93 that otherwise allows the cis orientation in AC_LA (Fig. S5). Further, flop moves αY190 closer to a salt bridge between αK145-αD200, and compacts the space between αW149 and N⁺ (Fig. 6). Overall, the simulated structures suggest that loop C flop increases binding energy both by repositioning αY190 near N⁺ and by allowing the cis-trans flip of the ligand. Rather than closing a door, the downward movement of loop C rearranges the furniture.

We interpret cis-trans rotation of the agonist about the N⁺-αW149 fulcrum (Fig. 5) with caution. Although the trans orientation is apparent in desensitized AChR cryo-EM structures, the cis orientation exists only in silico. Notwithstanding this limitation, in AC_LA→AC_HA the agonist’s flip into the binding pocket and the repositioning αY190 by loop C flop were observed consistently in all simulations with all agonists, and rationalize the increase binding energy in the hold transition. We speculate that the agonist flip triggers the rest of the hold restructuring but cannot rule out flop promoting flip or synchronous rearrangements. We have been unable to establish clearly the flip-flop temporal sequence.

For the agonists we examined, the somersault positions the tail deep into the binding cavity where it interacts via a structural water with the complementary subunit. It appears that agonists rotation increases binding free energy by displacing water and increasing vdW interactions and protein stability. A hydrophobic binding pocket is associated with functional activation in other ion channels (Furukawa et al., 2005; Sigel & Steinmann, 2012) and is observed in many proteins including β2-adrenergic receptor, opioid receptors, and HIV-1 protease (Huang et al., 2015; Rasmussen et al., 2011; Wlodawer & Erickson, 1993). However, these effects on the AChR pocket cannot be the only basis for the higher affinity AC_HA because several potent agonists lack a tail, for example tetramethylammonium. We hypothesize that both the reorganization of the aromatic groups in the pocket around N⁺ and binding site compaction/de-wetting contribute to stronger binding, in proportions that remain to be determined. Regardless, the agonist’s somersault shows that mobility is one feature that distinguishes a ligand from a covalently linked side chain (Indurthi & Auerbach, 2023). Not only can an unconnected agonist come and go, it can also reposition itself dramatically to induce a conformational change of binding pocket to affect a change in binding energy.

Movements analogous to flip and flop have been identified in other proteins. Multiple ligand positions in a binding pocket have been documented in 5-HT3 receptors (Basak et al., 2020), δ-opioid receptors (Shang et al., 2016), HIV-1 protease (Klei et al., 2007) and transthyretin (Klabunde et al., 2000). In enzymes, substrate motion has been observed in hexokinase (Bennett & Steitz, 1978) and alcohol dehydrogenase (Plapp, 2010). Regarding flop, substrate-induced loop displacements are known to reposition residues to enhance catalysis, for example in chymotrypsin (Ma et al., 2005) and triosephosphate isomerase (Brown & Kollman, 1987; Nickbarg et al., 1988). The agonist flip and loop C flop in AChRs are additional examples of ligand mobility and loop motions undergirding function.

The PCA plots (Fig. 3) show 3 structural populations, one of which (m2) occurs in time between m1 (AC_LA) and m3 (AC_HA). In m2, the agonist tail interacts with αY198, possibly to guide the clockwise rotation of the ligand and reduce the probability of a return to the cis orientation. ACh was the only agonist to rotate counterclockwise, perhaps because it lacks the appropriate functional group in its tail. The MM-PBSA calculations indicate that ΔG_m2 (that does not have an experimental correlate) is intermediate between ΔG_LA and ΔG_HA for ACh and CCh, but least stable for Ebt and Ebx (Table 1). We do not have an explanation for this difference. However, we observe m1→m2, m2→m1, m2→m3 and m3→m1 transitions in the trajectories (Video 1), so the depth of the m2 energy well (relative to the barriers that separate it from m1 and m2) contributes to agonist efficacy and efficiency.

MD simulations did not reveal an interaction between the αK145-αD200 salt bridge distance and the degree of loop C displacement (Cheng et al., 2006). However, this distance discriminate between η classes by decreasing AC_LA→AC_HA in high (ACh and CCh) but not in low (Ebt and Ebx) η classes, consistent with experimental results (Indurthi & Auerbach, 2023). All 4 agonists are high efficacy, so this result suggests that perturbation of this salt bridge is not a necessary requirement for the change from LA to HA. After loop C flop, the repositioned αY190 connects N⁺ to the salt bridge, suggesting that this interaction (rather than bridge tightening) could help trigger the gating transition (Beene et al., 2002; Gay & Yakel, 2007; Mukhtasimova et al., 2005; Padgett et al., 2007; Pless et al., 2011). The simulation did not reveal the conformational changes (forces) that drive the next step of the isomerization (ECD restructuring). However, the observation that for all agonists, interactions with the complementary subunit and pocket dehydration terminated (at 200 ns) the hold rearrangement lead us to speculate these may contribute.

Loop E (connecting β5 and β6 strands) has been implicated in ligand discrimination (Basak et al., 2020; Muroi et al., 2009; Pless & Lynch, 2009). The double mutation F106L+S108C in the β2 subunit of the α4β2 nicotinic AChR subtype significantly decreases the action of Ebt without altering that of ACh (Tarvin et al., 2017). Loop E residues are responsible for the binding energy difference of α-γ versus α-Δ interfaces because hydrogen bonds influence the position and participation of γW55 (Nayak et al., 2016), and for reducing instability of the α-Δ site consequent to loop C mutations (Purohit & Auerbach, 2013). MD simulations with more agonists and more binding sites are needed to elucidate the role of loop E in the AC_LA→AC_HA transition.

Agonists belong to efficiency (η) classes in which ΔG_LA/ΔG_HA is the same for all members. So far 5 η-classes have been identified experimentally in adult-type AChRs, indicating that are 5 pairs of AC_LA/AC_HA structures (Indurthi & Auerbach, 2023). Although in the simulations loop F movements and salt bridge distances were agonist dependent, with only 4 ligands we cannot associate these rearrangements with η classes.

The key requirement for allostery is that ligand binding energy be weaker to the resting versus active conformation of the target. Natural selection must generate both low and high affinity binding sites, so our attention must be on restructuring in both ^AC=AC_LA (catch) and AC_LA=AC_HA (hold). What prevents loop C flop in apo? What barrier(s) slows k_on to apo? Where is the agonist in the encounter complex? What forces promote the cis orientation? What is the role of water in agonist orientation, both cis and trans? What are the mechanisms and consequences of a wider salt bridge, apo versus AC_LA? What structural classes are associated with efficiency classes? The correlation in the free energy changes in catch and hold that entangles binding with gating implies a corresponding correlation in structural change. We anticipate that additional MD simulations of catch, hold, and post-hold domain rearrangements, for more agonists and with binding site mutations, will lead to a more complete understanding of how agonists turn on the AChR.

Materials and Methods

Hardware and Software

All computational studies were carried out on a Linux Ubuntu 20.04 based workstation with 126 GB RAM, 64 CPU cores, and two RTX A5000 24 GB GPU cards, as well as an Ubuntu 20.04 based server equipped with 132 GB RAM, 96 CPU cores, and two RTX A5000 24 GB GPU cards. The software utilized includes the Maestro release 2022-2 graphical user interface (GUI) of the Schrödinger software suite (BioLuminate, Schrödinger, LLC, NY, USA), Amber22 (Case et al., 2005), VMD 1.9.3 (Humphrey et al., 1996), ChemDraw Professional 15.1 (Purushotham et al., 2022), ChimeraX (Pettersen et al., 2021), and PyMOL (PyMOL Molecular Graphics System, Version 2.0, Schrödinger, LLC).

Protein preparation

Protein preparation was conducted using the Protein Preparation Wizard (PPW) (Purushotham et al., 2022) and OPLS4 force field. Bond orders were assigned to atoms based on the CCD database, with zero bond orders given to bonded metals. Missing hydrogen atoms were inserted, and disulfide bonds were placed appropriately. During preprocessing, structural refinement was performed with Maestro’s Prime utility, optimizing and adding any missing side chains and residues. Missing residues were formulated by combining SEQRS information from the PDB files and optimized using Prime’s ePLOP module. The ionization state of all heterogeneous groups, including ligands, bound metals, and charged amino acids, were determined using the EpiK tool at a pH range of 7.0 ± 2.0. Subsequent force field minimization allowed for water sampling and basic restraint minimization, during which only hydrogen atoms could be freely reduced.

The hydrogen bonding (H-bond) network underwent significant modification. The terminal chi angle for Asn, Gln, and His was sampled through 180-degree flips, altering their spatial H-bonding capabilities without affecting the fit to electron density. The neutral and protonated states of His, Asp, and Glu, as well as two His tautomers, were sampled. Hydrogens on hydroxyls and thiols were also analyzed to enhance the H-bond network. Finally, the ProtAssign method was employed in two modes: a "regular" mode, which completed in seconds, and an "exhaustive" mode that considered many more states and could take minutes or hours, depending on the H-bond network complexity.

Ligand Preparation

The two-dimensional SDF structures of carbamylcholine (CCh), acetylcholine (ACh), epibatidine (Ebt) and epiboxidine (Ebx) were retrieved from the PUBCHEM database. LigPrep (LigPrep, Maestro 11.2.014, Schrödinger LLC) was utilized to construct ligands, produce stereoisomers, and generate tautomers using the OPLS46 force field. Additionally, the Epik module was employed to desalt and protonate the ligands at pH 7.0. Default values were applied to the remaining parameters, and LigPrep was also used to prepare the ligands ACh, CCh, Ebt, and Ebx, with the option "retain specified chiralities" engaged.

In terms of ligand state assessment, Epik can calculate a penalty to quantify the energy cost of forming each state in solution based on the Hammett and Taft methods. This EpiK state penalty is entirely compatible with the GlideScore used for docking, as it is computed in kcal/mol units. This compatibility facilitates the examination of how the EpiK state penalty affects the GlideScore. The DockingScore in Glide represents the total of the GlideScore and the Epik state penalty and served as the basis for final ranking and enrichment evaluation. Furthermore, Epik offers a technique for handling metal binding states, involving a change in the pH range within which the state is generated. This approach underscores the comprehensive considerations made in preparing the ligands for docking and subsequent analyses.

Grid Generation and Molecular Docking

The receptor grid complex of 6UWZ.pdb was generated after removing alpha-bungarotoxin from the binding pocket using Glide (Glide, Maestro 11.2.014, Schrödinger LLC). A virtual grid box with dimensions 20Å × 20Å × 20 Å was formed, centered within the C and F loops in the α-Δ binding site of the energy-minimized resting state. This grid box was used to generate a receptor grid, and the default values were assigned for other parameters. Docking was performed using the ligand docking wizard built into the Glide module of Schrödinger, with the default scaling factor and partial charge cutoff set to 0.80 and 0.15, respectively. The docking procedure consisted of two stages, extra precision (XP) followed by post-processing with prime MM-GBSA. The best dock poses were selected based on dock score and MM-GBSA energy for further investigation.

Molecular Dynamics Simulations

MD simulations were conducted using the Amber22 software suite within the Amber module as detailed elsewhere (Singh et al., 2021; Srivastava et al., 2022). Ligand topology was prepared via the Leap module in AmberTools22, employing force field ff19SB (Tian et al., 2020) and the Generalized Amber Force Field (GAFF) (Vassetti et al., 2019). The AM1-BCC strategy and GAFF2 were used to assign atomic charges and additional parameters. Before simulation, the systems were minimized and equilibrated following a five-step minimization process, with each step involving 10,000 energy minimization steps. Systems were heated under NVT ensembles in two consecutive runs from 0–300 K for 1 ns each, followed by a 1 ns simulation at 300 K and 1 bar pressure under the NPT ensemble for equilibration. An additional 5 ns equilibration was performed before the 200 ns production MD for each system, conducted under periodic boundary conditions (PBCs) and utilizing an NPT ensemble with 300 K and 1 bar pressure. Long-range electrostatic interactions were treated using the particle mesh Ewald method (PME), including a direct space cutoff of 12.0 Å and van der Waals interactions. In production MD, a time step of 2.0 fs was employed, and the SHAKE algorithm was used to keep bond lengths at equilibrium. Isobaric (NPT) conditions were maintained using the Berendsen barostat, with temperature monitored by a Langevin thermostat. Coordinates from the production MD were recorded in the trajectory file every 20 ps, resulting in a minimum of 10,000 frames. Trajectories were analyzed using amber-tools CPPTRAJ (Roe & Cheatham, 2013), and further analysis and figure generation were completed with VMD, Pymol, Xmgrace and ChimeraX. To ensure reproducibility, MD trajectories were generated in triplicate.

Cluster Analysis

Cluster analysis of the ligand was performed using TCL scripts within the Visual Molecular Dynamics (VMD) software, as described elsewhere (Mittal, Tonk, et al., 2021). In brief, this script automates the identification of frames from the most stable conformations based on PCA minima. Using these frames as input, the algorithm employed RMSD values of ligand positions to cluster similar orientations. Specifically, clusters with an RMSD of ≤ 1Å were considered to share a similar ligand orientation. Frames that did not fit into these top clusters were automatically excluded by the script.

Principal Component Analysis (PCA) and Free Energy Landscape (FEL)

PCA identifies collective motions in atomic-level simulations of macromolecules (Mittal, Srivastava, et al., 2021; Srivastava et al., 2022). It emphasizes significant patterns by reducing noise in MD simulations, revealing essential dynamics underlying conformational changes. All equilibrated MD simulation trajectories were analyzed by PCA computations using the CPPTRAJ algorithm in Amber tools.

Prior to PCA, least-square fitting of the protein’s Cα atoms was carried out with respect to a reference structure (average protein coordinates) to eliminate rotational and translational degrees of freedom in the simulation box. A positional covariance matrix C (of dimensions 3N × 3N) was constructed using the Cα atom coordinates, and this matrix was diagonalized to produce eigenvectors (directions of motion) and eigenvalues (magnitudes). The elements of the covariance matrix C were obtained using:

where X_i and X_j are the Cartesian coordinates of the Cα atoms, and N is the number of Cα atoms. The ⟨⟩ notation represents the ensemble average of atomic positions in Cartesian space. The resulting principal components (eigenvectors) were sorted by the total motion they captured, generating 3N eigenvectors. The free energy landscapes (FELs; FIg.3) were created using the "g_sham" module in GROMACS, based on the probability distribution of the top two principal components (PCs). The FEL plot aids in visualizing possible conformations and corresponding energy levels. The free energy values (G) were calculated as

where k is the Boltzmann constant, T is the absolute temperature of the simulated system (300 K), N_i is the population of bin i, and N_max is the population of the most-populated bin.

Binding Free Energy Calculation

The binding free energy of the protein-ligand docked complex (ΔG_bind) was calculated in two stages, after docking (using MM-GBSA in the Prime module of Schrödinger) (Mittal, Tonk, et al., 2021; Purushotham et al., 2022) and after MD simulation on equilibrated MD trajectories (using MM-PBSA/-GBSA in AMBER22).

After docking, for the static docked pose of the protein-ligand complex the MM-GBSA binding energy calculation was performed using the OPLS4 force field and the variable-dielectric generalized Born (VSGB) model in the Prime MM-GBSA module in Schrödinger 2022. We used ΔG_{bind_dock} to rank protein-ligand docked complexes, selecting those with the highest negative values and dock scores as the best ligands for further studies. After MD Simulation, in the MD-equilibrated system binding energies were calculated through the MM-PBSA/-GBSA protocol, implemented in AMBER22. This was performed over the most stable clustered poses extracted from the equilibrated trajectory at minima m1, m2, and m3 on the free energy landscape (FEL) over a 200 ns time frame. AMBER utilizes the conventional g_mmpbsa module for MM-PBSA calculations using MM-PBSA.py script. These methods calculate the energies of electrostatic interactions, van der Waals interactions, polar solvation, and non-polar solvation in the equilibrated trajectory using a single trajectory protocol (STP). The binding energy of each protein-ligand complex (ΔG_bind) was calculated as

where ΔG_complex, ΔG_protein, and ΔG_ligand represent the absolute free energies of the complex, protein, and ligand systems, respectively. Additional terms, such as ΔH (enthalpy), T and ΔS (temperature and total solute entropy), ΔG_gas (total gas phase energy), and ΔG_solv (solvation energy), contribute to the calculation.

The Poisson-Boltzmann calculations were performed using an internal PBSA solver in Sander, while Generalized Born ESURF was calculated using ‘LCPO’ surface areas. The vibrational mode entropy contributions (TΔS) for protein-ligand interactions were calculated and averaged using normal-mode analysis in Amber. The entropy term provides insights into the disorder within the system, as well as how this disorder changes during the binding process. This value was added to the enthalpy term computed by the MM-PBSA method to determine the absolute binding free energy (ΔG_bind). These comprehensive calculations allowed for a precise evaluation of the interactions between the protein and ligand, providing valuable insights into the binding process and its energetic characteristics.

Pocket Volume Calculation

To investigate the variations in binding site volume across all four agonist-bound cases, pocket volume analysis was carried out using POVME3.0 (Wagner et al., 2017). From each of the four systems, frames were extracted at consistent intervals with a stride of 2.0, culminating in representative sets of 5000 protein structures for every system. Before commencing the volume calculation, the trajectory was aligned, and frames were selected from VMD, serving as the initial input for this method. Subsequently, an inclusion region was defined, encapsulating all binding pocket conformations within the trajectory. This region’s definition involved building a sphere upon which the inclusion grid-points were computed, using different atoms of the ligand for this calculation.