Unraveling of a Strongly Correlated Dynamical Network of Residues Controlling the Permeation of Potassium in KcsA Ion Channel

The complicated patterns of the single-channel currents in potassium ion channel KcsA are governed by the structural variability of the selectivity filter. A comparative analysis of the dynamics of the wild type KcsA channel and several of its mutants showing different conducting patterns was performed. A strongly correlated dynamical network of interacting residues is found to play a key role in regulating the state of the wild type channel. The network is centered on the aspartate D80 which plays the role of a hub by strong interacting via hydrogen bonds with residues E71, R64, R89, and W67. Residue D80 also affects the selectivity filter via its backbones. This network further compromises ions and water molecules located inside the channel that results in the mutual influence: the permeation depends on the configuration of residues in the network, and the dynamics of network’s residues depends on locations of ions and water molecules inside the selectivity filter. Some features of the network provide a further understanding of experimental results describing the KcsA activity. In particular, the necessity of anionic lipids to be present for functioning the channel is explained by the interaction between the lipids and the arginine residues R64 and R89 that prevents destabilizing the structure of the selectivity filter.

The simulations were performed using NAMD 2.8 and NAMD 2.9 [1] in the NpT ensemble with pressure 1.01 bar (Langevin piston implemented via Nosé-Hoover method [1]), temperature 310 K (Langevin thermostat [1]) and electrostatic forces were computed using the smooth particle-mesh Ewald method (SPME). [2] Metadynamics calculations were performed using the built-in colvar module. [3] A multiple-timestep algorithm was used [4,5] and different integration steps according to the scope of simulations: • 1 fs for unbiased simulations, bonded interactions computed every time step, non-bonded non-electrostatic forces computed every 2 time steps and electrostatic forces every 4 time steps; • 2 fs for biased simulations, bonded and non-bonded non-electrostatic forces computed every time step and electrostatic forces every 3 time steps; Using longer integration steps (2 fs) is a common choice for biased simulations, since it allows longer sampling. The CHARMM27 force field (FF) was used for the protein, with a modification in the Lennard-Jones term to represent the interaction between K + and the carbonyl oxygens of the protein, CHARMM36 for the lipids and TIP3P for water. [6][7][8] The system was prepared by embedding the X-ray structure (pdb code 1K4C, 6284 atoms; solved at 2 Å resolution [9]) into a membrane patch of 1-palmitoyl-2-oleoylphosphatidylcholine (POPC, 222 molecules and 29748 atoms) prepared using the VMD 1.8.6 membrane plug-in, modified in order to work with CHARMM36. The protein was solvated using the Solvate program by Grubmüller [10], and the VMD solvate plug-in (17740 water molecules, 53220 atoms). MD free energy sampling, and QM/MM calculations (quantum mechanical/molecular mechanical) suggest that, in the X-ray conformation, the E71/D80 pair is linked by a hydrogen bond (H-bond) where one residue is ionized, and the proton exchanges with a preference for E71. [11,12] Consistent with conventional MD force fields, this has been mimicked in this work by assigning the proton to E71 in every protein where the E71-D80 bridge is present (WT, R64A, Y82A, L81A). Water molecules inside the cavity that were not resolved in the X-ray structure were inserted using the DOWSER program. [13] The missing parameters for K + necessary for DOWSER were taken from the GROMOS87 FF for consistency. [14] The potassium concentration in the aqueous phase was 0.2 M, equal to the bulk concentration used in many experiments: [9,[15][16][17][18][19] 75 Cl − ions and 63 K + ions, with 2 K + in the selectivity filter and 1 K + in the cavity. The final system, showed in Fig. S1, consisted of 89390 atoms and a size of 95 × 95 × 100 Å.
Mutants were created in-silico from the crystal structure of the WT which allow more accurate comparisons, following a commonly-used approach for studying mutants features. [15,16,20] Relaxation of the WT and preparation of the mutants were performed in different steps before the production runs. (i) A 2 ns MD simulation in which the protein, water, ions and lipid head groups were frozen; this ensured melting of the lipid tails. (ii) A further 2 ns MD simulation in which the protein (only) was restrained with harmonic tethers. After this step, the membrane having relaxed around the protein, the in-silico mutants were created using the program psfgen [1] based on the CHARMM topologies. [7,21] A further 1 ns MD simulation in which all atoms were unrestrained was performed for each protein.

Well-tempered Metadynamics
Metadynamics (metaD) was originally developed as a flexible tool to explore new reaction pathways and accelerate the observation of rare events, as well as to reconstruct free energies. [22] It is based on a history-dependent potential added to the Hamiltonian of the system which biases the simulation toward sampling unexplored configurations. The bias potential is created as a sum of Gaussians centered along the trajectories projected onto a subset of degrees of freedom, usually referred to as collective variables (CVs). In the case of a comprehensive sampling, the detailed knowledge of the bias potential allows the unbiased Boltzmann distribution along the CVs to be retrieved, and thus allows calculation of the free energy hyper-surface. [23,24] In the well-tempered metadynamics (wt-metaD), the height of a Gaussians added is history-dependent, and this dependence is associated with a parameter having the dimension of a temperature ∆T.
Dynamics of L81/R64 were investigated using ∆T = 1200 K, the E71-D80 distance using ∆T = 1500 K. It is worth noting that, although important from a practical perspective, the computed FES will eventually converge to the real one irrespective of the particular choice of ∆T, for proper sampling. [24] The initial deposition rate ω = w 0 /τ G of the Gaussians affects relaxation over the orthogonal degrees of freedom. [24] The parameter w 0 appeared to be very important for the success of calculations and was chosen to be small (0.005 kcal/mol) in order to produce a smoother (quasi-adiabatic) filling of the biasing potential over the whole calculation, which is desirable to avoid the occurrence of non-physical events. The deposition time, τ G , was selected to be 200 fs and a few preliminary tests suggested that the final result is not sensitive to the choice of this parameter in the cases under investigation.

Permeation in the V76 flipped configurations of E71A and WT
The flipping of V76 has been reported in several published works [25][26][27][28][29] and it has been suggested as being able to generate non-conductive conformations responsible for the inactivation of the filter (C-type inactivation) or for the modal-gating (flicker mode). [26,29] Any conformation in which the G77-V76 peptide group is rotated with respect to the putative conductive X-ray structure has been referred to as "V76 flipped". The conductivity of the V76 flipped conformation for WT and E71A was tested by means of MD simulations to confirm this hypothesis.
Simulations were initiated from the last configurations obtained from simulations of WT and E71A in which V76 flipped appeared as meta-stable states with lifetime longer than 10 ns. Outward transitions of K + were induced by interchanging the positions of a K + ion from the bulk and a water molecule in the cavity. In this way, because of the closed state of inner gate which results in the cavity to be closed, the electrostatic repulsion of the two ions in the cavity led to promotion of the multi-ionic movements through the selectivity filter, associated with the permeation. Simulations were preceded by 10 5 steps of minimization (conjugate gradient and line search algorithm as implemented in NAMD [1]). Two simulations were performed for the WT with different initial velocities and one for E71A.
Both E71A and WT proved to be conductive and the flipped state of V76 to be reversible within a few ns when ion permeation occurs. It is likely that the creation of V76 flipped states, with longer lifetimes necessary to hinder K + current, with the longer lifetime being dependent on conformational changes in residues in the region behind the selectivity filter. Occurrence of a vacancy in the filter were observed during the permeation. . Notably, the value SC 80 ≈ 13.5 Å, which was very close to the value obtained for X-ray structure of the WT (≈ 13.3 Å), was found to be a metastable state for E71A, though completely stable (during the simulation) for the WT. This suggests the conformation of D80 found in the crystal structure is stabilized by the E71-D80 H-bond.

Time series for the transitions of E71A
Time series of the order parameters are reported in Fig S3 for the first 10 ns of the simulation (total length 24 ns) of E71A mutant. The order parameters are the following: i) the V76 ψ dihedral angle ψ 76 , ii) the position SC 80 of the D80 side chain (the distance between D80 C γ and A73 C α atoms); iii) the distance D80-R89 between D80 (C γ atom) and the nearest R89 side chain (C ζ ) iv) the ψ dihedral angle of D80 ψ 80 ; v) the z coordinates of the outermost ion (K1) bound to the selectivity filter, centered relative to the COM of the selectivity filter. Rapid fluctuations in the instantaneous values of the parameters obscured trends in some cases, especially for the D80 position SC 80 . So, only the smoothed trajectories are shown (moving averages procedure were calculated with the time window interval of 0.25 ns). The part of the simulation before the flipping of V76 (the first 1 ns have been omitted as relaxation) was subdivided into three sections : i) before the creation of D80-R89 H-bond, up to 2.5 ns (blue shadow area); ii) during the creation of D80-R89 H-bond, between 2.5 and 5.9 ns; and iii) after the creation of D80-R89 H-bond and before flipping of V76, between 5.9 and 9 ns (green shadow area). The time series show that the strengthening of the H-bond (D80-R89 distance reduces) is linked to D80 moving upward (SC 80 increases). Note that the change in the position SC 80 is more visible in Fig. S4 (discussed below). The creation of the H-bond is accompanied by the partial flipping of V76 (6.8 ns, a rotation by 50 degrees of ψ 76 ), soon followed (7 ns) by an inward movement of K1 (1.9 Å) to the position between sites S1 and S2, accompanied by an adjustment of the water molecule in S2 (not shown). At 9 ns the complete flipping of V76 was observed (red vertical line on the graphs in Fig. S3). Substituting the oxygen of the V76 backbone, this water molecule in site S2 re-established the square anti-prism geometry of the coordination sphere of K1, which then moved back to its original position in site S1. The final transition involved the D80 side chain (9.5 ns, ψ 80 changes from −43 to −22 • ) that  broke the H-bond and restored the initial uncorrelated motions of the D80 and R89 side chains (see Figure 2 in the main text). Fig. S4 shows the change in position of D80 following the creation of D80-R64 H-bond in the region before (blue) and after (light green). The mean value of the D80 position changes from 13.5 to 13.8 Å and two-sample Kolmogorov-Smirnov test indicates this difference is statistically significant with p-value < 2.2 · 10 −16 .

Comparison of the dynamics of residues
Fig. S5 reports the root mean square displacement (RMSD) of residues with respect to the reference X-ray crystal structure of WT. RMSDs computed for the simulations of WT, Y82A and R64A proteins. RMSD was estimated as time averaged deviation of the backbone atoms of a residue from the reference structure. Only the region surrounding the selectivity filter is considered (residues from 63 to 90). Note, the error bars in the figure are standard deviations of the RMSDs. Fig. S6 shows an example of a partial rotation of the D80 side chain around the χ1 dihedral angle. The rotation was observed several times during the simulation of Y82A mutant.

Simulation WT-R64D80
In the present section, additional material is reported for the illustration of changes in the pore region during the simulation commenced with residue R64 near D80 in all the four subunits (simulation WT-R64D80). Fig. S7 illustrate conformational changes of the protein during the first part of the simulation. Fig. S8 shows snapshots of simulation WT-R64D80 when H-bonds between D80 and R64, or R89 or both R64 and R89 residues were formed. The probability density of z-coordinate of ions K1 and K2 are shown in Fig. S9 for simulations starting with WT-R64D80 and WT X-ray configurations.

Details of the FES calculation of R64/L81 motions
In the main text, the well-tempered metadynamics calculation of energetics of R64/L81 motions is presented. The two-dimensional FES was computed as functions of the CVs: i) the distance of the R64 side chain (C ζ atom) from the COM of the selectivity filter (R64-SF); ii) the χ1 dihedral angle of the L81 residue (χ1 81 ). The calculation was obtained by considering only the configuration KwK0K of the K + in the pore, with a vacancy in S3. This configuration remained stable for the whole calculation and the FES shows a quick convergence.

Free energies for R89 motions
The FESs (Fig. S10) of the interaction between R69 and D80 was computed as a function of the distance, SF − R89, between the C ζ atom of the R89 guanidino group and the COM of the selectivity filter. The FES was obtained on the base of two different unbiased simulations starting with the WT X-ray structure and using data from all the four subunits. The total sampling interval is 224 ns.
Initial estimations of the FES revealed that the FES depends on a configuration of ions in the selectivity filter. This observation confirms the existence of the strong link between the position of the K + ions in the selectivity filter and the state of the surrounding region. From all possible, three different configurations were typical and, therefore, the FES was calculated for them. These configurations are S2/S4, S1/S3, and S0/S2/S4, and the corresponding FESs are shown in Fig. S10.
The FESs confirm that in a similar manner to R64, the side chain of R89 can create a H-bond with D80. Two minimums corresponding to the presence and absence of the D80-R89 H-bond are clearly recognizable in the FESs, respectively at SF − R89 ≈ 15.7 Å and SF − R89 ≈ 17.1 Å. The barrier for the bond formation is defined by the transition from the state with SF − R89 ≈ 17.1 Å to the state with SF − R89 ≈ 15.7 Å. This barrier depends on the ions configuration and it is relatively low between 0.25 and 0.7 kcal/mol. The opposite barrier for breaking the D80-R89 H-bond depends on the ion configurations as well. The barriers are comparable for configuration S0/S2/S4, but the bond breaking is less probable when the outermost sites are not occupied. Thus, the R89-D80 H-bond is more stable when an ion is not in S0. Such a configuration is observed in the case of a low concentration of K + in the bulk, which is the condition when the inactivation probability is large.
The dependence of the FESs on the ion occupancy further confirms the role of K + ions bound to the filter in the dynamics of the pore region. However, the low barriers between the states in the FESs also indicate that the destabilization effect of R89 on E71-D80-W67 linkages is weaker in comparison to  those delivered by R64 alone. So, R69 plays a facilitating role in the joint action with R64 for breaking E71-D80-W67 linkages.

Radial distribution functions
Fig. S11 reports the radial distribution functions (RDF) between Cl − ions and residues R64 and R89. The RDFs confirm the strong interaction between the arginines and negatively charged elements of the extracellular bulk. The functions were estimated on the base of an unbiased simulation of WT. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q
The various conformations of the filter content reported in the wt-metaD study, were selected based on extensive preliminary unbiased simulations. The simulations revealed that states with K + ions in S2/S4 and no ion in the cavity are more stable than others. In these states K + ions from the bulk can enter and leave site S0 on scales of tens of ns. Configurations with K + ions bound to S1/S2 had the shortest residence times. If an ion is present in the cavity then the configuration S1/S2 underwent transitions to S0/S2/S3 within 1 ns. In the absence of ion in the cavity, within few ps the configuration S1/S2 first is changed to S1/S3 followed by the transition to S2/S4.
Since configurations with ions in sites S2 and S4 are typically observed in the unbiased simulations, such ions occupancy was chosen in order to avoid additional restraints on ions that might affect the pore dynamics. A set of different K + configurations was selected for the free energy calculations according to the K + occupancy of S0 (S0/S2/S4 and S2/S4), and further divided depending on the presence of a vacancy in S3: KwKwK and wwKwK (no vacancies); KwK0K and wwK0K (a vacancy in S3). Note that the permeation mechanism which include vacancies have been reported. [30][31][32] The total 8 free energy calculations were performed, and the estimated FESs are reported in Fig. S12. The case with a vacancy in site S3 is discussed in the main text. Comparison between calculations performed with or without a vacancy in S3 confirms the similarity of the FESs and, therefore, the conclusions drawn in the main text. The differences in the FESs should be considered as being dependent on the flipping of V76, which is promoted by the presence of a water molecule in S2 or S3 [26]. The latter confirms the existence of a link between D80-region and V76 rearrangements. R64 close, no ion in S0 R64 close, K + in S0 Figure S12. The FESs were obtained by the wt-metaD calculation for to the distance, D80 − E71, between C γ atom of D80 and the H-bond donor oxygen atom of E71. The FES and distance D80 − E71 are shown in kcal/mol and in Å respectively.