This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Knecht, V. Articles by Grubmüller, H. Search for Related Content PubMed PubMed Citation Articles by Knecht, V. Articles by Grubmüller, H.

Mechanical Coupling via the Membrane Fusion SNARE Protein Syntaxin 1A: A Molecular Dynamics Study

Theoretical Molecular Biophysics Group, Max-Planck-Institute for Biophysical Chemistry, 37077 Göttingen, Germany

Correspondence: Address reprint requests to Helmut Grubmüller, Tel.: +49-551-201-1763; Fax: +49-551-201-1089; E-mail: hgrubmu{at}gwdg.de.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
SNARE trans complexes between membranes likely promote membrane fusion. For the t-SNARE syntaxin 1A involved in synaptic transmission, the secondary structure and bending stiffness of the five-residue juxtamembrane linker is assumed to determine the required mechanical energy transfer from the cytosolic core complex to the membrane. These properties have here been studied by molecular dynamics and annealing simulations for the wild-type and a C-terminal-prolongated mutant within a neutral and an acidic bilayer, suggesting linker stiffnesses above 1.7 but below 50 x 10^-3 kcal mol^-1 deg^-2. The transmembrane helix was found to be tilted by 15° and tightly anchored within the membrane with a stiffness of 4–5 kcal mol^-1 Å^-2. The linker turned out to be marginally helical and strongly influenced by its lipid environment. Charged lipids increased the helicity and H3 helix tilt stiffness. For the wild type, the linker was seen embedded deeply within the polar region of the bilayer, whereas the prolongation shifted the linker outward. This reduced its helicity and increased its average tilt, thereby presumably reducing fusion efficiency. Our results suggest that partially unstructured linkers provide considerable mechanical coupling; the energy transduced cooperatively by the linkers in a native fusion event is thus estimated to be 3–8 kcal/mol, implying a two-to-five orders of magnitude fusion rate increase.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
For exocytosis and for the selective transport of macromolecules between the various organelles of eukaryotic cells, the fusion of a transport vesicle membrane with a target membrane is an essential step (Alberts et al., 2002). A complex network of sequentially interacting proteins is involved in the tethering and docking of the vesicles, as well as the promotion and regulation of the fusion process (Misura et al., 2000; Mochida, 2000). For synaptic transmission, the highly conserved SNARE (soluble N-ethylmaleimide sensitive factor attachment protein receptor) proteins (Söllner et al., 1993) have been identified as central players and studied in detail (Hanson et al., 1997; Jahn and Südhof, 1999; Chen and Scheller, 2001), and specific sets of SNARE homologs were also found for a number of other transport pathways (Bennett and Scheller, 1993; Ferro-Novick and Jahn, 1994). Although the precise role of SNARE proteins in membrane fusion is still a subject of ongoing discussions (Ungermann et al., 1998; Tahara et al., 1998; Mayer, 1999; Brünger, 2001; Bruns and Jahn, 2002; Jahn and Grubmüller, 2002), in vitro fusion experiments suggest that SNAREs form a minimal fusion machinery (Weber et al., 1998; Nickel et al., 1999). According to the zipper model of membrane fusion (Hay and Scheller, 1997; Hanson et al., 1997; Bruns and Jahn, 2002), vesicle (v) and target (t) membrane-anchored SNAREs assemble in a zipperlike fashion into extended trans complexes (see Fig. 1). Because both ends are anchored within the membranes via transmembrane (TM) domains (tilted yellow and magenta rectangles), the membranes are pulled into close proximity, thus promoting fusion.

View larger version (21K): [in this window] [in a new window]   FIGURE 1  Model for SNARE-promoted membrane fusion. A trans complex is formed by synaptobrevin and SNAP-25 (blue), the H3 helix of syntaxin (cyan), a linker region (red), the transmembrane (TM) domains (tilted rectangles) of syntaxin (yellow), and synaptobrevin (magenta), which exerts mechanical force onto the vesicle (v) and the target (t) membrane and pulls them toward each other. The H3 helix tilt (pink) is studied in molecular dynamics simulations; the brown box marks the simulation system. The hydrophobic regions of the two membranes are colored green, their polar regions dark green.

It has been suggested (Weber et al., 1998) that formation of the trans complex induces significant strain within the membrane-proximate regions of synaptobrevin (magenta) and syntaxin (red, linker), respectively. Thus, and assuming that the TM helices remain perpendicular to the local membrane plane, energy provided by complex formation could, via a conformational coupling, be transmitted to the bilayers and thereby induce membrane bending as shown in Fig. 1. Evidently, the coupling efficiency critically depends on the mechanical bending stiffness of the membrane-proximate regions. In particular, the energy available for membrane juxtaposition is limited by the amount required to increase the H3 helix tilt angle (pink) from its equilibrium value to the one required for SNARE complex assembly.

For synaptobrevin, the -helical domain of the synaptic fusion complex x-ray structure (core domain) extends into the TM domain by two residues; this suggests that core and TM domain form a continuous -helix, which is expected to be stiff upon bending (Shaitan et al., 2002). Membrane bending and final merger allowed the synaptobrevin helix to gradually relax to an unbent state (McNew et al., 1999).

For syntaxin, the H3 helix is separated from the TM helix by a five-residue basic linker (red) of unknown structure. The only information comes from recent electron paramagnetic resonance (EPR) measurements, which provided some evidence for an unstructured linker (Kweon et al., 2002). Hence it is unclear if a similar mechanism holds here. In fact, inserting the five-residue sequence GPPKL containing two helix-breaking proline residues into the linker of syntaxin reduced the in vitro fusion rate by only 50% (McNew et al., 1999). Similarly, for a syntaxin homolog in yeast, a six-residue insertion between the expected core and the TM domain led to a comparable fusion rate decrease in vivo. The effect of the proline insertions on the linker stiffness is, however, unclear. Although disfavoring -helical conformations, it can also rigidify the protein backbone (Schultz and Schirmer, 1979). Moreover, a similar proline mutation of synaptobrevin had virtually no effect on fusion efficiency (McNew et al., 1999). These results suggest that either a stiff linker is not required for fusion, or sufficient conformational coupling can also be provided by a nonhelical linker.

To discriminate between these two possibilities, we have performed multinanosecond molecular dynamics (MD) simulations of the relevant region (Fig. 1, box; see also Fig. 6) comprising a 39-residue C-terminal part of syntaxin (cyan, red, and yellow) embedded within a lipid bilayer (green) and explicit solvent environment (in Fig. 6 water molecules are shown in blue). In particular, the following questions were addressed: 1), What is the native structure of the five-residue syntaxin linker? Is it -helical? 2), Can a nonhelical linker also provide substantial conformational coupling? 3), Can the bending stiffness be estimated from the simulations? 4), Which factors (e.g., lipid composition, environment, sequence) determine the mechanical properties of the linker? 5), Is the TM domain sufficiently perpendicular to the membrane plane, as required for conformational coupling? 6), Is the TM domain anchored well enough within the membrane to withstand the high intermembrane repulsion forces due to hydration and electrostatic interactions (Lipowsky, 1995), which tend to pull the TM domain out of the membranes?

View larger version (113K): [in this window] [in a new window]   FIGURE 6  Simulation system setup. Shown is the C-terminal segment of syntaxin, comprising an 11-residue segment of the H3 helix (Asp-250–Lys-260, cyan), the five-residue basic linker (Ala-261–Lys-265, red), and the 23-residue TM helix (Ile-266–Gly-288), embedded in a lipid bilayer (green, the green spheres indicate the phosphor, oxygen, and nitrogen atoms of the polar regions), with an explicit solvent environment (water molecules are shown in blue; ions are not shown). The magenta rectangle indicates the position of the 10-asparagine prolongation present in some of the simulations.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Room-temperature and annealing simulations were carried out for various systems. From these simulations, free-energy landscapes for relevant observables, such as helix tilts and vertical peptide position within the membrane, were obtained.

Room temperature simulations
Simulation systems
MD simulations were carried out for three peptides (compare Table 1): the wild-type (WT), a WT peptide with a truncated TM domain (WT^S), and mutant prolonged by 10 asparagine residues (N₁₀) described below.

Table 2 specifies the two lipid compositions that were chosen to study the influence of the lipid environment on the peptide properties. A neutral bilayer (POP⁰) was created with a binary lipid mixture of zwitterionic palmitoyl-oleoyl-phosphatidylcholine (POPC) and palmitoyl-oleoyl-phosphatidylethanolamine (POPE) molecules, and a more physiological acidic bilayer (POP-) with a ternary lipid mixture of POPC, POPE, and palmitoyl-oleoyl-phosphatidylserine (POPS) molecules. Table 3 summarizes all simulated systems that fall into three sets. Set A contains the four possible combinations of the peptides WT and N₁₀ with bilayers POP⁰ and POP^- and was used for the main production runs at T = 300 K. Set B was used for additional annealing runs. Set C was used to obtain first insight into membrane effects, e.g., by simulating a simple van der Waals membrane (VdW, see below). The WT systems comprised 20, 000 atoms, including 373 peptide atoms, and the N₁₀ systems comprised 22, 000 atoms, including 483 peptide atoms. All simulation systems with explicit lipid environment were set up in rectangular boxes with periodic boundary conditions and lateral dimensions a x b of 60 Å x 60 Å. The box sizes in the z direction (84–88 Å) were chosen such as to obtain a minimum distance of 7 Å between the peptide-membrane system and the boundary in z direction, yielding a water layer of 30 Å between the membranes.

Solvation methods
Systems were solvated by filling the periodic box with water molecules and removing all those water molecules that overlapped with the peptide or lipid atoms. For the explicit lipid bilayer, also all water molecules within the hydrophobic core of the bilayer were removed. Ions were added by replacing water molecules at energetically favorable positions (Coulomb criterion), or in shells of radius 10 Å around oppositely charged groups (Debye criterion). If not denoted otherwise, the type (Cl^- or Na⁺) and number of ions were chosen so as to compensate the net charge of the peptide-membrane system (charge compensation). In one case, excess ions were added and K⁺ rather than Na⁺ ions were used to mimic the intracellular conditions of typical mammalian cells (Alberts et al., 2002), and, in particular, to obtain a physiological concentration of [KCl] = 0.154 M, (intracellular salt conditions).

Pilot studies: absent and simplified lipid environment
To check whether the lipid environment is relevant for the linker structure and, therefore, should be included within the production simulations, two simulations of 2 ns length each without explicit lipid environment were carried out (Table 3, Set C). The first simulation involved no lipid bilayer at all, the second a simplified bilayer model. Accordingly, a peptide with a truncated WT^S could be used. The conformations of two linker residues (Ala-261 and Arg-262) obtained from these simulations were used to model a linker start structure that was improved in subsequent simulations with explicit lipid environment.

Fig. 2 helps the reader to keep track of the various simulations and their interdependencies both for the pilot runs as well as, in particular, for the subsequent equilibrium and annealing runs. To obtain a start structure (see Fig. 3, top left), Protein Data Bank (PDB) data (where available), and modeling were combined (see Fig. 2, green boxes). The H3 helix was taken from chain F of the x-ray structure (Sutton et al., 1998), PDB entry 1sfc (green arrows from PDB box). The TM domain was modeled as a right-handed -helix using Insight II (Molecular Simulations, Cambridge, UK). The linker was modeled with its dihedral angles chosen from the ß-region of the Ramachandran plot (Ramachandran et al., 1963) (for the side-chain dihedrals, the default settings of Insight II were used). Because in the pilot study the peptide was truncated at both ends, both termini were modeled neutral (NH₂ and COOH, respectively). Both modeling procedures are summarized by the respective green box and arrows.

View larger version (32K): [in this window] [in a new window]   FIGURE 2  Overview of simulation runs and their mutual interdependencies. Blue (or yellow) and red boxes denote room-temperature runs at T = 300 K (also tagged 300K), and annealing sets (ANNEAL), respectively; green fonts indicate simulations with intracellular salt conditions. Green arrows denote input from the Protein Data Bank or major modeling steps; black and gray arrows indicate that a final configuration has been used as a starting configuration for a subsequent simulation, either without (black) or with (gray) modification. For nomenclature and specification of the various simulation systems, see Table 3.

View larger version (35K): [in this window] [in a new window]   FIGURE 3  Starting peptide structures (left) and intermediate structures (right) of the two pilot simulations described in the text, without lipid environment (top), and with simplified lipid environment (bottom). In the top right figure, the spheres show the residue pair forming a hydrophobic contact, which stabilized the final structure of the peptide, i.e., the linker residue Ala-261 (orange) and the TM residue Ile-270 (yellow). In the bottom panel, the position-restrained water oxygen atoms that mimic the sterical effect of a bilayer are shown as red spheres. In the bottom right figure, the black circle marks the two linker residues that spontaneously folded from loop into -helical conformation.

The system WT^S/VdW (Fig. 3, bottom left) was set up in a similar manner. To mimic the reduced conformational space available for the peptide due to the presence of the lipid bilayer, the oxygen atoms of all water molecules in a 5-Å layer perpendicular to the TM helix (red spheres) were subjected to position restraints. Additionally, position restraints were imposed onto the backbone atoms of the TM helix. The linker conformation was changed such that the side chains of the linker and the H3 helix did not penetrate the membrane layer.

Production runs: explicit solvated lipid environment
In the simulation, WT^S/VdW/300 K, after 0.3 ns, two residues of the linker, namely Ala-261 and Arg-262, adopted an -helical conformation (Fig. 3, bottom right, black circle), which remained stable for 0.9 ns. Therefore the linker structure after 1.2 ns was taken as the start structure of the first simulation with explicit lipid environment, (Fig. 2, gray arrow).

The start structures of the H3 and the TM helix were obtained similarly as for WT^S, except that three more H3 residues and the full TM domain (including Gly-288) were included. The TM helix was modeled using MOLMOL (Koradi et al., 1996); the side-chain dihedrals were set to the default values provided by MOLMOL.

For the lipid environment, an equilibrated bilayer of 128 POPC molecules kindly provided by P. Tieleman (Tieleman, 1998, 2000) was used. For the peptide insertion, a cylindrical hole of radius 7 Å was created by removing those four lipid molecules from each layer, whose headgroups were located within the hole region. Subsequently, a weak repulsive and radially acting force centered at the vacated region was applied to the lipid molecules during a short MD run to remove lipid tail atoms from the hole region (Tieleman, 1998).

As can be seen in Fig. 6, severe steric conflicts between the linker of the peptide (red) and several lipid headgroups (green), which could not be resolved by minor changes of the backbone configuration of the three remaining nonhelical residues, prohibited complete peptide insertion. As a result, the TM helix protruded at its N-terminal side into the polar region of the bilayer by 5 Å, whereas, at the C-terminus, a hole within the hydrophobic core of the bilayer of about the same length remained. To avoid closure of that hole, a prosequence of 10 -helical asparagine residues, modeled in -helical conformation, was appended at the C-terminus (N₁₀ mutant) with the expectation that the peptide would relax and shift downward during the equilibration phase to fully bury the apolar TM residues in the hydrophobic core of the bilayer, such that the prosequence could then be removed. For the prosequence, the choice of asparagines was motivated by the expectation that their high polarity (Stryer, 1988) should drive them quickly out of the hydrophobic region of the bilayer. Furthermore, the neutral asparagines should disturb the system only slightly.

Two POPC molecules of each layer still overlapped with the peptide and were removed. To obtain the desired lipid composition of the neutral bilayer, 46 uniformly distributed lipid molecules were mutated into POPE molecules. The system was then solvated and subsequently equilibrated for 5 ns the superscripts denote 1–4 distance restraints for the TM helix (T) and two linker residues Ala-261 and Arg-262 (L).

From the final structure of this trajectory, the set of the four possible combinations of lipid environment and peptide type (Table 3, Set A) was obtained as shown in Fig. 2 as follows: The system N₁₀/POP⁰ was transformed into WT/POP⁰ by removing the prosequence, modeling a new C-terminus, and filling the resulting cavity with water molecules. The two POP⁰ systems (WT/POP⁰ and N₁₀/POP⁰) were transformed into respective POP^- systems by removing the solvent environment, mutating selected POPE into POPS molecules, resolvating the system, and adding ions for charge compensation according to the Debye criterion. For each of the four obtained systems, free dynamics production runs of 5 ns length each were carried out (four blue boxes) for further analysis and computation of free-energy landscapes. Two final structures (taken from WT/POP⁰/300 K and N₁₀/POP⁰/300K) served as starting structures for subsequent high temperature runs (red boxes).

Further room-temperature simulations were performed as shown in Fig. 2 (Table 3, Set B). From the equilibrated system WT/POP^-, a system with intracellular salt conditions was created by changing the Na⁺ into K⁺ ions and inserting excess K⁺ and Cl^- ions according to the Coulomb criterion. To check for a possible bias of the chosen linker start structure, an independent line of simulations was performed. These started with simulation (second box, left column) with the linker now modeled in fully -helical conformation (subscript H). Here, no steric hindrance occurred in the polar region, and the peptide could be readily inserted, so no membrane-stabilizing prosequence had to be appended. Otherwise, the simulation system was set up similar to To equilibrate the lipid environment, the peptide was first immobilized by position restraints on the backbone atoms of the TM domain (superscript P) and, in addition to helix H3 and the TM domain, also the linker was subjected to 1–4 distance restraints (D). In subsequent simulations, position and distance restraints (WT_H/POP⁰/300 K) were successively removed.

Simulation parameters
All simulations were performed with the GROMACS simulation package (Van der Spoel et al., 1995), version 2.1, using the GROMOS-87 force field (van Gunsteren and Berendsen, 1987) with modifications (van Buuren et al., 1993). Polar hydrogen atoms were simulated explicitly; nonpolar hydrogen atoms were described via compound atoms. The SPC water model (Hermans et al., 1984) was used. The POPC force field originated from Berger et al. (1997), parameters for the unsaturated carbons and POPE parameters from the GROMOS-87 force field; the POPS force field (Knecht, Ph.D., thesis, University Göttingen (unpublished)) was derived from the POPE parameters using serine as a template. Using LINCS (Hess et al., 1997) to constrain the bond lengths, and heavy hydrogen atoms (Feenstra et al., 1999), a time step of 4 fs could be chosen. For the trajectories of the pilot study, light hydrogen atoms and a time step of 2 fs were used.

To calculate the nonbonded forces, a twin-range cutoff was used on a charge group basis: In each step, short-range electrostatic and van der Waals forces were calculated for all pairs of a neighbor list with nearest-image distances r r_list = 10 Å; the neighbor list, as well as the long-range electrostatic force on each particle due to all particles in a shell r_list < r 18 Å, were updated every 10 steps (in the pilot studies, this was done for a shell r_list < r 15 Å). For the intracellular salt condition simulations, the particle mesh Ewald method (Darden et al., 1993; Kholmurodov et al., 2000) was used. For the neighbor list and for the calculation of the van der Waals and Coulomb forces, a cutoff radius of 9 Å was used. For the calculation of the reciprocal sum, the charges were assigned to a grid in real space with a lattice constant of 1.2 Å using cubic interpolation.

All systems were minimized by a steepest descent, with an initial step size of 0.1 Å until the step size (in nm) converged to the machine precision. Peptide, each lipid type, water molecules, and salt ions were coupled separately to a heat bath of 300 K using a Berendsen thermostat (Berendsen et al., 1984) with a coupling time constant of 0.1 ps. To retain the initial square shape of the bilayer patch, pressure coupling was applied to the xy plane and the z direction separately, with a reference pressure of 1 bar and a coupling constant of 1 ps. Snapshots of the system were recorded every 0.4 ps for further analysis.

Simulated annealing
In Fig. 2, red boxes indicate groups of simulated annealing runs (Kirkpatrick et al., 1983) (subsequently referred to as annealing sets, see Fig. 4), which were performed for different solvated peptide-membrane systems. Each annealing set consists of a high temperature (1500 K) sampling run (see Fig. 4, brown curve) and multiple annealing runs (blue) and was generated as follows. An equilibrated start system was stabilized by a) position restraints on the TM helix backbone and lipid tail atoms using a force constant of 359 kcal mol^-1 Å^-2, b) 1–4 distance restraints on the H3 and the TM helices (plus the prosequence if present), analog to the distance restraints described above, with a force constant of 837 kcal mol^-1 Å^-2 and c) backbone dihedral restraints. Subsequently, for the high temperature sampling run, the peptide, each lipid type, the water molecules, and the salt ions were separately coupled to a heat bath of 1500 K using a coupling time constant of 10 ps and simulated for 1 ns. To avoid expansion of the box due to vaporization of the water, pressure coupling was switched off. Because of the faster atomic motions at high temperature, a smaller integration time step of 2 fs had to be used. As indicated by the green circles, 27 structures were chosen from the sampling trajectory, each of which was cooled down to 300 K linearly in time within 400 ps (sloped blue lines), using a coupling constant of 0.1 ps. Subsequently, water molecules that had moved into the hydrophobic core of the bilayer were removed. The resulting 27 structures were equilibrated for 200 ps each (horizontal blue lines). In these simulations, all restraints except the 1–4 distance restraints on the H3 helix were removed, an integration time step of 4 fs was used, and pressure coupling was switched on. The 27 final structures thus obtained (red crosses) were used for further analysis.

View larger version (19K): [in this window] [in a new window]   FIGURE 4  Procedure employed for each annealing set. From a 1-ns high-temperature sampling run (brown curve), 27 snapshots (green circles) were selected, annealed for 400 ps, and subsequently equilibrated for 200 ps (blue curves). The final structures (red crosses) were kept for further analysis.

Analysis methods
Free-energy calculations
Free-energy landscapes were estimated for the four different observables defined in Fig. 5. These are the H3 helix tilt and the TM helix tilt _TM with respect to the bilayer normal (z axis), the position z_TM of the center of mass of the TM helix C-atoms relative to the average z-position of the phosphor atoms of the bilayer, and the position z_link of the center of mass of the linker C-atoms relative to the average z-position of the phosphor atoms of the proximate monolayer. The helix orientations were determined from the principal axis of the respective C-atoms, i.e., from Asp-250–Lys-260 for the H3 helix, and from residues Ile-266–Gly-288 for the TM helix.

View larger version (31K): [in this window] [in a new window]   FIGURE 5  Reaction coordinates used for free-energy calculations, H3 helix tilt angle , TM helix tilt TM, TM helix position zTM, and linker position zlink. The red rectangles indicate the polar regions, the green rectangle the hydrophobic core of the bilayer.

(1)

was obtained, where x_i is the center of bin i, k_B the Boltzmann constant, and T the temperature. For the room-temperature ensembles, G(x_i) was smoothed by a Gaussian filter of the width of five bins standard deviation. Additionally, to estimate the associated stiffnesses from both the room temperature and the annealing ensembles, harmonic fits,

(2)

to the free-energy landscapes were obtained by computing average (equilibrium) values

(3)

and variances

(4)

of the observable x. Effective stiffnesses k were obtained from these harmonic fits via

(5)

Free energy of helix formation
For each linker residue i (i = 261,...,265), the probability p_i and free energy G_i of helix formation was estimated as follows: Using the criteria of Kabsch and Sander for secondary structure elements (Kabsch and Sander, 1983), the number n_i of structures was determined for which one of the neighbors of the residue i was -helical, and the number f_i of structures in which the residue i was -helical itself. The free energy G_i of helix formation was estimated via

(6)

with

(7)

Confidence intervals for annealing results
For the annealing results, due to their relatively small sample size, confidence intervals with confidence level 1 - = 95% for the p_i were estimated from the Wilson score method (Wilson, 1927),

(8)

with = ^-1(1 - /2), and denoting the error function. The upper and lower bounds of G_i, were obtained from Eq. 1 by replacing p_i by its bounds . The confidence intervals of G_i were estimated similarly, here replacing n by n_i in Eq. 8. The confidence intervals [µ_-,µ₊] of µ were estimated from (E. Brunner, lecture note, Göttingen; C. Cenker, lecture note, Vienna)

(9)

where t_{n - 1,1 - /2} is the (1 - /2) quantile of Student's t-distribution with n - 1 degrees of freedom (Rohatgi and Saleh, 2001) (defined via (Rohatgi and Saleh, 2001)). The confidence intervals of were estimated as (E. Brunner, lecture notes; C. Cenker, lecture notes)

(10)

where _n-1,/2² is the /2 quantile of the _n-1² distribution (Rohatgi and Saleh, 2001).

Tilt energies
The energy necessary to tilt the H3 helix into its presumed orientation within the trans complex was estimated from

(11)

with = 80°. This angle was chosen somewhat smaller than 90° to account for the twist of the individual strands of the four-helix bundle with respect to the bundle axis. The parameters k and µ were obtained by harmonic fits to ensembles of simulated annealing final structures according to Eqs. 3 and 5. Confidence intervals [E_-,E₊] for a confidence level of 1 - = 95% were estimated from the 1 - percentile intervals obtained from a bootstrap (Efron and Tibshirani, 1993) run. Here, 5000 replica {^*}, each equal in size to the original sample {}, were generated by sampling with replacement from {} to obtain an empirical distribution of E^* estimates. To verify that convergence was reached, varying numbers of replica were considered.

Significance of observed differences
To determine whether observed differences in estimated equilibrium tilt angles or stiffnesses are statistically significant, a bootstrap hypothesis test on the equality of means (Efron and Tibshirani, 1993) was applied. For samples {x}, 30,000 replica {x^*}, were generated by sampling without replacement from where µ denotes the mean value of the merged sample; hence an empirical distribution of the test quantity

(12)

was obtained. The fraction of t that exceeds the observed value yields the desired significance. For the equilibrium tilt angles, the test was applied to respective sample pairs {}, (bootstrap mean test); for stiffnesses, it was applied to the absolute deviations from the mean, {| - |}, (bootstrap variance test).

For differences in the helicity of the individual linker residues, a modified bootstrap test for the equality of probabilities (bootstrap probability test) was devised: For two structure samples of sizes n₁ and n₂, let the studied residue be in helical conformation in k₁ and k₂ cases, respectively, yielding probability estimates p_i = k_i/n_i, i = 1,2. With the null hypothesis that the underlying probabilities are equal, p = (k₁ + k₂)/(n₁ + n₂) provides the best probability estimate. That estimate was subsequently used to generate bootstrap samples, thus yielding empirical distributions for k₁ and k₂ and, therefore, also for the test quantity t = |p₁ - p₂|. Hence the significance was estimated similarly as above.

Statistical independence of successive tilt angles
To check to what extent the H3 helix tilt angles _i and _i+1 of successive annealing start structures (Fig. 4, green circles) can be considered statistically independent, the (auto-) correlation coefficient r,

(13)

with and was calculated. The significance of the calculated correlation coefficients was assessed by a bootstrap correlation test for the null hypothesis of uncorrelated data. Bootstrap samples were generated by randomly permuting the _i, and the significance was estimated similarly as above from the obtained distribution of the test quantity t = |r|. A similar procedure was applied to the tilt angles of successive annealing final structures (Fig. 4, red crosses).

Hydrogen bonds
For the room temperature runs and for the final structures of each annealing set, the average number N_HB = n_HB of hydrogen bonds between the peptide (H3 and linker residues) and lipid headgroups was determined, as well as its fluctuation,

(14)

Here, n_HB denotes the the number of hydrogen bonds for each snapshot. Hydrogen bonds were counted if the distance between the hydrogen atom and the acceptor was below 2.5 Å and the angle between hydrogen atom, donor and acceptor was below 60°.

Additionally, for each set of n annealing final structures, the hydrogen bond probability p_i = k_i/n of each donor and acceptor i within the H3 helix and the linker was estimated by counting the number of structures k_i for which the respective hydrogen bond is detected. The significance of differing p_i between two annealing sets was estimated using the bootstrap probability test.

Figs. 3, 6, 14, and 15 were prepared using MOLMOL (Koradi et al., 1996).

View larger version (45K): [in this window] [in a new window]   FIGURE 14  Color-coded probabilities pHB of hydrogen bonds between the H3 helix (gray, top) or the linker region (light gray, mid), respectively, and adjacent lipid molecules (not shown). The structure is taken from the WT/POP0 annealing set (final structure) and is close to the average orientation. Part of the TM helix is shown in dark gray (bottom).

View larger version (28K): [in this window] [in a new window]   FIGURE 15  Color-coded significances (p-values) of hydrogen bond probability differences of WT/POP- (top) and N10/POP0 (bottom), compared to WT/POP0, denoting significant (p < 0.5) increase (green), confident (p 0) increase (cyan), significant decrease (red), and confident decrease (yellow), respectively, as determined by the bootstrap probability test. The structures are taken from the WT/POP- and the N10/POP0 annealing set, respectively, and are close to the respective average orientations. Atoms not involved in hydrogen bonding are shown in gray (H3 helix), light gray (linker), and dark gray (TM region).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

Room temperature simulations
The linker structure obtained in the pilot study was used as the start structure for the preparation of the explicit bilayer system ( in Fig. 2). Fig. 6 shows the obtained system consisting of three turns of the H3 helix (cyan), the linker (red), and the full TM domain (yellow), embedded in a neutral lipid bilayer (POP⁰, green) composed of POPC and POPE molecules with an explicit solvent environment (water molecules are shown in blue).

Peptide insertion proceeded via an intermediate prosequence of 10 asparagine residues (N₁₀) that was appended at the C-terminus of the TM helix (magenta rectangle). As described in the Methods section, that prosequence was initially appended for the purely technical reason to avoid closure of the hole within the hydrophobic core of the bilayer, which was required to enable peptide insertion. During the subsequent free dynamics, the z-position of the peptide was expected to relax from its biased initial value toward a stable equilibrium. Indeed, a fast relaxation motion was seen during the first 40 ps, where the peptide penetrated the bilayer by further 2 Å (data not shown), and a slower relaxation by 2.5 Å in the same direction during the subsequent 450 ps. Remarkably, and in contrast to the simple membrane model used in the pilot study, during that time, the TM-adjacent linker residue Lys-265 folded into a stable -helical conformation, presumably due to the increased hydrophobicity of the lipid headgroup region (Zhou and Schulten, 1995), which is described more realistically than in the pilot study. As will be analyzed in more detail below, this interaction is indeed quite complex. During the subsequent 1.4 ns, after a temporal backward movement, the peptide finally moved further toward the membrane center by 0.5 Å and reached a stable position, apparently independent of the chosen start position; for the next 3.2 ns, only a small drift is detectable, hence we assume the z-position of the TM helix to be close to equilibrium. The motion of the peptide brought the linker into the polar region of the bilayer. The TM helix, initially perpendicular to the bilayer, stabilized at tilts around 10°.

From the obtained simulation system, the start structures of the four main room temperature production runs were prepared (cf. Fig. 2; Table 3, Set A). For two of these (WT), the prosequence was removed, and 5 ns simulations for the neutral bilayer (POP⁰), composed of POPC and POPE, and a more physiological (Gennis, 1989) acidic bilayer (POP^-), composed of POPC, POPE, and POPS, were carried out (trajectories WT/POP⁰/300K and WT/POP^-/300K, respectively). For the other two simulations, the N₁₀ prosequence was kept.

Linker structure and stiffness
During the two WT simulations, two of the five linker residues remained in random coil configurations providing a hinge region between the H3 and the TM helix. Unexpectedly, already visual inspection of the trajectories suggested that the linker was quite stiff as opposed to, e.g., the highly flexible lipid tails. To quantify the stiffness, Fig. 7A shows the H3 helix tilt during the four simulations. For the neutral bilayer (POP⁰, blue curve), after 600 ps, the tilt angle stabilizes quickly at small values. For the acidic bilayer (POP^-, green curve), the tilt angle starts at higher values and shows a slow drift, until, after 3 ns, it appears to be nearly equilibrated. We have therefore chosen the interval between 3 ns and 5 ns (gray) for the free-energy estimates according to Eq. 1.

View larger version (35K): [in this window] [in a new window]   FIGURE 7  (A) Tilt angle of the H3 helix with respect to the bilayer normal during room-temperature free-dynamics runs of the wild-type (WT) and the prosequence mutant (N10), both embedded within uncharged (POP0) and charged (POP-) bilayer. The period that was used for the free-energy landscape estimate, G(), is shaded gray. (B) Smoothed free-energy landscape for the H3 helix tilt angle, derived from the four free-dynamics runs of 2-ns length each, according to Eq. 1. (C) Linker helicities determined from the four 2-ns free-dynamics runs. Shown are the relative frequencies phelix,i that each of the five linker residues is in helical configuration, given that one of its neighbors is helical.

Effects of prolongation
Unexpectedly, the tilt angle behaved quite differently already in the run, which led us to suspect that the prosequence has a pronounced impact on the mechanical properties of the linker—despite the considerable distance between prosequence and linker region. Therefore we also carried out 5-ns simulations of the mutant both for neutral and for the acidic bilayer (Fig. 7 A, red, magenta). Indeed, even though they start at small tilt angles, both systems show much larger angles and larger fluctuations already after 800 ps compared to the WT. Straightforward calculation of the respective free-energy landscapes is complicated, due to the notorious sampling problem. In particular, the system (red) jumps back to smaller angles after 4.2 ns and remains in that (metastable) state for the remaining simulation, such that, due to insufficient statistics, the computed free-energy landscapes (Fig. 7 B) have to be judged with care. In particular, the relative height of the two apparent minima cannot be estimated unless multiple transitions are observed. In the second part of the Results section, this problem will be approached more systematically by multiple simulated annealing runs.

Table 4 summarizes the results of the—in this case necessarily rough—harmonic fits to the ensembles of N₁₀ tilt angles. As can be seen, for the neutral bilayer, the average tilt angle of the linker of the mutant is about three times larger than that of the WT, whereas the stiffness decreased by a factor of five—mainly due to the fact that two metastable states are encountered. For the acidic bilayer, an even larger (by a factor of four with respect to the WT) average tilt is seen and the stiffness is about three times smaller. Although the individual fit parameters are likely not very accurate at the present stage of analysis, the qualitative difference between the mutant and the WT is evident and thus deserves closer inspection.

Also for the mutant, the level of helicity was studied as possible cause for the changed stiffness (Fig. 7 C, red, magenta). As can be seen, for both bilayer types the conformations of Arg-263, Lys-264, and Lys-265 are similar to the WT systems. In contrast, the quite stable helical conformation of Ala-261 for the WT systems is destabilized for both bilayer types, particularly for N₁₀/POP^- (magenta), which also showed the largest average tilt. Apparently, the helicity of Ala-261 and hence the stability of the backbone hydrogen bond between Ala-261 and Tyr-257 is a critical determinant for the spontaneous H3 helix tilt. The observed correlation is also seen in the time domain, e.g., for N₁₀/POP⁰. Here, the pronounced and abrupt decrease of the tilt angle from 40° to 15° at 4200 ps (Fig. 7 A, red curve) followed a transition of Ala-261 from an equilibrium between helical and nonhelical conformation to a stable helical one at 3070 ps where it remained for the rest of the trajectory. Also, at 3200 ps, Arg-262 folded into a stable helical conformation, which, for the WT, was shown to be crucial for the H3 helix tilt stiffness. Moreover, the observed sequence of events in this case rules out the possibility that helix formation is the consequence of tilt angle changes rather than their cause. Although an (unknown) common cause can not be strictly ruled out, we attribute for both mutant systems the increased average tilt and decreased tilt stiffness to the reduced linker helicity.

Linker embedding
The previous result that inclusion of a realistic membrane environment can change the helical content of the linker suggested a mechanism for the initially unexpected long-distance impact of the prolongation: By changing the z-position of the TM helix, the prolongation could expose the linker region to a more hydrophilic environment and thereby destabilize the backbone hydrogen bond between Ala-261 and Tyr-257, which, in turn, affects the mechanical properties of the linker. To check this hypothesis, for all four systems the vertical linker position z_link with respect to the average positions of the phosphor atoms of the proximate polar region (see Fig. 5) was monitored and the respective free-energy landscapes were determined using Eq. 1. As can be seen in Fig. 8 A, and also from the harmonic fits listed in Table 4 (z_link,0, k_z,link), all four systems exhibit narrow energy minima; accordingly, the linker position is well-defined. As expected, the position of the WT (blue, green) is generally shifted with respect to that of the mutant (red, magenta) in such a way that the linker becomes buried deeper within the membrane, which supports our hypothesis. The effect is particularly strong for the neutral bilayer with a shift of 5.4 Å and less pronounced for the charged bilayer with a shift of 2.7 Å. Interestingly, compared to the neutral bilayer, the lipid charges shift the linker positions of WT and mutant in opposite directions by 1 Å. As shown in Table 4, for N₁₀ the linker is located roughly at the z-position of the phosphate groups. For WT, the linker resides unexpectedly deep in the membrane, suggesting that the membrane environment strongly influences the linker properties in this case.

View larger version (27K): [in this window] [in a new window]   FIGURE 8  Smoothed free-energy landscapes (A) for the linker position zlink with respect to phosphate groups and (B) for the position zTM of the TM helix center with respect to membrane center (dashed line), both derived from the four free-dynamics runs of 2 ns length according to Eq. 1. In A, the polar region is indicated, here defined via the distribution of phosphate group positions, using a cutoff (dashed line) at two standard deviations around the average position (Zhou and Schulten, 1995).

Hydrogen bonds
Generally, the lipid charges are expected to strengthen the interactions between the lipids and the peptide. On the other hand, the prolongation, reducing the peptide penetration, may weaken the peptide-lipid interactions (and thereby increase tilt flexibility). To check these hypotheses, the average number of hydrogen bonds between linker and H3 domain to surrounding lipid molecules, as well as its fluctuations (Eq. 14), were determined for each system. As shown in Fig. 9, as anticipated, the POP^- systems show more frequent hydrogen bonding than the respective POP⁰ systems. This may contribute to the increased tilt stiffness observed for the POP^- systems. For POP^-, a small decrease in the number of hydrogen bonds for N₁₀, as compared to WT, can be observed. For POP⁰, unexpectedly, the average number of hydrogen bonds between peptide and lipids does not strongly differ between WT and N₁₀.

View larger version (13K): [in this window] [in a new window]   FIGURE 9  Average number NHB of hydrogen bonds between the peptide (H3 and linker residues) and lipid headgroups determined from production phases of free-dynamics runs and from ensembles of annealing final structures. The bars denote fluctuations (Eq. 14). Note that the error bars are much smaller than the fluctuations.

Table 4 lists the average positions z_TM,0 and stiffnesses k_z,TM obtained from harmonic fits for the TM positions. As can be seen from the table, as well as from Fig. 8 B, also the mutant, anchored in the neutral bilayer, is slightly shifted out of the symmetric position in the inverse direction by 1 Å. One likely reason is the sequence asymmetry between linker and prosequence and, in particular, the presence of the four positive charges of the linker. That the asymmetry is inverted for the acidic bilayer, which changes the electrostatic environment of the linker, further supports the relevance of the interaction between linker and lipid headgroup region discussed above. The stiffnesses are, overall, very large, corresponding to a force of 300 pN for a 1-Å deflection. Such stiffness is actually indispensable for the proposed membrane fusion mechanism, because it is required to keep the peptide tightly anchored within the bilayer, so that it can withstand the strong repulsion forces that occur when two membranes are pulled toward each other.

TM helix tilt and structure
The second important requirement, as mentioned in the introduction, is that the TM helix must be able to withstand tilting with respect to the membrane normal. To study if this actually is the case, free-energy landscapes also have been computed from the simulations for the TM tilt angle _TM (see Fig. 5). Within the equilibration phase of the simulation that angle departed from its initially chosen value of 0° toward values around 10° (data not shown). The minima of the free-energy landscapes (see Fig. 10) are significantly narrower than those of G(), indicating that _TM is more strongly restrained than . Table 4 lists the corresponding average tilts _TM,0 and the tilt stiffnesses k_,TM obtained from harmonic fits. Overall, the stiffnesses of the TM helix angle are larger than the ones of the linker angle, such that the mechanics should be mainly determined by the stiffness of the linker.

View larger version (19K): [in this window] [in a new window]   FIGURE 10  Free-energy landscapes of TM helix tilt TM with respect to the membrane normal from room temperature production runs.

Simulated annealing runs
Concerning the linker structure and tilt angle, the room-temperature simulations suffer from an apparent sampling problem, as can be seen, e.g., from the fact that only one transition is seen for N₁₀/POP⁰ (Fig. 7 A, red). Also, the physiological membrane fusion timescales are likely three orders of magnitude longer than the nanosecond simulation times (Cevc and Richardsen, 1999). Thus it is unclear as to what extent the linker structures seen in the simulations could be biased by the chosen start structure. Moreover, also the calculated average tilt angles and, in particular, the stiffnesses may be inaccurate due to insufficient sampling.

To reduce this problem and to improve the sampling, we have carried out high-temperature runs followed by multiple simulated annealing runs (Fig. 2, red boxes). For both the WT and the mutant, 1-ns high-temperature runs were carried out at 1500 K as described in the Methods section starting from the respective production run trajectories (black arrows). For the WT in the neutral bilayer, three different cooling rates were used to assess convergence. As also described in the Methods section, the TM helix and the lipid tails had to be restrained for the high-temperature and annealing runs to avoid degradation of the system. For the successive free-dynamics room-temperature runs, the restraints were removed. Inasmuch as z_TM and _TM are much better sampled and show only small fluctuations as compared to the H3 helix tilt angle (cf. Figs. 8 A and 10), we do not expect these restraints to severely limit the sampling of the linker structure and the linker tilt angle. In particular, results from recent measurements on the membrane immersion of the syntaxin linker are consistent with our results of the average linker position and in this regard support our assumptions (Kweon et al., 2002).

Fig. 11 A shows the results of a set of WT/POP⁰ annealing runs (cf. Fig. 4) compared with the respective room-temperature production trajectory (blue). Rather than the tilt angle , here the projection of the (normalized) H3 helix axis onto the xy bilayer plane has been plotted; the center corresponds to = 0°. As can be seen, the 1-ns 1500-K sampling run (brown dots) provided a much broader distribution of H3 helix orientations than the 2-ns room-temperature production phase. Furthermore, the region covered by the room-temperature run is only rarely visited within the sampling run, suggesting that the free energy of other regions is lower, at least at this high temperature. This finding can be explained by the fact that the distributions of H3 helix orientations are governed by respective free energy landscapes, G = H - TS, which, in addition to the enthalpic contribution H, also include temperature-dependent entropic contributions, -TS, which generally tend to widen the distribution for increased temperatures.