|
|
|
|
|
|
|
||
ARTICLE |
Correspondence to Sergei Sukharev: sukharev{at}umd.edu
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONREFERENCES |
|---|
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.
© 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/).
 |
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONREFERENCES |
|---|
) 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 |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONREFERENCES |
|---|
) 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 (http://www.jgp.org/cgi/content/full/jgp.200810000/DC1).
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONREFERENCES |
|---|
) 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.
|
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.
|
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).
|
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.
|
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.
|
|
; 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 (http://www.jgp.org/cgi/content/full/jgp.200810000/DC1).
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.
|
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 |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSRESULTSDISCUSSIONREFERENCES |
|---|
) 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
