INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Mechanosensitive channels, which gate in response to the application of tension to a phospholipid membrane, are implicated in a wide range of biological processes, such as touch, hearing, and circulation (Hamill and McBride, 1995). Despite this importance in multiple contexts, the mechanism(s) by which channel proteins can sense membrane tension and the structural changes associated with gating are not well understood. A more complete understanding of mechanosensitive signal transduction processes is increasingly desirable with the discovery of vertebrate mechanosensitive channels, such as the recently cloned vanilloid receptor-related osmotically activated channel (VR-OAC) (Liedtke et al., 2000). Recently, the 3.5-Å-resolution crystal structure of the Mycobacterium tuberculosis homolog of the bacterial mechanosensitive channel of large conductance (Tb-MscL) was reported by Rees and co-workers (Chang et al., 1998). This structure provides a unique opportunity to consider the structure-function relationships of a mechanosensitive channel on an atomic level. MscL, which is thought to protect bacterial cells from severe osmotic downshock (Levina et al., 1999; Nakamaru et al., 1999; Wood, 1999), is the best characterized mechanosensitive ion channel. The Escherichia coli homolog (Eco-MscL) was the first to be cloned (Sukharev et al., 1994) and has received the most attention (Batiza et al., 1999; Blount and Moe, 1999; Oakley et al., 1999; Spencer et al., 1999; Sukharev, 1999; Sukharev et al., 1997). One focus of these studies has been the determination of single site mutations that affect gating properties. In particular, several "gain-of-function" (GOF) mutations have been characterized in which a lower tension is needed to gate the channel (Blount et al., 1997, 1996b; Ou et al., 1998; Yoshimura et al., 1999). Recently, Tb-MscL has been subjected to more limited structure-function studies, a major goal of which has been to determine the extent to which the extensive studies on Eco-MscL can be translated to Tb-MscL (Maurer et al., 2000; Moe et al., 2000).

Although the Tb-MscL crystal structure has given increased insight into many aspects of MscL function, there are still some ambiguities in the structure. For example, the crystal structure revealed Tb-MscL to be a homopentamer with an -helical bundle formed by the intracellular C-terminal regions of each subunit (Fig. 1 A). This region has several acidic residues (Asp and Glu) that would be expected to repel one another electrostatically in a bundled structure at physiological pH. However, the crystallization protocol that succeeded in producing acceptable crystals involved a relatively low pH (3.6-3.8) (Chang et al., 1998; Spencer et al., 1999), which would protonate Asp and Glu, potentially reducing their repulsion and stabilizing the -helical bundle structure. It is thus interesting to consider the extent to which pH may have influenced this interesting structural element.

View larger version (52K): [in this window] [in a new window]   FIGURE 1   (A) Tb-MscL crystal structure with the five subunits shown in different colors. Also, the intersubunit H-bond between Arg 45 and Gln 51 residues in adjacent subunits and the point (Asp 108) where the structure is cleaved for embedded simulations are highlighted. (B) The embedded wild-type Tb-MscL simulation system after 100 ps of MD simulation, before restraints were gradually removed from protein C atoms. The system included 495 protein residues, 290 POPE molecules, and 17851 SPC water molecules for a total of 73,313 atoms. The protein is shown in white CPK with charged residues (Lys, Arg, Asp, Glu) highlighted in red. The lipid and water are shown in green and blue wireframe, respectively, with lipid phosphorous atoms in yellow CPK. (C) A single subunit of the Tb-MscL crystal structure with the regions used for embedded MD analyses shown in different colors: blue, TM regions; red, extracellular loop; green, C terminus. In addition, the two residues altered by computational mutations in embedded MD simulations, Val 21 and Gln 51, are shown as white space-filling models.

Here we report molecular dynamics (MD) simulations utilizing the Tb-MscL crystal structure to help develop structure-based explanations of these and other experimental observations. Since the early simulations of lipid bilayer systems (Egberts and Berendsen, 1988), MD calculations of proteins with explicit lipid and solvent molecules have become increasingly common. Many of these simulations have been reviewed recently (Forrest and Sansom, 2000). Some work has focused on MD simulations of ion channels, considering small channels such as gramicidin A (Chiu et al., 1999a,b; Tang et al., 1999; Woolf and Roux, 1996) or ion channel models often consisting of transmembrane (TM) helix bundles (Capener et al., 2000; Fischer et al., 2000; Forrest et al., 2000; Husslein et al., 1998; Law et al., 2000; Lin and Baumgaertner, 2000; Randa et al., 1999; Schweighofer and Pohorille, 2000; Tieleman et al., 1999a, 1998a). The determination of crystal structures for KcsA and OmpF also allowed simulations of more complete channels in both explicit lipid (Bernèche and Roux, 2000; Shrivastava and Sansom, 2000; Tieleman and Berendsen, 1998) and lipid-mimetic environments such as octane (Guidoni et al., 1999, 2000). These MD studies described the dynamics and conformations of protein structures, the movements and orientations of water and lipid molecules, and interactions between protein, lipid, and water throughout the simulations. Unlike KcsA and OmpF, the Tb-MscL channel has a domain that extends significantly outside the membrane. Roughly 15-20 C-terminal residues of each subunit in the crystal structure are intracellular, forming the helical bundle discussed above. This complicates the simulation, in that a much more extensive water layer must be present on the intracellular side.

This paper discusses two sets of MD calculations. The first set focused on determining the effect of pH on the structure of the C-terminal region of the channel. These simulations suggest that the observed C-terminal helical bundle conformation may have resulted from the low pH crystallization conditions. In addition, specific mutations that appear to promote the helix bundle conformation have been identified. The second set of simulations included one wild-type and two mutant (V21A and Q51E) trajectories of Tb-MscL embedded in an explicit palmitoyloleoylphosphatidylethanolamine (POPE) membrane, totaling over 4 ns in length. Previous work has indicated that both mechanosensitive responses in general (Berrier et al., 1989) and those due to MscL in particular (Blount et al., 1996a; Häse et al., 1997) are associated with the E. coli inner membrane. In the M. tuberculosis inner membrane, phosphatidylethanolamine (PE) lipids are more prevalent than other lipids for which parameters are readily available, such as phosphatidylcholine (PC) (Lee et al., 1996). POPE lipid was also used in a previous simulation of the E. coli OmpF porin (Tieleman and Berendsen, 1998). Although some recent simulations have begun to consider heterogeneous membranes, such as lipid/steroid (Pasenkiewicz-Gierula et al., 2000; Smondyrev and Berkowitz, 1999, 2000; Tu et al., 1998) or two lipid mixtures (Murzyn and Pasenkiewicz-Gierula, 1999), simulations of embedded proteins have thus far included only single lipids.

The present simulations provide insights into the general dynamics of the wild-type Tb-MscL channel and differences in the V21A and Q51E mutant channel dynamics potentially related to their observed phenotypes. Other analyses of the trajectories point to protein-lipid interactions that may increase our understanding of tension transduction between lipid and protein.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

C-terminal region simulation setup

Simulations of the C-terminal region used residues Tyr 94-Arg 118 from all five subunits of the Chang et al. (1998) crystal structure (1MSL) as a starting conformation. Uncharged N- and C-termini (-NH₂ and -COOH) were used, as neither end is a real protein terminus, and the additional charges could cause spurious interactions. Two simulations were performed, which we will refer to as pH 7 and low pH. For the pH 7 simulation, all ionizable residues (Asp, Glu, Lys, and Arg; no His are involved) in the protein were charged. Neutral forms of Asp and Glu were used in the low pH simulation. Additional simulations were performed with various Glu/Gln and Asp/Asn mutations. These were made to the original crystal structure by changing the second acidic O atom of the residue of interest to an amide N in all five subunits. All remaining ionizable residues were charged in these mutant simulations. These protein structures were then solvated by a box of simple point change (SPC) waters extending 1 nm from the protein in each direction. Ions were added randomly to these boxes: 10 Na⁺ and 10 Cl ions for pH 7, 20 Cl ions for low pH, 7 Na⁺ and 12 Cl for single mutations, and 5 Na⁺ and 15 Cl for the double mutant. These ions were added both to neutralize the charge in the low pH box and to maintain an approximately equal ionic strength between the simulations. The simulations included a total of ~27,000 atoms, and after equilibration the system boxes were ~6.8 nm × 6.8 nm × 6 nm.

These structures were minimized for 50 steps of steepest-descents minimization to reduce close contacts. Minimization was followed by a MD simulation that heated the system to 310 K over 20 ps; after heating, all simulations were performed at 310 K. This and all subsequent MD runs on these systems included restraints on the C atoms of the two N-terminal residues of each subunit to represent conformational restrictions imposed by the overall assembly of the channel.

Embedded Tb-MscL simulation setup

The Tb-MscL crystal structure (1MSL) solved by Chang et al. (1998) was used as a starting structure for all MscL simulations. Previous experimental work on Eco-MscL has shown that C-terminal cleavage at Ala 110 or later has no significant effect on channel function (Blount et al., 1996b). Thus, all simulations used Tb-MscL with a similar cleavage after Asp 108, based on sequence alignment to Eco-MscL, to reduce system size. The nine residues missing from the crystal structure N-terminus were not included as their conformation is unknown. Uncharged N- and C-termini were used. All potentially ionizable residues (no His are present in Tb-MscL) were charged, as they generally were exposed to the pore water or to bulk solvent after embedding in the membrane.

An equilibrated hydrated membrane of 340 POPE lipid molecules and 6728 SPC waters similar to that used for simulations of OmpF was obtained from D. P. Tieleman (Tieleman and Berendsen, 1998). Throughout, the membrane will be oriented in the xy plane with the z axis as the membrane normal. This membrane was equilibrated further for 500 ps of constant pressure and temperature (NPT) simulation at 310 K before embedding Tb-MscL. Embedding of Tb-MscL was performed by superimposing the channel structure onto the lipid and removing a minimal number of lipids that were in the pore region or had extensive overlap with MscL. The optimal positioning of the channel in the Z-direction was somewhat ambiguous; patterns of charged, polar, and aromatic residues were used to estimate the proper position. Control simulations with only the TM helices embedded in different positions implied that the positioning used was in an appropriate range (K. Philipson, D. Elmore, and D. Dougherty, unpublished data). This system was then solvated fully, including inside the pore, with SPC waters (Berendsen et al., 1981), which have been shown to reproduce experimental properties well in membrane simulations (Tieleman and Berendsen, 1996). Sufficiently large layers of water were placed above and below the membrane to prevent the intracellular and extracellular Tb-MscL domains from being within the electrostatic cutoff of one another through the box edges. This total system, which was initially ~9.6 nm × 9.5 nm × 10.5 nm, included 495 protein residues, 290 POPE lipids, and 17,851 SPC waters for a total of 73,313 atoms and is shown in Fig. 1 B. This initial system was used to begin the wild-type and Q51E simulations.

The wild-type simulation began with 150 steps of steepest-descents minimization on the embedded structure to reduce close contacts. The system was then heated to 310 K over 20 ps with restraints on all C atoms of the protein; after heating, all simulations were run at 310 K. Restraints were used for an additional 80 ps of simulation and were then released in gradual steps over the subsequent 165 ps. This was followed by 1335 ps of simulation without restraints, for a total trajectory of 1500 ps.

The Q51E mutation was made to the embedded structure by deleting the two amide protons of the Glu 51 side chains and changing their amide N atoms to unprotonated acidic O atoms. These changes were made for all five subunits of the channel. The Q51E mutation also required the addition of five sodium cations to neutralize the overall charge of the system box. The simulation followed the same protocol as the wild type, except that 50 steps of minimization were before and 100 steps were after the addition of Na⁺.

The V21A mutant trajectory began with the 700-ps structure from the wild-type trajectory. The V21A mutation was made by deleting the two Val 21 C methyl groups and converting the C to a methyl group for all five subunits. This mutated structure was then run at 310 K without restraints for 1100 ps. A schematic of the embedded simulations performed can be seen in Fig. 2.

View larger version (10K): [in this window] [in a new window]   FIGURE 2   Schematic of the embedded Tb-MscL MD trajectories performed in this study.

MD simulations

All minimizations and MD simulations were performed using GROMACS 2.0 (Berendsen et al., 1995). Standard GROMACS parameters were used for the C-terminal simulations and the embedded protein. The MD protocol used for membrane embedded simulations in this study followed those applied successfully for previous embedded full channel structure systems (Bernèche and Roux, 2000; Shrivastava and Sansom, 2000; Tieleman and Berendsen, 1998). Lipid parameters were taken from previous work by Berger et al. (1997), with extra parameters for the oleoyl double bond taken from the GROMOS force field. All MD runs used a time step of 2 fs along with the LINCS routine to constrain bond lengths (Hess et al., 1997) and SETTLE to constrain water geometries (Miyamoto and Kollman, 1992). The NPT ensemble with periodic boundary conditions was employed for all runs with pressure coupling constants (_p) of 1.0 ps (Berendsen et al., 1984). C-terminal runs used isotropic pressure coupling to 1 bar. Embedded simulations used anisotropic pressure coupling to 1 bar in each of the three directions individually. This ensemble should allow the embedded system to adjust its surface area appropriately after the embedding procedure. System temperatures were coupled separately for protein, lipid, solvent, and ions to a temperature bath with a coupling constant (_t) of 0.1 ps (Berendsen et al., 1984). For the C-terminal region simulations, Lennard-Jones and short-range Coulombic interactions were cut off at 1.0 nm, and long-range electrostatics were calculated with particle-mesh Ewald (PME) methods, using Fourier grid spacing of 0.10 nm and cubic interpolation. In the embedded systems, a twin-range cutoff was used for determining non-bonded interactions: 1.0 nm for Lennard-Jones and short-range Coulombic and 1.8 nm for long-range Coulombic interactions. Although there may be some advantage to using Ewald sum methods to calculate long-range electrostatics in lipid membrane systems (Feller et al., 1996), the cutoff method employed here is less computationally expensive and has been used previously in embedded channel systems using these lipid parameters (Capener et al., 2000; Fischer et al., 2000; Forrest et al., 2000; Randa et al., 1999; Shrivastava and Sansom, 2000; Tieleman and Berendsen, 1998; Tieleman et al., 1999a, 1998a).

Analyses of the trajectories were primarily performed with tools included in the GROMACS suite (Berendsen et al., 1995). Properties were averaged over the last 1000 ps of C-terminal and embedded trajectories unless otherwise noted. Additionally, the DSSP program was used to determine secondary structure (Kabsch and Sander, 1983), and the HOLE program was used to calculate pore diameter profiles (Smart et al., 1997). Channel regions used for analyses were TM1 (15-43), extracellular loop (44-69), TM2 (69-89), and C-terminal region (90-108), analogous to those defined in previous sequence analyses and highlighted in Fig. 1 C (Maurer et al., 2000). Hydrogen bond and salt bridge interactions are defined geometrically as interactions in which the distance between the hydrogen and acceptor atoms is <0.25 nm and the interaction angle is 60°. Amide N atoms were omitted as possible H-bond acceptors. Some figures were made using MOLSCRIPT (Kraulis, 1991) or RasMol (Sayle and Milner-White, 1995).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Effect of pH on the C-terminal helix bundle

The pH 7 and low pH trajectories were noticeably different from one another. The pH 7 simulation had a much larger C root mean squared (RMS) deviation from the crystal structure than the one at lower pH, showing the overall conformation was more preserved in the lower pH simulation (Fig. 3 A). In addition, the low pH case exhibited consistently lower C RMS fluctuations for the majority of the structure (Fig. 3 B). Although interpretation of RMS fluctuations for simulations on this timescale must be viewed with some caution, we anticipate that they will provide some qualitative insights into the relative stabilities of the trajectories.

View larger version (22K): [in this window] [in a new window]   FIGURE 3   (A) RMS deviation of C from the crystal structure for the pH 7, low pH, and E104Q/D108N C-terminal region simulations, calculated for each frame after fitting the trajectory C to the crystal structure.(B) RMS fluctuation of C around an average structure for C-terminal region simulations plotted per residue. The computed RMS fluctuation is averaged over the five chains. Single mutant fluctuations are omitted for clarity. Values for the E102Q and E116Q simulations were generally in the range of the pH 7 simulation, whereas E104Q and D108N showed lower fluctuations than pH 7 from the N-terminus to the mutation sites. (C) RMS deviation of C from the crystal structure for amide mutation trajectories, calculated as described in A.

In addition, the -helical secondary structure is much better preserved in the low pH simulation, although many of the -helices have been largely eroded by the end of the pH 7 trajectory. This can be seen by DSSP analysis per residue over the trajectories (Fig. 4) and from the final structures in Fig. 5. In the pH 7 simulation, the acidic groups appear to have repulsed each other, leading to the overall deformation of the crystal structure conformation. Such a repulsion is not seen with the neutralized residues in the low pH case, as can be seen by multiple Glu 104 and Asp 108 side chains pointing toward one another in the final structure.

View larger version (29K): [in this window] [in a new window]   FIGURE 4   Analysis of secondary structure for the C-terminal region simulations performed with the DSSP method for the pH 7, low pH, and E104Q/D108N simulations. Black portions represent residues that are -helical at a given point in the trajectories. Results for one representative subunit from each trajectory is shown.

View larger version (65K): [in this window] [in a new window]   FIGURE 5   Ribbon diagrams of the final frames (1500 ps) of the pH 7, low pH, and E104Q/D108N C-terminal region simulations. The analogous region of the crystal structure is also shown. In all structures, the acidic residues (or their mutated amide counterparts) are shown in CPK.

The present simulations suggest that the intracellular helix bundle of Tb-MscL was stabilized by the low pH of the crystallization medium, a possibility suggested by Rees and co-workers (Chang et al., 1998; Spencer et al., 1999). However, such a helix bundle might still be physiologically relevant, perhaps stabilized by binding to metal ions or polyamines. In a case where many acidic residues are constrained into close proximity they could exhibit a shifted pK_a, causing some or all to remain protonated under physiological conditions. However, such a pK_a shift seems unlikely in this situation as the residues C-terminal to the helix-bundle region were too disordered to be built into the crystal structure, and those N-terminal to the region are flexible enough to allow their rearrangement. In this apparent absence of structural constraints on the region, it would likely be energetically more favorable to have fully solvated ionized residues than to form a helical bundle with those residues buried and neutralized.

Additionally, ambiguity in the physiological relevance of the C-terminal helical bundle does not imply any uncertainty in the overall pentameric channel conformation observed in the crystal structure. In fact, the crystal structure conformation, particularly that of the TM helices, was well preserved in simulations of the channel described below with ionized side-chain protonation states. These channel simulations used a C-terminally cleaved structure because of evidence that much of the C-terminal region is not necessary for function in Eco-MscL as well as the uncertainty in the physiological C-terminal conformation.

Mutations to the C-terminal helix bundle

To determine whether diminishing electrostatic repulsions in the C-terminal region could stabilize the helix bundle at pH 7, we have selectively converted acidic residues (Glu and Asp) to their amide counterparts (Gln and Asn). These Asn and Gln side chains would be extremely similar sterically and electrostatically to the protonated Asp and Glu. Thus, trajectories for the C-terminal region with single mutations of all four acidic residues (E102Q, E104Q, D108N, and E116Q) were computed. For the most part these single mutations led to trajectories similar to that of the wild-type pH 7 simulation, as measured by C RMS deviation from the crystal structure (Fig. 3 C), C RMS fluctuation, and secondary structure conservation by DSSP (data not shown). However, visual inspection of the E104Q and D108N trajectories indicated that these mutations did more nearly conserve the helical bundle, albeit far less than seen in the low pH situation. As well, there was some reduction in C RMS fluctuation from the N-terminus to the mutation site in each of those simulations (data not shown).

These hints of partial stabilization with the E104Q and D108N mutations led us to perform an additional trajectory with an E104Q/D108N double mutant. As can be seen by C RMS deviation (Fig. 3 C), C RMS fluctuation (Fig. 3 B), and DSSP secondary structure determination (Fig. 4), the combined effect of the two mutations led to a trajectory equivalent to that of the low pH simulation. The final structure from the E104Q/D108N simulation is also strikingly similar to that from the low pH simulation (Fig. 5).

Thus, an E104Q/D108N double mutation may have similar effects to low pH conditions on the structure of the Tb-MscL C-terminal helix bundle. It would be interesting to experimentally investigate the effects of this mutation on both helix bundle propensity and overall structural stability. As well, determining how, if at all, such mutations alter Tb-MscL function may give insight into a possible functional role of the Tb-MscL C-terminal region.

Wild-type MscL simulation: protein structure

The average MD structure from the final 1 ns of the embedded wild-type Tb-MscL trajectory was generally close to that of the crystal structure, as can be seen by visual inspection (Fig. 6). This relative conservation of the crystal structure conformation with all side chains ionized suggests that neutral (low pH) residue ionization states would not alter the trajectory as drastically as in the C-terminal bundle simulations. In particular, the TM domain structure remained very similar to the crystal structure. However, a greater structural difference is seen in the extracellular loop and C-terminal regions. These trends of structural conservation are also apparent in the computed C RMS deviation after fitting to the crystal structure over the trajectory. The RMS was calculated both for all C and for C in specific regions of the protein (Fig. 7). The overall C RMS of ~0.33 nm is somewhat higher than that of previous channel simulations, which generally have reported values of 0.18-0.25 nm (Bernèche and Roux, 2000; Shrivastava and Sansom, 2000; Tieleman and Berendsen, 1998). However, the RMS deviation of the MscL TM helices (~0.17 nm) is quite comparable to previous simulations of both KcsA and OmpF, both of which have much less extensive extra-membrane regions. As well, the RMS deviation for Tb-MscL TM domains was similar to those seen for stable TM bundle calculations (Capener et al., 2000; Fischer et al., 2000; Forrest et al., 2000; Husslein et al., 1998; Law et al., 2000; Lin and Baumgaertner, 2000; Randa et al., 1999; Schweighofer and Pohorille, 2000; Tieleman et al., 1999a, 1998a). It is not particularly surprising that an increased deviation from the crystal structure was seen in the extracellular loop and C-terminal regions, as those regions had greater uncertainty in the crystal structure. In addition, repulsions between acidic residues in their negatively charged physiological protonation states caused large structural deformations in the C-terminal simulations discussed above. The computationally imposed cleavage at Asp 108 in the C-terminal region may also impact the RMS deviation in that region.

View larger version (65K): [in this window] [in a new window]   FIGURE 6   Superimposed C traces for the Tb-MscL crystal structure (light gray) and the average structure for the final 1000 ps of the wild-type simulation (dark gray). The crystal structure is cleaved at Asp 108 as was done in the simulated channel.

View larger version (33K): [in this window] [in a new window]   FIGURE 7   RMS deviation of C from the crystal structure for the wild-type Tb-MscL trajectory. RMS values were determined both for the full protein and for selected regions and were calculated for each frame after fitting the C for the region under consideration to the crystal structure.

The Tb-MscL extra-membrane regions are also more conformationally flexible than the TM domains, as can be seen by the substantially larger C RMS fluctuations in these regions (Fig. 8). The increased fluctuation is in agreement with the increased temperature factors for those regions in the crystal structure (Chang et al., 1998). Increases in fluctuation, although smaller in magnitude, were also seen for loop regions in the OmpF simulation (Tieleman and Berendsen, 1998).

View larger version (32K): [in this window] [in a new window]   FIGURE 8   RMS fluctuation of C around an average structure for the wild-type and Q51E trajectories plotted per residue. Fluctuations for the V21A mutant trajectory were very similar and are omitted for clarity. The computed RMS fluctuation is averaged over the five chains.

The relative structural stability of the 10 TM helices is also apparent in their conservation of -helical structure throughout the simulation, as measured with DSSP (data not shown). As well, no other regions of consistent secondary structure, such as -helices or -sheets, arise during the simulation.

Structural changes in mutant trajectories

The wild-type, V21A, and Q51E trajectories display similar C RMS deviations from the crystal structure (data not shown), and all three show similar C RMS fluctuations (Fig. 8 A). The C RMS deviations from the crystal structure for the TM, extracellular loop and C-terminal regions of the mutants likewise were very similar to those for the wild type. As well, secondary structures throughout the mutant trajectories were similar to those in wild type, with the 10 TM helices well preserved (data not shown). Thus, neither of the mutant trajectories showed a significantly larger difference from the crystal structure than wild type. As well, the equilibrium structures from all three trajectories were at least as similar to each as other as they were to the crystal structure (data not shown).

Although no significant global change in structure is seen for either mutant, the V21A mutant does display an increased pore radius in the general region around residue 21 (Fig. 9). Some increase was expected given the reduction in steric size of this putative plug residue from Val to Ala. This increase from steric reduction can be seen by comparing the radius of the 700-ps wild-type frame (the frame from which the V21A trajectory was begun) with and without the mutation. However, as the V21A trajectory equilibrates, the observed pore widening extends through a longer section of the channel than seen for the initial change from Val to Ala. It has been postulated that a V21A mutation may decrease the van der Waals interactions among pore-lining TM1 helices, weakening a hydrophobic lock that stabilizes the channel closed state (Moe et al., 2000; Yoshimura et al., 1999). The simulations support this view in that there is a widening of the pore, perhaps indicative of weakened helix-helix interactions, across a significant section of the channel. Note, however, that there is no evidence for an increased structural instability of the TM1 region in either comparisons of C RMS fluctuations for wild-type and V21A simulations (Fig. 8) or conservation of -helical secondary structure measured by DSSP (data not shown).

View larger version (15K): [in this window] [in a new window]   FIGURE 9   Pore radius profiles for the wild-type and V21A Tb-MscL simulations in the membrane region calculated with the HOLE program. Each profile given is the average of profiles taken every 100 ps over the final 500 ps of the simulations. The calculated pore radius of the initial V21A mutation frame (the 700-ps wild-type frame with the V21A mutation) is also shown to show the change in pore radius resulting from the Val to Ala steric change before equilibration of the mutant structure. Numbering on the z axis begins with the intracellular side of the simulation box at 0 nm. The average position of the plug residue 21 is denoted.

R45-Q51 hydrogen-bonding interaction

An inter-subunit hydrogen bond is apparent between the R45 and Q51 side chains in adjacent extracellular loop regions of the Tb-MscL crystal structure. The actual proximity of these residues has been verified through designed cross-linking reactions, and several mutations to these residues, such as Q51E, cause a distinct GOF phenotype, hinting at a possible physiological role of this interaction (Maurer et al., 2000). A goal of the Q51E simulations was to see whether a structural correlate of this functional effect could be found. To that end, the presence of hydrogen bonding (or in the case of Q51E, salt bridge formation) between these residues in the trajectories was considered.

The populations of the five possible inter-subunit interactions between residues 45 and 51 show very different patterns in the wild-type and Q51E Tb-MscL trajectories (Fig. 10). In the wild-type trajectory after ~250 ps, a maximum of two of the possible five hydrogen bonds are usually formed at any one time. Two of the residue pairs (pairs 3 and 4) drift apart after that first portion, and another (pair 1) is very rarely formed as the simulation continues. As well, one of the pairs that shows a reasonable interaction is fairly flickery (pair 5). Notably, one of the pairs that drifted apart early in the simulation (pair 3) does form a few isolated interactions later in the trajectory, showing that interactions are not prevented because of some major deformation in the loop structure. Similar trends of hydrogen bond occupancy are seen for the V21A mutant (data not shown).

View larger version (24K): [in this window] [in a new window]   FIGURE 10   Hydrogen bond/salt bridge interaction populations of the five possible R45-Q51 pairs in the wild-type Tb-MscL simulation (black) and of the five possible R45-E51 pairs in the Q51E simulation (gray). The average interaction lifetimes were calculated over the final 1000 ps of each simulation.

In contrast, all five possible salt bridges are formed early in the Q51E trajectory (Fig. 10). After that point, two of the interactions essentially break (pairs 2 and 3); two are almost always formed (pairs 4 and 5); and one other fluctuates (pair 1). The net result is that the 45···51 inter-subunit interaction is significantly more prevalent in the Q51E mutant than in the wild type. This is confirmed by the fact that the dwell time of salt bridges in the Q51E trajectory is considerably longer than that of the hydrogen bonds in the wild-type trajectory (Fig. 10). In other words, when a salt bridge is formed in Q51E, it remains intact longer on the average than a hydrogen bond that is formed in the wild type.

These results show that there is indeed a structural change in the protein that can be associated with the Q51E mutation. There is a definite change in the inter-subunit 45···51 interactions, which may be related to the observed GOF phenotype in Q51E mutant channels. This opens up the exciting possibility of using MD simulations to rationalize and ultimately predict the functional consequences of structural changes in an ion channel. Of course, simulations of other known GOF mutants at residues 45, 51, and elsewhere, along with additional mutagenesis studies, will be required before a firm structure/function relationship can be established.

On a more general note, our results may be considered a bit surprising. We have converted a hydrogen bond between a neutral acceptor (Q51) and a cationic donor (R45) to one between an anionic acceptor (E51) and a cationic donor (R45), a salt bridge. The quantitative contribution of salt bridges to protein stability remains controversial. The 45···51 interaction is in a region of the protein that is highly exposed to the aqueous medium (Fig. 1 A), and often a salt bridge is considered to be fairly ineffective in this environment (Dao-pin et al., 1991; Hendsch and Tidor, 1994). Other studies have shown that charge···charge side-chain interactions are slightly stronger and require less orientational specificity than charge···neutral interactions (Scholtz et al., 1993). As well, in this system there is the possibility of a double hydrogen bond in the R45···E51 interaction that is not possible with R45···Q51, and in fact over 65% of the R45···E51 interactions include such double interactions.

It might seem counterintuitive that strengthening an inter-subunit interaction in the closed state should make gating easier, as implied by a GOF mutation, but perhaps a strengthening of inter-subunit interactions in this region of the protein is a component of the structural change associated with gating. Of course, these are early studies, and we have much to learn about gating. In addition, protein-lipid interactions might change as a result of the mutation, as discussed in the following section.

Protein-lipid interactions

MscL is gated by membrane tension, and because purified MscL is fully functional when reconstituted alone into lipid vesicles it is thought that tension is transduced directly from the lipid molecules to the protein (Häse et al., 1995). A more complete description of protein-lipid interactions may give insight into this poorly understood tension transduction process. The total numbers of hydrogen bonds between the lipid and the extracellular loop, C-terminal region, and the whole protein are shown in Fig. 11 A. The total hydrogen bonding between protein and lipid steadily increases over the simulation, reaching a total of ~60 hydrogen bonds by the end of the run. This hydrogen bonding seems extensive, compared with the less than 10 hydrogen bonds between lipid and hexameric alamethicin channels (Tieleman et al., 1998a), which are known to exhibit some mechanosensitive properties (Opsahl and Webb, 1994). However, the POPC lipid in those simulations had a different headgroup than that used here, and alamethicin is comprised primarily of TM helices that do not cross the headgroup interface as extensively as Tb-MscL. The majority of the Tb-MscL hydrogen bonding occurs with the C-terminal region of the protein, which contains multiple polar and charged residues; relatively few hydrogen bonds were formed between POPE and the extracellular loop region.

View larger version (30K): [in this window] [in a new window]   FIGURE 11   (A) Total number of hydrogen bonds between POPE and protein in the wild-type Tb-MscL simulation. The numbers of hydrogen bonds involving only the extracellular loop or the C-terminal region of the channel are also shown. (B) The percentage of hydrogen bonds between POPE and protein that involve PE ammonium groups (-NH) in the wild-type Tb-MscL trajectory. These percentages are given for both hydrogen bonds to the full protein and for hydrogen bonds to only the extracellular loop or C-terminal regions. (C) The percentage of the hydrogen bonds and close contacts (less than 0.25 nm) between POPE and the extracellular loop that involve only the first eight residues of the loop.

About half of the total protein-lipid hydrogen bonds involve the ammonium (NH) of the PE headgroup (Fig. 11 B), and the percentage involving the ammonium group is even higher, around 70%, for the extracellular loop region, although there is considerable variability in this region. Thus, a large proportion of protein-lipid hydrogen bonds include a hydrogen donor from a very small part of the lipid molecules. It is also notable that 73% of the Coulombic interaction between protein and lipid in the wild-type trajectory arises from interactions including the PE ammonium group (Table 1). Other common lipid headgroups, such as phosphatidylcholine and phosphatidylglycerol, do not have a hydrogen-bond-donating group analogous to the PE ammonium and would therefore be expected to show different patterns of hydrogen bonding with embedded Tb-MscL. Such differences in hydrogen-bonding interactions could contribute to the differences in Tb-MscL gating tension observed for channels in membranes with different lipid compositions (Moe et al., 2000). Additional simulations using lipids with different headgroups would give further insight into lipid-protein hydrogen bonds and their effect on protein dynamics.

In addition to the alterations in inter-subunit hydrogen bonding discussed above, the Q51E mutation could also lead to a GOF phenotype by altering protein-lipid interactions. The Q51E trajectory showed no notable differences from the wild-type simulation in protein-lipid hydrogen bonding. However, some differences were apparent in average Coulombic interactions between the protein and lipid. Although the negatively charged Glu 51 mutant residues had an increased interaction with the positively charged PE ammonium groups, they showed a reduced electrostatic interaction with the overall headgroup due to unfavorable interactions with the phosphate group (Table 1). Simulations utilizing more sophisticated methods, such as PME, for calculating long-range interactions, would be necessary to more fully characterize any electrostatic differences such as these, and the significance of the difference for channel dynamics is unclear in the present simulations. Nonetheless, such changes in protein-lipid interactions could be plausible mechanisms by which loop mutations could lead to altered channel physiology.

The distribution of protein-lipid hydrogen bonds in the extracellular loop region is also notable. Around two-thirds of the hydrogen bonds between POPE and the extracellular loop involve the first eight residues (44-51), the first third, of the loop (Fig. 11 C). This distribution does not merely reflect an increased number of protein-lipid contacts for these particular loop residues, as only about one-third of the loop-lipid atomic contacts within the hydrogen-bonding distance cutoff occur in this section of the loop. The sequence of this portion of the Tb-MscL loop directly after TM1 is quite different from that in other non-mycobacterial MscL homologs, such as Eco-MscL. Moreover, a group of MscL homologs, including Eco-MscL, which have similar gating tensions in E. coli spheroplasts (Moe et al., 1998), and thus lower gating tensions than Tb-MscL (Moe et al., 2000), share a highly conserved region directly after TM1 (Maurer et al., 2000). The lipid hydrogen bonding of Tb-MscL channels mutated to include part or all of this conserved Eco-MscL sequence would be interesting to investigate. Along with the particularly high ratio of hydrogen bonds involving PE ammonium groups in the extracellular loop, the localization of Tb-MscL hydrogen bonds in the loop is particularly interesting in light of experimental implications of the loop region in attenuation of the Eco-MscL gating tension (Ajouz et al., 2000).

It should be noted that the positioning of the embedded protein in the membrane could affect protein-lipid interactions observed during these trajectories. However, the increased number of hydrogen-bonding interactions during the simulation implies that the interactions represent favorable conformations sampled during the MD run and not just close contacts in the first frames of the simulation. As noted in previous analyses of protein-lipid hydrogen bonds, increased simulation times would help ensure more complete sampling of relevant side-chain conformations (Forrest et al., 1999). Such extended simulation trajectories may be particularly interesting to consider for this system, as the lipid-protein hydrogen bonding equilibrated fairly slowly and was perhaps not fully equilibrated in these trajectories. In addition, it would be interesting to consider how applying tension to a system would affect the observed hydrogen bond distributions. Nonetheless, these initial analyses highlight intriguing trends in interactions between the Tb-MscL channel and its membrane environment that warrant further investigation.

Particularly strong interactions between Tb-MscL and lipid tails could also occur to aid in tension transduction. Thus, hydrocarbon tail ordering could be different for lipids bordering the channel and those in the bulk lipid. To probe for such ordering in the simulation, the deuterium order parameter (S_CD) was calculated for the saturated palmitoyl tails of the lipid. In this analysis, lipids that had at least one atom with an average minimum distance of 0.35 nm from the protein over the last 1 ns of the simulation were considered to be bordering the protein; all other lipids were defined as bulk. This cutoff was set to include all of the closest shell of lipids around the channel while omitting a maximal number of lipids outside of this shell (Fig. 12 A). The calculated S_CD values show that the channel has some disordering effect on the bordering lipids relative to bulk (Fig. 12 B). A similar but larger decrease in S_CD was also seen for OmpF, a non-mechanosensitive channel, in simulations with a POPE bilayer (Tieleman et al., 1998b). Other helical bundle channel models also have shown analogous trends in computed S_CD (Husslein et al., 1998; Tieleman et al., 1998b). As well, there was essentially no difference in the trans/gauche fraction for tail dihedrals between bordering and bulk lipids. Thus, although the lipid tails are energetically important in channel-lipid interactions (Table 1), it appears that Tb-MscL is not causing unusual lipid tail properties.

View larger version (40K): [in this window] [in a new window]   FIGURE 12   (A) Visual depiction of the cutoff between protein bordering lipids (dark gray) and bulk lipids (light gray). The protein is shown as a ribbon structure. The definition of bordering and bulk lipids is given in the text and selects all lipids in the first shell around the protein. (B) Deuterium order parameters (SCD) for protein bordering and bulk lipids in the wild-type Tb-MscL simulation.

Wild-type MscL simulation: pore water properties

Although the channel certainly does not sample an open state, or even an intermediate to an open state, during the simulations, consideration of waters in the pore region can give some insight into closed channel dynamics. As shown in Fig. 13, water in the closed Tb-MscL channel orients distinctly along the pore axis, with water O atoms pointing toward the intracellular region in the main pore and reversing their orientation in the wider region near the top of the membrane. This ordering most likely is driven by the helix dipoles of the TM1 helices, as the water oxygens point toward the helix N-terminus inside the pore and begin to reverse as they reach the extracellular side of the pore and pass the C-termini of the helices. The magnitude of dipole orientation is the greatest in the most constricted part of the channel, around residue 21, where the fewest waters would have to be ordered to obtain this effect. However, the magnitude of ordering is not greatly affected by the increase of pore radius in the V21A simulation (Fig. 13). Similar trends of water dipole ordering explained by helix dipoles are seen in simulations of alamethicin channels, the OmpF porin, and KcsA channels (Guidoni et al., 2000; Tieleman and Berendsen, 1998; Tieleman et al., 1999b). In addition, charged TM1 residues, such as Asp 36, could play a role in water ordering. Such a distinct water orientation in the MscL pore is perhaps surprising because it is not an ion-selective channel, but such electrostatic ordering would likely decrease in the presumably much wider open channel.

View larger version (26K): [in this window] [in a new window]   FIGURE 13   Average orientation of pore water dipoles along the membrane normal (z axis) over the last 1 ns of the wild-type and V21A simulations. Numbering on the z axis begins with the intracellular side of the simulation box at 0 nm. Positive dipole values result from the water O pointing toward the intracellular region (Z-position of zero); negative dipoles are directed in the opposite direction. The black vertical lines mark the average position of the phosphorous atoms of each lipid leaflet, and the average position of the plug residue 21 is denoted.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Calculations involving only the intracellular C-terminal region of Tb-MscL were used to investigate the effect of different side-chain protonation states on the stability of the C-terminal helix bundle. Conservation of the crystal structure conformation was greatly increased by neutralizing the acidic Asp and Glu residues in this region, which would occur under low pH conditions such as those used for crystallization. Several mutations to this region were made in attempts to replicate the effects of low pH under physiological pH protonation states. Although single mutations did not have much effect on the region, an E104Q/D108N double mutant afforded a helix bundle conformation equivalent to that of the low pH model. It will be interesting to experimentally consider whether this mutant does induce a C-terminal helix bundle structure. As well, any physiological effects of inducing helical stabilization are unknown, and their determination would greatly increase the current understanding of the role played by the C-terminal region in channel function. In addition, it is unclear whether other MscL homologs, such as Eco-MscL, can form similar bundled conformations in their C-terminal regions.

Simulations of the embedded channel showed the applicability of such computational methods to the Tb-MscL crystal structure. Wild-type and mutant trajectories showed an overall conservation of the crystal structure conformation. As expected, the TM regions show very little deviation from the crystal structure, whereas the extracellular loop and C-terminal regions are much more flexible. Simulations of two GOF mutants, V21A and Q51E, showed no large changes in the global structure of the protein. It is possible that in a longer-timescale simulation more significant changes would develop, but there is no hint of any such alterations in the over 1 ns of dynamics on each mutant presented here. Nonetheless, these mutations did lead to some interesting, more subtle changes in the trajectories. In particular, the Q51E mutation caused a dramatic alteration of inter-subunit hydrogen bonding between residues 45 and 51. This mutation also caused an alteration in protein-lipid electrostatic interactions. One or both of these changes could be a cause of the functional change associated with this GOF mutant, suggesting that further studies of mutations in this region could be interesting. The V21A mutation caused a widening of the constricted region of the pore. Although partially due to the steric reduction in the pore plug, the widening does extend past the close vicinity of residue 21. Although the widening was not accompanied by analogous structural deformation or instability in the pore region, it could be part of the perturbation leading to the lower gating tension of that mutant.

Several interesting aspects of protein-lipid interactions were apparent in the calculations. Half of the hydrogen bonds between protein and lipid, and two-thirds of those between extracellular loop and lipid, involved the ammonium moiety of the PE headgroup, which represents an extremely small portion of the entire POPE lipid. Additionally, this functionality would not be present in other common lipids, implying that hydrogen bonding between Tb-MscL and lipid would be very dependent on the types of lipids present. This hydrogen-bonding pattern may be related to differences in gating tension that have been observed for Tb-MscL in differing lipid environments, such as E. coli spheroplasts versus azolectin membranes (Moe et al., 2000). An analogous dependence on lipid composition also has been observed for Eco-MscL gating tensions (Sukharev et al., 1993), so an increased understanding of lipid-specific interactions may be of general interest for MscL homologs. It is also noteworthy that protein-lipid hydrogen bonding was highly localized to a particular region in the extracellular loop. In Tb-MscL, this region directly after TM1 has a very divergent sequence from Eco-MscL and many other non-mycobacterial MscL homologs (Maurer et al., 2000).

Thus, this series of simulations has given insight into the overall dynamics and conformational stability of wild-type Tb-MscL and two mutant forms of the channel. In addition, they suggest potential directions for future experimental and theoretical investigations of Tb-MscL, such as more complete considerations of protein-lipid interactions, inter-subunit interactions between residues 45 and 51, and structure-function relationships in the C-terminal region.