Molecular Dynamics Simulations of PeptideâSurface Interactions

Langmuir 2005, 21, 1629-1639

1629

Molecular Dynamics Simulations of Peptide-Surface Interactions Vivek P. Raut,† Madhuri A. Agashe,† Steven J. Stuart,‡ and Robert A. Latour*,† Department of Bioengineering and Department of Chemistry, Clemson University, Clemson, South Carolina 29634 Received September 1, 2004. In Final Form: November 14, 2004 Proteins, which are bioactive molecules, adsorb on implants placed in the body through complex and poorly understood mechanisms and directly influence biocompatibility. Molecular dynamics modeling using empirical force fields provides one of the most direct methods of theoretically analyzing the behavior of complex molecular systems and is well-suited for the simulation of protein adsorption behavior. To accurately simulate protein adsorption behavior, a force field must correctly represent the thermodynamic driving forces that govern peptide residue-surface interactions. However, since existing force fields were developed without specific consideration of protein-surface interactions, they may not accurately represent this type of molecular behavior. To address this concern, we developed a host-guest peptide adsorption model in the form of a G4-X-G4 peptide (G is glycine, X is a variable residue) to enable determination of the contributions to adsorption free energy of different X residues when adsorbed to functionalized Aualkanethiol self-assembled monolayers (SAMs). We have previously reported experimental results using surface plasmon resonance (SPR) spectroscopy to measure the free energy of peptide adsorption for this peptide model with X ) G and K (lysine) on OH and COOH functionalized SAMs. The objectives of the present research were the development and assessment of methods to calculate adsorption free energy using molecular dynamics simulations with the GROMACS force field for these same peptide adsorption systems, with an oligoethylene oxide (OEG) functionalized SAM surface also being considered. By comparing simulation results to the experimental results, the accuracy of the selected force field to represent the behavior of these molecular systems can be evaluated. From our simulations, the G4-G-G4 and G4-K-G4 peptides showed minimal to no adsorption to the OH SAM surfaces and the G4-K-G4 showed strong adsorption to the COOH SAM surface, which is in agreement with our SPR experiments. Contrary to our experimental results, however, the simulations predicted a relatively strong adsorption of G4-G-G4 peptide to the COOH SAM surface. In addition, both peptides were unexpectedly predicted to adsorb to the OEG surface. These findings demonstrate the need for GROMACS force field parameters to be rebalanced for the simulation of peptide adsorption behavior on SAM surfaces. The developed methods provide a direct means of assessing, modifying, and validating force field performance for the simulation of peptide and protein adsorption to surfaces, without which little confidence can be placed in the simulation results that are generated with these types of systems.

1. Introduction When a medical implant is placed in the body, surrounding proteins rapidly adsorb to the implant surface at the solid-liquid interface. Cells present in the body fluid, which then come in contact with the implant, do not interact with the actual implant material itself. Rather, they sense the implant by way of this adsorbed protein layer. Together with other interfacial processes, protein adsorption thus regulates the biological response to any implant material in contact with a biological fluid, and hence implant biocompatibility. This phenomenon is attributed to the specific biochemical signaling potential of a protein.1,2 For example, the adsorption of blood proteins, such as fibrinogen, can influence the adhesion of platelets or macrophages and ultimately lead to thrombus formation or fibrous encapsulation.2-5 Thus the study of protein-surface interactions is extremely im* Corresponding author. Address: Department of Bioengineering, 501 Rhodes Engineering Research Center, Clemson University, Clemson, SC 29634. E-mail: [email protected]. Phone: 864656-5552. Fax: 864-656-4466. † Department of Bioengineering. ‡ Department of Chemistry. (1) Anderson, J. M.; Bonfield, T. L.; Ziats, N. Int. J. Artif. Organs 1990, 13, 375-382. (2) Ratner, B. D.; Hoffman, A. S.; Schoen, F. J.; Lemons, J. E. In Biomaterials Science: An Introduction to Materials in Medicine; Academic Press: San Diego, 1996.

portant in the field of biomaterials and in many other areas of biotechnology.5-15 Because of the relevance of protein-surface interactions, much effort has gone into protein adsorption experiments and models over the past several decades14,16-18 with an ultimate aim to quantitatively (3) Jahangir, A. R.; McClung, W. G.; Cornelius, R. M.; McCloskey, C. B.; Brash, J. L.; Santerre, J. P. J. Biomed. Mater. Res. 2002, 60, 135-147. (4) Shen, M.; Martinson, L.; Wagner, M. S.; Castner, D. G.; Ratner, B. D.; Horbett, T. A. J. Biomater. Sci., Polym. Ed. 2002, 13, 367-390. (5) Ratner, B. D. J. Biomed. Mater. Res. 1993, 27, 837-850. (6) Texter, J.; Tirrell, M. AIChE J. 2001, 47, 1706-1710. (7) Nakanishi, K.; Sakiyama, T.; Imamura, K. J. Biosci. Bioeng. 2001, 91, 233-244. (8) Mulzer, S. R.; Brash, J. L. J. Biomed. Mater. Res. 1989, 23, 14831504. (9) Baurmeister, U.; Vienken, J.; Grassmann, A. Nephrol. Dial. Transplant. (Suppl.) 1991, 3, 17-21. (10) Linnola, R. J.; Werner, L.; Pandey, S. K.; Escobar-Gomez, M.; Znoiko, S. L.; Apple, D. J. J. Cataract Refractive Surg. 2000, 26, 17921806. (11) Widmer, M. R.; Heuberger, M.; Voros, J.; Spencer, N. D. Tribol. Lett. 2000, 10, 111-116. (12) Gura, T. A.; Wright, K. L.; Veis, A.; Webb, C. L. J. Biomed. Mater. Res. 1997, 35, 483-495. (13) Kasemo, B. Surf. Sci. 2002, 500, 656-677. (14) Hlady, V.; Buijs, J. Curr. Opin. Biotechnol. 1996, 7, 72-77. (15) Lange, R.; Bergbauer, M.; Szewzyk, U.; Reitner, J. Facies 2001, 45, 195-202. (16) Brash, J. L.; Horbett, T. A. In Proteins At Interfaces II: Fundamentals and Applications; American Chemical Society: Washington, DC, 1995; pp 1-23.

10.1021/la047807f CCC: $30.25 © 2005 American Chemical Society Published on Web 01/14/2005

1630

Langmuir, Vol. 21, No. 4, 2005

measure, predict, and understand the details of proteinsurface interactions. As described by Norde et al.19 and Sigal et al.,20 proteins typically adsorb strongly to hydrophobic surfaces and weakly to neutral hydrophilic surfaces. Charged surfaces are generally found to be more adsorptive for oppositely charged proteins, and the degree of adsorption is typically lower for similarly charged surfaces and proteins.1,19,20 While these trends are easily conceptualized, the numerous simultaneous interactions occurring between the functional groups of proteins, material surfaces, and surrounding body fluids are very complex in nature and the actual submolecular-level mechanisms and structural rearrangements involved in these reactions are not well understood. As a result, the postadsorption state of a protein (its conformation, orientation, and bioactivity) cannot currently be accurately predicted. Partly as a result of this problem, it has been widely accepted that nonspecific protein adsorption should be prevented altogether and peptide ligands, such as arginine-glycine-aspartic acid (RGD), should instead be used in an attempt to control cellular response at the biomaterials-biofluid interface. While we agree that uncontrolled protein adsorption should be avoided, we believe that great potential still exists for the use of surface chemistry to control adsorbed protein orientation and conformation and to then use the natural bioactivity of proteins to direct cellular response. This type of specific surface design, however, will be possible only if submolecular-level protein-surface interactions can be quantitatively understood and predicted. Molecular simulation provides one of the most direct methods to theoretically investigate molecular behavior in complex systems, such as in a system involving protein adsorption. Because of the size of the molecular systems involved, empirical force field based molecular modeling methods, such as Monte Carlo (MC) and molecular dynamics (MD), are required for the simulation of these types of processes.21 These methods employ a potential energy function (referred to as a force field) that calculates the overall potential energy of the system based on the summation of individual atom-atom pair interactions. The force field equation takes into account the contributions due to bonded interactions (bond stretching, bending, and torsion) as well as nonbonded interactions (van der Waals and electrostatic). These energy contributions are determined by a set of empirical parameters that are used by the force field to calculate the energy contribution for each type of interaction for the atom types that are defined in the molecular system being considered.21 To accurately simulate molecular behavior using an empirical force field, parameters of the force field must be balanced and tuned together to appropriately represent the behavior of the given molecular system being addressed. As a consequence of this, a force field designed for one type of application cannot be confidently applied to other applications without separate validation; this issue is known as force field transferability.21 For example, a problem of transferability was shown by Schuler et al.22 with the GROMOS96 force field, which was initially (17) Sadana, A. Chem. Rev. 1992, 92, 1799-1818. (18) Claesson, P. M.; Blomberg, E.; Froberg, J. C.; Nylander, T.; Arnebrant, T. Adv. Colloid Interface Sci. 1995, 57, 161-227. (19) Norde, W. In Adhesion and Adsorption of Polymers, Vol. 2; Plenum Press: New York, 1980; pp 801-826. (20) Sigal, G. B.; Mrksich, M.; Whitesides, G. M. J. Am. Chem. Soc. 1998, 120, 3464-3473. (21) Leach, A. R. In Molecular modeling: Principles and applications; Addison-Wesley Longman: Essex, U.K., 1996. (22) Schuler, L. D.; Daura, X.; van Gunsteren, W. F. J. Comput. Chem. 2001, 22, 1205-1218.

Raut et al.

parametrized to represent the behavior of proteins in aqueous solutions. When GROMOS96 was applied to represent the behavior of condensed phase hydrocarbons, it was found to result in substantial errors and required new parameters to be developed for this particular type of molecular system even though prior parametrization included similar types of alkane functional groups, but in a different molecular setting. Over the past few decades, several force fields have been developed and validated for the simulation of proteins in an aqueous solution (e.g., AMBER,23 CHARMM,24 GROMACS,25,26 and OPLS27). However, none of these force fields have been specifically parametrized with consideration of the adsorption behavior of proteins to synthetic material surfaces. Protein-surface interactions are dominated by nonbonded interactions between the functional groups presented by amino acids, various functional groups presented by material surfaces (e.g., polymers), and water. For a force field to accurately represent this type of system, it must be parametrized to properly balance all of these types of interactions with one another. As a further complication, however, when it comes to peptidesurface adsorption behavior, very little is known about the actual molecular behavior that a properly parametrized force field should reproduce. Therefore, before a given empirical force field can be used confidently in an MC or MD simulation to accurately represent the adsorption behavior of a protein to a surface, experimental methods must first be developed to allow validation of the simulations. These experiments should probe a fundamental characteristic of the adsorption behavior between specific peptide residues and surface functional groups for a system that can also be readily represented by molecular simulation. Simulation results with a given force field can then be compared with the experimental data in order to evaluate and validate or modify force fields for specific molecular systems. To directly address these needs, Vernekar and Latour previously developed an experimental method using surface plasmon resonance (SPR) spectroscopy28 to measure the adsorption free energy of individual peptide residues on functionalized alkanethiol self-assembled monolayer (SAM) surfaces on gold using a host-guest peptide system of the form GGGG-X-GGGG (G ) glycine, X ) variable residue). These methods were designed so that the effects of single midchain residue substitutions on the adsorption behavior of this small peptide on functionalized surfaces could be measured using a system that was readily amenable to MD simulation. These methods were applied to determine the free energy of adsorption (∆Gads) for the peptides G4-G-G4 and G4-K-G4 (G ) glycine, K ) lysine) on hydroxyl (OH) and carboxyl (COOH) functionalized SAM surfaces over a range of temperatures from 288 to 310 K. From these studies, no peptide adsorption (i.e., ∆Gads ) 0) was observed for either (23) Pearlman, D. A.; Case, D. A.; Caldwell, J. W.; Ross, W. R.; Cheatham, T. E.; DeBolt, S.; Ferguson, D.; Seibel, G.; Kollman, P. Comput. Phys. Commun. 1995, 91, 1-41. (24) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (25) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Commun. 1995, 91, 43-56. (26) van der Spoel, D.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Hess, B.; Feenstra, K. A.; Lindahl, E.; van Drunen, R.; Berendsen, H. J. C. In GROMACS User Manual, version 3.1; University of Groningen: Groningen, The Netherlands, 2003. (27) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. J. Phys. Chem. B 2001, 105, 6474-6487. (28) Vernekar, V. N.; Latour, R. A., Jr. Mater. Res. Innovations, in press.

Simulations of Peptide-Surface Interactions

Langmuir, Vol. 21, No. 4, 2005 1631

Figure 1. Molecular model of host-guest peptide design in the form of GGGG-X-GGGG. Two types of guest residues were evaluated: X ) lysine and glycine.

peptide on the OH SAM surface or for the G4-G-G4 peptide on the COOH SAM surface, while a strong adsorption response, e.g., ∆Gads ) -7.0 ( 0.1 kcal/mol (mean ( std deviation) at 300 K, was measured for the G4-K-G4 peptide on the COOH SAM surface. The objective of our present research was therefore to develop and assess MD simulation methods that could be used to calculate the adsorption free energy for these hostguest peptide-SAM systems for comparison with these experimental results, using GROMACS as an exemplary protein simulation force field. Peptide adsorption to an oligoethylene oxide (OEG) functionalized SAM surface was also simulated due to the expected nonadsorption behavior for this type of surface.29,30 The results of these simulations showed distinct differences in peptide adsorption behavior as a function of the guest peptide and surface functional group combinations, with the predicted trends being in general agreement with the experimental results of Vernekar and Latour.28 Several areas were also identified, however, where the simulation results show what is considered to be unrealistic adsorption behavior, which demonstrates the need for force field parameter modification if GROMACS is to be used to accurately simulate peptide or protein adsorption to these types of surfaces. 2. Materials 2.1. Computational Environment. The initial peptide models as well as the functionalized SAM surface models were designed using the InsightII [Accelrys Inc., San Diego, CA] and Web Lab Viewer Pro (version 3.20) software [Accelrys Inc.]. Molecular dynamics simulations were performed on a 26-node dual-processor 2.66 GHz Linux cluster with GROMACS software (version 3.1.3).26 Each simulation was run on a dedicated node, with a calculation time of about 1.5 days of CPU time per nanosecond of simulation. 2.2. Peptide Models. The peptide residues were designed as zwitterionic host-guest peptide models with a peptide sequence of GGGG-X-GGGG, where G ) glycine and X represents any residue type. This peptide design enables the middle peptide residue X (the guest residue) to impart specific characteristics to the peptide. Two types of guest peptide residues were modeled (29) Kreuzer, H. J.; Wang, R. L. C.; Grunze, M. J. Am. Chem. Soc. 2003, 125, 8384-8389. (30) Wang, R. L. C.; Kreuzer, H. J.; Grunze, M.; Pertisn, A. J. Phys. Chem. Chem. Phys. 2002, 2, 1721-1727.

Figure 2. Molecular structures of the functionalized alkane chains used to construct alkanethiol SAM surfaces: (A) OH functionalization, (B) COOH functionalization, (C) COOfunctionalization, (D) trans OEG functionalization, and (E) helical OEG functionalization (torsional angles set as follows: C-O-C-C ) 70°, C-C-O-C ) 70°, O-C-C-O ) 60° (refs 30 and 32)). to compare with our prior experimental studies: X ) lysine (K) and glycine (G). The structures of these peptides are shown in Figure 1. The glycine residue, with its simple -H side-group, is a relatively polar residue with a hydropathy index of -0.4,31 while the side-group of lysine {-[CH2]4-NH3+} is positively charged at a pH of 7.4. 2.3. SAM Surfaces. Alkanethiol SAM surfaces on gold [AuS(CH2)n-X] with X ) OH, COOH, and (O-CH2-CH2)2-OH) (OEG) functionalities were selected for the model surfaces for this study. The OH and COOH SAMs were selected to correspond to the surfaces used in the prior experimental studies by Vernekar and Latour.28 The OEG surface was used because of its expected nonadsorptive behavior for peptides and proteins.29,30 Two types of OEG structures were evaluated: an all-trans structure, which represents the minimum energy conformation of the OEG functional group in vacuum or air32 and a helical conformation, which is believed to represent the minimum energy structure of OEG in aqueous solution.29,30 Figure 2 shows molecular models of the structure of each type of alkanethiol chain used to construct the SAM surfaces. The lower portion of each chain has been truncated to provide a surface of approximately 10 Å thickness; (31) Kyte, J.; Doolittle, R. F. J. Mol. Biol. 1982, 157, 105-132. (32) Rigby, D.; Sun, H.; Eichinger, B. E. Polym. Int. 1997, 44, 311330.

1632


Raut et al.

Figure 3. Molecular models of alkanethiol SAM surfaces: (A) hydroxyl (-OH) terminated SAM, (B) COO-/COOH terminated SAM, (C) all-trans OEG terminated SAM, and (D) helical OEG terminated SAM. This structure is produced by repeating the oriented hydroxyl terminated alkanethiol chain in a hexagonal close-packed (x3 × x3)R30° array to form a surface of required dimensions. Table 1. GROMACS Atoms Types Employed to Define New Surface Topologies atom name

description

mass (amu)

O OM OA C C1 C2 C3 HO

carbonyl oxygen (CdO) carboxyl oxygen (CO- and C-O-C) hydroxyl oxygen (OH) bare carbon (CdO, C-N) aliphatic CH group aliphatic CH2 group aliphatic CH3 group hydroxyl hydrogen

15.9994 15.9994 15.9994 12.0110 12.0110 12.0110 12.0110 1.0080

this represents the SAM structure effectively while reducing the size of the system to be simulated. SAM surface models of each type were developed by replicating a single truncated alkanethiol chain in a hexagonal close-packed (x3 × x3)R30° array with 4.97 Å lattice spacing to form a 36 Å × 40 Å surface plane, with chain orientation set according to the well-established alkanethiol chain structure on a gold (111) plane33 that has been used in several of our prior SAM models.34-38 The COOH SAM surface was constructed with both COOH (protonated) and COO- (deprotonated) functional groups. One out of twenty (5%) of the COOH groups were deprotonated and evenly positioned throughout the SAM surface, as appropriate for a surface pKa of about 8.7 in a solution with pH of 7.4. This surface pKa was obtained from work by Creager and Clarke.39 Figure 3 shows molecular models of each of the four SAM surfaces: OH, COO-/COOH, trans OEG, and helical OEG. The GROMACS library defines the residue topologies of all the amino acids and some other atoms that are used to model larger proteins. However, GROMACS does not define the residue topologies of individual SAM chains. The topologies corresponding to the individual SAM chains therefore had to be defined in the residue topology parameter (.rtp) file of the GROMACS standard ffgmx2 force field.26 Table 1 summarizes the atom types used to define these molecules. For consistency with the GROMACS model, all aliphatic hydrogens were treated as united atoms, (33) Ulman, A.; Eilers, J. E.; Tillman, N. Langmuir 1989, 5, 11471152. (34) Latour, R. A.; Hench, L. L. Biomaterials 2002, 23, 4633-4648. (35) Latour, R. A. Curr. Opin. Solid State Mater. Sci. 1999, 3, 413417. (36) Basalyga, D. M.; Latour, R. A. J. Biomed. Mater. Res. 2003, 64A, 120-130. (37) Latour, R. A.; Rini, C. J. J. Biomed. Mater. Res. 2002, 60, 564577 (38) Wilson, K.; Stuart, S. J.; Garcia, A.; Latour, R. A. J. Biomed. Mater. Res. 2004, 69A, 686-698 (39) Creager, S. E.; Clarke, J. Langmuir 1994, 10, 3675-3683.

Table 2. Charge Distribution on the Terminal Surface Functionalities terminal group (X)

partial charges based on:

atom type

partial charge in |e|

charge group

OH

serine (SER)

C2 OA HO

0.150 -0.548 0.398

0 0 0

COO-

aspartic acid (unprotonated) (ASP)

C2 OM OM

0.270 -0.635 -0.635

0 0 0

COOH

aspartic acid protonated) (ASPH)

C2 O OA HO

0.530 -0.380 -0.548 0.398

0 0 0 0

C2 OM C2 C2 OM C2 C2 OA HO

0.200 -0.400 0.200 0.200 -0.400 0.200 0.263 -0.566 0.303

1 1 1 2 2 2 3 3 3

OEG

while the hydrogens associated with the functionalized surface groups were separately defined by modifying the hydrogen database (.hdb) file of the GROMACS standard ffgmx2 force field. The partial charges assigned to each atom of these SAM surface residues were based on residue topologies for the side-groups of amino acids defined in GROMACS with corresponding functional groups. However, the partial charges associated with the OEG functional groups do not have corresponding predefined functionality in GROMACS. In order maintain self-consistency with the GROMACS force field as much as possible, the partial charge distribution on the OEG SAM chain was therefore borrowed from the partial charges for this functional group from another class I protein force field (AMBER23) as determined using the InsightII software program [Accelrys Inc.]. Table 2 summarizes the partial charge distribution on the terminal functional groups of each type of SAM surface. As per GROMACS design rules, the united carbon atoms of the SAM surface (i.e., the CH2 groups of the alkane chains of the SAM surface that were not part of the OEG segments) were represented as having zero partial charge.

3. Methods 3.1. Assembly of Residue-Surface System. Before conducting the simulated adsorption studies, short MD


Figure 4. Typical molecular model constructed to depict a G4-K-G4 peptide over an OH SAM surface with explicitly represented water.

simulations were performed on each respective peptide in vacuum (300 K, 100 ps, 1 fs time step) in order to generate a random configuration. The randomly structured peptides were then placed over the alkanethiol SAM surface models with the distance between the CR atom of the guest peptide and the plane of the terminal groups of the SAM surface (defined as the surface separation distance, SSD) initially set to 7.5 Å. Three-dimensional periodic boundary conditions were defined, with the size of the periodic cell being 36 Å × 40 Å in the plane of the surface and 65 Å perpendicular to the surface. The system was then solvated with simple point charge (SPC) water molecules with Na+ and Cl- ions added to approximate a 0.150 M physiologic saline solution40 and to maintain overall system neutrality. Water molecules falling within a 2 Å radius from the center of heavy atoms (i.e., nonhydrogen atoms) were then deleted to establish the final molecular assembly, with the resulting models containing approximately 8900 atoms. Figure 4 illustrates the initial assembly of the complete peptide-surface system for a G4-K-G4 peptide over an OH SAM surface. 3.2. Simulation Specifications. Once created, each molecular system was then energy-minimized without constraints using the steepest descent integrator for 5000 steps with the initial step size of 0.1 Å (the minimization tolerance was set to 1000 kJ/(mol nm)). Molecular dynamics simulations were then performed on this energy-minimized system using the leap-frog algorithm in the NVT (canonical) ensemble. All bond lengths were constrained, which enabled a time step of 2 fs to be used for the MD simulations.26 These conditions were established based on preliminary studies that were conducted to compare the trajectories for constrained and unconstrained simulations with the time step set at 0.5, 1, and 2 fs. Similar peptide-surface interaction behavior was observed for each of these conditions, thus enabling bond length constraints and a 2 fs time step to be used confidently. Temperature coupling with a Berendsen thermostat was implemented with a time constant of 0.1 ps at a temperature of 300 K. No pressure coupling was (40) West, J. B. In Best and Taylor’s Physiological Basis of Medical Practice, 11th ed.; Williams and Wilkins: Baltimore, MD, 1985; pp 441-442.

Langmuir, Vol. 21, No. 4, 2005 1633

applied to the system. Initial velocities were generated according to the Maxwell distribution at 300 K. The cutoff distances used for both the Columbic and the van der Waals interactions were 17 Å. All simulations were run for 10 ns of simulated time, and the system configuration (atom coordinates, velocities, and energies) was saved every 1.0 ps. In all, eight different molecular systems were evaluated (G4-G-G4 and G4-K-G4 peptides; OH, COOH, trans OEG, and helical OEG SAM surfaces), with three independent simulations conducted for each system using different initial peptide orientations over the SAM surface and different initial velocities. The resulting trajectory files were then viewed and analyzed using VMD software.41 3.3. Data Analysis. The trajectories for each molecular system were analyzed to calculate the adsorption free energy. This was accomplished using the probability ratio method.42,43 For this analysis, the GROMACS tool “g-traj” was used to track SSD as a function of time, which was then transformed into a probability density distribution with respect to SSD. The probability density distributions were then analyzed to calculate the relative free energy as a function of SSD, from which the overall adsorption free energy for the system was calculated. The probability density distribution for each peptideSAM surface system was determined by first dividing the SSD range of the simulation into 0.2 Å intervals (∆xi). The frequency of a given peptide being in a particular interval of SSD (SSDi) was then calculated as the number of times (Fi) that the peptide was positioned in a given SSDi interval during the 10 ns simulation. The positional probability (Ai) of the peptide being in a given SSDi location was then calculated from these data as

Ai )

Fi

∑Fi

with

∑(Ai) ) 1

(1)

Ai was then divided by the width of the interval, ∆xi, to calculate the normalized probability density (Pi) for each SSDi, or

Pi )

Ai ∆xi

with

∑(Pi ∆xi) ) 1

(2)

The normalized probability density (NPD) plot over the SSD range for each MD simulation was then used to calculate the adsorption free energy as a function of SSD. Once the NPD values were determined, the relative free energy difference (∆Gi) between two SSDi positions was calculated as

[]

∆Gi ) Gi - G0 ) -RT ln

Pi P0

(3)

where R is the ideal gas constant constant, T is the absolute temperature, and Pi and P0 represent the probability density of the peptide being at SSDi and SSD0, respectively, where SSD0 represents a defined reference state for the given molecular system. Finally, the overall free energy of adsorption (∆Gads) for the system was calculated as the weighted sum of the relative free energies, or (41) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33-38. (42) Mezei, M. Mol. Simul. 1989, 3, 301-313. (43) Gunsteren, W. F. V.; Weiner, P. K.; Computer Simulation of Biomolecular Systems: Theoretical and Experimental Applications (Volume 1); ESCOM Science: Leiden, The Netherlands, 1989.

1634


∆Gads )

∫SSD Pi ∆Gi dxi ≈ ∑[Pi ∆Gi ∆xi] i

Raut et al.

(4)

Following the above-described methods, plots of SSD vs time, NPD vs SSD, and ∆Gi vs SSD were generated for each of three independent 10 ns simulations for each peptide-SAM surface system; the ∆Gi vs SSD plots were then integrated to calculate ∆Gads for each run, from which were obtained the mean and standard deviation of ∆Gads for each system. These values were then compared with the experimentally determined values of ∆Gads to assess the performance of the GROMACS force field to represent this type of molecular system. In addition, the trajectory data were graphically viewed to evaluate the specific functional group interactions predicted to occur between the peptides and the SAM surfaces that were responsible for the observed adsorption behavior. 4. Results 4.1. Peptide Adsorption on the OH SAM Surface. In the initial simulations with the OH SAM, the lower two carbon atoms of the SAM surface were fixed in order to maintain the orientation of each of its chains while enabling the rest of each chain to freely move in response to its environment. However, as shown in Figure 5, this resulted in the peptides adsorbing strongly to the surface, which should not occur based on our previous experimental results.28 Upon inspection of the functional group interactions responsible for this behavior, it was observed that the adsorption response was largely due to hydrophobic interactions between CH2 segments of the peptides and the CH2 groups of the SAM that were exposed when the OH groups on the SAM surface separated due to the motion of the SAM chains (see Figure 6). To prevent this behavior, the coordinate positions of all atoms of the SAM were then fixed except for the surface OH groups, which solved this problem, and the peptides no longer adsorbed to this surface (Figure 5). Subsequent simulations for all of the SAM surfaces were conducted in this same manner, with the atoms all fixed except for the top functional groups, which were free to respond to the surrounding atoms. This highlights one shortcoming of the GROMACS potential for SAM surfaces, in that the parameters for the united-atom CH2 groups result in overly mobile chains at the correct packing densities, which also generates incorrect surface heights and tilt angles without the imposition of additional constraints. Representative data plots for SSD vs time, the NPD distribution of SSD, and ∆Gi as a function of SSD are presented in Figure 7 for both the G4-G-G4 and G4-K-G4 peptides on the constrained OH SAM surface. Similar data analyses were applied for each of the three independent simulations for each of these peptides on the OH SAM surface, but only one exemplary set of data plots is shown for conciseness. In producing the NPD plots, it was necessary to discard the SSD vs time data for SSD > 35 Å because the use of the 3-D periodic boundary conditions caused the peptides to interact with the image of the bottom of the SAM surface (with its fixed CH3 functionality) when the peptide moved within the cutoff distance of the top of the simulation cell. To plot ∆Gi as a function of SSD, a reference state (SSD0) had to be defined to provide the value of P0 for the calculation of the free energy at each SSDi position relative to SSD0 (see eq 3). After inspecting the NPD vs SSD plots, it was decided to use the average value of Pi between the SSD of 25 and 35 Å for P0 because this represents a region in the simulation where the peptide is beyond the cutoff distance for nonbonded interactions with the SAM surface at both the

Figure 5. SSD vs time trajectory plot for the unconstrained OH SAM surface compared to the constrained OH SAM surface: (A) G4-G-G4 peptide and (B) G4-K-G4 peptide.

Figure 6. A snapshot from the trajectory file for the G4-K-G4 peptide on the unconstrained OH SAM surface showing the CH2 segment of a glycine residue (*) hydrophobically interacting with subsurface CH2 groups of the SAM that were exposed when the surface OH groups separated. Hydrogen bonds formed by the NH and CO of this glycine residue with the OH groups of the SAM surface also helped stabilize this interaction.

bottom and top of the simulation cell. This assumption is equivalent to setting the zero of the free energy at an SSD of 30 Å. Based on these analyses, the adsorption free energies (mean ( std deviation) for the G4-G-G4 and G4K-G4 peptides on the OH SAM surface were calculated to be 0.12 ( 0.08 and -0.09 ( 0.39 kcal/mol, respectively. As shown from the plots presented in Figure 7 and by the near zero value of adsorption free energy, the peptides did not adsorb on the OH SAM surfaces, indicating that water molecules formed more thermodynamically favorable interactions with the SAM surface and peptide functional groups than the SAM and the peptides did with each other. This is in agreement with the experimental studies of Vernekar and Latour, which also indicate that


Langmuir, Vol. 21, No. 4, 2005 1635

Figure 7. SSD vs time, NPD (Pi) vs SSD, and ∆Gi vs SSD for each of the peptides on an OH SAM: (A) G4-G-G4 peptide and (B) G4-K-G4 peptide.

neither the G4-G-G4 nor the G4-K-G4 peptide adsorbs to an OH SAM.28 From viewing the trajectory files, it was observed that each of these peptides translated away from the surface in an apparently random manner after a few initial interactions with the SAM surface in the form of temporary hydrogen bonds formed between the various charged and polar functional groups of the peptide and the OH functional groups of the SAM. These interactions, however, were readily displaced by the surrounding SPC water molecules, which formed more stable hydrogen bonds with the OH functional groups on the SAM surface. As a result of these more favorable interactions, a relatively stable water layer was maintained over the SAM surface with minimal interactions between the surface and the peptide functional groups. A representative configuration from the simulation of G4-G-G4 over the OH SAM surface is shown in Figure 8. 4.2. Peptide Adsorption on the COOH SAM Surface. Figure 9 shows the SSD versus time, NPD vs SSD, and ∆Gi vs SSD plots for the G4-G-G4 and G4-K-G4 peptides over the COOH SAM surface. As clearly shown, and unlike the OH SAM surface, both peptides adsorbed to the COOH SAM surface, with the positively charged G4-K-G4 peptide adsorbing most strongly. As shown in the SSD vs time trajectory plot, the SSD value for each system was maintained at about 3.5 Å for most of the simulation, with further approach being prevented by van der Waals forces of repulsion. As a result of the strong interactions between the peptides and this surface, and as shown in the NPD vs SSD plots (Figure 9), the 10 ns simulations did not provide adequate time to enable the peptide to sample the full range of SSD values between 0 and 35 Å over the SAM surface. Consequently, the probability density could not be calculated accurately in the unsampled areas. Thus, the probability distribution does not properly account for the relative probability distribu-

Figure 8. A snapshot from the trajectory file of the G4-G-G4 peptide over the OH SAM surface. The sodium ions (dark spheres) and chloride ions (light spheres) are shown with a space-filling model, while the water molecules are not shown for the sake of clarity. In this system, the OH functional groups on the SAM surface interacted more strongly with the SPC water molecules than with the peptide, thus resulting in negligible peptide adsorption to the SAM surface.

tion of the peptide for the whole range of SSD, and the value of P0 that is needed for calculating ∆Gi could not be determined in the same manner as for the OH SAM surface. To account for this problem, P0 was defined as the NPD value of the SSDi location furthest from the SAM

1636


Raut et al.

Figure 9. SSD vs time, NPD (Pi) vs SSD, and ∆Gi vs SSD for both peptides on the COOH SAM: (A) G4-G-G4 peptide and (B) G4-K-G4 peptide.

surface, and this value was then used to calculate and plot the normalized free energy as a function of SSD, as shown in Figure 9, with the resulting values of the adsorption free energy for the G4-G-G4 and G4-K-G4 peptides on the COOH SAM surface calculated as -2.10 ( 1.61 and -3.75 ( 0.65 kcal/mol (mean ( std deviation), respectively. It must be realized, however, that the true NPD ratios that would have been obtained if the system had been adequately sampled (i.e., with P0 representing the actual average NPD value for SSD between 25 and 35 Å) would reflect larger values of Pi /P0 and subsequently more negative values of ∆Gi and ∆Gads per eqs 3 and 4. Thus, from these simulations, it is more appropriate to express the adsorption free energy values for G4-G-G4 and G4-K-G4 peptides on the COOH SAM surface as an upper bound, with ∆Gads < -2.10 ( 1.61 kcal/mol and ∆Gads < -3.75 ( 0.65 kcal/mol, respectively. From viewing the trajectory files for these peptideSAM systems, it was observed that this adsorption behavior was a result of several contributions from different types of nonbonded interactions between the functional groups of the peptides and the SAM surface. The dominant effect appeared to be the electrostatic attractions between the positively charged NH3+ functional groups of the peptide (lysine side-chain for the G4K-G4 peptides and the N-terminus of both peptides) and the COO- functional group of the surface. These interactions, however, always occurred with two or three intervening water molecules between the NH3+ and COOfunctional groups. The NH3+ functional groups of the peptide were therefore not tightly held to the surface but underwent a fair amount of motion over the surface while maintaining their position in general proximity to the charged COO- functional groups of the SAM surface. Electrostatic repulsion was also generally observed between the C-terminus of the peptides and the surface,

with the C-terminus generally being positioned furthest away from the surface while the peptide was adsorbed. In addition to electrostatic effects, relatively stable hydrogen bonding interactions were also observed to occur between both charged-polar and polar-polar functional group combinations between the peptide and the SAM surface. Unexpectedly, hydrophobic interactions were also frequently observed between the CH2 groups in the peptide (from both the glycine and lysine side-groups) and the CH2 groups immediately below the COO-/COOH groups on the SAM surface that were exposed when the COO-/ COOH surface groups occasionally separated during the simulation despite the fact that all other atoms of the SAM were fixed in position. Thus, the strong adsorption of these peptides to the COOH SAM surface was due to a combined result of electrostatic, hydrogen bonding, and hydrophobic interactions. The disagreement with experiment,28 along with the unexpected occurrence of hydrophobic contacts, indicates that perhaps the GROMACS force field overestimates the strength of hydrophobic interactions. Figure 10 presents a representative configuration from the simulation of the G4-K-G4 on the COOH SAM surface that illustrates each of these types of interactions. 4.3. Peptide Adsorption on the OEG SAM Surfaces. Simulations of each peptide were performed on two types of OEG SAM surfaces: trans and helical. It was readily apparent from the simulations that the trans conformation was much more stable than the helical one, with the trans OEG segments retaining their conformation throughout the 10 ns simulations, while the helical OEG segments shifted into a trans conformation within about the first 100 ps of the simulation and then maintained the trans conformation for the remaining 9.9 ns. Because of this, the trajectory data were essentially similar for both types of starting conditions, and only the trans OEG results


Langmuir, Vol. 21, No. 4, 2005 1637

Figure 10. A snapshot from the simulation of the G4-K-G4 peptide over the COOH SAM surface. The single atoms are sodium ions (dark spheres) and chloride ions (light spheres) contained in the solvent, with the water molecules not shown for clarity. In this system, the COO-/COOH functional groups on the SAM surface interacted with the peptide through a combination of electrostatic, hydrogen bond, and hydrophobic interactions resulting in strong adsorption of the peptide to the SAM surface.

Figure 11. SSD vs time, NPD (Pi) vs SSD, and ∆Gi vs SSD for both peptides on the OEG SAM: (A) G4-G-G4 peptide and (B) G4-K-G4 peptide.

will therefore be presented. Accordingly, Figure 11 shows the SSD vs time, NPD vs SSD, and ∆Gi vs SSD plots for G4-G-G4 and G4-K-G4 over the trans OEG SAM surface. While it was anticipated that this surface would be very nonadsorptive like the OH SAM, as shown in Figure 11, surprisingly both peptides interacted more strongly with this surface than either one did with the OH SAM surface, although less strongly than with the COOH SAM surface. While the G4-G-G4 peptide did sample the entire range of SSD during the simulation, the G4-K-G4 did not sample SSD values greater than 10 Å, remaining closely adsorbed to the surfaces during the whole 10 ns of the simulation. Thus similar methods had to be applied to calculate the

∆Gi and ∆Gads as were used for the COOH SAM surface, with these values again under-representing the absolute magnitude of the free energy values. Accordingly, from the data generated by the simulations, the adsorption free energy values (mean ( std deviation) for the G4-G-G4 and G4-K-G4 peptides on the OEG SAM surface were calculated to be -1.37 ( 2.03 and < -2.73 ( 1.54 kcal/ mol, respectively. An inspection of the trajectory files for the adsorption of these peptides on the OEG surface showed that while the interactions involved hydrogen bonding between the charged and polar functional groups of the peptides with the OH functional groups of the OEG SAM, by far the

1638


Raut et al.

Figure 12. A snapshot from the simulation of G4-K-G4 over the OEG SAM surface. Note that the OEG chains are generally structured in a trans conformation. The peptide is strongly adsorbed to the OEG surface via a combination of hydrogen bonds with the OH groups of the SAM surface and hydrophobic interactions between CH2 groups of the peptide and the CH2-CH2 segments of the OEG chain that became exposed by the separation of the surface OH groups during the simulation.

dominant cause of adsorption was due to hydrophobic interaction between CH2 groups of the peptides and the topmost CH2-CH2 segments of the OEG chain that were exposed when the top OH functional groups of the OEG SAM separated during the simulation. It was also observed that the water molecules did not penetrate the intermolecular spacing between OEG segments during the simulation to form polar bonds with the inner oxygens present in the OEG segments, which has been proposed as the mechanism that stabilizes the experimentally determined helical conformation of the OEG chains.29,30 Figure 12 shows a representative configuration from simulation of the G4-K-G4 peptide over the OEG surface to illustrate the typical types of bonding that occurred between the peptides and the OEG surface. 5. Discussion The results of these simulations show areas where the simulated behavior is both in agreement and in conflict with experimental results. As presented above, to obtain a satisfactory simulation of peptide behavior over the OH SAM surface, it was necessary to constrain all the atoms of the SAM except for the surface OH group. Once this was done, the simulated behavior was in excellent agreement with the experimental results of Vernekar and Latour,28 with negligible interactions between the peptide and the surface and a negligible adsorption free energy. However, unexpected hydrophobic interactions between the peptides and the subsurface CH2 groups of the SAM still occurred with the COOH SAM surface even when it was constrained. These interactions contributed to the relatively strong adsorption of the G4-G-G4 peptide to this surface, which was contrary to the prior experimental studies, that found no adsorption for this system. Peptide adsorption was also unexpectedly observed for the OEG surfaces, again largely due to hydrophobic interactions. These observations underscore the problem of force field transferability. The GROMACS force field was primarily parametrized for the simulation of peptides and proteins in aqueous solution and was not specifically designed to model the behavior of alkanes when arranged to form a SAM structure or OEG chain segments. This is also

reflected by the fact that new GROMACS residue-type definitions had to be defined in order to represent the molecules needed to construct the SAM surface. The parameters for these molecular segments were defined based on the parameter set for similar-type functional groups present in the GROMACS library for amino acids or, for the case of OEG, another force field (i.e., AMBER). It is apparent that the resulting parametrization satisfactorily represents neither the behavior of the SAM surface nor the adsorption of peptides to these SAM surfaces. Further modification and verification work is therefore necessary before this force field can be confidently used to simulate peptide and protein adsorption behavior on these types of surfaces. A second problem with the presented simulation methods was clearly demonstrated with the G4-K-G4 peptide over the COOH SAM, in which the system clearly was not sufficiently sampled in order to accurately calculate adsorption free energy. Referring to eqs 3 and 4, to represent the experimentally measured adsorption free energy of -7.0 kcal/mol at 300 K, the probability ratio of the likelihood of the peptide being adsorbed to the surface vs being far away from the surface in a nonadsorbed state (i.e., Pi/P0) would have to average over 100 000. However, a 10 ns simulation with data saved every 1 ps only provides 10 000 data points. Thus it is not surprising that the simulations did not sample SSD positions far from the surface in the strongly bound systems. To solve this problem, either much longer simulation times are needed (on the order of microseconds) or biased sampling methods, such as umbrella sampling,21,44-47 must be applied. We are currently developing the latter approach to address this issue for future simulations. (44) Ryckaert, J. P.; Bellemans, A. Mol. Dyn. Liq. Alkanes, Faraday Discuss. 1978, 20, 95-106. (45) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187199. (46) Jorgensen, W. L.; Gao, J.; Ravimohan, C. J. Phys. Chem. 1985, 89, 3470-3473 (47) Jorgensen, W. L.; Buckner, J. K. J. Phys. Chem. 1987, 91, 60836085.


6. Conclusions As stated, the objectives of this research were to develop and assess methods to simulate the adsorption of a hostguest peptide on functionalized alkanethiol SAM surfaces on gold using a selected molecular mechanics force field and to calculate the free energy of peptide-surface adsorption from the trajectory data for comparison with experimentally determined adsorption free energies. The GROMACS molecular simulation program and force field were selected for our initial studies, and our results have identified several areas in need of further development, involving both the simulation methods and the parametrization of GROMACS for this specific application. These simulation results do not allow us to make firm statements regarding whether the discrepancies between simulation and experiment are due to problems with the parametrization of the SAM surfaces, the peptide-SAM surface interactions, or both. However, we can conclude that the parameter set used in these simulations did not provide a satisfactory representation of the peptide-SAM molecular system. Additional parameter fitting is required to accurately represent SAM surfaces in GROMACS, rather than transferring parameters from other parts of GROMACS or other force fields. The alkane-alkane interactions will need to produce a SAM structure that can be simulated dynamically, without requiring unphysical constraints. It is also possible that new parameters, or parameter scaling factors, will also be needed to adjust interactions between the peptides and the SAM at the solid-solution interface in order to obtain the right balance of hydrophobic and hydrophilic interactions. Interactions at this interface are distinct from those within

Langmuir, Vol. 21, No. 4, 2005 1639

the SAM surface and the aqueous peptide solution themselves. Due to physical effects such as differences in effective polarizability between the bulk and interfacial regions, the interfacial interactions do not necessarily follow directly from using standard combining rules with classical, nonpolarizable potentials. Finally, we also conclude from our studies that biased sampling methods, such as umbrella sampling, are needed in order to efficiently calculate the free energy of peptide adsorption. This is due to the very low probability of sampling the nonadsorbed state for strongly adsorbing peptide-SAM systems. Once these issues are resolved, the developed simulation methods for peptide adsorption to SAM surfaces, coupled with the previously developed SPR experimental methods, should provide a very effective approach for the evaluation, modification, and validation of existing force fields for the simulation of peptide and protein adsorption behavior to functionalized surfaces. Once validated, a protein adsorption force field should have great potential for use as a tool for surface design to predict and control adsorbed protein orientation, conformation, and bioactivity. Acknowledgment. We thank the National Science Foundation (NSF) (Award Number EPS-0296165), the State of South Carolina, and Clemson University for providing the funding support for this project. We also thank the NSF Center for Advanced Engineering Fibers and Films (CAEFF), Ms. Corey Ferrier, and Mr. Tim Shelling at Clemson University for computational resources and computer system support. LA047807F

Molecular Dynamics Simulations of PeptideâSurface Interactions

Recommend Documents

Molecular Dynamics Simulations of PeptideâSurface Interactions