Enhancing Protein Adsorption Simulations by Using Accelerated Molecular Dynamics

Abstract

The atomistic modeling of protein adsorption on surfaces is hampered by the different time scales of the simulation ( s) and experiment (up to hours), and the accordingly different ‘final’ adsorption conformations. We provide evidence that the method of accelerated molecular dynamics is an efficient tool to obtain equilibrated adsorption states. As a model system we study the adsorption of the protein BMP-2 on graphite in an explicit salt water environment. We demonstrate that due to the considerably improved sampling of conformational space, accelerated molecular dynamics allows to observe the complete unfolding and spreading of the protein on the hydrophobic graphite surface. This result is in agreement with the general finding of protein denaturation upon contact with hydrophobic surfaces.

About the Authors:

Christian Mücksch

Affiliation: Physics Department and Research Center OPTIMAS, University of Kaiserslautern, Kaiserslautern, Germany

Herbert M. Urbassek

* E-mail: [email protected]

Affiliation: Physics Department and Research Center OPTIMAS, University of Kaiserslautern, Kaiserslautern, Germany

Introduction

The interaction of proteins with different biomaterials plays a vital role in describing their biocompatibility. Depending on the surface the overall protein structure may vary from native to unfolded [1]. Especially on hydrophobic surfaces proteins tend to unfold [2]. Experiments on pyrolytic graphite have shown that proteins irrespective of their primary sequence, secondary structure and molecular weight unfold and form nanopatterns [3]. Molecular dynamics (MD) simulations on the other hand provide a direct method for theoretically analyzing the adsorption processes at an atomistic scale.

A major problem that arises when using MD to study protein adsorption is the limited time scale of only a few hundred ns. The experimental time scale when these processes are investigated usually ranges from s to hours [4]. It has been shown by Wei et al. that very long simulations are necessary [5] since major adsorption steps like dehydration take a long time and denaturing events even longer. Strategies to overcome this sampling problem are needed. A computationally efficient method using an implicit inviscid water environment has first been established for protein adsorption processes by Raffaini et al. [6], [7]. But Sun et al. [8] argue by comparing implicit and explicit water models that these models show different adsorption behavior with stronger hydrophobic interactions in the case of implicit water.

Based on the work by Voter [9] Hamelberg et al. established an efficient method called accelerated molecular dynamics (in the following abbreviated as: aMD) for biomolecule simulation, which allows to escape the nanosecond time scale limitations [10]. Here, a biasing potential is added to the natural dynamics in order to escape potential energy minima. The method has been applied to many problems like bilayer structural and dynamic properties [11] or the dynamics of bovine pancreatic trypsin inhibitor (BPTI) [12]. Note that it is not straightforward to assign a ‘real time’ to the simulational time due to the bias in the dynamics, especially for rather complicated systems [13].

In the present work, a comparison between classical and accelerated MD is used to show differences in the sampling of conformational changes during adsorption on a hydrophobic graphite surface. To our knowledge, it is the first time that the technique of aMD is used to study protein adsorption. We note that interest in adhesion studies focuses on the final structure and adhesion forces of the adsorbed molecule, since this is where comparison to experiment is possible. The actual adsorption dynamics are not so relevant; this feature makes adsorption a good candidate for aMD.

Our protein of interest is the bone growth factor BMP-2 (bone morphogenetic protein 2) which is used to improve osteointegration by functionalizing implant surfaces [14], [15]. So far, several simulational studies of adsorption of BMP-2 have been performed, but all with classical MD [16]–[19]. To compare the results with previous studies [18], [19] graphite is used as model surface, which has many applications as an implant material [20] in the pyrolytic form.

Methods

Simulation Details

All the simulations were carried out using NAMD 2.8 [21] with the CHARMM27 force field [22] and the TIP3P [23] water model. The protein structure of BMP-2 was obtained from the protein data base (ID: 3BMP [24]) and placed in two different starting orientations differing by a rotation around an axis parallel to the surface so that the lowest protein atom had a distance of 6 Å to the first graphite layer (see Fig. 1).

[Figure omitted. See PDF.]

Figure 1. Two different initial orientations of BMP-2 above the graphite surface (indicated by 2 gray lines).

Nonpolar residues are shown in blue, neutral residues in green and polar residues in red. The sodium ions are shown in orange and the chloride ions in cyan. The explicit water surrounding is visualized as a volumetric surface in the background.

https://doi.org/10.1371/journal.pone.0064883.g001

For all simulations periodic boundary conditions are used with Particle Mesh Ewald (PME) [25] for long range electrostatics. The system is comprised of 43,951 water molecules with a salt concentration of 0.1 M NaCl. Since the protein structure was obtained at a pH of 5.4 and the pI (isoelectric point) of the protein is 4.7 [26] the protein net charge of was neutralized for correct usage of PME by having 2 more sodium ions than chloride ions in total. In order to prevent adsorption on the periodic image of the surface the system was chosen sufficiently large with a 25 Å thick layer of water on top of the protein. The cutoff for the van-der-Waals interactions was chosen as 12 Å. To describe the graphite surface [2 layers of (0001)] by Lennard-Jones interactions the needed parameters were taken from Werder et al. [27] to reproduce the experimental contact angle of 86° on graphite [28]. Since these parameters were developed for a different force field and water model we reproduced a 1 ns contact angle simulation with 9000 TIP3P water molecules. This resulted in a contact angle of 85.6° (see Fig. 1) which is in very good agreement with the results by Werder et al. Throughout all simulations the carbon atoms of the surface were kept fixed.

[Figure omitted. See PDF.]

Figure 2. Simulation of water droplet on graphite surface and fit (black circle) through the droplet surface above 8 Å (red) to exclude the near wall region (blue) according to

[27]. The green line represents the tangent to the fitted droplet surface in order to determine the contact angle.

https://doi.org/10.1371/journal.pone.0064883.g002

After energy minimization all systems were equilibrated for 100 ps at a temperature of 310 K and constant pressure of 1 atm with the protein backbone atoms restrained in order to prevent adsorption before equilibration. Following the equilibration classical and accelerated MD simulations were carried out for 20 ns with a time step of 2.0 fs enabled by the SHAKE [29] algorithm to ensure rigid hydrogen atoms. Furthermore, orientation 2 was also subjected to a 100 ns simulation run using classical MD.

In order to sufficiently sample the conformational space of adsorption the implementation of accelerated molecular dynamics (aMD) [10] into NAMD [30] was used. In short, a boost potential is added to the real potential when this is below a chosen value [10]:(1)with the boost potential:

(2)Here a dual acceleration was carried out where the boost potential is applied both to the dihedral potential and to the total potential. This means that besides the protein conformational motion also the diffusivity is increased as the motion of the solvent molecules is also accelerated [31], [32]. The boost parameters were set based upon the average dihedral and potential energy calculated by performing a short classical molecular dynamics (cMD) simulation. They were set according to the recipe given by [31]:(3)(4)(5)(6)

We also let these simulations run for time steps, as done in the cMD simulations. As mentioned in the Introduction, the assignment of a real time corresponding to aMD time steps is nontrivial; for brevity we shall in the following also call these aMD simulations to run up to ‘20 ns’.

Analysis of the Results

The mass weighted radius of gyration is defined as:(7)with beeing the distance of atom with mass to the center of mass of the protein. In this work we focus on the component parallel to the surface, i. e., , which gives an estimate of the protein spreading on the surface.

Furthermore, principal-component analysis using the cartesian coordinates of the protein atoms was performed with the GROMACS [33] analysis tools in order to describe the conformational sampling of the protein adsorption.

The adsorption snapshots were rendered using VMD [34] and Tachyon [35].

Results and Discussion

During the first 100 ps both classical and accelerated simulations show a similar trend. The protein approaches the graphite surface driven by attractive van-der-Waals interactions. But very soon classical and accelerated simulations show quite a different adsorption behavior. In the case of aMD, conformational transitions start to occur instantly whereas cMD simulations show very little structural rearrangements after the protein got in contact with the surface. During the accelerated simulations not only the protein motions but also the solvent motions were affected. Since dehydration of the surface is an essential but time consuming step in the adsorption process [5], aMD is able to accelerate the simulated time scale by a considerable factor.

The final adsorption snapshots after 20 ns, which compare cMD and aMD simulations (Fig. 3), clearly show that accelerated MD has a great impact on the overall structure. A distinctive feature which can be observed is the parallel orientation to the graphite surface of most of the aromatic amino acids (His, Phe, Tyr, Tryp). This finding is in agreement with theoretical work by other authors who investigated amino acid surface interactions [36], [37].

[Figure omitted. See PDF.]

Figure 3. Snapshot of the BMP-2 molecule (orientation 1) after 20 ns. a) cMD; b) aMD with sideview (left) and topview (right).

https://doi.org/10.1371/journal.pone.0064883.g003

Fig. 4 sheds additional light on the improved performance of aMD by displaying the temporal evolution of the van-der-Waals interaction energies between the protein and the graphite surface. Since the surface is neither charged nor polarized in this model we can conclude that van-der-Waals interactions are governing the adsorption process. Note first that indeed the cMD simulation only binds the protein in a shallow metastable precursor state, from which it cannot escape in the 20 ns run; aMD, in contrast, allows for a steady decrease of the binding energy, as the protein continuously changes conformation and acquires deeper bound states. This figure also demonstrates that the number of time steps simulated in this study, appears sufficient to estimate the final adsorption energy to be around 800 kcal/mol. This is a factor of larger than the cMD simulations predict. Note finally that orientation 1 achieves a better binding (larger adsorption energy) than orientation 1, see our discussion of Fig. 4 below.

[Figure omitted. See PDF.]

Figure 4. Attractive van-der-Waals energies between BMP-2 and the graphite surface.

https://doi.org/10.1371/journal.pone.0064883.g004

The protein, in the case of aMD, has spread completely on the graphite surface and forms a flat peptide monolayer. This can be confirmed by the temporal evolution of the radius of gyration (Fig. 5). Here, orientation 2 has spread more, both for classical and accelerated MD. The reason hereto is that the molecule has more flexibility in this orientation for spreading out on the surface due to weaker van-der-Waals interactions formed in the early adsorption stage until time steps (see Fig. 4); this feature has already been discussed in our previous work [19].

[Figure omitted. See PDF.]

Figure 5. Evolution of the radius of gyration of BMP-2 (component parallel to the graphite surface).

https://doi.org/10.1371/journal.pone.0064883.g005

Interestingly, both orientations in the case of aMD show the same secondary structure content after adsorption as can be seen in Table 1. It seems that after 20 ns of accelerated MD simulation using the above mentioned boost parameters the final unfolded adsorption state has been found.

[Figure omitted. See PDF.]

Table 1. Secondary structure of BMP-2 before adsorption and at the end of the simulation calculated with dssp [43].

https://doi.org/10.1371/journal.pone.0064883.t001

When comparing aMD and cMD one can clearly recognize that the helical and -sheet content is only slightly reduced in classical simulations. Even during the 100 ns simulation there has not been a great reduction in the -sheet content while one can recognize fluctuations in the helical content comparing with the 20 ns cMD simulation of orientation 2. On the other hand, accelerated simulations show almost no remaining secondary structure except one small -helix. This finding of denaturation is in agreement with other adsorption studies on hydrophobic graphite surfaces [6], [7]. As can be seen in Fig. 3 b topologically distant protein strands show a roughly parallel arrangement which is believed to result from graphite's hydrophobicity, crystallinity, and smoothness [38] and the optimized intramolecular interactions [7]. This phenomenon of unfolded proteins with parallel strands was observed experimentally [3]; this agreement demonstrates that aMD provides reasonable adsorption structures.

Additionally, we performed cartesian principal component analysis to identify the conformational space during protein adsorption (Fig. 6). In this case, the trajectory was projected onto the first two eigenvectors of the protein atom covariance matrix which account for the largest internal motions of the protein. From these results we can clearly conclude that aMD provides considerably improved sampling of the conformational space. This result is in agreement with the findings of investigations using aMD in other applications [10], [12], [30]. Especially in the case of orientation 2 very large conformational motions could be identified. Furthermore, for aMD a clear unfolding pathway can be seen for both orientations where many conformational states were visited. Extending the classical adsorption simulation to 100 ns resulted only in a slight extension of the visited conformational space. Other results employing long classical simulations (300 ns) suggest that even after long simulation times, which seemingly equilibrated the system, sudden structural changes may occur [5].

[Figure omitted. See PDF.]

Figure 6. Conformational space obtained by projection of the trajectory onto the first two eigenvectors of the protein atom covariance matrix.

Each point represents a conformational state at a simulation time shown by the color code.

https://doi.org/10.1371/journal.pone.0064883.g006

Conclusions

Comparing our results obtained by classical and accelerated molecular dynamics simulations we conclude that sufficient sampling of the conformational space is key to understanding protein adsorption. Since adsorption phenomena take place in time scales extending from s to hours [4], the use of classical MD simulations of a few nanoseconds up to 100 ns in an explicit water surrounding is not adequate for studying and describing the entire adsorption process. Our approach uses an explicit water environment with standard temperature and pressure control. By modifying the potential energy landscape through accelerated MD, conformational changes can occur in the time range accessible to simulation. Especially during adsorption the protein can get stuck in metastable states that hinder further spreading on hydrophobic surfaces.

The advantages of this method are – among many others (see [10]) – (i) the potential for application in many different surfaces ranging from hydrophobic to hydrophilic, (ii) the realistic system setup (e. g. inclusion of ions, atmospheric pressure) with no need for very long energy minimizations, (iii) a faster relaxation of water molecules [32], and finally (iv) only minor additional computational effort (see Figure 1 in [30]). Beside many available advanced sampling methods, recently reviewed by Zuckerman [39], we would like to point out an approach by Wang et al. [40] which uses biased replica-exchange molecular dynamics (REMD) to study adsorption processes. Comparing biased-REMD to accelerated MD we emphasize that aMD, besides being easy to handle, is computationally more efficient since no different replicas have to be simulated; the latter feature limits biased-REMD to smaller systems. On the other hand, REMD will give exact dynamics with a known reweighting factor.

The easiness of use and computational efficiency make aMD an ideal candidate for studying adsorption processes in large systems. Other interesting sampling techniques such as conformational flooding [41], [42] might be used as well to address adsorption studies but appear not so straightforward as aMD. Conformational flooding for instance needs calculation of a ‘flooding matrix’ and possibly multiple flooding potentials. Finally, coarse-grained models that might provide an interesting alternative for calculating large systems over long time scales always bear the risk of producing artifacts.

Our results suggest that the monomeric form of BMP-2 does indeed unfold on a hydrophobic graphite surface as suggested by previous results [19]. Despite the different protein structure (monomer vs. bioactive dimer [24]) this is in clear contrast to recent results by Utesch et al. [18]. It appears logical that a classical 10 ns MD simulation, such as that performed in [18], cannot provide insight into structural rearrangements during adsorption since the conformational space is sampled only insufficiently. We note that even in cMD simulations extending over several hundred ns, it has been found that the final configuration is governed by sudden conformation jumps which may occur late in the simulation [5].

Even though the dynamics of the adsorption process have been altered, the final adsorption state is most relevant in adsorption studies, since it can be compared with experimental results. From what we know from experimental results [3], e.g. an unfolded protein structure or the formation of parallel protein strands, the adsorption state obtained by aMD provides a considerably more realistic prediction than the ‘unfinished’ adsorption state after 20 ns of classical simulation. The next step in characterizing protein adsorption processes will be the utilization of aMD for studying a variety of different protein-surface interactions including initial orientation and topology dependences. Comparison to experiment is for instance possible by an analysis of the secondary structure content of the adsorbed protein or by determining the forces necessary for protein desorption.

Acknowledgments

We appreciate the computational resources provided by the compute cluster ‘Elwetritsch’ of the University of Kaiserslautern.

Author Contributions

Conceived and designed the experiments: CM HMU. Performed the experiments: CM. Analyzed the data: CM. Wrote the paper: CM HMU.

Citation: Mücksch C, Urbassek HM (2013) Enhancing Protein Adsorption Simulations by Using Accelerated Molecular Dynamics. PLoS ONE 8(6): e64883. https://doi.org/10.1371/journal.pone.0064883

References

1. Ratner BD, Bryant SJ (2004) Biomaterials: Where we have been and where we are going. Annu Rev Biomed Eng 6: 41.

2. Nakanishi K, Sakiyama T, Imamura K (2001) On the adsorption of proteins on solid surfaces, a common but very complicated phenomenon. Journal of Bioscience and Bioengineering 91: 233.

3. Svaldo-Lanero T, Penco A, PratoM, Canepa M, Rolandi R, et al. (2008) Nanopatterning by protein unfolding. Soft Matter 4: 965–967.

4. Gray JJ (2004) The interaction of proteins with solid surfaces. Current Opinion in Structural Biology 14: 110–115.

5. Wei T, Carignano MA, Szleifer I (2011) Lysozyme adsorption on polyethylene surfaces: Why are long simulations needed? Langmuir 27: 12074–12081.

6. Raffaini G, Ganazzoli F (2003) Simulation study of the interaction of some albumin subdomains with a flat graphite surface. Langmuir 19: 3403.

7. Raffaini G, Ganazzoli F (2010) Protein adsorption on a hydrophbic surface: A molecular dynamics study of lysozyme on graphite. Langmuir 26: 5679–5689.

8. Sun Y, Dominy BN, Latour RA (2007) Comparison of solvation-effect methods for the simulation of peptide interactions with a hydrophobic surface. Journal of Computational Chemistry 28: 1883–1892.

9. Voter AF (1997) Hyperdynamics: Accelerated molecular dynamics of infrequent events. Phys Rev Lett 78: 3908–3911.

10. Hamelberg D, Mongan J, McCammon JA (2004) Accelerated molecular dynamics: A promising and efficient simulation method for biomolecules. The Journal of Chemical Physics 120: 11919–11929.

11. Wang Y, Markwick PRL, de Oliveira CAF, McCammon JA (2011) Enhanced lipid diffusion and mixing in accelerated molecular dynamics. Journal of Chemical Theory and Computation 7: 3199–3207.

12. Pierce LC, Salomon-Ferrer R, Augusto F de Oliveira C, McCammon JA, Walker RC (2012) Routine access to millisecond time scale events with accelerated molecular dynamics. Journal of Chemical Theory and Computation 8: 2997–3002.

13. Xin Y, Doshi U, Hamelberg D (2010) Examining the limits of time reweighting and kramers’ rate theory to obtain correct kinetics from accelerated molecular dynamics. The Journal of Chemical Physics 132: 224101.

14. Jennissen HP (2002) Accelerated and improved osteointegration of implants biocoated with bone morphogenetic protein 2 (bmp–2). Ann NY Acad Sci 961: 139.

15. Sachse A, Wagner A, Keller M, Wagner O, Wetzel WD, et al. (2005) Osteointegration of hydroxyapatite-titanium implants coated with nonglycosylated recombinant human bone morphogenetic protein-2 (bmp-2) in aged sheep. Bone 37: 699–710.

16. Dong X, Wang Q, Wu T, Pan H (2007) Understanding adsorption-desorption dynamics of bmp-2 on hydroxyapatite (001) surface. Biophysical Journal 93: 750–759.

17. Oliveira AF, Gemming S, Seifert G (2010) Molecular dynamics simulations of bmp-2 adsorption on a hydrophobic surface. molekulardynamische simulation des bmp-2-adsorption auf einer hydrophoben oberfl¨ache. Materialwissenschaft und Werkstofftechnik 41: 1048–1053.

18. Utesch T, Daminelli G, Mroginski MA (2011) Molecular dynamics simulations of the adsorption of bone morphogenetic protein-2 on surfaces with medical relevance. Langmuir 27: 13144–13153.

19. Mücksch C, Urbassek HM (2011) Adsorption of bmp-2 on a hydrophobic graphite surface: A molecular dynamics study. Chemical Physics Letters 510: 252–256.

20. Winter M, Chuinard C, Cikes A, Pelegri C, Bronsard N, et al. (2009) Surgical management of elbow dislocation associated with non–reparable fractures of the radial head. Chirurgie de la Main 28: 158–167.

21. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, et al. (2005) Scalable molecular dynamics with namd. Journal of Computational Chemistry 26: 1781–1802.

22. MacKerell AD, Bashford D, Bellott, Dunbrack RL, Evanseck JD, et al. (1998) All-atom empirical potential for molecular modeling and dynamics studies of proteins. The Journal of Physical Chemistry B 102: 3586–3616.

23. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML (1983) Comparison of simple potential functions for simulating liquid water. The Journal of Chemical Physics 79: 926–935.

24. Scheufler C, Sebald W, Hülsmeyer M (1999) Crystal structure of human bone morphogenetic protein–2 at 2.7 å resolution. Journal of Molecular Biology 287: 103–115.

25. Darden T, York D, Pedersen L (1993) Particle mesh ewald: An n · log(n) method for ewald sums in large systems. The Journal of Chemical Physics 98: 10089–10092.

26. Zuo A, Liang D, Sun P, Zhang R, Zhao R, et al.. (2008) Bmp-2 modified hydroxyapatite/chitosan composite scaffold for bone engineering. The 2nd International Conference on Bioinformatics and Biomedical Engineering : 888–892.

27. Werder T, Walther JH, Jaffe RL, Halicioglu T, Koumoutsakos P (2002) On the water–carbon interaction for use in molecular dynamics simulations of graphite and carbon nanotubes. Journal of the American Chemical Society 107: 1345.

28. Fowkes FM, Harkins WD (1940) The state of monolayers adsorbed at the interface solid–aqueous solution. Journal of the American Chemical Society 62: 3377–3386.

29. Ryckaert JP, Ciccotti G, Berendsen HJ (1977) Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. Journal of Computational Physics 23: 327–341.

30. Wang Y, Harrison CB, Schulten K, McCammon JA (2011) Implementation of accelerated molecular dynamics in namd. Comput Sci Discov 4: 015002.

31. Wereszczynski J, McCammon J (2012) Computational Drug Discovery and Design, volume 819 of Methods in Molecular Biology. Springer New York, 515–524 pp.

32. de Oliveira CAF, Hamelberg D, McCammon JA (2006) On the application of accelerated molecular dynamics to liquid water simulations. The Journal of Physical Chemistry B 110: 22695–22701.

33. Hess B, Kutzner C, van der Spoel D, Lindahl E (2008) Gromacs 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. Journal of Chemical Theory and Computation 4: 435–447.

34. Humphrey W, Dalke A, Schulten K (1996) VMD – Visual Molecular Dynamics. Journal of Molecular Graphics 14: 33–38.

35. Stone J (1998) An Efficient Library for Parallel Ray Tracing and Animation. Master’s thesis, Computer Science Department, University of Missouri-Rolla.

36. Grimme S (2008) Gibt es spezielle nicht–kovalente π–π–stapelwechselwirkungen wirklich? Angewandte Chemie 120: 3478.

37. Rajesh C, Majumder C, Mizuseki H, Kawazoe Y (2009) A theoretical study on the interaction of aromatic amino acids with graphene and single walled carbon nanotube. Journal of Chemical Physics 130: 124911–1.

38. Raffaini G, Ganazzoli F (2004) Surface ordering of proteins adsorbed on graphite. The Journal of Physical Chemistry B 108: 13850–13854.

39. Zuckerman DM (2011) Equilibrium sampling in biomolecular simulations. Annual Review of Biophysics 40: 41–62.

40. Wang F, Stuart S, Latour RA (2008) Calculation of adsorption free energy for solute-surface interactions using biased replica-exchange molecular dynamics. Biointerphases 3: 9–18.

41. Grubmüller H (1995) Predicting slow structural transitions in macromolecular systems: Conformational flooding. Phys Rev E 52: 2893–2906.

42. Lange OF, Schäfer LV, Grubmüller H (2006) Flooding in gromacs: Accelerated barrier crossings in molecular dynamics. Journal of Computational Chemistry 27: 1693–1702.

43. Kabsch W, Sander C (1983) Dictionary of protein secondary structure: pattern recognition of hydrogen-bonded and geometrical features. Biopolymers 22: 2577.