Introduction

Adult vertebrate skeletal muscle nicotinic acetylcholine receptors (AChRs) are allosteric proteins that mediate synaptic transmission between motor neurons and muscle fibers. These pentameric ligand-gated ion channels 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 orthosteric sites in the extracellular domain (ECD) increases the probability of a global change in protein 3D conformation that activates a distant (~60.0 Å) allosteric site (a gate) in the transmembrane domain (TMD), allowing cations to cross the membrane.

At the start of the channel-opening isomerization the neurotransmitter sites switch from having a low to a high affinity for agonists (LA→HA), and at the end the gate switches from closed to open (C→O) (Grosman et al., 2000; Purohit et al., 2013). Hence, the reaction scheme for the ‘gating’ allosteric transition is C_LA⇄O_HA. Without any bound agonists the probability of having the O_HA shape (P_O) is ~10⁻⁶ and baseline current is negligible (Jackson, 1986; Nayak et al., 2012; Purohit & Auerbach, 2009). With or without agonists, the receptor turns on and off spontaneously (influenced only by temperature). Agonists increase P_O (that is, preferentially stabilize O_HA) because they provide net favorable binding energy to their target when it is in the active, high-affinity conformation. Agonists also stabilize the transition state that separates C_LA from O_HA and therefore increase P_O almost exclusively by increasing the channel-opening rate constant.

The difference in binding free energies, HA minus LA (ΔG_HA-ΔG_LA), drives the receptor’s allosteric transition and determines relative agonist efficacy (maximum P_O). Another ligand attribute we use is efficiency that is a function of the ratio of binding free energies, (1-ΔG_LA/ΔG_HA) (Nayak et al., 2019). Efficacy (efficiency) is the amount (fraction) of binding free energy delivered to the gating apparatus. In a dose-response curve, efficacy is manifest as the high-concentration asymptote and efficiency as the correlation between this maximum and EC₅₀.

Adult-type muscle AChRs have 2 approximately equal and independent neurotransmitter sites (Nayak & Auerbach, 2017). The net LA→HA binding free energy change (ΔG) for two neurotransmitters (−10.2 kcal/mol) adds to that for the unfavorable unliganded isomerization (+8.3 kcal/mol at −100 mV) (Nayak et al., 2012) making the free energy change of the allosteric transition with 2 bound ACh favorable (−1.9 kcal/mol). ΔG is proportional to the logarithm of the corresponding equilibrium constant (by −RT, where R is the gas constant and T the absolute temperature, equal to −0.59 at 23 ^°C). Accordingly, after 2 neurotransmitters have bound, LA→HA rearrangements at the 2 orthosteric sites increases P from ~10⁻⁶ to ~0.96 (Figure 1A).

Figure 1:

The relationship between ΔG_LA and ΔG_HA binding free energies and the agonist properties affinity, efficacy and efficiency are considered elsewhere (Auerbach, 2024). Here, our goal was to associate the ΔG_LA→ΔG_HA change in binding free energy (4 agonists) with restructuring events at a neurotransmitter site. Molecular dynamics (MD) simulations were used to explore conformational changes in the LA→HA transition at the α−δ neurotransmitter site of the Torpedo AChR (6UVW; (Rahman et al., 2020)). In brief, binding free energies calculated in silico matched those estimated by using electrophysiology in vitro (Indurthi & Auerbach, 2023), allowing the identification of LA versus HA structures in the simulations. In the LA→HA transition for all agonists, three prominent rearrangements were rotation (flip) of the ligand, downward displacement (flop) of loop C, and compaction and stabilization of the binding cavity (fix). The initial LA and final HA conformations of the site were connected by a brief intermediate state. The simulations suggest that the movement of a weakly-bound agonist starts the change to strong binding and initiates the protein’s allosteric transition.

Results

Hold

Traditionally, receptor activation is divided into distinct stages called ‘binding’ (formation of a ligand-protein complex) and ‘gating’ (receptor isomerization) (Figure 1). Here, we are concerned with the connection between these stages as molecular movements at a neurotransmitter site.

In AChRs, agonist binding is a 2-step process. First, the agonist (A) diffuses to the target and forms an encounter complex (A+C⇄A-C), after which a local ‘catch’ rearrangement establishes the LA complex (A-C⇄AC_LA) (Jadey & Auerbach, 2012). Gating, too, involves multiple steps (Gupta et al., 2017; Purohit et al., 2013). First, the site undergoes a local ‘hold’ rearrangement that establishes the HA complex (AC_LA⇄AC_HA), after which other protein domains restructure without changing agonist affinity, eventually opening the gate (Figure 1B).

The only agonist-dependent free energy changes in activation occur in the two local site rearrangements, ΔG_LA (in catch) and ΔG_HA−ΔG_LA (in hold). The terminal state in the ‘binding’ stage of activation is AC_LA and the initial state in the ‘gating’ stage is AC_HA. Our interest here is hold, AC_LA⇄AC_HA, the local rearrangement that connects these two states. The free energy changes associated with the other events in the activation sequence, namely the formation of the encounter complex and downstream domain restructurings, are approximately the same for all agonists and so are not relevant to this study.

Activation free energy changes incurred outside catch and hold, namely from the chemical potential associated with encounter complex formation (Phillips, 2020) and from rearrangements in other domains (and water), together represent an energy offset that is approximately agonist independent. With a few exceptions, energy coupling between the agonist binding pocket and distant residues (>15.0 Å) typically is weak (Gupta et al., 2017), so only the agonist’s energy changes in catch and hold determine the pharmacology attributes of affinity (ΔG_LA or ΔG_HA), efficacy (ΔG_HA−ΔG_LA) and efficiency (1-ΔG_LA/ΔG_HA). In AChRs, ΔG_LA and ΔG_HA have been measured experimentally by using electrophysiology for 23 agonists and 53 neurotransmitter site mutations (Indurthi & Auerbach, 2023).

Apo, resting-C with a bound toxin, and high-affinity desensitized AChR structures have been solved by using cryo-EM (Rahman et al., 2022; Rahman et al., 2020; Zarkadas et al., 2022). Desensitized AChRs have a similar (perhaps identical) high affinity as AO_HA (Auerbach, 2020; Nayak & Auerbach, 2017). Below, we describe the conformation dynamics that connect the end states of the hold rearrangement, AC_LA→AC_HA (Figure 1B, bottom), neither of which has been captured as a structure. However, as described below, they can be identified in MD simulations that use resting-C as a start.

Docking

We began by removing α-bungarotoxin from the high-resolution (2.69 Å) structure 6UWZ (Rahman et al., 2020) (Figure 1C) and docking agonists into the empty pocket. The docking simulations had strong binding scores ranging from −5.82 to −9.18 kcal/mol (Figure 2A), suggesting that the interactions were stable. However, unlike the agonist orientation apparent in cryo-EM structures of desensitized (HA) AChRs (Noviello et al., 2021; Rahman et al., 2022; Rahman et al., 2020; Walsh et al., 2018; Zarkadas et al., 2022) (Figure 1D), the docked agonist’s tail pointed towards (cis) rather than away (trans) from the α subunit (Figure 2B, top). This was the case for the top 3 docking positions for all 4 ligands.

Figure 2.

Because the cis orientation was not observed in any cryo-EM structure, we examined the top 200 docking scores for CCh. All poses were cis, with the tail pointing towards the α subunit and away from the main binding cavity. Apparently, in the apo structure the trans orientation is energetically unfavorable. Hence, to start our analysis of the hold rearrangement we hypothesized that the agonist orientations in the end states are cis in AC_LA versus trans in AC_HA.

The LA→HA transition

We performed 200 ns MD simulations for each of the docked complexes, in triplicate (with different seed values). Video 1 shows a simulation with CCh on 2 times scales, expanded for the initial 200 ns and condensed for the full 1.0 μs extension. The 200 ns simulations for all three replicates for all agonists showed similar root-mean-square-deviation (RMSD) patterns with differences consistently below 1Å (Figure 2C, Figure 2-Source Data 1) suggesting that the trajectories were sufficiently 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 showed flexibility around the pocket mainly in loops C and F (but not E) (Figure 2-figure supplement 1). The fluctuations at the base of the ECD were anticipated because the TMD that offers stability here was absent. As described below, the match between experimental binding energies (whole receptors) and those calculated in the simulations (ECD dimers) is further evidence that the ECD-TMD interfacial region does not influence agonist binding energy.

The ECD dimer in the absence of a bound agonist (apo) exhibited a similar RMSF and RMSD pattern for a 200 ns simulation, to that observed with ligand present (CCh) (Figure 2-figure supplement 2A, Figure 2-figure supplement 2B). The displacement of loop C apparent with agonists (see below) was also present but less pronounced in apo (Figure 2-figure supplement 2C).

In addition to the fluctuations in loops C and F during the simulations, the ligand orientation inverted, cis→trans. For all ligands, at the start of the simulations the agonist’s tail points towards α, but by the end points towards δ (Figure 2D). The agonist re-orientation represents a nearly complete somersault (a ‘flip’) about the agonist’s main cationic center (N⁺) that remains close to αW149 throughout (see Figure 5). The final MD configuration (with CCh) aligns well with the CCh-bound cryo-EM desensitized structure (7QL6.PDB) (RMSD < 0.5 Å) (Figure 2-figure supplement 3). Video 1 shows that the final trans orientation with CCh remains constant for at least 1 μs, indicating that structure had settled by 200 ns, with critical rearrangements conserved.

To test whether rearrangements at the binding pocket in the simulations enforced the trans configuration, we removed the bound CCh from the final 200 ns MD structure and re-docked the agonists. The preferred poses for all 4 ligands were trans (Figure 2A, bottom). In addition, binding energies from docking scores were higher for the 200 ns structures compared to the initial binding site structures, consistent with a LA-to-HA transition (Figure 2B). 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 (Figure 2-figure supplement 4).

Together, these results suggest that ECD dimers are adequate to study local rearrangements of the α–δ AChR neurotransmitter binding site. Further, they provide support for the hypotheses that in the simulations the site restructures from AC_LA to AC_HA, and that the preferred agonist pose is cis in LA versus trans in HA structures.

Principal component analysis (PCA)

To identify the most prominent conformational states of the receptor during the simulations, PCA was carried out based on the most significant fluctuations (principal components). Considering the first two principal components (PC-1 and PC-2) that captured the most pronounced C_α displacements (Figure 3-Source Data 1), backbone fluctuations revealed three energy minima (m1, m2, and m3) that reflect the stabilities of loops C and F (Figure 3). To assess consistency across independent simulations, we calculated the inner products of the top 10 principal components from each run. The resulting heatmaps reveal the degree of similarity between runs, with the inner product values of PC-1 and PC-2 reflecting considerable overlap (50-90%) in the principal dynamic modes captured in each simulation (Figure 3-figure supplement 1).

Figure 3.

For all 4 agonists, trajectories started in m1 and ended in m3 (Figure 3-figure supplement 2). Accordingly, we hypothesized that m1 corresponds to AC_LA and that m3 corresponds to AC_HA. For all agonists, the m2 population had a greater variance (was less stable) compared to m1 and m3. This, along with its temporal position, suggests that m2 represents an unstable intermediate configuration within the hold transition, between AC_LA and AC_HA. For all agonists and trajectories, m3 had the least variance (was the most stable).

The pronounced fluctuations in loop C are evident in residue-wise RMSF plots and align well with the contributions of the first two principal components (PC-1 and PC-2) (Figure 2-figure supplement 2D). The apo and CCh occupied proteins adopt similar conformational spaces during the 200 ns simulation (Figure 2-figure supplement 2E). However, presence of the agonist appears to increase conformational diversity, generating a broader distribution in the PCA plot. The bottom of the m1, m2 and m3 energy wells in Figure 3 represent the most stable configurations of each population and serve as reference points.

Computed versus experimental binding energies

Although the well bottoms in Figure 3 represent the most stable overall conformations, they do not directly convey information regarding agonist stability or orientation. To incorporate this information, we performed a cluster analysis of the ligand configuration using frames from the bottom of each PCA well as input. The top three clusters, each having RMSD of ≤ 1.0 Å, shared a similar ligand orientation (Figure 3-figure supplement 2) and were selected to compute binding free energies (Figure 4, Figure 4-Source Data 1, Figure 4-Source Data 2). The remaining frames were disregarded. The fraction of frames from each minimum accepted for free energy calculations ranged from 20% (ACh, m1) to 71% (Ebt, m3). On average, the fraction of frames selected was ~50%, showing the dynamic nature of the ligand in the pocket (Video 1).

Figure 4.

The above procedure defined populations of structures that we hypothesized represent the end states of hold, AC_LA (m1) and AC_HA (m3), plus a previously unidentified intermediate state (m2). We used PBSA (Poisson-Boltzmann Surface Area) computations of the selected structures for each agonist in each cluster population and compared binding free energies calculated in silico with those measured previously in vitro (Indurthi & Auerbach, 2023). There are limitations with regard to PBSA accurately predicting absolute binding free energies (Genheden & Ryde, 2015; Hou et al., 2011) that also depend on parameterization of the ligand (Oostenbrink et al., 2004). However, because we knew the actual binding free energies from electrophysiology experiments we used the PBSA calculations only to test our hypothesis that m1 is AC_LA and m3 is AC_HA.

Calculated PBSA and experimental ΔG values are compared in Figure 4A left (Figure 4-Source Data 1). For ACh and CCh there was good agreement between ΔG_m1 and ΔG_LA and between ΔG_m3 and ΔG_HA. For these agonists, in silico values overestimated experimental ones by ~8% and ~25%. The agreement was not good for the other 2 agonists, with calculated values overestimating experimental ones by ~45% (Ebt) and ~130% (Ebt). However, in all cases ΔG_m3 was more favorable than ΔG_m1, and the % overestimation was approximately the same for m1 versus m3. Despite the mismatches with the PBSA calculations, the comparison is consistent with the hypothesis.

Further support for the assignments comes from comparing calculated vs experimental efficiencies that depend on the binding free energy ratio, ΔG_LA/ΔG_HA, rather than absolute ΔGs. As described elsewhere, efficiency is the agonist’s free energy change in hold relative to its total free energy change, (ΔG_HA-ΔG_LA)/ΔG_HA or (1 −ΔG_LA/ΔG_HA) (Nayak et al., 2019). Figure 4A right (Figure 4-Source Data 1) shows that efficiency values calculated in silico agree almost perfectly with those measured experimentally. This result not only further supports the hypothesis that m1 represents AC_LA and m3 represents AC_HA, but also indicates that efficiency can be computed from structure alone. We conclude that MD simulations of a single pdb file successfully model the end states of the hold transition and that the trajectories likely reproduce the conformational dynamics in AC_LA→AC_HA.

All agonists show an increase in van der Waals (VdW) interaction energy in this rearrangement (Figure 4B). The calculated increase (PBSA calculations) was most pronounced with Ebt and least for Ebx, with CCh and ACh falling in between (Figure 4-Source Data 2). AC_LA→AC_HA restructuring of the binding pocket generates a smaller, more-tightly packed pocket (Tripathy et al., 2019), in part because of the increased VdW interactions (Video 1). This is shown in Figure 4C as a reduction in binding pocket volume and in Figure 4D as a decrease in the number of water molecules. The VdW contacts that form concomitant to agonist flip and loop C flop establish a more compact, hydrophobic and stable local environment (Video 1).

Flip-Flop-Fix

In AC_LA→AC_HA, a cis→trans rotation of the agonist (flip) was prominent in all simulations (Figure 5). All agonists underwent a >130° rotation about the N⁺-αW149 fulcrum. Also, in all instances the flip occurred near the start of the hold rearrangement.

Figure 5.

A second consistent and prominent restructuring event in AC_LA→AC_HA was the downward displacement (flop) of loop C toward the pocket center (C_α of αW149). This well-known ‘clamshell closure’ motion is conserved in related binding proteins and thought to play a role in AChR activation by agonists (Basak et al., 2020; Hansen et al., 2005; Hibbs et al., 2009). In the simulations, loop C flop was measured as the displacement of its tip (C_α of αC192) in m3 relative to its apo position. The distance traveled varied with the agonist (Figure 6-table supplement 1) but we do not discern a pattern in the extent of loop C flop with regard to agonist affinity, efficacy or efficiency.

In hold, residues that form the pocket cavity move, presumably to accommodate the reoriented ligand. In AC_LA, αY190, αY93, and δW57 are spaced apart, providing an open slot that accommodates the agonist’s tail in the AC_LA, cis orientation (Figure 5-figure supplement 1). However, in AC_HA loop C flop repositions αY190 closer to αY93, effectively filling the gap to create a unified surface.

The αW149, αY190, αY198 and αY93 rings surround the N⁺ (Figure 1D). In adult-type mouse AChRs, Ala substitutions of these 4 aromatic residues have different consequences and result in ACh LA binding free energy losses of 3.0, 2.7, 2.2 and 1.1 kcal/mol, respectively (Purohit et al., 2014). δW57A does not change ACh binding energy (Nayak & Auerbach, 2013). α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 distance between αW149 and N⁺ decreases only slightly, AC →AC (Figure 6, Figure 6-table supplement 1). This cation-π bond appears to act as a slippery anchor (a hinge) that maintains stability of the N⁺ position even as the tail rotates about it.

Figure 6:

In AChBP (perhaps C_LA), the hydroxyl group of αY93 forms a hydrogen bond with the ligand head group (Celie et al., 2004; Hansen et al., 2005). In the AChR AC_LA structure, the cis orientation appears to be stabilized by this hydroxyl as it is close to a functional group on the agonist’s tail. Likewise, 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 may guide the flip and provide transient stability. The absence of a secondary nitrogen in ACh possibly explains its alternate re-orientation route, as only this ligand was observed to flip counterclockwise (Figure 5). It may also explain why the m2 basin in the PCA plots for ACh is especially broad and shallow (Figure 3).

The AC_LA→AC_HA transition compacts, dehydrates and stabilizes the binding cavity (Figure 4; Video 1). 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) (Figure 6-figure supplement 1). As reported elsewhere (Blum et al., 2010; Rahman et al., 2020; Zarkadas et al., 2022), 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 amine of δL121 in the complementary δ subunit backbone via a structural water (Figure 6). The agonist’s flip occurs at the start of m1→m2, and the establishment of this water-mediated H-bond completes in m2→m3 (Figure 6) (Video 1). This water-mediated bonding and the increase in VdW interactions appear to both be crucial for stabilizing the loops and agonist in the HA conformation of the neurotransmitter site.

Salt-bridge

Near the pocket, the αK145-αD200 salt bridge (Figure 1D) has been suggested to play an important role in receptor activation (Beene et al., 2002; Gay & Yakel, 2007; Mukhtasimova et al., 2005; Padgett et al., 2007; Pless et al., 2011). The hold rearrangements of αK145 are complex and agonist dependent. In AC_HA, this residue i) makes an apparent contact with αY93, ii) contacts the αY190 hydroxyl group (Figure 6-table supplement 1) and, iii) approaches loop F of the complementary subunit.

Interestingly, 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 and 0.46) (Indurthi & Auerbach, 2023). This agonist dependence in function has a parallel in structure. In the hold simulations, the αK145-αD200 distance (Å) increases from only 4.4 in apo to 6.1 or 5.1 in AC_LA (CCh or ACh). This spreading may relate to the agonist’s tail being inserted into the αY93 slot in the cis orientation. In AC_LA→AC_HA, for some agonists the αK145-αD200 separation (Å) shortens substantially to 2.8 or 2.6, making it 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 αK145-αD200 and αK145-αY190 distances can primarily be attributed to αK145 (Figure 6-figure supplement 2).

Discussion

After diffusing to the target, the only agonist-dependent rearrangements (energy changes) in receptor activation are in catch and hold, A-C⇄AC_LA⇄AC_HA (Figure 1B, bottom) (Auerbach, 2024; Jadey & Auerbach, 2012). By comparing binding energies calculated in silico with those measured in vitro we were able to identify in the MD trajectories the hold end states, AC_LA and AC_HA, plus an intermediate configuration not detected in electrophysiology experiments. Three conspicuous hold rearrangements discussed below are: i) a cis→trans rotation of the agonist (flip), ii) an up→down displacement of loop C (flop), and iii) reductions in pocket volume, water and fluctuations (fix).

There was reasonably good agreement between calculated and experimental binding free energies for ACh and CCh, but the match was excellent when free energy ratios (efficiencies) rather than absolute values were compared (Figure 4A). Absolute in silico values overestimated experimental ones by approximately the same factor (different for each ligand) that cancelled in the ratio. We do not know the origin of this factor but speculate that it is caused by errors in ligand parameterization. Regardless, the results show that agonist efficiency can be a useful metric for associating structure with function, and that efficiency can be estimated from structure alone.

Below we consider the structural elements associated with flip, flop and fix, namely the agonist (cis→trans), loop C (up→down) and the pocket (large→small, et al). MD simulations of the α–δ AChR neurotransmitter site suggest the following approximate and jittery generic sequence of rearrangements in AC_LA→AC_HA (Video 1). The stable states are bold.

1. ACLA (m1): agonist cis, loop C up, pocket wide/wet/wobbly,

2. …agonist leaves cis and starts to rotate about the N+-αW149 fulcrum,

3. …loop C starts to move down (towards the fulcrum).

4. Intermediate (m2): agonist half-rotated, tail stabilized by αY198,

5. …agonist rotates to trans (little change at the fulcrum), tail not secured in the pocket,

6. loop C is fully down.

7. ACHA (m3): agonist trans, tail secured by VdW contacts and H-bonds, pocket small/dry/rigid.

Agonist

With all agonists and in all simulations, the rotation of the agonist about the N⁺-αW149 fulcrum begins at the start of the AC_LA→AC_HA transition. Movements of bound ligands have been identified in other proteins. In AMPA-type glutamate receptors, the neurotransmitter glutamate binds in either of 2 poses, crystallographic or inverted (Yu et al., 2018). In one MD simulation these interconverted, suggesting the bound ligand is not rigidly fixed. Although the alternative poses were not associated with binding energies, the binding cleft was more closed in the crystallographic (presumably HA) versus inverted (possibly LA) orientation. The re-orientation of a bound glutamate in AMPA receptors is possibly related to the flip of the agonist in AChRs.

Multiple ligand positions in binding pockets have also been reported in Glycine receptors (Yu et al., 2014), 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, bound substrate movement is apparent in hexokinase (Bennett & Steitz, 1978) and alcohol dehydrogenase (Plapp, 2010). It is not unusual for movement of a bound ligand to trigger an otherwise unfavorable local rearrangement that drives protein function (an ‘induced fit’; (Richard, 2022). Not only can an untethered ligand come and go, it can also reorient itself dramatically to change binding energy and, hence, the local environment that couples to the full allosteric transition.

In AChRs, such ligand reorientation cannot pertain to all agonists. Some with high potency are perfectly symmetrical (decamethonium) and some lack a tail (tetramethylammonium). Although we have not yet identified the conformational dynamics that determine ΔG_HA, factors other than the flip certainly contribute, and in proportions that are likely to be agonist dependent (see below).

The PCA plots (Figure 3) show 3 structural populations, one of which (m2) occurs in time between m1 (AC_LA) and m3 (AC_HA) (Figure 3-figure supplement 2). In m2, the agonist’s tail (CCh, Ebt and Ebx) interacts with αY198, possibly to guide the clockwise rotation of the ligand and reduce the probability of a return to the cis orientation. The 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 (Figure 4-Source Data 1). We do not have an explanation for this difference.

Loop C

A number of observations place focus on the role of the downward displacement of loop C in receptor activation. In ACh binding proteins, its extent correlates with agonist potency (Basak et al., 2020; Du et al., 2015; Polovinkin et al., 2018) even if 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). Simulations of the apo AChR indicate that in the absence of agonists loop C has flopped partially, narrowing the aqueous connection between the bath and the pocket (Figure 2-figure supplement 2C).

Given the flop in loop C, it is puzzling that electrophysiology measurements show the agonist association rate constant (k_on) is substantially faster to C_HA (narrow entry, loop C down) than to C_LA (wide entry, loop C up), as well as more correlated with agonist potency (Grosman & Auerbach, 2001; Nayak & Auerbach, 2017). Docking results offer a clue to the apparent contradiction that narrowing an entrance speeds association up to the diffusional limit (4 x10⁹ M⁻¹s⁻¹). All 4 agonists dock in the cis orientation to C_LA (wide) versus in the trans orientation to C_HA (narrow) (Figure 2B). Hence, despite the wider aperture, trans is disfavored and lower affinity cis is preferred. We hypothesize that binding in the trans orientation is prohibited in the apo pocket, and that loop C flop (and other rearrangements) both narrow the entrance and remove this prohibition. Accordingly, k_on with loop C up (LA) is slower than with loop C down (HA). We do not know the structural element(s) that prohibits trans with a wide pocket. Any or all of the constellation of rearrangements in unliganded C_LA→C_HA could contribute.

As an aside, the inward displacement of loop C toward the pocket center resembles the ‘clamshell closure’ apparent at other receptor binding sites (Chen et al., 2017; Hansen et al., 2005; Hibbs & Gouaux, 2011). However, in light of the k_on experimental results, using ‘closure’ (or ‘capping’) to describe this motion is misleading because it implies putting a lid on the pocket that would reduce, not increase, the agonist association rate constant. For this reason, we use functionally neutral words ‘flop’ and ‘up/down’ for loop C displacement and position.

AChRs can turn on without a loop C. Truncation of loop C eliminates activation by ACh but has almost no effect on unliganded activation (Purohit & Auerbach, 2013). Interestingly, leaving only loop C residues αY190 and αY198 in place saves activation by ACh. Hence, repositioning of one (or both) of these in up→down is sufficient to maintain activation by the neurotransmitter. Other agonists were not investigated.

αY190 is our interest. An Ala substitution here results in a large loss of ACh binding energy with about half coming from the removal of the hydroxyl group (Purohit et al., 2014). In contrast, the mutations αY198F and αY93F have essentially no effect on ACh binding energy. These results suggest that repositioning αY190 consequent to loop C flop is important for activation by ACh (Figure 6). The αY190 ring fills the slot that opens next to αY93 when the agonist exits cis, preventing a return (Figure 5-figure supplement 1). However, the distance between the αY190 ring and N⁺ does not change significantly so it is unlikely that the deployment of αY190 increases binding energy via an interaction with the ring (Figure 6-table supplement 1).

Perturbation of the salt bridge residue αD200 by the αY190 hydroxyl has been proposed to trigger the AChR allosteric transition (Beene et al., 2002; Gay & Yakel, 2007; Mukhtasimova et al., 2005; Padgett et al., 2007; Pless et al., 2011). In support of this idea, the mutation αK145A reduces ACh (and choline) binding energy substantially but has no effect when the αY190 hydroxyl is absent (Bruhova & Auerbach, 2017). In the simulations, the downward displacement of loop C positions the αY190 hydroxyl to within ~5.0 Å of αD200 and αK145 (Figure 6-table supplement 1). We sought, but did not find, structural water in this gap. Although there was no correlation between these separations and the degree of loop C displacement, the gap was smaller with high- (ACh and CCh) versus low-efficiency agonists (Ebt and Ebx). In summary, some evidence supports the above proposal, but so far, the simulations have not illuminated a mechanism. Substrate-induced loop displacements are known to deploy residues to enhance catalysis, for example in chymotrypsin (Ma et al., 2005) and triosephosphate isomerase (Brown & Kollman, 1987; Nickbarg et al., 1988). Repositioning αY190 in the AChR AC_LA→AC_HA transition could be analogous.

Although attention on loop C has focused on the role of the flop, it is also important to consider the functional consequences of the up position. The structural elements of the up (wide cleft) configuration that prohibit agonist docking into the trans orientation have not been identified. A key requirement for allostery is that the ligand bind more strongly to one of the alternative conformations of the orthosteric site, so attention must also be given to the conformational dynamics of the ‘catch rearrangement, A-C⇄AC_LA. In both GABA_A (Jones et al., 2001; Mortensen et al., 2010) and ACh receptors (Jadey & Auerbach, 2012) the LA agonist association rate constant is correlated with potency rather than diffusion constant, indicating that the formation of AC_LA requires an energy (structure) change.

Pocket

The αW149-N⁺ hinge remains relatively unchanged in the hold transition (Figure 6, Figure 6-table supplement 1). Hence, we hypothesize the bulk of the extra favorable agonist binding energy in AC_HA is not generated here but somewhere else, cis→trans. For the agonists we examined, the rotation places the tail into the main binding cavity (Figure 6-figure supplement 1) (Video 1). The trans orientation is secured by VdW interactions (possibly caused by further downward loop C displacement and dehydration) and a hydrogen bond with a structural water bonded to the loop E backbone (in the complementary subunit). The net result to the overall pocket is a reduction in volume (Figure 4C) and increase in stability of loops C and F (Figure 3).

In the simulations, these two structural changes, distributed over the large surface area of the cavity, were robust and consistent for all 4 agonists. It is possible that the increase in binding free energy in hold derives mainly from restructuring of the pocket rather than from a change in the N⁺-aromatic interactions or, indeed, any one local ligand-protein interaction. Compact, hydrophobic binding pockets are associated with functional activation in other ion channels (Furukawa et al., 2005; Sigel & Steinmann, 2012) and observed in other proteins including β2-adrenergic receptors, opioid receptors, and HIV-1 protease (Huang et al., 2015; Rasmussen et al., 2011; Wlodawer & Erickson, 1993).

After hold, the next step In AChR gating is the rearrangement of the extracellular domain (Auerbach, 2024; Purohit et al., 2013). In related receptors, this domain becomes more compact and stable (a rearrangement called ‘unblooming’) (Sauguet et al., 2014; Taly et al., 2009), structural changes that echo those in hold at the AChR neurotransmitter site. It is therefore possible that the AC_LA→AC_HA rearrangement of the neurotransmitter site nucleates the subsequent restructuring of the extracellular domain. Extending this speculation, the structural correlates of energy change in gating might at first spread in contact area (from N⁺, to pocket, to extracellular domain) and then refocus (from transmembrane domain, to gate region, to pore water).

In the complementary δ subunit (Figure 1D), loop E (β5-β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). In addition, loop E residues are responsible for the binding energy differences α−γ versus α−δ (Nayak et al., 2016) and reduce site instability caused by loop C mutations (Vij et al., 2015). Although our simulations do not shed light on these matters, the temporal sequence flip-flop-fix suggests that perturbation of the loop E backbone could also play a role in downstream activation of the receptor.

Agonists belong to efficiency classes in which ΔG_LA/ΔG_HA is the same for all members. So far 5 such classes have been identified experimentally in adult-type AChRs, indicating that also are 5 pairs of AC_LA/AC_HA structural classes (Indurthi & Auerbach, 2023). Although in the simulations loops C and F movements, salt bridge distances were agonist dependent, with only 4 ligands we cannot associate any these rearrangements with efficiency classes. Additional simulations using more agonists and with binding site mutations may lead to better understanding of the dynamics that underpin molecular recognition in catch and protein animation in hold, and perhaps reveal the structural basis of their linkage (Auerbach, 2024).

Implications

Several aspects of this work may have general implications. The LA and HA end states of hold emerged from simulations of a single protein structure and were identified by comparing experimental with calculated binding energies. If general, this approach could be useful because cryo-EM structures of ligands bound with LA are lacking. Agonist efficiency was calculated accurately from structure alone. This, too, could be useful because efficiency allows affinity and efficacy to be calculated from each other (Indurthi & Auerbach, 2021). In AChRs the up→down movement of loop C flop does not impose a lid on the pocket. Rather, the up position serves to enforce catch (ligand recognition) by blocking directly binding in the trans orientation, and the down position serves to stabilize the AC_HA pocket and, perhaps, to deploy αY190. It is possible that the ‘clamshell closures’ apparent in other receptors also serve to activate for reasons other than imposing a steric barrier. Finally, it appears that the movement of the agonist, bound with low affinity, is the initial event in gating. If confirmed, we think it both fitting and beautiful that the motion of a structural element added to the system from the outside jumpstarts the allosteric transition.

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. Because agonist energy changes are approximately independent of energy changes in the protein outside the binding pocket, for the simulations we removed the transmembrane domain and simulated only α−δ subunit extracellular domain dimers. The starting structure for the hold transition was resting C (PDB: 6UWZ) after toxin removal and is essentially the same as the apo structure (Zarkadas et al., 2022) (Figure 1C). The 4 agonists we examined were acetylcholine (ACh), carbamylcholine (CCh), epibatidine (Ebt) and epiboxidine (Ebx) (Figure 2A). These are of similar size (MW 187, 183, 209 and 178, respectively) but have different efficiencies (0.50, 0.52, 0.42, and 0.46, respectively) (Indurthi & Auerbach, 2023).

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 residues 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. These agonists were chosen because they have a similar size and represent different efficiency classes (Indurthi & Auerbach, 2023). 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 AM1BCC 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 using 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; Figure 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

Experimental binding free energies were estimated by using electrophysiology, as described elsewhere (Indurthi & Auerbach, 2023). Briefly, the constants K_HA and L₂ were estimated by using electrophysiology (dose-response curves, for example), and K_HA was calculated using the thermodynamic cycle with L₀ known a priori. The logs of K_LA and K_HA are proportional to ΔG_LA and ΔG_HA. 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). 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).

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 was calculated as ΔG_bind=ΔG_complex−ΔG_protein−ΔG_ligand 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.