Conduction and Gating Properties of the TRAAK Channel from Molecular Dynamics Simulations with Different Force Fields.

In recent years, the K2P family of potassium channels has been the subject of intense research activity. Owing to the complex function and regulation of this family of ion channels, it is common practice to complement experimental findings with the atomistic description provided by computational approaches such as molecular dynamics (MD) simulations, especially, in light of the unprecedented timescales accessible at present. However, despite recent substantial improvements, the accuracy of MD simulations is still undermined by the intrinsic limitations of force fields. Here, we systematically assessed the performance of the most popular force fields employed to study ion channels at timescales that are orders of magnitude greater than the ones accessible when these energy functions were first developed. Using 32 μs of trajectories, we investigated the dynamics of a member of the K2P ion channel family, the TRAAK channel, using two established force fields in simulations of biological systems: AMBER and CHARMM. We found that while results are comparable on the nanosecond timescales, significant inconsistencies arise at microsecond timescales.

Introduction

Potassium channels are integral membrane proteins responsible for the efficient and selective conduction of K+ ions across the membrane. Depending on their overall architecture and functionality, three main classes can be distinguished: voltage-gated (Kv), inward rectifier (Kir), and two-pore-domain (K2P) K+-channels.[1] Even though the Kv and Kir are the most studied families, a growing interest in K2P channels has recently arisen.[2] Unlike the vast majority of homotetrameric K+-channels, members of the K2P family are homodimers but they are nonetheless arranged in a 4-fold-like symmetry around the pore axis of the protein (see Figure a).[3] Another of their peculiar features is the presence of an extracellular domain, known as the “cap”, whose function and specific role in ion conduction, if any, is yet to be established.[4] At the extracellular entrance of the pore, each subunit carries two P-loop regions (P1 and P2) with distinct signature sequence motifs. Like in canonical K+-channels, these are arranged to form a narrow constriction, the selectivity filter (SF), where a series of binding sites for potassium ions are hosted (S0 to S4 starting from the extracellular side, see Figure b).[5] Each site consists of eight coordinating oxygen atoms of the backbone, except for S4, which is made by backbone and sidechain oxygen atoms of threonine residues, and S0, which contributes only with four coordinating oxygens from the protein.[5,6] K2P channels are responsible for setting the resting potential of the membrane close to the electrochemical equilibrium for K+ ions and counterbalancing membrane depolarization.[7] Their activity is influenced by several physicochemical stimuli, including temperature, pH, membrane stretch, and even changes in the membrane potential, just to mention a few.[7,8] This multifactorial regulation is peculiar to this family of channels and its complexity has only recently started to be fully acknowledged. From this standpoint, molecular dynamics (MD) simulations and related methods are becoming essential tools to complement experimental observations and to characterize the structural and molecular mechanisms implicated in the functionality of K2P channels.[9−14]

Figure 1

(a) Representative snapshot of the simulation system where the protein is shown in cartoon representation, and the lipid molecules and the potassium ions are represented with van der Waals spheres. The two subunits of the channel are displayed in white and gray colors. Water molecules and bulk ions are omitted for the sake of clarity. (b) Close-up view of the SF in the two initial configurations considered in this study (only three out of the four chains are shown). The residue name of the amino acids that contribute to the SF from the two P-loop regions of each subunit (P1 and P2) are also shown. Ion binding sites are labeled S0 to S4.

MD simulations have traditionally played a key role in the field of ion channel research as aspects related to conduction and selectivity are governed by detailed atomic-level features that cannot be directly observed by experiments.[15] Since the pioneering studies on the prototypical KcsA K+-channel, MD simulations have been widely used for this purpose,[16−24] and they have proved instrumental in portraying two alternative mechanisms of ion conduction differing by the presence[17,19,20] or absence[25,26] of interposed water molecules (KWK and KK mechanisms, respectively, see Figure b). Despite their popularity, MD simulation methods have historically suffered from two main drawbacks. The first is related to the time-consuming nature of the calculations involved which limits the observation of events to the timescales that can be accessed through the available computational resources. Owing to recent software and hardware improvements,[27] together with advances in methods and protocols suited to enhance sampling efficiency,[28,29] multi-microsecond-long simulations can be routinely achieved in high-performance computational facilities at present. To date, MD simulations have characterized at a fully atomistic level several functional aspects of K+-channels, including inactivation and gating processes.[9,12,30−32] The second drawback of MD simulations concerns their reliability, which is dictated by both the quality of the underlying potential energy function and the level of realism used to describe the model system under investigation. The latter issue is related to the choice of fixed protonation states for titrable residues (e.g., Asp, Glu, His) as well as oversimplified descriptions of lipid composition for membrane systems. Concerning force fields, while they have reached a remarkable maturity in the field of biomolecular simulations,[33] they are still severely challenged by the complexity of specific interactions like those taking place in the SF of ion channels. This is particularly relevant when it comes to simulating ion conduction where the detailed balance of attractive and repulsive electrostatic interactions, together with the lack of explicit treatment of polarization effects, push empirical force fields to their limits. Furthermore, state-of-the-art simulations, those reaching microsecond timescales, disclose artifacts or suboptimal features of classical force fields that were previously unnoticed simply by virtue of the shorter timescales attainable.

AMBER and CHARMM are two distinct families of force fields that have gained great popularity in the field of MD simulations of ion channels. While they provide similar results when simulations in the nanosecond timescales are considered, significant differences emerge in the microsecond regime. In a recent study, we have demonstrated that the structure of the SF in the KcsA K+-channel is particularly sensitive to the choice of the force field, displaying remarkable stability with AMBER and collapsing in simulations performed with CHARMM.[34] To understand if such a behavior is specific to KcsA or if, in contrast, it reflects a more general trend, here, we extended the comparison of the two force fields to the simulations of the TRAAK channel.

TRAAK (TWIK-related arachidonic acid activated K+-channel, KNCK4, K2P4)[35] is a member of the TREK subfamily of K2P channels that have recently been the subject of intense experimental and computational investigation.[11,12,36,37] We opted to study the dynamics and ion conduction of this channel since the SF bears significant differences compared to the sequence of canonical K+-channels (T75V76G77Y78G79 in KcsA). In detail, these include: (i) different sequences in P1 and P2, (ii) V76 is replaced by isoleucine in P1, and (iii) Y78 is replaced by phenylalanine in P2. In addition, TRAAK is a mechanosensitive channel gated by membrane stretch, the positive curvature of the lipid bilayer, and polyunsaturated fatty acids.[38,39] Gating upon membrane stretch is a characteristic shared by other members of the TREK subfamily of K2P channels. X-ray structures show that the TM4 helices of these channels can bend around the position of a conserved glycine (Gly268 in TRAAK) giving rise to two distinct conformational states, most likely involved in mechanosensing. These states are referred to as “up” and “down” depending on the position adopted by the TM4 helix of one subunit in comparison to the TM2 of the other chain of the channel. In the up state, TM4 and TM2 interact, while in the down state, TM4 is significantly displaced from the other subunit, resulting in a lateral fenestration that opens on the TM region of the channel directly facing the lipid environment. MD simulations have been used to investigate the gating mechanism of this family of channels, and it is therefore important to perform a comparative study of the behavior of the transmembrane segments directly involved in the gating process with the two force fields. Analysis of microsecond-long trajectories of ion conduction carried out with starting configurations matching both the KWK and KK mechanisms, for a total simulation time of 32 μs, showed significant differences between AMBER and CHARMM simulations. Taken as a whole, these results confirm the need for reconsidering the performances of widespread empirical force fields in light of currently achievable timescales.

Methods

System Setup and MD Simulations

The model of the TRAAK channel was based on the Protein Data Bank entry 4WFE.[38] The entire transmembrane domain of the channel was considered from residue Arg2 to residue Arg284. The protonation and tautomeric state of residues at pH 7.0 was assigned with PROPKA.[40] The channel was inserted into a homogeneous 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane using the CHARMM-GUI webserver[41,42] according to the information stored on the Orientations of Proteins in Membranes (OPM) database.[43] The system was solvated using TIP3P water molecules,[44] and a 200 mM concentration of KCl was added to the simulation box. Potassium ions were manually placed at the SF sites S0, S2, and S4. For simulations starting from the KWK configuration, water molecules were also added at sites S1 and S3. Two force field models belonging to the CHARMM and AMBER families were used: CHARMM36m[45] and ff14SB,[46] respectively. Thus, based on the initial configuration of the SF, two subsets of MD simulations were conceived for each force field, hereafter referred to as AMBER-KK, AMBER-KWK, CHARMM-KK, and CHARMM-KWK.

For the CHARMM set of simulations, the standard Lennard-Jones parameters included in the CHARMM36m distribution were used[45] and van der Waals interactions were switched off between 10 and 12 Å. For the AMBER set of simulations, the Joung and Cheatham parameters for ions were used,[47] and van der Waals interactions were truncated at 9 Å together with the standard AMBER scaling for 1–4 interactions. Note that the different setups reflect the most commonly used combinations of force fields and parameters adopted by the community, that is the CHARMM force field parameters for protein, lipids, and ions and AMBER for proteins and lipids together with Joung and Cheatham parameters for ions. This choice leads to the use of a different cutoff for nonbonded interactions in both sets of simulations in consistence with those employed in the development of the corresponding force field. The same setups have recently been used for comparing the AMBER and CHARMM force fields in the study of the KcsA channel.[34] In all the simulations, long-ranged electrostatic interactions were treated with the particle mesh Ewald method using a grid spacing of 1.0 Å,[48] and the SETTLE algorithm was used to restrain bonds involving hydrogen atoms to use a 2 fs timestep.[49] The standard CHARMM-GUI protocol whereby the system is heated up and any initial restraints are smoothly released was adopted. The temperature was controlled using a Langevin thermostat with a damping coefficient of 1 ps–1, while the pressure was maintained at the constant value of 1 atm by coupling the systems to a Nosé-Hoover Langevin piston, with a damping timescale of 25 ps and a period of 50 ps.[50] The membrane potential was simulated by applying a constant electric field acting along the whole length (L) of the z-axis of the simulation box. In particular, an electric field, E, corresponding to two different transmembrane potentials, ΔV, of +100 and +200 mV, was applied to all of the atoms of the system[51,52]Production runs were performed in the NVT ensemble at the target temperature of 310 K as previously described.[52,53] The NAMD-2.12 MD engine was used for all of the simulations with a timestep of 2 fs.[54] For each set of simulations, two copies were performed at equilibrium conditions and three copies were run under +100 and +200 mV. Every run lasted 1 μs of sampling, amounting to a total simulation time of 32 μs. Data were aggregated unless otherwise stated. A summary of the simulations considered in this study is presented in Table .

Table 1

Summary of the Simulations Considered in This Study and the Relevant Features of the Different Simulation Runs

force field	initial SF configuration of ions	membrane potential (mV)	replicas	prevalent conformational statea	lipids in cavityb	conduction events (#)
AMBER	KK	0	#1	up/up	no	0
AMBER	KK	0	#2	up/up	no	0
conductance: –
AMBER	KK	100	#1	up/up	no	2
AMBER	KK	100	#2	up/up	no	4
AMBER	KK	100	#3	up/up	no	2
conductance: 4.3 ± 1.9 pS
AMBER	KK	200	#1	up/up	no	3
AMBER	KK	200	#2	up/up	no	11
AMBER	KK	200	#3	up/up	no	8
conductance: 5.9 ± 3.2 pS
AMBER	KWK	0	#1	up/up	no	0
AMBER	KWK	0	#2	up/up	no	0
conductance: –
AMBER	KWK	100	#1	up/up	no	1c
AMBER	KWK	100	#2	up/up	no	0
AMBER	KWK	100	#3	up/up	no	2c
conductance: 1.6 ± 1.6 pS
AMBER	KWK	200	#1	up/up	no	0
AMBER	KWK	200	#2	up/up	no	0
AMBER	KWK	200	#3	up/down	yes	0
conductance: –
CHARMM	KK	0	#1	up/down	yes	0
CHARMM	KK	0	#2	up/down	yes	0
conductance: –
CHARMM	KK	100	#1	up/up	no	0
CHARMM	KK	100	#2	down/up	yes	0
CHARMM	KK	100	#3	up/up	no	0
conductance: –
CHARMM	KK	200	#1	up/up	no	0
CHARMM	KK	200	#2	down/up	yes	11
CHARMM	KK	200	#3	up/up	no	0
conductance: 2.9 ± 5.1 pS
CHARMM	KWK	0	#1	up/down	yes	0
CHARMM	KWK	0	#2	up/down	yes	0
conductance: –
CHARMM	KWK	100	#1	up/up	no	0
CHARMM	KWK	100	#2	up/up	no	0
CHARMM	KWK	100	#3	up/up	no	0
conductance: –
CHARMM	KWK	200	#1	down/up	yes	0
CHARMM	KWK	200	#2	down/up	yes	0
CHARMM	KWK	200	#3	up/down	yes	0
conductance: –
Total Simulation Time: 32 μs

The duration of each simulation was 1 μs amounting to 32 μs. Experimental conductance: 65.4 ± 3.2 pS at +100 mV under symmetrical 150 mM KCl.[70]

The distance threshold between up and down states is considered to be 8 Å. A conformational state is regarded as prevalent in a single run if it persists by more than 50% of the total simulation time.

A lipid is considered to be inside the cavity if the headgroup is found at a distance lower than 10 Å in the channel and for more than 50% of the total simulation time.

The conduction events were sampled after the exit of water molecules from the selectivity filter.

Analysis of Trajectories

Analysis of the SF stability, backbone dihedral angles, and permeation events was performed with the python libraries available from the MDAnalysis toolkit.[55,56] The distance of lipid atoms from the channel axis was computed with the driver functionality of PLUMED-2.4.[28] The ggplot2[57] package running under the R 4.0.2 programming environment[58] was used to carry out statistical analysis on the distribution of atomic distances and to plot results.

Two dimensionality reduction methods were used to visualize and compare the conformation of the SF in the different set of simulations: principal component analysis (PCA) and sketch-map.[59,60] PCA was performed using the dimRed package running under the R 4.0.2 programming environment.[57] Specifically, the values of the backbone dihedral angles (φ,ψ) of all of the residues of the SF and the χ1 angles of Thr129 and Thr238 were converted into sin- and cos-transformed variables prior to carrying out PCA.[61,62] Conversely, the dihedral angle values of the same set of residues were directly used for the sketch-map analysis. Sketch-map is a nonlinear multidimensional scaling (MDS) method that seeks to preserve middle-ranged distances of the input data structure by minimizing the following stress function[59]where R and r are the distances between two points evaluated in the high- and low-dimensional (low-d) space, respectively, while F and f are two sigmoid functions of the formandThe effect of the sigmoid functions is to shrink and expand the space for all of the distances below and above the characteristic length determined by σ, respectively. The exponents A, B and a, b control the rate at which the distances are switched from 0 to 1 in the high- and low-dimensional (low-d) space, respectively.[59] To reduce the computational overhead, the sketch-map allows building the low-d embedding using a relevant subset of points from the entire data set (landmark points). Then, the remaining points are projected on the low-d space through a robust out-of-sample procedure.[59,60] In this work, we used a total number of 1000 landmark points, and the sigmoid functions were tuned with the following parameters: σ = 2.5, A = B = 12, a = 1, and b = 2. Trial sketch-map analyses were performed with different combinations of parameters. As expected, the choice of σ was found to be the most critical for achieving an effective dimensionality reduction. The value of 2.5 for σ turned out to be the best compromise to project distinct conformations of the filter in farther regions of the low-d space while preserving a good resolution detail on similar conformations. A density-based cluster analysis (DBSCAN) was then performed on the projections of the trajectories on the low-d space using the dbscan package running under the R 4.0.2 programming environment.[57] For every subset of simulations, the minimum number of points was set to 3 while the epsilon parameter was equal to 0.8 except for the AMBER-KK series where a value of 0.1 was used. The different parameters used in the latter subset of simulations was necessary to cope with the different data structure showing a higher density of points. The cluster analysis was limited to the most informative region of the sketch-map space, which is the one comprised between −60 and 60 units in both dimensions.

Potential sites for the interaction of nonannular lipids and other hydrophobic molecules with the TRAAK channel have been described in structural biology studies.[38,39] The movement of lipids toward the cavity of the channel was also analyzed, and the projection of the distance in the membrane plane of selected atoms of a single lipid molecule on the pore axis was considered. In particular, these distances were computed using the XYDISTANCES multicollective variable function available in PLUMED-2.4.[28] The opening of the fenestration was monitored by evaluating the distance between the TM2 and TM4 helices of the two subunits (A and B) of the channel (TM2A-TM4B and TM2B-TM4A distances). Specifically, the separation between the helices was determined by computing the distance between the Cα-carbon of Pro155 (TM2) and Ile279 (TM4) as a reference.

Results

The TRAAK channel in a fully conductive state (up/up configuration of the TM4 helices) was embedded in a lipid membrane and simulated under conditions of symmetrical ion concentration and a positive electrostatic gradient. Two initial configurations of ions and water in the SF were considered according to the KWK and KK models of ion conduction (Figure b).[17,19,25,26] For each of these initial states, three replicas were considered both with the CHARMM36m[45] and the AMBER ff14SB[46] force fields (see also Table ). Furthermore, for each setup, two additional replicas were simulated in equilibrium conditions (ΔV = 0 mV) as a control. In the following, we will outline the major similarities and differences from the analysis of the conformations of the selectivity filter, ion conduction, and the behavior of the transmembrane helices of TRAAK that are directly involved in channel gating in the trajectories obtained with both force fields considered.

Conformational Preference of the SF

In Table , a list of the MD simulations performed in this study together with the initial simulation conditions and the number of conduction events recorded on a microsecond timescale for each set of simulations are reported. As expected, no conduction events were recorded under equilibrium conditions (ΔV = 0 mV). However, it is remarkable how ion conduction could only be observed in some of the subsets of the simulations performed under electrostatic gradients, namely, AMBER-KK (+100 and +200 mV), AMBER-KWK (+100 mV, but only after the exit of water molecules from sites S2 and S3), and CHARMM-KK (+200 mV). In the case of the CHARMM-KWK subset, ion conduction was never recorded even in the presence of an applied electric field. Moreover, except for AMBER-KK at +200 mV, the number of conduction events was very small. These results not only suggest a dependence of ion conduction on the initial presence of water molecules inside the SF but also a force field-specific behavior that should be carefully analyzed in-depth.

In Figure , the root-mean-square deviation (RMSD) of the backbone atoms of the SF is reported for the four subsets of simulations. The data shows that the initial structure of the SF is better preserved when the simulations are performed with the AMBER force field and are initiated from the KK rather than the KWK configuration. It is also found that the RMSD maximum is more sharply peaked in the simulations with the AMBER force field compared to the simulations with CHARMM. The most striking behavior is observed for CHARMM-KWK where three maxima are found, implying a substantial flexibility of the SF that is not detected in the other subsets of simulations.

Figure 2

Comparison of the root-mean-square deviation (RMSD) of the SF backbone atoms evaluated in the four sets of simulations. Blue, green, orange, and red correspond to the AMBER-KK, AMBER-KWK, CHARMM-KK, and CHARMM-KWK sets, respectively.

To obtain a more informative picture of the conformational space spanned by the SF in the different subsets of simulations, we resorted to dimensionality reduction techniques such as the well-established PCA and the more recent sketch-map approach.[59] In all of the cases, dimensionality reduction was performed on the dihedral angles of the SF backbone plus the χ1 torsions of Thr129/238 for a total of 44 degrees of freedom as a whole (see the Methods section for further details). In this way, the 44-D conformational space of the SF was reduced to a 2-D map that could be easily displayed and further analyzed aiming at identifying major conformational states and recurrent conformations among the distinct subsets. In Figure S1, the low-d space obtained through PCA is shown. Both with the PCA and sketch-map, the conformations explored by the SF in the AMBER-KK set were restricted to a small portion of the map (blue points) reflecting smaller fluctuations compared to those experienced by the other subsets. The opposite behavior was observed in the CHARMM-KWK set in which the projections of the MD trajectories spanned a substantial area of the low-d space. This behavior can be better appreciated in the sketch-map plot shown in Figure . A sketch-map is a nonlinear MDS method that attempts to preserve the middle-ranged distances of the high-dimensional data structure, under the assumption that these encode the most informative features to describe the process under investigation.[59,60] Thus, in the low-d embedding, points closer than a certain threshold distance (σ = 2.5 in this work) are collapsed, while the remaining ones are separated further away providing an intuitive picture of the whole conformational space sampled by the system. As a result, similar conformational states are grouped together leading to the typical island-like appearance of a sketch-map plot like the one reported in Figure . More precisely, regions that are denser in points correspond to distinct metastable states of the SF that were either specific or shared among different subsets of simulations. To better compare the states, a density-based cluster analysis was performed on the projection of the trajectories of each subset of simulations onto the low-d space obtained by the sketch-map. Representative configurations of recurrent conformational states of the SF identified by cluster analysis are shown in the insets of Figure . More details regarding the results of cluster analysis can be found in Figure S2 and Table S1 of the Supporting Information, SI.

Figure 3

Sketch-map representation of the conformational space sampled by the SF during all of the simulations considered in this study. The recurrent conformational states observed are labeled from I (native/conductive state) to VII. Representative configurations of each conformational state are shown and mapped onto the low-d space with solid colored circles. Blue, green, orange, and red correspond to the AMBER-KK, AMBER-KWK, CHARMM-KK, and CHARMM-KWK sets, respectively. The size of the points is proportional to the populations of the main states using a logarithmic scale. Chains A and B are colored in white and gray, respectively. Flipped carbonyls are highlighted in yellow.

The most striking feature of the sketch-map plot reported in Figure is that the native and conductive conformation of the SF (state I) is located at the center of the low-d space, while all of the other conformational states around it correspond to conformations where one or more carbonyl groups are flipped away from the channel axis. As reported in Table , the canonical conductive state of the SF prevails over other conformations in the AMBER-KK set. In AMBER-KK simulations, only 2.2% of the configurations deviated from the canonical conductive state, by means of carbonyl flipping of glycine residues in S0. States II and II′ represent two symmetric conformations of the SF characterized by the flipping of Gly240 in chains A and B of the channel, respectively, affecting the boundary between S4 and S3 sites. These states were observed in the AMBER-KWK, CHARMM-KK, and CHARMM-KWK sets but with different probabilities (see Table ). In particular, the total probability of visiting states II-II′ is as low as 1.0% in the case of the CHARMM-KK set and as high as 47.1% in the CHARMM-KWK simulations, yielding the conformation of the SF with one carbonyl flipping at the S4–S3 boundary as the predominant one in the latter subset of simulations. Other conformational states of the SF involving a single carbonyl flipping that can be observed in the low-d space are states III and IV. These are significantly populated only in the AMBER-KWK and CHARMM-KK sets but with similar probability ratios (Table ) and they involve the carbonyl flipping at the S3–S2 and S1–S0 boundaries, respectively. All of the remaining states identified by the sketch-map represent multiflipped conformations of the SF. For example, in state V, the SF is flipped at both the S4–S3 and S3–S2 boundaries, and while this conformation is observed in all subsets of simulations except the AMBER-KK set, only in the CHARMM-KWK set it significantly contributes to the pool of sampled conformations (4.8%, see Table ). Conversely, state VI involves one carbonyl flipping at S2–S1 and two symmetric flipping at the S1–S0 boundaries. State VI is observed with similar probabilities only in the AMBER-KWK and CHARMM-KK sets of simulations. Finally, state VII is an example of a multiflipped conformation of the SF that involves all of the boundaries between sites that can be subjected to carbonyl flipping, which are S3–S2, S2–S1, and S1–S0, and it is significantly populated only in the CHARMM-KK set. The states described above represent recurrent conformations that account for a substantial portion of the conformational space of the SF in all of the subsets of simulations except for the CHARMM-KWK subset, in which almost half of the conformations are specific to this simulation condition (44.7%, Table ). This situation is manifested in the sketch-map plot of Figure in which several spots are only populated by conformations that come from CHARMM-KWK trajectories (see the red dots in Figure ).

Table 2

Recurrent Conformational States of the SF Mapped onto the Sketch-Map Spacea

		relative population (%)
		AMBER		CHARMM
conformational state	carbonyl flipping boundaries	KK	KWK	KK	KWK
I	NA	97.8	43.2	52.3	2.7
II	S4–S3		11.0	0.7	35.5
II′	S4–S3		4.6	0.3	11.6
III	S3–S2		16.6	7.1
IV	S2–S1		8.5	4.0
V	S4–S3 and S3–S2		0.4	0.2	4.8
VI	S3–S2 and S1–S0		13.3	10.3
VII	S0–S1, S1–S2, S2–S3			7.6	0.7
others	NA	2.2	2.4	17.5	44.7

Full results of cluster analysis are reported in Table S1.

Stability of the TM Region of the Channel and Interactions with Lipid Molecules

Like the other members of the TREK subfamily of K2P channels, TRAAK is gated by membrane stretch and the positive curvature of the membrane. Even though some contradictory findings have been reported,[38,63] it is widely accepted that the active state of the channel is represented by the up/up conformation, with the TM4 helices of both subunits kinked upward (Figure a), while the down state is regarded as a low conductance one. According to a proposed gating hypothesis, the low conductance associated with the down state can be explained by impairment of ion conduction due to the presence of lipid tails reaching the cavity through the lateral fenestrations.[37] Upon membrane stretch, the hindrance imposed by the lipids would be released, and the up/up state with sealed fenestrations would be adopted. More recently, an allosteric effect directly linking the “down to up” transition of TM4 with the activation of the SF has been proposed[64−66] and investigated by several groups using different computational approaches.[67−69] In this model, the SF would enter into a nonconductive state when the down conformation of TM4 is adopted, resembling the C-type inactivation process observed in Kv channels.[66]

Figure 4

(a) Representative configuration of the up and down states of the protein at one of the interfaces between the two subunits of the channel. The fenestration that appears in the down state at the TM2-TM4 interface. (b) Scatter plots of the distance between the TM2 and TM4 helices of different subunits obtained from the AMBER and CHARMM sets of simulations are highlighted. Blue and orange circles correspond to the maxima of the AMBER and CHARMM distributions, respectively. A and B refer to the two subunits of the channel.

In Figure b, a scatter plot is reported showing the distance between the TM4 and TM2 helices calculated for both subunits (i.e., TM4A-TM2B and TM4B-TM2A) in the sets of simulations performed with the AMBER and CHARMM force fields (see also the individual plots for each simulation run in Figure S3). When the two helices interact, the distance between the two reference carbon atoms (see Methods for further details) approaches a value close to 4 Å (up state), while distances above 8 Å are an indication of a substantial separation between the helices with the concurrent opening of the fenestration (down state). All of the simulations were initiated from a fully conductive up/up state, and in both sets of simulations, “up to down” transitions were recorded (one and five transitions in the AMBER and CHARMM simulations, respectively, see Table ). Notably no backward transitions to the up/up state were ever sampled in this study. As the plot in Figure b shows, in the AMBER set of simulations, there is a more pronounced tendency to preserve the up/up state compared to the CHARMM set. Regardless of the number of complete up to down transition events, which are more frequently sampled with the CHARMM force field, the maximum in the distribution of the distances is found at about 4 Å in the AMBER simulations. This value increases by more than 1 Å in the CHARMM simulations which illustrated that the initial conditions are better preserved with the AMBER force field.

The conformational state adopted by TRAAK in MD simulation is influenced by the behavior of lipid molecules at the interfaces between the two subunits of the channel. In Figure a, the probability of observing a series of reference atoms of a lipid molecule at a given distance from the channel axis in the two sets of simulations is shown. A representative lipid molecule with the atoms used for the analysis is shown in Figure b. The degree of penetration of lipid molecules inside the pore of the channel was measured by considering the boundary between the cavity of the channel and the membrane environment at about 10 Å. We observed a substantial penetration of lipids in the simulations performed with the CHARMM force field (see the lower tail of the distribution for atoms P, C21, C36, C218, and C316), while this behavior is much less pronounced in the simulations with the AMBER force field. Values of the distance between the P atoms of the lipids and the pore axis lower than 10 Å observed in the AMBER set of simulations are associated with configurations where the POPC headgroup only slightly faces the cavity from the intracellular side of the membrane (see Figure a). Conversely, in the CHARMM simulations, lipids enter the cavity in an almost complete way, remaining anchored to the membrane through the Sn-1 chain (see the distribution of atom C316 in Figures a and 6b) and occasionally, pointing the headgroup toward the SF. In spite of this, lipid tails are frequently observed in the cavity even without invoking full lipid excursions. However, as previously reported, this occurrence is unlikely to have any functional implication in ion conduction.[36] The aggregated data shown in Figure a is also reported as individual plots for each simulation run in Figures S4 and S5. From these plots, it can be observed that most of the lipid occlusion events sampled in the CHARMM simulations showed the tendency of a single POPC molecule to enter the cavity entirely (see Figures S5 and S6).

Figure 5

(a) Presence of lipid atoms in the cavity of the TRAAK channel in the AMBER and CHARMM sets of simulations. (b) POPC molecule in stick representation highlighting in yellow the atoms used in the analysis presented in (a) with their corresponding labels.

Figure 6

Comparison between typical configurations of lipids facing the internal cavity of the channel obtained with the (a) AMBER and (b) CHARMM force fields. The orientation of the protein is identical in both snapshots. Lipids are depicted in van der Waals representation in blue and orange, respectively. The protein is depicted in gray and the up state is shown in (a), while the down state is displayed in (b).

Discussion and Conclusions

The microsecond-long MD simulations performed in this study under different initial conditions facilitate exploring the effects of either vacancies or water molecules inside the SF and the choice of the force field in the simulation of ion conduction in the TRAAK channel. The most important observation emerging from our data is that the combination of the AMBER force field with the absence of water molecules inside the SF as the initial configuration is a requirement to obtain a significant number of conduction events at both +100 and +200 mV. This occurrence is strongly related to the fact that only in this subset of simulations, the native conductive state of the SF was preserved in the vast majority of the MD runs (more than 95% of the total simulation time, see Table ). Conversely, when simulations performed with the same force field are initiated from the KWK configuration, the probability of finding the filter in the native conductive configuration drops at about 43% of the total simulation time and so does the total number of conduction events (from 8 to 3 and from 22 to 0 at +100 and +200 mV, respectively). Interestingly, the other significantly populated states observed in this set of simulations mostly involve carbonyl flipping at the S4–S3 boundary (states II and II′, 15.6% in total) or the S3–S2 boundary (state III, 16.5%), even though other multiflipped states involving these coordination states can be found (see Table ). The same flipping is observed in CHARMM-KWK where states II and II′ account for about 47% of all of the conformations visited by the filter, and in CHARMM-KK but with much lower probabilities. This is an indication that the presence of water molecules triggers the conformational change that brings the SF away from the native conductive state. Once flipping at the S4–S3 boundary occurs with a water molecule in S3, it is difficult to recover the conductive state in a microsecond regime. Moreover, we highlight that all of the conduction events observed in this study took place according to the KK mechanism of ion conduction. Notably, the conduction events recorded for AMBER-KWK at +100 mV were observed after the expulsion of the intercalating water molecule from the S3-binding site. The same behavior was not observed in the CHARMM-KWK subset of simulations, where either the corresponding water molecule remained in the SF or the SF became too unstable to allow ion conduction. The observations described above are consistent with the results of a recent study performed on the KcsA channel.[34]

Concerning the comparison between simulations performed under the same initial conditions but employing different force fields, the most striking result is that the probability of finding the SF in the native conformation is reduced by half going from the AMBER to the CHARMM force field, when the KK configuration is chosen to start with, and conduction is virtually suppressed if the starting configuration is the KWK one. This is a clear indication that AMBER and CHARMM force fields behave differently on a microsecond timescale, highlighting an incipient instability of the SF when the CHARMM force field is considered that can be critical for describing the dynamics of the system under investigation. These findings are in line with the results obtained from the comparison of the two force fields in preserving the native state of the SF for the KcsA channel.[34] However, while in the case of KcsA the CHARMM force field leads to a full collapse of the SF on the microsecond timescales regardless of the choice of Lennard-Jones parameters for the ions, in the present work the collapsed state of the SF was actually never sampled. This observation suggests that the choice of the force field might influence the structural properties of the simulated channels in a subtler way than it was observed in the other study, also showing a case-dependent behavior. Even though the main goal of this work was to highlight the different states of the SF sampled under different initial conditions, we also attempted to estimate the conductance of the channel using the number of ion translocation events recorded in the independent MD runs (Table ). The estimated conductance in the AMBER-KK subset of simulations was found to be 4.3 ± 1.9 pS at +100 mV. Thus, the experimental conductance of 65.4 ± 3.2 pS at +100 mV under symmetrical 150 mM KCl is underestimated by more than 1 order of magnitude.[70] Again, the discrepancies between the experimental and estimated values of the conductance are in line with those obtained for the KcsA channel from simulations performed with the AMBER force field and under the same transmembrane potentials.[34]

Since TRAAK is activated by polymodular stimuli including membrane stretch, in this work, we have also investigated the behavior of the TM helices directly involved in the gating functionality of the channel. Our data clearly suggest that the initial conductive state of the channel (up/up state) is better preserved in microsecond-long simulations employing the AMBER force field. Indeed, while up to down transitions are found in both sets of simulations, these are more frequently observed in the CHARMM simulations. In addition, even when configurations consistent with the up/up state are considered, these appear to be intrinsically less structured in the CHARMM trajectories compared to AMBER ones. This is evident in the broader distributions of distances between the transmembrane helices TM2 and TM4, shown in Figure b. It is noticeable that without imposing any membrane stretch to keep the channel in an active conformation, the up to down conformational transition is expected to occur. However, the rate at which this event is observed and its manifestation in several of the CHARMM trajectories highlight a clear difference of this force field in describing the conformational dynamics of the TM4 helix in comparison to the AMBER force field. Results obtained with the AMBER force field are in line with those previously reported by Harrigan et al. in simulations of the closely related TREK-2 channel using the same AMBER force field.[67] These authors described recurrent and reversible activation events of the channel taking place on a multi-microsecond timescale, the down conformation being the preferred state under equilibrium pressure conditions.[67] The channel activates upon membrane stretch, so in equilibrium conditions (without stretch), the closed state, corresponding to the down conformation, should be the most likely conformational state.

In the present study, the lower structural stability of the TM region of the channel experienced in the simulations using the CHARMM force field seems to be correlated with an increased probability of lipids reaching the cavity from the open fenestrations present between TM2 and TM4 in the down state. Notably, in the case of CHARMM, these events involve the insertion of almost an entire lipid molecule into the cavity, with the headgroup directly facing the bottom of the S4 site of the SF. Since the presence of lipid tails in the cavity of this subfamily of K2P channels has been suggested to be implicated in the mechanosensing gating mechanism, our results suggest that care must be taken when interpreting data sets from MD simulations. On a final note, we would like to emphasize that the above observations are based on the assumption of a causative relationship between the stability of the protein and the presence of lipids inside the cavity of the channel. However, the two sets of simulations (AMBER and CHARMM) also differ in the force field parameters describing the lipid molecules. Therefore, while it is unlikely that lipid penetration and conformational changes of the TM4 helix are not correlated, this hypothesis cannot be ruled out, nor that the diverging behavior of the TM helices of the protein can be an effect of different lipid parametrization.

Currently, MD simulations are widely employed to study biomolecular systems at a fully atomistic level. The ever-increasing availability of high-performance computational facilities is narrowing the gap between the timescales able to be sampled by simulations and experiments. This, in turn, expands the scope of the computational investigations considered and allows us to benchmark the force field developed in the 1980s using picosecond timescales for the study of processes occurring on much longer timescales. In this work, we compared results from the AMBER and CHARMM force fields with the aim to test their ability to describe several functional properties of the TRAAK channel, including the stability of the SF, the occurrence of ion permeation, and the preferential states adopted by the transmembrane regions of the protein directly involved in gating. While these two force fields behave similarly on the nanosecond regime, results presented in this study demonstrate that they render a strikingly dissimilar picture of the channel dynamics on longer timescales. Even though simulations performed with the AMBER force field under applied membrane potentials of magnitudes comparable to experimental ones provide a more realistic description of ion conduction for the TRAAK channel over the CHARMM force field, the computed values of conductance are still underestimated by an order of magnitude, leaving room for further development in the field.