Gating of the Mechanosensitive Channel Protein Mscl: Protein-Induced Membrane Elasticity

1

Supplementary Material

Gating of the mechanosensitive channel protein MscL: the interplay of membrane and protein

Jonggu Jeon and Gregory A. Voth

Nonequilibrium molecular dynamics (NEMD)simulation

In the NEMD simulation, the applied strain is homogeneously reflected in the equations of motion of the individual atoms as well as in the system size. The NEMD algorithm in the authors’ previous study(1, 2)was integrated with the pressurecontrol of Hoover type (3). However, the current study uses the GROMACS program which does not provide a correct implementation of the Hooveror Parrinello-Rahman barostat. Thus, a differentequation of motion is employed here for use with the Berendsen barostat as follows

(1)

where mi, ri, pi are the mass, position, and momentum of particle i,rCM is the center of mass of the system in thecentral simulation box,ei (i=1,2,3) is the unit vector in each Cartesian direction, is the applied strain rate, is the system compressibility in the z direction, is the barostat relaxation time, and are thezzcomponents of the current and target pressures, respectively, and Lx,y,z are the box size in each Cartesian direction. The terms involving the pressure arise from the Berendsen barostat. Note that the momentum equation is not modified by the applied strain to be consistent with the Berendsen barostat, which does not scale velocities. The NEMD with constant is available in the GROMACS program since version 3.3b. Thus, Eq. 1 with stepwise strain (see Eq. 2 of paper) can be implemented with modest modification of the “deform” routine in GROMACS. In practice, the stepwise strain was approximated by the following function

(2)

witha = 248.58 ps-1 and b = 0.025 ps. With these parameters, Eq. 2 closely resembles the stepwise strain with a rise time of 2b = 0.05 ps.

Despite the differences in the MD package and pressure control algorithm as detailed above, the calculated relaxation function of the POPE membrane after area expansion (Table 5 in the paper) is very similar to that of the DMPC membrane in our previous study (Table 1 of Ref. 1). Specifically, both systems show (i) sub-ps responses with amplitude larger than 10 kbar, (ii) slower and weaker responses with tens of ps relaxation time, and (iii) much weaker responses with relaxation times in the hundreds of ps to ns range. These different modes have been attributed to (i) the inter- and intramolecular vibrations of water and lipids, (ii) acyl chain isomerization, (iii) rotations of lipid head groups and the entire lipid about its long axis, respectively, in Ref. 1. We note that the static relaxation function from the previous DMPC study (–0.704 kbar: Table 1 of Ref. 1) carries uncertainty of as large as 1.7 kbar. However, thanks to the longer simulation timeand the more reliable determination of the equilibrium lateral pressure from separate NPAT simulations on the same system (see Computational Methods in the paper), the current study yields a more reasonable of 1.517 kbar (Table 5 in the paper).

Preparation and equilibration of the POPE lipid bilayer

A single POPE lipid structure obtained from Tieleman and Berendsen (4) was replicated to make a bilayer composed of 392 (=14×14×2) lipids with the area of 12544 Å2 (112 Å in the x and y directions) and the thickness of 48 Å measured between the locations of phosphorus atoms. 7281 SPC water molecules were then added to the space right above and below this lattice structure, giving the initial system thickness of 70 Å. The water molecules were not allowed to fill the space between the phosphorus atoms of the two monolayers. This initial system (called PW0; see Fig. S1) was first energy-minimized with the steepest descent method to remove bad contacts between atoms andthen equilibrated for 0.5 ns under the constant NPT condition (at 1 bar and 308 K with the Berendsen semiisotropic barostat and the Nose-Hoover thermostat, respectively). After additional 0.08 ns of NPT dynamics with Parrinello-Rahman barostat with anisotropic scaling at 1 bar, the system shrank to 105.8, 106.2, 59.4 Å in the x,y,z directions, respectively. Additional water molecules were added to this system to make the thickness 70 Å again. This system (referred to as PW1 hereafter) had 392 POPE lipids and 10160 water molecules for a total number of atoms of 50864 and the system size 105.8×106.2×70 Å3.

After system PW1 was equilibrated for ~0.7 ns (the Parrinello-Rahman barostat at 1 bar, the PME electrostatics with the real-space cutoffof 10 Å, and the Lennard-Jones (LJ) cutoffof 10 Å were used), the two monolayers were observed to separate from each other, as reported in the paper. This indicates insufficient attraction between hydrocarbon tails of the two monolayers. Thus, the simulation was run again at the 0.6-ns point withincreased to 15 Å. With this change, the integrity of the bilayer structure was maintained during the next 17.3 ns of simulation (17.9 ns since the preparation of system PW1).

The GROMACS lipid force field used in this study (5) was originally developed on a pentadecane system with of 10 Å and the electrostatic interaction cut off at 18 Å. We have separately tested these parameters on thePOPE system and it resulted in the unstable bilayer. Thus, the instability of the POPE system with = 10 Å could be attributed to the presence of an unsaturated C-C bond in POPE which tends to result in less packing of hydrocarbon tails. Apparently, the parameters developed on saturated chains are not transferable to an unsaturated system. Considering this and the additional computational cost associated with the use of large , we tried = 12 Å which preserved the bilayer integrity during a 9-ns long simulation starting with an equilibrated system at = 15 Å. Consequently, all production computations were performed with = 12 Å. Specifically, four POPE configurations after 7.3, 8.3, 9.3, and 10.3 ns of equilibration with = 15 Å were taken and further equilibrated at = 12 Å for 1 ns. These four configurations, hereafter referred to as PWE1 to PWE4, were used in the NEMD and NPAT simulations of the POPE system (see Methods in the paper). Separately, the MscL-POPE system was prepared with the 4.3-ns configuration of the POPE bilayer (system PWM0).

The area per lipid and the membrane thickness of the systems PWE1-PWE4 are summarized in Table S1. For comparison, the four trajectories leading to PWE1-PWE4 were continued,each for 4 more ns, and averages from these longer simulations (20 ns in total) are also reported in Table S1. The thickness was determined as the average distance along the membrane normal between phosphorus atoms in the two monolayers. Fig. S2 shows the deuterium order parameter SCD of the palmitoyl chains of the POPE bilayer. Since the united atom model was used, |SCD|was determined from the orientation of C-C bonds assuming tetrahedral configuration (6). The calculated area per lipid of POPE bilayer is somewhat smaller than a previous simulation result of 53.3 Å by Gullingsrud and Schulten (7) but these authors expected the value to decrease on longer equilibration. Table S1 and Figure S2 show good agreement between single configuration values from PWE1-PWE4 and the averages over 20 ns trajectory, indicating a good equilibration of the starting structures for the NEMD simulation.

Table S1. Properties of the POPE bilayers used in the NEMD simulations. The entries PWE1-PWE4 are single configuration properties and the entry “Average” is the average over all four trajectories that generated PWE1-PWE4, each 5-ns long.

Area per lipid (Å2) / Thickness (Å)
PWE1 / 49.7 / 44.5
PWE2 / 49.2 / 45.1
PWE3 / 49.4 / 44.9
PWE4 / 49.3 / 45.0
Average / 49.5 / 44.7

Figure S1. Equilibration scheme of the POPE system. Key configurations defined in text are enclosed in rectangular boxes. See text for details.

Figure S2. Deuterium order parameter |SCD| of the palmitoyl chains of the POPE bilayer systems. The carbon numbering starts at the carbonyl carbon of the acyl group.

Figure S3. Changes in (a) the lateral pressure and (b) lateral dimension during simulations with negative lateral pressure. Trajectories NSB1 and NSB2 have target lateral pressure of -100 bar. NSB1 employs barostat relaxation time of 50 ps while NSB2 uses 1 ps. NSB2 was started with 5 ns configuration of NSB1. See Table 1 of main text for further details.

Figure S4. Changesin the radius of gyration of MscL components during the simulations with negative lateral pressure. The components in each plot are defined in Fig.1 of main text. Trajectories NS-NSF have target pressure bar. NSB1 and NSB2 have bar. The trajectories are described in Table 1. Note that NSB2 starts at 5 ns of NSB1.

Figure S5.Pore radius profiles of MscL. (a) Radii of the crystal structure (2OAR) andtwoconfigurations at the end of equilibrium simulations F4 and SF2 under zero tension. (b) Radii of the crystal structure andthree configurations at the end of trajectories NS, NSE, NSF. NS-NSF were run with = -1000 bar. Trajectories are defined in Table 1of the main text. Note that simulations SF2 and NS-NSF were carried out with truncated MscL and thus no pore is defined below ~70 Å of the vertical axis. The HOLE program(8) was used in the calculation.

Figure S6.Changesin the radius of gyration of MscL components during the simulations with radial bias force. The components in each plot are defined in Fig. 1 of main text. Trajectories RS-RSF have target pressure bar. The trajectories are described in Table 3 of main text.

References

1.Jeon, J., and G. A. Voth. 2005. The Dynamic Stress Responses to Area Change in Planar Lipid Bilayer Membranes. Biophys. J. 88:1104-1119.

2.Ayton, G., A. M. Smondyrev, S. G. Bardenhagen, P. McMurtry, and G. A. Voth. 2002. Calculating the Bulk Modulus for a Lipid Bilayer with Nonequilibrium Molecular Dynamics Simulation. Biophys. J. 82:1226--1238.

3.Melchionna, S., G. Ciccotti, and B. L. Holian. 1993. Hoover NPT dynamics for systems varying in shape and size. Mol. Phys. 78:533-544.

4.Tieleman, D. P., and H. J. C. Berendsen. 1998. A Molecular Dynamics Study of the Pores Formed by Escherichia coli OmpF Porin in a Fully Hydrated Palmitoyloleoylphosphatidylcholine Bilayer. Biophys. J. 74:2786-2801.

5.Berger, O., O. Edholm, and F. Jahnig. 1997. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J. 72:2002-2013.

6.Egberts, E., and H. J. C. Berendsen. 1988. Molecular dynamics simulation of a smectic liquid crystal with atomic detail. J. Chem. Phys. 89:3718.

7.Gullingsrud, J., and K. Schulten. 2004. Lipid Bilayer Pressure Profiles and Mechanosensitive Channel Gating. Biophys. J. 86:3496-3509.

8.Smart, O. S., J. G. Neduvelila, X. Wanga, B. A. Wallacea, and M. S. P. Sansom. 1996. HOLE: A program for the analysis of the pore dimensions of ion channel structural models. J. Mol. Graph. 14:354-360.