Our analysis of the CRAC channel began with the equilibrium simulation of the closed state (PDB ID: 4HKR). The functional state of an ion channel is customarily predicted (Rao et al., 2019) on the grounds of the radius profile. Supplementary Figure S1 shows the radius profile of the crystal structure and the profile of the radius averaged along the equilibrium simulations. Both profiles are everywhere larger than the radius of a single water molecule showing the absence of steric obstruction. From a purely geometric point of view the channel is thus expected to be water-permeable. This prediction is confirmed by the space filling representation of the pore shown in Supplementary Figure S2. According to the color code employed by the HOLE program (Smart et al., 1996), green regions are wide enough to accommodate no more than a single water molecule, while blue regions can accommodate two or more water molecules. Supplementary Figure S2 shows a bottleneck at the level of the Glu-ring and other two at the level of the basic region. In these points, the side chains project towards the center of the pore reducing its radius but there is still enough space for a single chain of water molecules. Except for these three points the pore can everywhere be occupied by two or more water molecules.

This analysis is inconsistent with the electrophysiological characterization of this structure (Hou et al., 2012) clearly showing it corresponds to the closed state. This inconsistency can be explained if the radius profile is compared with the hydrophobicity profile and the Potential of Mean Force of water permeation (Figure 2). Figure 2 shows that, even if in principle the pore is wide enough to host one or more water molecules (Figure 2A), the chemistry of the wall determines a peak of the hydrophobicity profile (Figure 2B) corresponding to the hydrophobic region of the pore. The hydrophobicity peak also corresponds to the maximum of the Potential of Mean Force for water permeation (Figure 2C). Taken together, these results suggest that the pore is functionally occluded even in the absence of steric block, which is the signature of hydrophobic gating.

Hydrophobic gating is a phenomenon of evaporation in conditions of nanoconfinement (Roth et al., 2008). If the water-wall interactions are weaker than water-water interactions, a bubble may form, preventing the flow of liquid water and ions. This pattern is clearly manifested in our equilibrium simulation. Figure 3 shows the time evolution of the number of water molecules in the whole pore region (Figure 3A) and in the hydrophobic region (Figure 3B). The water count in the whole pore slightly increases during the first half of the simulation and eventually it stabilizes around an average value of 200 water molecules. Conversely, if we restrict our attention to the hydrophobic region, this district appears to be almost always devoid of water safe for a few occasional intrusions of 1 or 2 water molecules at the boundary with the Glu-ring and the basic region. Indeed Figure 3D, corresponding to the final frame of the simulation, shows that the hydrophobic region of the pore is dry. The bubble appears early, already during the constrained stage of the equilibration, and it is stably maintained throughout the whole simulation of 100 ns. As a result, the radial-axial PMF of water permeation reported in Figure 3C, highlights a high free energy zone (in red) in the hydrophobic region where water has very little likelihood to reside.

The main purpose of this work is to study the gating transition steering the CRAC channel from the open to the closed state and back using Targeted Molecular Dynamics simulations (Schlitter et al., 1994). Since this technique requires knowledge of the end states of the transition, it is extremely important to computationally characterize the recently resolved structures of two putatively open states, mutants H206A (PDB ID: 6BBF) (Hou et al., 2018) and P288L (PDB ID: 6AKI) (Liu et al., 2019). This analysis will allow us to choose the most appropriate conformation to use as the open state in TMD simulations.

Figure 4 shows the radius profile of the crystal structures of the closed and open states (Figure 4A) as well as the radius profile averaged during the equilibrium simulation of the three systems (Figure 4B). The radius profiles of the experimental structures are very similar at the level of the Glu-ring, but in the hydrophobic and basic regions the pore of 6BBF appears to be much larger than that of the other two structures. Interestingly enough, even if 6AKI was predicted (Liu et al., 2019) to represent the open state, its radius profile is much more similar to that of the closed state (4HKR) than to that of the other putatively open state (6BBF). During the equilibrium simulation, the radius profile of 6AKI tended to drift towards that of 6BBF until it settled on a curve intermediate between those of the closed state (4HKR) and of the other putatively open state (6BBF). As shown in Figure 4C, the different pore geometries affect the Potential of Mean Force of water permeation. As expected, the closed state is characterized by a high free energy barrier at the level of the hydrophobic region where the vacuum bubble prevents the flow of liquid water. Conversely, during the equilibrium simulation the wide pore of 6BBF is completely water filled. Since water density inside the pore is comparable to that in the bulk, the PMF profile is almost flat. Finally, the PMF profile of 6AKI also shows a barrier, smaller than that of the closed state but higher than what expected for an open state. The origin of this barrier can be understood monitoring the number of water molecules in the hydrophobic region of the three systems (Figure 4D). The plot clearly shows that, during the first 40 ns of the simulation, the hydrophobic region of 6AKI is also occupied by a bubble, which, during the second half of the trajectory, becomes water-filled. The process is also shown in Supplementary Movie SM1. This event, similar to a spontaneous transition from the closed to the open state, suggests the opportunity of a deeper investigation on the structural stability of 6AKI.

Figure 5 reports the time evolution of the RMSD distance between 6AKI and the experimental structures of 4HKR and 6BBF during the equilibrium simulation of mutant P288L. In this calculation the two structures have been aligned using all backbone atoms while the actual RMSD calculation has been performed using all backbone atoms (Figure 5A), the backbone atoms of helices TM4 (Figure 5B) and the backbone atoms of the pore helices (Figure 5C). Figure 5A, shows that, during the equilibrium simulation, 6AKI becomes increasingly similar to both 4HKR and 6BBF but, unexpectedly for a putatively open structure, it is more similar to 4HKR than to 6BBF. Figure 5B shows that the increasing similarity between 6AKI and 4HKR is due to some structural rearrangement at the level of the outermost ring of helices. Indeed, as also shown in Supplementary Movie SM1, at the beginning of the simulation, TM4 helices of 6AKI are in an extended conformation, but during the run they spontaneously bend adopting a structural arrangement similar to that of the closed state. Finally, Figure 5C shows that no major structural rearrangement occurs at the level of the pore helices, even if, also in this district, 6AKI appears to be more similar to 4HKR than to 6BBF. The disruption of the bubble must thus be ascribed to some fine tuning of the pore conformation. In order to assess the variability across runs of the 6AKI-4HKR and 6AKI-6BBF RMSD, we performed a block average analysis (Frenkel and Smit, 2002) of the equilibrium simulation of 6AKI. Details can be found in the Supplementary Material section Analysis of RMSD variability where we also comment on the seemingly large RMSD values appearing in Figure 5.

The structural instability exhibited by 6AKI suggests that this structure is unlikely to represent the open state and more probably corresponds to an intermediate or a conformation on the slope of the free energy barrier which separates the closed state (4HKR) and the open state (6BBF). However, if we represent the free energy profile as a function of a single collective variable, a structural transition that brings 6AKI closer to 4HKR would automatically move it away from 6BBF and vice versa (Supplementary Figure S3A). This would be inconsistent with what observed in the RMSD analysis, i.e., that, during the simulations, 6AKI tends to get closer to both 6BBF and 4HKR. This seeming contradiction can be solved by considering that free energy is more realistically a function of several collective variables (Supplementary Figure S3B).

In order to reproduce the conformational transition from the closed to the open state and vice versa we ran Targeted Molecular Dynamics simulations (Schlitter et al., 1994). Due to the structural instability of mutant P288L (PDB: 6AKI) (Liu et al., 2019) we chose as open state the structure of mutant H206A (PDB: 6BBF) (Hou et al., 2018), while the structure of the wild type solved by Hou et al (PDB: 4HKR) (Hou et al., 2012) was considered as representative of the closed state. Specifically, we ran two TMD simulations: a 100 ns simulation from the closed to the open state and a 500 ns simulation from the open to the closed state. As can be noted, these simulations are extremely long as compared to the typical length of TMD simulations reported in the literature (1 ns or less) (Schlitter et al., 1994; Xiao et al., 2015; Liang et al., 2017; Li et al., 2013). The length of our simulations was not required to induce the structural transition but rather for the destruction or formation of the bubble in the hydrophobic region of the pore. This behaviour is not completely unexpected considering that in TMD the biasing potential is applied to the protein structure while water molecules redistribute spontaneously as a result of the changing chemical environment around them. In addition, the formation/destruction of the bubble is typically a first order transition, which, once triggered by structural changes, can be irreversible on the simulation timescale (Giacomello and Roth, 2020); therefore very slow changes in the biasing potential are needed to avoid artifact. In this regard, it is also notable that the O → C simulation had to be 5 times longer than the C → O simulation. This suggests that the energy barrier in the O → C transition is higher than that in the C → O transition implying a lower free energy of the open state as compared to the closed configuration. The water content of the pore is monitored in Figure 6. If we consider the water count in the whole pore (Figures 6A,B), we can note that the number of water molecules increases or decreases in a rather gradual way. However, if we focus on the hydrophobic region only (Figures 6C,D), we observe an abrupt transition where the number of water molecules directly jumps from zero to the final value or vice versa. This behaviour suggests that bubble formation or destruction is a rare event possibly related to the overcoming of a high energy barrier. Moreover, this highly cooperative transition was also predicted in the thermodynamic model by Roth et al. (2008) where gating occurs in a narrow window of pore radii or hydrophobicities. Finally, the fact that bubble formation and destruction occurs at the end of the simulation both in the opening and closing transition, suggests the process to be strongly hysteretic, possibly due to the presence of hidden collective variables. This issue will be the object of a future investigation using techniques successfully applied in both model nanopores (Tinti et al., 2017; Camisasca et al., 2020) and biological ion channels (Costa et al., 2021). A more detailed description of the water and ion distribution in the pore is provided in the Supplementary Material section Water and ion distribution during TMD simulations.

A long-standing controversy in CRAC channel literature concerns the gating mechanism, that, according to some authors (Yamashita et al., 2017), is due to a rotation of pore helices; however, this event was not confirmed in a number of subsequent experimental (Hou et al., 2018; Liu et al., 2019) and computational works (Frischauf et al., 2019). Section Rotation of pore helices and Supplementary Figure S14 in the Supplementary Material show that the pore helices of CRAC channel have an intrinsic tendency to rotation that is present even during equilibrium simulations; in addition, the rotations seen in the experimental structures, reported in Supplementary Table S4, show good agreement with those observed in simulations. In the following we characterise rotation occurring during biased simulations. Figures 7A,B show that rotation of pore helices does occur during the TMD simulations. However, the rotation angles measured at different levels along the pore axis suggest that rotation does not occur as a rigid body rotation but rather as a twisting. It is also noteworthy that the rotation angles of the Glu-ring and hydrophobic region are much smaller than those of the basic region, presumably because the latter is located close to the N-terminus where chain constraints are expected to be weaker.

The value of the rotation angle of the hydrophobic region averaged along the last 20 ns of the C → O simulation is ∼11 ^o while the value averaged during the last 100 ns of the O → C simulation is ∼16 ^o . These values are in reasonably good agreement with the average helix rotation of ∼15 ^o at residue F171 measured by Yamashita et al. (2017). It is important to note that, while our rotations were measured during a conformational transition from the open to the closed state and back, the rotation angles reported by Yamashita et al. (2017) are due to frequent spontaneous conformational fluctuations from the crystallographic structure that they measured during equilibrium simulations. This may be due to the fact that, even at equilibrium, the system transiently visits conformations relevant for its functional dynamics.

Finally, in the C → O transition it is possible to identify a change in slope of the rotation angles around 70 ns, just little before bubble destruction. This suggests that the rotation of pore helices might be the trigger for bubble destruction. However, the time coincidence between the change in rotation angles and bubble formation is less evident in the closing transition.

After the rotation of pore helices is ascertained, it is important to investigate the functional role of this event. According to the thermodynamic model by Roth et al. (2008), hydrophobic gating is under the control of two key parameters: the hydrophobicity of pore wall and the pore radius. We therefore checked whether either of these parameters was affected by pore helix rotation. Figures 7C,D show the time evolution of the Solvent Accessible Surface Area (SASA) in the most important regions of the pore: the Glu-ring, the hydrophobic region, and the basic region. It can be noted that, both in the opening and closing transition, the SASA plot shows a change in slope approximately at the same time of the change in slope of rotation angles. However, while the SASA change of the basic region and E-ring favour the increase of water content expected during opening and the decrease expected during closing, the SASA change of the hydrophobic region tends to oppose the wetting expected during opening and the de-wetting of the closing process. These results therefore suggest that the rotation of pore helices is not finalized to the creation of a thermodynamically favourable environment for bubble formation or destruction.

We now explore the possibility that the rotation of pore helices might cause a change of the pore radius. In Supplementary Figure S8, we consider two possible scenarios. In the first scenario, the rotation of pore helices is the main cause of the increase of the pore radius through a displacement of the bulky side chains of hydrophobic residues. In this case, if we compute a PMF as a function of rotation angle and pore radius, the low free energy region should be slanted. On the other hand, if we repeat the calculation after excluding the side chains from the calculation of pore radius, the low free energy region should be flat. In the second scenario, helix rotation and pore radius increase are simultaneous events but rotation is not the main cause of the growth of the pore radius. In this case, either performing the calculation with or without side chains, the low free energy region of the PMF should be sloped.

This strategy can be employed to analyze the PMFs derived from our TMD simulations and illustrated in Figure 8. In the opening simulations, both computing the PMF with (Figure 8A) and without side chains (Figure 8B), the minimal free energy region is inclined. This suggests that helix rotation and radius enlargement occur simultaneously but there is no causal link between the two phenomena. However, if we consider the closing simulation, we note that, when the PMF is computed using side chains (Figure 8C), the increase of the rotation angle is accompanied by a decrease in the pore radius while this effect is suppressed when removing the side chains (Figure 8D). The conclusion of these results is that rotation of pore helices might play some role in the change of pore radius. However, this driving force is probably not the only one and the change of pore radius also depends on helix displacement.

In order to reconstruct the gating mechanism, for each frame of the TMD trajectories we built a contact map (two residues are considered to be in contact if at least two heavy atoms of the side chains are closer than 5.0 Å). In this way it is possible to determine the sequence of contact breakdown (Figure 9 and Supplementary Figure S9) and formation (Figure 10 and Supplementary Figure S10). The black squares in the maps represent contacts present in the initial conformation and not broken during the simulation. The coloured symbols represent contacts formed or broken at different times during the simulation. In particular, we only considered contacts formed or broken in at least 3 of the subunits and we averaged their formation or breakdown times. In the contact maps, the elements on the diagonal represent contacts inside the same alpha helix. The elements orthogonal to the diagonal represent contacts between different alpha helices of the same subunit. Finally the islands of off-diagonal elements represent inter-subunit contacts. The maps show that all the broken and formed contacts are either intra-subunit or between neighbouring subunits n and n + 1.

During the opening simulation, the number of broken contacts is much higher than the number of formed contacts, consistent with the formation of a looser structure. As shown in Figure 9, the contact breakdown originates at the coiled coil between TM4-ext helices of neighbouring subunits and propagates towards the extracellular side through the TM3-TM4b interface. The contact breakdown then propagates towards the inner rings of helices. In particular we can observe the breakdown of contacts in the upper part of the TM2-TM3 and TM1-TM2 interface and a breakdown of contacts spanning the whole length of neighbouring TM1 helices. The few formed contacts (Figure 10A) are scattered throughout the contact map, but indeed they are very localized in space (Figure 10B). Most of them, indeed, involve the lower part of helices TM2(n)-TM4b(n + 1) and TM1(n)-TM3(n). This opening mechanism is consistent with the fact that activator STIM1 is known [Hou et al. (2012) and references therein] to break down the TM4e(n)-TM4e(n + 1) coiled-coil to form a new one. The conformational change is then expected to allosterically propagate (Yeung et al., 2018) to the innermost circle of helices where gating occurs.

The closing process, depicted in Supplementary Figure S9 and Supplementary Figure S10, is basically the reverse of the opening one. In this case the number of formed contacs is much higher than that of the broken ones. As illustrated in Supplementary Figure S10, contact formation starts in the extracellular side of helices TM3-TM4b and propagates downwards to restore the coiled coil between neighbouring TM4e helices. Contact formation also propagates towards inner helices. In both the opening and closing processes, the last events involve H206 on helix TM2. As shown in Figure 9B, during opening H206 breaks its hydrogen bond with L202 (TM2) and breaks the hydrophobic contacts with L168 (TM1) to form new hydrophobic contacts with A166 (TM1). Other broken contacts are those with F259 and G255 of TM3 (Figure 10B). The fact that the contacts involving H206 are formed and broken at the very end of the process, suggests that the whole mechanism is finalized to these events highlighting the importance of H206. This is consistent with the experimental evidence (Yeung et al., 2018) that most mutations of H206 create constitutively open mutants. Indeed our open structure was taken from the backmutated H206A mutant (Hou et al., 2018).

The contact analysis has highlighted two main events during the gating process: the massive formation or breakdown of contacts between helices TM3 and TM4b and the final formation and breakdown of several contacts involving H206. We will now analyze more in depth these two events starting with the formation/destruction of contacts in the TM3/TM4b interface.

In the C → O transition, after the massive breakdown of contacts, helices TM3 and TM4b move apart (Figure 11B) and they face a different destiny. Helix TM3, in particular its lower part, moves centrifugally. The displacement of helix TM3 makes room for a corresponding centrifugal displacement of the lower part of helix TM1: this is the opening of the basic region of the pore. These motions are described in Figure 11A where we show the displacement of helices TM1 and TM3 from the center of mass of the pore. Both helices move centrifugally, but since TM1 takes a longer step, the TM1-TM3 distance decreases (Figure 11B) justifying the formation of contacts that we observed in the contact analysis (Figure 10).

We now turn to the fate of TM4b. In Figure 11C we show the displacement of TM4b(n + 1), TM3(n), and TM2(n) with respect to the initial position of the centers of mass. The displacements have been projected on the vectors connecting the COMs of the three helices so that a positive displacement corresponds to an approach while negative displacements correspond to a distancing. Our plots show that, while TM4b(n + 1) moves closer to TM2(n) and TM3(n), also TM2(n) and TM3(n) move towards TM4b(n + 1) so that there is a decrease of distances TM4b(n + 1)-TM3(n) and TM4b(n + 1)-TM2(n) as illustrated in Figure 11D. This pattern is in agreement with the TM4b(n + 1)-TM2(n) contact formation observed in the contact analysis (Figure 10).

As shown in Supplementary Figure S11, what happens during the closing transition is the reverse of what occurs during opening. The lower parts of helices TM1 and TM3 move towards the center of the pore (Supplementary Figure S11A) and since TM1 takes a longer step, the distance between the two helices increases (Supplementary Figure S11B) and we have breakdown of TM1-TM3 contacts (Supplementary Figure S9). We also note that, as helices TM4b(n + 1), TM2(n) and TM3(n) move apart from each other (Supplementary Figure S11C), their distance increases (Supplementary Figure S11D) and we have the breakdown of TM4b(n + 1)-TM2(n) contacts highlighted in the contact analysis (Supplementary Figure S9).

The structural mechanism that underlies the distance plots in Figure 11 is depicted in Figure 12. Panels (A) and (B) of Figure 12 show a rear view of a single subunit at time zero and at 100 ns during the opening process. At time zero, in the closed state, helix TM4 has a double bend, one between TM4a and TM4b and one between TM4b and TM4-ext. In this particular arrangement it can be observed that helices TM4b and TM4-ext are behind helix TM3 that therefore is prevented from moving back. At time 100 ns, helix TM4 is fully extended. In this conformation, the space behind helix TM3 has been cleared and, since helix TM4b is no longer in the way, helix TM3 can move backwards. Panels (C) and (D) of Figure 12 show a side view of the same process. What is notable here is that the extension of TM4b allows the backward movement of the lower part of TM3. This motion, in turn, clears the space behind TM1 allowing a backward movement of the lower part of this helix that opens the basic region of the pore. This process is dynamically illustrated in Supplementary Movie S2.

As already discussed, the contact analysis has highlighted two key events. After analyzing the implications of the TM3-TM4b contact breakdown, we now turn to the events involving H206. In Figure 11 and in Supplementary Figure S11 we have observed a significant rearrangement of helices TM1, TM2, TM3. In particular, we have noted that, during the opening process, helix TM1 moves centrifugally while helices TM2 and TM3 of subunit n move towards helix TM4b of subunit n + 1. Figure 13A shows that these events cause an increase in the distances TM1-TM2, TM1-TM3, and TM2-TM3. Since the side-chain of H206 leans in the space between these three helices, their distancing causes a decrease in the atom density around H206, see Figure 13B showing the number of atom contacts of H206 as defined by a distance cutoff of 7.5 Å. Finally, Figure 13C shows that as soon as the atom density around H206 decreases, the side chain of this residue swings due to a rotation around dihedral angle χ ₁. Supplementary Figure S12 shows that the reverse process occurs during closing. In this case there is a decrease of the distances, TM1-TM2, TM1-TM3, and TM2-TM3 (Supplementary Figure S12A). This causes an increase in the atom density around H206 (Supplementary Figure S12B) so that the side chain of this residue is pushed back (Supplementary Figure S12C) to the conformation typical of the closed state.

Figure 14 provides a structural view of the conformational transition involving H206. For the sake of clarity only a short stretch of helices TM1 and TM2 is shown. Figures 14A,B offer a view from the front while Figures 14C,D from the rear. In the closed state, as predicted by Yeung et al. (2018), H206 is hydrogen bonded with L202 also located on TM2. H206 is also close to S165, but differently from Frischauf et al. (2019), no hydrogen bond can be seen. Finally, in the closed state, H206 forms hydrophobic contacts with L168 on helix TM1. When opening occurs, H206 rotates around the χ ₁ dihedral breaking the hydrogen bond with L202 and swinging from the contact with L168 to a hydrophobic contact with A166 (also located on TM1). It is important to note that this rearrangement of H206 creates a free space behind the hydrophobic region of TM1. A dynamic view of this process is provided by Supplementary Movie S3.

The conformation of the closed state of CRAC channel was taken from the Protein Data Bank (PDB ID: 4HKR) and the missing loops 181–190 and 220–235 were modelled with the MODELLER 9.21 software (Sali and Blundell, 1993). This structure contains mutations C224S, C283T introduced to prevent non specific disulfide formation and P276R, P277R (in the hypervariable TM3-TM4 loop) introduced to produce well ordered crystals (Hou et al., 2012). In order to better compare our computational results with the experimental data in the literature, these mutations were retained. The structures of the constitutively open mutants H206A and P288L were also taken from the Protein Data Bank (PDB ID: 6BBF and 6AKI respectively). In order to prepare the systems for Targeted Molecular Dynamics simulations (Schlitter et al., 1994) that require the same number of atoms in the initial and final conformations, short N-terminal and C-terminal extensions have been added to both 6BBF (residues 144–147 and 328–334) and 6AKI (residues 144–149 and 331–334) to match the sequence length of 4HKR, using the MODELLER software (Sali and Blundell, 1993). The structures of 6BBF and 6AKI have been back-mutated to the wild type sequence (mutation A206H in 6BBF and L288P in 6AKI). In order to match the sequence of 4HKR mutations P276R and P277R have also been introduced. Protonation states have been assigned using the WHATIF server (Vriend, 1990). A number of open issues pertaining to the experimental structures of the closed state (PDB: 4HKR) and the putatively open states (PDB: 6BBF, 6AKI) is discussed in Supplementary Material section Open issues on CRAC experimental structures. The three channels were embedded in a lipid bilayer comprising 533 molecules of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) and solvated by 62,927 water molecules using the CHARMM Membrane Builder (Wu et al., 2014). We added 184 sodium ions and 172 chloride ions so as to neutralize the charge of the channel and reach a final concentration (NaCl) = 0.15 M. Overall, the system comprised 279,195 atoms and the simulation box had dimensions 145 × 145 × 130 Å.

All simulations were performed with the NAMD 2.11b2 suite of programs (Phillips et al., 2005) using the ff15ipq force field (Debiec et al., 2016) for the protein, the Lipid17 force-field (Dickson et al., 2014) for the phospholipids and the SPC/E water model (Takemura and Kitao, 2012).

The three systems first underwent 10,000 steps of conjugate gradient minimization. During equilibration harmonic restraints were applied to nonhydrogen atoms of the protein backbone and side-chains as well as to the phospholipid heads. A harmonic restraint was also applied to the dihedral angle formed by carbons 8–11 of oleoyl acid and to the improper dihedral C₁-C₃-C₂-O₂ involving the three carbons of the glycerol unit and the hydroxyl oxygen linked to its central carbon. The equilibration was organized in 12 stages whereby the constraints were gradually released. The values of the force constants used in the 12 stages can be found in Supplementary Table S1.

The production run, following the restrained equilibration, was carried out in the isothermal isobaric ensemble for 100 ns. Pressure was kept at 1 atm by the Nosé-Hoover Langevin piston method, and the temperature was kept at 300 K by coupling to a Langevin thermostat with a damping coefficient of 1 ps⁻¹. Long-range electrostatic interactions were evaluated with the smooth particle mesh Ewald algorithm. For the short-range nonbonded interactions, we used a cutoff of 12 Å with a switching function at 10.0 Å. The integration time step was 2 fs, and the bonds between hydrogen and heavy atoms were fixed to eliminate the most rapid oscillatory motions.

Pore radius profiles have been computed using the HOLE program (Smart et al., 1996). The Potential of Mean Force (PMF) of water and ions as a function of the axial position and of the distance from the pore axis was computed using equation (Im and Roux, 2002): PMF(z, r) = − k _B T log(ρ(z, r)/ρ _bulk ) where k _B is Boltzmann constant, T is the absolute temperature, ρ _bulk is the water or ion density in the bulk and ρ(z, r) is the density in position (z, r). In order to avoid a divergence of the PMF in ranges of z devoid of samples, we arbitrarily assign to these bins the maximal value of the free energy of the profile, i.e. the value of the free energy corresponding to the bin with the smallest non-zero water density. The assignment of this constant value to all the bins devoid of samples leads to the emergence of an artificially flat region near the maximum of the PMF as can be observed, e.g,. in Figure 3C; for this region the physical expectation is to find water molecules in the gas state. The electrostatic potential inside the pore has been computed using the PMEPot plugin (Aksimentiev and Schulten, 2005) of VMD as detailed in the Supplementary Methods section of the Supplementary Information. This section also provides information on our protocols to compute rotation angles, solvent accessible area and hydrophobicity profile.

Targeted Molecular Dynamics (TMD) (Schlitter et al., 1994) is an enhanced sampling technique used to drive a system between given end states. The algorithm applies a harmonic biasing potential to the RMSD between the current conformation and the target conformation: U T M D = 1 2 k N [ R M S D ( t ) − ρ ( t ) ] 2

where k is the force constant, N is the number of atoms where the biasing force will be applied and ρ(t) is the target RMSD that during the simulation is linearly decreased from the RMSD between initial and final structure to zero. In our case the system was steered from the final frame of the equilibrium simulation of the closed state (PDB: 4HKR) to the final frame of the equilibrium simulation of the open state (PDB: 6BBF) and vice versa. We applied a biasing force constant of 100,000 kcal/mol/Ų and we chose a simulation lenght of 100 ns for the C → O transition and 500 ns for the O → C transition. Since the initial distance between closed and open state is 13.94 Å, this choice ensured a sufficiently slow transition with a variation of the RMSD of ∼0.14 Å/ns in the C → O and of ∼0.03 Å/ns in the O → C transition. This extremely slow protocol should allow the relaxation of the degrees of freedom that were not biased during the simulation. The good superposition between the time-evolutions of the current and target RMSD (Supplementary Figure S13) is an indication that the TMD was run slowly enough to allow the system to follow the target RMSD schedule. Finally, in order to avoid unfolding events an harmonic constraint was applied on the α-helical content of TM4-ext helices (force constant 50,000 kcal/mol/Ų).