Structural Stability of Wild Type and Mutated α-Keratin Fragments

Mutated r-Keratin. Fragments: Molecular Dynamics and Free Energy Calculations ... heptad periodicity within the rod domain is interrupted at sever...
0 downloads 0 Views 355KB Size
Biomacromolecules 2004, 5, 2165-2175

2165

Structural Stability of Wild Type and Mutated r-Keratin Fragments: Molecular Dynamics and Free Energy Calculations Cristian Danciulescu, Birgitta Nick,* and Franz-Josef Wortmann German Wool Research Institute, Veltmanplatz 8, 52062 Aachen, Germany Received April 8, 2004; Revised Manuscript Received June 24, 2004

The objective of this study is to investigate the influence of point mutations on the structural stability of coiled coil fragments of the human hair intermediate filament by molecular dynamics simulations and free energy calculations. Mutations in the helix termination motif of human hair keratin gene hHb6 seem to be connected to the hereditary hair dystrophy Monilethrix. The most common mutations reported are Glu413Lys and Glu413Asp, located at the C-terminal end of the coiled coil 2B rod domain of the IF. According to our simulations, significant conformational changes of the side chains at the mutation and neighboring sites occur due to the Glu413Lys mutation. Furthermore, the differences in electrostatic interactions cause a large change in free energy during transformation of Glu413 to Lys calculated by the thermodynamic integration approach. It is speculated that the structural rearrangement necessary to adapt the interactions in the mutated coiled coil leads to changes in the IF assembly or its stability. The second mutation, Glu413Asp, only leads to a small value of the calculated free energy difference that is within the error limits of the simulations. Thus, it has to be concluded that this mutation does not affect the coiled coil stability. Introduction R-Keratin fibers such as hair and wool show a complex morphological structure.1 They exhibit a high content of fibrous proteins, grouped together in filaments and fibrils oriented parallel to the fiber axis. The main structural components of the keratin fiber are the cuticle and the cortex. The cuticle is a surface layer that covers the core of the fiber. It consists of layers of cells that overlap each other, resulting in a scale structure, with the cuticle edges pointing toward the tip of the fiber. The cortex consists of elongated cells of approximately 100 µm length. They have an irregular crosssection of a few micrometers, separated from one another by the cell membrane complex. Cortical cells consist of long filaments oriented parallel to the axis of the fiber, traditionally called microfibrils. These microfibrils are in current semantics referred to as intermediate filaments (IF) due to their compositional and structural similarity with the IFs of the cytoskeleton of the animal cells. In the cortical cell, the intermediate filaments are embedded in a sulfur-rich matrix and grouped together in units of about 0.5 µm in diameter, called macrofibrils. The intermediate filament structure can be described through several hierarchical levels of increasing order, from the primary structure to the quaternary structure. The models of the IF molecular architecture have a lower level of confidence for the higher levels of structural organization. Thus, the tertiary and the quaternary structures of the intermediate filament remain poorly described in detail.2 The IF monomeric units have a three-domain structure that consists of a central R-helical “rod” domain, flanked by * To whom correspondence should be addressed. Phone: +49 241 4469 118. Fax: +49 241 4469 100. E-mail: [email protected].

Figure 1. Schematic representation of the intermediate filament dimer (adapted from ref 2).

“head” and “tail” domains. The rod domain reveals a heptad repeat pattern, which indicates an R-helical structure. The heptad periodicity within the rod domain is interrupted at several positions (Figure 1), resulting in four consecutive R-helical segments differing in length, namely 1A, 1B, 2A, and 2B. The helical segments are connected by short nonhelical linkers L1, L12 and L2. The initial step in IF assembly is the formation of a dimer (Figure 1) by the parallel alignment of two compatible IF monomeric units, type I (acidic) and type II (neutral to basic), resulting in a left-handed coiled coil structure that is 45-47 nm long. The next level of aggregation is the protofilament (tetramer) that is formed by the antiparallel alignment of two coiled coil dimers in one of the basic modes named A11, A22, and A12, based on which portions of the rod domains overlap,2 i.e., A11, an antiparallel arrangement of dimers with their 1B segments largely overlapping, A22, an antiparallel arrangement of dimers with their 2B segments largely overlapping, and A12, an antiparallel arrangement of dimers with their 1B and 2B segments largely overlapping. It is likely that two protofilaments wind around a common axis in a right-handed manner and form the protofibril. Finally, four protofibrils (32 chains in cross-section) would likely coil around one another in a left-handed manner to form the intermediate filament.3 Recently, several studies of IFs have shown that mutations in the genes encoding the various chains often result in

10.1021/bm049788u CCC: $27.50 © 2004 American Chemical Society Published on Web 09/10/2004

2166

Biomacromolecules, Vol. 5, No. 6, 2004

Danciulescu et al.

Figure 3. Schematic illustration of the parallel arrangement of two R-helical chains in a coiled coil structure (axial projection).

Figure 2. Polarized light micrograph of hair from a monilethrix patient (adapted from ref 4).

diseases.4-7 Most of the mutations that have been observed involve single nucleotide changes, resulting in amino acid substitutions in the keratin chain. Frequently observed mutations in the helix termination motif of two type II (basic) human hair keratin genes, hHb1 and hHb6, encoding the R-helical 2B domain seem to be linked to an hereditary hair dystrophy, known as Monilethrix.4-7 Clinically, monilethrix is characterized by a dystrophic alopecia of varying size and rough follicular papules. Hairs from affected individuals show a beaded structure of alternating elliptical nodes and constrictions (Figure 2). These internode regions exhibit a high propensity to weathering and fracture. Ultrastructural hair studies of monilethrix patients suggest local edematous degeneration in matrix cells of the hair bulb, together with cytolysis and disrupted keratin filament packing in the cortical cells.6 Monilethrix patients exhibit a point mutation in the hHb6 gene encoding a type II hair keratin that leads to the substitution of a highly conserved glutamic acid with lysine (Glu413Lys). Another point mutation in the same glutamic acid codon of hHb6, which results in an aspartic acid substitution (Glu413Asp) was also reported.4 Clinical information from monilethrix families suggests that the Glu413Lys mutation is generally associated with a more severe phenotype of the disease, whereas the Glu413Asp mutation is typical for a milder phenotype.4,5 Other studies with synthetic peptides also indicate that Glu f Asp substitution results in only mild alteration of the filament structural stability.8 More general studies have indicated that many proteins that form two-chain coiled coils possess a 13-residue “trigger” motif that confers special stability for the formation of a coiled coil.9,10 The motif as currently defined is absent in the IF chains, but regions of some similarity occur toward the ends of the 1A, 1B, and 2B rod domain segments. Wu et al. showed that there is a region in each of the 1B and 2B segments that is essential for the stability of coiled coils and thus may serve as the trigger motif.11 The highly conserved region located at the end of the 2B segment is called IF “consensus” motif and consists of the following residues: YRKLLEGEE.12 Based on the atomic structure of this consensus motif, there are structural characteristics that are supposed to contribute significantly

to the IF assembly. In vitro data clearly prove that this amino acid sequence directs the proper formation of tetramers, and controls the number of subunits per filament cross-section.12 Some of the residues from this consensus motif may have a dual role in the intermediate filament structure: (a) they may participate in the trigger formation to confer special stability to the two-chain coiled coil and (b) they may participate in stabilization of the two-dimer assembly. The role of each amino acid from the end of the 2B domain in the IF assembly (dimer or tetramer) is not yet completely understood. The YRKLLEGEE motif may vary considerably for various types of IF. For example, phakinin and filensin from vertebrates that form IF-type filaments exhibit significantly different sequences: YHALLDREE and YHRIIENEG, respectively.13,14 In contrast, in cytokeratins, comparatively minor amino acid replacements in the consensus motif appear to have a drastic effect on IF structure. Thus, mutations occurring in the consensus motif of keratin intermediate filaments may cause severe diseases. The glutamic acid (E) at the 8th site of the YRKLLEGEE sequence occupies a d position in the heptad repeat of the R-helix. So this residue is directed into the coiled coil interface (Figure 3). The heptad positions a and d are occupied usually by nonpolar residues that build up a hydrophobic interface between the helical strands. The role of this glutamic acid in the stability of IF structure has not been investigated yet. It is the objective of this study to investigate the influence of the specific mutations involved in the monilethrix disease (Glu413Lys and Glu413Asp located on hHb6) on the structural stability of affected monomeric and dimeric segments from the end of the 2B domain. For this purpose, molecular dynamics (MD) simulations of the selected protein segments are performed in order to obtain information about their stability and the features of their dynamical behavior, such as the flexibility of the chains and the mobility of the amino acid side chains in the regions of interest. Special attention is paid to the dynamics of Glu413 and Lys413 side chains and their electrostatic interactions with neighboring charged residues. The relative thermodynamic stability of the IF fragments is estimated by calculating relative free energy differences between the wild type and the mutated forms (Glu413Lys and Glu413Asp). Previous studies have proven that molecular dynamics simulation is a powerful tool for investigating the stability of the secondary structure of the R-keratin intermediate

Wild Type and Mutated R-Keratin Fragments

Biomacromolecules, Vol. 5, No. 6, 2004 2167

Table 1. Amino Acid Sequences of the Investigated IF Fragmentsa amino acid sequence segment name

a

b

c

d

e

f

g

K(+) 396 I 403 L 410 L 417

Q 390 L 397 A 404 E(-) 411 C 418

E(-) 391 G 398 T 405 G 412 E(-) 419

V 392 L 399 Y 406 E(-) 413 G 420

M 393 D(-) 400 R(+) 407 E(-) 414 V 421

N 394 I 401 R(+) 408 Q 415 G 422

S 395 E(-) 402 L 409 R(+) 416

Y 389 K(+) 396 I 403 L 410 L 417

Q 390 L 397 A 404 E(-) 411

E(-) 391 G 398 T 405 G 412

V 392 L 399 Y 406 E(-) 413

M 393 D(-) 400 R(+) 407 E(-) 414

N 394 I 401 R(+) 408 Q 415

E(-) 388 S 395 E(-) 402 L 409 R(+) 416

Y 339 R(+) 346 I 353 L 360 L 367

Q 340 A 347 N 354 E(-) 361

V 341 R(+) 348 T 355 S 362

L 342 L 349 Y 356 E(-) 363

L 343 E(-) 350 R(+) 357 D(-) 364

D(-) 344 C 351 S 358 C 365

E(-) 338 V 345 E(-) 352 L 359 K(+) 366

single helix (hHb6)

coiled coil

chain A (hHb6)

chain B (hHa4)

a The residues are arranged into seven columns, indicating the heptad positions a to g. The columns a and d that form the interface between the two strands of the coiled coil are in boldface. The last five residues of single helix model belong to the tail domain of the IF monomer. The charge assigned to each ionizable amino acid in the MD simulations is displayed in brackets.

filament of wool. Choosing an ideal R-helical conformation as the start conformation for each IF segment, Knopp et al. carried out MD simulations in vacuum for 500 ps at systematically increasing temperatures (ranging from 280 to 440 K). The stability of the helical segments was estimated from their ability to stay folded as the temperature is raised.15-19 Similar approaches were applied by Daggett and Levitt in their study of protein denaturation20,21 and by Hirst and Brooks in the investigation of the stabilities of helical fragments of myoglobin.22 The study of the kinetic stability of IF fragments by molecular dynamics methods was further extended by Nick et al.23,24 and Kastanja et al.25 to mutated forms of IF fragments. These studies point out the potential of the molecular dynamics methods for investigating the relative stability of secondary and tertiary structures of keratin IF fragments. Materials and Methods Selected IF Fragments. Two R-keratin IF fragments from the end of the helical 2B domain, namely a single helix of 33 residues and a coiled coil segment of 2 × 30 residues, were selected. The single helix model is a segment of the hHb6 IF monomeric unit that corresponds to residues 390 to 422 of the primary structure.26,27 The chain A of the coiled coil is also a hHb6 segment and consists of the amino acid residues 388 to 417.26,27 The chain B of the coiled coil is a segment of the hHa4 IF monomeric unit, and consists of the amino acid residues 338-367.28 The hHa4 is supposed to be a plausible partner of hHb6 in a coiled coil dimer.29 The secondary and tertiary structures of the IF segments were assumed in conformity with the description of ref 29. The amino acid sequences of the wild-type segments are listed in Table 1 in one-letter code. The amino acid sequence of single helix segment (33 residues) differs slightly from that of the chain A (30 residues) of the coiled coil, namely, the single helix model

contains the last 28 residues of the 2B domain and the first 5 residues from the tail domain of the IF monomer. As already mentioned, position 413 (involved in the mutations) is situated four residues from the end of the helical segment 2B. Due to the fact that the ends of a single helix in solution have the tendency of fraying,15,20 five more residues were added to the helical segment in order to locate the mutation in a more stable region of the helix. This allows a proper investigation of the influence of the mutation on the dynamical behavior and the stability of the helical segment. The single helix segment was generated by building up the amino acid sequence using Swiss PDB-Viewer30 and imposing an ideal right-handed R-helical conformation onto the peptide chain. The three-dimensional structure of the coiled coil segment was obtained from the GCN4 leucine zipper31 that also exhibits a left-handed coiled coil of two right-handed R-helices of similar twist parameters. The amino acid side chains of GCN4 were replaced with the side chains corresponding to the IF monomeric units hHb6 and hHa4. Afterward the structure was refined by several energy minimization steps toward a suitable starting conformation.32 In addition to the wild type segments, mutated forms were obtained by replacing the glutamic acid at the mutation site 413 with both experimentally reported mutations, i.e., Lysine and Aspartic acid. The single helices and the coiled coil segments resulted from substituting Glu with Lys and Asp, as well as the wild-type segments, were used to perform MD simulations and free energy calculations, using the Gromacs 3.1.4 software package.33 MD Simulation Setup. The MD simulations were carried out using periodic boundary conditions. To this, rectangular simulation boxes were defined with the dimensions of 4.5 × 4.2 × 4.8 nm3 for the single helix and 5.0 × 5.0 × 5.5 nm3 for the coiled coil. Each protein fragment was placed in the center of the box; subsequently the boxes were filled with water. To do this, the protein segment was overlaid by equilibrated simple point charge (SPC)34 water boxes. All

2168

Biomacromolecules, Vol. 5, No. 6, 2004

Danciulescu et al.

Figure 4. Thermodynamic cycle with physical (horizontal) and nonphysical (vertical) paths for determining the relative free energy of association, ∆∆G, from experiments and from computer simulations.

water molecules with the oxygen atom lying within 0.23 nm of a non-hydrogen protein atom were removed. In this way, the single helix system contains about 2940 water molecules and the coiled coil system contains about 4230 water molecules. A steepest-descent energy minimization was performed in order to remove local strains in the peptide and to relax the shells of water molecules around it. The interactions of the molecular systems were described by the GROMOS96 43A1 force field.35 This force field treats the aliphatic hydrogen atoms as united atoms.35 The ionizable amino acids (Glu, Asp, Lys, Arg) were considered in their charged state. Thus, the single helix carries a total charge of -3 and the coiled coil segment carries a total charge of -6 in the wild type form (Glu413). In the mutated form (Lys413), the single helix carries a total charge of -1 and the coiled coil carries a total charge of -4. The nonbonded interactions are evaluated using two different methods: (a) The twin-range cutoff method (with a short-range cutoff of 1.0 nm and a long-range cutoff of 1.4 nm): the van der Waals and electrostatic interactions within the short cutoff distance (1.0 nm) are evaluated at every time step by taking into account the pairs of atoms from a neighbor (pair) list generated with a cutoff radius of 1.0 nm at every 10 time steps. A long cutoff distance of 1.4 nm is set only for the electrostatic interactions. The electrostatic interactions between atoms separated by a distance ranging between 1.0 and 1.4 nm are updated every 10 time steps, being kept unchanged between the updates of the pair list (10 steps). (b) Electrostatic interactions are evaluated using the Particle-Mesh Ewald (PME) algorithm.36 This method provides a better treatment of the electrostatic interactions since every atom interacts with all the other atoms in the simulation box and with all of their images in an infinite array of periodic cells. The MD parameters for PME method were set to the default values, namely: cutoff distance for the neighbor list ) 0.9 nm, van der Waals cutoff ) 0.9 nm, Coulomb type ) PME, Coulomb cutoff ) 0.9 nm. The fast Fourier transform grid dimensions are controlled with a grid spacing of 0.12 nm and an interpolation order of 4 (corresponding to cubic interpolation). The molecular dynamics simulations of the IF systems were performed using the leapfrog integration algorithm with a time step of 2 fs. This time step was possible due to the use of bond constraints. The bonds between the atoms of the protein were constrained to ideal values using the LINCS algorithm.37 The geometry (bonds and angles) of the water

molecules was constrained using the SETTLE algorithm.38 The initial velocities of the atoms were taken from a Maxwell-Boltzmann distribution at the temperature of interest. The temperature and the pressure were maintained to the desired values by coupling to temperature and pressure baths.39 The temperature of the peptide and the solvent were independently coupled to a Berendsen thermostat using a relaxation time of 0.1 ps. The pressure of the system was weakly coupled to the pressure bath with isotropic scaling and a relaxation time of 0.5 ps. For each system, multiple MD simulations, usually three, were performed by assigning different initial velocities to the starting conformation. Free Energy Simulations. The molecular dynamics free energy simulation method is well suited to explore the energetic consequences of amino acid substitutions on protein structure using detailed atomic models. The method is based on the fact that the free energy is a state function. Thus, the sum of the free energy changes for a series of reversible transitions in a closed path is zero. Free energy calculations in conjunction with molecular dynamics have been used to study protein-ligand or protein-protein interactions40-43 as well as protein stability.44-46 In this study, free energy calculations are used to investigate the thermodynamic stability of wild type and mutated coiled coils. Employing a standard thermodynamic cycle to investigate the protein stability, free energy calculations have to be performed for both the folded and denatured states. However, inherent difficulties arise when attempting to assume a structure for the denatured state. It has already been shown that the calculated free energy difference ∆∆G depends significantly on the structure of the denatured state.47 Especially when the mutation involves charged residues, it is unlikely that an extended peptide model45,48,49 is appropriate for the denatured state. Consequently, a different approach is followed here. Instead of calculating the free energy difference due to the mutation in folded and denatured conformation, the relative free energy of association of R-helical segments to form wild type and mutated coiled coils is calculated. This approach involves the use of well-defined structures, i.e., single R-helices and coiled coils. The thermodynamics of the association of helical segments can be investigated by employing the thermodynamic cycle shown in Figure 4. ∆Gassoc.1 and ∆Gassoc.2 represent the free energies of associations of the wild-type helical segments and the mutated segments, respectively. ∆G1 and ∆G2 are the free energy changes due to mutating Glu f Lys or Glu f Asp

Wild Type and Mutated R-Keratin Fragments

Biomacromolecules, Vol. 5, No. 6, 2004 2169

in an isolated helix and in the coiled coil segment, respectively. The free energies of association, ∆Gassoc.1 and ∆Gassoc.2, of the wild type and the mutated proteins can in principle be determined experimentally; ∆G1 and ∆G2 are accessible only by computer simulation. The simulation results can be related to the experimental data via the thermodynamic cycle. Thus, it can be written ∆∆G ) ∆Gassoc.2 - ∆Gassoc.1 ) ∆G2 - ∆G1

(2)

To calculate the free energy differences associated to the interconversion of residue Glu413 to Lys and to Asp the thermodynamic integration approach is used.35 Thus, the Hamiltonian is made a function of a coupling parameter λ that varies progressively from 0 to 1, to change the system interactions from the initial state (A) to the final state (B). The difference of free energy between the molecular systems A and B is then calculated using ∆G ) G(λB) - G(λA) )

∫01〈



∂H(λ) dλ ∂λ λ

Figure 5. Schematic representation of the hybrid residue E2K used for simulating the mutation Glu to Lys. DN, DH, and DO (circled) represent dummy atoms.

(3)

The simulations are done at a few discrete points λi along the pathway and the integral is calculated numerically. The transformation of glutamic acid into lysine or aspartic acid is associated with a change in the number of interacting atoms. Dummy, noninteracting atoms must be added to the molecule to model the amino acid transformation at a constant number of atoms during the simulation. A dummy atom is an atom for which the nonbonded interactions, namely the Lennard-Jones and electrostatic interaction parameters with all other atoms have been set to zero. Thus, a normal atom can be “created” in the system by gradually mutating a dummy atom into the requested normal atom. To remove a normal atom from the system, the atom should be gradually mutated into a dummy atom. The introduction (or removal) of an atom into a system causes a rapid variation of the potential energy at the beginning (or end) of the MD simulation50 because of the occurrence of close contacts (r = 0) between two atoms. The singularities in the potentials are removed by using soft-core potentials.51,52 To perform the mutation Glu f Lys by free energy simulations, the hybrid residue Glu-Lys (referred to as residue E2K) was designed and included in the GROMACS topology file. The hybrid was built up by connecting a group of 4 dummy atoms to a carboxyl oxygen of the glutamic acid residue (Figure 5). The dummy atoms DN(DH)3 have the spatial arrangement of an amino group NH3+. During the transformation Glu f Lys, one of the carboxyl oxygens becomes the dummy atom DO at the end state and the other carboxyl oxygen is mapped onto the Lys-C (the united atom CH2). The C atom of the COO- group is mapped onto the Lys-Cδ atom (united atom CH2). The group DN(DH)3 becomes the ammonium group (NH3+) of Lys at the end of the transformation. Figure 6 shows the equivalent hybrid residue E2D for the mutation Glu f Asp. A hybrid topology is designed in the manner described before, i.e., by introducing dummy atoms where necessary and mapping the atoms of Glu onto those of Asp.

Figure 6. Schematic representation of the hybrid residue E2D used for simulating the mutation Glu to Asp. DO represents a dummy atom.

Because the same bond lengths and force constants are used for dummy atoms as for their corresponding normal atoms, no free energy contribution arises from these terms. The dihedral and the improper dihedral angles are treated differently, namely the force constants of the harmonic potentials that describe these angles are set to zero for the dummy atoms. However, the free energy contribution arising from the change in the dihedral angles is small and it cancels in the thermodynamic cycle. To model the free energy change due to the mutation, the glutamic acid located at the mutation site (413) is replaced with the hybrid residue, either E2K or E2D, that goes through the alchemical transformation. The backbone atoms, N, CR, C, are positionally restrained during the free energy simulations, to preserve the starting conformation during the entire simulation. No restraints are applied to the side chains of the amino acid residues. Additionally, MD simulations without position restraints are also performed, to verify the influence of the restraints on the free energy calculations. All MD simulations are performed over 200 ps at a temperature of 300 K. The atom-atom interactions are described using the GROMOS96 force field.35 The nonbonded interactions are evaluated using the twin-range cutoff method. To check the influence of the cutoffs on the calculations, two different sets of cutoffs are used, namely 1.0/1.4 nm and 1.4/1.8 nm. The singularities in the potentials due to creation/annihilation of atoms in the system are

2170

Biomacromolecules, Vol. 5, No. 6, 2004

Danciulescu et al.

removed by using soft-core potentials. The soft-core parameter R is assigned a value of 1.5 and the radius of interaction σ a predefined value of 0.3 nm in cases where C6 or C12 is zero. For each system (single helix and coiled coil) separate simulations are performed at a number of λ points (minimum 12) from λ ) 0 to 1. At each λ point, the system was first equilibrated for 100 ps and then data were collected for the further 100 ps. The averages of the derivatives, dG/dλ, at each λ were integrated using the trapezoidal method to obtain the free energy difference, ∆G. Results and Discussion MD Simulations with Coiled Coil Structures. Molecular dynamics simulations on the coiled coil structures are carried out at 300 K using the twin-range cutoff method (1.0/1.4 nm) and PME method for calculating the nonbonded interactions. The end conformations of the coiled coil segments after 1 ns MD simulations using the twin-range cutoff method are shown in Figure 7. No significant difference in stability can be quantified from these snapshots. To this, the RMS deviation of the backbone atoms is calculated from each trajectory. For both runs the RMSD values are below 0.2 nm during the first 300-400 ps. Afterward, they slowly increase with the simulation time, for both wild type and mutated segments. Thus, the average RMSD value corresponding to the end conformations of the wild-type coiled coil obtained from different MD trajectories at 300 K is 0.32 ( 0.07 nm. The corresponding RMSD value for the mutated coiled coil is 0.35 ( 0.15 nm. The analysis of the RMS fluctuations of the CR atoms puts into evidence a higher flexibility of the terminal residues of the helices. By visual inspection of the MD trajectories, one can see bending of the helices toward the chain ends or fraying of the marginal residues. Since the system under investigation is highly charged equivalent MD simulations were carried out using the PME method. This method is expected to give an more accurate description of the electrostatic interactions of the system. Both, the wild-type and the mutated coiled coil segments exhibit a remarkable stability of conformation at 300 K during an 1 ns MD simulation (Figure 7). The RMS deviations of the backbone atoms from the starting structure are low during the whole simulation, fluctuating around the average value of 0.1 nm. These low RMSD values indicate a stable structure over the whole simulation time. Again, no difference in the structural stability of the wild type and mutated coiled coils can be detected from the RMSD analysis. Electrostatic Interactions at the Coiled Coil Interface. To understand the underlying effects of the interactions stabilizing the coiled coil structure, the interactions at the coiled coil interface are investigated in more detail. The electrostatic interactions between the residues of chain A and the residues of chain B of the coiled coil were analyzed in detail using MD trajectories at 300 K obtained with the PME method. The pairs of interacting charged residues were determined by plotting the distance between all pairs of

Figure 7. End conformations of wild type (a, a′) and mutated (b, b′) coiled coil from 1 ns MD simulations at 300 K using twin-range cutoff (a, b) and PME (a′, b′) method. The side chains of Glu413 (wild type) and Lys413 (mutated) are schematically represented.

charged residues (positive and negative) as a function of time. To limit the number of interacting pairs and consider only those that have the main energetic contribution, a cutoff distance of 1.0 nm was set and the groups that were never closer than this distance were not inspected further. The selected pairs of charged amino acids, the position of each residue in the heptad repeat and the relative position of the interacting residues in the coiled chains, as well as their interacting distances are presented in Table 2 (wildtype coiled coil) and Table 3 (mutated coiled coil). In these tables, the heptad positions as well as the relative residue positions of chain A and B are denoted without and with the ′ character, respectively. The analysis of the electrostatic interactions in the coiled coil segments is focused on the qualitative changes in the potential energy induced by the mutation Glu(-) f Lys(+). Examining Table 2 one observes that Glu413 is not involved in attractive interactions with residues from chain B. On the contrary, Glu413 has repulsive interactions with

Wild Type and Mutated R-Keratin Fragments

Biomacromolecules, Vol. 5, No. 6, 2004 2171

Table 2. Electrostatic Interactions in the Wild-Type Coiled Coil Segmenta pairs of interacting amino acids (chain A-chain B) negative-negative Asp400-Glu352 Glu413-Glu363 Glu413-Asp364 Glu414-Glu363 positive-positive Lys396-Arg346 Lys396-Arg348 Arg407-Arg348 negative-positive Glu388-Arg346 Glu391-Arg346 Glu402-Arg357 Arg407-Glu352 Glu414-Lys366 Arg416-Asp364

position in the coiled coil (chain A-chain B)

average distance between charged groups r (nm) ( s.d.

e-g′ (n - n′+2) d-d′ (n - n′) d-e′ (n - n′+1) e-d′ (n - n′-1)

0.95 ( 0.12 0.75 ( 0.15 0.60 ( 0.11 0.64 ( 0.13

a-a′ (n - n′) a-c′ (n - n′+2) e - c′ (n - n′-9)

0.88 ( 0.07 0.92 ( 0.12 0.99 ( 0.23

g-a′ (n - n′+8) c-a′ (n - n′+5) g-e′ (n - n′+5) e-g′ (n - n′-5) e-g′ (n - n′+2) g-e′ (n - n′-2)

1.00 ( 0.24 0.76 ( 0.18 0.54 ( 0.21 0.52 ( 0.14 0.89 ( 0.18 0.89 ( 0.24

a The heptad positions of the interacting residues (a through g) and their relative positions in the chains (n - n′( x) are displayed in the second column. The character ′ indicates chain B. The interactions of amino acid occupying the mutation site with other amino acids are in bold type.

Table 3. Electrostatic Interactions in the Mutated Coiled Coil Segment (Lys413)a pairs of interacting amino acids (chain A-chain B) negative-negative Asp400-Glu352 Glu414-Glu363 positive-positive Lys396-Arg346 Lys396-Arg348 Arg407-Arg348 Lys413-Arg357 negative-positive Glu391-Arg346 Glu402-Arg357 Arg407-Glu352 Lys413-Glu361 Lys413-Glu363 Lys413-Asp364 Glu414-Lys366 Arg416-Asp364

position in the coiled coil (chain A-chain B)

average distance between charged groups r (nm) ( s.d.

e-g′ (n - n′+2) e-d′ (n - n′-1)

1.00 ( 0.14 0.63 ( 0.17

a-a′ (n - n′) a-c′ (n - n′+2) e-c′ (n - n′-9) d-e′ (n - n′-6)

0.89 ( 0.07 1.02 ( 0.15 0.97 ( 0.19 1.07 ( 0.13

c-a′ (n - n′+5) g-e′ (n - n′+5) e-g′ (n - n′-5) d-b′ (n - n′-2) d-d′ (n - n′) d-e′ (n - n′+2) e-g′ (n - n′+2) g-e′ (n - n′-2)

0.59 ( 0.17 0.43 ( 0.13 0.49 ( 0.15 0.77 ( 0.17 0.76 ( 0.08 0.29 ( 0.07 0.88 ( 0.27 0.83 ( 0.25

a The heptad positions of the interacting residues (a through g) and their relative positions in the chains (n - n′( x) are displayed in the second column. The character ′ indicates chain B. The interactions of amino acid occupying the mutation site with other amino acids are in bold type.

Asp364 (average interaction distance 0.6 nm) and Glu363 (average interaction distance 0.75 nm). These interaction distances around 0.6 to 0.75 nm are relatively short distances for electrostatic interactions which act on much longer ranges. Consequently, the repulsive interactions due to the presence of Glu413 are expected to locally (at the C-terminal end) induce a destabilization of the coiled coil structure. An analysis of the hydrogen bonds showed that Glu413 has no hydrogen bonds with polar groups from the chain B of the coiled coil. The side chain of Glu413 is partially buried in the interface of the coiled coil, but its carboxyl group is

Figure 8. Distances Glu413-Asp364 (wild-type coiled coil) and Lys413-Asp364 (mutated coiled coil), respectively, as a function of simulation time. The distance is calculated as the minimum distance between the residues, i.e., the distance between the closest two atoms at any given time.

exposed to the solvent, participating in hydrogen bonds with the surrounding water molecules. The substitution of Glu413 with Lys significantly changes the electrostatic interactions concerning the ratio between the attractive and the repulsive interactions at the C-terminus of the coiled coil (Table 3). Lys413 participates in attractive electrostatic interactions with the adjacent charged residues of chain B. The main interaction is Lys413-Asp364 which has an average distance of 0.29 nm. This indicates that Lys413 and Asp364 form a salt bridge during a large part of the simulation time. Examining the distances between Glu/Lys413 and Asp364 (d-e′) as a function of time as presented in Figure 8, a significant difference in the side chains distance is observed. The side chains of Lys and Asp point toward each other for most of the time and form a salt bridge, while the side chains of Glu413 and Asp364 point in opposite directions due to the repulsive interaction. There are two other attractive interactions in the mutated coiled coil segment, Lys413-Glu361 and Lys413-Glu363, that locally influence the electrostatic energy, most likely by stabilizing the structure. The other side chain - side chain interactions in the mutated coiled coil exhibit similar average interacting distances as in the wild-type coiled coil (see Tables 2 and 3). Thermodynamics of the Interactions in the Coiled Coil Segments Studied by Free Energy Calculations. Prior to the free energy calculations, MD simulation with the single helix model were carried out. These simulations show the R-helical structures to be stable over the first 200 ps of the MD simulation time at 300 K which is an essential prerequisite for the free energy calculations. Furthermore, any effects introduced by replacing the normal Glu413 residue with the designed hybrid residues E2K and E2D are checked during MD simulation. The data show that the hybrid residues are able to reproduce the dynamics and the spatial arrangement/orientation of the residues Glu413, Lys413, and Asp413 in the coiled coil structure. Free Energy Calculations on Glu413Lys. The estimation of the free energy difference due to the mutation Glu413Lys

2172

Biomacromolecules, Vol. 5, No. 6, 2004

Danciulescu et al.

Figure 9. Ensemble average 〈dG/dλ〉 as a function of λ, obtained during the transformation Glu(-) f Lys(+) in the single helix and coiled coil (1st run, cutoffs 1.0/1.4 nm). The error bars corresponding to the standard deviation of dG/dλ at each λ value are shown for the change of λ in the single helix. The standard deviations of dG/dλ are in the same range for the coiled coil; for reasons of clarity the error bars were not plotted on this graph.

Figure 10. Potential energy of the interactions of residue E2K to single helix and to water as a function of the coupling parameter λ.

Table 4. Calculated Free Energy Differences Corresponding to Mutating Glu(-) f Lys(+) in the Single Helix and in the Coiled Coila

method 1.0/1.4 nm (1st run) 1.0/1.4 nm (2nd run) 1.0/1.4 nm, no restraints 1.4/1.8 nm 1.4/1.8 nm, no restraints

∆ G1 (single helix) (kJ mol-1) -116 -109 -106 -110 ( 5

∆G2 (coiled coil) (kJ mol-1)

∆∆G (calc.) (kJ mol-1)

-169 -186 -140 -174 -152 -164 ( 18

-53 -77 -28 -68 -46 -54 ( 19

a

The double difference ∆∆G represents the relative free energy of association of wild-type and mutated helices (chain A and chain B). In the case of the unrestrained systems ∆∆G is calculated using the mean of the restraint single helix ∆G1 values, i.e., -112 and -106, respectively.

is based on MD simulations at various intermediate states of the coupling parameter λ. The average of the derivative of the potential energy at the various λ-dependent states with respect to λ as well as the error bars corresponding to the fluctuations of dG/dλ at each λ value are plotted in Figure 9 for the single helix and the coiled coil system. It can be seen in Figure 9 that the derivatives dG/dλ are very large. This may lead to large statistical errors in the estimation of ∆G. However, the initial equilibration of the system, followed by collecting of data over 100 ps keeps the statistical error within reasonable limits. The nearly linear shape of the curves dG/dλ as a function of λ clearly indicates that the free energy change is determined by electrostatic interactions. This is in agreement with expectations since the charge of the side chain changes from -1 to +1 when residue E2K changes from glutamic acid to lysin. The final results of the free energy calculations obtained after integration of the derivative of ∆G with respect to λ for the mutation Glu413Lys in the single helix segment and in the coiled coil are shown in Table 4 for various simulation conditions applied. Considering the large values in Figure 9 the values for ∆G1 and ∆G2 are reproducible and quite independent of the cutoff. Removing the restraints during

Figure 11. Potential energy of the interactions of residue E2K with chain A, chain B, and water as a function of the coupling parameter λ.

the coiled coil simulation does increase the ∆G2 value up to 20%, but this does not change the final result. Unrestrained simulations of the single helix were not performed due to the high flexibility of the helical chain ends. The fraying of the marginal residues would significantly influence the interactions at the mutation site (located near the C-terminus) and would probably lead to erroneous free energy differences. Nevertheless, the limitation in performing unrestrained simulations of the single helix increase the uncertainty concerning the influence of the mutation in the isolated R-helix. However, a significant deviation of the calculated ∆G1 from the value shown in Table 4 is not expected. As one can see, mutating Glu to Lys results in all cases in a negative free energy difference for both systems (single helix and coiled coil), with larger negative values corresponding to the coiled coil system. The calculated ∆∆G of -54 ( 19 kJ mol-1 indicates a more favorable association of the helical segments involving the mutated form (Lys413). To understand the contributions to the free energy changes, the energies involved in the interactions of residue E2K with the rest of the peptide and with the water are analyzed for the single helix and coiled coil in Figures 10 and 11, respectively. In the case of the single helix, the interaction energy of residue E2K with the rest of the helix shows only a slight decrease with λ, i.e., mutating Glu to Lys does not

Wild Type and Mutated R-Keratin Fragments

affect significantly the interaction energies inside the single helix. The interactions of residue E2K with water are dominated by electrostatic energy. For λ ) 0 and 1 corresponding to a charge of -1 and +1, respectively, the potential energy is about -500 kJ mol-1. At λ ) 0.5, the charge is nearly zero, and consequently, the potential energy is almost zero. The separation of the various contributions to the free energy for the coiled coil system is shown in Figure 11. In contrast to the intramolecular interaction in the single helix system, the interactions of residue E2K within the coiled coil show large changes in the potential energy as a function of λ. Thus, the interaction of residue E2K with chain A shows an increase of the energy from -300 kJ mol-1 to about 10 kJ mol-1. The increase of the energy indicates a repulsive interaction of residue E2K with residues of chain A. This can be explained by the reduced mobility of the side chain of residue E2K that cannot freely rearrange and, consequently, experiences repulsive interactions with the neighboring side chains when it changes from a negative charge (Glu) to a positive charge (Lys). The interactions of residue E2K with chain B show a decrease of the potential energy from about +220 to -350 kJ mol-1 when λ changes from 0 (Glu) to 1 (Lys). These values indicate that Glu413 has a strong repulsive interaction with chain B, whereas Lys413 has an attractive interaction. Thus, it is clear that the main contribution to the free energy difference due to mutating Glu f Lys in the coiled coil comes from the interaction of Glu/Lys413 with the residues of chain B. The change in electrostatic interaction due to the mutation was already observed in the analysis of the interaction distance of the charged groups of chain A and B. Obviously, the attractive electrostatics dominate the interactions of the coiled coil and give the main contribution to the free energy difference. It should be noticed that the interaction of residue E2K with water results in a slight increase of the potential energy from -420 kJ mol-1 at λ ) 0 (Glu) to -350 kJ mol-1 at λ ) 1 (Lys). Considering the orientation of the side chains as previously discussed, this can be explained by the fact that the carboxyl group of Glu413 is totally exposed to the water while the amino group of Lys413 is involved in a salt bridge with Asp364. Consequently, it is less exposed to the solvent. Free Energy Calculations on Glu413Asp. Although the side chains of Glu and Asp differ only by a methylene group, the hybrid residue E2D (Figure 6) involves three dummy atoms. Annihilation and creation of these three heavy atoms is expected to cause large values of the derivative of the free energy with respect to λ at the beginning and at the end of the transformation. The calculated curve of the derivative of the potential energy with respect to λ is shown in Figure 12. It can be seen that extreme values at λ ) 0 and 1 are avoided due to the use of soft core potentials. The comparison to Figure 9 directly indicates completely different interactions in the system. During the Glu f Asp mutation electrostatics do not play a significant role, the mutation is characterized by Lennard-Jones interactions. However, during integration the large contributions from annihilation and creation of the atoms almost cancel. The

Biomacromolecules, Vol. 5, No. 6, 2004 2173

Figure 12. Ensemble average 〈dG/dλ〉 as a function of λ, obtained during the transformation Glu(-) f Asp(-) in the single helix and in the coiled coil (cutoffs 1.0/1.4 nm). Table 5. Calculated Free Energy Differences Corresponding to Mutating Glu(-) f Asp(-) in Single Helix and the Coiled Coila Twin-range cutoff

∆G1 (single helix) (kJ mol-1)

∆G2 (coiled coil) (kJ mol-1)

∆∆G (calc.) (kJ mol-1)

1.0/1.4 nm

0.8

-0.6

-1.3

a The double difference, ∆∆G, represents the relative free energy of association of wild type and mutated helices (chain A and chain B).

final results of the free energy calculations of the mutation Glu413Asp in the single helix and in the coiled coil are shown in Table 5. Mutating Glu413 to Asp results in small free energy changes in both single helix and coiled coil. In the single helix the free energy change due to the mutation is about 0.8 kJ mol-1 and in the coiled coil it is about -0.6 kJ mol-1 resulting in a relative free energy difference of -1.3 kJ mol-1. The error limits for this calculation are expected to be larger than the calculated value for ∆∆G of - 1.3 kJ mol-1. It can be concluded that the result of the free energy calculation on mutating Glu413 to Asp does not give evidence for a significant influence of this mutation on the stability of the coiled coil structure. Conclusions In this study, the conformational stability of different R-keratin intermediate filament fragments was investigated by means of molecular dynamics and free energy simulations. The IF fragments include both wild type and mutated forms (Glu413Lys and Glu413Asp) of the human hair monomeric IF unit that seem to be connected with a hereditary hair dystrophy, known as Monilethrix. Both approaches, the MD simulation and the free energy calculations, indicate that the mutation Glu413Lys stabilizes the coiled coil structure. The analysis of the electrostatic interactions and of the side chain dihedrals shows significant differences between the behavior of Glu413 and Lys413 in the coiled coil structure. The most relevant for this study is the interaction of Glu413 and Lys413 with the adjacent residue Asp364. The distance analysis of the relevant interaction between chain A and chain B namely Glu413-

2174

Biomacromolecules, Vol. 5, No. 6, 2004

Asp364 reveals strong repulsive interactions with an average interaction distance of 0.6 nm. In contrast, the equivalent analysis of the mutated coiled coil structure shows that Lys413 and Asp364 form a salt bridge during the whole simulation time exhibiting average interaction distance of 0.29 nm. These data suggest that Lys413 would stabilize the coiled coil structure, in comparison to the wild-type structure. Additionally, it has to be noted that due to the interaction with Asp364, the side chains of Glu413 and Lys413 have different orientations in the coiled coil structure. The glutamic acid side chain points away from the interface in order to minimize the repulsive interaction, while the lysine side chain points toward chain B because of attractive interactions to Asp364. The thermodynamics of these effects were further investigated by free energy calculations for the wild type and mutated IF fragments. The differences in electrostatic interactions cause a large change in free energy during transformation of Glu413 to Lys calculated by the thermodynamic integration approach. The association of the mutated chain A to chain B is favored by about 54 kJ mol-1 in comparison to the wild type coiled coil structure. The statistical errors associated with our free energy calculations might be slightly higher then presented due to the limitation in performing unrestrained simulations of the isolated helical segments. The structural rearrangement that is observed seem to be necessary to adapt the interactions in the mutated coiled coil. This may lead to changes in the IF assembly or its stability. The second mutation Glu413Asp does not lead to strong changes in the electrostatic interactions of the system. It is dominated by changes in the Lennard-Jones interactions caused by annihilation and creation of atoms during the amino acid transformation. The calculated difference in free energy ∆∆G ) -1.3 kJ mol-1 is within the error limits of the simulations. Therefore, it has to be concluded that this mutation does not affect the coiled coil stability. Based on these results, speculations about the effects of the mutations on the IF structure, probably connected to the clinical expression of Monilethrix, are possible. It is generally accepted that the IF assembly starts with a gradual association of the dimers by several distinct modes of lateral interactions that leads to the formation of tetramers. In vitro data12 indicate that the YRKLLEGE413E motif is crucial for the formation of authentic tetrameric complexes and also for the control of the filament width during assembly. The results of the present study that was carried out on an IF segment exhibiting a similar sequence indicate a rearrangement of the side chains when Glu413 is replaced by Lys. Moreover, due to the large value of the calculated free energy difference induced by the mutation, one may presume that larger structural arrangements might take place at the C-terminus of 2B segment. Due to the orientation of the side chains in the wild type dimer, the side chain of Glu413 may be involved in lateral interactions with neighboring IF dimers. According to our simulations, the conformational changes of the side chains at the mutation site may alter the higher order assembly of keratin IF due to the Glu413Lys mutation. However, the molecular basis of the structural changes in keratin IF caused by Monilethrix can only be explained when

Danciulescu et al.

detailed data at atomic resolution of the dimer-dimer associations of the IF become available. Acknowledgment. The authors gratefully acknowledge generous amounts of computer time provided by the Rechenzentrum of the RWTH Aachen University. We also greatly appreciate valuable discussion with PD Dr. J. Schweizer from the German Cancer Research Center, Heidelberg. References and Notes (1) Feughelman, M. Mechanical properties and structure of alpha-keratin fibres; UNSW Press: Sydney, 1997. (2) Parry, D. A. D.; Steinert, P. M. Q. ReV. Biophys. 1999, 32, 99. (3) Watts, N. R.; Jones, L. N.; Cheng, N.; Wall, J. S.; Parry, D. A. D.; Steven, A. C. J. Struct. Biol. 2002, 137, 109. (4) Winter, H.; Rogers, M. A.; Langbein, L.; Stevens, H. P.; Leigh, I. M.; Labreze, C.; Roul, S.; Taieb, A.; Krieg, T.; Schweizer, J. Nat. Genet. 1997, 16, 372. (5) Horev, L.; Glaser, B.; Metzker, A.; Ben-Amitai, D.; Vardy, D.; Zlotogorski, A. Human Heredity 2000, 50, 325. (6) Winter, H.; Rogers, M. A.; Gebhardt, M.; Wollina, U.; Boxall, L.; Chitayat, D.; Babul-Hirji, R.; Stevens, H. P.; Zlotogorski, A.; Schweizer, J. Human Genet. 1997, 101, 165. (7) Winter, H.; Labreze, C.; Chapalain, V.; Surleve-Bazeille, J. E.; Mercier, M.; Rogers, M. A.; Taieb, A.; Schweizer, J. J. InVest. Dermatol. 1998, 111, 169. (8) Steinert, P. M.; Yang, J. M.; Bale, S. J.; Compton, J. G. Biochem. Biophys. Res. Commun. 1993, 197, 840. (9) Kammerer, R. A.; Schulthess, T.; Landwehr, R.; Lustig, A.; Engel, J.; Aebi, U.; Steinmetz, M. O. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 13419. (10) Steinmetz, M. O.; Stock, A.; Schulthess, T.; Landwehr, R.; Lustig, A.; Faix, J.; Gerisch, G.; Aebi, U.; Kammerer, R. A. EMBO J. 1998, 17, 1883. (11) Wu, K. C.; Bryan, J. T.; Morasso, M. I.; Jang, S.-I.; Lee, J.-H.; Yang, Jun-Mo; Marekov, L. N.; Parry, D. A. D.; Steinert, P. M. Mol. Biol. Cell 2000, 11, 3539. (12) Herrmann, H.; Strelkov, S. V.; Feja, B.; Rogers, K. R.; Brettel, M.; Lustig, A.; Haner, M.; Parry, D. A. D.; Steinert, P. M.; Burkhard, P.; Aebi, U. J. Mol. Biol. 2000, 298, 817. (13) Gounari, F.; Merdes, A.; Quinlan, R.; Hess, J.; FitzGerald, P. G.; Ouzounis, C. A.; Georgatos, S. D. J. Cell Biol. 1993, 121, 847. (14) Merdes, A.; Gounari, F.; Georgatos, S. D. J. Cell Biol. 1993, 123, 1507. (15) Knopp, B. Dissertation RWTH, Aachen, 1995. (16) Knopp, B.; Jung, B.; Wortmann, F.-J.; Hoecker, H. Macromol. Symp. 1994, 81, 377. (17) Knopp, B.; Jung, B.; Wortmann, F.-J. Macromol. Symp. 1996, 102, 175. (18) Knopp, B.; Jung, B.; Wortmann, F.-J. Macromol. Theory Simul. 1996, 5, 947. (19) Knopp, B.; Jung, B.; Wortmann, F.-J. Macromol. Theory Simul. 1997, 6, 1. (20) Daggett, V.; Levitt, M. J. Mol. Biol. 1992, 223, 1121. (21) Daggett, V.; Levitt, M. J. Mol. Biol. 1993, 232, 600. (22) Hirst, J. D.; Brooks, C. L., III. Biochemistry 1995, 34, 7614. (23) Nick, B.; Kastanja, I. N.; Schweizer, J.; Wortmann, F.-J. Abstracts of Papers, 223rd ACS National Meeting, Orlando, FL, 2002. (24) Nick, B. Habilitation RWTH, Aachen, 2002. (25) Kastanja, I.; Nick, B.; Schweizer, J.; Wortmann, F.-J. DWI Rep. 2000, 123, 496. (26) Rogers, M. A.; Langbein, L.; Praetzel, S.; Moll, I.; Krieg, T.; Winter, H.; Schweizer, J. Differentiation 1997, 61, 187. (27) Rogers, M. A.; Winter, H.; Langbein, L.; Wolf, C.; Schweizer, J. J. InVest. Dermatol. 2000, 114, 464. (28) Rogers, M. A.; Winter, H.; Wolf, C.; Heck, M.; Schweizer, J. J. Biol. Chem. 1998, 273, 26683. (29) Korge, B. P.; Hamm, H.; Jury, C. S.; Traupe, H.; Irvine, A. D.; Healy, E.; Birch-Machin, M.; Rees, J. L.; Messenger, A. G.; Holmes, S. C.; Parry, D. A. D.; Munro, C. S. J. InVest. Dermatol. 1999, 113, 607. (30) Guex, N.; Peitsch, M. C. Electrophoresis 1997, 18, 2714. (31) O’Shea, E. K.; Klemm, J. D.; Kim, P. S.; Alber, T. Science 1991, 254, 539. (32) Kastanja, I. Personal communications. (33) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Model. 2001, 7, 306.

Wild Type and Mutated R-Keratin Fragments (34) Berendsen, H. J.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. In Intermolecular forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (35) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger, P. H.; Kru¨ger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular simulation: the GROMOS96 manual and user guide, Hochschulverlag AG an der ETH Zu¨rich. (36) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (37) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463. (38) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952. (39) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (40) Wong, C. F.; McCammon, J. A. J. Am. Chem. Soc. 1986, 108, 3830. (41) Bash, P. A.; Singh, U. C.; Brown, F. K.; Langridge, R.; Kollman, P. A. Science 1987, 235, 574. (42) Kuczera, K.; Gao, J.; Tidor, B.; Karplus, M. Proc. Natl. Acad. Sci. U.S.A. 1990, 87, 8481. (43) Simonson, T.; Brunger, A. T. Biochemistry 1992, 31, 8661.

Biomacromolecules, Vol. 5, No. 6, 2004 2175 (44) Dang, L. X.; Merz, K. M., Jr.; Kollman, P. A. J. Am. Chem. Soc. 1989, 111, 8505. (45) Tidor, B.; Karplus, M. Biochemistry 1991, 30, 3217. (46) Prevost M.; Wodak S J.; Tidor B.; Karplus M. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 10880. (47) Sun Y. C.; Veenstra D. L.; Kollman P. A. Protein Eng. 1996, 9, 273. (48) Tanimura, R.; Saito, M. Mol. Simul. 1996, 16, 75. (49) Sugita, Y.; Kitao, A. Biophys. J. 1998, 75, 2178. (50) van Gunsteren, W. F.; Beutler, T. C.; Fraternali, F.; King, P. M.; Mark, A. E.; Smith, P. E. In Computer Simulation of Biomolecular Systems, Volume 2; van Gunsteren, W. F., Weiner, P. K., Wilkinson, A. J., Eds.; Escom Science Publishers: Leiden, 1993. (51) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Chem. Phys. Lett. 1994, 222, 529. (52) Zacharias, M.; Straatsma, T. P.; McCammon, J. A. J. Chem. Phys. 1994, 100, 9025.

BM049788U