Mapping allosteric linkage to channel gating by extracellular domains in the human epithelial sodium channel

The epithelial sodium channel (ENaC) mediates sodium absorption in lung, kidney, and colon epithelia. Channels in the ENaC/degenerin family possess an extracellular region that senses physicochemical changes in the extracellular milieu and allosterically regulates the channel opening. Proteolytic cleavage activates the ENaC opening, by the removal of specific segments in the finger domains of the α- and γ ENaC-subunits. Cleavage causes perturbations in the extracellular region that propagate to the channel gate. However, it is not known how the channel structure mediates the propagation of activation signals through the extracellular sensing domains. Here, to identify the structure–function determinants that mediate allosteric ENaC activation, we performed MD simulations, thiol modification of residues substituted by cysteine, and voltage-clamp electrophysiology recordings. Our simulations of an ENaC heterotetramer, α1βα2γ, in the proteolytically cleaved and uncleaved states revealed structural pathways in the α-subunit that are responsible for ENaC proteolytic activation. To validate these findings, we performed site-directed mutagenesis to introduce cysteine substitutions in the extracellular domains of the α-, β-, and γ ENaC-subunits. Insertion of a cysteine at the α-subunit Glu557 site, predicted to stabilize a closed state of ENaC, inhibited ENaC basal activity and retarded the kinetics of proteolytic activation by 2-fold. Our results suggest that the lower palm domain of αENaC is essential for ENaC activation. In conclusion, our integrated computational and experimental approach suggests key structure–function determinants for ENaC proteolytic activation and points toward a mechanistic model for the allosteric communication in the extracellular domains of the ENaC/degenerin family channels.

Members of the epithelial sodium channel (ENaC) 3 /degenerin ion channel family share functions that can be backtracked to a common ancestor. This ancestral receptor was likely essential to the maintenance of electrolyte balance in eukaryotic cells as they adapted to an electrolyte-rich environment (1). The ion channels in this family control and regulate the flux of cations entering cells, thus the channels control the cell volume and the electric properties of the cell membrane. Maintaining proper cell volume is essential for cellular and organismal survival (2). Therefore, the ENaC/degenerin ion channel family is crucial to electrolyte homeostasis in organs such as the kidneys, lungs, and brain (3). Mutations in channel subunits can cause severe electrolyte disorders, such as hereditary hypertension Liddle syndrome and neuronal degeneration in Caenorhabditis elegans (1,(3)(4)(5)(6)(7)(8). Although the channels in the ENaC/degenerin family share considerable homology, each channel has adopted a unique biochemical mechanism to regulate its function.
All of the channels contain an extracellular degenerin site that is essential for proper channel function. The domains in the extracellular region sense and process several extracellular clues, such as changes in the concentrations of electrolytes, neuropeptides, and proteases (9). It remains unknown, however, how the detection of extracellular signals modifies the operation of the distal gating mechanism. The lowest homology among the channel subunits resides in the finger domain, specifically the hypervariable region (HVR). The HVR has diverged to adopt specific folds that enable each channel to sense and respond to specific extracellular signals (10). The HVRs in the ␣and ␥ ENaC-subunits are partially disordered, and they undergo site-specific channel activation by proteolytic cleavage that releases inhibitory peptides (11). Because perturbations in the finger domains are essential to couple the activation signals to the channel gate, we hypothesized that key extracellular residues and inter-residue interactions transduce HVR conformational changes to the gate region.
The ENaC/degenerin ion channel family utilizes a remarkable allosteric machinery that enables the large extracellular domains to act as receptors to sense specific changes in the environment (12,13), such as the release of neuropeptides in the extracellular milieu (1). Transduction of extracellular signals into conformational changes at the gating pore implies dynamical coupling in the protein network (14). Better understanding of how extracellular signals are processed in the ENaC/degenerin family has been challenging due to limited structural data (15). The X-ray structures of the chicken acid sensing ion channel (cASIC) in multiple functional states have provided a platform for further structure-function studies of the allosteric mechanisms that regulate the function of the ENaC/degenerin ion channels (10, 16 -21). Results of previous computational work provided insights for developing mechanistic models that describe the inter-domain motions in ASIC subunits resulting from extracellular proton binding (19,(22)(23)(24)(25). Particular residues in the ASIC palm domain have been found to modulate the kinetics of desensitization upon stimulation by proton binding. Other investigators have used cysteine mutagenesis and chemical modification to identify functionally important extracellular sites (18). Here, we developed an integrated computational and experimental strategy to define a structural basis for the allosteric activation of ENaC by proteolytic cleavage. We used a graph-based approach to analyze the allosteric pathways that are derived from inter-residue dynamical coupling obtained from molecular dynamics simulations. We related the changes in the structural dynamics of ENaC upon proteolytic cleavage to the kinetics of proteolytic activation.
We have combined molecular dynamics simulations, chemical modifications of residues substituted by cysteine, and electrophysiology recordings to investigate the structure-function determinants of ENaC proteolytic activation. We identified residues in the ␣-subunit that modulate the kinetics of proteolytic activation and uncovered allosteric pathways that dynamically link the HVR region in the ␣ENaC finger domain to the channel gate. Our molecular dynamics simulations of ENaC heterotetramer ␣ 1 ␤␣ 2 ␥ predict the conformational changes and intersubunit coupling that mediate ENaC activation. In summary, our findings offer a dynamic model of the structural basis for the allosteric sensing and processing of extracellular activation signals by the ENaC/degenerin family.

Cleavage-triggered inter-subunit conformational changes alter intra-subunit dynamics
Because ENaC cleavage occurs mainly in the ␣and ␥-subunits of the protein, we focused our analysis on the interactions at the ␣␥ interface. By analyzing multiple discrete molecular dynamics (DMD) trajectories (26,27), we found that cleavages in the ␣and ␥-subunits relax the finger domains of the subunits, because the residues flanking the cleavage sites fold into a set of thermodynamically favorable conformations (Movie S1). All of the sampled conformations fell into an ensemble of structures that were stabilized by H-bonded interactions formed by the finger domain of the cleaved subunit (␣ or ␥) and the knuckle domain of the interfacing subunit (␣, ␤, or ␥). We observed interactions between the knuckle domain of ␣ENaC and the finger domain of ␥ENaC that occur subsequent to cleavage. Upon relaxation of the flanking residues in the ␥-subunit, the positively charged residues at the furin site, 135 RKRR 138 , formed stabilizing contacts with neighboring polar residues. Mainly, the observed inter-subunit interactions are mediated by residues in the ␣-knuckle domain, such as ␣Asp 527 and ␣Glu 531 . These acidic residues formed H-bond contacts with the lysines and arginines in the ␥-furin site (Fig.  1A).
Upon cleavage, the residues flanking the furin cleavage site in ␥ENaC undergo major relative displacements, up to ϳ4 Å (Fig.  1C). The second major conformational change occurs in the knuckle domain of the ␣-subunit. As the furin site in the ␥-subunit folds into a stable state, the polar residues in the knuckle domain of the ␣-subunit undergo displacements that are quantitated by an increase in the root mean square fluctuation (r.m.s. fluctuation) of the C␣ atoms (Fig. 1, A and C). These observations suggest that inter-subunit dynamical coupling stabilizes intermediate structures that can promote activation of ENaC upon proteolytic cleavage. To understand the allosteric impact of the observed conformational changes on the intra-subunit dynamics, we constructed a graph, residue interactions networks (RINs), representing the protein structure.
We used the contact frequency data derived from DMD simulations to assign weights to the constructed graph edges (13,14,28). We then performed simulations of random walks on the constructed networks using the MapEquation method that utilizes the clustering algorithm Infomap to identify communities of inter-connected residues (see "Materials and methods"). The communities found using the MapEquation method ( Fig. 2) resemble closely the domain architecture suggested by the crystal structure of the cASIC (17). The knuckle and the palm domains are part of the largest module, community I, with Ala 549 as its hub. The wrist domain is shared between community I and community VI (Fig. 2, A and C). The boundaries between communities may facilitate the exchange of dynamical coupling across the network through the correlated inter-residue fluctuations. The set of inter-community links can serve as a guide for understanding the communication pathways in ␣ENaC.
We used the fluctuation-based cross-correlation data derived from DMD simulations to compare the inter-residue dynamical coupling between the cleaved and uncleaved RINs. To identify mechanisms of information transfer upon cleavage, we have analyzed two information transfer metrics: the spatial distribution of unweighted shortest paths and the flow of correlations through paths in the network. From DMD simulations, we observed local conformational changes at the cleavage sites of the ␣and ␥-subunits (Fig. 1). We hypothesized that the conformational changes accompanying proteolytic cleavage cause a reconfiguration of the RIN. To test this hypothesis, we aligned the two networks and, indeed, we found that upon cleavage, the uncleaved RIN gained and lost some edges due to changes in inter-residue contacts. Moreover, cleavage allowed for additional communication through an additional set of paths (Fig. 3B). The differences in the distribution of contacts in the protein structure are reflected in the network parameters and configuration (Fig. 3, A and B). We used the network analysis algorithm in Cytoscape to compute the network parame-

Allosteric coupling in ENaC extracellular domains
ters in the cleaved and uncleaved RINs. The characteristic path length in unweighted networks reflects how information can be optimally transferred through the network (28). By analyzing the average path lengths between the source and target residues (see "Materials and methods") in the cleaved and uncleaved networks, we found that cleavage changes the network configuration resulting in a reduced average path length. Consequently, the dynamical perturbations due to cleavage may pro-

Allosteric coupling in ENaC extracellular domains
vide more efficient information transfer in the allosteric network of ␣ENaC.

Allosteric pathways that can propagate activating signal in ␣ENaC
To choose a source node or residue, which can sense the allosteric communication triggered by cleavage, we analyzed the pairwise correlations between the residues close to the cleavage site. As expected, these residues exhibit a significant change in the cross-correlation upon cleavage. Similarly, a large cluster of residues in the knuckle domain is also dynamically perturbed upon cleavage (Fig. 2B). Additionally, using the MapEquation (see "Materials and methods") method, we calculated the flow of random walks through the RIN. We found that Ser 270 has the highest visit frequency among the residues with maximal perturbations upon cleavage. The random walk analysis suggests that Ser 270 contributes strongly to the communication within its community. Thus, we used Ser 270 as the source of signal transmission. Based on the crystal structure of the cASIC, the largest conformational change of activation occurs in the wrist domain (29). We chose the target residues to be where the largest conformational changes occur in the ENaC/ degenerin family. Therefore, we analyzed the allosteric paths that communicate changes in the dynamical coupling between Ser 270 and residues that lie at the inter-domain interface between the wrist and palm domains.
Considering that cleavage can alter not only the spatial distribution of inter-residue contacts, but also the magnitude of the dynamical coupling between residues, we analyzed the weighted RINs of ␣ENaC, where each edge was weighted according to the inter-residue pairwise correlations obtained from DMD simulations. For these networks, the optimal paths were ranked with respect to the flow of correlations along the paths, i.e. how inter-residue correlated motions propagate along each path (see "Materials and methods"). We used PathLinker to search for and rank paths according to the correlation-based flow. By analyzing the spatial distribution of the paths sampled by PathLinker, we found that both the cleaved and uncleaved networks share one main cluster of optimal paths (Fig. 4A). The optimal paths pass through communities V, III, and I. Although the clusters of paths are identical in terms of the communities they pass through, the flow of the top optimal paths is higher in the cleaved state than in the uncleaved state (Fig. 4B). The increase in flow of the ensemble of optimal paths upon activation by proteolytic cleavage suggests that the reconfiguration of the allosteric network is mediated in part by the stabilization of interactions along the optimal paths.  Right, schematic representation of the ␣-subunit. The cleavage fragment is colored black. One cluster of optimal paths is represented on the structural model using PyMOL. A gray sphere indicates each residue on the path. Edges connecting residues on each path are modeled as cylinders. Edges are colored yellow in the optimal paths. The target site in the wrist domain is colored orange. B, a plot of the natural algorithm of flow for the cleaved (gray) and uncleaved (blue) networks. The horizontal axis shows the index of each sampled path from the PathLinker algorithm. The indices range from 1 to 100,000, so that a path with index k ϭ 1 has the highest flow, and path with index k ϭ 100,000 has the lowest flow in the top 100,000 paths obtained from the PathLinker algorithm.

Allosteric coupling in ENaC extracellular domains Lower palm domain in ␣ENaC is essential for ENaC proteolytic activation
Because we found dynamical differences between the uncleaved and cleaved networks in the ␣-subunit, we hypothesized that there are sites in ␣ENaC that regulate the structural dynamics needed for activating the channel. Moreover, our DMD simulations predicted a functional role for interactions at the ␣␥ inter-subunit interfaces. To uncover dynamically sensitive sites that can propagate the activation signal, we used sitedirected mutagenesis to substitute residues with cysteines and chemically modified these cysteines during electrophysiology recordings. We used thiol modification to probe the functional impact of chemical perturbation of specific sites on the protein structure. We selected sites that lie at the inter-subunit interfaces of our model for ENaC heterotetrameric complex. This approach allowed us to monitor, in real-time, the effect of sitespecific chemical perturbations of ENaC subunits on the channel activity. First, we recorded the basal activity of ENaC by monitoring the amiloride-sensitive currents. Then we perfused the oocytes expressing cysteine-engineered ENaC subunits with a cysteine-reactive MTS reagent, such as MTSET (Fig.  5A). As a negative control, we recorded ENaC currents from oocytes expressing WT channels. As expected, the function of WT channels was insensitive to MTSET perturbations (Fig.  5A). However, MTSET modification of a site in a loop in the lower palm domain (Fig. 5B), ␣Glu 557 , which belongs to the largest community in the protein network (Fig. 2C) of the ␣-subunit, elicited a marked inhibition of basal activity (Fig. 5A). We found similar results for ␣Leu 558 , but not for ␣Lys 556 (Fig. 5C). The observed inhibitory response suggests that ␣Glu 557 and ␣Leu 558 may be important for the stabilization of a closed state of ENaC.
To determine whether the MTSET-triggered inhibition interferes with the proteolytic activation of ENaC, we perfused the oocytes with 10 g/ml of chymotrypsin to proteolytically stimulate the channel. To isolate post-cleavage events that maybe related to the rate-determining conformational changes leading to activation of near-silent channels, we applied a brief protease treatment (30 s) followed by excess aprotinin (Fig. 6A). To ensure that aprotinin effectively stopped continuing chymotrypsin cleavage of ENaC, we confirmed that co-perfusion of chymotrypsin and aprotinin for 2.5 min led to no stimulation. Therefore, the cleavage events that activate ENaC occur within the first 30 s of chymotrypsin exposure. We infer that the slow stimulation of ENaC by chymotrypsin monitors events set in motion by cleavage, such as the molecular rearrangements that activate the channel.  MTSET does not impact neither the basal activity nor the rate of proteolytic activation. MTSET treatment of ␣E557C slows down the rate of proteolytic activation, asterisks indicate statistically significance differences between ␣E557C groups, p Ͻ 0.0001 using unpaired two-tail t test. Data shown are from 8 to 9 oocytes repeated in 2-3 batches.

Allosteric coupling in ENaC extracellular domains
We hypothesized that the exponential activation rate of ENaC currents in a single oocyte is an ensemble record produced by thousands of channels undergoing sequential conformational changes that culminate in channel opening. Surprisingly, we found that the MTSET treatment slows the rate of proteolytic activation of ␣Cys 557 by ϳ2-fold, but had no effect on stimulation of WT ENaC (Fig. 6A). Absent MTSET, the time constant of activation generated by brief protease stimulation, , is ϳ50 s for WT ENaC and the Cys mutant ␣Cys 557 /␤/␥. Applying MTSET to oocytes expressing ␣Cys 557 /␤/␥ increases to ϳ100 s (Fig. 6B). We found that the fold-increase in current above basal activity was comparable between MTSET-treated and untreated groups. This slowed activation due to covalent modification of ␣Cys 557 suggests that this residue is dynamically coupled to both the cleavage site and the channel gate. Interestingly, this site belongs to the ensemble of optimal paths predicted from DMD simulations (Fig. 4A), suggesting that the optimal paths may be important for the stabilization of nearsilent channels, which remain mostly closed until an activating cleavage event occurs. The presented data suggests that ␣Glu 557 is critical for propagating rate-determining conformational changes that are required for the proteolytic activation of ENaC.
The results from DMD simulations and electrophysiology experiments suggested a novel functional role for the lower palm domain of the ␣-subunit in ENaC proteolytic activation. Specifically, the loop containing residues ␣Glu 557 and ␣Leu 558 is functionally critical for ENaC proteolytic activation. We tested whether the homologous loops in the ␤and ␥-subunits can exhibit a similar functional response upon MTSET perturbation. Interestingly, neither ␣/␤Q498C/␥ nor ␣/␤/␥D511C channels were affected by MTSET treatment in the whole cell voltage-clamp experiments (Fig. 5C). MTSET treatment did not affect basal activity of ENaC currents for the ␤ and ␥ mutant channels, but ␣E557C/␤/␥ and ␣L558C/␤/␥ channels showed inhibitory responses to MTSET treatment. Our structural model of ENaC suggests that the homologous sites are solvent exposed. In addition to the loop in the lower palm domain of the ␣-subunit, we have tested an additional site in the palm domain of the ␣-subunit. ␣D148C/␤/␥ channels were irresponsive to MTSET perturbation with an I-MTSET of 1 (Fig. 5C). These data suggest that among all the tested sites in the ␣-, ␤-, and ␥-subunits, ␣Glu 557 and ␣Leu 558 , when mutated to cysteines, are the only sites that are responsive to MTSET treatment. Moreover, the rate of proteolytic activation is reduced for ␣Glu 557 upon dynamical perturbation of these tested sites.

Discussion
We propose a structural model for the activation of ENaC by proteolytic cleavage in the ␣and ␥-subunits. We identify key allosteric pathways and residues that facilitate the propagation of proteolytic events to the channel gate. We used a combination of computational and experimental approaches to uncover key structure-function determinants for ENaC proteolytic activation. Based on the data generated from DMD simulations of the ENaC heterotetramer ␣ 1 ␤␣ 2 ␥ in the cleaved and the uncleaved states, we built residue interaction networks to identify allosteric pathways that link the cleavage state of ENaC to the gate region of the channel. We substituted amino acid residues with cysteines and used thiol modification agents to identify among the substituted residues those that modulate the rate-determining conformational changes of channel activation. By integrating the predictions from molecular dynamics simulations and the experimental results from electrophysiology recordings, we derived a mechanism for ENaC proteolytic activation. Our model specifies that the lower palm domain of the ␣-subunit of ENaC is critical for propagating the proteolysis-triggered conformational changes, initiated by the cleavage of the prostasin site in ␥ENaC, to the gate region. To our knowledge, this work presents the first evidence for the role of the lower palm domain of ␣ENaC in modulating the kinetics of channel gating.
To determine the structure-function determinants of allosteric activation of ENaC, we analyzed the differences in dynamics in the allosteric networks of ␣ENaC upon cleavage. Network analysis of protein structures using graph-based approaches can give insight to allosteric sites and mechanisms that can regulate and diversify protein function. Recently, engineering studies of mechanically coupled networks have used allosteryinspired approaches to find strategies to create new functionalities in networks (30,31). Several structure-function network parameters have been previously used to characterize information flow in networks. Shortest path analysis has been instrumental in understanding the mechanisms of long-range communications in allosteric systems such as ligand-activated proteins (32). del Sol et al. (33) used shortest path analysis of protein structures from several protein families. The investigators found that residues that maintain the connectivity in protein networks are conserved. Our network analysis of structures of ␣ENaC in the uncleaved and cleaved states revealed a decrease in the average shortest path length in the ␣ENaC network. Changes in the average shortest path length suggest a rewiring in the paths that can communicate the activation signal to the channel pore region.
The inter-subunit conformational changes induced by proteolytic activation of ENaC cause changes in the intra-subunit dynamics of ␣ENaC. To explain the altered inter-residue interactions upon cleavage, we analyzed the inter-residue dynamic coupling in paths that connect the cleavage site to the pore region. A set of dominant paths spanning the ␤-ball, knuckle, and palm domains exhibited higher inter-residue correlations upon cleavage. del Sol et al. (32) have suggested that the dynamical properties that regulate the allosteric networks of ligandactivated proteins are tuned by an ensemble of pre-existing pathways. Our results from analyzing the flow of correlated fluctuations in the allosteric network of ␣ENaC further support this conformational ensemble model. Consistent with previous work by Berman et al. (34), our molecular dynamics simulations predict formation of stabilizing inter-residue contacts between the knuckle domain of the ␣-subunit and the finger domain of the ␥-subunit. The positively charged residues at the ␥-cleavage site, RKRR, undergo conformational changes upon cleavage and form H-bonds with the acidic residues in the ␣-knuckle domain. Interestingly, similar polar contacts form between the knuckle domains and the finger domains of additional inter-subunit interfaces, such as ␣␤ and ␥␣. We have

Allosteric coupling in ENaC extracellular domains
found that perturbations in the knuckle domain dynamics are important for the activation of ENaC. Moreover, disease-associated substitutions in the ␣-subunit knuckle domain result in changes in the inter-subunit structural dynamics of ENaC (35). Here, we found a loop in the lower palm domain of ␣ENaC to be critical in linking the cleavage sites to the channel gate. Interestingly, introduction of cysteines at the homologous residues in the ␤and ␥-subunits did not generate a functional response upon MTSET perturbations. These results suggest the lower palm domain in ␣ENaC may have unique dynamic properties to regulate ENaC gating. Moreover, the lack of inhibition response in the ␤and ␥-subunits may indicate the asymmetric nature of the structural dynamics in ENaC subunits. In summary, we present an allosteric model of ENaC proteolytic activation (Fig. 7). To identify critical residues that regulate the kinetics of proteolytic activation of ENaC, we combined computational modeling with two-electrode voltage-clamp electrophysiology recordings of chemically modified channels. We found mechanisms with which the allosteric network in ␣ENaC can be dynamically perturbed by proteolytic cleavage. This work explains the conformational transitions and interactions that allow ENaC to cross the thermodynamic barriers to conduct sodium. The changes in structural dynamics by cleavage perturb the dynamic coupling in the residue interaction networks. The slow kinetics of ENaC proteolytic activation point to the remarkable allosteric machinery required for the channel function. Future studies should address the conformational transitions that result from specific cleavage events. How dynamically and functionally reactive is the protein network to perturbations by single cleavage events? What is the extent of asymmetry in dynamics in the different ENaC subunits, ␣, ␤, and ␥? Fluctuations in allosteric networks can be a driving force to allow for conformational transitions that activate ion channels in the ENaC/degenerin family. Our work advances the structural and dynamic understanding of the microscopic variables that tune the kinetics of the conformational transitions required to activate ion channels in the ENaC/degenerin family.

Electrophysiology recordings
Xenopus oocytes were harvested and maintained in 85 mM NaCl, 1 mM KCl, 2.4 mM NaHCO 3 , 0.82 mM MgSO 4 , 1.88 mM CaCl 2 , 0.33 mM Ca(NO 3 ) 2 , 16.3 mM HEPES, pH 7.4, and 10 M amiloride. The protocols used to obtain oocytes from Xenopus laevis frogs are approved by the University of North Carolina IACUC (protocol ID 17-261.0). Oocytes were injected with 0.3 ng of each cRNA encoding ENaC subunits ␣, ␤, and ␥. Twoelectrode voltage-clamp recordings were performed 24 h after injection. ENaC currents were recorded using Geneclamp amplifier (Axon Instruments) in a constant perfusion system. Oocytes were voltage clamped at Ϫ100 mV. Using a Digidata 1200 A/D converter (Axon Instruments) and AxoScope software, currents were digitized and recorded. We analyzed recorded traces using the pCLAMP software. Recordings were initiated in the presence of 10 M amiloride. Following amiloride wash, basal currents were recorded in frog Ringer's solution (2.5 mM KCl, 1.88 mM CaCl 2 , 120 mM NaCl, and 10 mM HEPES, pH 7.35) to obtain amiloride-sensitive currents I Na . Experimental groups of 5 to 6 oocytes were used, repeated in 2-4 separate oocytes batches.
We analyzed the kinetics of proteolytic activation by fitting the recorded trace of protease-mediated ENaC current to a single exponential function. We have shown before that the rate of exponential activation in recorded ENaC currents is related to the proteolytic susceptibility of activation (36). We normalized the time constant for proteolytic activation that is obtained from the MTSET-treated oocytes to the obtained from oocytes with no MTSET treatment. We refer to this normalized rate constant as normalized -MTSET. A -MTSET of 1 indicates no MTSET effect on the proteolytic activation kinetics, and -MTSET Ͼ1 indicates slowing down of proteolytic activation rate. Similarly, we calculated the fraction of currents responsive to MTSET as I-MTSET. I-MTSET Ͻ1 indicates an inhibitory effect of MTSET on basal activity.

Modeling ENaC proteolytic cleavage
We performed DMD simulations (26,27,(37)(38)(39) of a heterotetramer model of ENaC (35) to catalogue structural changes accompanying ENaC proteolytic cleavage. From multiple DMD trajectories, we analyzed the backbone fluctuations of ␣ENaC in the fully cleaved and uncleaved states. In different DMD trajectories, we observed cleavage-dependent conformational changes. To simulate the complete proteolytic activation of ENaC, we removed the inhibitory tracts in the ␣and ␥-sub-

Allosteric coupling in ENaC extracellular domains
units. Cleavages in both ␣and ␥-subunits are required for full activation of ENaC. Thus, we analyzed the structural dynamics at the ␣␥ interface to uncover possible inter-subunit dynamic coupling resulting from release of the cleavage fragments. Proteolytic cleavage occurs in the hypervariable region in the finger domains of the ␣and ␥-subunits. Thus, we considered the interacting residues surrounding the cleavage site in the finger domains as the sources of signal propagation for proteolytic activation.

Discrete molecular dynamics simulations
The starting structure of the ENaC heterotetramer was prepared and minimized as described previously (35). Briefly, relaxation and equilibration simulations were performed at a low temperature to minimize the potential energy of the system. Seven independent DMD trajectories were simulated at temperature (T) ϭ 0.4 DMD reduced units (kcal/(mol of k B )) for 10 6 DMD steps (37), which is approximately equivalent to 50 ns of simulation time for each trajectory. Harmonic constraints were applied on all backbone atoms in the transmembrane helices. A soft potential with a force constant of 0.1 kcal mol Ϫ1 Å Ϫ2 was applied to restrain the motion of the transmembrane helices to mimic lipid-protein interactions.

DMD simulations analysis
To analyze the structural dynamics in the extracellular regions, we obtained the equilibrated structures from each DMD trajectory after discarding the first 25 ns of equilibration. To analyze the dynamical coupling between residues in contact, we used in-house scripts to calculate the r.m.s. fluctuations of C␣ atoms and the cross-correlations between the C␣ fluctuations. Based on the data obtained from the trajectories we constructed a graph G(V, E), where V nodes correspond to the C␣ atoms of each residue in the ␣-subunit. The nodes in the graph are linked to each other by a set of edges or links E; each edge e ij connects a residue pair v ij that are in contact. A contact between two residues was assigned if the C␣-C␣ distance between two residues of less than or equal to 7.5 Å was maintained for at least half of the simulation time for each trajectory (13,40). We obtained the absolute value for the cross-correlation values and averaged them across the different trajectories. The average cross-correlation values were used to define the weights for edges (E) that connect the nodes or vertices (V) in the graph G(V, E). This graph was used to describe the RIN. We visualized the RINs using Cytoscape.
We used the last 25 ns from DMD trajectories to analyze networks corresponding to cleaved and uncleaved proteins. We computed the change in the absolute value of the ␣ENaC C␣ cross-correlation values upon cleavages in the ␣and ␥-subunits. We used the common nodes and edges in the uncleaved and the cleaved states to map the nodes that are most perturbed upon proteolytic activation. As these two states correspond to inactive or closed and active or open states of the ENaC, respectively, comparing properties of two corresponding networks may uncover allosteric determinants for the activation in the ENaC heterotetramer.

RIN analysis
To find out how the observed differences in fluctuations can propagate from the cleavage site to the pore region, we analyzed the structural origins of information transfer for each network. To identify the loopless optimal paths that connect the cleavage site in the ␣-subunit to the wrist domain, we employed the algorithm PathLinker to find the k highest-scoring paths (28). The details of the algorithm are described in Ref. 29. Briefly, PathLinker utilizes a modified and optimized version of Yen's algorithm that iteratively uses the Dijkstra's algorithm as a subroutine. Two types of rankings were used during optimal path search: the total number of nodes along the path for the unweighted graph and the amount of flow through each path. Flow was defined as a probabilistic measure to predict the likelihood that fluctuations will propagate through a certain path from the source node to the target node. Flow ⍀ k was quantified as the multiplication of all the edge weights on each path, ⍀ k ϭ ⌸ iϭN iϭ1 c i,iϩ1 , whereas c i,iϩ1 is the average cross-correlation value between two consecutive residues c i and c i,iϩ1 , on path with index k. For unweighted shortest-path analysis, we used the Pesca application in the Cytoscape program to find the tree of paths with the minimal number of edges (41). Optimal paths were visualized using PyMOL to identify whether the paths clustered into groups on the protein structure. Unweighted shortest paths were visualized using Cytoscape (42). We used the network analysis algorithm in Cytoscape to compute the network parameters in the cleaved and uncleaved RINs. We analyzed the distribution of the shortest-paths that connect residues that are in direct contact. To analyze the distribution of the flows between dynamically coupled residues in each network we obtained the top 100,000 optimal paths in each network using the algorithm PathLinker. The algorithm calculates the top-ranking k paths according to the flow of correlations from a source node to a target node. We define each path by two structure-function descriptors: (i) the residues and communities that the path passes through, and (ii) the amount of flow of correlated fluctuations through the path. The flow in a path describes the relative probability of the fluctuations propagating from the cleavage site to the pore region.

Community clustering analysis
We clustered the protein regions into communities of interconnected residues using the Mapequation method, which utilizes the Infomap algorithm (42,43). This method samples the protein network using random walks. The random walker defines the size of the community and the boundaries between communities according to the visit frequency and the number of links or edges each residue has with its neighbors. The trajectories of the random walks resemble the flow of information in the network (43). We mapped the communities identified by this method on the structural model of ␣ENaC. We used the intercommunity boundaries to guide the analysis of the paths found by the PathLinker algorithm (42). Residues in each community are more dynamically coupled to their neighbors than to residues in other communities. The size of each community is characterized by the flow of a random walker within the module relative to the visit frequencies in the rest of the modules in Allosteric coupling in ENaC extracellular domains the network (43). The inter-community link resembles how freely information is transferred between adjacent communities. We rank each community based on its total flow.