Molecular Dynamics Simulation of Protein Adsorption at Fluid Interfaces

Sep 14, 2010 - School of Life Sciences and International Centre for Brewing and Distilling, Heriot-Watt University,. Riccarton, Edinburgh, EH14 4AS, U...
1 downloads 0 Views 2MB Size
Download PDF
Biomacromolecules 2010, 11, 2781–2787

2781

Molecular Dynamics Simulation of Protein Adsorption at Fluid Interfaces: A Comparison of All-Atom and Coarse-Grained Models Stephen R. Euston* School of Life Sciences and International Centre for Brewing and Distilling, Heriot-Watt University, Riccarton, Edinburgh, EH14 4AS, United Kingdom Received July 27, 2010; Revised Manuscript Received August 24, 2010

The adsorption of LTP at the decane-water interface was modeled using all-atom and coarse-grained (CG) molecular dynamics simulations. The CG model (300 ns simulation, 1200 ns scaled time) generates equilibrium adsorbed conformations in about 12 h, whereas the equivalent 1200 ns simulation would take about 300 days for the all-atom model. In both models the LTP molecule adsorbs with R-helical regions parallel to the interface with an average tilt angle normal to the interface of 73° for the all-atom model and 62° for the CG model. In the all-atom model, the secondary structure of the LTP is conserved upon adsorption. A considerable proportion of the N-terminal loop of LTP can be found in the decane phase for the all-atom model, whereas in the CG model the protein only penetrates as far as the mixed water-decane interfacial region. This difference may arise due to the different schemes used to parametrize force field parameters in the two models.

Introduction The adsorption of proteins at interfaces is a phenomenon observed in many areas of biology. It is of particular technological interest in a number of industries, including the biomedical1 and food and beverage industries.2 In the latter, proteins are commonly used to stabilize fluid interfaces (e.g., oil-water and air-water interfaces) in, for example, food emulsions and foams.2 In this application, the proteins adsorb to the interface between water and the second nonaqueous phase, lower the interfacial tension, thus, facilitating interface formation, and create a viscous steric stabilizing layer that slows the rate of breakdown of the emulsion or foam structure. The conformation that the protein adopts at the interface is an important factor in the stability of the interface. The adsorbed conformation of a protein is determined by many things, some of which are related to the proteins native structure (presence of secondary and tertiary structure and the distribution of amino acids in the sequence), while others relate to the type of interface (e.g., hydrophobicity and presence of charged groups). For fluid interfaces in food systems, the oil-water or air-water interfaces at which the protein adsorbs are not charged, and the main driving forces for adsorption appear to be entropic in nature. These include dehydration of the oil surface by adsorbing amino acid residues and possibly entropic contributions associated with secondary and tertiary structural changes in the protein at the interface. Studying the conformation of adsorbed protein is not a trivial matter. It is possible to follow changes in secondary structure in adsorbed molecules, but the tertiary structure is more difficult because many of the experimental methods that are used for protein native structure determination (e.g., X-ray diffraction and NMR) either cannot be used or are more difficult to apply to interfaces. Consequently, detailed conformational information on individual adsorbed proteins at the tertiary structure level is lacking, with much of the available data confined to studies of the overall properties of a composite layer * To whom correspondence should be addressed. E-mail: s.r.euston@ hw.ac.uk.

obtained using indirect methods. One method that may be useful in studying adsorbed protein conformation is computer simulation.3,4 We have recently used a molecular dynamics (MD) simulation to probe the structural changes that occur in nonspecific barley lipid transfer protein (LTP) when it adsorbs to either a water-vacuum interface or a decane-water interface.5,6 Using a conventional MD simulation, we found that LTP adsorbed to both the water-vacuum and decane-water interfaces within the restricted time scale achieved, between 20 and 40 ns. The LTP molecule adsorbed to the water-vacuum interface but did not penetrate into the vacuum-space and did not undergo a significant conformational change from its native state.5 At the decane-water interface, however, adsorption was accompanied by a significant penetration of the molecule into the decane phase and by evidence of significant structural change during the simulation.5 These preliminary MD simulations were only run for a short time scale, up to 40 ns. Experimental studies suggest that any conformational change that does occur in adsorbed proteins will take place over much longer time spans.7 Therefore, it is not surprising that we did not see substantial conformational change in our simulations. There have been several recent advances in molecular dynamics methodology that may allow longer time scale dynamics of adsorbed proteins to be probed.8–13 These methods look to speed up simulations by overcoming energy barriers that limit the sampling of phase space by preventing transitions between minima on the molecules energy landscape. Often these methods require the use of biasing potentials or alternative coordinate systems which require careful selection for successful implementation. It is not always obvious how these are selected, and this can limit the application of accelerated methods. An alternative methodology that promises to be useful in simulating protein adsorption dynamics and one that is far easier to implement is coarsegraining.14 In this method a simplified representation of molecules is introduced where groups of atoms are represented as composite beads. For coarse-grained (CG) schemes such as the popular Martini force-field, introduced by Marrink and co-workers,14–16 there is a 4:1 mapping of atoms/CG beads. A

10.1021/bm100857k  2010 American Chemical Society Published on Web 09/14/2010

2782

Biomacromolecules, Vol. 11, No. 10, 2010

similar treatment can be applied to water molecules. This leads to an immediate 4-fold reduction in the number of atoms to be simulated in any simulation and, thus, a faster simulation. For the Martini approach, this is combined with the use of softinteraction potentials, which allow much larger integration time steps to be used, typically 20-30 fs compared to 2 fs for conventional MD with GROMACS, which allow a further speed-up of the simulation run. This methodology has been used successfully to study proteins embedded in lipid membranes,15 which suggests it will be applicable to proteins adsorbed at other fluid interfaces. In this study we compare an all-atom MD simulation of LTP adsorption to a decane surface carried out over several 100 s of ns, with a CG model for LTP at a decane surface. We are particularly interested in finding out if the two methodologies give comparable results. Our rationale for this is that, whereas a typical all-atom simulation for the size of system we are using runs at approximately 1.0 ns per day per processor, we can achieve around 150 ns per day per processor for a comparably sized CG simulation. Obviously, if the coarse-grained simulations give similar results, then this would allow us to access the long time scale dynamics of adsorbed protein molecules. In addition, barley LTP has technological significance to the brewing industry. It is a major contributor to foam formation and stability in beer,17,18 and thus, studies of the adsorption of LTP will have relevance to this important quality factor in alcoholic beverages. Simulation Methodology. Molecular dynamics simulations were carried out using the Gromacs MD package version 4.0.5.19 The adsorption of LTP at the decane-water interface was studied using an all-atom description of the LTP conformation and a coarse-grained (CG) description using the 4:1 mapping protocol of the Martini force field, as proposed by Marrink and co-workers.14–16 The structure for the barley LTP was that proposed by Heinemann et al.20 For the all-atom MD simulations, the methodology used (with some modifications) has been described elsewhere.5,6 Briefly, a simulation box was set up with dimensions of 6.5 × 6.5 × 5 nm and filled with 4923 SPC water molecules and an LTP molecule at the center. The box was expanded by 2 nm at each z-face to form a vacuum space, and this was filled with 408 decane molecules. The number of water and decane molecules was chosen by trial and error to give a system after equilibration with a water density close to 1000 g/L and a decane density close to 720 g/L. The system was equilibrated for 2 ns using an NPT ensemble (P ) 1 bar, T ) 300 K, Berendsen barostat and thermostat21) and for a further 2 ns using a NVT ensemble (T ) 300 K, Berendsen thermostat), both with the positions of the heavy atoms of the LTP molecule position restrained. Finally, a production run of at least 300 ns was carried out within the NVT ensemble at 300 K using a Berendsen thermostat. For the all-atom simulations, the nonbonded interactions in the system were described by the GROMOS98 force field (G43a1). Electrostatic interactions were calculated using the particle mesh Ewald summation method (PME)22,23 with a switching function to decay the Lennard-Jones interactions to zero after a cutoff of 0.9 nm. At separations greater than 0.9 nm, the interactions decay exponentially to zero at 1.2 nm. Linear center of mass motion was removed from all atoms over the course of the simulation. For coarse-grained simulations, the LTP starting conformation was converted to a CG representation of the structure using the atom2cg awk script and the seq2itp perl script, both of which were downloaded from the Martini website. In its current form, the Martini force field requires that the secondary structure of

Euston

the protein be defined at the start of the simulation and remain fixed throughout. The secondary structure of the native LTP was determined from the starting structure using the DSSP software.24 This data was used to define the topology of the coarse-grained LTP, along with the positions of the four disulfide bonds, in the seq2itp script. The coarse-graining of the protein structure in Martini involves a 4 to 1 mapping of atoms. In other words, four atoms in the protein structure are replaced with a single bead of composite mass. For simplicity, only four general bead types are defined, corresponding to polar, nonpolar, apolar, and charged beads. Each type has a number of subtypes that have differing degrees of polar, nonpolar, apolar, and charged character. This allows for a more realistic representation of the chemistry of the defined molecule. The force field parameters for CG beads have been derived by optimization based on the reproduction of partition-free energies between polar and apolar phases of a large set of common compounds. The Martini force field also incorporates a coarse-grained water model and a coarse-grained decane model. For water there is also a 4:1 mapping, with each CG water being equivalent to four all-atom water molecules, while for decane a three bead representation (3.333:1 mapping) was found to give a better fit to the experimental partition free energies used to parametrize the model. The Martini force field also uses soft potentials for the nonbonded interactions forces. This allows much larger time steps to be used (30 fs was used in this work for the CG simulations compare to 2 fs for the all atom simulations). The combination of smaller system size and longer time steps means that the CG simulation runs much faster than the all atom simulation, albeit at the expense of structural complexity. To set up the CG simulation, a box of 7.5 × 7.5 × 7.5 nm was defined and filled with CG water and a single LTP molecule at the center. The box was expanded along the z-axis to a dimension of 7.5 × 7.5 × 15 nm, which formed two vacuum spaces of 3.75 nm in the z-dimension. The vacuum spaces were filled with CG decane. The numbers of CG water and CG decane molecules were chosen by trial and error to give a final density of water and decane close to 1000 and 720 g/L, respectively, after equilibration. The system was equilibrated for 30 ns using an NPT ensemble (P ) 1 bar, T ) 300 K, Berendsen barostat and thermostat) and for a further 30 ns using a NVT ensemble (T ) 300 K, Berendsen thermostat), both with the positions of the CG beads of the LTP molecule position restrained. For the CG simulations, the nonbonded interactions in the system were described by the Martini force field. Electrostatic interactions were calculated using a shift function with a cutoff of 1.2 nm. Linear center of mass motion was removed from all atoms over the course of the simulation. The justification for the different summation methods for the CG system is that the electrostatic interactions in the Martini force field are not very accurate in the first place, as the screening in the system is set to be uniform across the system with a screening constant of 15, and the force field has been parametrized with short-range shifted electrostatic interactions. Production runs were carried out for 300 ns (approximately equivalent to 1200 ns CG time). The structures of the simulated adsorbed and nonadsorbed LTP conformations were characterized in a number of ways. The root-mean-square deviation (rmsd) and radius of gyration (Rg) were used as an indicator of changes in the tertiary structure in the conformations as a function of time. Secondary structure changes were followed by calculating the number of amino acid residues existing in helices, β-sheet, or in random coils as a function of simulation time using DSSP.24 Density profiles of water, decane, and protein across the simulation cell were

Molecular Dynamics Simulation of Protein Adsorption

Biomacromolecules, Vol. 11, No. 10, 2010

2783

Figure 1. Radius of gyration (Rg) as a function of simulation time for (a) all-atom LTP in a water box, (b) all-atom LTP at a decane-water interface, (c) CG LTP in a water box, and (d) CG LTP at a decane-water interface.

Figure 2. Root mean square displacement (rmsd) as a function of simulation time for (a) all-atom LTP in a water box, (b) all-atom LTP at a decane-water interface, (c) CG LTP in a water box, and (d) CG LTP at a decane-water interface.

calculated using analysis algorithms from the GROMACS software.19 An estimate of the difference in protein conformational entropy between adsorbed and nonadsorbed all-atom LTP was calculated using the method outlined below. The Rg was used as a conformational marker, and a histogram of the Rg versus probability of that Rg (pi(Rg)) occurring was assumed to give a realistic estimate of the conformational probability density for the adsorbed and nonadsorbed LTP molecules. If we assign the Rg to one of n equally spaced bins, then the conformational entropy S is calculated from the Gibbs entropy equation

stress that interpretation of the time in CG simulations should be done with care. The mean Rg averaged over the last 100 ns for the all-atom LTP in a water box (Figure 1) is 1.17 ( 0.01 nm. This is slightly lower than the experimental Rg of 1.23 ( 0.01 nm reported by Da Silva et al.27 for tobacco LTP1 but comparable to the molecular dynamics simulation of Lai et al.28 who report a simulated Rg of 1.17 nm for rice LTP. The sequence homology between rice LTP and barley LTP is about 63%, that is, 57 of the 91 residues are conserved. This includes all eight cysteines, most of the hydrophobic residues in the binding pocket, prolines and glycines in the turns and loops, and two tyrosines at the C- and N-terminal ends.29 In addition the structures of barley20 and rice,30 LTP show a similar fold that differs only in the details.31 It is therefore not surprising that simulated Rg values are the same. Tobacco LTP shares a similar tertiary structure fold to barley and rice LTP27 but has structural differences that alter the way in which it binds lipids in comparison to cereal LTPs. In particular, tobacco LTP contains a hydrophobic cluster in its structure that is absent from cereal LTPs, and this limits binding in the hydrophobic cavity to just one lipid, whereas cereal LTPs are able to bind one or two lipid molecules.27 This may partly explain the difference between the simulated Rg for barley and rice LTP and the experimental value for tobacco LTP. The first thing we notice about the rmsd curves (Figure 2) is that for both LTP models in water it takes a long time for the rmsd to settle down to equilibrium. If we look at the rmsd for CG LTP in a water box (Figure 2) we see that the CG conformation reaches a relatively constant value after about 800 ns of CG time. For the all-atom conformations in water, we can see that the rmsd is higher than for the CG conformation and appears to reach equilibrium by 300 ns. Change in LTP conformation during the bulk water simulations is most probably due to the fact that the systems were equilibrated (NVT and NPT simulations) with the positions of the LTP constrained, and once these were relaxed, it took some time for the LTP conformation to equilibrate. LTP is known to be a very stable protein, so it is perhaps not a surprise that conformational relaxation takes place over a long time scale. LTP is reported to have a denaturation temperature of 100 °C and to be able to withstand extensive boiling during the brewing process with only partial denaturation.32 This suggests a highly stable structure that is resistant to conformational change. The differences observed between the rmsd of the all-atom and CG conformations almost certainly arise due to the nature of the model. The apparent longer time that it takes for equilibration of the CG structure in water may be due to the extra restraints

n

S ) -kb

∑ pi(Rg)ln(pi(Rg))

(1)

i)1

Obviously the calculated conformational entropy is only an estimate because we are assuming that the same Rg value indicates the same unique conformation, which is not necessarily true. However, for our purposes, the estimate should be accurate enough to assess whether an adsorbed simulated LTP molecule is able to access more conformations and is therefore more flexible than LTP in water. The hydrophobicity along the LTP sequence was calculated using the amino acid hydrophobicity weighting scale proposed by Kyte and Doolittle,25 as implemented on the ExPASy proteomics server26 (http://expasy.org/tools/protscale.html).

Results and Discussion Some idea of the changes in tertiary structure of the LTP models, both at interfaces and in bulk water can be deduced by looking at the root-mean-square deviation (rmsd) of the atomic positions and the radius of gyration (Rg) for both models. Figures 1 and 2 show the radius of gyration (Rg) and root-mean-square deviation (rmsd) of LTP in solution and adsorbing to the decane-water interface for both all-atom and CG models. A direct comparison of rmsd and Rg between our simulated allatom and CG models is not valid due to the differing representations used for the chains. However, comparison between adsorbed and nonadsorbed conformations for the same model is valid. In keeping with the suggested 4-fold scaling in time in Martini simulations15 (based on a comparison of the experimental self-diffusion coefficient of water and that calculated from CG simulations) each CG time point is multiplied by 4 so that the full 300 ns CG simulation is assumed to be equivalent to a 1200 ns all-atom simulation. This scaled time will be referred to as CG time. Monticelli et al.,15 however, do

2784

Biomacromolecules, Vol. 11, No. 10, 2010

that are placed on the secondary structure in the CG model. The secondary structure is held fixed in the CG model, while in the all-atom model it is free to rearrange and fluctuate under the forces that the individual atoms experience. The Rg for both LTP models in water (Figure 1) remains relatively constant throughout the course of the simulations, even though there are noticeable changes in the rmsd. This would suggest that the changes in rmsd in Figure 2 are not due to changes in molecular shape (i.e., unfolding) but are due to internal rearrangements of the molecule, which take place on a longer time scale. Even though there are observable differences between the ways in which the different LTP models behave in a water box, it is clear from Figures 1 and 2 that the adsorbed molecules adopt a different, more open conformation to that adopted by the native structure. This occurs in both the CG and all-atom models, although the unfolding of the CG conformation is clearly more extensive where there is a definite increase in Rg as the LTP molecule adsorbs and unfolds at the decane interface. If the Rg is decomposed into its x, y, and z components (relative to the principle axes), it is found that there is a significant decrease in Rgz normal to the interface and a significant increase in the Rgx and Rgy parallel to the interface, that is, the molecule flattens at the surface. It is also noticeable from Figure 1 that the CG conformation in the decane-water system has a lower starting Rg than the CG conformation in a water box. This could be due to a perturbing effect of being confined between the two decane slabs. It was noticed that during the equilibration steps prior to the production run that the LTP molecule did move closer to the decane interface even though position restraints were applied. This led to a conformational change that gave a lower overall Rg value for the staring conformation at the decane-water interface. It is not considered that this will affect the final adsorbed conformations for the CG LTP at the decane-water interface since conformational changes of this type are likely to occur as a natural consequence of an adsorbing protein approaching an interface. The fact that this difference in starting Rg was not seen for the all-atom model is likely to be a consequence of the slower dynamics in the all-atom model, that is, the equilibration period was not long enough for the conformational change to occur. To put these differences in perspective, the starting Rg for the CG LTP in water is about 5% higher than the Rg for the CG LTP in the decane-water system, and the final Rg is 5% higher for the adsorbed conformation. So these differences are relatively small, but are significant compared to the differences between the replicates for the CG simulations. The low degree of unfolding, that is, close to a 10% increase in Rg on adsorption for the CG LTP molecule is consistent with the high conformational stability of the molecule. It is more difficult to estimate the degree of unfolding for the all-atom model because this Rg probably has not reached a steady state. From the data available, the change in Rg from the starting conformation to the final adsorbed conformation is about 9%, which is comparable to the value obtained with the CG model. Barley LTP has a structure characterized by four R-helices (labeled A-D from the Nterminal end) connected by loops of random coil structure.20 Three of the helices, A-C form an up-down-up motif, where helix A takes up a central position with B and C on either side. The overall structure is stabilized by four disulfide bonds which connect the helices. Cys27 and Cys13 connect helix A and B, Cys3-Cys50 connects A to C, and Cys28-Cys73 connects B to D. At the C terminal end residues Asn74-Tyr91 adopt a random coil extended configuration. This contains an extended loop structure between Val75-Ile90, which forms hydrophobic con-

Euston

tacts with all four helices and is stabilized by a disulfide bond between Cys48-Cys87. This structure imparts a high degree of conformational stability on the LTP molecule, and it is not perhaps surprising that the conformational change upon adsorption is low. LTP interaction with phospholipid monolayers has been studied experimentally by Subirade and co-workers33,34 for a nonspecific wheat and maize LTP using attenuated total reflection infrared spectroscopy (ATR). On binding of both LTP’s to the phospholipid interface there was evidence of a significant change in the ATR spectrum, particularly in the amide I region, which indicates conformational change. The LTP conformation at the interface also appears to be surface pressure dependent. Changes in secondary structure during adsorption of proteins have been observed before, although this is highly dependent on the protein molecule and on the interface type. At the oil-water interface the milk protein β-lactoglobulin (β-lac) undergoes a surface-induced unfolding of its β-barrel structure which then reforms intermolecularly.35 The small R-helix region present in β-lac on the other hand increases in proportion particularly when the molecule is adsorbed to a hydrophobic surface.36–38 In contrast, Lad et al.39 report little change in secondary structure when β-lac adsorbs to the air-water interface. They also report no change in secondary structure for bovine serum albumin at the air-water interface, but lysozyme adopted a non-native adsorbed conformation containing a significant proportion of antiparallel β-sheet.39 Subirade et al.33,34 also probed the change in secondary structure on adsorption of LTP to the phospholipid interface. At high surface pressure ATR spectral changes suggested an increase in R-helical content, which was not observed at low surface pressure. Overall, the LTP did not undergo a surface induced unfolding of the helical regions. Subirade et al.33,34 concluded that the helical regions were stabilized upon adsorption. This is an important point with regards to our simulations since it justifies the restraining of the secondary structure that is used in the CG model. To check that the secondary structure of the adsorbed and nonadsorbed LTP in the all-atom model was also conserved we calculated the secondary structure using the DSSP program.24 This is plotted in Figures 3a for simulated all-atom LTP in water and in Figure 3b for the simulated allatom LTP at the decane-water interface. There are some fluctuations in the helical regions of the LTP molecule for both the adsorbed and nonadsorbed molecule, but in general adsorption does not lead to significant unfolding of the R-helical regions of the molecule. There is some loss/fluctuation in helix A for both the LTP in water and at the decane surface. In both systems helix A appears to be split into two smaller subhelices, one of which unfolds for part of the simulation. Both helices B and C are well conserved for the full 300 ns of the simulation. Helix D appears to be more stable in the decane-water system. This helix is found embedded in the decane phase during most of the simulation run. The observation of increased stability for helix D is consistent with the experimental conclusion that the R-helical structure of LTP is stabilized at a lipid interface. The C-terminal loop also shows differences in structure between the conformations in a water box and adsorbed at the decane-water interface. In water there is a tendency for the LTP molecule to form short β-sheet (particularly after about 100 ns), turn, and bend structures. When the LTP is adsorbed at the decane surface formation of β-sheet is not observed in the C-terminal loop, and the proportion of bends and turns is also reduced. In the adsorbed conformation the C-terminal loop

Molecular Dynamics Simulation of Protein Adsorption

Biomacromolecules, Vol. 11, No. 10, 2010

2785

Figure 3. Secondary structure assignment calculated using DSSP24 as a function of simulation time for (a) all-atom LTP in a water box and (b) all-atom LTP at a decane-water interface: (white) coil; (red) β-sheet; (black) β-bridge; (green) bend; (yellow) turn; (blue) R-helix; (purple) 5-helix; (brown) 3-helix.

Figure 4. Simulated adsorbed conformation for all-atom LTP adsorbed at the decane-water interface. (a) View looking parallel to the decanewater interface; (b) view looking perpendicular to the decane-water interface. The four helices and the region from residue 55-91 are shown in different colors. (red) helix A; (blue) helix B; (green) helix D; (lt. blue) residues 55-91 (excluding helix D and including C-terminal extended loop). Decane molecules are colored olive green. Water molecules have been omitted for clarity.

is found penetrating into the decane phase of the simulation box (Figure 4). Obviously this favors formation of a random coil structure over ordered secondary structure, possibly because the driving force for hydrogen bond formation is reduced in the nonaqueous environment. The C-terminal loop is known to form hydrophobic contacts with all four helices, which stabilizes the LTP structure.20 In a hydrophobic environment these are less likely to form and the C-terminal extended loop will become more disordered and flexible. There are both simulation40 and experimental41,42 evidence to support the penetration of adsorbed protein molecules into an oil phase. Anderson et al.40 used a Monte Carlo lattice model to study globular proteinlike molecules adsorbing at an oil-water interface. They found that the conformational entropy of the molecules increased upon adsorption, which indicated the chains were entering the oil-phase and that this allowed them to become more flexible and to adopt a larger number of conformations. This finding is supported by the experimental results of Murray et al.41,42 who have found that globular protein are more expanded at the oil-water interface compared to the air-water interface, probably because the oil is a better solvent for the hydrophobic regions of the molecule, thus allowing a greater penetration into the oil-phase.

We have calculated the difference in conformational entropy between adsorbed and nonadsorbed LTP using eq 1 for both all-atom and CG LTP and have found that there is an increase of +0.45kb for the all-atom model and +0.64kb for the CG model, that is, in both models the LTP molecule is able to occupy more conformations at the interface. The surface pressure has also been observed to influence the orientation of the helices. The orientation of the helices was determined experimentally by Subirade et al.33,34 from the dichroism spectra of the amide I adsorption band. At low surface pressure the helices tend to sit almost parallel to the interface with an estimated angle of 60-70° to the normal of the interface. At high surface pressure the results suggest that the helices are in a more random orientation. We have determined the angle of tilt averaged over all helices in both all-atom and CG simulations during the course of the simulations. Both models show a tendency for the helices to orient parallel to the decane interface with an average tilt angle for the all-atom LTP (calculated over the last 100 ns of the simulation) being 73° and the CG model slightly lower at 62°. These values agree well with those reported by Subirade et al.33,34 This similarity between experimental and simulation results is encouraging. It

2786

Biomacromolecules, Vol. 11, No. 10, 2010

Euston

Figure 6. Hydrophobicity as a function of position along primary sequence of barely LTP.20 The hydrophobicity was calculated using the ExPASy proteomics server26 applying the weighting scale proposed by Kyte and Doolittle.25

Figure 5. Density profiles normal to the decane-water interface (zdirection) for (A) all-atom LTP and (B) CG LTP: (a) density profile for LTP; (b) density profile for water; (c) density profile for decane.

should be remembered, however, that the experimental results are measured in the presence of other adsorbed protein molecules, which have may have a perturbative effect on the adsorbed protein structure. The simulation results are for an isolated protein, unperturbed by surrounding molecules. The experimental methods used by Subirade et al.34 cannot distinguish between the orientations of the individual helices. They have suggested34 that based on the findings that the helices are not exactly parallel in solution43 it is possible that helix D may lie perpendicular to the interface, while helices A-C are parallel. From our simulations helices A-D have individual tilt angles of 62, 76, 79, and 73° to the normal of the interface in the all atom simulation, and 74, 68, 43, and 62° in the CG simulation, which suggests that all helical regions prefer to lie parallel to the interface. The orientation of helices parallel to the surface is evident when looking at the final (300 ns) conformations for the adsorbed LTP, and in particular in the all-atom model (Figure 4). It is also evident from Figure 4 that in the all-atom model the LTP molecule penetrates a significant distance into the decane layer. For the CG LTP molecules this is not as obvious. In fact if we look at the density profile along the z-axis for the system (Figure 5) the penetration of the all-atom LTP into the decane layer is confirmed, while it would appear that the CG model only penetrates as far into the decane-layer as do the water molecules, that is, the protein residues are confined to the aqueous phase and the mixed-water decane region at the interface, but do not enter the pure decane phase. Subirade et al.34 note that their infrared spectroscopy results indicate that at low surface pressure at least part of the LTP molecule penetrates into the phospholipid layer. It was not possible for Subirade et al.34 to deduce which regions of the LTP molecule penetrate the oil phase. In our all-atom simulation we can see that the C-terminal region of the molecule, from residue 55-91, which includes helix D, makes up the major part of the molecule in the decane phase (Figure 4). In the CG model the same region, although it does not penetrate into the decane phase, is in contact with the surface. A plot of the hydrophobicity of the LTP

molecule, based on criteria proposed by Kyte and Doolittle25 shows that the C-terminal end of the molecule from about residue 55 is on average more hydrophobic than the region from the N-terminal up to residue 55 (Figure 6). This combined with the greater flexibility of the terminal loop makes it more likely for this region to penetrate the oil phase. Concerning the lack of penetration of the terminal loop region into the decane phase in the CG simulations, it should be noted that coarse-graining of a protein structure by its very nature leads to a reduction in the conformational entropy of a simulated protein molecule. In the Martini force field this is compensated by reducing the enthalpy contribution to the free energy used to parametrize the force field, which leads to a bias in the entropy/enthalpy balance at the CG level.14,15 It is possible that this accounts for the differing behavior of the terminal loop in the all-atom and CG models and that the partitioning of the C-terminal loop into the decane phase is entropy driven, a factor that will be reduced in the coarse-grained model. Evidence for the disordering effect of the decane phase on the C-terminal end of the all-atom LTP conformation can be seen in the DSSP secondary structure plots of Figure 3. In Figure 3a there is clear evidence of intermittent secondary structure formation (small regions of β-sheet, turn and bends) at the C-terminal loop for the all-atom LTP simulated in water. This is much reduced when the LTP is adsorbed to a decane interface, and the C-terminal loop penetrates into the decane and adopts a more random coil structure. It is clear from the results presented in this study that molecular dynamics simulation is capable of generating adsorbed conformations of LTP molecules that are consistent with the experimental studies on LTP adsorbed conformation at phospholipid interfaces. This is true for the all-atom model and to a certain degree for the CG models. The rmsd and Rg clearly show that the LTP molecule unfolds and flattens at the decane-water interface. In particular, both models reproduce the preferred parallel to the interface alignment of the R-helical regions of the LTP molecule. There are two major differences between the all-atom and the CG adsorbed conformations for LTP. In the CG model the LTP model seems to form an equilibrium unfolded conformation at the decane surface, where the rmsd and Rg reach a constant value for much of the second half of the simulation. In the all-atom simulation the rmsd and Rg still appear to be increasing after 400 ns. This confirms the expectation that the conformational dynamics and unfolding will be speeded up in the CG simulations. Taking in to account the anticipated increase in time scale for the CG simulation, a 300 ns CG simulation would be roughly equivalent to a 1200 ns all-atom simulation. A visual comparison of the rmsd of the all-atom and CG simulations would suggest that the pathway

Molecular Dynamics Simulation of Protein Adsorption

of unfolding is similar over the first 400 ns. We might expect further unfolding of the all-atom LTP conformation to occur if the simulation were to run for longer. The second difference between the adsorbed conformations generated in the two models is that the CG LTP C-terminal extended loop does not penetrate into the decane phase in the same way that it does in the all-atom model. This is most likely due to the differences in the way in which the force fields are parametrized and to the loss of conformational entropy in the coarse-graining process. The major advantage of the CG model is the greatly reduced CPU time needed. A 300 ns CG simulation (1200 ns CG time) takes about 12 h on four processors, while a 300 ns all-atom simulation takes over 3 months. The latter is clearly less unfolded than the CG simulation, and to generate a 1200 ns all-atom simulation, which we would hope would have a similar degree of unfolding to the CG simulation, would take about a year of CPU time. Clearly simulation techniques such as CG allow access to unfolding time scales that are difficult to achieve using all-atom models, and will allow us to probe the dynamics of adsorbed protein conformational changes. A potential problem is that we have seen clear differences between the all-atom and CG adsorbed conformations, which is surprising given the success that the CG model has had in describing dynamics of membrane embedded proteins. This should be weighed against the fact that the CG and all-atom models also show some obvious similarities such as the parallel alignment of helices to the decane-water interface. It is possible that the parametrization of the Martini force field, while optimized for membrane proteins is not optimal for surface adsorbed proteins and that improvements could be made to the similarity between all-atom and CG models by revisiting the CG force field parametrization for the adsorbed models. However, given that the original parametrization of the Martini force field is based on partitioning of the individual amino acids between water and a hydrocarbon phase (butane or cyclohexane) it is difficult to see how this improvement could be achieved.

Conclusions CG modeling of adsorbed protein conformation would appear to be a useful way of accessing long-time scale dynamics of protein adsorption. It has the advantage over other accelerated MD methods that it is very easy to set up and run and requires no biasing potentials to force the simulated protein to sample more of conformation space. It should be borne in mind that the coarse-graining procedure may produce artifacts (possibly through changes in the conformational entropy of the adsorbed conformation), which lead to differences between all-atom and CG adsorbed conformations. Further study of the adsorption of other protein molecules at an oil-water interface is required to see whether this is a general phenomenon or is specific to LTP. Further work is in progress for other globular proteins such as the bovine milk protein β-lactoglobulin to test this. Acknowledgment. The author acknowledges Heriot-Watt University Computer Centre for use of the HWU central HPC facility.

References and Notes (1) Chen, H.; Yuan, L.; Song, W.; Wu, Z.; Li, D. Prog. Polym. Sci. 2008, 33, 1059–1087. (2) McClements, D. J. In Modern Biopolymer Science; Kasapsis, S., Norton, I., Ubbink, J., Eds.; Elsevier Academic Press: London, 2009; pp 129-166.

Biomacromolecules, Vol. 11, No. 10, 2010

2787

(3) Euston, S. R. Curr. Opin. Colloid Interface Sci. 2004, 9, 321–327. (4) Euston, S. R.; Costello, G.; Naser, Md., A.; Nicolosi, M. In Understanding and Controlling the Microstructure of Complex Foods; McClements, D. J., Ed.; CRC Press: Cambridge, 2007; pp 334-379. (5) Euston, S. R.; Hughes, P.; Naser, Md A.; Westacott, R. Biomacromolecules 2008, 9, 1443–1453. (6) Euston, S. R.; Hughes, P.; Naser, Md A.; Westacott, R. Biomacromolecules 2008, 9, 3024–3032. (7) Beverung, C. J.; Radke, C. J.; Blanch, H. W. Biophys. Chem. 1999, 81, 59–80. (8) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562– 12566. (9) Spiwok, V.; Lipovova, P.; Kralova, B. J. Phys. Chem. B 2007, 111, 3073–3076. (10) Kubitzki, M. B.; de Groot, B. L. Biophys. J. 2007, 92, 4262–4270. (11) Hayward, S.; de Groot, B. L. Methods Mol. Biol. 2008, 443, 89–106. (12) Perez, D.; Uberuaga, B. P.; Shim, Y.; Amar, J. G.; Voter, A. F. Annu. Rep. Comput. Chem. 2009, 5, 79–98. (13) Abrams, C. F.; Vaden-Eijden, E. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 4961–4966. (14) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. J. Phys. Chem. B 2007, 111, 7812–7824. (15) Monticelli, L.; Kandasamy, S. K.; Periole, X.; Larson, R. G.; Tieleman, D. P.; Marrink, S. J. J. Chem. Theory Comput. 2008, 4, 819–834. (16) Lopez, C. A.; Rzepiela, A. J.; de Vries, A. H.; Dijkhuizen, L.; Hunenberger, P. H.; Marrink, S. J. J. Chem. Theory Comput. 2009, 5, 3195–3210. (17) Sørensen, S. B.; Bech, L. M.; Muldbjerg, M.; Beenfeldt, T.; Breddam, K. MBAA Tech. Q. 1993, 30, 136–145. (18) Evans, D. E.; Sheehan, M. C. J. Am. Soc. Brew. Chem. 2002, 60, 47–57. (19) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435–447. (20) Heinemann, B.; Andersen, K. V.; Nielsen, P. R.; Bech, L. M.; Poulsen, F. M. Protein Sci. 1996, 5, 13–23. (21) Berendsen, H. J. C.; Postma, J. P. M.; Di Nola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684–3690. (22) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089– 10092. (23) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. J. Chem. Phys. 1995, 103, 8577–8593. (24) Kabsch, W.; Sander, C. Biopolymers 1984, 22, 2577–2637. (25) Kyte, J.; Doolittle, R. F. J. Mol. Biol. 1982, 157, 105–132. (26) Gasteiger, E.; Hoogland, C.; Gattiker, A.; Duvaud, S.; Wilkins, M. R.; Appel, R. D.; Bairoch, A. In The Proteomics Protocols Handbook; Walker, J. M., Ed.; Humana Press: Totrowa, NJ, 2005; pp 571-607. (27) Da Silva, P.; Landoa, C.; Beltoise, R.; Ponchet, M.; Vovelle, F. Proteins: Struct., Funct., Bioinf. 2006, 64, 124–132. (28) Lai, Y.-T.; Cheng, C.-S.; Liu, Y.-N.; Liu, Y.-J.; Lyu, P.-C. Proteins: Struct., Funct., Bioinf. 2008, 72, 1189–1198. (29) Cheng, H.-C.; Cheng, P.-T.; Peng, P.; Lyu, P.-C.; Sun, Y.-J. Protein Sci. 2004, 13, 2304–2315. (30) Lee, J. Y.; Min, K.; Cha, H.; Shin, D. H.; Hwang, K. Y.; Suh, S. W. J. Mol. Biol. 1998, 276, 437–448. (31) Poznanski, J.; Sodanao, P.; Suh, W. S.; Lee, Y. J.; Ptak, M.; Vovelle, F. Eur. J. Biochem. 1999, 259, 692–708. (32) Lindorff-Larsen, K.; Winther, J. R. FEBS Lett. 2001, 488, 145–148. (33) Subirade, M.; Marion, D.; Pe´zolet, M. Thin Solid Films 1996, 284285, 326–329. (34) Subirade, M.; Salesse, C.; Marion, D.; Pe´zolet, M. Biophys. J. 1995, 69, 974–988. (35) Lefevre, T.; Subirade, M. J. Colloid Interface Sci. 2003, 263, 59–67. (36) Cornell, D. G. J. Colloid Interface Sci. 1982, 88, 536–545. (37) Cornell, D. G.; Patterson, D. L. J. Agric. Food Chem. 1989, 37, 1455– 1459. (38) Dufour, E.; Dalgalarrondo, M.; Adam, L. J. Colloid Interface Sci. 1998, 207, 264–272. (39) Lad, M. D.; Birembaut, F.; Matthew, J. M.; Frazier, R. A.; Green, R. J. Phys. Chem. Chem. Phys. 2006, 8, 2179–2186. (40) Anderson, R. E.; Pande, V. S.; Radke, C. J. J. Chem. Phys. 2000, 112, 9167–9185. (41) Murray, B. S.; Nelson, P. V. Langmuir 1996, 12, 5973–5976. (42) Murray, B. S. Prog. Colloid Polym. Sci. 1997, 103, 41–50. (43) Gincel, E.; Simorre, J. P.; Vovelle, F.; Marion, D.; Ptak, M. Eur. J. Biochem. 1994, 226, 413–422.

BM100857K