Structural correlates of affinity in fetal versus adult endplate nicotinic receptors

Adult-type nicotinic acetylcholine receptors (AChRs) mediate signalling at mature neuromuscular junctions and fetal-type AChRs are necessary for proper synapse development. Each AChR has two neurotransmitter binding sites located at the interface of a principal and a complementary subunit. Although all agonist binding sites have the same core of five aromatic amino acids, the fetal site has ∼30-fold higher affinity for the neurotransmitter ACh. Here we use molecular dynamics simulations of adult versus fetal homology models to identify complementary-subunit residues near the core that influence affinity, and use single-channel electrophysiology to corroborate the results. Four residues in combination determine adult versus fetal affinity. Simulations suggest that at lower-affinity sites, one of these unsettles the core directly and the others (in loop E) increase backbone flexibility to unlock a key, complementary tryptophan from the core. Swapping only four amino acids is necessary and sufficient to exchange function between adult and fetal AChRs.

E ndplate AChRs are heteropentamers that have two a(1) subunits and one each of b, d and either g or e. Each receptor has two functional neurotransmitter binding sites located in the extracellular domain at subunit interfaces, either ad þ ag (fetal) or ad þ ae (adult) (Fig. 1a, inset). Agonist affinities in mouse AChRs have been measured for these sites, both separately and as pairs [1][2][3][4] . The resting equilibrium dissociation constant for ACh (K d ACh ) is B30-fold lower at ag compared with ad or ae. The affinity of the fetal, ag-site is similar to that of the Lymnaea stagnalis acetylcholine binding protein (AChBP) 5 .
The fetal, g-subunit is required for the proper maturation of the neuromuscular synapse [6][7][8] . In mice, g-null mutations are lethal 9 and in humans g-subunit mutations cause both lethal and non-lethal (Escobar) types of multiple pterygium syndromes [10][11][12] . The reason(s) for the g-subunit requirement in synaptogenesis is not known, but possibilities include the higher agonist affinity, smaller single-channel conductance, longer open-channel lifetime, smaller gating-voltage dependence, lower Ca þ þ permeability and lower probability of opening constitutively of fetal AChRs 1,13,14 . Physiologically, the higher affinity of the fetal agonist binding site for both ACh and choline will lead to larger cellular responses 1,15 at the low concentrations of these agonists that prevail at developing neuromuscular synapses 16 .
There are B250 side chain differences between the adult eand fetal g-subunits. Mutational correlates of differences between the fetal and adult AChR function have been investigated previously. The mice expressing AChRs having an e-g-subunit chimera joined at a conserved glycine in loop E, g(r113) þ e(Z114), have a fetal-like open-channel lifetime and adult-like conductance 7 . Synapses expressing these AChRs undergo normal endplate differentiation but show altered innervation patterns. In mouse AChRs, the M3-M4 linker has been shown to influence fetal versus adult open-channel lifetime 17 . In Xenopus AChRs, residues at positions 6 0 /7 0 in M2 (N/I in g versus S/V in e) were implicated in setting fetal versus adult conductance 18 , and in rat AChRs, swapping the M2-20 0 amino acid from K to Q results in a partial exchange of conductance 19 . None of the above studies investigated the structural basis of fetal versus adult affinity, which is the topic we address here. At all three kinds of agonist site (ad, ae and ag), ACh is stabilized in the binding pocket by a core of five aromatic residues (Fig. 1a). On the basis of their individual contributions to affinity, these can be divided into two working groups, an 'aromatic triad' and a 'special pair' 20 . The triad is aW149 (indole ring), aY198 (benzene ring) and aY190 (benzene ring and hydroxyl). The four functional groups of these three amino acids each make a similar contribution to ACh binding energy at all the three agonist sites 15,21 . In adult AChRs, these groups provide almost all of the neurotransmitter binding energy, which is equal to þ 0.59 lnK d (23°C). Together, these groups stabilize ACh by approximately À 5.1 kcal M À 1 per site ( À 10.2 kcal M À 1 for ad þ ae sites combined). The action of the triad is led by aY190, because the ring and -OH each contribute approximately À 2 kcal M À 1 at ad, ae and ag.
A second working part of the agonist site apparatus is the special pair, aY93 (benzene ring) and e/gW55 (dW57; indole ring). At ae and ad, these two groups make only a small contribution to ACh affinity, but at the fetal ag site, the aromatic rings contribute significantly to bring the single-site-total to approximately À 7.1 kcal M À 1 and the ad þ ag total to À 12.2 kcal M À 1 (ref. 15). Moreover, at ag, the effects of alanine substitutions of this pair are not independent, as removing either ring interferes substantially with the other's ability to stabilize ACh. The energy contribution of this tryptophan differs massively between the fetal and adult sites. At ag, the gW55A mutation decreases affinity by B2,000-fold, whereas the homologous substitution at ae reduces it by only B13-fold and at ad it actually increases affinity by B2-fold.
Because the aromatic core has the same composition at all the three kinds of agonist sites, the difference in affinity between the fetal and adult AChRs can be attributed to residues in the complementary, non-a-subunit (d, e or g). Our approach to finding these amino acids was to use molecular dynamics (MD) simulations of homology models based on AChBP to identify complementary-subunit amino acids near the core that influence the contribution of the special pair to affinity. Then, we exchanged those side chains in vitro (d and e2g) and estimated affinity from single-channel currents of AChRs expressed in cells.
The results indicate that four residues in combination determine adult versus fetal resting affinity, three in loop E (111-113; g-subunit numbers) and one in the b5-b5 0 linker (104). In the b5 0 -b6 hairpin region, mutations of the human g-subunit (V108, S111, P112 and P121) cause multiple pterygium and Escobar syndromes 12 . The simulations suggest that these amino acids influence core properties by changing the structure and dynamics of the b5-b5 0 linker and the complementary b-sheet, to effect the action of the special pair and the affinity for the agonist. Swapping the four side chains is both necessary and sufficient to exchange fetal versus adult affinities and open-channel lifetimes.
The results are presented in five sections. First, we build and test homology models for the fetal and adult agonist sites. Second, we use the models to identify residues in the complementary subunit that influence ACh affinity in silico. Third, we use electrophysiology to measure in vitro affinities of AChRs having the identified residues swapped, fetal2adult. Fourth, we report the effects of point mutations on affinity, from the identified group and of a critical, complementary-side tryptophan. Fifth, we analyse the simulated structures and develop hypotheses for the mechanisms that undergird the fetal versus adult affinity difference.

Results
Tuning the model. Our approach was to use MD simulations to guide mutagenesis, and electrophysiology experiments to measure in vitro affinities. We used an empirical and computationally inexpensive method of estimating agonist binding energy in silico (see Methods section). By adjusting the free parameters for van der Waals and electrostatic contributions in the calculations, simulated and experimental affinities could be correlated (Fig. 1b). This was possible, in part, because in AChRs the entropy component of affinity is small 22,23 .
The circles in Fig. 1b are simulated versus experimental ACh binding energies for ag and ad agonist sites with mutations of core residues. The free parameters in the binding energy calculation (Equation 1, Methods section) were obtained only from the alanine mutations. Two non-alanine substitutions fell on the same regression line. The arrows indicate that the large and opposite effect of an alanine substitution of the complementary tryptophan on affinity at ag versus ad that is apparent in vitro was reproduced in silico. Also, in both simulations and experiments the ag WT site showed B40% more favourable binding energy for ACh and tetramethylammonium (TMA) compared with the ad WT site (Fig. 1c). These results suggest that brief simulations and approximate binding energy estimates can be used to screen for amino acids that influence affinity.
C4 in silico. The first goal was to identify residues in the complementary, gversus d-subunit that are responsible for the 30fold higher affinity for ACh in ag versus ad. We used two criteria to select candidate residues to be exchanged in silico. We chose side chains that in the model were within 20 Å of the quaternary ammonium (QA) group of ACh and were different between g and d but homologous between e and d. A total of 10 amino acids satisfied both of these criteria (C10; Fig. 2a).
In the first set of simulations, all the C10 side chains were swapped, g-d and d-g, and ACh binding energies were calculated. The results were that the C10-mutated ad site provided about the same binding energy as the ag WT site, and the C10-mutated ag site provided about the same as the ad WT site (Fig. 2b). This suggested that the mutation(s) we were seeking were within the C10 group.
We then winnowed the mutation list. First, each member of the C10 group was dropped one at a time. The result of these nineswap simulations was that the affinities no longer reversed when one of 4 of the original 10 amino acids was omitted. When g(L104, S111, P112 or D113) or d(Y106, Y113, D114 or S115) side chains remained as in the parent subunit, the full affinity reversal observed for C10 was incomplete. We call this set of four complementary subunit side chains C4.
Next, we simulated affinities using models in which only the C4 amino acids were exchanged as a group, to make constructs ag C4d and ad C4g (regular text is parent and superscript is target). Swapping all C4 side chains was sufficient to make the ag site have an affinity like ad WT and the ad site to have an affinity like ag WT . In silico, exchanging the C4 set of side chains was sufficient to exchange affinities, stably (Fig. 2b).
In the next set of simulations, each of the six other residues of the C10 starting set was added to C4, one at a time. None of these made a significant difference in the magnitude of the affinity swap. Finally, we dropped each mutation from the C4 set, one at a time, and estimated ACh affinity. Unlike in the above simulations, the trajectories of these three-mutant swaps were not stable (both fluctuating and divergent root mean square deviation; Supplementary Fig. 1). In silico, C4 was the minimum construct required for exchanging stably agand ad-binding energies.
C4 in vitro. We next made corresponding experimental measurements of ACh affinity (K d ) from single-channel currents recorded from AChRs expressed in HEK cells (Fig. 3). The C4 residues were investigated at ag-, adand ae-binding interfaces. In these electrophysiology experiments, the constructs were ad C4g First, the in vitro measurements were made using AChRs that had only one functional binding site, its companion being knocked out by a mutation(s) 24 . The results were that the C4 e-g and d-g exchanges resulted in a nearly complete swap of affinity (Fig. 3a). The ae C4g and ad C4g constructs each had the doseresponse profile and affinity of the target, ag WT site. The K d ACh values for ae C4g and ad C4g are shown in Supplementary Table 1.
The behaviour of ag C4d was more complex (Fig. 3b). This construct alternated between two distinct modes of single-channel activity. One mode retained the high-affinity Left-to-right: swapping all C10 residues almost completely exchanges affinity; dropping any one of four residues from C10 reduces the magnitude of the affinity change (C4; violet spheres in a); swapping all C4 side chains together is sufficient to exchange affinity; adding residues does not improve significantly the effect of the C4 swap.
(Supplementary Table 1) and longer open-channel lifetime characteristic of the parent (ag WT ) and the other had the lowaffinity and briefer open-channel lifetime characteristic of the target (ad WT ). The prevalence of each mode was approximately equal (isoenergetic), with a switching time constant of B20 ms. The C4 g-d affinity exchange in silico was stable, but in vitro the swap in binding and gating functions was bimodal.
In addition to decreasing the ACh equilibrium dissociation constant (increasing affinity), the ae/d C4g mutations generated a fetal-like, slower channel-closing rate constant in the presence of the neurotransmitter ACh. The effects of these C4g mutations in the absence of agonists or in the presence of saturating [ACh] are shown in Fig. 4a,b, with the results summarized in Fig. 5a and Supplementary Table 1. The C4 substitutions were effective in exchanging ACh affinity. We also measured affinity for the partial agonists carbamylcholine (CCh), TMA and choline and these exchanged affinity partially with C4 mutations. On average, the low-high affinity C4 swaps were B90% effective and the high-low affinity swaps were B70% effective.
We also made in vitro measurements from AChRs having only three of the C4 side-chain exchanges. The kinetic properties of these C3 receptors were complex and we did not attempt to estimate affinity.
The two WT AChR agonist sites act nearly independently in so far as binding energies measured from two sites combined are approximately the sums of single-site energies 15 . In the next set of experiments, we examined in vitro two different mutated, adult AChRs (both sites functional), ae C4g þ ad WT and ae C4g þ ad C4g . The prediction was that the first would have a fetal-type affinity and the second would be a 'super' AChR having the high, fetaltype affinity at both agonist sites. These expectations were confirmed, approximately (Fig. 5b,c).
Point mutations. The effects of point mutations of each C4 residue at each agonist site are shown in Fig. 5d. In experiments, none of these single exchanges had a large (41 kcal M À 1 ) effect on affinity or gating without agonists. With regard to affinity, the sum of the d-g single-residue C4 swaps was À 2.8 kcal M À 1 , which is close to that for the ad C4g combination, but the sum of the g-d single-residue swaps was only B þ 0.5 kcal M À 1 . With regard to the allosteric constant, none of the C4 swap sums were sufficient to account for the B14-fold smaller unliganded gating equilibrium constant ( þ 1.6 kcal M À 1 ) apparent in fetal versus adult AChRs at À 100 mV (ref. 1).
The mouse AChR e-g chimera was joined at a conserved glycine in loop E 7 . We were unable to record single-channel currents from AChRs having a point substitution here (gG114 to A, L or P) because these receptors either failed to express or open.
The action of the complementary, special pair tryptophan is the main distinction between ag versus ad sites 15,25,26 . The huge difference in the effect of an alanine substitution at W55 in g versus d in vitro was also apparent in silico (Fig. 1b). To further test both the model and the ability of C4 to swap the properties of the aromatic core, we compared simulated and experimental ACh binding energies for the C4 constructs in the presence of the point mutation g/eW55A or dW57A. If C4 was completely responsible for the different contributions of the special pair, we expected that an alanine substitution would have either no or little effect at ag C4d/e (as in ad/e WT ), but would cause a large decrease in affinity at ad/e C4g (as in ag WT ). The results agreed with these predictions, but partially ( Fig. 6). We did not test ag C4d ( þ dW55A) in electrophysiology experiments because ag C4d gave rise to complicated, heterogeneous responses without the tryptophan mutation (Fig. 3).
Analyses of structures. The modelling and electrophysiology results suggest that the b-sheet of the complementary surface of the agonist site can fold stably into, or isomerize between (Fig. 3b), either of two approximately isoenergetic configurations to produce a high, g-like or a low, d-like resting affinity (HA or LA, not to be confused with active-state vs resting-state affinities). Further, they suggest that the C4 side chains bias which of these super-secondary structures is adopted. Because there are no atomic-resolution structures of fetal or adult AChR agonist sites with agonist bound, we analysed structures generated in the MD simulations to explore possible difference between the alternate conformations.
The two regions of interest were the core and the complementary b-sheet.  estimated from C4-swapped constructs were similar to those from the parent subunits of corresponding affinity. Representative snapshots of the core from ag C4d and ad C4g trajectories are shown in Fig. 9a. In LA versus HA structures, the b5 0 linker residue Y104 is close and face on with aW149, and b2 residue W55 is distant and not orthogonal to either aW149 or aY93. The structural differences of the a-subunit aromatic triad were less pronounced, but in the LA structures these side chains a-subunit, white; complementary subunit, tan; core aromatics, green; C4 amino acids, violet; ACh nitrogen, blue sphere; aC of conserved loop E glycine, yellow. Bottom: high-resolution view of core in the parent (white) and C4-swapped (green) constructs. In the higher affinity (HA) structures, the binding pocket is more compact and W55 is closer to agonist (quaternary ammonium group, large blue sphere). were all further from the agonist's quaternary amine (QA). Overall, the simulations indicate that the HA core is more compact, organized and stable.
In AChBP there is a structural water within the core that is H-bonded to ACh and the backbones of b5 0 and b6 23,27,28 (Fig. 1a). In the simulations, both this water and ACh were dynamic (Fig. 9c). The water-b6 H-bond had the same propensity in both HA and LA structures, but the H-bonds with ACh and b5 were less prevalent in LA trajectories (Supplementary Table 2). In the simulations, the N-C-C-O dihedral angle (t2) of the neurotransmitter was bimodal, being either À 60 o or þ 60 o (Figs 8a and 9a). The þ 60 o configuration was approximately five times more common in LA trajectories.
We also estimated structural parameters for the complementary surface, which is made up of the b5 0 linker, the b5 0 -b6 hairpin (that includes loop E) and strand b2 (Fig. 1a).  Table 2). With LA amino acids (YDS), the 111-115 amide-oxygen H-bond is present so that the first residue is part of the hairpin rather than the loop. Also, in HA structures, the S111 side chain has an appropriate distance and geometry to H-bond with the loop E backbone (D113 or G114). The other difference in H-bonds is near the base of the hairpin, where the 109(N)-117(O) bond is less prevalent in LA structures. Figure 9b (bottom) shows that in the simulations, the overall hairpin is more upright (by B8 o ) and less twisted (by B9 o ; Fig. 8a) in LA structures.
There were fewer H-bonds between the W55 backbone and neighbouring residues in LA structures ( Fig. 9c and Supplementary Table 2 Finally, we used root-mean-square fluctuation (RMSF) analyses to estimate the dynamics of the complementary b-sheet backbone (Fig. 8b). The results are summarized in Fig. 9c. In LA structures, the RMSF of the b5 0 linker and loop E (that is, all C4 residues) was significantly smaller than in HA structures. This   Table 2). In summary, in the simulations, lower affinity was associated with a straighter hairpin, a lessdynamic b5 0 linker, a shorter and less-dynamic loop E, and a more-dynamic backbone near the core (base of hairpin, structural water and W55).

Discussion
It was possible to swap in vitro the affinities of the fetal and adult agonist sites based on mutations identified in silico. Apparently,  Figure 9 | Structures suggest a mechanism for fetal versus adult affinity. (a) Representative snapshots of C4-swapped cores (LA, low affinity ad WT /ag C4d and HA, high affinity ag WT /ad C4g ). Left, the LA core is spread, gW55 is far from ACh, C4 residue gL104Y (purple) is face-on to aW149 and ACh is bent. Right, the HA core is compact, dW57-aW149-aY93 are orthogonal, dY106L points away from the agonist and ACh is straight. (b) H-bonds of the b5 0 -b6 hairpin. Top, dotted lines are H-bonds. pink, bonds that correlate with affinity (are B2 Â more prevalent in simulations of LA or HA sites); black, bonds that are equally prevalent in LA and HA trajectories (Supplementary Table 2 From the two sets of gating rate constants in the bimodal construct ag C4d (Fig. 3b), we estimate that f ¼ 0.57, a value that is characteristic of many residues in the a-subunit transmembrane domain 33 . We do not know whether this low f-value can be attributed in general to complementary-side residues or to a difference between gand e-subunits. An indepth exploration of amino acid f-values on the complementary surface and in the ag core will be taken up, elsewhere.
The C4 mutation set does not account completely for differences in fetal versus adult AChR affinity. First, simulations showed a stable, high-to-low affinity exchange in ag C4d but in electrophysiology experiments, the behaviour of this construct was bimodal (Fig. 3b). Apparently, with this swap, the energy barrier separating the alternative, complementary b-sheet folds was low enough to allow a reversible isomerization on the B20 ms time scale, which is far beyond the simulation time frame. Second, the effect of C4 swaps was only partial with regard to the affinity of W55A (both in silico and in vitro; Fig. 6). Third, the affinities of the C4 constructs for the partial agonists choline and TMA were not completely swapped (Supplementary Table 1). These results suggest that the C4 group does not fully account for all differences in fetal versus adult binding properties. The addition of mutations to C4, either from the original C10 set or elsewhere in the extracellular domain, might improve the fetal versus adult match in function. Also, some of the differences between simulated and experimental results may arise from longrange interactions that were not modelled.
Although the C4 combination was sufficient to exchange affinities, we cannot be sure that this set of substitutions is unique in this regard. There are 1,023 possible mutant combinations for the starting set of 10, of which we examined only 84 in silico. It is possible that other combinations could also generate the affinity swap. Also, the original selection criteria could have left out important residues from the starting set. For example, gE57/ eG57/dD59 in the complementary b2 strand were not included in the starting set because the eand d-side-chains are not homologous. Nonetheless, the ability of the C4 combination to swap affinity and lifetime suggests that the C4 residues are an important basis for the functional differences between fetal versus adult endplate AChRs.
Recently, it was found that in a4b2 neuronal nicotinic AChRs, three complementary-subunit residues determine differential agonist potency of a4-a4 versus a4-b2 binding sites 34 . The corresponding positions in the endplate AChR g-subunit are L109, Y117 (in C10) and L119, none of which belong to the affinity-changing C4 group. This lack of correspondence may reflect a difference between these neuronal versus endplate AChRs, or may be traced to the benchmark of potency versus affinity. Nonetheless, the approach of combining simple and rapid simulations with in vitro energy estimates could be useful, in general, for revealing affinity-influencing amino acids in zones surrounding the core of ligand-binding sites in other ligand-gated ion channels.
In the simulations, the C4 mutations caused changes in structure that paralleled those that distinguished LA, ad WT and HA, ag WT sites (Fig. 8a). One of the C4 residues is in the b5 0 linker, adjacent to the core, and the other three are in loop E and far from the core. Hence, the effects of the C4 side chains appear to take place by both local and non-local mechanisms, to generate a core that is more compact in HA versus LA structures (Fig. 9a).
The primary features that distinguish the complementary subunit in LA versus HA structures were as follows. (i) The b5 0 linker is less dynamic, with the aromatic Y104 side chain close and face-on to aW149 and suggestive of a direct, p À p interaction. (ii) Loop E is one residue shorter. (iii) There are fewer H-bonds near the core, at the base of the hairpin with ACh and near W55 in b2. (iv) The W55 side chain is displaced from the core and not orthogonal to aW149 or aY93. (v) The loop E backbone is less dynamic and the b2 backbone near W55 is more dynamic. Overall, these differences suggest that three processes are required for C4 residues to generate a low affinity-unsettling the core by a local effect of the linker (104) side chain with aW149, loosening the b-sheet near the core by reducing the probability of H-bonds, and increasing the dynamics of the W55 backbone to unlock the special pair from the core.
We hypothesize that the dynamics and orientation of 104 are both determined by the aromatic versus aliphatic character of the side chain, which in LA is Y and in HA is L. The apparent pstacking of the dY104 and aW149 aromatic rings (Fig. 9a, left) likely affects core architecture directly. This interaction may influence that between aW149 and ACh (B0.7 kcal M À 1 less favourable at adult sites 15 ) and could also be a reason for the lower dynamics of the b5 0 linker. The unsettling of the core by the linker Y is, however, not sufficient to lower affinity completely because the effect of the point, L2Y mutation is relatively small (Fig. 5d), perhaps because gW55 remains locked in place.
Recently, we showed that the ad C4g mutations reduce modal changes in affinity induced by loop C mutations at the ad site 35 . It is possible that increased compactness of the site, along with the loss of a p À p interaction with aW149 stabilizes the aromatic triad and affinity.
In LA constructs the b5 0 -b6 hairpin is straighter, and loop E is shorter and less dynamic. It is possible that these differences affect interactions between core residues and the agonist directly, but more likely they influence affinity indirectly by changing the H-bond network of the b-sheet to increase backbone dynamics near the core. The simulations suggest that in LA structures, the bottom half of the b-hairpin and near W55, locations that have fewer H-bonds, are more dynamic (Fig. 9c). Our hypothesis is that the loop E C4 residues unlock indirectly the W55 backbone and allow this side chain to move away from the core, reducing its interaction with its partner aY93 and diminishing the contribution of the special pair towards affinity.
The simulations also showed differences in the ACh molecule between LA and HA constructs. In LA structures, the ACh t2 dihedral angle is þ 60 o rather than À 60 o , perhaps because the H-bond with water is less probable as a consequence of the greater backbone dynamics nearby. However, the extent to which the differences in agonist orientation and H-bonding relate to fetal versus adult affinity is uncertain. The agonist sites show a similar affinity difference for TMA (Fig. 1c, bottom), but this ligand does not have a rotatable bond and is too distant from the water to form a H-bond. We suspect that the core water is a factor that sets the higher affinity for ACh versus TMA, but may not be critical in determining the relative affinity of fetal versus adult sites.
The results suggest that both the structure and dynamics of complementary b-sheet are influenced by the C4 amino acids as a group, and that these residues are the bases, in part, for adult versus fetal endplate receptor function. At the agonist-site core, the orientation of the 104 side chain and dynamics of the W55 backbone together appear to generate the B30-fold affinity difference for the neurotransmitter. The simulations suggest that alternative conformations of the complementary b-sheet influence the architecture and affinity of the core. Further analyses of loop E interactions with loops A, B and C in the a-subunit (that hold other core aromatics) and with nearby residues in b2-b3 linker (the MIR) may further illuminate affinity mechanisms in AChRs.

Methods
Electrophysiology. AChRs were expressed in HEK cells by transient transfection of mouse endplate AChR subunits. Single-channel currents were recorded in the cell-attached patch configuration (23°C). The bath solution was (mM): 142 KCl, 5.4 NaCl, 1.8 CaCl 2 , 1.7 MgCl 2 , 10 HEPES/KOH, pH 7.4 and the pipette solution was: 137 NaCl, 0.9 CaCl 2 , 2.7 KCl, 1.5 KH 2 PO 4 , 0.5 MgCl 2 and 8.1 Na 2 HPO 4 , pH 7.3. To estimate the fully liganded gating equilibrium constant, a saturating concentration of agonist (100 mM; 4100 Â K d ) was added to the pipette. The agonists were acetylcholine (ACh), tetrametylammonium (TMA), carbamylcholine (CCh) and choline. To place the gating rate constants into a range suitable for kinetic analysis, we added background mutations that changed the allosteric constant or changed the membrane voltage. Neither of these perturbations had any effect on affinity 36,37 .
Experimental estimates of affinity and binding energy. Resting-state binding energy (in kcal M À 1 ) is equal to þ 0.59 lnK d (23°C). K d was estimated from single-channel current interval durations by using QUB software 38 . In one approach (Fig. 3a), K d was estimated from the ratio of dissociation/association rate constants, which in turn were estimated by fitting interval durations across multiple [agonist] using a A þ C2AC2AO scheme (one-site AChRs; A is the agonist, C is closed and O is open). In another approach (Figs 4 and 5d; Supplementary Table 1), K d was estimated from a ratio of gating equilibrium constants, as follows. In both WT and mutant AChRs, K d EE 0 /E 1 , where E is the C2O gating equilibrium constant and the subscript indicates the number of bounds agonist (E 0 is the allosteric constant) 39 . E 1 was estimated as the gating equilibrium constant for a single site at high (agonist; full saturation, 45 Â K d ). E 0 was estimated from the gating equilibrium constant in the absence of any agonists using a constitutively active background 40,41 . In the cross-concentration approach, the membrane was held at À 100 mV, and in experiments at high (agonist) the membrane was depolarized to þ 70 mV to reduce channel block by the agonist. Membrane potential has no effect on K d (ref. 37). When both approaches were used for the same construct, the K d estimates agreed.
Protein engineering. The mutations were incorporated into AChR subunits using the QuickChange site-directed mutagenesis kit (Agilent Technologies, CA) and were verified by nucleotide sequencing. Many AChR mutations away from the agonist sites only influence E 0 (ref. 33), which is 7.4 or 0.52 Â 10 À 7 in WT adult or fetal AChRs at À 100 mV (refs 40,41). We measured E 0 for every mutant construct by adding background mutations that increased it by known extents, but had no effect on affinity 36 . In selecting the backgrounds, we chose those that were energetically independent, so that the aggregate E 0 was the product of the WT value and the individual fold changes. The values shown in Fig. 5d have been corrected for the backgrounds.
To study AChRs having just one functional binding site, we added mutations to the e-, g-, d-subunits that eliminate binding at the mutated site 24 . To make agor ae-only AChRs, we added dP123R, and to make ad-only, we added e/gP121R, sometimes in combination with gW55R. These mutations also change E 0 , which was measured for each knock-out construct.
Homology model. The homology model of the extracellular domain of the fetaltype AChR was built based on the Aplysia californica ACh binding protein (AChBP) bound to epibatidine (pdb ID: 2BYQ) by using MODELLER 42 . The sequences of AChBP and AChR subunits were aligned using CLUSTALX 43 (Fig. 10a). The AChR subunits share B20% sequence identity with Aplysia AChBP. AChR subunits were modelled simultaneously so that spatial reciprocity was maintained at the interfaces. Residues 128 and 142 in the a-subunit and the corresponding residues in the other subunits were constrained as disulfide bonds. A protocol of conjugate gradient optimization, simulated annealing and molecular dynamics were used to refine the structure.
First, 100 structural models were generated. MODELLER has various assessment methods and objective functions to test the validity of a homology model (such as molpdf and DOPE scores), but these are not recommended to be used for multi-chain proteins 44 . Therefore, we used PROCHECK scores based on G-factor 45 to rank the models. The model with the best G-factor score was chosen for docking and simulations 46 . The top five models based on G-factor were similar in both structure (backbone root-mean-square deviation (RMSD) o0.7 Å; the equilibrated structure RMSD from the simulations was 1.2 Å) and scores. The selected model also ranked high in the MODELLER scores (first in molpdf and fifth in DOPE score), bad contacts (third) and Ramachandran criteria (third). Further minimization of the selected model reduced the bad contacts to zero.
Ligand docking. ACh and TMA were docked at the ad and ag binding sites using the Lamarckian genetic algorithm in AUTODOCK 47 . ACh and the core aromatic residues were allowed to be flexible, and 30 Å cubic search grid was used at the expected binding site with 0.375 Å grid spacing. Docked structures were analysed and selected on the basis of lowest energy and RMSD clustering. The CHARMM force field parameters for ACh and TMA were obtained from the Charmm Generalized Small Molecule Force Field webserver (CGenFF) 48,49 .
MD simulations. Point mutations were introduced in silico using the VMD mutator plug-in 50 . AChRs with ligands or mutations were optimized and equilibrated by using energy minimization and MD simulation. The system was solvated in a water box with TIP3P water model 51 and the box boundary was extended at least to 10 Å from the periphery of the protein in each dimension. Na þ and Cl À ions were added to neutralize the system and bring it to an ionic concentration of 150 mM NaCl.
Molecular dynamics simulations were run using NAMD 52 version 2.8, with the CHARMM27 force field 53 . First, a 20,000-step minimization was done using the steepest descent method, and with gradual release of restraints on the protein backbone. The system was heated to 300 K over 100 ps, a 500 ps equilibration run was performed in the NVT ensemble and then 20 ns MD simulations were performed in the NPT ensemble at a temperature of 300 K and pressure of 1 atm using the Nosé-Hoover method 54 . Following minimization, harmonic constraints (force constant ¼ 1 kcal M À 1 Å À 2 ) were applied on the Ca atoms of residues, which were 425 Å away from the ligand. These restraints maintained the global backbone conformation of the model while allowing relaxation of all side chains.
Periodic boundary conditions were applied. A 10 Å switching distance and a 12 Å cutoff distance were used for non-bonded interactions. The particle mesh Ewald method 55 was used to calculate long-range electrostatic interactions. The SHAKE algorithm 56 was used to constrain bond lengths of hydrogen-containing bonds, which allows a time step of 2 fs for MD simulations. Four MD simulation trajectories were obtained for each system. The coordinates of the systems were saved every 1 ps during MD simulations for later analyses. The protein RMSD became stable by the first 10 ns (Fig. 10b). All the analyses and binding energy calculations were done on snapshots extracted every 20 ps over the last 10 ns of each trajectory. The ensemble for each system therefore contained 2,000 snapshots of the system.

Calculation of affinity.
A total ligand-protein binding energy (DE) was estimated as described elsewhere 57 . In brief, this energy was calculated using a continuum solvent model from an ensemble of 2,000 snapshots: where E vdW is the van der Waals contribution and DE elec and is the electrostatic contribution calculated using the Poisson À Boltzmann (PB) method (PBEQ module of CHARMM 58 in a salt concentration of 140 mM). The molecular dielectric surface was defined by a probe radius of 1.4 Å. Previously optimized atomic Born radii 59 for the 20 amino acids were used to estimate the electrostatic free energy in explicit water molecules. The dielectric constant of the protein interior and the aqueous environment were set to 4 and 80, respectively. The coefficients a and b in equation (1) were evaluated by comparing simulated and experimental DE values for alanine mutations and optimizing the coefficients to reproduce the experimental energies (Fig. 1b). We calculated in silico binding energies for the WT (with ACh/TMA) and Ala mutations at each of the core aromatic residues (a93, a149, a190, a198 and d55), as described above. Values of a and b were scanned within a range of 0-2 to minimize the RMSE (root-mean-square error) between simulated and experimental energies for total 13 sample points: where n is each sample point and N ¼ 13.
The overall fitting quality was tested by the coefficient of multiple determination R 2 . The result was a ¼ 0. 22 15 . The small a is possibly due to the loss of protein-water VdW contacts in the AChR, whereas the larger b (41) may be due to cation-p interactions between the quaternary ammonium group of the ligand and the binding site aromatics, which is unaccounted for by the classical force field. With these values of a and b, the RMSE between simulated and experimental results were at a minimum (0.34 kcal M À 1 ; R 2 ¼ 0.92), with a correlation coefficient of 0.97. To test the generality of these parameters, we plotted simulated and experimental affinities for two additional, non-alanine mutations (dY106L and gL104Y) in Fig. 1b. These fell on the same regression line. The approximate binding energy calculated from MD simulations using the above method (DE) and the free energy estimated from single-channel currents ( þ 0.59 ln(K d )) are not equivalent. However, the entropy contribution to ligand binding energy has been experimentally shown to be negligible in case of muscletype AChRs 22 , and here we show that changes in these quantities can be compared for the purpose of engineering biological AChRs. In the text, we use 'affinity' for both the in silico and in vitro binding energy estimates.
Structure analyses. The geometric centres of the aromatic rings of interest and the ACh quaternary amine (QA) nitrogen were used as reference points. Structural analyses were done using the last 10 ns of each trajectory for the hetro-pentamer simulations. We measured all the distances and angles using VMD.
The volume of the ligand binding pocket was calculated by joining the centroids of the aromatic rings to form two adjoining tetrahedrons. The volume of each tetrahedron was calculated using the three-simplex determinant method from the coordinates of the vertices, as described elsewhere 15 . The b-hairpin twist was measured from the angle between the central axial vectors of b5' and b6. The central axis was defined by the least squares linear regression fit of the coordinates of the backbone atoms of residues 107-110 in b5' and residues 115-118 in b6 (g-subunit numbers). The regression fit was calculated by singular value decomposition of the coordinates.
To identify a hydrogen bond between two atoms (that is, acceptor and donor) a donor-acceptor distance of o3.5 Å and a H-donor-acceptor angle of r30°were used as criteria 60 . We used VMD to identify and calculate the occupancies of all H-bonds within the last 10-ns ensemble.
RMSD and RMSF analyses. To assess the conformational stability of the MD simulations, we calculated the RMSD of all backbone atoms in the aplus complementary subunit relative to the starting structures. The RMSD plots showed a small amount of drift over time in some of the trajectories (Fig. 10b, left). We determined that most of this drift came from the complementary-side loop F, which was highly flexible and disorganized. We repeated the RMSD analyses after removing a part of loop F from the calculation (residues 166-183 in g and 164-177 in d), and both the drift and variance in the RMSD were reduced significantly (Fig. 10b, centre). The plots show that the system stabilized within B10 ns. This region of loop F had been modelled ab-initio into the homology model because there is no corresponding region in AChBP. To test the role of loop F in affinity, we partitioned the binding energy (Equation 1) into contributions from individual residues using CHARMM. None of the loop F residues that were excluded from the RMSD calculations were in the top 5% (nB25) of those contributing to affinity, that comprise B98% of the total binding energy. The choice of 5% corresponds to a P value of 0.05.
As a further test of stability, we calculated a running average on the RMSD (without loop F) using a rolling window of 5 ns (Fig. 10b, right inset). The slopes of the regression lines fitted to the last 10 ns were small. The overall mean absolute slope value was 0.002 þ 0.001 Å ns À 1 ( þ s.d) for all the four constructs combined (n ¼ 16 trajectories). In addition, we calculated the average RMSD for the four trajectories for each construct (Fig. 10b, right; s.d. calculated from 5-ns bins). The drifts in these averages during the last 10 ns (B0.02 Å) were smaller than the s.d.s. of the fluctuations (B0.07 Å). These tests establish that all of the systems became stable within 10 ns and could be used for affinity estimation.
To compare the flexibility of ligand binding interface between the two states, we performed RMSF analysis on the backbone atoms based on the last 10 ns of the MD simulations. The RMSF of each residue was calculated with respect to the average structure of the ensemble using VMD version-1.9. This gives a residue-wise average of all possible fluctuations in the trajectories.