Mechanosensitive channels of small conductance (MscS) are ubiquitous turgor pressure regulators found in many walled cells and some intracellular organelles. Escherichia coli MscS acting as a tension-activated osmolyte release valve shows a nonsaturable conductance (1.2 nS in a 39 mS/cm electrolyte) and weak preference for anions. Pursuing the transition pathways in this channel, we applied the extrapolated motion protocol (cycles of displacements, minimizations, and short simulations) to the previously generated compact resting conformation of MscS. We observed tilting and straightening of the kinked pore-forming TM3 helices during the barrel expansion. Extended all-atom simulations confirmed the stability of the open conformation in the bilayer. A 53° spontaneous axial rotation of TM3s observed after equilibration increased the width and polarity of the pore allowing for stable voltage-independent hydration and presence of both cations and anions throughout the pore. The resultant open state, characterized by a pore 1.6 nm wide, satisfied the experimental conductance and in-plane expansion. Applied transmembrane electric field (±100 to ±200 mV) in simulations produced a flow of both K+ and Cl−, with Cl− current dominating at higher voltages. Electroosmotic water flux strongly correlated with the chloride current (∼8 waters per Cl−). The selectivity and rectification were in agreement with the experimental measurements performed in the same range of voltages. Among the charged residues surrounding the pore, only K169 was found to contribute noticeably in the rectification. We conclude that (a) the barrel expansion involving tilting, straightening, and rotation of TM3s provides the geometry and electrostatics that accounts for the conductive properties of the open pore; (b) the observed regimen of ion passage through the pore is similar to electrodiffusion, thus macroscopic estimations closely approximate the experimental and molecular dynamics-simulated conductances; (c) increased interaction of the opposing ionic fluxes at higher voltages may result in selectivities stronger than measured near the reversal potential.
Introduction
Mechanosensitive channels MscL and MscS act as turgor pressure release valves in bacteria (Levina et al., 1999) and both undergo large-scale conformational transitions driven directly by tension in the cytoplasmic membrane (Sukharev et al., 2005, 2007; Perozo, 2006; Blount et al., 2007). These channels are convenient systems for basic biophysical studies as they generate robust current responses in patch-clamp experiments in native or reconstituted membranes and their crystal structures are available (Chang et al., 1998; Bass et al., 2002; Steinbacher et al., 2007). A similar mechanism of direct activation by stress in the lipid bilayer has been suggested for several types of eukaryotic mechanosensitive channels (Patel and Honore, 2001; Zhou et al., 2005).
It is commonly accepted that in order to detect mechanical stimuli, the receptor must comply with external forces driving the molecule's transition into an active state (Corey and Hudspeth, 1983; Sachs and Morris, 1998; Markin and Sachs, 2004; Wiggins and Phillips, 2005). The work produced by the force biases the relative energies of distinct conformational states in the receptor, which is manifested as changes in transition rates and state occupancies with tension. According to generalized models of stretch-activated channel gating, the spatial scale and free energy of the transition can be estimated from open probabilities (Po) measured as a function of tension (Sukharev et al., 1999; Chiang et al., 2004). The crystal structure of MscL from Mycobacterium tuberculosis representing the closed state of the channel (Chang et al., 1998) and the derived homology model of Escherichia coli MscL (Sukharev et al., 2001b) along with experimentally estimated spatial parameters gave first structural grounds that allowed us to predict the iris-like character of helical movements during gating transitions in these proteins (Sukharev et al., 2001a), which was confirmed by disulfide cross-linking (Betanzos et al., 2002) and EPR (electron paramagnetic resonance) studies (Perozo et al., 2002). The studies have also demonstrated that the thermodynamic spatial parameter of the transition (in-plane area change) can be reliably used as a constraint in the modeling of the gating process.
In contrast to MscL, which is largely prokaryotic, mscS-like genes have been identified in essentially every group of walled organisms, both prokaryotes and eukaryotes (McLaggan et al., 2002), making MscS a more general model. MscS orthologues have been shown to play critical roles in the maintenance of chloroplasts in unicellular alga (Nakayama et al., 2007) and higher plants (Haswell and Meyerowitz, 2006). E. coli MscS, the best studied representative of this family, acts as an osmolyte release valve regulating turgor pressure. It exhibits a complex functional cycle involving closed, open, desensitized (mode-shifted) and completely inactivated states (Akitake et al., 2005, 2007). Its crystal structure, initially deemed to represent the open state (Bass et al., 2002), was found to be nonconductive at physiological voltages in simulations (Anishkin and Sukharev, 2004; Sotomayor and Schulten, 2004; Spronk et al., 2006). When embedded into a bilayer without restraints in all-atom molecular dynamics (MD) simulations, this structure was unstable, displaying a quick asymmetrical collapse of the pore (Sotomayor and Schulten, 2004), which posed a number of questions as to what functional state it might represent or resemble. It has been suggested that at least the resting conformation of the channel must be more compact than the splayed crystal structure (Miller et al., 2003).
There were several attempts to envision MscS opening in MD simulations by applying membrane tension, steering forces or high electric fields on the crystal structure (Sotomayor and Schulten, 2004; Sotomayor et al., 2006; Spronk et al., 2006). The expanded crystal structure often featured a fully hydrated pore; however, it always displayed a strong anionic selectivity in both explicit MD (Spronk et al., 2006; Sotomayor et al., 2007) and Brownian dynamics MC (Sotomayor et al., 2006; Vora et al., 2006) simulations. In many instances these preexpanded models appeared to be stable and conductive only under conditions of high transmembrane voltage. Because in reality MscS has a stable long-lived open state, weak selectivity, and nearly linear conductance around zero voltage, these first models unlikely represented a fully open state. They rather pointed out that more profound rearrangements in the barrel should occur in order to stabilize the conductive pathway and change the electrostatics in the lumen.
In the preceding paper (Anishkin et al., 2008) we presumed that the sources of instability of the crystal structure in the previous all-atom MD simulations (Sotomayor and Schulten, 2004; Spronk et al., 2006) were the splayed conformations of the lipid-facing helices (TM1 and TM2), the absence of unresolved N-terminal segments, and improper positioning of the protein complex relative to the bilayer midplane resulting in a dislocated lateral pressure profile acting on the protein wall. To alleviate the stress in the protein caused by interactions with lipids, we transformed the crystal structure by aligning the TM1–TM2 helical pairs in a more parallel fashion along the central (TM3) shaft and by adding Rosetta-predicted N-terminal segments. The transformations were done with the fast “extrapolated motion” protocol in vacuum, and only after substantial rearrangement the energy-minimized model was equilibrated in all-atom MD. The models revealed that the most compact packing of transmembrane helices can be achieved when the characteristic TM3 kink is moved two helical turns down, from its crystal position near G113 to more conserved G121. As a result, the new structure well fitted the hydrophobic profile of the lipid bilayer, which was positioned more periplasmic when compared with the previous assignments (Sotomayor and Schulten, 2004). This compact resting conformation was stable in MD simulations in the explicit fully hydrated lipid bilayer, remaining nonconductive. The mutual positions of the peripheral (TM2) and core (TM3) helices in the resting state were then experimentally tested using disulfide cross-links (Anishkin et al., 2008). The recent assessment of the resting MscS using EPR combined with computational modeling suggested a comparable conformation with the helical packing more tight than in the crystal structure and N termini extensively interacting with the headgroups in the periplasmic leaflet of the membrane (Vasquez et al., 2008).
In the present paper, we attempt to predict the open conformation of MscS by first modifying the closed structure in vacuum using the computationally inexpensive extrapolated motion technique (Anishkin et al., 2008). From a set of extrapolated trajectories we choose candidate conformations for the open state satisfying the experimental parameters of conductance and lateral expansion, and then refine these models in extended all-atom simulations. The extrapolated trajectories reveal a consistent picture of the MscS barrel expansion, whose limits are structurally defined by straightening and tilting of TM3s. The refinement resulted in a pore-stabilizing axial rotation of TM3 helices that follows expansion. We present analysis of the pore electrostatics in this new open conformation, macroscopic estimations, and explicit MD simulations of its conductance and compare them with experimental conductances measured under similar conditions. The model with straightened TM3s well satisfies the conductive properties of the pore, where the electrolyte behaves essentially in a bulk-like fashion. The conductance simulations also indicate substantial coupling of ion and water fluxes predicting for the first time a dependence of channel selectivity on voltage.
Materials And Methods
Molecular Simulations
All simulations were performed using NAMD2 (Phillips et al., 2005) with the CHARMM27 forcefield (MacKerell et al., 1998), Particle-mesh Ewald (Darden et al., 1993) method for long-range electrostatics estimation, 10 Å cutoff for short-range electrostatic and van der Waals forces, and Langevin thermostat set at 310 K. VMD (Humphrey et al., 1996) was used for visualization, molecular modifications, and analysis using embedded Tcl language, providing flexible and convenient environment for analysis of trajectories. Scripts were written for calculations of protein in-plane areas, tilts of helices, pore diameters and asymmetry, macroscopic estimations of conductance, computation of ion concentrations in predefined segments of the pore, detection of permeation events, as well as for computation of ionic and coupled water fluxes. The extrapolated motion protocol (Anishkin et al., 2008) was also implemented as a Tcl script run with VMD and NAMD. The previously described compact resting state model, generated from the original crystal structure (PDB ID 1MXM), was used as a starting conformation. This original crystal structure and the refined version (PDB ID 2OAU) are equally usable for explorations of the transmembrane part of the protein since substantial refinements involving registry shifts took place only in the cytoplasmic “cage” domain. The compact model was expanded in repeated extrapolation trials in vacuum, and multiple conformations were assessed for compatibility with the experimental parameters of the open state (estimated conductivity and in-plane expansion) and conflict-free transitions to and from the closed state. The candidate open-state models were then equilibrated and simulated an all-atom setting. Upon equilibration, the electrostatics of the model and its conductance under physiological voltages were computed.
Predictions of the Transition Pathway Using the Extrapolated Motion Protocol
Gating transitions in MscS occur in microsecond time scale (Shapovalov and Lester, 2004). To overcome limitations of MD simulations in an all-atom system where large-scale and slow macromolecular motions are prohibitively expensive to compute, we explored such motions of the MscS barrel in vacuum using sequences of fast extrapolation-minimization cycles. The method has been described previously (Akitake et al., 2007; Anishkin et al., 2008), here we only outline its principle and provide specific parameters used in MscS expansion simulations in brief. An initial small displacement of a specific domain or all protein atoms in arbitrary directions is the first step that is followed by energy minimization and a short relaxing MD simulation. The new positions of atoms are now “corrected” by these relaxing steps in both absolute value and direction, indicating the local pathway preferred by the system itself, which can be different from the initial arbitrary displacement. In the next cycle, displacements for each atom are calculated through the linear extrapolation of coordinates based on the previous and current (minimized) positions. After the second displacement, the energy is minimized and a relaxing MD simulation is repeated again. A sequence of such cycles produces a smooth curvilinear trajectory of self-permitted protein conformations, which can be repeated many times. Finding an energetic obstacle, the protein may stop its motion and reverse the direction. The range of protein motion covered by a trajectory can be changed by introducing “amplification coefficient” (usually chosen between 1.0 and 1.2), which scales the absolute value of the coordinate extrapolation vector and thus helps to overcome local barriers. Importantly, the initial assignment of thermal velocities is random, and Langevin temperature control at each relaxing MD step imposes random forces on the molecule. For these reasons, trajectories initiated from the same configuration are not identical. While requiring sufficient statistics, this randomness permits more extensive exploration of the conformational space.
Previously, exploring compact resting conformations of the MscS transmembrane barrel, we used 0.1–1 Å displacements in the radial direction to initiate expansion or compression. This time, we simulated the initial conformation for 1 ps, and the accumulated drift of the system, averaged among the seven identical subunits, was used as the initial displacement. Each extrapolation cycle involved 100 steps of energy minimization, a 1-ps relaxing MD simulation at 310 K, and another energy minimization (100 steps) with sevenfold symmetry restraints. The amplification coefficient was set to either 1.00, 1.05, or 1.10, the later allowed to explore a broader range of expansions.
Estimations of the In-Plane Areas, Pore Cross Sections, and Resistance
The barrel was divided into thin (0.5 Å) sections normal to the z axis (pore symmetry axis) and the area of each section was computed within its solvent-accessible boundary. The solvent-accessible surface was created with a 1.4-Å probe radius using the CHARMM27 VDW radii for the protein. The radius of a circle equal in its area to the area of the slice was taken as the effective protein radius at that particular z level. Smooth radii profiles are presented as the twofold nearest-neighbor average of r over z slices. In the graphs presented below the z coordinate corresponding to the midplane of the membrane was set as zero. Correspondingly, effective in-plane radii for the periplasmic and cytoplasmic rims of the protein were calculated for the ranges 12.5–17.5 Å and −17.5 to −12.5 Å, which approximately corresponds to the levels of the carbonyl oxygens in an unperturbed bilayer. The same method was used for calculations of effective pore radii of the outer hydrophobic chamber (−13 < z < −7 Å) and pore constriction (−23 < z < −17 Å). The computations of cross-sectional areas for the pore allowed us to estimate total pore resistance in a continuum conductor approximation (Hall, 1975; Hille, 1992). The ion-accessible area in each slice was determined with a 2.83-Å probe radius, which corresponds to the average ion radius (between K+ and Cl−) with half of its first hydration shell. The access resistances were determined from the effective radii of the cytoplasmic and periplasmic pore entrances (Hall, 1975). Using the same approach, resistance of the side windows of the cytoplasmic cage was estimated based on their cross-sectional area profiles. The total resistance of the cytoplasmic cage was represented by seven parallel resistances of individual windows connected in series with the transmembrane pore. Addition of the cytoplasmic cage was estimated to diminish the conductance of the open barrel by ∼38%. In this approximation, the resistance of the cage is independent of voltage. Experimental studies (Koprowski and Kubalski, 2003) indicate that in the course of the gating cycle, the cage might experience significant conformational changes, but the exact dynamics of this domain is not well understood for any of these states. Thus, at the present stage any calculations of the cage resistance, implicit or explicit, should be considered as estimations only.
Analysis of Expansion Trajectories
To visualize the trajectories of the gate-bearing TM3 barrel and characterize the distribution of outcomes at every iteration step, helical positions were scored in terms of their distance from the pore axis and tilt angle relative to the axis. The backbone of TM3a helix (residues V96 to N112) was fitted with a line, and the helical tilt was defined as the angle between the line and the pore axis. The average surface-to-surface distance at the level of the gate (residues L105 and L109) calculated with CHARMM27 atomic radii was taken as the constriction diameter.
Molecular Dynamics Simulations
The open-state model satisfying the experimental parameters was truncated to reduce the size of the system. Only the transmembrane domain with the adjacent portion of the cytoplasmic cage (residues 1–130, 137–140, 147–154, and 162–175 of each subunit) was included. During the simulations, the α carbons of the terminal residues (points where the cytoplasmic domain was truncated) were softly restrained. The barrel was equilibrated in a fully hydrated POPC bilayer (220 lipids) with TIP3P water (Jorgensen et al., 1983) in a flexible hexagonal periodic box at 1 atm pressure (Langevin piston method). It was then simulated in an unrestrained mode (for 35 ns in total) with periodic (every 4 ns) symmetry-driven simulated annealings (1 ns each). In most trials, the membrane tension was set at 10 dyne/cm (Gullingsrud et al., 2001; Gullingsrud and Schulten, 2004), and additional simulations were performed at 50 dyn/cm. During the symmetry annealing stages all protein atoms were driven toward their continuously updated sevenfold symmetric average positions using gradually stiffening harmonic restraints. The rest of the system (lipids, water, and ions) was unrestrained. In all of the unrestrained and annealing simulations the number of potassium and chloride ions corresponded to a 200 mM salt. The total system size was ∼100,000 atoms. Additional technical details of the simulations can be found in the previous paper describing the generation of the resting state model (Anishkin et al., 2008).
Explicit Computation of Conductance and Electrostatic Maps
Estimations of conductance in explicit all-atom MD simulations were done essentially as described in Aksimentiev and Schulten (2005). All conductance calculations used a model of the open state obtained in extrapolated-motion simulations, which was truncated and subjected to unrestrained relaxation in the explicit lipid bilayer (4 ns at 10 dyne/cm) followed by 1-ns symmetry-driven simulated annealing. In this structure, TM3 helices have already experienced axial rotation resulting in both partial retraction of the L105 and L109 side chains from the pore lumen and exposure of glycines that increase the diameter and polarity of the pore lining. In our conductance simulations, the protein backbone was restrained softly (0.01 kcal/mol/Å2) to allow for conformational freedom yet prevent large-scale drift, whereas the salt concentration in the aqueous phase surrounding the channel was increased to 1500 mM. All other parameters were as described above for simulations in the explicit POPC bilayer under 10 dyne/cm. Electrostatic maps were calculated using the “PME electrostatics” plug-in over a 4-ns trajectory obtained without an electric bias across the membrane as in Sotomayor et al. (2007). Ion movement was followed at ±50, ±100, ±150, and ±200 mV across the entire simulation cell of 106 Å, with the potential on the cytoplasmic side assigned as zero, for the duration of 8–16 ns at each voltage. Ion crossing events in the pore were scored using custom-written Tcl scripts in VMD.
To provide a “benchmark” for electrolyte conductivity of the bulk solution, we simulated a cube (80 Å side) filled with 1500 mM KCl solution (15,500 waters, 443 Cl−, and 443 K+ ions, ∼47,000 atoms in total) in the same conditions as the channel above, except that 1 atm pressure was applied isotropically. The box was simulated for 8 ns twice, with the same external electric field as in simulations of MscS to provide +100 and −100 mV potential difference across the simulation cell. Conductive characteristics were estimated with the same algorithms as for the channel. The simulated average conductivity of 1500 mM KCl solution was 109.7 mS/cm, 1.42 times lower than the experimentally measured 156.3 mS/cm for the same solution. However, the simulated conductance was compared with the experimental data recorded not in pure KCl, but rather in KCl containing buffer (below). The measured conductivity of the recording buffer with 1.5 M KCl was 110.4 mS/m, close to the simulated bulk conductivity. Thus, to compare with experimental values, a scaling factor of 1.006 was applied to the simulated conductance.
Patch-Clamp Experiments
The D3N, D8N, R88Q, and K169Q mutations in MscS were introduced using a QuickChange kit (Invitrogen). Wild-type (WT) MscS and the mutants were expressed in MJF465 cells (Levina et al., 1999) from the pB10b vector (Okada et al., 2002). Preparation of giant E. coli spheroplasts and patch-clamping procedures were conducted as described previously (Akitake et al., 2005; Anishkin et al., 2008). Channel recording was conducted on excised inside-out patches in symmetrical buffers containing 200 mM KCl, 10 mM CaCl2, 40 mM MgCl2, and 10 mM HEPES (pH 7.2). In some experiments, KCl concentration was increased to 0.5 or 1.5 M. Mechanical stimuli were delivered using a high-speed pressure clamp apparatus (HSPC-1, ALA Scientific).
To measure the I-V curves in a broad range of voltages two protocols were used. In the first “pressure ramp” protocol, 1-s ramps to subsaturating pressure were applied to facilitate the observation of the openings of a few “early” channels at a given voltage. Repeated several times, these ramps allowed us to capture full opening events that are rare at voltages beyond 100 mV in either direction. I-V curves were then created from single-channel currents for WT and mutants in the range from −150 to +150 mV.
The second “voltage pulse” protocol was used to determine the I-V relationships from the integral current of a small channel population. For that purpose, we used noninduced spheroplasts expressing ∼10 channel channels per patch through the natural “leakage” of the promotor. In that protocol, a saturating pressure was attained along a ramp in 1.0 s and held for additional 2.25 s. At the start of the pressure ramp, the voltage across the patch is stepped from 0 to +5 mV. On attaining the saturation pressure, the patch is stimulated with voltage pulses of 50 ms duration separated by 500 ms. Each subsequent pulse was higher in magnitude by a predetermined step to collect the data in the range from −200 to +200 mV. The potential of +5 mV maintained between the pulses allowed us to monitor the stability of the baseline. Recordings from patches destabilized by pressure and voltage were discarded. The voltage pulses of low to moderate amplitude produced current spikes that had two components, the fast capacitive decay and the constant current. At higher amplitudes, voltage pulses produced slowly decaying inactivation currents in addition to capacitive spikes. By fitting the slowly decaying component of the current with exponents, we determined the initial integral current before inactivation or slippage into subconductive states took place. Both protocols are illustrated in Fig. 7.
Online Supplemental Material
The PDB model of the open state presented in Fig. 3 B and a movie illustrating the simulation of ion permeation through the open pore are available as online supplemental material ().
Results
Fig. 1 shows the helical arrangement of the transmembrane barrel of MscS in its compact resting conformation (Fig. 1 B) and expanded open-like state (Fig. 1 C), both derived from the original crystal structure (pdb 1MXM, Fig. 1 A). To create the resting state model, in the preceding work (Anishkin et al., 2008) we performed initial compaction of the peripheral helices followed by a partial expansion of the entire barrel, attachment of the de novo predicted N-terminal segments (unresolved in the crystal structure), and cycles of expansion–compaction. The compact model was extensively tested in all-atom MD equilibration with lipids and water and found to be stable and nonconductive in the physiological range of voltages. The model was then refined in a 1-ns symmetry-driven simulated annealing (Anishkin et al., 2008). The key features of the resting state differentiating it from the crystal structure include an almost parallel packing of the peripheral helices and the presence of polar N termini, all critically improving the protein match to the surrounding bilayer. When compared with previous assignments (Sotomayor and Schulten, 2004; Spronk et al., 2006), the position of the bilayer relative to the protein is shifted ∼8 Å toward the periplasmic end of the barrel, which puts positive charges on TM1 and TM2 (R46, R54, R74) not in the middle of the hydrocarbon as was suggested previously (Bass et al., 2002), but next to the polar headgroups of lipids. The attachment of TM1–TM2 pairs to the central TM3 shaft restored buried contacts between the lipid-facing domain and the channel gate. In contrast to the crystal conformation bearing a characteristic kink near G113, the TM3 helices in the new resting state bend near G121 (Fig. 1 B), having the kink shifted two helical turns down and thus permitting for a conflict-free packing of TM1–TM2 loops above the splayed TM3b domains.
Extrapolated Motion Simulations
Using this compact resting conformation as the initial state, we performed repeated extrapolated expansions of the transmembrane barrel. The motion was initiated by spontaneous thermal fluctuations accumulated over 1 ps of unrestrained simulation and averaged among seven subunits. The following energy minimization, short (1 ps) relaxing MD simulation and symmetry-restrained minimization produced new coordinates of a slightly expanded structure (average radial expansion of ∼0.1 Å), and this sequence comprised the first iteration cycle. By extrapolating the coordinate displacement vector in the same direction either exactly or with a small amplification coefficient (g) and then repeating the minimization, relaxation, and symmetrization steps, we complete the second iteration. With an amplification coefficient chosen between 1.0 and 1.1, the average increase of the effective radius of the outer boundary of the protein with each cycle varied between 0.01 and 0.17 Å, while the pore radius increased by ∼0.0008–0.07 Å . It took 60–80 iterations to produce a quasi-continuous trajectory toward an expanded state of the transmembrane barrel with a gap-free wall and a diameter of the pore lumen approaching 16–18 Å. We observed that extrapolated transitions initiated by thermal fluctuations produce trajectories similar to the extrapolations initiated by small radial displacement “kicks” (Anishkin et al., 2008). The choice of amplification coefficient was found to be more critical for the character of trajectories than the way the initial displacement was achieved. Fig. 1 (C and F) shows a typical expanded conformation representing the consensus of several states satisfying the experimental conductance and in-plane area change estimated from the slope of measured dose–response curves on tension. This expanded model obtained in extrapolations in vacuum was chosen as a candidate open state for further refinement in unrestrained all-atom MD simulations (below).
Before explaining the refinement, we should clarify how a consensus of 50 repeated opening extrapolations was found. We observed that separate trials started from the same conformation rarely produced identical trajectories. The protocol includes a stochastic component since at every cycle the structure is subjected to a short relaxing MD simulation at 310 K initiated with random initial velocity assignment and maintained through the Langevin dynamics temperature control. Because the extrapolation of coordinates utilizes the previous and the current position for every atom, each new position depends on the variations of atomic positions at every step. Despite keeping the extrapolated conformations symmetrical through the imposed sevenfold symmetry restraints at the end of each cycle, the expansion trajectories varied substantially. We also found that dissimilar trajectories starting from different initial conditions often converge in the same region of the expansion (helical tilt versus pore diameter) diagram and produce very similar expanded states, which suggest clusters of favorable conformations that can be searched for and achieved in several different ways.
Fig. 2 illustrates two series of extrapolated expansions performed with amplification coefficients of 1.00 and 1.05, presented in the coordinates of helical tilt versus constriction diameter (see legend for definitions). The initial conformation of the barrel (Fig. 2 A) was the same, and each trial started with random thermal fluctuations which, when amplified, eventually lead to the opening of the channel. The position of this initial state characterized with ∼6.5 Å pore and 24° tilts of TM3s is denoted with square in Fig. 2 C. When “kicked” with a thermal fluctuation and extrapolated with g = 1.00, the barrel typically experienced limited expansion (blue trajectories), roaming around the initial state in the range of diameters between 5 and 10 Å and tilts between 20° and 30°. A smaller group of conformations assumed higher helical tilts (∼35°), however, without being able to expand beyond 10 Å. The states characterized with smaller, slower growing tilts (20°–30°) made more frequent excursions to expanded conformations with the pores up to 15–17 Å in diameter. This range agrees well with the earlier electrophysiological measurements (14–18 Å, Sukharev, 2002) and our previous estimations of required pore expansion (15–16 Å, see Anishkin and Sukharev, 2004). With an amplification coefficient of 1.05 (trajectories shown in red), the barrel expanded more decisively, occupying the ranges of tilts between 12° and 45° and diameters up to 18 Å, where the expansion in some of the trajectories stopped and reversed to contraction. Expanding beyond this point, the structure usually became unstable and tended to expand infinitely as illustrated by the trajectories continuing to the right.
Note that the expansion trajectories represent self-permitted pathways for the protein motion but do not include energy terms for hydration and protein–lipid interactions and thus cannot predict limits for expansion reliably. To find a meaningful limit for the barrel expansion we compared the estimated conductance and in-plane expansion area for extrapolated conformations with the experimentally determined values. When stimulated with triangular ramps of pressure, WT MscS displays a substantial hysteresis with the slopes of Po–tension dependencies corresponding to 18 and 9 nm2 on the activation (upward) and closure (downward) branches of the trace. To avoid this ambiguity we used symmetric activation curves for the A98S mutant with no hysteresis and both slopes corresponding to ∼12 nm2 (unpublished data).
Fig. 2 D shows the pore conductance (G) estimated in Hall's approximation (see Materials and methods) plotted against the total in-plane expansion (ΔA) for nine frames from one representative trajectory. The conductance–expansion trajectory crosses the rectangle that represents the range of experimental parameters determined with ∼10% accuracy (ΔA = 12 ± 1.2 nm2 and G = 1 ± 0.1 nS). The expansion trajectory passes through the middle of the expansion–conductance box in Fig. 2 D, which corresponds to a conformation of the TM3 barrel characterized by 30° tilts and pore diameters near 15 Å. Out of ∼20 extrapolated conformations that satisfied both conductance and expansion, 10 had very similar tilt and diameter, clustering tightly within the box in Fig. 2 C. These conformations, shown by thin lines in Fig. 2 B, deviate no more than 2.3 Å (rmsd) from the average position of the helical backbone. An extrapolated conformation closest to the average position (1.2 Å rmsd) was chosen as a “consensus” (shown by thick wire in Fig. 2 B). Examination of the entire transmembrane domain indicated that the consensus model (presented in Fig. 1 D as a candidate open-state model) was characterized with a gap-free packing of the TM1–TM2 helical pairs along the central TM3 barrel and the lowest distortion of the secondary structure compared with the crystals. The increase in the TM3a tilt angle between the closed and candidate open state is close to 10°.
How accurate are structural predictions obtained with the extrapolated motion? At the end of each extrapolation cycle the minimized structure fully satisfies near-equilibrium bond lengths and angles and excludes steric clashes of VdW surfaces according to the CHARMM force field. Since the magnitude of displacements introduced at each step is generally small (<1 Å), the transformation preserves the secondary structure and the repeated cycles of minimizations and short relaxing simulations produce pseudo-continuous transitions along the energy valleys defined by the protein. Computed in vacuum, these energies do not take into account the barriers and wells that could be imposed by the aqueous or lipid environments. Therefore, even though the candidate open-state conformation may satisfy the experimental conductance and barrel expansion, and the chosen segment of the extrapolated trajectory may well approximate the opening transition, the purposeful incompleteness of the extrapolation system permitting computational speed still requires refinements and tests in all-atom equilibrium simulations in an explicit medium.
All-Atom Simulations of the Open Barrel in the Fully Hydrated POPC Bilayer
The model representing the consensus of extrapolated open-like conformations (Fig. 1 C) was embedded in the fully hydrated lipid bilayer, minimized, equilibrated, and simulated under membrane tension of 10 dyn/cm in unrestrained mode for the total of 20 ns. Fig. 3 shows the lipid-embedded open state, presented alongside the closed state previously equilibrated under similar conditions (Anishkin et al., 2008). In the open state, the peripheral and inner helices maintained tight association throughout the simulation and the polarity pattern of the lipid-facing surface matches that of the bilayer. Observation of the system subjected to moderate tensions (10 dyn/cm) for the first 4 ns indicated an initial evolution toward stabilization of the water-filled conformation. The main result was a gradual rotation of TM3 by 53° that largely removed L105 and L109 side chains from the pore and packed them against the side chains of L111 and L115 (Fig. 3, D and E). Simultaneously Gly101 and Gly104 on all seven TM3a helices became exposed into the lumen, thus making the pore wider and more stable due to favorable hydration. Most of the rotation occurred in the TM3a segment of the helix as a result of additional straightening of TM3 around G113. The “latching” of the pairs of aliphatic side chains appears to stabilize the open conformation of the barrel, better defining the size of the pore. Indeed, if the barrel was narrower, the side chains of L105 and L109 would clash and would have to turn back into the pore, thus forcing the helix to rotate back into its initial position, increasing exposure of the hydrophobic surface to water and retracting the polar glycines from the lumen. On the other hand, in an overly expanded barrel, the leucine pairs would be disjoined, the decreased VdW interactions between the aliphatic chains would also increase exposure of their hydrophobic surfaces, and the TM3 helical assembly would produce solvent-permeable gaps. Thus the predicted leucine–leucine interactions impose limits on the TM3 barrel expansion on the cytoplasmic side; however the analysis of multiple extrapolated trajectories indicated that the amount of expansion at the more flexible periplasmic rim can vary in a wider range without compromising the integrity of the wall. We also observed that during the entire 20 ns of unrestrained simulations of the open model at 10 dyne/cm, the previously reconstructed TM2–TM3 interhelical contacts (Anishkin et al., 2008) remained stable. This buried contact formed primarily by the residues V65, F68, L69, and L72 on TM2 apposing L111 and L115 on TM3 appeared to be strong enough to transmit force from the peripheral helices to the core and drive the gate open. Preliminary data from experiments and steered MD simulations suggests that the TM2–TM3 interaction is critical for channel activation and can be interrupted during inactivation (unpublished data).
Simulated under low (10 dyne/cm) tension, the open barrel had a tendency to gradually close. After ∼10 ns we observed the beginning of closure/inactivation associated with protrusion of L105 and L109 side chains into the pore, and the beginning of kink formation near G113. These tendencies are in good correspondence with the results that G113 kinks are characteristic of the inactivated state (Akitake et al., 2007). The process of closing was slow on the nanosecond timescale with an average rate of the decrease in pore radius ∼0.05 Å/ns, so the projection of the observed time course estimated that an arrival to a 6-Å (fully closed) pore would take ∼100 ns. The reason for closing under such tensions can be an overcompressed state of the POPC bilayer in CHARMM (Gullingsrud and Schulten, 2004). Additional simulations of the transmembrane part of MscS conducted at 50 dyne/cm for 16 ns revealed a stably open barrel showing no tendency to close (unpublished data). At this tension, the protein area has stabilized only after 8 ns, predicting the expansion area of ∼16 nm2 compared with the unstressed resting state.
Electrostatic Maps and Explicit Simulations of Conductance through the Transmembrane Barrel
The presence of ions throughout the pore in the expanded conformation indicated a conductive state of the channel. To obtain sufficient statistics of ion permeation events, we increased the concentration of KCl in the simulation cell to 1.5 M. An equilibration of the open barrel with high salt at zero transmembrane voltage for 4 ns resulted in an almost uniform distribution of both K+ and Cl− ions throughout the pore with the exception of the ring of K169 residues where the concentration of Cl− was considerably higher (Fig. 4 B). Upon gate opening, the K169 ring became the narrowest part of the pore (z position = −37 Å). Ions were also present at lower density in the gate and outer chamber (−30 < z < −5 Å), likely due to unfavorable dielectric environment inside the hydrophobic part of the pore. A marked domination of Cl− in the constriction and slight enrichment in the outer chamber, as well as slightly higher concentration of K+ near the outer vestibule (z = +15A) can be due to the presence of fixed charges, K169, R88, D8, D3, and E2 (Fig. 4 A). Note that the profile of the ionic concentrations was calculated for a narrow (4 by 4 Å) column positioned at the pore axis for uniform treatment of the constriction region in the closed and open structures. The periplasmic vestibule of the pore bearing the negative charges is considerably wider (∼30 Å in diameter). In the explicit medium, the negative charges on the pore wall are effectively screened by water and ions within 1 nm distance from the surface. Therefore, the ions scored in the column amount for a very modest increase of the K+ concentration while the bare protein features moderately negative potential in that region. The exact electrostatics of the N-terminal region in our models should be taken with caution since the positions of E2, D3, and D8 were derived not from the crystal structure but predicted de novo by Rosetta (Chivian et al., 2005) and adjusted in extrapolated motions and equilibrium simulations (Anishkin et al., 2008). In the course of pore expansion, R88s facing the pore lumen in the original crystal structure retracted into the protein wall, leaving the electrostatic potential in the lumen at this level close to zero. Fig. 4 illustrates the distribution of the electrostatic potential around the bare protein in vacuum with all charged residues in their default ionization state (Fig. 4 A) and the potential in the complete all-atom simulation cell with lipids, water and ions after equilibration (Fig. 4 C). Without the medium, the channel essentially represents a giant dipole with positive potential concentrated near K169 and a smaller negative potential near the N terminus bearing acidic residues. Previous analysis of MscS electrostatics in the crystal structure and preexpanded conformations (Sotomayor et al., 2006) indicated the ring of R88 as the region of highest positive potential imparting strong anionic selectivity to the current. In our model R88 does not contribute to the pore electrostatics substantially.
The presence of lipid, water, and ions largely equalizes the electrical potential in the aqueous solution, strongly positive potential remains only inside the membrane and protein (Fig. 4, compare A and C). Note that the color scale in C represents much narrower range of potentials than in A. Moderate negative potential corresponds to the areas of increased Cl− concentration near K169, and in the center of the pore around the midplane of the membrane, where local deviations from the complete charge balance in the electrolyte can be seen. Comparisons of the separate electrostatic maps for lipids, ions, and water (unpublished data) indicated that most of the screening of the structural charges is done by ions, and only a small fraction of the field is compensated by water dipoles and the lipids.
Application of moderate, near-physiological voltages (±10 to ± 200 mV, denoted as periplasm vs. cytoplasm) resulted in a steady flow of ions through the channel. At voltages <100 mV, the onset of steady flow was delayed by 2–3 ns (unpublished data), whereas this transient lag period was shorter at higher voltages. The two top panels in Fig. 5 illustrate fluctuations of ion density in the pore constriction in the course of 16-ns simulations at ±100 mV. At any given moment, numbers of K+ and Cl− ions in this region generally match well, satisfying average electroneutrality (Poisson-Boltzmann behavior); however, the total numbers of ions fluctuate with a variance of 30% of the mean, illustrating the thermally driven dynamics of electrolyte density on a nanosecond scale.
The integrals of ion passage events counted over time are represented by near-linear curves (Fig. 5, C and D) reflecting steady permeation rates. At both polarities, the fraction of charge (Q) carried by chloride ions exceeds that of potassium. The difference is more pronounced at −100 mV. Shorter (8 ns) simulations performed at +200 and −200 mV (bottom panels) gave qualitatively similar pictures with a substantially more suppressed K+ current at −200 mV. We noticed that the statistics of permeation events collected at voltages <50 mV was prone to large fluctuations and the values of transferred charge per unit time could not converge in the course of 8 ns, indicating that these voltages require longer simulations. Cumulative data obtained at ±100, ±150, and ±200 mV are presented in the form of I-V curves in Fig. 6 A along with selected experimental data. The inflection shows a slight rectification in the region from −100 to −150 mV, which is then overcome at −200 mV. The slope of the I-V curve near 0 mV for the transmembrane barrel predicts 4.3 nS, which was simulated with 1.5 M KCl solution, producing simulated specific conductivity of 109.6 mS/cm in “benchmark” runs. When adjusted for the presence of the cage contributing 38% of total resistance (see Materials and methods), the simulated I-V curve essentially coincides with the experimental I-V curve measured in symmetrical buffer of 110.4 mS/cm (see below).
We should specify that estimations of cage resistance were done on the crystal conformation and were assumed to be voltage independent, although it is known that the cage is dynamic (Koprowski and Kubalski, 2003; Grajkowski et al., 2005) and the side windows can potentially be wider. Our estimation of resistance is somewhat higher than the experimentally observed 15–20% difference in conductances between WT MscS and the Δ266–286 MscS truncation mutant (Edwards et al., 2007). This mutant lacks a β-barrel at the bottom of the cage which apparently destabilizes the cage (Schumann et al., 2004) and makes it more conductive. The configuration of conductive pathways in the truncated cage, however, is unknown.
Analyzing the ion permeation trajectories we also scored the numbers of water molecules carried with ions. We observed that water follows the dominant flux of chloride ions with water/Cl− ratio of ∼8:1, which is slightly higher than the hydration number for chloride. In our benchmark simulation of the bulk 1500 mM KCl solution, 6.9 waters were in the first hydration shell around Cl−. The observations of individual ionic trajectories in the pore indicated that although the instantaneous concentrations of K+ and Cl− closely match each other (everywhere except in the constriction), the cations are often “swept” in an opposite direction by the counterflux of Cl− and water, thus the fraction of K+ current is disproportionately lower (Fig. 6 B). The cations are often swept in an opposite direction by the counterflux of Cl− and water, thus the fraction of K+ current is disproportionately lower (Fig. 6 B). These simulations show for the first time the microscopic picture of frictional interaction between two oppositely directed ionic fluxes in a weakly selective pore. A movie illustrating cation sweeping events at −100 mV can be found in the online supplemental material ().
Experimental Measurements of I-V Curves in a Wide Range of Membrane Potentials
The above simulations required transmembrane potentials between 100 and 200 mV to obtain reliable statistics of permeation events. Experimental conductivity measurements on MscS under similar conditions are challenging because at pipette voltages beyond +80 mV the channels tend to slip into subconductive states. Beyond −40 mV MscS does not exhibit full conductance either but inactivates through a series of substates (Vasquez and Perozo, 2004; Akitake et al., 2005). To overcome these limitations, we have applied two different approaches to experimental recording of I-V curves in a wider range of voltage, all illustrated in Fig. 7. The first approach was recording of the opening events in response to 1-s ramps of pressure (Fig. 7 A). The pressure at the end of the ramp was chosen at a subsaturating level just to activate a few channels and avoid excessive mechanical stress of the patch. Repeated ramp stimulations provided ample statistics of full opening events, even though many of them showed only substates. This approach worked reliably in the range of ±100 mV for WT and all tested mutants, and for some of them it was possible to record full openings at voltages up to ±150 mV. I-V curves for single channels measured in 0.2 and 1.5 M KCl buffers are shown in Fig. 7 B. The currents are normalized to the specific bulk conductance of electrolyte and, as a result, presented in units of nV·cm versus voltage. The stretches of curves between −20 and +40 mV essentially coincide, showing that the channel conductance at low voltages is strictly proportional to the bulk conductance of the electrolyte, suggesting that the overall concentrations and mobilities for ions in the pore are linearly related to the bulk values. At more positive (hyperpolarizing) pipette voltages, the normalized current at 1.5 M KCl visibly bends down from that measured at 0.2 M salt, suggesting that the increased voltage changes the conductance of the pore, likely by skewing the steady-state concentrations of charge carriers. At depolarizing (negative pipette) voltages, both curves partially level off, showing moderate rectification, more pronounced at the lower salt concentration.
To extend the I-V curves beyond 100 mV, we applied short pulses of progressively increasing voltage on a preactivated population of channels. For the purpose of avoiding high patch currents, we used uninduced spheroplasts exhibiting <10 channels per patch. To increase the stability of the fully open state, we also examined them at saturating pressures. The typical pressure stimulus, train of voltage pulses, and current trace are presented in Fig. 7 C. The right side of the panel shows the magnified top part of the current response to the last (200 mV) pulse. The current displays a nearly exponential decaying kinetics with the rates increasing with voltage. Exponential fits of these responses allowed us to assess the initial value of the current right after the voltage onset. These values normalized to the number of channels in the patch are plotted in Fig. 7 D for three different concentrations of KCl. The I-V curve obtained in the highest concentration (1.5 M) is closest to being linear, showing only a small inflection. The curve obtained in the lowest concentration (0.2 M) is appreciably rectifying, whereas the curve obtained at 0.5 M shows intermediate behavior. The rectification, however, is overcome at more negative pipette voltages. The raw data obtained with the single-channel technique in the narrower voltage range are plotted on the same graph to show correspondence between the two methods. The data for 1.5 M KCl obtained from the single-channel and ensemble currents are presented in Fig. 6 in comparison with the explicit MD simulation data. A good correspondence between the two datasets was achieved after the simulated curve was readjusted to account for the cage resistance. According to our in silico estimations, removal of the cage might increase the conductivity of MscS by ∼61%. This qualitatively agrees with the recent experimental data obtained in the Δ268–288 truncation mutant (Edwards et al., 2007).
The weaker rectification at higher salt points to the self-screening effect and electrostatic nature of the nonlinearity of the I-V curve. Indeed, a number of fixed charges are present near the permeation path of the transmembrane pore. These are E2, D3, D8, and R88 in the outer vestibule and the prominent positively charged K169 ring near the cytoplasmic end of the transmembrane pore (Fig. 4). The importance of R88 and K169 has been recognized in the previous analysis of the MscS pore and cage electrostatics (Sotomayor et al., 2006, 2007). We tested the charge-neutralizing D3N, D8N, R88Q mutations in the outer vestibule and none of them changed the character of rectification significantly compared with WT (unpublished data). However, the K169Q substitution visibly reduced the asymmetry of the I-V curve compared with WT (Fig. 7 E). For WT, the ratio of conductances measured at positive and negative pipette potentials (between −100 and +100 mV) is 2.9, whereas for K169Q this ratio decreased to 1.8 (i.e., conductance increased by 65%). Apparently, the removal of positive charges from the cytoplasmic vestibule of the pore permits more potassium ions to cross the membrane into the periplasm at negative pipette potentials, whereas in WT there are mostly chlorides. When negative voltage is applied, the chlorides are driven away to the cytoplasm, leaving less charge carriers at the cytoplasmic entrance. The neutralization of K169, however, did not straighten the I-V curve completely, indicating that charges on the cage domain may contribute to the rectification more.
This mechanism of rectification has been analyzed previously by Sotomayor et al. (2006), who concluded that the preexpanded crystal conformation should rectify in the opposite direction compared with the experimentally observed open-state conductance. Note that the addition of the charged N-terminal domain to our model and the outward movement of R88 associated with tilting and rotation of TM3s produced a distribution of electric potential around the pore (Fig. 4) very different from what was reported for the unmodified crystal structure (Sotomayor et al., 2006). The pore does not seem to contribute strongly in rectification and the fact that most of the rectification remained after the neutralization of K169 points to the critical role of the cage domain in defining the directionality of the conductance and selectivity. The selectivity measurements in MscS pore mutants recently published by Edwards et al. (2007) also showed that R88 is not a part of the selectivity mechanism. Consistent with the expanded state of the open pore, the R88S, T93R, G101D, and G113R/D mutations did not affect selectivity substantially, suggesting that opening moves these residues sufficiently far from the center of the pore. The G113R substitution placing extra positive charge in the cytoplasmic vestibule near K169 increased rectification, but the Δ266–286 deletion mutation, which opens a hole at the bottom of the cage (and apparently destabilizes the cage; Schumann et al., 2004), essentially removed rectification (Edwards et al., 2007). The previous computational analysis of cage electrostatics also revealed strong separation of anions and cations between the upper and lower hemispheres (Sotomayor et al., 2006), and the data are consistent with our conclusions of a limited role of the charges immediately surrounding the pore in setting the electrostatics for ion permeation. The intriguing questions of potential involvement of the electroosmotic water flux observed in simulations (Fig. 6 C) in MscS conduction and regulation, and the disappearance of rectification at higher depolarizing voltages (Fig. 7 D) should be addressed in the near future in combined experiments and simulations.
Discussion
First simulations of the transmembrane segment of MscS (Anishkin and Sukharev, 2004) or entire crystal structure (Sotomayor and Schulten, 2004) revealed essentially dehydrated state of the pore constriction and no ion permeation events until the transmembrane voltage was increased to 500–1,500 mV (Spronk et al., 2006; Vasquez et al., 2008). This “dry” state of the narrow crystallographic pore made the macroscopic conductance estimations with HOLE (Smart et al., 1996) or Monte-Carlo simulations (van der Straaten et al., 2005), both implying a continuous aqueous path and near-bulk diffusion coefficients, hardly applicable. It was estimated that the diameter of the constriction should be increased to ∼16 Å to reach a stably hydrated pore satisfying the experimentally observed conductance (Anishkin and Sukharev, 2004). The first attempts to enlarge the pore by applying membrane tension to the splayed peripheral helices just increased the hydration of the constriction and stabilized the barrel from asymmetrical collapse (Sotomayor and Schulten, 2004). Further studies using tension, steering forces, or high transmembrane voltage produced several expanded models (Sotomayor and Schulten, 2004; Sotomayor et al., 2006, 2007; Spronk et al., 2006). Some of these models reached pore diameter of 16 Å and were estimated to conduct at levels close to the experimentally observed, albeit at much higher voltages. Subjected to lateral pressure of lipids from all sides, the peripheral helices slightly decreased their tilts and the TM2–TM3 gaps partially closed. In several of the reported pre-expanded states, the characteristic kink at G113 was observed to straighten in some of the subunits, implicating this type of motion in the gating mechanism (Sotomayor et al., 2006; Spronk et al., 2006).
Either explicit (Spronk et al., 2006; Sotomayor et al., 2007) or Monte-Carlo (Sotomayor et al., 2006; Vora et al., 2006) simulations of expanded states invariably indicated strong Cl− selectivity, with PCl–/PK+ permeability ratios varying between 20 and 100. Experimentally measured PCl–/PK+ ratios were always between 1.2 and 2 (Sukharev et al., 1993; Li et al., 2002; Sukharev, 2002; Edwards et al., 2007). Although the process of relaxation of these expanded conformations at low voltages has not been reported, it appeared that these models were stable and conductive only under high membrane voltage (>500 mV), which implied strongly nonlinear I-V curves. This was in disagreement with the experimentally observed long-lived and flicker-free opening events as well as nearly linear conductance at voltages around zero (Sukharev, 2002; Shapovalov and Lester, 2004; Edwards et al., 2007). These properties of the open state obviously required more profound rearrangements of the crystal structure that could change the electrostatic potential in the lumen and stabilize the conductive pathway against collapse.
To search for permitted conformations, we developed extrapolated motion protocol, based on small extrapolated displacements followed by minimization and relaxation steps. Implementing this method, we used the flexibility of the Tcl scripting language and the power of the NAMD-VMD suite (Humphrey et al., 1996; Phillips et al., 2005). It allowed us to overcome many of the inherent scale limitations of traditional MD (reproducing canonic thermodynamic ensembles) and to explore self-permitted pathways for protein motion. This iterative protocol previously helped us to generate a compact model of the resting state (Anishkin et al., 2008). An alternative model of resting MscS independently deduced from an extensive set of EPR assignments has similar parallel orientations of peripheral helices and the N termini well exposed to the outer leaflet of the lipid bilayer (Vasquez et al., 2008).
The extrapolated motion protocol also permitted us to envision the pathway for barrel expansion. Because the protocol is computationally inexpensive, we were able to generate more than 50 expansion trajectories. The presence of a “random component” during short periods of MD created scatter in these simulations, allowing for more effective exploration of the conformational space. Bundling of the trajectories under such circumstances pointed to a common expansion pathway. The frequently visited regions on the pore diameter–helical tilt diagram (Fig. 2 C) identified plausible closed and open conformations. The consensus model satisfying the experimental parameters was then subjected to further refinement. As a modeling approach, the extrapolated motion protocol appears to be less “subjective” than modeling “by hand.” The fact that it is capable of generating a variety of similar conformations is more consistent with the modern view that each of the functional states, closed or open, when subjected to thermal fluctuations is in fact represented by a set of different conformations characterized by similar geometry of the gate.
To make the extrapolated motion procedure more efficient, transformations were done in vacuum, purposely disregarding the energetic contributions of water and lipids. Due to the incompleteness of the system, we opted to use experimental data to find the limits for the transition and then refine models in all-atom simulations. The experimentally determined in-plane area change (∼12 nm2) and open-state conductance (1 nS in a 39 mS/cm buffer) defined the functionally meaningful limits for extrapolated expansion (Fig. 2 D). We found a set of similar expanded conformations that satisfied both the conductance and area changes at the same time, suggesting that the degree of lateral expansion, thickness of the channel wall, and the packing of the peripheral helices were chosen realistically. This coincidence allowed us to choose a segment of the extrapolated trajectory approximating the transition. At the end of expansion, the TM3 helices invariably straightened, and this structural feature defined the geometry of the open pore. The concerted straightening of helices acting as “collapsible struts” under tension appeared to be the major factor stabilizing the open state. Indeed, the mutations that increase the helical propensity of TM3s at two “hinge” sites, G121 and G113, have recently been shown to increase the stability of the open state, whereas increased flexibility at these sites increased the rates of desensitization (G121) or inactivation (G113) (Akitake et al., 2007).
We should remember that our current models for the open and closed states of MscS have been developed from the originally deposited 1MXM crystal structure of MscS (Bass et al., 2002). This structure was recently replaced by the refined 2OAU structure (Steinbacher et al., 2007). If the new structure was used as a starting point for our modeling of the gating transition, it would probably have a negligibly small effect on the outcome. Indeed, the revised crystallographic model has several registry shifts in the cytoplasmic cage, whereas our simulated transitions that represent gating occurred primarily in the transmembrane part. The most noticeable change in the transmembrane part of the 2OAU structure compared with the previous one is the relocation of the R88 side chain further away from the pore axis, which decreased the positive potential in the pore at that level. However, a similar retraction of R88, albeit of larger amplitude, occurred spontaneously in the course of our 1MXM crystal structure extrapolated transitions, leading to even more substantial lowering of the potential in the pore. Other minor changes in the refined atomic coordinates of the transmembrane domain are not expected to bias our models consistently since the conformational changes experienced by the crystal structure in extrapolated-motion transformations are of much larger scale.
The candidate model for the open state derived from multiple extrapolation trials was stable in extended all-atom MD simulations. During the unrestrained MD equilibration, a rotation of helices was observed that made the pore more hydrophilic (Fig. 3, C–E). L105 and L109 side chains swung out of the pore, joining L111 and L115 on the adjacent helix. At the same time, glycines 101, 104, and 108 became exposed to the pore lumen, considerably increasing the polarity of the pore lining. This “hydrophilic” stabilization of the open state is highly consistent with the previous results by Edwards et al. (2005) that alanine substitutions for these glycines make the channel “stiffer,” i.e., opened at higher tension, whereas more polar A106G (and A106S) substitutions lead to easy opened “gain-of-function” channels. It appears that the open state of MscS is doubly stabilized through the favorable exposure of glycines and association of leucines.
Both the extrapolated and MD-equilibrated structures demonstrated that straightening, tilting, and outward motion of TM3s define the scale of the transition and the geometry of the open pore. Supplemented with the Rosetta-predicted N-terminal segment, the stably hydrated pore predicts a slight positive potential in the lumen. The time-averaged electrostatic maps predict that the open pore at low potentials should have weak anionic selectivity consistent with all experimental observations. We performed explicit simulations of open MscS conductance in a manner previously reported by Aksimentiev and Schulten (2005) and found that the regimen of ion transport is close to that in the bulk. At 1.5 M salt concentration and voltages between 100 and 200 mV, the fraction of current carried by Cl− ions was higher than that of K (Fig. 5) not only because more Cl− is in the pore, but also because Cl− carries more water and the counterflux of K ions is swept by the electroosmotic water stream created by the Cl− current. To our knowledge, our explicit conductance simulations for the first time illustrated interaction of two ionic fluxes inside an essentially nonselective water-filled pore.
After an adjustment of the simulated channel conductance for the presence of the cytoplasmic cage, we observed a surprisingly good correspondence between the simulated conductance, the inflected shape of I-V curves, and the experimentally measured characteristics under similar ionic conditions (1.5 M KCl) (Fig. 6). A stronger rectification was experimentally observed in low-salt (0.2 M) buffers (Fig. 7), but simulations under such conditions would require considerably longer times to conclude on conductance and selectivity and thus were not feasible. We should also emphasize that in the current configuration our simulation system may not be able to reproduce the rectification and selectivity precisely because the cage domain is absent.
How realistic are explicit conductance simulations in general? It is known that bulk properties of physiological electrolytes can be well reproduced in MD simulations; however, the accuracy varies nonmonotonously with salt concentration (Aksimentiev and Schulten, 2005; Sotomayor et al., 2007). Correspondence to experimental data is very good at low concentrations of KCl, at intermediate concentrations (∼100 mM KCl) MD overestimates conductivity by ∼25%, and then simulated specific conductivity declines and approaches the experimental values again at ∼1000 mM concentrations. To enhance the computed and recorded ionic currents, our simulations and experiments were performed at even higher concentration of KCl (1.5 M). The specific conductivity of 109.6 mS/cm obtained in our benchmark simulations underestimates the experimental value of 156.3 mS/cm for 1.5 M KCl solution, thus continuing the previously reported tendency to decline. However, our experimental data with MscS were obtained not in pure KCl, but in a buffer (see Materials and methods) of 110.4 mS/cm–specific conductivity. This value is fairly close to the simulated conductivity, therefore a rather insignificant scaling factor of 1.006 was applied to simulated data before comparison with experimental values. The simulated ratio of the chloride to potassium currents contributing to the conductivity of the bulk electrolyte was 1.06—very close to the experimental ratio of equivalent conductivities (Lide, 2007), suggesting that the error in simulated Cl−/K+ current ratio in MscS introduced by the force field inaccuracies is very small.
Having plausible models for the resting and opens states of MscS, one may envision a transition driven by membrane tension applied to the peripheral helices only. The tight associations between TM1–TM2 pairs and TM3, especially in their cytoplasmic regions, appear to be sufficiently strong to provide a mechanical link between the lipid-facing periphery and the gate. Since both models equilibrated at low tension have essentially cylindrical outer wall facing the lipid (Fig. 3), the transition driven by membrane tension may appear as a uniform or slightly hourglass-like expansion. At this moment we cannot exclude that a more realistic pathway for opening may proceed through a funnel-like conformation. Indeed, the periplasmic part of the barrel appears to be “softer” and may comply more easily with tension, whereas the dehydrated gate, acting as a hydrophobic sphincter located below the level of phospholipid headgroups may impart stiffness to the cytoplasmic end. Cone-shaped intermediate states may expose hydrophilic groups in the pore and thus promote hydration of the gate, resulting in full opening and transformation of the barrel into a more cylindrical shape depicted in Fig. 3 B.
If the straight conformation of theTM3 helices defines the open state, how does the channel close? As was shown previously (Akitake et al., 2007), the balance between lateral tension, hydration forces (sphincter-like action of the gate-keeping leucine rings), and TM3 flexibility at two separate glycines, G113 and G121, defines the stability of the straight conformation and correspondingly the overall life time of the open state. TM3s thus act as collapsible “struts” making the frame of the open pore. While the exit rate from the open state directly depends on the flexibility of TM3s near G121 promoting reformation of the hydrophobic gate, the probability of inactivation relies on helix bending near G113 (Akitake et al., 2007). The distance between the C ends of TM3s resting on the roof of the cytoplasmic cage thus become an important parameter, and therefore the dynamics of the cytoplasmic cage is predicted to have a strong regulatory role (Koprowski and Kubalski, 2003; Edwards et al., 2004; Schumann et al., 2004; Grajkowski et al., 2005).
Although in the extrapolated-motion simulations all domains of the MscS protein were free to move, the rigidly folded cage generally experienced much smaller changes compared with the transmembrane domain. In the course of opening transition, an average displacement in the region of the cage proximal to the transmembrane domain was only ∼3 Å. For the subsequent all-atom simulations of the transmembrane region in the explicit bilayer, the cage was largely removed and the truncation points were harmonically restrained at the modeled positions to prevent artifactual drift. Note that these all-atom simulations were designed to relax and explore certain stable states in the explicit medium rather than force any transitions. We presume that in the course of the gating cycle, the cage may experience conformational transitions of a scale much larger than 3 Å, but these processes may be slower than the dynamics of the transmembrane barrel. Further studies and better understanding of the cage dynamics are especially important in the light of experimental studies (Koprowski and Kubalski, 2003; Edwards et al., 2004; Schumann et al., 2004; Grajkowski et al., 2005), suggesting that certain modes of cage motions may be decisive in regulating the gating of the pore.
While closure, desensitization, and inactivation will be subjects of future research, the beginning of this process was already observed in our extended all-atom simulations under low tension. We saw a tendency for slow closing manifested as turning helices back, reformation of kinks, and relaxation of the height of the barrel. The return of aliphatic side chains of L105 and L109 back to the lumen appears to be the beginning of gate reformation that should lead to a pinching off the column of water and initiation of dewetting (Anishkin and Sukharev, 2004). We should note that the CHARMM forcefield produces a noticeably overcompressed POPC bilayer and might require 30–40 dyne/cm (Gullingsrud and Schulten, 2004) just to maintain the area per lipid corresponding to the experimental values (Klauda et al., 2006; Kucerka et al., 2006). Thus tension set at 10 dyn/cm in our equilibrium simulations may actually exert a considerable positive compressing action at which the open conformation cannot be maintained infinitely long. We are currently exploring more stable open conformations observed under 50 dyne/cm, which may become next modeled approximations for the open MscS structure.
© 2008 Anishkin et al. This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.jgp.org/misc/terms.shtml). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 3.0 Unported license, as described at http://creativecommons.org/licenses/by-nc-sa/3.0/).
Abbreviations used in this paper: EPR, electron paramagnetic resonance; MD, molecular dynamics; MscS, mechanosensitive channels of small conductance; WT, wild-type.
Acknowledgments
The authors thank Dr. H. Robert Guy for extensive discussions and the set of first models of open MscS that he provided. These models stimulated the progress of this work. The authors also thank Miriam Boer and Naili Liu for assistance with MscS mutants, spheroplasts, and preparation of the manuscript.
The project is supported by National Institutes of Health grant R01GM075225 to S. Sukharev.
Lawrence G. Palmer served as editor.