| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
* Department of Physics, and Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana-Champaign, Illinois
Correspondence: Address reprint requests to Klaus Schulten, E-mail: kschulte{at}ks.uiuc.edu.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTS AND DISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTS AND DISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
–5). In bacteria, mechanosensitive channels have been proposed to act as safety valves preventing the cell from bursting upon osmotic shock (2,6–8). Two bacterial mechanosensitive channels have been cloned and crystallized: the mechanosensitive channel of large conductance, MscL (Mycobacterium tuberculosis) (9,10) and the mechanosensitive channel of small conductance MscS (Escherichia coli) (7,11). While the gating mechanism for MscL is fairly well understood and multiple experimental and theoretical studies have been performed on it (12–19), the voltage-modulated and stress-dependent gating mechanism of MscS remains controversial (11,20–25).
The crystal structure of MscS displays three transmembrane helices (labeled TM1, TM2, TM3A, and TM3B) forming a pore of radius 
3.3 Å, originally suggested to represent the open state of the channel (11) with a conductance of 1 nS and a slight anionic selectivity as determined in patch-clamp experiments (2,23,26). However, subsequent computational studies suggested that the crystal structure is not completely open (20). In addition, an experimental characterization of MscS in spheroplasts (23) revealed voltage-independent activation and strong voltage-dependent inactivation (first reported in (27)). Moreover, the channel activity was found to depend strongly on the rate of pressure applied and the crystal conformation was suggested to represent an inactive state (23). A conductance substate for MscS has been recently discovered through an optimized patch-clamp setup (28) and two different closed conformations of MscS have been proposed as well. Indeed, molecular dynamics simulations of the entire MscS structure revealed an asymmetric closure when the channel was relaxed in a POPC membrane for several nanoseconds (21) (also reported in (29)). On the other hand, a symmetric closed model was proposed based on Monte Carlo simulations and site-directed mutagenesis studies that characterize the role of the glycine-rich TM3A helix in MscS gating (22).
In addition to the transmembrane pore, MscS features a large cytoplasmic domain of unclear function with seven side openings and one distal pore. Recent studies have indirectly implicated this balloon-shaped domain in gating (30–32), desensitization and stability (33,34), as well as transport and selectivity (21). Specific interactions between the cytoplasmic domain and the transmembrane domain were implicated in MscS gating (21).
While the various studies mentioned above have shed light on architectural and functional properties of MscS, many questions remain unanswered and new ones have arisen: What is the actual physiological role of MscS, i.e., is MscS only a safety valve? Is the published structure significantly distorted due to crystallization conditions? Does the crystal structure show a so-called conductance substate, not a fully open state? What residues of MscS are relevant for its gating? What is the role of MscS's large cytoplasmic domain?
Molecular dynamics (MD) simulations provide a necessary complement to the experimental effort toward an answer to these questions. However, realistic simulations of a complete channel embedded in a membrane environment present a formidable problem (35–38). Relevant cellular dimensions range from 10 nm to 1 mm, while specific channel behaviors, such as switching or selectivity, are also controlled by sub-nanometer functional subunits of the protein. In addition, strong short-range ion-ion and ion-protein interactions require ion dynamics to be resolved on a femtosecond timescale while reliable estimates of current require simulations lasting several 100 ns.
Although large-scale MD simulations provide detailed insight into the mechanisms of channel structure, gating, and selectivity (see (39) and references within), reliable estimates of channel conductance are barely possible with resources currently available (24,40). Continuum models based on Poisson-Nernst-Planck (Drift-Diffusion) theory have recently been proposed as a computationally efficient alternative means for computing channel conductances (41–54); however, the results of these models must be interpreted with caution, since the continuum paradigm deals with average densities rather than discrete particles of finite size, a premise that is problematic in channels only a few Angstroms wide (41,53). To extend simulation times well beyond the present limits of MD simulations, while still retaining the essential particle nature of the ions, a coarse-grained particle approach is needed. The Biology Monte Carlo (BioMOCA) ion channel simulator has been developed for this purpose (55), based on the Boltzmann transport Monte Carlo (56) and particle-particle-particle-mesh methodologies (57).
In this work we study the electrostatic properties of the crystal and various open conformations of MscS using a combined approach in which conformations obtained from MD are used as starting points for BioMOCA simulations, thus reaching timescales of hundreds of nanoseconds. The simulations revealed that ions permeate MscS through side openings of the cytoplasmic domain rather than the distal cytoplasmic pore; time-averaged electrostatic potential maps and ion densities reflect a strong localized separation of anions and cations inside and around the cytoplasmic domain and the transmembrane pore. We found that ionic currents computed for MscS conformations similar to those of the available crystal structure (11) were too small, while currents computed for wider conformations reached through steered MD simulations reproduce the experimentally determined values. However, ionic currents were found to be mainly driven by Cl– ions, indicating a high selectivity for anions over cations, in relative disagreement with present experimental results that suggest slight anionic selectivity. Rectification at high voltages (>50 mV) was found to be in disagreement with experiments as well.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTS AND DISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
)) and a coarse-grained method (BioMOCA (55)). In this section we present first the protocols of all-atom molecular dynamics simulations used to generate MscS conformations. Simplified models of channels with different charge distributions are described as well. Then, details and parameters used in BioMOCA simulations of these conformations and models are described along with the diagnostics utilized.
MscS conformations from molecular dynamics simulations
Four different conformations of MscS obtained from molecular dynamics simulations were utilized for BioMOCA simulations (see below). In all molecular dynamics simulations, the crystal structure of MscS (Protein Data Bank code 1MXM (11)) was embedded in a membrane bilayer of 299 POPC lipids solvated with over 50,000 explicit water molecules (TIP3P) and 200 mM of KCl. The simulated system included altogether 224,340 atoms. The first conformation, referred to as crystal, was obtained after 4.5 ns of dynamics with the backbone atoms of the protein restrained to their original positions in the crystal (k = 1 Kcal/mol/Å2), but all other atoms free to equilibrate otherwise (Fig. 1, A and B; and Supplementary Material's Fig. 7 A). This crystal conformation features a pore slightly narrower than that of the crystal itself (rmin = 2.8 Å compared to 3.3 Å; all pore radii were computed using HOLE (59) and default AMBER van der Waals radii for consistency with previous work, see (21) for details). The second conformation, labeled open1, was obtained after 4.6 ns of dynamics in which restraints to backbone atoms were gradually released while a surface tension of 20 dyne/cm was applied to the system (Fig. 1 C and Supplementary Material's Fig. 7 B). Details of the corresponding molecular dynamics simulations leading to crystal and open1 conformations can be found in Sotomayor and Schulten (21).
|
)). During the first MD simulation, a surface tension of 50 dyne/cm was applied to the system in the plane of the membrane (x,y-plane; see (13,21)) along with a constant force of 10 pN acting on each of the C
atoms of residues Tyr27, Val29, Arg46, Arg54, Arg59, Arg74, and Arg128. Forces were applied radially, pointing outwards from the center of the channel, in the x,y-plane, for 12 ns (Supplementary Material's Fig. 7, E and F), after which the forces were increased to 40 pN and the system was simulated for an additional 4 ns, leading to conformation open2 (Fig. 1 E and Supplementary Material's Fig. 7 C). The position of the center of the channel, used to define the direction of the forces applied, was updated at simulation times 4.5, 6.2, 9, 10.8, 12, 13.7, and 15 ns. The magnitude of the forces is comparable to those used previously to open MscL under similar conditions (18).
The second MD simulation, leading to conformation open3 (Fig. 1 F), consisted of 1 ns of dynamics with an applied surface tension of 20 dyne/cm, followed by 1 ns of dynamics with an applied surface tension of 40 dyne/cm. Then, the surface tension was turned off, a constant volume protocol was assumed, and a force of 10 pN was applied to each of the C atoms of residues 96–113 (TM3B helix) over 1 ns (Supplementary Material's Fig. 7, G and H). Forces were then increased to 20 pN and the simulation continued for an additional 2 ns. Forces were applied radially in the x,y-plane pointing outwards from the center of the channel to induce contacts between TM1–TM2 and TM3, as loss of these contacts has been suggested to cause inactivation (23). All external forces were applied through the NAMD/Tcl interface (58).
Conformations crystal, open1, open2, and open3 were taken from the last frame of the corresponding simulations described above and mapped to a rectilinear grid (see below) including protein and lipids, but not solvent and ions. Explicit lipids provide a more realistic, irregular membrane, leaving some space between crevices formed at the cytoplasmic side of the protein for ions to visit. Mobile water molecules within the bilayer are known to contribute to the membrane potential (60), but it is not clear how to best include them in the present context. A fourth conformation labeled crystalNC was obtained from crystal by excluding the cytoplasmic domain of MscS (residues 129–280). The mapping of each of the five conformations mentioned above (crystal, open1, open2, open3 and crystalNC) was performed using the standard set of partial charges of CHARMM (61) along with a standard set of radii for atoms used in VMD (see Supplementary Material's Table 3). We decided not to use the standard set of CHARMM radii, since they have not been designed for continuum simulations (see (51,62–64) for examples in which atomic radii have been optimized for continuum simulations). Indeed, the combined sets of CHARMM charges and VMD radii are very similar to those of parameters for solvation energy (i.e., PARSE (65)), previously used for estimation of currents with continuum models. However, to assess the effect of this choice, a sixth conformation, labeled crystalCR, was constructed identical to crystal but using the standard set of CHARMM radii during the mapping to a rectilinear grid.
Schematic models
Four schematic models, referred to as mod1, mod2, mod3, and mod4, were built containing a slab of membrane of thickness 50 Å with a cylindrical pore of radius 5 Å. Each model featured three rings of seven atoms each (ratom = 1.5 Å, unless otherwise stated). The first two rings were located at the cytoplasmic and periplasmic ends of the pore. The third ring was located 15 Å away from the membrane in the cytoplasmic bath. All rings were set neutral for mod1, whereas the rings at both ends of the channel were set positive and the third ring was set neutral for mod2. Model mod3 featured positively charged rings at the ends of the channel while its third ring was negatively charged. The fourth model, mod4, designed as a control model, included two negatively charged rings at the ends of the channel and a third, neutral, ring.
Three additional schematic models (mod5, mod6, and mod7) with pores of radius 10 Å were prepared for BioMOCA simulations as well. Model mod5 was decorated with three neutral rings of seven atoms each, as described above. Model mod6 featured positively charged rings at the ends of the channel and a third, neutral, ring located 15 Å away from the membrane in the cytoplasmic bath. Finally, model mod7 featured only two, neutral, rings at the ends of the channel.
Coarse-grained, particle-based simulations
Calculations of the electrostatic field around MscS and the current through its pore were carried out using the program BioMOCA; since the program has been described extensively elsewhere (55), we provide only a brief description here. In contrast to MD simulations, where the motion of every particle is explicitly resolved, in BioMOCA the protein, lipid, and water molecules are treated as static continuum dielectric media and only the trajectories of the mobile ions in solution are computed. This reduces the number of particles from O(105) to O(102), allowing much longer simulation times than what can currently be achieved using MD.
The molecular structures were taken from snapshots of molecular dynamics simulations or schematic models as described above, and mapped onto a rectilinear grid using the scheme implemented in the adaptive Poisson Boltzmann Solver (66), thus partitioning the simulation domain into three contiguous regions: protein, lipid, and electrolyte. For all of the simulations discussed here, a uniform grid with a mesh spacing of 1.5 Å was used. Probe radii of 1.38 Å and 4 Å were utilized to determine the protein and lipid boundaries, respectively.
Ion trajectories were computed in real space as sequences of free flights which are randomly interrupted by scattering events. This is equivalent to solving (with the Monte Carlo methodology) the Boltzmann Transport equation, widely used in solid-state physics to describe transport of electrons (57,67). The flight times tf between scattering events were generated statistically from the total scattering rate for all scattering processes 
according to
 | (1) |
depends on the particle momentum 
which in turn will depend on the ion's local environment, and constitutes the sum of individual scattering rates for various scattering processes (in the context of solid-state physics, any impurities, crystal defects, or thermal vibrations of the crystal lattice would be classified as different scattering processes). When an ion reaches the end of its free flight, another random number is used to select which type of scattering event took place and the final state of the ion is selected according to the particular scattering model implemented. In regions of restricted volume, such as the interior of the channel pore, ion-water scattering processes are likely to be significantly different from bulk electrolyte processes and additional scattering processes may prevail. BioMOCA has been, therefore, equipped to handle multiple scattering processes, and space- and momentum-dependent scattering rates. However, in this study we have assumed a single scattering mechanism that thermalizes the ions, with a constant scattering rate 
i for each species, linked to an ion's mass mi, and diffusion coefficient Di in water,
 | (2) |
), yielding scattering rates of K = 3.26e13 s–1 and Cl = 3.46e13 s–1. At the end of each free flight the ion's velocity was reselected randomly from a Maxwellian distribution. The assumption of a single constant scattering rate yields ion transport that is equivalent to Brownian dynamics.
The use of a bulk-like, uniform diffusion coefficient for K+ and Cl– will likely result in the overestimation of computed currents for the narrower structures. Indeed, previous studies using all-atom MD simulations and hydrodynamic theory have shown that ion mobility is highly reduced in narrow channels (51,54,69). This effect seems to be similar for both K+ and Cl– ions and, therefore, should not affect selectivity (70). While the reduction in mobility should not be dramatic for a fully conductive channel with a large conductance like MscS, future studies should use position-dependent diffusion coefficients computed from all-atom MD simulations for MscS's transmembrane pore and all openings of its cytoplasmic domain.
Ion dynamics and field solution
In between scattering events, the ions moved in the electrolyte according to Newton's laws of motion in response to the local field. Trajectories were synchronously propagated in time and space by integrating the equations of motion using the second-order accurate leapfrog scheme (57). High ion-water scattering rates necessitated a short trajectory integration timestep of 10 fs. The total force on the ion was the sum of a pairwise ion-ion repulsive component, which represents ionic core repulsion, and an electrostatic component, which was obtained self-consistently by solving Poisson's equation,
 | (3) |
is the local electrostatic potential, 
and 
are, respectively, the densities of mobile ions and permanent fixed charge on the protein and lipids, and 
is the local dielectric constant discussed below. Solving Eq. 3 over the entire domain subject to specific Dirichlet boundary conditions provided a simple way to include an applied bias and the effects of image charges induced at dielectric boundaries. The mobile charge density was recomputed from the positions of all ions in the system every 100 ion trajectory time steps, i.e., every 1 ps, as well as each time an ion entered or left the system, i.e., every time the total charge in the system changed. Each ion's charge was mapped to a finite rectilinear grid using the cloud-in-cell scheme (57) to weight the charge to the eight nearest grid points. The same procedure was used to map the fixed partial charges residing on the protein and lipid to the grid at the start of the simulation. Equation 3 was discretized using finite differences and solved using conjugate gradient (71) methods, subject to specific boundary conditions described below. Once the potential was known, the electrostatic field 
at each grid point was obtained using second-order finite differences and the field at the eight nearest grid points was then weighted back to the ion using the same interpolation kernel as was used to map the charge density. Discretization of Poisson's equation leads to an unavoidable truncation of the short-range component of the electrostatic force, which was corrected using the particle-particle-particle-mesh scheme (55,57). A pairwise force which prevents ions from coalescing unphysically was also included as mentioned above. For computational efficiency we implemented a truncated-shifted form of the Lennard-Jones 6-12 potential that retains only the short-range repulsive component (55). We note, however, that the parameters describing pairwise interactions do not take into account the effect of missing water molecules with a finite volume (hydration shells), and have not been validated against known activity coefficients.
The protein and lipid regions were deemed inaccessible to ions by assigning a finite radius, equal to the Pauling ionic radius (72), to each ion. For K+ and Cl– ions these ionic radii are 1.33 Å and 1.81 Å, respectively. If any point within the spherical ion overlapped the protein or membrane boundary during a time-step, the ion was returned to its position at the beginning of the time-step and reflected diffusively. This process was repeated until the ion's final position lay outside the protein/membrane region.
Dielectric coefficients
Perhaps the most challenging aspect of coarse-grained ion channel simulations is the assignment of appropriate values for the protein, lipid, and electrolyte dielectric coefficients, which determine the strength of the electrostatic interaction between two charged species as well as the self-force felt by an ion as it approaches a dielectric boundary. The protein environment responds to external fields in several ways, each with its own relaxation timescale. For models that include the permanent protein dipoles explicitly, but treat the protein and water reorganization implicitly, the dielectric coefficient is hard to define, particularly when ion motion takes place on the same timescale as the protein's response to its presence. In such cases it has been suggested that the value of the protein dielectric coefficient depends on the nature of the interaction: for charge-dipole interactions a value between 4 and 6 is recommended, while for charge-charge interactions the protein dielectric coefficient can be >10 (73). For this work we kept 
l = 2 in the lipid region and used p = 5 in the protein region, unless stated otherwise. The electrolyte bath regions were assigned the dielectric coefficient of bulk water; e.g., w 
80. Note, however, that in most ion channels the pore is very narrow and is often lined with highly charged residues and/or strong permanent dipoles. As a result, water molecules can be highly ordered within the channel pore, restricting their response to external fields (74,75) and the dielectric coefficient ch is likely to be <<80 and anisotropic. Although the transmembrane pore of MscS is wide enough to accommodate a few files of water molecules, all-atom MD simulations of MscS have shown intermittent permeation of water through the pore, or so-called dewetting transitions (20,21), which may also affect the dielectric properties of this region.
Ideally, the dielectric coefficient should decrease gradually from the bulk-like bath regions to the channel pore, but computing electrostatic forces acting on ions in regions surrounding a sharp transition in dielectric coefficient poses additional numerical complications. For the simulations presented here we have adopted the same approach as several other groups (51,64,76–78) and used ch = 80 inside the channel pore as well as in the electrolyte bath regions. However, we note that future studies should further investigate the value of the dielectric constant within the cytoplasmic domain of MscS and along its transmembrane pore.
Boundary conditions
Experimentally, the electrical properties of a single ion channel can be measured by inserting the channel into a lipid bilayer (membrane) separating two baths containing solutions of specific concentrations (79). Electrodes are immersed in the baths to maintain a constant bias voltage across the membrane. To adequately represent these contact regions without the computational overhead of including extremely large baths, we assumed that beyond a Debye length from the membrane surface the average electrostatic potential and ion densities do not vary appreciably, an assumption supported by the results of recent continuum simulations (80). For salt concentrations considered here, this distance is of the order of 10 Å. Therefore, we imposed Dirichlet boundary conditions on the potential at the two domain boundary planes that lie parallel to the membrane, which are located 15 Å from the periplasmic side of the lipid and 15 Å from the distal end of the cytoplasmic domain of the protein. Homogeneous, i.e., zero-field, Neumann boundary conditions were imposed at all other domain boundaries.
Electrolyte baths of fixed concentration were modeled on either side of the membrane as follows: Ions that left the domain through the Dirichlet boundary planes were removed from the simulation. If the ion population for each species in either bath fell below the value representing the desired density, new ions were injected with a Maxwellian velocity to maintain the given density. The newly injected ions were placed randomly within a 4 Å (6.25 Å for simplified models) slab from the Dirichlet boundary planes. It should be noted, however, that ions were never artificially removed from these buffer regions. Ions that attempted to cross the Neumann boundary planes were reflected elastically.
Calculation of potential maps, ion concentrations, and currents
Ionic currents were computed by counting the number of ions that crossed the membrane during the simulation and normalizing by the simulation duration (see (40,81) for alternative methods). For smaller channels like gramicidin this is an unambiguous means of estimating stationary current. For MscS, however, the number of ions that cross the membrane may not be an accurate measure of the actual experimental current since MscS's large cytoplasmic domain may lead to substantial, partial transport across the simulation volume. Nevertheless, trajectory analysis indicated that ions passed through the transmembrane pore relatively quickly (O(100 ps)) and entered or left the cytoplasmic domain through its side openings within the simulation timescale used, i.e., partial transport was rare.
In addition to ionic currents, time-averaged three-dimensional maps of the electrostatic potential and ion concentrations were computed. Every time the Poisson equation was solved, the instantaneous ion density was sampled and added to a running average. One-dimensional average ion concentration 
plots were computed from the final concentration map by integrating the ion density n(i, j, k) in the plane parallel to the Dirichlet, i.e., x,y-boundary planes, and normalizing by the total ion-accessible volume in that plane. For this purpose we employed the expression
 | (4) |
 | (5) |
3 = (1.5 Å)3 and C = 6.022 x 1026. |
RESULTS AND DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTS AND DISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
Electrostatic and transport properties of MscS in the crystal conformation
Four independent BioMOCA simulations of the crystal conformation of MscS were carried out (100 ns each) for symmetric baths of 0.1 M KCl and zero bias potential (0 voltage, see Table 1). Representative snapshots of the time-averaged electrostatic potential as well as the density of negative and positive ions are shown in Fig. 2. The distribution of ions shows a clear pattern with negative ions found mainly in and around the transmembrane pore of MscS and around the distal C-termini (Fig. 2, B and E). On the other hand, positive ions localize in and around the distal zone of the cytoplasmic domain (Fig. 2, A and E). The ionic density is particularly high around some regions of the cytoplasmic domain (Asp195,197,199,262; Arg128,131; and Lys169) that may serve as binding spots for charged species (see Fig. 2 A for positive ions and Fig. 2 B for negative ions). The averaged electrostatic potential follows the corresponding trend: positive at the transmembrane region of MscS and the distal C-termini, negative at the distal zone of the cytoplasmic domain (Fig. 3 A). Spontaneous diffusion of ions through the transmembrane pore (crossings in Table 1) was monitored during these four simulations. The net current was approximately zero in all cases, as expected since the applied voltage was set to zero. However, only Cl– ions crossed the transmembrane pore, whereas positive ions did not cross in either direction, indicating strong selectivity for anions over cations.
|
|
|
,23,26), but the electrostatic pattern and distribution of ions matches MD results reported earlier (21). Average ionic concentration and equilibrium occupancy along the axis of the channel (z) show features consistent with this pattern (Fig. 4, A and C), namely, higher concentration of negative ions at the ends of the transmembrane domain and higher concentrations of positive ions at the distal zone of the cytoplasmic domain. Note that the distribution of ions seems to be rather independent of the biasing potential.
|
All ionic currents monitored in the BioMOCA simulations of MscS in its crystal and alternative conformation crystalNC were smaller than experimental values, particularly at negative potentials (see Fig. 5 A). While experimental ionic currents of MscS at –100 mV are close to –100 pA, the ionic currents observed in BioMOCA simulations are –20 pA. Interestingly, ionic currents observed for positive potentials are in better agreement with the experimental values. These currents are in the range of previous estimates (70–250 pS (20,21)). Note, however, that the crystal conformation is slightly narrower than the crystal structure itself (see Methods). The simulations also indicate an extremely selective transmembrane pore. BioMOCA simulations of MscS's crystal structure using CHARMM radii instead of VMD radii (crystalCR, see Methods and Supplementary Material's Table 3) resulted in even smaller currents and higher selectivity (Table 1), indicating the relevance of the chosen set of radii for the computed MscS's transport properties.
|
p = 20 for the protein region. The corresponding currents (Table 1) are similar to those obtained with p = 5 (larger in magnitude but still smaller than what is expected for the fully conductive conformation of the channel).
Electrostatic and transport properties of MscS in open conformations
Further BioMOCA simulations were performed on three candidates for the MscS open conformation: open1 (rmin = 3.75 Å), open2 (rmin = 5.8 Å), and open3 (rmin = 7.85 Å, see Methods and Fig. 1). It is interesting to note that all the open conformations featured a less pronounced kink at Gly113, leading to more straight TM3 helices than those observed in the crystal conformation, thereby suggesting a possible gating mechanism (see Supplementary Material's Fig. 7). Since further experimental tests are needed to judge these open conformations, we refrain from further speculation.
BioMOCA simulations of the open1 conformation revealed an electrostatic pattern similar to that observed for the crystal conformation (data not shown). Indeed, average ion concentration and occupancy (Fig. 4, B and D) along the axis of the pore are just slightly different from those observed for the crystal conformation. Preliminary studies on MscS homologs show that this pattern seems to be conserved throughout different species as well (L. Kunz, M. Sotomayor, and K. Schulten, unpublished). Currents monitored during simulations without biasing potential were close to zero, as expected, whereas application of a 50 mV biasing potential led to a net ionic current of 37 pA for two simulations lasting 200 and 606 ns (Table 1). When the biasing potential was set to –50 mV the ionic current was found to be –22.4 pA (200 ns); ionic currents in simulations where the biasing potential was set to 100 and –100 mV were found to be 75.3 and –59.3 pA, respectively. The currents are compatible with a conductance slightly lower than 1 nS. Additional simulations performed using symmetric baths of 0.2 M KCl and biasing potentials of 0 and ±50 mV showed similar results with slightly larger currents and no clear rectification (Table 1). Although the net currents observed for open1 are in reasonable agreement with those observed experimentally (slightly smaller, particularly for negative potentials), all currents originated from Cl– ions. Therefore, we further tested selectivity by performing two additional simulations using asymmetric baths of 0.3/0.1 or 0.1/0.3 M KCl. Asymmetric conditions are used in experiments to infer the selectivity of ion channels. A nonselective channel would let both positive and negative ions go down the concentration gradient, thereby generating a zero net current. The ionic currents over 200 ns of simulation under asymmetric conditions were found to be 13.6 and –21.6 pA, mainly driven by Cl– ions, corroborating the strong selectivity observed in the previous simulations (and indicating that at least at low voltages the channel in the open1 conformation may favor conduction of ions from the cytoplasm to the periplasm). In contrast, experimental results (2,26) suggest only slight selectivity. Trajectories of ions were visually inspected for one simulation of the open1 conformation with a biasing potential of 100 mV. Positive and negative ions were found to enter and leave the cytoplasmic domain of MscS through its side openings as previously observed in MD simulations (21), but mainly negative ions passed through MscS's transmembrane pore.
Simulations of the open2 conformation (see Fig. 1 E) performed using symmetric baths of 0.2 M KCl and external potentials of 0 and ±100 mV exhibited ionic currents larger than those observed in experiments for positive potentials (see Table 1). Although the number of positive ions crossing the transmembrane domain of MscS was larger than for crystal and open1 conformations, the pore was found again to be highly selective. On the other hand, ionic currents at negative potentials (likely overestimated as discussed below; see Fig. 5 B) are in reasonable agreement with experiments.
BioMOCA simulations of the open3 conformation, which features the widest transmembrane pore (rmin = 7.85 Å, see Fig. 1, D and F), revealed electrostatic patterns similar to those described above, but ionic currents were found to be larger than those calculated for all the other conformations and those reported in experiments at positive potentials. Although the number of positive charges crossing the transmembrane domain is significantly larger than for conformations crystal, open1, and open2, relative currents continue to show strong selectivity, likely because of the presence of Lys169 at the cytoplasmic end of the transmembrane domain (see Fig. 1 F), and Arg88 located at the periplasmic end of the pore, as well as other charged residues around the cytoplasmic side openings.
Schematic models, selectivity, and rectification
The strong selectivity observed in the simulations of different conformations of MscS led us to design schematic models that could shed light on the origin of this phenomenon. The schematic models consist of an artificial membrane 50 Å thick ( = 2), a cylindrical pore across the center of the membrane with radius 5 Å, and in some cases rings of seven unitary charges placed so as to mimic MscS's charge distribution (see Methods and Fig. 6 A). Snapshots of the time-averaged electrostatic potential and ionic distributions obtained through BioMOCA simulations of schematic models mod1 and mod3 are shown in Fig. 6, B and C, respectively. As expected, ionic distributions are uniform for the all-neutral mod1 model, while negative and positive ions are highly concentrated around the decorating charges placed at the ends of the pore in the mod3 model. Ionic currents for each of the simplified models are listed in Table 2. Permeation events over 200 ns were only observed for models mod2, mod3, and mod4, indicating that conduction of ions was triggered by the presence of charges at the ends of the pore. Furthermore, the channel was switched from being completely selective for anions to completely selective for cations by changing the sign of the charges placed at the ends of the pore (see Table 2). A similar behavior has been observed in modified carbon nanotubes simulated explicitly over few nanoseconds (82). A slight difference between the currents driven by positive and negative charges (mod2 and mod4) was observed and is likely caused by the difference in size and diffusion coefficient for K+ and Cl– ions. The presence of a third, negatively charged, ring on the cytoplasmic side in the mod3 model seems to have little impact on the net ionic current computed. Wider models (rmin = 10 Å) mod5, mod6, and mod7 exhibited a similar behavior (see Table 2). While the all-neutral models mod5 and mod7 showed low conductance and low selectivity, the positively charged rings of mod6 enhanced conductance threefold and induced a high selectivity for anions.
|
|
The schematic models served as control systems as well. Using Hille's equation (72) one can estimate an upper limit to the resistance Rchannel of a channel of radius a and length l bathed in a solution of resistivity 
as
 | (6) |
= 80 
cm, l = 50 Å and a = 10 Å (mod7), one obtains Rchannel = 1.67 G, equivalent to 598 pS. A pore subject to a 100 mV voltage drop across the membrane should then exhibit a current of 60 pA, which is larger than the value observed for mod7 (36 pA). Similarly, a pore of radius a = 5 Å subject to a voltage drop of 100 mV would exhibit a current of 17 pA. The small value may explain why ionic currents were not observed for simulations of mod1 over 200 ns. Likely, the pore would exhibit a current smaller than the limit set by Hille's equation, and an even smaller value might be expected due to obstruction effects caused by atoms decorating the ends and cytoplasmic entrance of mod1. Such a small current is probably within the limits of what can be observed during 200 ns of simulation (one ion crossing equivalent to 1 pA). However, control simulations performed on mod1 using a smaller radius for decorating atoms and a smaller mesh spacing for the grid did not exhibit permeation events either.
Sources of errors
Among multiple sources of error and approximations of the models presented here, four issues need to be addressed in future simulations and modeling of MscS; a brief discussion about each of them follows.
The first problem that needs to be addressed is the effect of missing residues in the crystallographic structure of MscS. The first 26 residues of the N-termini and the last six residues in the C-termini were not resolved in the crystal structure (11) and, therefore, not included in present and previous (21) simulations. While there is some experimental evidence that the missing domains are not fundamental for MscS function (33,34), details of the electrostatic potential at the C- and N- termini may be influenced by missing charged residues (three negatively charged per subunit at the N-termini; two negatively charged and two positively charged residues per subunit at the C-termini) and the charges of the dissociated termini themselves.
The second issue is related to probable deformations of the MscS structure caused by crystallization conditions. The crystal was solved using a detergent solution that may not represent the hydrophobic membrane environment well. Moreover, the large cytoplasmic domain of MscS participates in multiple lattice contacts. Indeed, each molecule in the unit cell is interacting with two or three other molecules (see Supplementary Material's Fig. 8). While the transmembrane domains do not participate in these lattice contacts, TM1–TM2 loops of different molecules are within 4 Å between each other in some cases (see Supplementary Material's Fig. 8 D), which may explain the unusual tilt of MscS's transmembrane helices. In addition, lattice contacts between different domains of the cytoplasmic domain can directly affect the conformation of TM3 and, thereby, affect the architecture of the pore. Indeed, previous MD simulations showed a displacement of the cytoplasmic domain and the TM3 domain toward the membrane (21). Experiments have also shown that co-solvents can induce changes in the cytoplasmic domain and, thereby, affect MscS gating (32).
The third problem is related to the intrinsic limitations of the methodology used in the present investigation. The coarse-grained description utilized allowed us to reach timescales of 100 ns or more, but the spatial accuracy of the model is limited. How to best represent the membrane bilayer as a continuum dielectric material (including explicit lipids/hydration water molecules) remains an open question. Furthermore, the dynamics of the protein is not taken into account and, therefore, conformational changes induced by external fields or hydrodynamics effects are ignored. Moreover, ionic currents presented here represent an upper bound since bulk diffusion constants for ions were used for all regions, even in the narrow transmembrane pore where diffusion of ions may be affected by dewetting transitions and dehydration of ions (not included either in the present model). Smaller diffusion constants would reduce the magnitude of the measured currents without affecting selectivity (70). In addition, changes in the dielectric coefficients and protonation states used for the protein may affect the electrostatic potential and even the strong selectivity observed. Ideally, all-atom simulations could deal with these methodological issues; however, with present computer technology, such in situ simulations involving 225,000 atoms permits only several simulations of tens of nanoseconds (40), compared to multiple simulations of altogether several microseconds presented here.
Finally, the location and composition of the membrane need to be carefully analyzed. Molecular dynamics simulations, theoretical modeling, and experiments (20,21,33) suggest that the position of the membrane is likely to be slightly different from that suggested in Bass et al. (11), with residues 22–26 immersed in the bilayer. Moreover, the overall positive potential of the transmembrane domain may attract anionic lipids present in E. coli. How anionic lipids interact with MscS is unknown, and what the effect of these lipids on selectivity might be, remains to be elucidated. Other bacterial channels are known to depend on anionic lipids to function (83).
|
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTS AND DISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
), whereas the structure revealed a heptamer (11). The structure was suggested to represent a fully open conformation of the channel (11), but further studies make it more likely that the crystal conformation represents instead a nonconductive state (20,23). However, distances measured from cross-linking of site-specific cysteines in a probable closed state do not correspond well with those observed in the structure (30), suggesting large conformational changes upon gating or inactivation (23). Moreover, the size of the pore, when fully open, remains experimentally undetermined and the in-plane area expansion of the channel was first predicted to be 8.4 ± 0.4 nm2 (26), and then estimated to be 18 nm2 (23). While some studies have reported inactivation (7,23,27,31,32,34), others do not mention it (2,26), and the original voltage dependence for activation of MscS (2) has been challenged in more recent studies (23,27). In addition, the cytoplasmic domain was first characterized as essential for the stability of the protein (33), but now is being suggested to be nonessential and only responsible for increased stability and activity (34). Finally, the magnitude of ionic currents at relatively small voltages (>±40 mV) becomes variable due to an abundance of subconducting states (23) and further evidence of a subconductance state has been recently reported (28) using recording techniques with 3-µs resolution. Indeed, the conductances reported for MscS vary from 970 pS (650 pS for depolarizing voltages) (2), and 540 pS (270 pS) (26), to 1250 pS (22), likely due to differences in buffer conductivity.
In the context of the experimental results described above, the simulations presented here, although far from conclusive and mainly of qualitative significance, suggest that the crystal structure does not represent the fully open state of MscS. Slightly wider channels are found in our simulations to reproduce the experimentally determined currents. However, the simulations reveal an unexpected rectification at high voltages and a strong selectivity for anions over cations, which is in relative disagreement with several experimental studies showing only slight anionic selectivity. Several issues may cause this disagreement in selectivity, namely, missing residues and charged N-termini, protonation states of charged residues, or large conformational changes of the structure upon application of tension or voltage. While the crystal structure could represent an inactive state in which the TM3 helices are detached from helices TM1–TM2, which interact with the membrane (20,23), contacts of TM1–TM2 loops with the end of TM3 helices at the cytoplasm (observed in the present and previous studies (21)) may couple TM3 with external helices and membrane. Indeed, forces applied to the end of TM3 helices are required to widen the channel (data not shown), for which TM3 exhibits a less pronounced kink.
The role of the cytoplasmic domain remains unknown, but the simulations presented here indicate that it has a strong influence on the distribution of ions and perhaps on the selectivity of the channel, while its side openings can readily conduct ions. An interesting possibility is that the cytoplasmic domain could act as a molecular sieve (11). Multiple and functional side openings, like the ones observed in this study, seem to give further support to this hypothesis. Upon osmotic shock and exocytosis, these openings would preclude the blockage of the main transmembrane pore by large solutes, and would provide alternative exit pathways for small solutes. Since the main anion in bacterial cells is glutamate (84) and not Cl– (utilized in patch-clamp experiments), it would be interesting to investigate glutamate's interactions and possible pathways through MscS. Another possible role of the cytoplasmic domain may be related to its influence on gating. Further simulations of the entire (including missing residues and charged termini) MscS channel using different co-solvents may clarify how the cytoplasmic domain dynamics and interaction with intracellular elements could gate the channel and, thereby, provide a sensor for cytoplasmic crowding.
|
SUPPLEMENTARY MATERIAL |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTS AND DISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
|
ACKNOWLEDGEMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTS AND DISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
). This work was supported by the National Institutes of Health (NIH grant No. P41 RR05969, NIH grant No. 1 RO1 GM067887, and NIH grant No. R01-GM073655), by the National Science Foundation (Network for Computational Nanotechnology grant No. EEC-0228390), and by the National Center for Supercomputing Applications. The authors also acknowledge computer time provided by the National Science Foundation (National Resource Allocation committee grant No. MCA93S028).
|
FOOTNOTES |
|---|
Submitted on December 20, 2005; accepted for publication February 27, 2006.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTS AND DISCUSSIONCONCLUSIONSSUPPLEMENTARY MATERIALACKNOWLEDGEMENTSREFERENCES |
|---|
2. Martinac, B., M. Buechner, A. H. Delcour, J. Adler, and C. Kung. 1987. Pressure-sensitive ion channel in Escherichia coli. Proc. Natl. Acad. Sci. USA. 84:2297–2301.
3. Sachs, F., and C. Morris. 1998. Mechanosensitive ion channels in non-specialized cells. In Reviews of Physiology and Biochemistry and Pharmacology. M. Blaustein, R. Greger, H. Grunicke, R. Jahn, L. Mendell, A. Miyajima, D. Pette, G. Schultz, and M. Schweiger, editors. Springer-Verlag, Berlin, Germany. 1–78.
4. Hamill, O. P., and B. Martinac. 2001. Molecular basis of mechanotransduction in living cells. Physiol. Rev. 81:685–740.
5. Sukharev, S., and D. P. Corey. 2004. Mechanosensitive channels: multiplicity of families and gating paradigms. Sci. STKE. 219:Re4.
6. Berrier, C., A. Coulombe, I. Szabo, M. Zoratti, and A. Ghazi. 1992. Gadolinium ion inhibits loss of metabolites induced by osmotic shock and large stretch-activated channels in bacteria. Eur. J. Biochem. 206:559–565.[Medline]
7. Levina, N., S. Totemeyer, N. R. Stokes, P. Louis, M. A. Jones, and I. Booth. 1999. Protection of Escherichia coli cells against extreme turgor by activation of MscS and MscL mechanosensitive channels: identific
) and conformation open1 (
) in 100 mM KCl symmetric solution and (B) for conformation open2 (
), (26
) and (23
). Experimental curves have not been normalized to the specific buffer conductivity of different experimental setups and, therefore, provide only an approximate reference.