Simulation of Adhesion Forces and Energies of Peptides on

Alessandro Motta , Marie-Pierre Gaigeot , and Dominique Costa. The Journal of Physical Chemistry C 2012 116 (44), 23418-23427. Abstract | Full Text HT...
0 downloads 0 Views 3MB Size
Download PDF
pubs.acs.org/Langmuir © 2010 American Chemical Society

Simulation of Adhesion Forces and Energies of Peptides on Titanium Dioxide Surfaces Susan K€oppen†,‡ and Walter Langel*,† † ‡

Institut f€ ur Biochemie, Universit€ at Greifswald, Felix-Hausdorff-Strasse 4, 17487 Greifswald, Germany, and Hybrid Materials Interfaces Group, Faculty of Production Engineering and Bremen Center for Computational Materials Science, University of Bremen, TAB-Building, Am Fallturm 1, 28359 Bremen, Germany Received April 27, 2010. Revised Manuscript Received August 25, 2010

Force field molecular dynamics simulations on a decapeptide in contact with a rutile (100) surface in aqueous solution are reported. The peptide sequence is part of R1-collagen. Force-distance curves yield discrete peaks to the rupture of charged lysine and glutamate side chains from the surface according to the model of contact points. The rupture forces are 0.2-2.2 nN, and the values strongly depend on the charges of surface hydroxyl groups. Adhesion energies are evaluated from the areas of the rupture peaks. For proton charges of 0.4 and 1, adhesion energies between 40 and 190 kJ/mol were found being comparable to recent ab initio molecular dynamics results. Flips in the torsional angles of the peptide are observed during restrained desorption. The partial charges of hydroxyproline are revised, and the polarization of the Cβ-Cγ bond as well as the ring pucker conformations are taken into account. It is shown that transition from collagen to helix fold is more likely for hydroxyproline than for proline and might be relevant for differences in the properties of POG (proline, hydroxyproline, glycine) and PPG collagen sequences.

Introduction Recently, interfaces between inorganic surfaces and biomolecules attracted great interest,1 playing an important role for the biofunctionalization of implants, for example. Titanium based implants may be covered by specific proteins for improving the biocompatibility. Here, a detailed understanding of the interaction between the passivation layer and the coating is essential for its optimization. The surfaces of passivated Ti are commonly modeled both as rutile ((110) or (100)) and anatase ((101) or (001)).2 We focus on rutile (100) which tends to be the most reactive one and adopt a model for force field molecular dynamics simulations on the perfect partially hydroxylated surface from.3 Collagen coatings on titanium considerably enhance bone growth and thus favor osseointegration.4,5 The adhesion of the protein on the oxidized surface is driven by local interactions between side chains of those amino acids which are charged at physiological pH values and single surface charges such as Ti and O surface atoms or singly and doubly coordinated hydroxyl groups.6 The carboxyl groups of glutamic and aspartic acid have stable distances of 1.5 A˚ between the carbonyl oxygen and the surface proton. The side chain of lysine has an N-O distance of 2.8 A˚ with short living alternating hydrogen bonds between the amino protons and the surface hydroxyl oxygen atoms. Consequently, short peptide sequences with about 10-12 residues may already show specific binding properties to inorganic *To whom correspondence should be addressed. E-mail: [email protected]. (1) Morra, M. Eur. Cells Mater. 2006, 13, 1–16. (2) Diebold, U. Surf. Sci. Rep. 2003, 48, 53–239. (3) K€oppen, S.; Langel, W. Surf. Sci. 2006, 600, 2150–2160. (4) Teng, S.; Lee, E.; Park, C.; Choi, W.; Shin, D.; Kim, H. J. Mater. Sci.: Mater. Med. 2008, 20, 2553–2561. (5) Schliephake, H.; Aref, A.; Scharnweber, D.; Bierbaum, S.; Sewing, A. Clin. Oral Implant Res. 2009, 21, 31–37. (6) K€oppen, S.; Ohler, B.; Langel, W. Z. Phys. Chem. 2007, 231, 3–21. (7) Oren, E. E.; Tamerler, C.; Shahin, D.; Hnilova, M.; Seker, U. O. S.; Sarikaya, M.; Samudrala, R. Bioinformatics 2007, 24, 2817–2823.

15248 DOI: 10.1021/la101687m

surfaces,7 and adhesion properties of collagen may be studied by modeling typical short sequences in atomistic detail rather than by taking into account the extended strand itself. A typical feature of collagen sequences is the partial hydroxylation of proline (three letter code: PRO) to hydroxyproline (HYP) resulting in a significant increase of the melting temperature.8 As stated by Park et al.,9 proline and hydroxyproline have two different ring puckering conformations, which are distinguished by two torsion angles. The backbone angle Φ (CR-C-N-CR) is around -75° and -58° for the endo and exo forms, respectively, and the respective values of ξ in the side chain (N-Cδ-Cγ-Oδ) are 140° and 85°. Both in the monomer as well as in collagen-like, peptides the endo conformation is energetically favored for PRO whereas exo HYP has a lower energy than that of its endo form. The preference of a different pucker conformation in HYP might be one reason for the natural hydroxylation of proline in collagen.9 A straightforward experimental approach for quantifying the adhesion on surfaces is atomic force microscopy (AFM).10 Single biomolecules are grafted on the tip by silanization, and forcedistance curves are recorded.11 On nonpolar surfaces, no discrete bonds are pulled apart, but the hydrophobic effect results in plateaus in the force-distance curves being due to lateral friction.12,13 Bond breaking between molecules and polar surfaces results in peaks of sawtooth shape, whose maxima correspond to the rupture forces. Experimental values range from 1 nN for smaller systems with only few contacts such as bovine serum (8) Shoulders, M. D.; Raines, R. T. Annu. Rev. Biochem. 2009, 78, 929–958. (9) Park, S.; Radmer, R. J.; Klein, T. E.; Pande, V. S. J. Comput. Chem. 2005, 25, 1713–1717. (10) Butt, H. J.; Berger, R.; Bonaccurso, E.; Chen, Y.; Wang, J. Adv. Colloid Interface Sci. 2007, 143, 91–104. (11) Friedsam, C.; Gaub, H. E.; Netz, R. R. Biointerphases 2006, 1, 1–22. (12) Horinek, D.; Serr, A.; Geisler, M.; Pirzer, T.; Slotta, U.; Lud, S. Q.; Garrido, J. A.; Scheibel, T.; Hugel, T.; Netz, R. R. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 2842–2847. (13) Serr, A.; Horinek, D.; Netz, R. R. J. Am. Chem. Soc. 2008, 130, 12408– 12413.

Published on Web 09/13/2010

Langmuir 2010, 26(19), 15248–15256

K€ oppen and Langel

Article

Figure 1. (A) Structure and atom names of the (from left to right) protonated, stoichiometric, and hydroxylated titanium dioxide units. Units have the same state of hydroxylation on bottom and top.3 (B) Structure and atom names of the heavy atoms of the hydroxyproline residue (ball and stick).The two amino acids connected to the residue in the peptide are also indicated (lines).

albumin (BSA) on enamel14 to 1.2 μN for collagen on titanium dioxide.15 The data sensitively depend on the structure of the adsorbed molecule, and the addition of a single functional group can have a dramatic effect. Dihydroxyphenyl alanine (DOPA) on titanium dioxide has a rupture force of 0.8 nN, whereas removal of one OH group results in a reduction to 0.1 nN for tyrosine.16 Restrained molecular dynamics simulations as inspired by these single molecule AFM experiments yielded rupture forces of peptides in several systems.11-13,17-21 Peptides have been desorbed by fixing restraints at single CR atoms20,21 and pulling off the molecule from the surface by shortening the restraint with a constant velocity. The displacement rate of the set point is crucial for the effect on the protein. Typical experimental velocities are as low as 0.4-40 μm/s.10 This is too costly for molecular dynamics simulations where rates as high as 0.1-50 m/s have to be employed.17,18,21 The dependence of the forces on pulling rates between 0.1 and 10 m/s has been investigated for spider silk on hydrophobic surfaces, and stable values of F = 54 ( 15 pN could only be found below 0.1 m/s.12 In ref 20, fibronectin has been stretched with 50 m/s inducing denaturation, which was accompanied by flips of some of the Φ and Ψ angles from (-16°, -67°) to (-113°, -5°). The adhesion energy is an essential parameter for adsorption configurations of proteins on solid surfaces. Energies of localized contact points are calculated for single amino acids and small hydrocarbons on rutile and anatase by constrained desorption using first principles molecular dynamics.22-24 Bonding was stronger on hydroxylated TiO2 than on stoichiometric TiO2.23 As long as the interaction is mainly of electrostatic origin, with some enhancement by hydrogen bonding, force field simulations should be adequate for yielding microscopic structural information and reproducing adhesion energies in larger systems. To the best of our knowledge, the energetics of the force distance curves from force field simulations have not been discussed so far, and we present a method for extracting energies for single surface (14) Schwender, N.; Huber, K.; Al Marrawi, F.; Hannig, M.; Ziegler, Ch. Appl. Surf. Sci. 2002, 252, 12. (15) J€ahnchen, G.; Mertig, M.; Pompe, W. Surf. Interface Anal. 1999, 27, 444– 449. (16) Lee, H.; Scherer, N. F.; Messersmith, P. B. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 13999–14003. (17) Heymann, B.; Grubm€uller, H. Chem. Phys. Lett. 1999, 305, 213–219. (18) Heymann, B.; Grubm€uller, H. Chem. Phys. Lett. 1999, 303, 1–9. (19) Heymann, B.; Grubm€uller, H. Phys. Rev. Lett. 2000, 84, 6135–6139. (20) Krammer, A.; Lu, H.; Schulten, K.; Vogel, V. Proc. Natl. Acad. Sci. U.S.A. 1999, 96, 1451–1456. (21) Lu, H.; Isralewitz, B.; Krammer, A.; Vogel, V.; Schulten, K. Biophys. J. 1998, 75, 662–671. (22) K€oppen, S.; Langel, W. Phys. Chem. Chem. Phys. 2008, 10, 1907–1915. (23) K€oppen, S.; Bronkalla, O.; Langel, W. J. Phys. Chem. C 2008, 112, 13600– 13606. (24) Langel, W.; Menken, L. Surf. Sci. 2003, 538, 1–9.

Langmuir 2010, 26(19), 15248–15256

contacts from force field simulations in aqueous solution. Essentially, sawtooth peaks for bond rupture in force distance curves are integrated numerically. The aim is to calibrate molecular dynamics (MD) simulations in comparison to density functional theory (DFT) work and to establish models for the contact of larger proteins to inorganic surfaces, which can be cross checked with experiments, such as single molecule atomic force microscopy.

Methods The generation of the perfect rutile (100) surfaces is described in detail in ref 3. A TiO2 slab is built up of 10 10 Ti5O10 units (Figure 1A), which consist of 5 titanium and 10 oxygen atoms and have an area of 4.595  2.95 A˚2 each and a height of 11 A˚. The atomic positions are fixed and the slab is positioned parallel to the xy-plane covering the full cross section (45.95  29.5 A˚2) of the simulation cell. Under periodic boundary conditions, a contiguous infinitely extended surface is mimicked. In the bulk, each titanium atom is octahedral coordinated by six oxygen atoms, which are in turn threefold coordinated to titanium. On the surface, in plane titanium and bridging oxygen atoms are fivefold and twofold coordinated, respectively. According to two pH dependent equilibriums,2 the open valences are in part saturated by protonation of the O and hydroxylation of Ti atoms. The position of the hydroxyl O corresponded to that of the missing sixth bulk oxygen.3 Proton positions are taken from Car-Parrinello molecular dynamics (CPMD) calculations24 and are frozen during the force field simulation (Figure 1A). All surfaces are adjusted to the physiological pH value of 7.4, where the ratio of protonated to neutral to hydroxylated units is 17:51:32.3 The electrostatic interaction of the peptides with the surface, which is prevailing, fully depends on the parametrization of the partial charges of the surface atoms. Charges of the bulk atoms were fixed to 1.152 (Ti) and -0.576 e (O),25 but three different models have been applied for the hydroxyl groups: (1) Formal charges of -2 for oxygen and þ1 for hydrogen, as in earlier work.3 They clearly describe an upper limit of the surface polarity generating a highly hydrophilic surface.26 (2) Partial charges from the fit of the electrostatic potential (ESP)27-29 to DFT calculation (Table 1) using the Lautrec code30 with the Car-Parrinello (CP) method.31 The PW91 exchange correlation generalized gradient (25) Swamy, V.; Gale, J. D. Phys. Rev. B 2000, 62, 5406–5413. (26) Ohler, B.; Langel, W. J. Phys. Chem. C 2009, 113, 10189–10197. (27) Cox, S. R.; Williams, D. E. J. Comput. Chem. 1981, 2, 304–323. (28) Chandra Singh, U.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129–145. (29) Woods, R. J.; Khalia, M.; Pell, W.; Moffat, S. H.; Smith, V. H. J. Comput. Chem. 1989, 3, 297–310. (30) De Vita, A. D.; Canning, A.; Galli, G.; Gygi, F.; Mauri, F.; Car, R. EPFL Supercomput. Rev. 1994, 6, 22. (31) Car, R.; Parrinello, M. Phys. Rev. Lett. 1985, 55, 2471–2474.

DOI: 10.1021/la101687m

15249

Article

K€ oppen and Langel

Table 1. Partial Charges in Units of the Elementary Charge e for Hydroxylated Rutile (Atom Names as in Figure 1A) as from (1) Formal Charges, from (2) ESP Calculations, and from the (3) TIP3P Water Model As Described in the Texta protonated

stoichiometric

hydroxylated

TiO2 All models -0.576 1.152

bulk O1,2,3,4,5,7,8,10 bulk Ti1,3,4 Formal, TIP3P

-0.576 1.152

surface O6,9 surface Ti2,5 ESP surface O6,9 surface Ti2,5

-0.576 1.252

-0.576 1.152

-0.676 1.152

Hydroxylation Formal surface OH surface H

þ1

-

-2/þ1 -

-

-0.87/0.47 -

ESP surface OH surface H

þ0.4 TIP3P

surface OH -0.834/0.417 surface H þ0.417 a For the ESP model the values for the surface Ti and O atoms were slightly adjusted providing integer total charges for the protonated and hydroxylated slabs.

(3)

approximation (GGA) functional32 and the projected augmented wave (PAW) method33 for charge evaluation are used as described in ref 34. The systems with a 2  2 surface consists of four units as described above. Two of them are protonated or hydroylated, respectively, and the others are stoichiometric. The minimizations of the electronic states and the geometry optimizations have been performed under periodic boundary conditions. This approach yielded charges not only for the surface hydroxyl groups but also for the bulk TiO2. The Ti and O charges in one unit were strongly scattering (0.9-1.7 for Ti), and for consistency with our previous work we used the charges from ref 25 instead. Partial charges of -0.834 and þ0.417 for O and H, respectively, as in the transferable intermolecular water potential with three points (TIP3P). The third approach assumes that no charge transfer occurs during the water dissociation on the surface. Its charge values are very similar to the DFT results (Table 1).

A fragment of a collagen type I triple helix6 with the sequence N-GLY101-PRO102-HYP103-GLY104-GLU105-LYS106GLY107-PRO108-HYP109-GLY110-C was positioned on the TiO2 slab. The 10 residues of the decapeptide have been numbered 101 to 110 and contain two -Xaa-Yaa-GLY- triplets. The Xaa and the Yaa positions are occupied by proline (single letter code: P) and hydroxyproline (O), respectively, resulting in POG triplets with glycine (G). The contacts to the surface are found at a (32) Perdew, J. P.; Wang, Y. Phys. Rev. B 1992, 45, 13244–13249. (33) Bl€ochl, P. E. Phys. Rev. B 1994, 50, 17953–17979. (34) Schneider, J.; Ciacchi, L. C. Surf. Sci. 2010, 604, 1105–1115.

15250 DOI: 10.1021/la101687m

Table 2. Data from AMBER9435 and from RESP Calculation in This Work for Proline and Hydroxyproline, Respectivelya atom

proline (PRO)

hydroxyproline (HYP)

N -0.255 -0.221 C 0.590 0.656 O -0.575 -0.478 CR -0.027 -0.073 HR 0.064 0.059 Cβ -0.007 -0.384 Hβ1 0.025 0.126 Hβ2 0.025 0.126 Cγ 0.019 0.335 Oγ1 -0.617 Hγ1 0.023 0.383 Hγ2 0.023 0.065 Cδ 0.019 -0.101 Hδ1 0.039 0.062 Hδ2 0.039 0.062 a For atom names, see Figure 1B; hydrogen atoms are indexed according to the heavy atoms to which they are connected.

charged EKG triplet sequence in the middle of the decapeptide with the residues 105 and 106 consisting of glutamic acid (E) and lysine (K). The HYP residue is derived from the predefined PRO in the AMBER force field35 by adding a hydroxyl group to the Cγ atom (Figure 1B). The charges are recalculated according to the standard restrained electrostatic potential (RESP) procedure using the Gaussian98 program36 (Table 2). The Cβ-Cγ bond is strongly polarized by the hydroxyl group, and the Cβ atom has a significant negative charge, which is in accord with chemical intuition. The partial charges for hydroxyproline from refs 9 and 37 were not adopted here, since they do not reflect this polarization. The cell with an initial height of 150 A˚ has been filled up with about 5000 TIP3P water molecules. The lengths of bonds involving hydrogen are fixed by a SHAKE algorithm.35 Filling of the simulation cell by standard procedures resulted in incomplete packing. Standard density is attained by reducing the cell height in steps of 0.1 A˚ and redistributing the water molecules by short energy minimization and runs at constant particle number, volume, and temperature of 300 K (NVT), keeping the protein and the titanium dioxide fixed.3 NpT runs at constant pressure (p = pressure) would in general result in an inclination of the z-axis of the cell with respect to the orientation of the TiO2 slab destroying the contiguous oxide surfaces. The water orientation, the ion concentration, and the diffusion coefficient of water and ions is affected up to 15 A˚ away from the surface. Subtracting the thicknesses of the TiO2 slab and of the two water layers on both sides of it from the final total cell height of 125.2 A˚, the depth of bulk water comes out to at least 80 A˚. Simulations have been performed using the program sander in the AMBER package. After cell initialization as described above, NVE production runs at constant energy (E = energy) with 5  105 steps corresponding to 0.5 ns were started for relaxation and formation of stable contact points. Three-dimensional plots and (35) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollmann, P. A. J. Am. Chem. Soc. 1995, 128, 5189–5207. (36) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; € Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Daniels, A. D.; Farkas, O.; Gaussian09, revision G98; Gaussian, Inc., Wallingford, CT, 2009. (37) Klein, T. E.; Huang, C. C. Biopolymers 1999, 49, 177–193.

Langmuir 2010, 26(19), 15248–15256

K€ oppen and Langel

Article

Figure 2. (A) Scheme of restrained desorption: A CR atom within the adsorbed peptide is connected to a fixed atom by restraint with a force constant of krest = 14 N/m. The distance s0 is shortened with a pulling rate of 2 m/s, and the restraint applies a force to the CR atom of Frest = krestΔs according to Hookes0 s law. (B) (top) The restraint stores energy itself, Erest, due to difference between set point and actual value and induces elastic deformation of the molecule with the energy, Eelast,i, both indicated as springs. At the leading edge of a rupture peak, Erest increases due to the increasing difference between s0 and actual distance s. Elastic energy may be stored in the adsorbed molecule when pulling one CR atom and fixing another atom to the surface. (bottom) After desorption, both restraint and elastic energies are released. For desorption, the restraint energy Erest must at least be equal to the adhesion energy Eadh. distance over time plots were generated with the VMD package,38 and quoted distances are averages over NVE runs with error bars due to thermal fluctuations. Adhesion forces and energies have been measured by induced desorption in NVT runs at 300 K. A harmonic restraint with a force constant krest has been defined between one CR of the peptide and a fixed water molecule in the top part of the simulation cell (Figure 2A). By reducing the set point s0 for this restraint with a typical velocity of s0 = 2 m/s = 1 A˚/50 ps, the molecule is lifted from the surface into the water on top of it. In general, the actual value s of the restrained distance is back lagging s0, and the force on the molecule due to the restraint is given by Hooke’s law: F ¼ krest ðs - s0 Þ ¼ krest Δs

ð1Þ

From such trajectories, force-distance curves of F over the actual height s of the respective CR atom are obtained. Forces were read out as maxima in the force distance curves after smoothing by convolution with a Gaussian of 0.11 A˚ full width at halfmaximum (fwhm). R Pulling on the restraint results in energy storage of Erest = Frest ds, which is the integral of the applied force over the distance, by which the CR is lifted during desorption. In addition to that, the peptide may be deformed under stress and take up the elastic energy Eelast,i. Obviously, at least the adhesion energy Eadh must be transferred to the adsorbate by the restraint for desorption, but the elastic energy of the adsorbate molecule after desorption, Eelast,f, may be smaller than the initial one due to relaxation after stress release (Figure 2B). Consequently, the energy balance reads Eadh þ Eelast, i ¼ Erest þ Eelast, f w Eadh ¼ Erest þ ðEelast, f - Eelast, i Þ e Erest That is, the adhesion energy is a lower limit to the observed restraint energies. The force distance curves are recorded with a resolution of 0.02 A˚, and the respective sawtoothed peaks in the raw data were integrated numerically. The area of well separated peaks does not depend significantly on the choice of the integration limits, with the average force outside the peaks being very close to zero. (38) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33–38.

Langmuir 2010, 26(19), 15248–15256

Reasonable variations of the limits resulted in deviations of less than 1 kJ/mol. Overlapping peaks with two separate maxima are split into two integrals, taking the minimum between the peaks as upper and lower integration limit for the first and second one, respectively. If there is no valley between two peak components, only the sum of both contacts was evaluated. For checking the reproducibility of the quoted energies, several parameters have been scanned: The choice of the force constant of the restraint determines the accessible range of forces. Thermal fluctuation around the restrained distance is induced by collisions of the peptide with water molecules. This random motion superposes to the directed displacement of the molecule resulting in noise with an amplitude Δsn on the function of Δs over the simulation time t. Consequently, only forces can be measured, which are larger than Fn = krestΔsn resulting in values of Δs at least comparable to Δsn. The theoretical average noise amplitude at temperature T is given by   pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi krest krest Fn 2 kB T w Fn ¼ krest kB T ð2Þ ðΔsn Þ2 ¼ ¼ 2 2 2 krest where kB is the Boltzmann constant. Noise amplitudes observed in the trajectories are well consistent with the estimate given in eq 2. Obviously, recording small forces affords a low krest. On the other hand, Δs should be small compared to the heights of the simulation cell of about 90 A˚, that is, not larger than 5 A˚, for the maximum forces occurring, which gives an upper limit for krest. We used krest = 14 N/m, which allows forces to be evaluated between 0.25 and 5 nN. A further crucial parameter is the pulling velocity. Computational costs impose a lower limit, but high velocities may lead to erroneous results and even denaturation. We have traced the secondary structure during desorption by Ramachandran plots of single amino acids along the trajectory. Even without denaturation, rupture forces in discrete water might be overestimated at high velocities by an additional Stokes force. For a single spherical amino acid with a molar volume of 100 cm3 and a radius of 0.34 nm in water (10-3 Pa 3 s), the force is 14 pN at 2 m/s, which is lower by about 1 order of magnitude than the total forces calculated here. Finally, force-distance curves were recorded for different pulling velocities to check the convergence. DOI: 10.1021/la101687m

15251

Article

K€ oppen and Langel

Dependence of the results on the starting geometry was checked by using different end points of the restraint. During the simulations, runs pulling at each of the CR atoms were performed. The other end point was a water molecule with a height of about 40 A˚ above the surface. In general, the molecule was centered in the x-y-direction in the cell. As a check of reproducibility, the water end point was shifted by about (20 A˚ in the x-direction in some runs, modifying the direction in which the molecule is pulled away from the surface.

Results and Discussion Contact Points of the Decapeptide. Simulations on the decapeptide resulted in four stable contact points for all charge models (Figure 3). The negatively and positively charged GLU 105 and LYS 106 side chains attach to protonated bridging oxygen atoms and singly coordinated hydroxyl groups on the surface, respectively, as observed before in refs 6 and 23. Similar stable configurations have been found in 3 ns MD simulations of three different collagen fragments on the (100) rutile surface, and we consider the observed surface contacts as typical. Hydrogen bonds of glutamic acid to surface hydroxyls with formal charges were very strong having lengths around 1.4 A˚. This value increased to 1.8-2.0 A˚ for the more realistic charge models. The influence of charge reduction was much smaller on the lysine contact, with the N-O distance slightly increasing from 2.8 A˚ to 2.9-3.0 A˚. In addition to that, contacts between the backbone oxygen of GLY 104 and between the charged amino group of the N-terminal GLY 101 and the surface were stable here. These two types of contact points will be less important in larger folded peptides. Above a certain number of residues, the interaction of the terminal groups with the surface does not play much role as compared to the adhesion of the side chains, and folding hinders contacts between the backbone and surface for steric reasons. Force-Distance Curves. All force-distance curves for the decapeptide had well-defined peaks with sawtooth shape. The rupture forces on highly polar surface (1) and the respective signal/ noise ratios are high (Figure 4A). On the two surfaces with reduced hydroxyl charges, some signals were close to the noise limit. If only one single contact is broken and the final state of the peptide is sufficiently relaxed, an isolated peak is obtained (see GLY 101 in Figure 4A), and its area reproduces the respective adhesion energy of this contact point. In several cases, peaks from neighboring contacts overlap. When pulling at GLY 107, for example, only the GLY 101 (N) signal is clearly separated, whereas rupture peaks for LYS 106, GLY 104 (O), and GLU 105 overlap, and only the sum of the energies for the latter three contacts is evaluated. (Figure 4A). We first used a force constant of 7 N/m for the restraint which was close to the value of 4.2 N/m in ref 20. In the resulting forcedistance curves, the maxima did not clearly correlate with the rupture of the surface contacts (Figure 4B). Instead, some of the peaks split into a small one at the point of rupture and a second one shifted by typically 100 ps. This effect may be rationalized by assuming a system consisting of two force constants for a restraint spring between the pulling point and the moving mass and for a bond spring between moving mass and surface (Figure 2). Before rupture both springs are under stress. When the bond to the surface breaks, this spring suddenly relaxes, and the restraint is first relaxed and then stressed again like a yo-yo. We speculate that such an effect is at the origin of the observed peak splitting. Obviously, it would depend on the sizes of the moving mass and the respective strings in each individual system. The problem was indeed cured here by doubling the force constant. 15252 DOI: 10.1021/la101687m

Figure 3. Decapeptide on the rutile (100) surface. The four residues forming stable contact points are indicated by a ball and stick drawing: GLU105 with its CO2- group and GLY 104 with the backbone O pointing toward protonated bridging atoms, and LYS 106 and N-terminal GLY 101 with positive NH3þ groups in contact with singly coordinated surface hydroxyl groups. Water molecules are omitted for clarity.

In ref 12, very small pulling rates had to be employed for measuring forces on a hydrophobic surface. We used a significantly higher rate of 2 m/s and had to check how far this might affect the simulated energies and forces. For two systems (restraint at GLU105 and PRO 108), the peptide was pulled from the surface with four different velocities of 2, 1, 0.5, and 0.1 m/s, and fully convergent force distance-curves were obtained (Figure 4B) for the surface model (1). The lower forces for surface models (2) are consistent but affected by noise. Possibly the low forces of about 50 pN as found in ref 12 for hydrophobic surfaces were much more affected by a high pulling rate than the data above 0.2 nN as recorded here for electrostatic interactions with the surface. The dependence of the force distance curves and adsorption energies on the pulling geometry was checked by starting with identical configurations but fixing the restraint at different CR atoms in the decapeptide chain. Figure 4C shows the two forcedistance curves for restraints on the CR atoms of the HYP 103 and HYP 109 close to the C- and N-terminals, respectively, using model (2). The order of the rupture of LYS 106, GLU 105, and GLY 101 is reversed by pulling at the opposite end of the peptide. In most cases, contact points which are closer to the respective CR atom break first. A long strand between pulling point and contact has a small effective force constant and is significantly stretched before rupture is attained. Varying the H2O molecule, at which the restraint is fixed, had a small impact on the results, modifying the forces and energies by only (20%. Forces. Forces for the formal charge model (1) are higher by about a factor of 2-7 (Table 3) than forces for reduced charges, (2) (Table 4) and (3) (Table 5). The forces for the system with the TIP3P derived charges are slightly smaller than those for the DFT model, in consistence with the charge values (Table 1). For the contact of glutamic acid, the forces of 2.2 nN (surface (1)) and 0.3 nN (surfaces (2) and (3)) only weakly depend on the restrained CR. In some cases, the carboxyl group rotated by 180° exchanging the O atom in contact to the surface while pulling delaying desorption by about 50 ps. This process was not observed in stable adsorption geometries, and we attribute it to a reaction of the glutamic acid side chain to mechanical stress. Lysine data show significant dependence on the pulling direction with respect to the surface in all cases. When restraining a CR near the N-terminal, forces are significantly lower than those for pulling close to the C-terminal (Tables 3-5). As stated in ref 24, Langmuir 2010, 26(19), 15248–15256

K€ oppen and Langel

Article

Figure 4. Force-distance curves showing rupture peaks for the contact points described in the text. GLU and LYS indicate desorption of the GLU 105 and LYS 106 side chains. GLY(N) and GLY(O) correspond to N-terminal GLY 101 and backbone O in GLY 104. Curves are rescaled as indicated and shifted in the y-direction for clarity. Relative shifts of the contact points in different traces are indicated by lines. (A) Influence of the hydroxyl charges. Restraints with 14 N/m were applied to GLY107. Curves refer to (1) formal charges, (2) DFT calculation, and (3) TIP3P charges (from top to bottom). For (1), the raw data before (black) and after smoothing by convolution with a Gaussian of fwhm = 0.11 A˚ (red) are plotted; for all other data in this figure, only the smoothed curves are shown. (B) Traces for different pulling rates of 2 m/s (black), 1 m/s (red), 0.5 m/s (green), and 0.1 m/s (blue) clearly converge, indicating that the pulling rate of 2 m/s as applied in most calculations in this work is appropriate for the systems studied here. Restraints with 14 N/m were applied to GLU105. From top to bottom: Charge model (1), force constant 14 N/m; model (2), force constant 14 N/m; model (2), but 7 N/m. The spurious peak as seen for 7 N/m is indicated by an arrow. (C) Variation of the restraint CR modifies the force distance curves and even reverses the order of the peaks. Model (2), 14 N/m, restraints at HYP 103 near N-terminal (top) and at HYP 109 close to C-terminal (bottom). Table 3. Rupture Forces and Adhesion Energies of Four Contact Points (GLY 101, GLY 104, GLU 105, LYS 106) of the Decapeptide on the Partially Hydroxylated Surface during Restrained Desorption Pulling at the Cr Atoms of Different Residues on Surface Model (1) with Formal Charges for Surface Hydroxyl Groupsa GLY 101 (N) restraint

GLY 104 (O)

GLU 105

force/nN energy/kJ/mol force/nN energy/kJ/mol force/nN energy/kJ/mol

LYS 106 force/nN

energy/kJ/mol

PRO108 0.7 120 f 2.25 220 2.1 315 GLY107 1.3 225 2.25 92 2 142 2.6 370 LYS106 1.25 240 f 2.25 277 2.5 301 GLU105 0.9 75 1.6 108 2.3 213 1.5 170 GLY104 1.5 140 1 85 2.4 365 1.5 175 HYP103 2.9 329 1.9 160 2 190 1.25 186 PRO102 1.4 170 1.1 86 2.35 208 1.25 162 GLY101 2.1 147 1.1 106 2.1 297 1.3 172 lower limits 1.1 ( 0.2 120 ( 30 1.1 ( 0.1 90 ( 10 2.2 ( 0.2 190 ( 50 2.4 ( 0.2 (C-term); 1.4 ( 0.2 (N-term) 170 ( 20 a Entries in bold contain presumed significant contributions of elastic energy in the peptide chain. Lower limits of forces and adhesion energies are assumed to describe desorption only (see text). Arrows indicate overlap of the respective second rupture peak on this entry.

Table 4. Rupture Forces and Adhesion Energies of Four Contact Points (GLY 101, GLY 104, GLU 105, LYS 106) of the Decapeptide on the Partially Hydroxylated Surface Model (2) during Restrained Desorption Pulling at the Cr Atomsa restraint

GLY 101 (N)

GLY 104 (O)

GLU 105

LYS 106

GLY110 0.65 175 0.3 30 f 0.6 HYP109 0.55 97 f f 0.4 PRO108 0.5 113 f 0.4 70 0.5 GLY107 0.5 109 f f 0.7 LYS106 0.6 150 f f 0.7 GLU105 0.5 95 f 0.6 63 0.4 GLY104 0.7 137 f 0.5 63 0.4 HYP103 0.8 97 f 0.2 30 0.35 PRO102 0.4 85 r 0.3 56 0.3 GLY101 0.5 80 f f 0.3 lower limits 0.5 ( 0.1 90 ( 10 0.3 30 0.3 ( 0.1 60 ( 10 0.6 þ 0.2 (C-term); 0.4 ( 0.1 (N-term) a Charges on surface model (2) are derived from DFT calculations (entries in bold; see Table 3).

the threefold symmetry of the amino group and the two- or fourfold symmetries of the (100) rutile sites do not match. The amino group is attached to a singly coordinated hydroxyl group embedded in two rows of bridging oxygen atoms whose Ti-O bonds are not directly perpendicular to the surface but slightly inclined (Figure 3). During the simulation time, the NH3þ group is reorienting along the C-N-axis, tending to form short-lived hydrogen bonds to surface oxygen atoms.6 When restraining CR atoms close to the C-terminal end, lysine detaches in one step. While pulling at the other end of the peptide, the amino group of Langmuir 2010, 26(19), 15248–15256

90 105 55 140 137 52 56 50 53 115 50 ( 10

lysine is elongated from the single coordinated hydroxyl group but finds an intermediate contact on the adjacent row of bridging atoms. The final desorption takes place from there, which results in smaller maximum forces and wider peaks with two close maxima. The rupture forces of 1.1, 0.5, and 0.4 nN for the N-terminal amino group of the peptide when pulling at CR atoms close to it are similar to those for the lysine side chain contact. A major exception is the result from restraint at HYP 103 on the surface (1) with F = 2.9 nN. The trajectory indicates that in this DOI: 10.1021/la101687m

15253

Article

K€ oppen and Langel

Table 5. Rupture Forces and Adhesion Energies of Four Contact Points (GLY 101, GLY 104, GLU 105, LYS 106) of the Decapeptide on the Partially Hydroxylated Surface Model (3) during Restrained Desorption Pulling at the Cr Atomsa restraint

GLY 101 (N)

GLY 104 (O)

GLU 105

GLY 110 0.4 90 0.2 20 0.35 HYP 109 0.25 50 r 0.35 PRO 108 0.53 103 0.25 22 0.37 GLY 107 0.4 70 f f LYS 106 0.58 82 f 0.37 GLU 105 0.42 76 0.2 10 f GLY 104 0.63 75 f 0.35 HYP 103 0.25 26 0.28 27 0.23 PRO 102 0.48 53 r 0.35 GLY 101 0.4 30 r 0.42 lower limits 0.4 ( 0.1 65 ( 20 0.2 ( 0.1 20 ( 10 0.3 ( 0.1 Ab initio MD23 a Charges on surface model (3) are derived from TIP3P water model.35

configuration considerable elastic energy is stored in the backbone deformation of the peptide and suddenly released during rupture of the N-terminal. Energies. Energies from numerical integration of rupture peaks for nonoverlapping peaks are compiled in Tables 3-5. The values cluster at a lower limit, which is quoted as adhesion energy. The energies for the formal charge model are three to four times higher than those for the two other surface models, being of the order of magnitude of covalent bonding. For DFT or TIP3P charges, the adhesion energies of the N-terminal are significantly higher than those of the lysine side chain (90 and 65 kJ/mol vs 50 and 55 kJ/mol). This does not correspond to the forces, which are equal for both N-containing groups, but is consistent with the predefined partial charges. They are -0.60 and 0.45e for N and H, respectively, for the N-terminal in glycine and -0.25 and 0.29e in the lysine side chain.31 Only few of the exceptionally high energy values for model (1) with high charges (Table 3, bold entries) are due to overlap of rupture of different contact points, since in many cases sharp welldefined peaks were obtained. Additional elastic deformation of the desorbing molecule is likely to take place, and the energy introduced by pulling at the CR possibly is also stored in other degrees of freedom than the bond to the surface. The situation is different for the weaker peaks on the surfaces with reduced charges, and most energy values, which are significantly higher than the respective lower limit, include contributions of a second contact. This is not the case for the two energies for GLY 101 with restraints at GLY 110 and LYS 106, and they are attributed to significant deformation of the peptide during desorption. It is clearly seen that the force field reproduces the preferential adsorption of charged side chains. The energies for GLU and LYS are 60 and 50 kJ/mol, respectively, for surface model (2) and 40 and 55 kJ/mol for surface model (3). The small adhesion energies of 30 and 20 kJ/mol for the GLY 104 on surface (2) and (3) are subject to large error bars, since they are close to stronger peaks of glutamic acid or lysine contacts in most cases. Flip of the HYP Backbone. Denaturation of proteins by rapid stretching may start with flips of torsional angles in the peptide backbone. Several flips of the N-CR-C-N backbone torsion angle ψ between R- and β- areas are found in the Ramachandran plots (Figure 5A). ψ angles of the HYP residues 103 and 109 in Yaa position reversibly switch between values around 170° and -35° (Figure 5B), frequently even in both residues during the same simulation. This corresponds to a local transition between initial collagen and R-helical structure (Figure 5C). As a result, the HYP backbone oxygen and the nitrogen of the next GLY reorient, forming a new hydrogen bond between the N-H of glycine and the ring nitrogen atom of the hydroxyproline. The 15254 DOI: 10.1021/la101687m

46 25 40 68 42 23 54 68 40 ( 10 110

LYS 106 0.4 0.52 0.3 0.45 0.6 0.52 0.45 0.3 0.3 0.26 0.4 ( 0.1 (C-term);0.3 ( 0.1 (N-term)

74 60 40 80 66 60 50 54 44 30 55 ( 20 160

process as induced by pulling at the peptide is not a simple stretching of the backbone, since the N-N distance across the HYP and GLY is even reduced from 3.7 to 2.8 A˚ , but may rather be considered as a transition from one local energy minimum to another one triggered by the external perturbation. These flips occur at all pulling rates between 0.1 and 2 m/s having no significant correlation to the velocity. The intervals of about 5-20 A˚ between the flips are related to the shortening of the restraint which probably induces static structural deformations (Figure 5B). If a flip takes place simultaneously to detaching of a side chain from the surface, the energy for overcoming the barriers between these two states may contribute to the observed peak area (Tables 3 and 4, bold entries). When detaching, for example, LYS 106 by pulling at its CR, the flip occurs after 140 ps, whereas the desorption peak extends from start to 210 ps. In the relaxation of the molecule under stress, several degrees of freedom may be involved, and no discrete energy values can be assigned to single flips. In our simulations, only the HYP at the Yaa position flipped, whereas the torsional angles ψ for the neighboring Xaa PRO were stable (Figure 5A). This seems to be rather due to the structure of the HYP and PRO residues than to differences in the environment of the Xaa and the Yaa positions. When replacing HYP by PRO (PPG), the flip occurred only in one PPG triplet rather than in both POG and was often delayed with respect to the equivalent trajectory with HYP. Φ and ξ values in our simulations were similar to those in ref 9, making clear assignments of the HYP and PRO residues to the exo/endo configurations possible (Figure 5B). The hydroxyprolines switched from the exo to the endo form responding to external stress during shortening the restraint. In contrast to that, we did not find stable exo intervals for proline, for which the endo/ exo energy difference is 3 times as high as that for the hydroxyproline molecule. Monomer values are 6 versus -2 kJ/mol, and energy differences might be larger in a more rigid longer peptide. In parallel to Φ and ξ, we monitored here the flips of the Ψ torsion angles which were fixed in the simulations of Park et al.9 There is an obvious correlation between endo/exo transition and Ψ flip (Figure 5B), and we speculate that in HYP the latter induces a transition to the less stable endo form. In proline, the reverse transition from endo to exo is unlikely due to the high energy difference, and this might be related to the observation that Ψ flips are less frequent here. In one case, we found a reversible decrease of the number of hydrogen bonds of the peptide to the solvent during desorption at 2 and 1 m/s which was not seen at the lower pulling rates (Figure 6). This effect depends on the pulling velocity and thus is not correlated to the Ψ flips of hydroxyproline. In spite of the different behavior of the hydrogen bonds at different velocities, Langmuir 2010, 26(19), 15248–15256

K€ oppen and Langel

Article

Figure 6. Integrated RDF of water hydrogen atoms within 2.25 A˚ around the side chain hydroxyl O of HYP as a function of time when restraining GLU 105 during desorption. Consecutive traces are shifted by 2 units on the y-axis. At 2 and 1 m/s, the hydrogen bonds are opened reversibly during the desorption process (arrows). This is not visible for the lower pulling rates (coloring as in Figure 4).

Conclusion

Figure 5. Flip of the torsion angle ψ in hydroxyproline. (A) Ramachandran plots for HYP and PRO residues (as indicated in the bottom right corners) during desorption by restraining GLU 105. For the HYP 109 (POG sequence) data, points both in the collagen and close to the R-helix ranges are seen indicating flips between both folds. PRO 108 shows stable torsion angles of ψ = 140-190°. (B) ψ of HYP 109 as a function of time when restraining GLU 105 during desorption. The traces correspond to different pulling rates and time intervals for each frame, respectively: (black) 2 m/s and 1 ps, (red) 1 m/s and 2 ps, (green) 0.5 m/s and 4 ps, as well as (blue) 0.1 m/s and 20 ps. The traces are shifted by 270° on the y-axis for each consecutive pulling rate, and two states at 170° and -35° indicating local collagen and R-helix structures are clearly seen in each case. Flips occur on similar intervals between 5 and 20 A˚ with respect to the shortening of the restraint. The two bottom traces show Φ, and ξ (see text). The dashed red lines indicate (1) conformational change to exo conformation induced by Ψ flip from collagen to R-helix range, (2) back to endo conformation because of mechanical stress during shortening the restraint, and (3) second conformational change to exo configuration accompanied by flip of Ψ. C Steric view of the HYP-GLY unit for collagen-like (left) and R-helix fold (right). The torsion angle ψ switches from more than 130° to less than 0° (dashed line). During the flip, the positions of the HYP backbone oxygen and the GLY nitrogen are exchanged (green arrows). The R-helical structure after the flip is stabilized by a hydrogen bond of the N-H of glycine to the ring nitrogen of hydroxyproline.

the force-distance curves in this example are highly consistent, indicating that the hydrogen bonding to the solvent has no significant influence on the force-distance curves. Langmuir 2010, 26(19), 15248–15256

Force-distance curves from classical MD calculations had been applied before for determining rupture forces from surfaces. Here, consistent adhesion energies were evaluated from the areas of rupture peaks. Thereby, adsorption reactions induced by electrostatic interaction can be quantified by force field simulations. The results depend sensitively on the parametrization of surface charges. Variation of the charges for surface hydroxyls in a feasible range results in energies ranging from 40 to 190 and 60 to 170 kJ/mol for the contacts of the glutamic acid and lysine side chains to the surface, respectively. These values are similar to 110 and 160 kJ/mol from a recent first principles calculation.23 The low energies from the force field calculation with chemically feasible hydroxyl charges might indicate that the first principle method includes significant electronic contributions to the energies, which are not taken into account in a force field simulation. The use of high formal charges for the surface hydroxyls would then compensate for this deficiency by overestimating the purely electrostatic interaction. This incites more effort on the parametrization of the TiO2 substrate beyond the stable description of the bulk crystal only. In principle, data sets from both classical and first principles molecular dynamics are subject to some uncertainty. In the present work, the energies crucially depend on the partial charges, whereas in DFT calculations sampling over a large number of initial configurations is very costly. Thus, the overall consistence of force field and CPMD data shows that both recently developed methods of extracting adhesion energies on surfaces in aqueous solution from molecular dynamics runs yield valid results. In general, the desorption energy found by induced desorption is higher than the bond energy between the adsorbed molecule and surface, since the adsorbate takes up additional elastic energy. This depends on the configuration and especially on the CR atom, on which the restraint is applied. The effect is an artifact of the method but at least coincides qualitatively with the theory of dissociation reactions. In unrestrained large molecules, redistribution of energy among internal degrees of freedom will result in a higher apparent desorption energy and slower kinetics than those in small systems with the same bond strength to the surface. In all runs, stable contact points of the decapeptide at the surfaces were found. The electrostatic interaction with the polar TiO2 surface (about 0.2-3 nN) is by about a factor of 10 higher than that with hydrophobic systems as studied in ref 12. Consequently, consistent force-distance curves were obtained at rather high pulling rates in the range of 0.1-2 m/s. Pulling on single atoms by restraints has been shown earlier to induce denaturation under certain conditions. Flips of the ψ angle as induced by structural deformation due to restraints may be considered as elementary steps for folding or unfolding. DOI: 10.1021/la101687m

15255

Article

Hydroxyproline is shown to be more susceptible to switch from the collagen to a local helix structure than proline. This effect might contribute to the differences in the properties of hydroxylated POG and nonhydroxylated PPG collagens. Acknowledgment. We gratefully acknowledge assistance of Stefan Horn, Malte Passvogel, and Helge Grosch (Universit€at

15256 DOI: 10.1021/la101687m

K€ oppen and Langel

Greifswald) during some of the classical calculations, Julian Schneider (Universit€at Bremen) for evaluating partial charges of the surface hydroxyl groups, Prof. Lucio Colombi Ciacchi, Mohammad Koleini (Universit€at Bremen), Bastian Ohler, and Armin Marx (Universit€at Greifswald) for very useful discussions, as well as financial support by the DFG (LA700/8-1 and KO3811/ 1-1) and the COST D41 action.

Langmuir 2010, 26(19), 15248–15256