INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Mechanosensitive (MS) channels play an important physiological role in living cells of diverse phylogenetic origin. They are ubiquitous in prokaryotes and have recently been characterized in archaebacteria (Le Dain et al., 1998) as well as mammals (Patel et al., 1998; Maingret et al., 1999). In eukaryotes, MS channels play a role in such important biological functions as hearing, touch, and cardiovascular regulation (Corey and Hudspeth, 1983). In bacteria, response to the osmolality of their environment is essential for maintaining viability of the cell. In Escherichia coli, three MS channels have been identified, and one of these, MscL, has been cloned (Sukharev et al., 1994). Several studies (Blount et al., 1997; Ou et al., 1998; Ajouz et al., 1998) have confirmed the importance of this channel for osmoregulation of the bacterial cell. A bacterial cell exposed to osmotic downshock experiences an increase in cell membrane tension, which can lead to cell lysis unless the osmotic gradient can be relieved. In these circumstances, MS channels gate to allow K⁺ and other osmoprotectants to be excreted from the cell. In eukaryotes, the stimulation of exocytosis by mechanical strain is thought to be mediated by stretch-activated channels (Xu et al., 1996; Weber et al., 2000).

The MscL protein exhibits a high degree of primary sequence conservation within a group of bacteria that includes E. coli, on which most physiology experiments have been performed, as well as Mycobaterium tuberculosis, from which the crystal structure was obtained. The determination of the crystal structure of MscL (Chang et al., 1998) revealed a protein with a homopentameric structure, approximately 50 Å wide in the plane of the membrane and 85 Å tall. Each 151-residue subunit consists of two transmembrane helices, labeled TM1 and TM2, and a cytoplasmic helix that extends some 35 Å below the membrane. The TMl helices are arranged so as to block diffusion through the channel at their N-terminal ends; this region of the protein also exhibits very high sequence conservation. A loop region between TM1 and TM2 extends into the pore, which may also contribute to the conductance of the channel. Excision of the cytoplasmic domains has been found to have little effect on the gating properties of the channel (Ajouz et al., 2000).

In a series of patch-clamp experiments, Sukharev et al. (1999) characterized the response of MscL from E. coli to tension applied to reconstituted membranes. Conductance measurements indicated that MscL forms a pore at least 30 Å across. Single channel measurements revealed the existence of at least five subconductance states; the rate-limiting step in gating was found to be the transition between the closed state and the first subconductance state, with a barrier of 38 k_BT when no tension is applied. This first transition is also the only tension-sensitive part of the gating process; the free energies of the higher conducting states appear to be insensitive to tension. The authors conclude from their measurements that during the initial application of tension, the in-plane area of E. coli MscL increases greatly without a concomitant increase in conductance; subsequent internal rearrangements give rise to the subconductance states seen in patch-clamp measurements.

Yoshimura et al. (1999) studied the effect of mutating Gly22 to all other natural amino acids. It was found that both the growth rate of mutants as well as the threshold pressure of activation varied directly with the hydrophobicity of the substitution, with particularly acute growth inhibition from acidic substitutions.

Ou et al. (1998) randomly mutagenized mscL in living bacteria and screened for mutants with leaky channels and hampered growth; it was found that most of the mutations in this group could be mapped to between residues 13 and 30, corresponding to one side of the TM1 helix. The authors infer from their data that this face of the helix moves from a hydrophobic to a hydrophilic region during gating, then back to a hydrophobic region. Interestingly, it was also observed that the mutation V23G results in a severe gain-of- function phenotype, though glycine is hydrophobic and smaller than valine. Thus, hydrophobic as well as specific packing interactions must play a role in MscL gating.

Knowledge of the dynamics of MscL, e.g., of its thermal fluctuations, within a lipid bilayer environment would help to evaluate models that have been suggested for the gating of MscL by membrane tension, and may even suggest new mechanisms. Molecular dynamics simulations can give a detailed picture of the dynamics of MscL on the time scale of a few nanoseconds. The effect of surface tension can also be incorporated into molecular dynamics simulations. Examination of specific protein-lipid interactions will shed light on how MscL is able to gate by membrane tension alone. In such simulations, one can rationalize the results of mutagenesis experiments based on the environment observed around the residues. The overall rigidity or flexibility of the protein may give us some insight into the gating mechanism of MscL. Finally, the response of the protein to surface tension may give one clues as to how the protein responds to mechanical stress.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Initial coordinates for the protein were taken from the crystal structure of Chang et al. (1998) (PDB entry lmsl). Residues following Glu104 were excised from the structure; these residues correspond to the C-terminus cytoplasmic helices and have been shown to be nonessential for channel gating and function (Ajouz et al., 2000). Residues 1-9, which were disordered in the crystal structure, were not modeled. The remaining structure (residues 10-104) was modeled using X-PLOR (Brünger, 1988) to place atoms that were absent from the crystal structure and to remove bad contacts.

The membrane used to provide the lipid environment for the protein was constructed from a palmitoyl-oleoyl-phosphatidylcholine (POPC) membrane taken from Heller and coworkers' 1993 molecular dynamics simulation (Heller et al., 1993). The system consisted of an equilibrated rectangular bilayer with 100 POPC lipids in each leaflet. Beginning from this structure, we performed a series of modeling steps in order to prepare a membrane suitable for our protein system. First, a narrow strip of lipids and water along the short side of the membrane was cut off and manually positioned on the long side of the membrane, making the membrane more square. After energy minimization to remove bad contacts, this new structure was equilibrated at 1 atm and 340 K for 1 ns, using the same methodology as described below, i.e., full electrostatics, periodic boundary conditions, and a flexible unit cell. This equilibrated membrane was too small to adequately contain the MscL structure. We constructed a larger membrane patch by replicating the equilibrated structure four times. The completely built protein structure was manually positioned in the membrane; overlapping water and lipid molecules were subsequently removed. After insertion, the protein and a square patch of membrane and water measuring 88 Å on each side was cut out of the large replicated system.

A pre-equilibrated water box was overlaid on the protein-lipid system in order to hydrate completely the aqueous part of the protein structure. Five chloride ions replaced five water molecules to bring the entire system to charge neutrality. The complete system was comprised of 7370 protein atoms, 195 POPC lipids, 7387 water molecules, and 5 chloride ions for a total of 55,666 atoms.

Molecular dynamics simulations were carried out using the program NAMD2 (Kale et al., 1999), with v.26 of the CHARMM force field (MacKerell et al., 1998) for proteins (MacKerell et al., 1992) and lipids (Schlenkrich et al., 1996). Bonds to all hydrogen atoms were kept rigid using SHAKE (Ryckaert et al., 1977), permitting a time step of 2 fs. The system was simulated in periodic boundary conditions, with full electrostatics computed using the particle mesh Ewald (PME) method (Darden et al., 1993) with a grid spacing on the order of 1 Å or less.

The system was energy minimized using the Powell algorithm, then heated for 2 ps under Langevin dynamics at a temperature of 310 K and with a damping coefficient  of 10 ps^l. The system was then equilibrated for 1 ns at constant pressure and temperature. Pressure was maintained at 1 atm using the Langevin piston method (Feller et al., 1995), with a piston period of 200 fs, a damping time constant of 100 fs, and piston temperature of 310 K. Temperature coupling was enforced by velocity reassignment every 2 ps.

Finally, the system was simulated for an additional 2 ns in the NpT ensemble using Langevin dynamics/Langevin piston at a temperature of 310 K and a damping coefficient of 10 ps¹. Coordinates during this phase of the simulation were saved every picosecond.

Simulations were also performed with the same protein structure as above, but with no membrane or water. The Langevin piston method (Feller et al., 1995) was used to control the applied surface tension. Target values for the three Cartesian components of the pressure tensor were set using the formula (Chiu et al., 1995; Feller and Pastor, 1999)
Here P_x, P_y, and P_z are the x, y, and z diagonal components of the pressure tensor and L_z is the size of the unit cell along the z axis, which is perpendicular to the plane of the membrane. The positions and momenta of the atoms in the protein were thereby coupled to changes in the area of the protein through the equation of motion of the Langevin piston. Simulations were performed with periodic boundary conditions, and with a cutoff of 14 Å for non-bonded interactions. Use of PME for the simulation of the bare protein was not warranted, given the qualitative nature of the model; moreover, the rapidly changing size of the dimensions of the unit cell would have made the accuracy of the electrostatics calculation degrade over time.

Analysis of trajectories and energetics was performed using X-PLOR (Brünger, 1992), VMD (Humphrey et al., 1996), and Matlab (MathWorks, Natick, MA).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The system at the beginning of equilibration is shown in Fig. 1. At the beginning of the simulation, the protein was well immersed in the membrane, with no significant gaps between membrane and protein. The membrane patch was sufficiently large to accommodate MscL, whose function depends sensitively on its lipid environment. Throughout this article, we refer to TM1 as residues 15 to 43 of the protein, and TM2 as residues 69 to 89, following the previous nomenclature (Chang et al., 1998).

View larger version (97K): [in this window] [in a new window]   FIGURE 1   (a) Top view and (b) side view of MscL embedded in a POPC membrane at the start of the simulation.

Stability and fluctuations

Stability of the simulated protein can be assessed by analyzing the deviation of the structure from the initial crystal structure. Fig. 2 shows the root mean square deviation (RMSD) of the C atoms of the protein during equilibration and analysis. Preparation and heating of the protein are not shown in the figure; hence the RMS begins at ~2 Å at t = 0. As expected, the transmembrane segments of the protein are considerably more stable than the protein as a whole, which contains loop regions on each side of the membrane that extend into the solvent. It is evident from Fig. 2 that the transmembrane helices are quite stable by the time Langevin dynamics is begun at t = 1 ns. The RMSD during the dynamics period, ~2.5 Å, is comparable to the value obtained in other simulations of transmembrane helices (Forrest et al., 2000; Randa et al., 1999) as well as to that obtained for KcsA (Berneche and Roux, 2000; Shrivastava and Sansom, 2000).

View larger version (14K): [in this window] [in a new window]   FIGURE 2   RMSD of C atoms during the simulation relative to the crystal structure. Solid line: all C atoms; dotted line: only C atoms in the transmembrane helices (residues 15-43 and 69-89). Overall center-of-mass translation of all transmembrane C atoms was removed before calculation of the RMSD.

Fluctuations about the mean positions indicate which parts of the system are most mobile. Fig. 3 shows snapshots of a trace of the backbone atoms during the last 2 ns of simulation. It is evident from the figure that the loop regions of the protein, corresponding to residues 46-68, are the most mobile part of the system, whereas the transmembrane helices are quite stable, both in terms of orientation and in terms of secondary structure.

View larger version (68K): [in this window] [in a new window]   FIGURE 3   Backbone trace at 200-ps intervals during the last 2 ns of simulation.

Fig. 4 shows the RMS fluctuations for the C atom of each helix residue for TMl and TM2 separately. The data show that the least mobile part of the protein corresponds to the first 4 to 5 residues of TM1, which pinch together to form a non-leaky occlusion. The TM2 helices show no such pattern; one of the TM2 helices moves as a rigid body away from its initial position, accounting for its large RMS values. The immobility of this part of the protein is in qualitative agreement with corresponding electron spin resonance experiments (G. E. Perozo and B. Martinac, personal communication).

View larger version (26K): [in this window] [in a new window]   FIGURE 4   RMS fluctuations for C atoms in (a) TMl and (b) TM2.

Analysis of key residues

It is possible that secondary structure formation in the extracellular loop regions of the protein could affect the conductance of the channel. For this reason we examined hydrogen bond formation between backbone atoms in this region of the protein. Three pairs of residues were involved in hydrogen bond formation in at least two of the five subunits: Ile59-Ile67, Ile61-Ile65, and Ile67-Va148. In the case of Ile61-Ile65, the hydrogen bond formed between the backbone oxygen atom of lle61 in one subunit and to the amide hydrogen of the neighboring subunit. These hydrogen bonds exhibited varying degrees of stability, as shown in Fig. 5.

View larger version (40K): [in this window] [in a new window]   FIGURE 5   Backbone hydrogen bonding in the loop region of MscL; the first residue of the pair is the acceptor, the second is the donor. Distances shown are between the oxygen of the acceptor and the hydrogen of the donor. (a) Ile59-Ile67, (b) Ile67-Ile48, (c) Ile61-Ile65; in the last case, the two residues belong to neighboring subunits. Data shown are running averages of 20 data points, sampled at 1-ps intervals.

We focused particular attention on a set of hydrogen bonds formed between Arg45 and Gln51 of neighboring subunits, as shown in Fig. 6. Interaction between these residues has been shown (Maurer et al., 2000) to have a strong effect on MscL gating; when cysteine cross-links are introduced via mutagenesis at these sites, a gain-of-function mutant results. Fig. 6 shows the degree to which these two residues remained in close proximity during the 3-ns simulation. Two pairs of residues interact strongly during the entire simulation period, while a third appears to be drifting toward such an interaction and may have formed a long-lasting hydrogen bond as well. The loop regions in the other two subunits have unfolded too much for the bonds to have a chance to form.

View larger version (30K): [in this window] [in a new window]   FIGURE 6   Hydrogen bonding between Arg45 and Gln51 of neighboring subunits. (a) Distance between glutamine oxygen and sidechain oxygen of Arg45 (see Fig. 7). (b) Distance between glutamine oxygen and amide hydrogen of Arg45. Data shown are running averages of 20 data points, sampled at 1-ps intervals.

We also examined interactions between transmembrane helices, both within subunits and in neighboring subunits. We found only one set of interactions that appeared with any stability in more than one subunit: the interaction between the side chain of Lys33 and the side chains of Ser74 and Asn78 of the neighboring subunit. Fig. 7 shows the distance between these residues during the 3-ns simulation, and a representation of the hydrogen-bonding network formed by these three residues. It is clear from Fig. 7 that at all times, both residues interact strongly with one or more hydrogen atoms of Lys33.

View larger version (29K): [in this window] [in a new window]   FIGURE 7   Inter-subunit hydrogen bonds between Lys33 and (a) Ser74 and (b) Asn78. Data shown are running averages of 20 data points, sampled at 1-ps intervals.

Penetration of water into the pore region of the protein is a key component of this analysis, in that the hydrophobicity of residues lining the pore has been shown to be critical for proper gating of the channel. Fig. 8 shows the extent to which the water molecules penetrated the pore at the end of the simulation. From the extracellular side of the pore, water molecules reach only to Thr25. Residues deeper in the gating region, including Ala20 and Va121, are not exposed to water. The interfacial region near the intracellular side of the protein contains a number of water molecules, but these do not appear to have affected the stability of the pore.

View larger version (84K): [in this window] [in a new window]   FIGURE 8   Water penetration in the pore at the conclusion of the simulation. Solid surface represents protein residues 15 to 40 and 69 to 89. Water molecules within 3 Å of TMl residues are shown. Ala20 is shown in blue, Val21 in tan, and Thr25 in purple. Phosphate atoms in the intracellular leaflet of the membrane are shown for reference as gray spheres.

Simulation of bare protein under surface tension

Though a realistic simulation of MscL must include the membrane and surrounding water, we can investigate the mechanics of the protein itself without these external media. To this end, we conducted a series of simulations of the same protein structure as in the membrane simulations, but with no membrane or water present. The simulations were conducted at constant surface tension and zero normal pressure.

Several factors influenced our choice of surface tension for these simulations. The value must be large enough to provoke a conformational change in the protein before the protein becomes unstable due to its unnatural environment; however, the surface tension must not be so large as to stretch the protein excessively. We ran a series of simulations with a normal pressure of 0 atm and a surface tension in the range of 10-200 dyn/cm. Previous membrane simulations (Chiu et al., 1995; Feller and Pastor, 1999) have suggested that a surface tension of 10-50 dyn/cm gives the best agreement with the measured lipid density of the membrane. It should be emphasized here that our intention was to induce a non-equilibrium conformation change in the protein, while biasing as little as possible the pathway taken by the protein. Under all conditions studied, the protein refolded into an open conformation with minimal loss of secondary structure. The most visible difference between the simulations was the rate at which the protein refolded; this rate was nearly inversely proportional to the applied surface tension.

We describe here one representative simulation carried out with a surface tension of 60 dyn/cm. Analysis was performed for the first 115 ps, after which the rescaling introduced by the constant pressure method caused unphysical large changes in the protein structure.

Fig. 9 shows the radius of the MscL pore during the applied surface tension simulation, computed using the program HOLE (Smart et al., 1993). In the closed state of the channel, and in the snapshots at 50 and 100 ps, there were two primary points of constriction in the channel. Val103 and Glu104 formed the narrowest constriction, and Val21 and Thr25 formed a second constriction at the end of the extracellular pore. A third point of constriction was formed by the extracellular loops comprising residues 44-68.

View larger version (18K): [in this window] [in a new window]   FIGURE 9   Radius of MscL pore as a function of the position along the pore axis, as calculated by the program HOLE, at four points in the simulation of the bare protein. Solid line: t = 0 ps; dotted line: t = 50 ps; dashed line: t = 100 ps; dot-dashed line: t = 113 ps. Protein coordinates at times later than t = 0 were fitted to the coordinates at t = 0 using the positions of the C atoms in TM1 and TM2.

During the first 100 ps of simulation with applied surface tension, the extracellular loops retracted from the center of the pore, resulting in the expansion of the pore radius by about 4 Å at z = 10-30 Å. The total in-plane area of the protein, as well as the angle formed by transmembrane helices with the membrane normal, remained essentially unchanged. During the next 13 ps, a dramatic shift in the tertiary structure of the protein took place. Both TM1 and TM2 helices tilted downward, producing a shortening of the total length of the pore. Val21 and Thr25 moved apart to allow water molecules to diffuse through this part of the channel. The total in-plane area of the protein increased from ~2300 Ų to 4100 Å². Val103 and Glu104 were pulled up toward what would be the interior of the membrane as the TM2 helices tilted downward; this is seen in the shift of the minimum pore radius in Fig. 9. Though these residues were sufficiently far apart to allow diffusion of water and ions, they still diminished significantly the total conductance of the channel. Fig. 10 illustrates the conformational changes associated with the gating during this simulation.

View larger version (55K): [in this window] [in a new window]   FIGURE 10   MscL pore structure at (a) t = 0 and (b) t = 113 ps. The solid blue region in each picture corresponds to pore radius as calculated by the program HOLE. Pink tubes correspond to TM1 residues, white tubes represent the extracellular loop region (residues 44-68), and purple tubes correspond to residues 89-104. In both pictures, residues Val21, Thr25, Val103, and Glu104 are shown in space-filling representation; (a) also shows Ile14 near the bottom of the figure and Gly63 near the top. TM2 helices are omitted for clarity.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We have described the first molecular dynamics simulations of a representative of an important and ubiquitous class of proteins, the mechanosensitive ion channels. These simulations were only recently made possible by the determination of the crystal structure of MscL. Our simulations proceeded along two avenues of investigation: first, we thought to explain the nature of the protein in the closed state when placed in a realistic membrane and water environment; second, we asked how the protein would react to application of surface tension. Though it might still prove possible to investigate both questions with a single simulation, for reasons of computational efficiency we sought to divide our efforts between two different approaches.

Results of the simulations of the protein in the full membrane-water system reveal a protein that is quite stable in the closed state. This is to be expected from patch-clamp data (Sukharev et al., 1999), which reveal a channel that is absolutely nonleaky until significant tension is applied. Large-scale changes in the shape of the protein could not be expected during the progress of a 3-ns simulation; it is possible, therefore, that a much longer simulation could reveal a somewhat different closed state. We believe that we have described the essential features of this protein on the time scale of several nanoseconds, and find encouraging correspondence with experiments. Fluctuations on the scale of individual residues were found to be in good agreement with corresponding measurements from electron spin resonance experiments, confirming the validity of our protein model. Water penetration in the pore was found to extend only to hydrophilic residues, i.e., only as far as Thr25. This result lends support to proposed mechanisms of MscL gating that postulate a change in the solvent environment of hydrophobic residues in the constricted region of the protein during gating.

Our simulations of the bare protein using an applied surface tension to induce conformational change provided remarkably consistent results: the protein retained its secondary structure while radically reforming its tertiary structure to form a large pore. Retention of secondary structure was an important validity check, since the native lipid environment would not have allowed alternative hydrogen bonds to form. The observation that the transmembrane helices flattened out corresponds well with recent measurements made of the effect of membrane thickness on MscL gating (Kloda and Martinac, 2001). In these measurements, it was found that when MscL was placed in a thinner membrane, it remained in its open state for a longer period. This would seem to suggest that the open conformation of MscL is flatter than the closed structure. The simulation of the bare protein also suggested a role for particular sets of residues in the gating process. Val21 and Thr25 formed a tight constriction at the end of the extracellular pore, which was abolished only after significant tilting of the transmembrane helices. Val103 and Glu104 retained a constricted arrangement even after the rest of the pore expanded. It remains to be seen whether these observations remain valid in a fully hydrated environment; the absence of water in the pore could have contributed to a decrease in stability of the constriction formed by Val21. The interactions between MscL and the surrounding bilayer are also quite complex, and we are currently pursuing simulations in which tension is applied to a complete protein-membrane-water system.

Molecular dynamics simulations of membrane protein systems are becoming increasingly common even as the systems' sizes continue to grow (Lin and Baumgaertner, 2000; Berneche and Roux, 2000; Shrivastava and Sansom, 2000). Our simulations were of comparable duration and of slightly larger size than other related studies; this reflects in part our desire to construct a membrane large enough to encapsulate the protein fully within a realistic environment without the effects of boundary conditions. There is accumulating evidence (Chiu et al., 1995; Feller and Pastor, 1999) that surface tension plays an important role in obtaining the most realistic membrane structure. Though a number of investigators have simulated helices and proteins embedded in membranes without the use of a nonzero surface tension and obtained apparently reasonable results, we expect to see more simulations in which surface tension is controlled. Due to the large size of thermodynamic fluctuations in two-dimensional systems such as membranes, and due to the slow time scales in which lipid diffusion takes place, a substantial computational effort must be made in order to pursue these investigations.

In the case of MscL, we have even more reason to look closely at the role of surface tension. Our simulations of MscL under the influence of surface tension are to our knowledge the first application of this type of external force to study conformational changes in proteins. Our work on MscL will continue with constant surface tension simulations of the protein embedded in the membrane-water system.