Mechanism of pore opening in the calcium-activated chloride channel TMEM16A

The anion channel TMEM16A is activated by intracellular Ca2+ in a highly cooperative process. By combining electrophysiology and autocorrelation analysis, we investigated the mechanism of channel activation and the concurrent rearrangement of the gate in the narrow part of the pore. Features in the fluctuation characteristics of steady-state current indicate the sampling of intermediate conformations that are successively occupied during gating. The initial step is related to conformational changes induced by Ca2+ binding, which is ensued by rearrangements that open the pore. Mutations in the gate shift the equilibrium of transitions in a manner consistent with a progressive destabilization of this region during pore opening. We come up with a mechanism of channel activation where the binding of Ca2+ induces conformational changes in the protein that, in a sequential manner, propagate from the binding site and couple to the gate in the narrow pore to allow ion permeation.

L igand-dependent gating is a tightly regulated process in ion channels whose open probability increases in response to the binding of an agonist to a specific site of the protein.
Depending on the class of channels, the chemical nature of such ligands may range from ions to small molecules or even proteins. Ligand binding generally triggers a conformational change in the protein around the binding site, which is propagated to the pore region to open a gate that impedes ion conduction in the closed state of the channel. Gating is frequently found to be a cooperative process that proceeds via a defined set of intermediate states from the closed to the open conformation of the protein.
Whereas these intermediates are transient and thus often escape structural characterization, they were in many cases successfully characterized by kinetic analysis of single-channel recordings [1][2][3] .
TMEM16A is a ligand-dependent anion-selective channel that is activated by an increase in the intracellular Ca 2+ concentration [4][5][6] . The protein harbors two pores that are each contained within a single subunit of the homodimeric protein [7][8][9] . Both pores act independently and are activated by the binding of two Ca 2+ ions to a conserved site that is located within each subunit in proximity to the ion permeation path 10,11 . The location of the site within the membrane confers voltage dependence to the binding step and the proximity of bound Ca 2+ to the pore shapes the electrostatic potential for conduction 12 . Besides the change of pore electrostatics, the binding of Ca 2+ also triggers a conformational change in a membrane-spanning helix of the protein (α6) that coordinates the bound divalent cations after its rearrangement 8 . The movement of α6 is in turn coupled to the release of a steric gate in the narrow neck of an hourglass-shaped pore. An accompanying study 13 has identified this gate to be composed of three hydrophobic residues that are located at the intracellular end of the narrow neck. Whereas major factors controlling ion flow have been identified in previous studies and a general model for activation was proposed 14 , the detailed sequence of events and the existence of intermediates in the gating process, which together define the activation mechanism, are still elusive. Although this process would ideally be characterized by single-channel analysis, such studies are prohibited by the low conductance of TMEM16A 10,15 .
Here, we combine macroscopic electrophysiological measurements and autocorrelation analysis to investigate the mechanism of channel activation. We show that the fluctuation characteristics of steady-state current in TMEM16A is consistent with the sampling of intermediate conformations that are successively occupied during gating. Our results suggest a mechanism of channel activation where the binding of Ca 2+ induces conformational changes in the protein that, in a sequential manner, propagate from the binding site and couple to the gate in the narrow pore to allow ion permeation.

Results
Autocorrelation analysis. Gating in ion channels is defined as the transition between non-conducting and conducting states, which is usually associated with conformational changes of the protein.
Whereas gating events can be readily observed in single-channel recordings for channels of large conductance, they are hidden in the noise for those with low conduction rate such as TMEM16A 10,15 . Nonetheless, stochastic transitions in an ensemble of independent single channels give rise to fluctuations around the mean of the macroscopic steady-state current 16,17 ( Supplementary Fig. 1a). The statistical properties of these fluctuations may be quantified from its power spectrum [18][19][20][21] , which is the Fourier transform of the autocorrelation function of the current 22 . For single-channel fluctuations, the autocorrelation function is characterized by an exponential decay with time constants corresponding to the system's relaxation times 19,20 . Time intervals shorter than the relaxation times are expected to yield higher correlation, as it is more likely that the channel remains in or resamples the open state, while the correlation becomes zero at much longer times as any co-occurrence of opening events separated by these time intervals is essentially random.
We illustrate these properties by calculating the autocorrelation function and power spectrum from simulated single-channel trajectories using various gating models ( Supplementary Fig. 1). In all cases, the autocorrelation function consists of a sum of exponentials and the power spectrum, being the Fourier transform of the former, is a sum of Lorentzian components ( Supplementary  Fig. 1b, c). The number of such components defines the minimum number of transitions between states in the underlying mechanism and their respective amplitudes and corner frequencies are determined by the corresponding rate constants. Because the power spectrum of a macroscopic current is a multiple of the spectrum of a single-channel record for independent and identical channels, analysis of the former is equivalent to analyzing its single-channel counterpart in Fourier space, from which model parameters may be estimated (Supplementary Figs. 2 and 3 and "Methods"). The effect of non-stationarity ( Supplementary Fig. 4) and an examination of the parameter estimation process are discussed in the methods. In this study, we applied this analysis to macroscopic currents of TMEM16A to investigate ligand binding and the subsequent transition into transient and stable conformational states. We also explored the relationship of the observed states to structural determinants during channel activation.
Mechanism of Ca 2+ activation. Recent cryo-EM structures of TMEM16A have revealed a conformational change involving the rearrangement of a pore-lining helix (α6) upon Ca 2+ binding in a process that precedes pore opening 8 . To gain mechanistic insight into this process and into the sequence of events following Ca 2+ binding, we used the described kinetic analysis to identify additional conformational transitions that escaped structural characterization. For that purpose, we obtained a family of power spectra by recording steady-state currents over a range of Ca 2+ concentrations, as shown in Fig. 1a, b. Evident from the trajectories in the time domain is that slow fluctuations (manifested in Lorentzian components with lower corner frequencies) vanish and fast fluctuations (those with higher corner frequencies) become more prominent as the Ca 2+ concentration increases (Fig. 1a, b). The ligand dependence of the slow fluctuations and their inverse correlation with saturation hint at events related to ligand association and dissociation. A shift in the spectral frequencies as Ca 2+ concentration is elevated is consistent with the hastening of the activation response time in concentration 8,23,24 and voltage jump experiments 25 , which are governed by the same set of time constants. At saturating Ca 2+ concentrations, the spectrum is reduced to reflect fluctuations that are caused solely by gating transitions as Ca 2+ remains bound with very high probability (Fig. 1a-c). The presence of three Lorentzian components under saturating conditions (as judged by Akaike's information criterion (AIC) of models with a different number of components 26 , Supplementary  Fig. 5) indicates the sampling of at least four conformational states when the channel is fully occupied by Ca 2+ . In conjunction with known structures, where the Ca 2+ -bound conformation might resemble a ligand-activated pre-open state, our data thus point towards the presence of intermediate gating steps that relay α6 activation to the opening of the pore.
The above results motivated us to construct a kinetic model where activation of a single TMEM16A pore proceeds in three steps, traversing different conformational intermediates. In these states α6 is either mobile and disengaged as found in the Ca 2+ -free structure of the channel ('α6-loose'), or rigid and in contact with the binding site as displayed in the Ca 2+ -bound structure ('α6tight'). From the 'α6-tight' conformation, the protein transits into another partly activated but still non-conductive state ('pre-open') before reaching the conductive open state ( Supplementary Fig. 2). Three different Ca 2+ occupancies (0, 1, and 2 Ca 2+ ) give rise to 12 hypothetical states including three discrete open states reflecting the influence of bound Ca 2+ on conduction, which was previously observed for mutants showing pronounced basal activity 8,12 . All states are related by an allosteric Monod-Wyman-Changeux (MWC) mechanism 27,28 . In light of the low basal activity of wild type (WT) and the high cooperativity of activation, a reduced version of the model where only seven states are explicitly included was used in our subsequent kinetic analysis ( Fig. 1d and Supplementary Fig. 2). Of the resulting 14 kinetic parameters, 10 (k 01 , k 12 , k 21 , k 23 , k 32 , k 42 , k 35 , k 53 , k 54 , k 56 ) were determined independently by fitting the data (Fig. 1e). The rate k 65 describing the binding step of the first Ca 2+ was assumed to be diffusionlimited. The values of k 24 and k 45 are defined by microscopic reversibility and the empirically estimated limiting Ca 2+ binding affinity obtained from an accompanying manuscript 13 . k 10 , the final opening rate constant was determined using knowledge of the maximum open probability of the channel (P o max ) obtained from non-stationary noise analysis ( Supplementary Fig. 6). The resulting equilibrium constants relating to the connected states (L 10 , L 21 , L 32 , L 42 , L 53 , L 54 , L 65 ) were obtained from the ratio of corresponding rates.
By simultaneously considering the concentration-response relation, we performed a global fit of the entire set of power spectra to the described model where Ca 2+ binding transitions are explicitly included (Fig. 1b, c). The wealth of data obtained from six different Ca 2+ concentrations ensured a unique fit of the kinetic parameters ( Fig. 1e and Supplementary Table 1). Consistent with classical ligand activation models, the mechanism involves a cycle that couples affinity increment and transition efficacy (Fig. 1d) and can be viewed as the mechanistic counterpart of the conformational transition between apo and Ca 2+ -bound TMEM16A observed in the structures, which is also governed by ligand binding. The affinity of binding of the second Ca 2+ (L 53 ) increases from~1 µM in the resting state to~40 nM (L 42 ) in the next closed state (Fig. 1d), consistent with the origin of this conformational transition being α6 activation. The rate of binding of the second Ca 2+ to the singly bound resting state (C 5 , rate constant k 53 ) is approximately three orders of magnitude slower than the binding of the first Ca 2+ to the apo state (C 6 , described by the diffusion-limited rate constant k 65 ), but becomes two orders of magnitude faster in the next closed state (C 4 , k 42 ) (Fig. 1e and Supplementary Table 1), which, again, is consistent with the described mechanism of α6 activation where additional hydrophilic and negatively charged residues on α6 get into direct contact with the bound Ca 2+ (ref. 8 ). Such electrostatically assisted association is also evident in Poisson-Boltzmann calculations, where the electrostatic potential of the lower binding site becomes more negative when α6 assumes an activated conformation ( Fig. 2a, b). The close relationship between the electrostatic potential at the binding site and the estimated association rate constants of the respective configurations ( Fig. 2c) emphasizes the correspondence of the proposed mechanism to structural states. A thermodynamic consequence of this process is the 34-fold increase in the efficacy of α6 activation when two Ca 2+ are bound (L 32 ) compared to that of the mono-liganded state (L 54 ) ( Fig. 1d and  Based on the estimated parameters, we examined the most probable sequence of conformational changes upon Ca 2+ binding (Fig. 2d). As there is no detectable basal activity and the incorporation of a single apo closed state is sufficient to account for the data, the binding of the first Ca 2+ most likely occurs in the state where α6 is in its resting conformation (C 6 ). The route taken by the channel is then dependent on Ca 2+ concentration due to the Ca 2+ dependence of the L 53 transition. In the physiological range (at sub-micromolar concentrations), the transition of α6 into a tight conformation occurs mainly through the singly bound state (C 5 to C 4 ) until the Ca 2+ concentration is high enough to allow the second Ca 2+ binding step to 'compete' for the singly bound state (C 5 to C 3 ) (Fig. 2d). Once two Ca 2+ ions are bound and α6 is in its activated conformation (C 2 ), the channel can then undergo the 'gating steps', which involve a pre-open state (C 1 ) and a final transition that opens the pore (L 10 ).
Rearrangement of the gate region during activation. Equipped with the described framework for the kinetic analysis of macroscopic recordings, we sought to understand the mechanistic origin of the gating steps (i.e., the steps connecting C 3 , C 2 , C 1 , and O 0 ). To this end, we combined spectral and double-mutant cycle analysis [29][30][31][32] involving residues at the inner gate (consisting of Ile 550, Ile 551, and Ile 641) (Fig. 3, Supplementary Fig. 7 and Supplementary Table 2) that give rise to constitutive activity when mutated to alanine as described in an accompanying paper 13 . In non-stationary noise analysis, we observed that only the alanine mutation of Ile 641, which forms the core of the gate, gives rise to a substantial increase in P o max in the Ca 2+ -bound state, while I550A displays a modest increase and I551A a moderate decrease in P o max (Fig. 3a, c). In the case of I641A and I550A, this provides additional evidence for the stabilization of the open pore conformation by the respective mutations. Analysis of its power spectrum at a saturating Ca 2+ concentration suggests that the macroscopic phenotype of I641A originates from an increase in transition efficacy in the initial and the final steps (manifested in the increase in L 32 and L 10 ) that correspond to α6 rearrangement and pore opening ( Fig. 3d-f), suggesting a destabilization of the gate during these two transitions. In contrast, no substantial change in the respective equilibrium constants was observed for I550A and I551A (Fig. 3d-f), suggesting that these two residues may affect all of the connected states to a similar extent.
In the absence of Ca 2+ , the corresponding power spectra of gate mutants reveal comparable transitions as observed at a saturating Ca 2+ concentration ( Supplementary Fig. 8), which Although much higher than for WT, both the efficacy and the kinetics of the pre-opening (L 21 ) and the opening steps (L 10 ) are consistently lower in the apo than in the fully Ca 2+ -bound state, leading to a moderate reduction in the P o max (Supplementary  Tables 2 and 3). This observation further confirms the role of the bound Ca 2+ ions in influencing the energetics of the gating transitions even in mutants with considerable basal activity. This enhancement acts concomitantly with the release of an electrostatic gate that impedes anion conduction in the open state of the apo channel 12 .
To examine how the truncation of sidechains of Ile 550 and Ile 551 influences gating, we analyzed the respective mutations in double-mutant cycles 29,32 . In this analysis, interdependent contributions of a mutation originating from interactions with another residue can be quantified by their coupling energy (G coupling ). G coupling denotes the difference between the energetic effects of a mutation introduced on the wild-type protein and in the background of the other mutation. If the two perturbations are independent, its magnitude is zero as the respective backgrounds do not have an impact on the effect of the mutation. In contrast, G coupling would deviate from zero if both residues interact functionally. Conversely, energetic contributions of a mutation that are independent of a particular pairwise interaction (expressed as G indep ) may be inferred from the effect of the same mutation introduced on the background of a second mutation where the sidechain of the interacting residue is truncated (see "Methods").
An exchange of interactions at the α4-α6 interface stabilizes the open state. In the Ile 550/Ile 641 and Ile 551/Ile 641 mutant cycles, the reversal of the I641A effect on the opening step by mutations I550A and I551A suggests that, besides their role in the stabilization of the closed gate by interacting with Ile 641, the two residues on α4 may also be involved in a stabilization of the open state by interactions that are independent of Ile 641. Because Ile 551 and Gln 649 are closely apposed in the Ca 2+ -bound structure (Fig. 4a), we investigated their potential interaction in stabilizing the open pore in a double-mutant cycle. Since the mutation Q649A consists of perturbations that change both the polarity and volume of the sidechain, we paired I551A with either Q649A or Q649L to disentangle the two effects (Fig. 4).
When introduced on its own, the mutation of Gln 649 to Ala exerts little effect on the gating transitions whereas its mutation to Leu increases the efficacy of the first and last step (Fig. 4b-d). The introduction of I551A on the Q649A background reverses the dampening effect of I551A on the opening step (L 10 ), which now increases the efficacy of opening (Fig. 4d). This results in a considerable negative coupling energy (Fig. 4e), suggesting that Ile 551 and Gln 649 interact to stabilize the open pore, a notion that is also supported by elevated conduction barriers in these mutants (presumably by a more collapsed pore geometry) as reflected in the outward rectification of current (see below). Examination of the stepwise mutant cycles reveals that the coupling persists when the partial charges of Gln 649 were removed on the Q649L background (Fig. 4e), suggesting that the interactions are predominantly mediated by the volume instead of the polarity of the sidechain and the open pore may thus in part be stabilized by van der Waals forces. Representation is as in Fig. 3a. b-d Equilibrium constants for mutants at the interface for transitions b, L 32 , c, L 21, and d, L 10 . Bars indicate the best-fit values of the averaged data shown in Supplementary Fig. 7d (Q649A, n = 7; Q649L, n = 7; I551A/Q649A, n = 8; I551A/Q649L, n = 7). Errors are 95% confidence intervals. e Coupling energy (G coupling ) and f independent energetic contribution (G indep ) of the indicated residues. Bars indicate quantities calculated using e, Eqs. 31-32, 35 and f, Eq. 34 from the best-fit values shown in b-d. Errors correspond to standard errors. Asterisks indicate significant deviation from zero in a two-sided one-sample t-test (e I551/Q649: ***p = 2e−7 and ***p = 3e−5; I551/L649: ***p = 5e−6 and *p = 0.023; f ***p = 1e−7, ***p = 6e−4, and ***p = 0.001). Consistent with the previously described weak interaction between Ile 551 and Ile 641 that stabilizes the closed pore in the final gating step (L 10 , Fig. 3g), we show here that Ile 551 stabilizes the closed pore independently of Gln 649 (as reflected in the positive G indep for the same step for I551A, Fig. 4f). Together, these observations explain how I551A reverses the effect of I641A on pore opening (L 10 ) in the double mutant. In the described process, Gln 649 competes with Ile 641 for interaction with Ile 551, with the former interaction stabilizing the open (O 0 ) and the latter the pre-open state (C 1 ). This suggests a mechanism where Ile 551 dissociates from Ile 641 and in turn interacts with Gln 649 as the gate opens.
Arrangement of the gate in the open pore. Next, we investigated the interactions of the same residues when the pore is open. For quantification, we extracted the energetic effects of alanine mutants of the three gating residues and Q649A on ion conduction and casted these in double-mutant cycles (Fig. 5a). Energies were obtained from a fit of I-V relations to a threebarrier model that was introduced previously 7,12 . The fit yields two kinetic parameters (σ β and σ h ) that can each be converted into an energy difference of the inner and central barrier relative to the outer barrier (Fig. 5b), which allows the approximate localization of effects of mutations on the ion conduction path.
Our data indicate a slight enhancement of the relative rates of anion diffusion at the inner entrance of the neck (σ β ) and inside the narrow pore (σ h ) for I550A and I641A due to a decrease of both barriers (Fig. 5c, d). This behavior is also reflected in the increase of the single-channel conductance in both mutants as observed in non-stationary noise analysis ( Fig. 3c and Supplementary Fig. 6). Reversed effects are found in mutants Q649A and I551A, where the relative rate across the inner barrier  (Figs. 3c, 5c, d and Supplementary Fig. 6). For the I550A/I641A and I551A/I641A pairs, the coupling energy for the inner barrier is small and positive, while that for the middle barrier is not significantly different from zero (Fig. 5e). The low amplitudes of the coupling energies are consistent with a widening of the gate region in the open pore where obstruction to conduction is minimal. For the I551A/Q649A pair, however, the coupling energy is slightly negative for the inner and strongly positive for the central barrier (Fig. 5e), suggesting that both residues interact in the open state to keep the neck in a conductive conformation. Examination of the stepwise mutant cycle on the Q649L background, which removes the effects of partial charges and which consists of only contributions from steric interactions, reveals considerable negative coupling energy for the central barrier (Fig. 5e). This indicates that the positive polarity of the coupling in the overall cycle is likely due to a compensatory effect in the I551A/Q649L cycle. Therefore, while the coupling between Ile 551 and Gln 649 is small in stabilizing the open pore conformation at the inner entrance of the neck, it is required to maintain a widened conformation inside the neck. This observation is consistent with an interaction between these two residues that stabilizes the open state as inferred from kinetic analysis (Fig. 4).
Timing of motion of the gate region. Our analysis revealed the successive rearrangements of the gate region leading to pore opening during channel activation. To gain further insight into the location of energy maxima connecting the different states, we analyzed the gating kinetics of a series of mutants in the gate region collectively to obtain the relationship between the forward (k f : k 32 , k 21 , k 10 ) or backward rate constants (k b : k 23 , k 12 , k 01 ) and their corresponding forward equilibrium constants (K eq : L 32 , L 21 , L 10 ) (Fig. 6). The relative timing of a local conformational change in a global transition can be characterized by its phi value (0-1), which is reflected in the slope when k f is plotted against K eq on a logarithmic scale [33][34][35] . A phi value close to one corresponds to a case where a mutation affects the transition state and the final state to a similar extent. This scenario is consistent with the involvement of the site of mutation in a structural rearrangement that occurs early on the reaction coordinate. Conversely, a value close to zero would suggest a rearrangement that occurs relatively late during the global transition. During the gating of TMEM16A, the locations of respective transition states of the gate region differ for distinct transitions. Whereas both k f and K eq are affected to a similar extent in the series of mutants for the intermediate transition (L 21 ), as reflected in a phi value close to 1, k f is virtually insensitive to changes in K eq for transitions involving the destabilization of the gate region (L 32 and L 10 ) which translates to a phi value of zero (Fig. 6). In all cases, the spread of K eq values obtained for different mutants further emphasizes the energetic involvement of the gate region in each of the global transitions. These results suggest that during transitions in which the gate is destabilized (i.e., L 32 and L 10 ), the associated rearrangement occurs late and is likely amongst the final motions involved in the widening of the pore.

Discussion
In this study, we have used electrophysiology in combination with a detailed kinetic analysis to unravel the gating mechanism of TMEM16A. By employing stationary noise analysis of macroscopic recordings at increasing Ca 2+ concentrations, we have identified distinct states of the channel that are successively occupied during activation. Transitions that dominate at saturating Ca 2+ concentrations where binding events are scarce are consistent with the sampling of at least four distinct protein conformations (Fig. 7a). The strongly Ca 2+ -dependent step corresponds to an early transition where the helix α6 is either loose as defined in the structure of TMEM16A obtained under Ca 2+ -free conditions or tight as observed in the Ca 2+ -bound form of the protein. The conformational rearrangement of α6 is coupled to the gate at the intracellular entry to the narrow neck, and leads to its successive destabilization via another pre-open intermediate to finally reach a conductive conformation.
Our kinetic analysis suggests that successive Ca 2+ binding to TMEM16A proceeds via distinct trajectories depending on its concentration (Fig. 7a). Whereas the first Ca 2+ binds to a state where α6 is in its loose conformation, the binding of the second Ca 2+ would proceed in the same α6-loose conformation only at unphysiologically high Ca 2+ concentrations, while it would occur more frequently in the state where α6 has assumed an activated conformation at physiological Ca 2+ concentrations. The mutual stabilization of the two events increases the affinity for the binding of the second Ca 2+ , underlying the positive cooperativity of channel activation.
Following the binding of two Ca 2+ ions and α6 activation, the protein further progresses towards pore opening. The presence of pre-open intermediates during this transition suggests that the channel undergoes successive conformational rearrangements before entering a conducting state. This could be reflected in the Ca 2+ -bound structure of TMEM16A, where such a pre-open intermediate might have been stabilized in a detergent environment and the absence of bound PI(4,5)P 2 , the latter of which prevents current rundown in excised patch recordings 36,37 . In this structure, α6 has already undergone its ligand-induced rearrangement whereas the narrowest part of the pore might not have expanded sufficiently to accommodate permeant anions. Although still not conductive, the partial widening of the inner pore in this structure signifies a destabilization of the gate region as a likely consequence of the tightened conformation of α6. Such coupled movement is also manifested in our kinetic analysis where Ile 641 was found to stabilize the relaxed conformation of Fig. 6 Rate-equilibrium free-energy relations of the gating steps. Plots show the experimental relation between the forward (blue) and backward (red) rate constants and the corresponding equilibrium constants for the transitions relating the gating intermediates and end points (from C 3 to O 0 ). Each data point corresponds to the parameter(s) obtained for a mutant in the gate region shown in Supplementary Fig. 7d. Errors correspond to 95% confidence intervals (respective error bars are not visible when they are smaller than the symbols). Solid lines are fits to the pair of rate-equilibrium free energy relations for each transition (Eq. 36). The estimated value of phi is presented as best-fit ± standard error.
α6 (α6-loose) by coupling with Ile 550 and Ile 551 in the gate region (Fig. 3). Consistent with this assumption, the disruption of any of the three gate residues in the apo state is likely to have a reciprocal energetic consequence on the α6 conformation, as observed in our kinetic model and also in the cryo-EM structure of the activating mutant I551A described in our accompanying study 13 .
The spectra of mutants with basal activity obtained in the absence of Ca 2+ reveal constitutive sampling of equivalent transitions as found in the Ca 2+ -bound state, underlining that the general mechanism of pore opening does not strictly require the presence of the ligand as expected for a mechanism that is based on an MWC model 27,28 . Despite the similarity of the accessible states, local structural differences of the open states in the presence and absence of Ca 2+ might be expected given the strong electrostatic influence of the bound Ca 2+ ions on the protein conformation. Such differences have been observed in the apo structure of I551A described in an accompanying study 13 , where α6 adopts a partially activated conformation in which the transition from an αto a π-helix conformation cannot be completed due to the energetic penalty of the transition and the electrostatic repulsion from the vacant binding site. Thus, although activation and pore opening can proceed without a transition into the strained π-helix conformation, the decrease in the efficacy of both the pre-opening (L 21 ) and opening steps (L 10 ) in the apo state suggests that this might have an energetic consequence on gating.
An opening equilibrium constant in WT that is several fold larger than one suggests that the partially destabilized gate found in the α6-tight and pre-open intermediates is intrinsically less stable compared to the fully open conformation (Fig. 1d). This might be attributable to the combined effect of weakened interresidue interactions in the destabilized state and the release of ordered water around the hydrophobic cage in the dewetted region upon opening. On the basis of coupling energies obtained from mutant cycles (Figs. 3-6), we come up with a structural interpretation of the events leading to the opening of the gate (Fig. 7b). Immediately preceding pore opening, retained interactions between Ile 550 and 551 on α4 with Ile 641 on the opposing α6 still stabilize the partially widened gate region, thereby maintaining a pore diameter that is sufficiently small to favor spontaneous dewetting, a condition prohibitive for ion conduction. Pore opening results as these interactions become disrupted and is accompanied by the formation of an alternative arrangement between α4 and α6 where Ile 551 engages in an interaction with the neighboring Gln 649 (on α6) to stabilize the open pore. This sequence of events is compatible with the at times counterintuitive characteristics of Ile 551, which, by toggling respectively between interaction with Ile 641 or Gln 649 on α6, may act as a switch that stabilizes either the pre-open or the open pore conformations.
Further insight into the relative location of transition states connecting the gating intermediates and end states was obtained from phi value analysis, which relates changes in the rate and corresponding equilibrium constants. Whereas a phi value close to one in the passage from the 'α6-tight' to the 'pre-open' conformation suggests a transition state that is reached early on the reaction coordinate, phi values close to zero for the first and last gating steps are consistent with rate-limiting steps towards the end of the process (Fig. 7a). For the last step, this suggests that the rearrangement of the gate is amongst the final events in a delocalized motion that opens the pore (Fig. 7b). Thus, although the gate region remains in a closed-like conformation until the end of the final gating step, other parts surrounding the pore might have already approached an open-like state earlier in the gating process. In this respect, a more extended rearrangement elsewhere in the pore region is likely, as evident from systematic mutagenesis that we report in an accompanying manuscript 13 where residues located towards the extracellular end of the neck also impact gating and contribute to a widening of the central part of the pore in the open conformation.
In conclusion, by establishing a general method to measure steady-state kinetics for channels whose low unitary conductance prohibits single-channel recording, we have provided a first detailed mechanistic view of the gating process in TMEM16A. While our current study has focused on the gate region, the expansion of the mutational analysis towards other residues of the pore might further define the molecular motions during gating. Given the conserved architecture of the family, our findings will be relevant for studying activation in other TMEM16 ion channels and lipid scramblases, and provide a functional template Fig. 7 Activation mechanism. a Cartoon depiction of the activation mechanism. The states are as in Fig. 1d. Inset (top) shows a schematic representation of the energy maxima in the gate region between connected states. The relative location of the transition states was obtained from phi value analysis. The two alternative pathways for activation at physiological and at high Ca 2+ concentrations are indicated. b Schematic depiction of the relationship between residues of the gate region in different steps during activation. The view is from the extracellular side. Gate residues located on α4 (Ile 550 and Ile 551) retain their interaction with the gate residue Ile 641 on α6 in the transition between the three closed conformations although the interaction between Ile 551 and Ile 641 successively weakens until it is replaced by an interaction with Gln 649 upon transition into the open state where the gate has dissociated and the pore has dilated to permit ion conduction.
for the therapeutic targeting of TMEM16A in pathological conditions such as cystic fibrosis 38 .

Methods
Molecular biology and cell culture. HEK293T cells (ATCC CRL-1573) were maintained in Dulbecco's modified Eagle's medium (DMEM; Sigma-Aldrich) supplemented with 10 U ml -1 penicillin, 0.1 mg ml -1 streptomycin (Sigma-Aldrich), 2 mM L-glutamine (Sigma-Aldrich), and 10% FBS (Sigma-Aldrich) in a humidified atmosphere containing 5% CO 2 at 37°C. HEK293T cells were transfected with 3 μg DNA per 6 cm Petri dish using the calcium phosphate coprecipitation method and were used within 24-96 h after transfection. Mutations were introduced with a modified QuikChange method 39 using the a,c variant of mouse TMEM16A (UniProt identifier: Q8BHY3-1) as the template and were verified by sequencing. Primers are listed in Supplementary Table 4.
Electrophysiology. Recordings were performed on inside-out patches excised from HEK293T cells expressing the construct of interest. Transfected cells were identified via the fluorescence of the Venus tag. Patch pipettes were pulled from borosilicate glass capillaries (O.D. 1.5 mm, I.D. 0.86 mm; Sutter Instrument) and were fire-polished with a microforge (Narishige) before use. Pipette resistance was typically 3-8 MΩ when filled with the recording solutions detailed below. Seal resistance was typically 4 GΩ or higher. Voltage-clamp recordings were made using Axopatch 200B, Digidata 1550, and Clampex 10.6 (Molecular Devices). Analog signals were filtered with the in-built 4-pole Bessel filter at 10 kHz and were digitized at 20 kHz. Solution exchange was achieved using a gravity-fed system through a theta glass pipette mounted on an ultra-fast piezo-driven stepper (Siskiyou). Liquid junction potential was found to be consistently negligible given the ionic composition of the solutions and was therefore not corrected. All recordings were performed at 20°C.
A symmetrical ionic condition was used throughout. Stock solution with Ca 2 + -EGTA contained 150 mM NaCl, 5.99 mM Ca(OH) 2 , 5 mM EGTA, and 10 mM HEPES at pH 7.40. Stock solution with EGTA contained 150 mM NaCl, 5 mM EGTA, and 10 mM HEPES at pH 7.40. Free Ca 2+ concentrations were adjusted by mixing the stock solutions at the required ratios calculated using the WEBMAXC program (http://web.stanford.edu/~cpatton/webmaxcS.htm). Patch pipettes were filled with the stock solution with Ca 2+ -EGTA, which has a free Ca 2+ concentration of 1 mM.
Analysis of current-voltage (I-V) relations. I-V data were fitted to a minimal permeation model that accounts for the most fundamental biophysical behavior of mTMEM16A as described previously 7,12 , where I is the current, n is the number of barriers, c i and c o are the intracellular and extracellular concentrations of the charge carrier, z is the valence of Cl − , V is the membrane voltage, and R, T, and F have their usual thermodynamic meanings. A = β 0 v is a proportionality factor where β 0 is the value of β when V = 0 and v is a proportionality coefficient that has a dimension of volume. σ h and σ β are respectively the rate of barrier crossing at the middle and the innermost barriers relative to that at the outermost barrier (β). The best-fit values of σ β and σ h at zero and saturating Ca 2+ concentrations were used to calculate ΔE α(σβ) and ΔE a(σh) , the difference between the activation energy at the innermost barrier and the middle barrier relative to that of the outermost respectively, using Non-stationary noise analysis. Non-stationary noise analysis was performed as described previously 10,40 . The current and variance were sampled by repeatedly activating and deactivating the channel using regularly spaced concentration jumps. The variance of such 50-100 aligned successive and kinetically identical currents at each time point during deactivation was calculated by computing the mean of the squared successive difference 41 , which mitigates the effect of nonstationarity at each isochrone and therefore allows the estimation of the variance in the presence of current rundown. Assuming the presence of a dominating conducting level, the data were fitted to where σ 2 N is the variance for N channels, i is the unitary current, I is the mean current, and the subscript bg denotes background. The introduction of I bg and σ 2 bg allows the accommodation of x-and y-translations respectively. To merge data from different patches, individual σ 2 bg -and I bg -corrected σ 2 N -I plots were brought to the same scale by normalization in both x and y directions according to a patch-specific parameter iN, the maximum achievable I for each patch if the Po was 1. Because each data pair ( I j ; σ 2 N j ) is unique, they were sorted according to the I values and were averaged using a Gaussian moving average filter, which allowed the central tendency of the normalized parabolas from different patches to be estimated. The averaged data were re-fitted to Eq. 3 without the σ 2 bg and I bg terms. This procedure allows the estimated Po to be directly read from the merged σ 2 N -I plots.
Computing the power spectrum. The power spectrum of a steady-state current recorded at +80 mV was obtained via Fast Fourier transform (FFT). For analysis, currents were recorded over a continuous period of typically 50 s, or 100 s for cases where sampling of lower frequencies was required. Recordings were filtered at 10 kHz and sampled at 20 kHz as described above. To reduce spectral leakage, a Hamming window was applied to mitigate edge discontinuities before Fourier transform 42 . The magnitude of the spectral components (S) is given by the sum of squares of the amplitudes derived from the real (a real ) and imaginary parts (a imag ), which was scaled according to Clampfit 10.6 (Molecular devices) using to yield P, the one-sided spectral density in A 2 Hz , where N is the transform length, f is the sampling frequency in Hz, and ω is a scale factor for the window function (f w ). The background spectrum of each patch, recorded at 0 mV where the current reverses, was subtracted from the raw spectrum to obtain the power spectrum specific to the channel. This resulting spectrum typically consists of three major components, which are the 1/f-like component at low frequency 18,43 , the Lorentzian components corresponding to molecular fluctuations, and an effectively constant term likely consisting of very fast components whose corner frequencies are not resolved within the bandwidth of the spectrum. In order to subtract the 1/flike and constant components, the total spectrum was fitted to to obtain an empirical description, where f is frequency, n is an exponent describing the decay, a 0 and a i are respectively the amplitude of the 1/f-like and the Lorentzian components, f ci is the corner frequency, and c is a constant. The number of components was assessed using the AIC calculated for models with different number of transitions 26 . It was found that at least three Lorentzian components were required to account for the spectrum for cases at saturating Ca 2+ concentrations.
Mechanism and parameter estimation. We hypothesized that TMEM16A gating can be described by the following mechanism, where C and O correspond to the closed and open states respectively, and the subscripts denote the number assigned to the states. The matrix notation of this mechanism 44 is where the subscripts indicate the transition described by the rate constant k in s −1 , for example, k 01 corresponds to the rate constant of the transition from state 0 to 1.
In the case where the Ca 2+ binding steps were included, the following mechanism, with the superscripts denoting the number of Ca 2+ bound, was used and the corresponding matrix is given by  2   6  6  6  6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  7  7  7  5 ð8Þ where x is the molar concentration of Ca 2+ . In both cases, the below relation was used to decrease the number of free parameters where L is the forward equilibrium constant with the subscript indicating the transition, and P o was supplied as an experimental estimate from non-stationary noise analysis at saturating Ca 2+ concentrations (P OCa→∞ ). At zero Ca 2+ , the open probability (P O0Ca ) was calculated by rearranging where I is the mean current, i the unitary current, and the subscripts indicate at zero Ca 2+ (0Ca) and at saturation (Ca ! 1), which were derived from experimental values. In the case of the full mechanism, microscopic reversibility 45 and the knowledge of the highest Ca 2+ binding affinity (K d(a2) ), which was empirically estimated to be 3.6 × 10 −8 M in an accompanying manuscript 13 , were also used and the on rate of the first Ca 2+ binding step was assigned a diffusion-limited value (1 × 10 11 ) 46 . Following the methods of Colquhoun and Hawkes 44 , the equilibrium occupancy of states was calculated from where Y 0 is the initial occupancy and V was obtained from the Eigen decomposition of Q Q ¼ VΛV À1 are the Eigenvalue and Eigenvector matrices respectively. The corresponding spectral matrices are given by The general form of the single-sided power spectrum due to Markovian fluctuations is given by 19 where N is the number of conducting units, V is the membrane potential when the current is linear and reverses at 0 mV, is the steady-state occupancy of open states 1 to k, is a submatrix of the spectral matrix and o 1 o k ; o 1 o k denote the upper left elements, is the conductance of the states arranged in a matrix form, and is a unit vector of length corresponding to the number of open states. Because the amplitude of the power spectrum concerns the number of channels and their conductance, which are variables not related to the mechanism, we fitted the experimental power spectra using a normalized form where G(0) is a constant corresponding to the power at a very low frequency.
For spectra obtained at saturating concentrations, model parameters (θ) were estimated by minimizing the sum of squares between Eq. 16 (4-state Q-matrix, concentration-independent) and the experimental spectra (y(f j )).
In the case where a family of power spectra was fitted, the sum of squares between Eq. 16 (7-state Q-matrix, concentration-dependent) and the experimental spectra obtained at the indicated Ca 2+ concentrations (y(x i , f j )) and those between Eq. 12 and the experimental open probability (p(x i )) were used to calculate the total sum of squares where z is a scaling factor. The variance of the best-fit parameters was obtained from the diagonal elements of the variance-covariance matrix 47,48 where H and J are the Hessian and Jacobian matrices at the least squares estimates respectively, the superscript T indicates transpose, and n d and n p are the number of data points and parameters respectively. The square root of the variance was used to approximate the standard deviation error, from which the 95% confidence interval was computed.
Validation of kinetic analysis via simulations. As discussed in the results section, the power spectrum provides an alternative means to analyze microscopic kinetics for channels where single-channel recording is challenging. Although mechanistic information can be inferred, direct fitting of the power spectrum is not possible without additional experimental information as the number of rate constants to be estimated is always larger than the number of observables (amplitudes and corner frequencies of the Lorentzian components) even for linear mechanisms. This limitation may be overcome by using the open probability (Po), which can be estimated independently from non-stationary noise analysis, as an additional constraint.
To validate this approach, we analyzed the accuracy and precision of parameter estimation when the independent experimental information is combined. As described below, we randomly sampled the power spectra of simulated singlechannel trajectories and Po from a Gaussian distribution with a spread similar to the experimental standard error, and estimated parameters from these synthetic data. The distributions and correlations of these estimates are shown in Supplementary Fig. 3. As expected from an adequately determined system, the 95% confidence intervals are finite and the estimates are centered at the true value of the parameter (Supplementary Fig. 3a). Moreover, although the estimates of some of the rate constants are somewhat correlated, the estimates of the equilibrium constants show little if any correlation ( Supplementary Fig. 3b, c). These results indicate that both the rate and equilibrium constants estimated from the power spectrum can be interpreted meaningfully for linear mechanisms.
In addition, as most experimental records are not truly stationary, we analyzed the effect of non-stationarity and/or current rundown, typical of TMEM16A current recorded from excised patches, on simulated power spectra. In our examples, we modeled channel rundown as a stochastic first-order decay in the number of channels (Supplementary Fig. 4a) and, as an alternative mechanism, a slow entry into a non-conducting state from the open state ( Supplementary  Fig. 4d). In both cases, the resulting power spectrum can be described by the superposition of the spectrum describing the mechanism and that describing the NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-20788-8 ARTICLE NATURE COMMUNICATIONS | (2021) 12:786 | https://doi.org/10.1038/s41467-020-20788-8 | www.nature.com/naturecommunications autocorrelation of the survival trajectory of the number of activatable channels ( Supplementary Fig. 4b, c, e, f). This is likely a consequence of the separation of time scales, meaning that the concurrent current decay does not affect the determination of model parameters when it is slow compared to the time scale of interest and that its power can be subtracted from the experimental spectrum to yield the mechanism-dependent components.
Analysis of the gating pathway. Similar to Colquhoun and Lape 49 , timeindependent probabilities for a transition out of a specified state were calculated from the optimized model parameters according to Specific to the mechanism that we proposed, two alternate routes of Ca 2+ binding can be taken during activation. One route consists of two consecutive binding events before the channel undergoes any gating transitions (6→5→3→2→1→0), which we term non-coupled binding, and the other is coupled to a conformational change (6→5→4→2→1→0), termed coupled binding. The probabilities for these routes are P nonÀcoupled binding ¼ π 65 π 53 π 32 π 21 π 10 P coupled binding ¼ π 65 π 54 π 42 π 21 π 10 ð28Þ and are dependent on Ca 2+ concentration. The relative usage of these routes was calculated using for the range of Ca 2+ concentrations displayed in Fig. 2d. The preferred route of deactivation was calculated as above, except that the corresponding sequence of transitions was flipped.
Double-mutant cycle analysis. The free energy of transition (ΔG) was calculated from the forward equilibrium constant using where R and T have their usual thermodynamic meanings, L is the forward equilibrium constant and the subscript indicates the transition from state i to j. A double-mutant cycle 29 The coupling energy between X and Y (ΔΔΔG XY ) was calculated using either the X or Y mutations The perturbation caused by the mutation X ! 0 on the wild-type background (X, Y) may be thought of as a combination of effects interdependent (ΔG XY ij ) as well as independent (ΔG X ij ) of the other residue Y in the cycle, hence Assuming that the inter-dependent component (ΔG XY ij ) is largely abolished when the same mutation X ! 0 is analyzed on a background when Y is mutated (X, 0), this perturbation may be interpreted as and the coupling energy as Consequently, the interdependent component (ΔG XY ij ) might be extracted from the coupling energy (ΔΔΔG XY ij ), and the independent effect (ΔG X ij ) may be approximated by ΔΔG ð0ÀX;0Þ ij . The standard error (σ) of the parameter estimates for each subtraction was propagated as described in the Statistics section. Deviation of ΔΔΔG XY ij from zero was detected using a two-sided one-sample t-test with a significance level of 0.05.
Rate-equilibrium free-energy relation analysis. The rate-equilibrium free-energy relation 34,51 consists of the following pair of relations that describe the effect of a series of perturbations on the rate constants (k f and k b ) as a fraction (ϕ and ϕ − 1) of their effect on the forward equilibrium constant L. ϕ can adopt values between 0 and 1. k i is the rate constant when L = 1. The parameters ϕ and k i were estimated by minimizing the total sum of squares for the set of equations for each transition.
Simulation of single-channel records. Trajectories of single-channel opening and closing (X t ) were generated using an approximation of a stochastic simulation algorithm 52 . The initial state sampled at t = 0 was assigned according to the deterministic equilibrium occupancy of the states (P 1 P n ) (Eq. 12). A random number (r) was drawn from the uniform distribution in the unit interval and the following acceptance criteria were used to determine which state is sampled, If r < P i then X t¼0 ¼ state i If P i ≤ r<P i þ P j then X t¼0 ¼ state j If P i þ P j ≤ r < P i þ P j þ P k then X t¼0 ¼ state k and so forth. Once initialized, the transition to the next state was assigned according to the transition probability in the infinitesimal time interval (dt) 44 , A random number (r) was again drawn from the uniform distribution in the unit interval and the following acceptance criteria were used to determine which state is sampled next, If r < p ij dt ð Þ then X tþdt ¼ state j given that Þ then X tþdt ¼ state l given that X t ¼ state i and so forth. This procedure was then repeated until the desired length of the trajectory is reached.
The above algorithm is a finite approximation to our choice of the infinitesimal interval dt. We used this approximation, instead of the exact stochastic simulation algorithm that directly samples the dwell-time distribution 52 , because the resulting trajectory is intended to serve as an input for FFT that requires equally spaced data points. Nonetheless, when a short time interval was used as in our case, the trajectories generated using the above algorithm were found to provide a relatively accurate description of the system's behavior, as we observed good agreement between their power spectra with the corresponding deterministic solutions.
Simulation and analysis of rundown. The trajectory of rundown was simulated using an algorithm similar to the above with the following modifications. We modeled channel rundown as a stochastic first-order decay of the number of conducting units with a decay probability in the infinitesimal time interval (dt), An array consisting of N single-channel trajectories was used to describe a channel population and the sum of these trajectories gives the ensemble current at steady-state. The number of trajectories was recorded as N and the ensemble current was calculated using A random number (r) was drawn from the uniform distribution in the unit interval and whether a decay occurs in the next time step was assigned according to the following acceptance criterion. If r < p decay N t ; dt ð Þthen N tþdt ¼ N t À 1, and one of the single-channel trajectories was removed according to a random integer drawn from the discrete uniform distribution in the interval (0, N t ). The ensemble current at t + dt was calculated by taking the sum of the surviving single channels whether or not a decay has occurred. This procedure was then repeated until the desired length of the trajectory is reached. As an alternative mechanism, we also modeled channel rundown as a slow entry into a non-conducting state from the open state. The trajectory was simulated as in the simulation of single-channel records, but the initial occupancy of states was assigned according to those calculated for only the closed and open states to mimic rapid equilibration before the current decays. In both cases, the resulting macroscopic current is characterized by a slow decay when k decay > 0, which reflects the decrease in the number of activatable channels over time.
The power spectra of the simulated ensemble current with or without rundown, from the same set of single-channel trajectories, and the survival trajectory were obtained via FFT. The one-sided deterministic power spectrum of the survival trajectory, modeled as an exponential decay, was calculated using where f is frequency, N and P o are defined as above, and T is the duration of the simulated trajectory.
Poisson-Boltzmann calculations. The electrostatic potential along a path connecting the lower Ca 2+ binding site and the cytosolic space was calculated by solving the linearized Poisson-Boltzmann equation in CHARMM 53,54 on a 240 Å × 240 Å × 260 Å grid (1 Å grid spacing) followed by focusing on a 160 Å × 160 Å × 160 Å grid (0.5 Å grid spacing). Partial protein charges were derived from the CHARMM36 all-hydrogen atom force field 55 . Hydrogen positions were generated in CHARMM. The membrane boundaries and dielectric properties of the system were as described previously 12 . The Ca 2+ -free and -bound structures of mouse TMEM16A (PDB: 5OYG and 5OYB respectively) were used to represent the Ca 2+ -sensing helix α6 in its resting and activated conformations respectively. Three configurations of sub-maximally bound states-resting α6 without any Ca 2+ bound, resting α6 with the upper Ca 2+ bound, and activated α6 with the upper Ca 2+ bound-were used to model the ligand-binding intermediates described in the full mechanism that was shown to adequately account for both Ca 2+ binding and gating in TMEM16A at steady-state.
Monte Carlo estimation of confidence intervals. N single-channel records were simulated using the parameters estimated for the wild-type channel at a sampling frequency of 20 kHz. 1.5% of the records were randomly drawn from the discrete uniform distribution in the interval (0, N) to mimic random sampling during an experiment. N was typically 1000. The averaged power spectrum of the randomly sampled records was calculated and model parameters were estimated by minimizing the sum of squares. To account for uncertainty in the P O estimated from nonstationary noise analysis, a random number drawn from the Gaussian distribution centered at the P O of the wild-type channel with a standard deviation of 0.02 (which typically results in values ranging from~0.77 to~0.88 and centered at 0.82) was used as an input for each fit. This sampling procedure was then repeated 1000 times to estimate the uncertainty of each rate and equilibrium constant and their correlation.
Statistics. Data were collected from individual cells obtained from different transfections. Data analysis were performed using Clampfit 10.6 (Molecular Devices), Excel (Microsoft), NumPy (https://www.numpy.org), and SciPy (https:// www.scipy.org). For numerical operations and simulations, NumPy and SciPy were used. Parameter optimization was performed with the least_squares function in SciPy using the described sum of squares objective function, which also computes the Jacobian matrix that was used to estimate the 95% confidence intervals. Standard error uncertainties were propagated using 56 σ ðaþb or aÀbÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi σ 2 a þ σ 2 b q σ ðab or a=bÞ f ða; bÞ j j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi σ a a j j