A mechanism for acetylcholine receptor gating based on structure, coupling, phi, and flip

Neuromuscular acetylcholine (ACh) receptors (AChRs) have two transmitter-binding sites (TBSs) in the extracellular domain (ECD) and a gate in the transmembrane domain (TMD). AChRs activate constitutively, but agonists bind with a higher affinity to the active-O(pen) versus resting-C(losed) conformation to increase the probability of being open (P_O) over the basal level (Monod et al., 1965; Karlin, 1967; Jackson, 1989; Auerbach, 2012). The complete C↔O allosteric transition (“gating”) is global, so rearrangements at the TBSs (that undergird affinity) and the gate (that undergird conductance) are linked.

AChRs have a molecular mass of ∼300 kD and share a common architecture with other pentameric ligand-gated ion channels (pLGICs; Fig. 1 A, left; Miyazawa et al., 1999; Sine, 2012; Taly et al., 2014; Cecchini and Changeux, 2015; Jaiteh et al., 2016). The ECD is mostly β-barrel (β1–10, in each subunit) with each TBS located midlevel at a subunit interface. The TMD is mostly α-helix (M1–4) with a hydrophobic gate located near the equator and formed by M2 from all five subunits. X-ray structures of other pLGICs show that in C→O there is a rotation between these two domains, a compaction of the ECD, a rearrangement of TMD helices that includes a radial tilt of M2 and straightening of M1, and dilation of the gate (Hibbs and Gouaux, 2011; Sauguet et al., 2013, 2014; Althoff et al., 2014). In AChRs, neurotransmitter affinity is determined mainly by a group of aromatic residues in loops at the TBS.

Actions at the ECD–TMD interface are important in the gating isomerization, where loops connect β-strands, the ECD–TMD, and the M2–M3 helices (M2M3; Fig. 1 B, bottom; Bouzat, 2012; Bertozzi et al., 2016; Morales-Perez et al., 2016). In mouse muscle AChRs, mutations of αM2M3 residues influence substantially the gating equilibrium constant and P_O (Grosman et al., 2000a; Bafna et al., 2008). Unnatural amino acid substitutions of a conserved M2M3 proline in 5HT3_A receptors led to the proposal that a cis-trans isomerization here opens the pore (Lummis et al., 2005), and simulations of other pLGICs support the idea that this proline is a “handle” that pulls open M2 (Zhu and Hummer, 2010).

There have been several suggestions regarding the mechanism of TBS–gate communication. In one, the TBS and the β1–β2 loop (loop2) at the base of the ECD are linked mechanically (by a rigid-body motion), so that agonist occupancy perturbs the β2 strand, then loop2, M2M3, the top of M2 (“pin-into-socket”) and the gate (Miyazawa et al., 2003; Kash et al., 2004; Lee et al., 2008). In another, agonist occupancy perturbs loop C and the β10 strand, then an interfacial salt bridge near the base of β10, loop2, M2, and the gate (Lee and Sine, 2005). A third proposal derives from phi values, which provide a snapshot of the relative free energy of an amino acid at the gating transition state (TS) on a scale from 1 (O-like) to 0 (C-like; see Materials and methods). In AChRs, diliganded gating α subunit phi values decrease in steps between the TBS and the gate (Fig. 1 A, right), which led to the proposal that in C→O the protein passes through a series of brief intermediate states, as a conformational cascade (Grosman et al., 2000b; Auerbach, 2005). A similar phi pattern has been observed in CFTR (Sorum et al., 2015). Structures of GLIC having an O-like ECD and a C-like TMD (Prevost et al., 2012; Sauguet et al., 2014) and an approximate microsecond shut component in high-resolution electrophysiology recordings (“flip” or “primed”; Lape et al., 2008; Mukhtasimova et al., 2016) are more-direct evidence for an intermediate state or states in pLGIC gating.

Basal activity of WT receptors is nearly zero and neurotransmitter application causes rapid and nearly complete activation, so it is understandable that in all of the above proposals the opening process starts with a perturbation at the TBS that propagates to the gate. However, some observations challenge this orthodoxy. WT AChRs expressed in muscle cells open and close without agonists (Jackson, 1984, 1986). The unliganded P_O of WT AChRs is small, but many mutations throughout the protein can increase this value substantially (as much as an agonist), apparently by making the same energy changes and with approximately the same phi values as when agonists are present (Purohit and Auerbach, 2009; Jadey et al., 2011). Furthermore, most mutations of TBS residues have little or no effect on the unliganded gating equilibrium constant, so this region is not special when it does not contain an agonist. Removal of loop C (β9–β10) or a key, complementary tryptophan on β2 does not affect unliganded gating (Purohit and Auerbach, 2013; Nayak and Auerbach, 2015). Finally, the top-down gradient in phi is not perfectly uniform. With agonists, some residues in αM2M3 have phi values as high as in the TBS (Bafna et al., 2008), and without agonists, gating phi values of some TBS residues become smaller and similar to those in the rest of the ECD (Purohit and Auerbach, 2010). These results question the assumptions that the AChR opening process starts at the TBS and that long-range, TBS–gate communication occurs via a discrete, mechanical pathway.

To investigate the AChR gating process, we used double mutant cycle analysis (MCA) of single-channel currents to measure interaction energies between αM2M3 with the TBS and proximal loops and linear free energy relationships (phi values) to sequence amino acid gating rearrangements. The phi pattern was used to compute energy landscapes of the gating TS ensemble (TSE), and simulations were used to reproduce the approximate microsecond flip/primed shut component. We propose a modified version of the conformational cascade model, in which C→O starts at αM2M3 and the TBS, to promote a twist of the ECD, a tilt of the TMD helices, a dilation of the gate, and the formation of a gate “bubble” that collapses to initiate ion conduction.

The paper is organized as follows. In the Materials and methods, in addition to describing experimental protocols, we provide summaries of the cyclic activation mechanism and phi value analysis. After a brief consideration of structures, we describe remote and local interactions with residues in the αM2M3 linker. We then show phi analyses that indicate αM2M3 is the first region of the protein to reach the TSE in unliganded gating and use simulations to show that flip reflects sojourns in the TSE mainly associated with gate bubble formation. We describe methods to compute a TSE energy profile from phi and combine phi and flip measurements to explore energy landscapes and structure-based models for AChR gating with and without agonists. Animations based on the models are shown at different time scales in Video 1.

Mutations were made to mouse muscle AChR cDNAs using the QuikChange site-directed mutagenesis kit (Agilent Technologies) and were verified by nucleotide sequencing. HEK 293 cells were maintained in Dulbecco’s modified Eagle’s medium supplemented with 10% (vol/vol) fetal bovine serum and 1% (vol/vol) penicillin-streptomycin, pH 7.4. The cell cultures were transfected with a mixture of cDNAs encoding WT or mutant subunits by calcium phosphate precipitation, by incubating for ∼20 h with 3.5–5.0 µg DNA per 35-mm culture dish. The AChR subunit cDNAs were added at a ratio of 2:1:1:1 (α/β/δ/ε or γ); cDNA for GFP (0.3 µg/µl) was also added to the mixture to facilitate the selection of AChR-expressing cells. The media was changed, and electrophysiology recordings commenced within 16–18 h.

Single-channel currents were recorded in the cell-attached patch configuration at 23°C. Pipettes were pulled from borosilicate capillaries and coated with Sylgard (Dow Corning). For the coupling measurements, an agonist concentration at least five times higher than the resting equilibrium dissociation constant was added to the pipette (≥1 mM ACh or ≥20 mM choline). To reduce channel block by the agonist, the pipette potential was −100 mV, which corresponds to a membrane potential (V_m) of ∼70 mV. The bath and pipette solutions were Dulbecco’s phosphate-buffered saline containing (mM) 137 NaCl, 0.9 CaCl₂, 2.7 KCl, 1.5 KH₂PO₄, 0.5 MgCl₂, and 8.1 Na₂HPO₄, pH 7.4.

Data acquisition and kinetic analyses were performed by using QUB software (Nicolai and Sachs, 2013). Currents were low-pass filtered at 20 kHz and digitized at a sampling frequency of 50 kHz. At high agonist concentrations, openings occurred in clusters that represent C↔O gating activity of individual AChRs with the long, silent intervals between clusters reflecting periods when all of the AChRs in the patch were desensitized. Clusters were selected by eye, and the intra-cluster currents were idealized into noise-free intervals after digitally low-pass filtering at 12 kHz and by using the segmental k-means algorithm. Rate constants were estimated from the idealized interval durations by using a maximum log-likelihood algorithm and a two state, C↔O model after imposing an approximate dead time correction of 25 µs. When necessary, an additional shut state was added to the model and connected to O to eliminate from the rate constant estimation gaps associated with sojourns in an approximate millisecond, apparently desensitized state (Salamone et al., 1999; Elenes et al., 2006). Without the inclusion of this state, the opening rate and gating equilibrium constants would be underestimated. The diliganded gating equilibrium constant was calculated as the ratio of the gating rate constants. Our time resolution was such that the approximate microsecond shut intermediate component (flip/primed) was not detected but rather was incorporated into C↔O TSE, represented by the arrow.

Methods for analyzing modal activity are described in detail elsewhere (Vij et al., 2015). In brief, clusters were selected by eye, and P_O was calculated from intra-cluster interval durations. In these experiments, which did not involve rate constant estimation, idealization was after low-pass filtering at 5 kHz and by using a half-amplitude criterion. A k-means clustering algorithm was used to determine the number of P_O populations and to segregate the clusters. An adult WT background was used for αP265A and αP265A+δW57A (100 µM ACh, V_m = −100 mV). For the loop C mutants, coupling was measured for each mode using [ACh] = 100 mM (to saturate the TBS) and V_m = 70 mV (to reduce channel block). To prolong the open duration and facilitate rate constant estimation, a background mutation was added to these constructs, βT456F for αY190F or αC418Y for αY198A. Without the αP265A mutation, there were no modes with either the backgrounds or the loop C mutations on these backgrounds.

AChRs activate by a closed allosteric cycle (no external energy) that involves C↔O gating with and without bound agonists and agonist binding to both C and O conformations (Fig. S1). In adult mouse AChRs, the two binding sites are equivalent and independent for ACh and choline (Nayak et al., 2014), so

E₂ is the diliganded gating equilibrium constant, E₀ is the unliganded gating equilibrium constant, and K_d/J_d is the ratio of equilibrium dissociation constants to C versus O (the coupling constant). In WT adult mouse AChRs at −100 mV (23°C), E₀ is 7.4 × 10⁻⁷, K_d^ACh = 170 µM, E₂^ACh = 25, and the coupling constant is ∼5,800 (Auerbach, 2010). The log of an equilibrium constant is proportional to the free energy difference between the end states. Taking the log of Eq. 1,

where ΔG₂ is the diliganded O minus C energy difference, ΔG₀ is the unliganded O minus C energy difference, and ΔΔG_cc is the log of the coupling constant and is the per-site difference in agonist binding free energy (bound minus free) to O versus C. The presence of intermediate states between C and O does not affect these energy differences. In WT adult mouse AChRs, at −100 mV, ΔG₀ is 8.3 kcal/mol and ΔΔG_cc is −5.1 kcal/mol for ACh and −3.3 kcal/mol for choline. For the TBS residues we studied, the effect of a mutation on ΔΔG_cc was much larger than on ΔG₀, so the pairwise interaction values involving these amino acids are mainly with regard to ΔΔG_cc.

In about half of the experiments, a background mutation was added to adjust E₀ so that the emergent interval durations with agonists were in a readily measurable range (0.1–10 ms; Jadey et al., 2011). The background mutations only influenced E₀ (ΔG₀) and had no effect on the coupling constant. ACh and choline provide different energies from the low to high affinity change (ΔΔG_cc values), which when added to the engineered ΔG₀ produced diliganded gating equilibrium constants (ΔG₂ values) that placed P_O into the optimal range for analysis (Eq. 2). By correcting the apparent rate constants for the effects of depolarization and the background mutation or mutations, gating equilibrium and rate constants were normalized to a standard condition (WT, −100 mV, 23°C). This approach eliminates many of the technical difficulties associated with idealization and modeling of events close to the time resolution of the patch clamp and extends greatly the range over which gating constants can be estimated. The backgrounds are shown in Tables S1 and S2.

The main method we used to estimate the interaction between side chains is MCA. A pair of residues was selected and mutated. The fold change in the diliganded gating equilibrium constant (E₂^mut/E₂^WT) was measured separately for each mutant. Multiplying the natural logarithm of the fold change by −RT (where R is the gas constant and T is the absolute temperature; RT = 0.59 at 23°C) gives the effect of the mutation as a free energy change (kcal/mol). All MCA experiments were with fully liganded AChRs.

Next, the fold change in E₂ was measured for the mutation pair. A deviation of the sum of the individual energy changes from that of the pair indicates an interaction ("coupling") between the two mutants. The coupling values are expressed as pair minus sum. For αP265 coupling, only the predominant, highest P_O mode was used, except for with αY198A and αY190F, for which two modes were measured. The coupling values were all from mouse endplate AChRs and either are new (n = 79) or were published previously (n = 98; Cadugan and Auerbach, 2010; Purohit and Auerbach, 2010; Jadey et al., 2011; Jha et al., 2012; Gupta et al., 2013).

We estimate that we are able to measure an approximately twofold change in E₂ (∼0.4 kcal/mol). Assuming independent errors, by the law of propagation of errors (Taylor, 1997), the combined error for a mutation pair is √(0.4² + 0.4²), or ∼0.57 kcal/mol. Another way to estimate the error in the coupling energy estimate is from the SD of experimental values for residues separated by >15 Å (excluding αM2M3), which was 0.56 kcal/mol. The agreement indicates that the errors in the equilibrium constant estimates are independent. Below, we set the threshold for an interaction being significant at ±0.6 kcal/mol.

An independent method of estimating interactions between the agonist and a mutation is by comparing mutational fold changes in E₂ and E₀ (Eq. 1). The two TBSs are approximately equivalent for the agonists we used, ACh and choline (Jha and Auerbach, 2010; Nayak et al., 2014), so √(E₂/E₀) is equal to the coupling constant. A deviation in ΔΔG_cc ≥0.6 kcal/mol compared with the WT value indicates an interaction with the TBS.

A mutation that changes the gating equilibrium constant causes the relative free energy (structure) of O versus C to be different than in the WT. This happens because the substituted side chain is more (less) stable in the local environment compared with the WT. So far, we have not found an AChR mutation away from the TBS that changes significantly the ΔΔG_cc of any agonist or that has an effect on ΔG₂ that is agonist dependent. This suggests that in AChRs, the agonist and separated mutations have mainly independent and short-range effects. Consequently, gating energy (structure) changes can be localized to the mutated amino acid and put into sequence by using phi value analysis, as follows.

Phi is the slope of a log-log plot of the full C→O (opening) rate constant versus the full C↔O (gating) equilibrium constant for a series of mutations of one residue. This slope gives the relative energetic character of that amino acid at the gating TS on a scale from 1 (O-like) to 0 (C-like). The opening rate constant is proportional to the free energy difference between the gating TS and C (ΔG_TS − ΔG_C), and the gating equilibrium constant is proportional to the free energy difference between O and C (ΔG_O − ΔG_C). Phi is the ratio of these free energy differences. Phi = 1 (the two energy differences in the previous sentence are equal) indicates that ΔG_TS and ΔG_O changed to the same extent with all mutations of that amino acid, and phi = 0 indicates that ΔG_TS and ΔG_C changed to the same extent. Energy is related to structure, so a residue with a high phi is structurally more O-like at the TS than one with a low phi. Phi gives the relative position in the forward isomerization at which a side chain changes energy (structure) from C to O (see supplemental text). Methods for determining and phi values and population means (the color assignments in Fig. 1 A, right) are described in detail elsewhere (Purohit et al., 2013).

Unliganded gating phi values of αM2M3 amino acids were estimated after adding background mutations that increased constitutive activity (Table S2). The level of activity of αP265 mutants was extremely low, and many patches were required before observing single-channel clusters, perhaps because of low expression.

AChR diliganded gating phi values are modal (mean values of ∼1, 0.8, 0.6, 0.3, and 0) and follow a longitudinal spatial gradient (Purohit et al., 2013). There are several interpretations regarding the physical meaning of a fractional phi value. (1) In the “multiple intermediates” view, each phi population represents a short-lived gating intermediate state (see supplemental text; Auerbach, 2005). The five AChR phi populations thus reflect a sequence of transitions between six states: C, four short-lived intermediates (which together comprise the TSE), and O. Higher phi residues change energy (structure) earlier in this reaction chain. (2) In the “multiple paths” view, there are two gating reaction trajectories connecting C and O, with TSs at phi = 1 and 0. A fractional phi value thus reflects the relative probability of taking either pathway (Ternström et al., 1999; Oliveberg, 2001). Phi still gives the position in the reaction where the mutated amino acid reaches the TS, but possibly with no metastable intermediate states between C and O. (3) In the “multiple perturbation” view, an amino acid can experience an energy perturbation at more than one point in a single reaction trajectory (Jha et al., 2009). For instance, a TBS residue might be bumped in the low→high affinity change and again when the ECD compacts, or an M2M3 linker residue might be perturbed at the onset of opening and again when the TMD helices rearrange. The emergent phi of an amino acid is thus a weighted mean of the energy perturbations at multiple positions in the reaction trajectory. In the case of a reaction with one metastable intermediate state separated from the end states by barriers at phi = 1 and 0, equal energy perturbations in both transitions would yield phi ∼0.5. These three interpretations are not mutually exclusive.

Of these three interpretations, only the multiple intermediate view predicts intermediate states and modal phi values that are organized spatially. The multiple path and multiple perturbation views (that only have barriers at phi = 1 and 0) both predict a continuum of fractional phi values and say nothing regarding spatial organization. We therefore use the multiple intermediate interpretation, with four metastable TSE states (T_1–4) separated by five microscopic energy barriers (the arrows). The relative energies of the barriers in the TSE were calculated from the diliganded gating phi values as described in the supplemental text and elsewhere (Zhou et al., 2005).

According to this interpretation, phi values reveal the relative heights of the TSE energy barriers in the reaction chain but nothing regarding the energy wells. To estimate these, we incorporated knowledge regarding flip/primed, assuming that this signal arises from sojourns in one or more of the TSE states (T_1–4, in brackets in Scheme 1). We simulated single-channel currents using QUB and Scheme 1 assuming all of the TSE states are nonconducting and using different values of the unconstrained rate constants while maintaining correspondence with the phi values. The objective was to match the apparent lifetime of flip/primed. We used parameters that mimic a high-resolution patch clamp experiment (sampling frequency of 1 MHz and low-pass filtering at 30 kHz) with idealization by using a half-amplitude criterion. After simulation using Scheme 1, the idealized interval durations were fitted using QUB using Scheme 2 after invoking a dead time of 8 µs, where sojourns in C′ correspond to those in the T states of Scheme 1 and the flip/primed shut component. In experiments with human adult AChRs, the lifetime of C′ is ∼9 µs, with left and right exit rate constants of ∼82 and 19 ms⁻¹ (Mukhtasimova et al., 2016).

In the Scheme 1 simulations, the rate constants were changed to adjust the depths of the TSE wells while keeping the relative barrier heights (set by the phi values) constant. For each simulation, we calculated an energy landscape from the rate constants (k) from the relationship ΔG^‡ = −0.59ln(k/A*) (in kcal/mol), using A* = 2 × 10⁶ s⁻¹ (Jones et al., 1992; Chakrapani and Auerbach, 2005). Changing A* alters the absolute, but not relative, heights of the barriers and therefore does not affect phi values. As described in the supplemental text, we also calculated a committor (‡), which is the point in the trajectory where there is a 50% probability of rapidly entering C or O (Bolhuis et al., 2002; Zhu and Hummer, 2012), and a transmission coefficient, defined as the inverse of the mean number of times the system crosses ‡ before reaching the other end state.

Distances were measured using GluCl as the structural template (PDB accession code 3RIF) using VMD (Humphrey et al., 1996). In the MCA analyses, separations were calculated as distances between α carbons (αC). The secondary structural assignment, loop versus helix, was by the measurement of the backbone dihedral angles and hydrogen bonds using the DSSP program (Kabsch and Sander, 1983; Touw et al., 2015). The conformational energy of the M2M3 backbone was calculated as the combination of bonded energy terms (bond, angle, Urey-Bradley, dihedral, and improper energy) using the CHARMM27 force field with CMAP corrections (MacKerell et al., 1998; Brooks et al., 2009). The O and C structures were GluCl 3RIF and 4TNV (residues 265–276; Hibbs and Gouaux, 2011; Althoff et al., 2014), and GLIC 4HFI and 4NPQ (244–255; Sauguet et al., 2013, 2014).

The supplemental text describes the method for computing the TSE energy landscape from phi values. Fig. S1 shows the allosteric activation cycle, and Fig. S2 shows energy landscapes for simulations with stabilized intermediate states. Video 1 describes the model and shows cartoon structures, energy landscapes, and single-channel currents for simulations of gating at different time scales. Tables S1 and S2 report the background mutations used to make the measurements. Tables S3 and S4 give the simulation parameters and results.

Classically, the global C↔O allosteric isomerization is called "gating," with "opening" as the forward process and "closing" as the backward process. The detection of a shut gating intermediate (C↔C′↔O; Scheme 2) has suggested a change in this nomenclature, with the first step called "flipping" (or "priming") and the second called "gating." Below, we present evidence that there are four intermediate shut states between C and O and name each of the five connecting steps according to a correlated structural change. In this paper, we will continue the classical use of "gating" to refer to the full transition between the stable end states of the reaction. For example, here "closing" is O→C′→C (where C′ represents all intermediate shut states) and is not necessarily the same as a step decrease in single-channel current that could reflect O↔C′↔O.

Each TBS is at an α subunit interface, ∼30 Å from the ECD–TMD interface and ∼50 Å from the gate (Fig. 1). The principal side (α subunit, in heteromeric AChRs) of each TBS is formed by loops A (β4–5), B (β7–8), and C (β9–10). At the two adult sites (αδ and αε), three aromatic amino acids provide most of the ACh binding energy, αY190 (loop C) > αW149 (loop B) ≈ αY198 (loop C), but at the fetal αγ site, the energetically coupled pair αY93 (loop A) and γW55 also contribute (Nayak and Auerbach, 2013). Mutations of two TBS glycines and two complementary side prolines also influence resting affinity (Sine et al., 1995; Ohno et al., 1996; Purohit and Auerbach, 2011; Gupta et al., 2013; Jadey et al., 2013; Nayak et al., 2016).

Fig. 1 A (right) shows a map of the relative free energy (structure) of residues in the α subunit at the diliganded gating TS (Purohit et al., 2013). The highest, phi ∼1 residues (purple) are at the TBS and in αM2M3 (αI260, αP265, and αS268). The ECD is mostly phi ∼0.8 (blue), and the TMD is mostly phi ∼0.6 (green), with a discrete border near a proline-induced π-helix in M1 (Purohit et al., 2015). The lowest, phi ∼0.3 residues (red) are at the gate and lipid-facing residues in M3 (not depicted). At the TS, the TBS and some residues in αM2M3 are O-like, the gate is C-like, and the ECD and TMD are intermediate. The phi map indicates that in the C→O conformational change, the α subunit TBS+αM2M3, ECD, transmembrane helices, and gate reach the TS in that order.

In mouse AChRs, the high-phi and large-effect amino acids in the TMD align with the αM2M3 region of α₄β₂ AChRs (Fig. 1, B and C; Morales-Perez et al., 2016). In the related pLGICs GluCl and GLIC, in O versus C structures, the M2M3 linker is displaced upward relative to M1 (Fig. 1 B, bottom inset). In both of these pLGICs, the calculated backbone conformational energy in C is 4 and 2 kcal/mol less stable than in O. We also calculated the nonbonded interactions of these residues with their local environment, and for both receptors, the C configuration was 5.2 kcal/mol less stable. These calculations did not consider side chain energies but suggest that the M2M3 backbone may be locally more stable in O versus C.

Interaction free energies were estimated by using MCA as the free energy difference between the effects of a mutation pair on the gating equilibrium constant minus the sum for the two mutations separately (Fig. 2 A). Below, we use “coupling” specifically to refer to the MCA interaction energy, “significant” to characterize an absolute interaction 0.6–1.2 kcal/mol, and “strong” for one that is >1.2 kcal/mol (see Materials and methods).

We used MCA of diliganded AChRs to estimate coupling between mutations of four αM2M3 residues (to A or G) and five TBS aromatic residues (Fig. 2, B and C). αS268A interacts strongly with αY190F and significantly with αW149A, two mutations that have large (∼2 kcal/mol) effects on ACh affinity in adult AChRs. αP265A showed strong interactions with all five aromatics. A glycine mutation of αL263 interacts significantly with the TBS aromatics. Overall, 74% of the αM2M3–aromatic pairs were coupled significantly, with all but one producing a smaller effect than predicted from independence.

Fig. 2 C (inset) shows coupling values for these αM2M3 amino acids and prolines or glycines at the TBS. 11 of the 14 residue pairs showed a significant degree of absolute coupling. αP265A was coupled significantly with αG147, εP121, and γP112, a residue in loop E that influences affinity. There was a strong interaction between αS268P and αY190P. For the glycine and aromatic TBS substitutions, the coupling values pertain to mutations in two α subunits, whereas for the proline mutants only one non–α subunit was substituted. Hence, the mean absolute coupling with the lone P121 (0.8 kcal/mol; n = 7) may indicate a stronger per-residue interaction than for the double glycines (0.9 kcal/mol; n = 7) and aromatics (1.1 kcal/mol; n = 22).

All of the coupling energies are summarized as a function of separation in Fig. 3 A. αM2M3 and the TBS interact more strongly than expected given their ∼30-Å separation. Also, many closely apposed residues do not interact significantly. Of the four αM2M3 residues we examined, αP265 and αL263 were the most interactive. At the TBS, strong interactions were divided approximately equally between the five aromatics.

In WT AChRs, kinetic properties of clusters are essentially homogeneous, but loop C mutations (in particular at αP197) generate discrete kinetic modes arising from differences in affinity at the αδ TBS (Vij et al., 2015). Fig. 3 B shows that modes are also apparent with the αP265A mutation in M2M3. As was the case with αP197 mutations, the αP265A modes disappear when δW57A is added, suggesting that these, too, arise from affinity differences at αδ. However, the αP265A and αP197A modes are not identical. αP265A modes persist with the addition of either αY190F or αY198A, whereas with αP197, the addition of αY190F but not αY198A reduces modal activity. Also, the addition of δW57A to αP265 causes a large reduction in P_O but has a smaller effect with αP197A. Regardless, the modal activity of αP265A independently supports the existence of long-range interactions between αM2M3 and the TBS. The mechanism for this connection is not known, but given the significant MCA coupling between αP265A and both εP121 and γP112, this could possibly involve long-distance backbone interactions between αM2M3 and the β-sheet of the complementary subunit at the TBS (Nayak et al., 2016).

Fig. 4 A shows coupling values for alanine mutations in αM2M3 (except αA270G) and nearby regions at the ECD–TMD interface. Regarding αloop2, strongly coupled pairs were Q48–I260/V261 (in the M2 helix) and E45–P265. Q48 and E262 do not interact, and neither of the phi ∼1 M2M3 residues (P265 and S268) interact with V46. Only a few αloop7 pairs were measured, with the only significant coupling being between V132–S268. Regarding εloop9, previous experiments showed there is strong coupling (−1.7 kcal/mol) between εG183I–αP265G (Jha et al., 2012). The new measurements show significant coupling between εG183 and T267–S268, but not with S266. Regarding εpreM1, εR219 interacts strongly with both S266 and T267, but εL221 was not coupled.

We also examined four alanine mutant pairs within αM2M3. αI260–αP265 are coupled by −0.6 kcal/mol, but the interactions between αI260–αS268, αS268–αK276, and αP272–αK276 were insignificant.

Overall, 10 of 29 local pairs showed strong or significant MCA coupling. For the 10 residue pairs that were coupled significantly, the absolute interaction free energy was 1.5 ± 0.2 kcal/mol. There were strong intra-subunit interactions with loop2 and cross-subunit interactions with εloop9 and εpreM1. αP265 was the only M2M3 amino acid coupled to three different structural elements (αM2, αloop2, and εloop9).

Fig. 4 B shows cross-subunit coupling at the diliganded fetal (αγ) TBS. Previously, it was found that in adult-type AChRs, εP121A interacts strongly with αG147A (2.6 kcal/mol) when ACh is present (Gupta et al., 2013). At the fetal site, γP121A interacts strongly with the three aromatic amino acids that provide most of the ACh binding energy from the α subunit (αY190 and αY198 in loop C and αW149 in loop B) and significantly with αG153. γP121A did not interact with its neighbor γP120A, αG147A, or the pair of residues that provide substantial binding energy at the fetal site and are themselves coupled (αY93 and γW55).

The MCA results indicate that some long-range interactions between αM2M3 and the TBS are nearly as strong as short-range ones at the domain and TBS interfaces. In the next set of experiments, we compared phi values at αM2M3 and the TBS, with versus without agonists.

Fig. 5 A shows unliganded gating rate-equilibrium linear free energy relationships (REFERs; slope is phi) for four αM2M3 residues. Without agonists, these amino acids have a mean phi = 0.89 ± 0.05 (mean ± SEM), which is the same as with agonists (0.90 ± 0.03; Fig. 5 B). However, TBS amino acid phi values without agonists (0.79 ± 0.03) are smaller than with agonists (0.95 ± 0.02; Purohit and Auerbach, 2010). The agonist, which we think of as a removable TBS side chain, has a phi value of 0.93 (Grosman et al., 2000b). These results indicate that at the unliganded gating TS, αM2M3 is more O-like than the TBS.

We attempted to measure diliganded gating phi values for residues in the β10 strand of the α subunit. However, the effects of 14 mutations on the gating equilibrium constant were all small (<0.5 kcal/mol), and phi values could not be estimated (Purohit et al., 2013). The mutations were αT202A/C, αY203A/C, αH204A/C/Y, αF205A/C/T, and αV206A/C/I/W.

In AChR diliganded gating, amino acid phi values are distributed as five Gaussians with means of ∼0.95, 0.79, 0.58, 0.33, and ∼0.06 (Fig. 6 A). The phi populations also show a clear spatial organization, with the first four mapping to the α subunit TBS+M2M3, ECD, TMD, and gate, respectively (Fig. 1 A), and the last mapping mainly to gate residues in non–α subunits. (Lipid-facing residues in αM3 also have a phi ∼0.3.) We used Scheme 1 to model the intermediate states of the gating TSE, with the ratios of the exit rate constants from each of these (r values) calculated from phi (supplemental text) so that 6 of the 10 rate constants were constrained. When the overall C↔O equilibrium constant is 1, Scheme 1 becomes Scheme 3. The rate constants are ms⁻¹, and the bracketed states (T_1–4) comprise the intermediates of the TSE. Each arrow in Scheme 3 is associated with a microscopic barrier crossing representing a structural change that can be localized spatially according to the phi map.

We hypothesize that the αM2M3 and TBS rearrangements are synchronized because of coupling, so that the first step (C↔T₁) incorporates gating movements in both regions. The subsequent transitions in Scheme 3 are rearrangements of the ECD (T₁↔T₂), the TMD (T₂↔T₃), the gate (T₃↔T₄), and water entry into the gate region (T₄↔O). An energy landscape corresponding to Scheme 3 is shown in Fig. 6 B (right). In the opening process, the relative TSE barrier heights (calculated from the phi values) show an increase associated with ECD movement, followed by a nearly flat region. The committor (‡; the position in the TSE that leads rapidly to C or O with equal probabilities) is near T₃ (Fig. 6 B). Because of the upward tilt of the TSE barriers in opening, when C→T₁ is slow relative to exits from the TSE, a mean of ∼20 exits from C must be attempted before O is achieved (supplemental text). Because of the downward tilt of the TSE in the reverse direction, closing is more efficient with only approximately three exit attempts from O required for a complete transit to C.

Next, we adjusted the unconstrained rate constants of Scheme 3 while maintaining consistency with the phi distribution. The goal was to set the TSE well depths so that the lifetime in T₁–T₄ was consistent with flip/primed. In experiments, this duration was ∼9 µs and arose from left and right exit rate constants from C′ (Scheme 2) of 82 and 19 ms⁻¹ (Mukhtasimova et al., 2016). Fig. 6 B (left) shows example simulated single-channel currents and, below, corresponding open and shut interval duration histograms.

First, we used a single value k for all of the unconstrained rate constants (in Scheme 3: k₃, k₅, k₇, and k₉ all equal to 300 ms⁻¹). This kinetic scheme predicts a broad distribution of four brief-shut components. Fig. 6 C shows the results of filtering and fitting by Scheme 2. In facsimile high-resolution patch-clamp recordings, the four components merge to become the tail of a single exponential with a lifetime of τ^flip ∼6.6 µs. Fig. 6 B shows the corresponding energy landscape and that ‡ is right of center.

Fig. 6 D (and Table S3) shows the simulated TSE and flip time constants and the ratio of the C′ exit rate constants estimated by using Scheme 2 as a function of k. Values of k <200 ms⁻¹ generate τ^flip values that are too long-lived to be consistent with the experimental results.

In the next experiments, we used a large value of k (500 ms⁻¹) to make most intermediate-state sojourns extremely brief and made just one or two of the TSE wells deeper (Fig. S2 and Table S4). These schemes approximate a single intermediate shut state while maintaining consistency with the phi pattern. It was possible to reproduce an ∼7-µs shut interval by lowering any one of the wells by 0.5–1 kcal/mol. Stabilizing T₁ and T₂ together, each by 0.5 kcal/mol, also produced a flip-like event. Apparently, there are many possible combinations of k values and energy wells that can generate an approximate microsecond intermediate shut event, but we did not attempt to simulate these exhaustively.

In the simulations, the lifetime of each intermediate state can be calculated as 1/(k_left + k_right). Using Scheme 3 and with a uniform k = 300 ms⁻¹, these lifetimes are 0.8, 1.4, 1.5, and 1.4 µs for T₁–T₄, respectively. With a uniform k = 500 ms⁻¹, these are 0.5, 0.9, 0.9, and 0.9 µs. The lifetimes of the TSE microstates in landscapes having 1 or 2 wells stabilized are shown in Fig. S2.

The coupling values (Fig. 3 A) show that (a) except for αM2M3–TBS pairs, most significant interactions are apparent with αC separations of less than ∼15 Å, and (b) many residues that are this close (or closer) to each other do not interact significantly. Evidently, for many amino acids, the C versus O free energy change is not influenced by those of other residues. In the AChR allosteric transition, side chains appear simply to “resettle” to influence the gating equilibrium constant and P_O by whatever changes happen to their local (within ∼15 Å) environment.

αM2M3–TBS coupling and the modal activity of αP265A indicate that these separated regions interact long range. The results show only a fuzzy “cloud” of communication with the TBS rather than a specific set of interacting amino acids. We suspect that the strengths of the αM2M3–TBS interactions are sufficient to account for their common phi values when agonists are present because they are on par with short-range coupling magnitudes. We hypothesize that long-distance interactions synchronize the gating rearrangements and energy changes at αM2M3 and the diliganded TBS.

With agonists, and for the TBS mutations we investigated, the main energy perturbation was with regard to the coupling constant (ΔΔG_cc; Eq. 2). Although we did not measure M2M3–TBS communication in unliganded gating, evidence suggests that it is minimal. Mutations that induce a high constitutive P_O also enforce a high affinity (Purohit and Auerbach, 2013), so the unliganded TBS appears to undergo a low→high affinity rearrangement without agonists, with an associated energy change that is part of ΔG₀. Without agonists, aromatic TBS residues have phi values that are similar to the rest of the ECD (Fig. 5 B). This indicates that their gating energy changes occur not with the affinity rearrangement (phi ∼1) but rather with ECD rotation/compaction (phi ∼0.8). Also, phi values of αM2M3 and TBS residues are similar in liganded but different in unliganded gating. We speculate that without agonists, M2M3 and the TBS do not communicate.

The coupling mechanism has been suggested to involve a mechanical link between the TBS and M2M3, such as the displacement of a β strand or a domino-like cascade of side chain/domain bumping. The structural elements of a rigid-body movement (or a cascade) should all reach the TS together and therefore have the same phi. However, with agonists, the TBS and αM2M3 are separated by a sea of lower-phi amino acids (Fig. 1 A), and with or without agonists, phi values of residues in αloop2 (base of β2), the salt bridge (base of β10), and αloop7 (base of β7) are characteristic of the rest of the ECD rather than of the TBS and αM2M3. The current phi map does not support a mechanical link, but phi values for many ECD residues have not yet been determined (including in β9 and β10), and a high-phi, rigid-body connection remains a possibility.

Proteins are not static objects—bonds fluctuate, side chains adopt rotamers, and domain linkages are flexible. One way to move energy (between M2M3 and the TBS) and leapfrog an intervening region (the ECD β-barrel) is by correlated changes in backbone dynamics (Taly et al., 2005; Bu and Callaway, 2011; Changeux, 2014; Garret et al., 2014). Long-range communication, including allostery, can be mediated by changes in protein “breathing” motions rather than by changes in mean atomic positions (Popovych et al., 2006; Tsai et al., 2008; Zheng and Thirumalai, 2009; Ma et al., 2011). The structural changes in gating at the AChR ECD–TMD domain interface are not known, but it is conceivable that the upward displacement of the M2M3 linker apparent in GluCl and GLIC could be sufficient to influence the dynamics of both M2M3 and TBS loops, to link these separated regions energetically without perturbing the intervening ECD. For example, a loop↔helix transition is correlated with change in structure and function at a distant site in kinesin (Kikkawa et al., 2001), myosin (Skolnick, 1987; Harrington et al., 1988), GPCRs (Kjeldgaard et al., 1996; Sprang, 1997), and AcrR (Manjasetty et al., 2016). Certainly, M2M3–TBS coupling in pLGICs could be by both a rigid-body link and backbone dynamics, as these mechanisms are not mutually exclusive.

The short-range interactions between αM2M3 and other loops at the domain interface are complex (Fig. 4 A). Side chains in the linker interact with those both in the same subunit (loop2 and loop7) and the adjacent subunit (εpreM1 and εloop9). Consistent with the pin-into-socket hypothesis, the N-terminal αM2 positions αI260 and αV261 interact with αQ48 in loop2. Interestingly, amino acids at the top of M2 (αE262 and αL263) also interact long range with the TBS. Within αM2M3, the only significant interaction was between αI260 and αP265, both of which are phi ∼1 amino acids.

Overall, it appears that αP265 plays a central role in orchestrating gating rearrangements at the domain interface. This proline is highly conserved (Jaiteh et al., 2016), has a phi ∼1 with and without agonists, induces affinity modes, is coupled to most TBS amino acids, and interacts with the ECD and TMD of the same and the adjacent subunit. These results are consistent with the previous suggestion that a change in the M2M3 proline omega dihedral angle initiates the global opening transition, although perhaps by a partial rather than a complete cis-trans isomerization. However, this proline was our focus from the outset, and more coupling measurements elsewhere at the interface are needed to confirm whether or not it is uniquely important.

The strong coupling between αQ48A and the top of M2 supports the notion that gating movements of loop2 and M2 are linked. The αQ48A mutation has a “catalytic” effect on gating (Chakrapani et al., 2004), which suggests that the loop2–M2 interaction mainly influences the energy barrier for ECD rotation/compaction rather than the stability of the active-state configuration. It is possible that in the opening process, an extension of the M2 helix (Fig. 1 C) facilitates ECD twisting by a direct interaction with loop2.

Although we did not find many αM2M3 residues coupled with αloop7, we note that only a few pairs were examined and not with the absolutely conserved proline here because all mutations of αP136 in AChRs eliminate expression or function.

γP121 on the complementary side of the TBS interacts significantly and long range with two high-phi αM2M3 residues, αS268 and αP265 (Fig. 2). At the fetal, liganded TBS, this proline also interacts locally and strongly with residues in loops B and C (Fig. 4 B). Hence, part of the energy connection between αP265 and the TBS may involve γP121. It has not escaped our notice that many of the long-distance interactions involve prolines, for instance γP112 (loopE), αP197 (loopC), γP121 (β5′), and αP265 (M2M3) and that αP221 (M1) forms the discrete ECD–TMD phi border (Purohit et al., 2015). It is possible that in AChRs these (and other) prolines are nodes in a network for long-distance communication via the backbone.

Neither phi nor flip provide complete descriptions of the intermediate states that separate C from O, but together, they begin to illuminate the gating TSE. Phi suggests the number of intermediates, the sequence and locations of the rearrangements, and the relative energy barriers associated with each, but it does not provide temporal information. Flip/primed shows intermediate states directly along with lifetime information, but it is at the extreme edge of the patch clamp’s capabilities and therefore requires two assumptions, that the brief shut intervals all have exactly zero current (but see Cymes and Grosman [2012]) and arise from a single exponential component (a single intermediate state) even if all of the detected events are longer than that time constant. The results suggest that the modal phi distribution and the flip/primed component are complementary and lead us to propose that the TSE of AChR gating is a single reaction trajectory with four intermediate states that each have a lifetime of ∼1 µs.

We are at present unable to derive a unique energy landscape for the AChR gating TSE. Flip/primed could reflect sojourns in the ensemble of all four TSE intermediates, a combination of a few, or perhaps just one, as proposed. It is easy to imagine a far more complex TSE landscape. For instance, the ECD may not rotate and compact synchronously, or the TMD helices might rearrange in a specific sequence. Revealing these microscopic features is now beyond the resolution of phi and flip analyses.

The results indicate that (a) many mutational energy changes have only short-range effects (Fig. 3 A), (b) phi values are modal (Fig. 6 A), and (c) with each phi population mapping to a discrete structural domain (Fig. 1 A). Also, a reaction chain with four brief intermediate states can account for both (d) the phi distribution (supplemental text) and (e) the flip/primed shut component (Fig. 6). Together, these observations suggest the following model for AChR gating (Fig. 7 and Video 1). (Again, C and O represent stable structures and gating refers to the full C↔O isomerization.) We use AChR α subunit diliganded phi values to locate, sequence, and define an energy landscape for the gating rearrangements, flip/primed to provide temporal information about the intermediates, MCA values to describe interactions as being local versus nonlocal, and other pLGICs to associate the on–off configurations with intermediate structures.

The model has five moving protein parts that switch between a locally off or locally on configuration: M2M3, TBS, ECD, TMD, and gate. In resting-C, all component structures are off, and in active-O, all are on. Between these two structures, there is a mixture. For clarity, each reversible on↔off local structural transition is named as follows: “click” (upward displacement of M2M3), “hold” (rearrangements at the TBS that underpin the affinity change), “twist” (rotation and compaction of the ECD), “tilt” (rearrangements of the TMD helices including straightening of M1 and radial bending of M2), and “dilate” (unpacking of the hydrophobic gate to form a bubble). The lipid, too, is perturbed in the dilate step (Purohit et al., 2013). The nonprotein components of the model are the agonist, ions, and water. “Wetting” of the gate region (collapse of the gate bubble) initiates fusion of the extra- and intracellular water compartments and ion conduction. We assume that M2M3–TBS coupling synchronizes click-and-hold so there are only four short-lived TSE intermediates (preopen states), all assumed to be nonconducting, which in some combination give rise to flip/primed. The energy barriers separating the TSE intermediates are small, so there is a thermal, random walk across the barrier region. As parts of a Brownian machine, the protein components are a linear actuator (M2M3), a latch (TBS) with a pin (agonist), a screw cap (ECD), a valve stem (TMD), and a plug (gate/bubble).

In the first step of the diliganded opening process, M2M3 and the TBS move (phi ∼0.95). These regions are coupled, so the agonist sites switch to the high-affinity conformation (hold) at approximately the same point in the reaction as the upward displacement of M2M3 (click). The second step is ECD twisting (rotation/compaction; phi ∼0.8). The relative heights of the exit barriers from T₁ indicate that this twisting is more difficult than reversing click-and-hold, so the receptor returns rapidly (<1 µs) to C a mean of ∼3.2 times before the ECD twist is achieved (Fig. 6 B). It is possible that twist and hold stabilize each other by virtue of favorable and local interactions at the TBS. Continuing, the ECD twist, when it eventually does happen, perturbs the domain interface to favor TMD tilting (phi ∼0.6) by short-range interactions that involve ECD β strands, interfacial loops, and TMD helices. Tilting causes gate dilation (phi ∼0.3), perhaps by a direct mechanical linkage with the αM2–13′ amino acid that, based on its relatively high phi ∼0.5, has been suggested to hold together the hydrophobic cluster (Purohit et al., 2013). We speculate that unpacking the gate creates a bubble (Sauguet et al., 2016) that in the final step of opening (phi ∼0) collapses and fills with water to start ion flow (Roth et al., 2008; Zhu and Hummer, 2010). The closing process is the reverse of these steps and begins with a spontaneous separation of the water into intra- and extracellular compartments, without the need for agonist dissociation (Calimet et al., 2013). In the model, gating appears as a staggering and reversible conformational cascade (Video 1).

The simulations suggest that the approximate microsecond shut intervals observed in high-resolution electrophysiology recordings correspond to sojourns in the TSE, entered from O. In the model, these events mainly reflect bubble reformation. We cannot distinguish whether flip is a pause in the reaction chain in just one TSE micro-well (Scheme 2) or whether this component reflects sojourns in an ensemble of intermediate states (Scheme 3). A time resolution of <1 µs might distinguish these possibilities and determine whether the brief shut component is a single exponential or the tail of several closely spaced exponentials.

With two bound ACh molecules, the apparent opening and closing rate constants (C↔O, full transits across the TSE) in WT mouse adult AChRs are ∼50 and 2.5 ms⁻¹, to yield shut- and open-channel lifetimes of 20 and 400 µs (−100 mV, 23°C, ∼25-µs time resolution). Using Scheme 3 with 300 ms⁻¹ for the unconstrained rate constants to match phi and flip, and assuming 1,000 ms⁻¹ for the C→T₁ rate constant (close to the maximum; Chakrapani and Auerbach, 2005), the extended kinetic model with all states fully liganded (gating with 2 ACh) is shown in Scheme 4. The rate constants are in ms⁻¹. In Scheme 4, the lifetime of diliganded C is only 1 µs, whereas the experimental shut duration associated with the C→O transition is 20 µs. This discrepancy is because, on average, approximately three exits from C and sojourns in the TSE (lifetime ∼6 µs; Fig. 6 B) are required to achieve O (Eq. S7 in the supplemental text). The experimental lifetime of diliganded O is ∼150 µs, which in the model is the time required for the gate bubble to reform. This also happens, on average, approximately three times per full passage back to C, to give an apparent open lifetime of ∼500 µs. There is no single rate-limiting step in opening with ACh, but ECD twisting appears to be the most difficult. In contrast, in the model, bubble formation is rate limiting for channel closing. Bubble reformation requires ∼150 µs, which is at least 100 times longer than the time required for bubble collapse. Apparently, the bubble is unstable compared with contiguous water.

The apparent opening and closing rate constants for unliganded gating of WT AChRs are ∼10⁻³s⁻¹/10⁴ s⁻¹ (Purohit and Auerbach, 2009; Jha and Auerbach, 2010). There is no phi map or flip measurement for unliganded gating, so we cannot calculate this TSE. However, we know that in a REFER the unliganded data point falls on the same line as that for agonists (Fig. 2 in Auerbach [2010]) and that many residues in the TMD have similar phi values with and without agonists. These results indicate that the unliganded TSE profile is similar to that for diliganded gating, except for “early” barriers close to C.

We speculate that unliganded WT AChRs have a small opening rate constant because with only water at the binding sites, the energy barriers for TBS hold and ECD twist both are large. Scheme 5 divides these energy penalties equally and reproduces approximately the apparent unliganded gating rate constants (gating without agonists). The increase in the T₁↔T₂ barrier (Fig. 6 E) without a concomitant one for C↔T₁ predicts a different phi pattern without ligands. However, the calculated phi values of Scheme 5 are 0.99, 0.83, 0.61, 0.35, and 0.0 (supplemental text), which are indistinguishable from the measured diliganded ones (Fig. 6 A). We hypothesize that the unliganded opening rate constant is tiny because both hold and twist are difficult, so that ∼18,000 exit attempts from C are required before O is finally reached.

Both with and without agonist, the apparent gating rate constants measured in a standard patch clamp recording do not give accurate estimates of the exit rates from C and O. Either ∼3 (with ACh at both TBS), ∼20 (with choline), or ∼18,000 (with water) exit attempts from C are required to reach O, to make the apparent opening rate constant underestimate the exit rate from C by these factors. The model predicts that during silent periods in electrophysiology experiments, AChRs are dynamic structurally and undergo partial gating rearrangements without generating a current signal (Video 1).

The kinetic and structural models for diliganded AChR gating are consistent with local energy changes, phi values, locally closed GLIC structures, flip/primed, a central role for the M2M3 proline, and suggestions that a gate bubble is a metastable barrier to ion flow. They are also consistent with mechanical interactions at the domain interface that have been suggested to be involved in opening, including those involving the β2 and β10 strands. However, in our model, these strains are mainly local and are initiated by ECD twisting rather than being linked directly to actions at the TBS.

We thank Marlene Shero, Mary Merrit, and Mary Teeling for technical assistance, Tapan Kumar Nayak for helpful discussions, and Aashutosh Vihani for data shown in Fig. 3 A.

This work is supported by National Institutes of Health grants NS-023513 and NS-064969.

The authors declare no competing financial interests.

Author contributions: S. Gupta, S. Chakraborty, and R. Vij performed research and analyzed data; S. Chakraborty derived the equations in the supplemental text; A. Auerbach designed experiments, simulated currents, and wrote the paper.

Sharona E. Gordon served as editor.