Predicting the Structure–Activity Relationship of Hydroxyapatite

Article pubs.acs.org/Langmuir

Predicting the Structure−Activity Relationship of HydroxyapatiteBinding Peptides by Enhanced-Sampling Molecular Simulation Weilong Zhao,† Zhijun Xu,† Qiang Cui,∥ and Nita Sahai*,†,‡,§ Department of Polymer Science, ‡Department of Geosciences, and §Integrated Bioscience Program, University of Akron, Akron, Ohio 44325-3909, United States ∥ Department of Chemistry and Theoretical Chemistry Institute, University of WisconsinMadison, Madison, Wisconsin 53706-1322, United States Downloaded via DURHAM UNIV on January 3, 2019 at 04:55:10 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: Understanding the molecular structural and energetic basis of the interactions between peptides and inorganic surfaces is critical to their applications in tissue engineering and biomimetic material synthesis. Despite recent experimental progresses in the identification and functionalization of hydroxyapatite (HAP)-binding peptides, the molecular mechanisms of their interactions with HAP surfaces are yet to be explored. In particular, the traditional method of molecular dynamics (MD) simulation suffers from insufficient sampling at the peptide−inorganic interface that renders the molecularlevel observation dubious. Here we demonstrate that an integrated approach combining bioinformatics, MD, and metadynamics provides a powerful tool for investigating the structure−activity relationship of HAP-binding peptides. Four low charge density peptides, previously identified by phage display, have been considered. As revealed by bioinformatics and MD, the binding conformation of the peptides is controlled by both the sequence and the amino acid composition. It was found that formation of hydrogen bonds between lysine residue and phosphate ions on the surface dictates the binding of positively charged peptide to HAP. The binding affinities of the peptides to the surface are estimated by free energy calculation using parallel-tempering metadynamics, and the results compare favorably to measurements reported in previous experimental studies. The calculation suggests that the charge density of the peptide primarily controls the binding affinity to the surface, while the backbone secondary structure that may restrain side chain orientation toward the surface plays a minor role. We also report that the application of enhanced-sampling metadynamics effects a major advantage over the steered MD method by significantly improving the reliability of binding free energy calculation. In general, our novel integration of diverse sampling techniques should contribute to the rational design of surface-recognition peptides in biomedical applications. regeneration of damaged tissue.6,7 Compared to recombinant proteins, peptides have several significant advantages as the functionalizing biological factors, including being less expensive, easier to purify, and less vulnerable to denaturation.8 Regarding connective tissue applications, peptide sequences mimicking natural proteins in the ECM of bone and teeth, such as fibronectin, statherin, osteopontin, and bone sialoprotein, are viable candidates of functionalizing motifs.9−12 Alternatively, a combinatorial screening method, such as phage display, has been employed to identify peptides that have specific binding to pure HAP or bone-like apatite mineral.13−15 By combining these peptides with other signaling molecules, their preferential adsorption to HAP can be exploited to design the material’s control over cell adhesion, proliferation, and differentiation.14−18 However, despite the discovery of a variety of natural and genetic-engineered HAP-binding peptides, the

1. INTRODUCTION Regeneration and restoration of the structure and function of pathological connective tissues, such as bone, teeth, and cartilage, require the biomaterials to elicit the formation of extracellular matrix (ECM). The specific cellular responses that lead to ECM formation in vivo can be manipulated by establishing cell-biomaterial communication at the molecular level.1 For example, surface modification of implantable materials by biomimetic soluble signals is a promising approach to enable the implants to gain a high level of control over the cellular function and host tissue integration.2,3 Calcium phosphate (CaP) ceramics have been widely used in the applications of dental and orthopedic implants, as CaP is the major inorganic component of teeth and bone and displays notable biocompatibility with the target tissue.4 The prominent CaP polymorph formed during odontogenesis or osteogenesis is hydroxyapatite (HAP, ideal stoichiometric form Ca5(PO4)3(OH)).5 Functionalizing the surface of HAP with biological factors will enhance its bioactivity associated with HAP-cell and HAP-ECM recognition, thus promoting the © 2016 American Chemical Society

Received: April 25, 2016 Revised: June 10, 2016 Published: June 21, 2016 7009

DOI: 10.1021/acs.langmuir.6b01582 Langmuir 2016, 32, 7009−7022

Article

Langmuir Table 1. Sequence and Physicochemical Properties of the Four Peptides Studied Here

From ref 27, determined as the equilibrium adsorption amount on biomimetic apatite thin film. bBasic, acidic, and neutral polar residues are heighted in blue, red, and green text, respectively. cReference 14. dDe novo sequence from this study. eReference 27.

a

focus on the molecular structure of binding, whereas the underlying energetics, such as the free energy profile of HAP− peptide interactions, remained marginally understood. To complement classical MD simulations, Gray and co-workers have developed a novel, bioinformatics-based molecular simulation tool, RosettaSurface, for modeling protein−HAP interactions.41,42 Combined with constraints from ssNMR data, RosettaSurface was successfully employed to obtain a highresolution (∼3 Å) picture of the entire statherin and a portion of amelogenin bound to HAP.43,44 While providing energyminimized structures, RosettaSurface calculates system energy (“score”) by relying on existing crystal structure of reference sequences. Hence, only relative interaction energies are obtained for conformations corresponding to a particular sequence, which makes it difficult to compare binding free energies between different peptide sequences. In order to link the molecular-level picture with energetic aspect of HAP−peptide interaction, it is imperative to assemble a simulation pipeline that leads to appropriate extrapolation to experimentally measured binding affinity. In essence, the theoretical characterizations of this value in previous efforts were based on calculating the potential of mean force (PMF) as a function of the binding molecule’s vertical distance to the surface of the inorganic substrate.38,45,46 By comparing the difference in the PMF profile between the bound and unbound states, the binding affinity can in principle be recovered. Yet in practice, due to the strong adsorption of the molecule at the bound state, the sampling efficiency of the binding−unbinding transition may be compromised, bringing about large statistical error of PMF at both end states. Recent development of several advanced biomolecular sampling approaches, which are centered on the application of metadynamics simulation, has made the reliable calculation of PMF of biomolecular adsorption technically viable.47−49 However, the improvement of the sampling efficiency by metadynamics does require a careful consideration on the issues of statistical convergence and computational cost, which become the bottleneck in the case of sampling the multidimension degree of freedoms implicated in the interaction between large biomolecules and charged surfaces.50 In particular, this problem has not been adequately addressed for computing the binding of peptides to inorganic crystals; this motivated us to integrate metadynamics with the enhanced-sampling simulation protocol discussed below. To illustrate the protocol, the molecular binding mechanisms of several 12-mer peptides identified by phage display14,27 (Table 1) at the HAP (100) surface were studied in this work. In contrast to previous experimental and computational studies

molecular level and energetic basis of the binding, in terms of the structure−activity relationships of the peptides, remains elusive. A critical process of normal odontogenesis or osteogenesis is the mineralization of ECM through the crystallization of HAP.5,19 Insoluble collagen and soluble noncollagenous proteins (NCPs) compose ECM, and their interactions with inorganic ions or other CaP polymorphs regulate the nucleation and growth of the final HAP phase.20−23 In particular, it has been shown that proteins and peptides can mediate the transition of amorphous and crystalline CaP polymorph precursors to HAP in vitro.24−26 Hence, peptides bearing a high binding affinity to HAP, when introduced in vivo, have the potential to control the crystallization of HAP through biomimetic pathway and further improve osteo- or odontotissue regeneration.4,27 When functionalized with fluorescent molecular probes, the HAP-binding specificity of the peptides can also be exploited to monitor the mineralization process in situ, leading to a variety of diagnostic and therapeutic applications.28 Therefore, a comprehensive understanding of molecular structures at the inorganic−organic interface is essential to the rationalization of designing HAP-recognition peptides. Useful insights into the structure−activity relationships at the molecular interface between HAP and proteins/peptides have been obtained by means of solid-state nuclear magnetic resonance (ssNMR), atomic force microscopy (AFM), quartz crystal microbalance (QCM), and molecular simulation.29−33 Taking advantage of advanced ssNMR techniques, several investigators have resolved the structural basis of natural proteins involved in odontogenesis, including statherin and amelogenin, in atomistic details.11,34,35 In these studies, acidic and basic amino acid residues are shown to be important for surface binding. However, the energetic basis of the peptide− HAP interaction, with respect to relative affinities, cannot be determined by ssNMR. Conversely, AFM and QCM have provided valuable data on the binding affinity of peptides to HAP surface,36,37 although these techniques do not provide information about the molecular level structures. Molecular simulation presents a promising technique to connect the molecular-level structures with energetics on HAP−biomolecule interaction.32 Molecular dynamics (MD) has been widely used to study the adsorption of biomolecules onto HAP.38−40 While yielding helpful information about the interface, the correct interpretation of MD simulation heavily depends on factors such as accuracy of the force fields, sufficiency of the sampled conformational space, and truthfulness of the surface models. More importantly, previous MD studies tended to 7010

DOI: 10.1021/acs.langmuir.6b01582 Langmuir 2016, 32, 7009−7022

Article

Langmuir

2.2. Simulation Details. 2.2.1. RosettaSurface. The sequence of peptides of interest and their properties are listed in Table 1. The initial extended chains of peptides were built from the primary sequence. The HAP (100) surface (Figure 1a) was prepared as

focusing on natural sequences in NCPs, we selected these genetic-engineered peptides, which have few acidic residues and bear low charge density.9 Intriguingly, these peptides were observed to bind to HAP with similar affinities compared with natural sequences.18,27,29 Understanding their binding mechanisms to HAP is significant for exploring their potentials as HAP surface-functionalizing agents and mineralization probes. In particular, our aim is to elucidate the relationships between the primary sequence and the charge distribution with the binding conformations and affinities of the peptides. To this end, we used the sequences initially discovered by Segvich et al., namely, peptide VTKHLNQISQSY (VTK) and its scrambled variant, IYQSKHTLSNQY (VTK_s).14 These two peptides have similar binding affinities to HAP,14,27 while our preliminary RosettaSurface calculations suggested that they adopted helix and random coil conformations on the HAP (100) surface, respectively. To examine the effect of charge distribution, the neutral glutamine (Gln) residue on the seventh position of VTK was mutated to a charged glutamate (Glu) residue to create a new peptide, VTKHLNEISQSY (VTK_7E, Table 1). This point mutation was designed to provide an extra binding site to the surface. The fourth peptide was built by phosphorylating VTK (pVTK, Table 1) to further elucidate the effect of primary amino acid composition on the binding conformation and affinity.27 The (100) surface of HAP was chosen as the substrate because it is the most dominant face identified in both biogenic crystals and synthetic samples.51,52 RosettaSurface was combined with classical MD to establish the most stable binding structures of peptides at the HAP surface, from which the molecular details of HAP−peptide interactions were identified. To address the sampling challenge for computing binding PMF, parallel-tempering metadynamics was applied to calculate the binding free energies of the peptides to the HAP surface. We demonstrated that this method effected a significant improvement on free energy estimation over the widely used steered MD approach. The enhanced sampling of binding process by metadynamics corroborated the reliability of the simulation protocol in identifying the molecular pathway and the energetic basis of surface−peptide interaction. Results suggested that the binding affinities of the peptides are determined predominantly by the density and the location of charged amino acids while the binding conformations have a minor effect. The computed affinity difference between the various peptides is consistent with experimental measurement of peptide adsorption isotherms.27 Based on the microscopic and energetic data, the structure−activity relationship of the HAP-binding peptides and its implications in biological mineralization are discussed.

Figure 1. Simulation setup: (a) Side and top views of the HAP (100) surface. Color designation: red, oxygen; tan, phosphate; cyan, calcium; white, hydrogen. (b) Schematic of simulation pipeline. The lowest energy conformation of HAP−peptide was extracted from RosettaSurface and transferred to MD. Water molecules were added, followed by 40 ns MD simulation. Then, PT-Meta-WTE and SMD were performed on the configurations taken from the production run of the preceding MD step. The free energy profile along the z-direction, ΔG(z), is calculated using both methods to evaluate the computationally predicted binding affinity. described earlier38 and was used throughout the present work. The simulation was performed using the RosettaSurface module in Rosetta3.4 software, with an implicit solvation model.41 40 000 structures (also called “decoys”) were generated for both the solution and the adsorbed states for each 12-mer peptide. Structure prediction was performed without any preset conformational constraint, such as experimentally measured backbone dihedral angle or hydrogen bond information. Energy values generated by RosettaSurface were ranked to find the lowest energy conformation in the adsorbed state, which was used subsequently as the initial configuration of MD simulation. For the structural analysis, 1000 lowest energy conformations were extracted, and the secondary structural properties were calculated by the Dictionary of Secondary Structure of Protein (DSSP) software.58 2.2.2. Classical MD. Simulations were performed using GROMACS 4.5.5.59 The configuration of HAP−peptide was obtained from RosettaSurface as described above, and then TIP3P water molecules were filled into the simulation box by the GROMACS genbox tool. Proper numbers of sodium or chloride ions were added not only to neutralize the system but also to set the ionic strength of bulk solution phase to 0.1 M. Energy minimization of the configurations was carried out using the steepest descent algorithm. The systems were then equilibrated under the NVT ensemble for 100 ps, followed by 500 ps NPT simulation. After equilibration, 40 ns MD production run was performed under the NVT ensemble. For NVT and NPT simulations,

2. METHODS 2.1. Force Fields. The Lennard-Jones interatomic interactions, the bond, and the angle potentials for HAP were derived from Hauptmann et al.,53 as described in our previous publications.38 For water molecules, the TIP3P model was used.54 CHARMM27 force fields55,56 were employed to describe the intramolecular and intermolecular interactions of the peptides. The Lorentz−Berthelot mixing rule was used to obtain the cross-terms of interatomic interactions. In our previous work, the properties of HAP−water−organic interface obtained using the current force fields have been benchmarked against experimental data and density functional theory (DFT) calculations with a combination of explicit and implicit solvent models.38,57 Therefore, the accuracy of our force field scheme is likely adequate for the purpose of understanding relative binding affinities of the peptides in Table 1. 7011

DOI: 10.1021/acs.langmuir.6b01582 Langmuir 2016, 32, 7009−7022

Article

Langmuir

Figure 2. Structural analysis of results from RosettaSurface in which implicit solvent is used: (a) Configurations of global energy minima. Peptide backbones are represented in blue cartoon. Amino acid side chains that are closer than 8 Å to the surface are highlighted using licorice representation. Color designation: silver, carbon; red, oxygen; blue, nitrogen; tan, phosphate; cyan, calcium; white, hydrogen. (b) Percentage of amino acid residues belonging to different peptide secondary structrures. RC: random coil. (c) Root-mean-square-fluctuation (RMSF) of molecular conformation plotted against residue index number. Nosè−Hoover temperature coupling60,61 and Parrinello−Rahman pressure coupling62 were used, respectively. In particular, semiisotropic pressure coupling was applied in the NPT ensemble to obtain realistic density of water in the interface and bulk regions. The particle mesh Ewald (PME) method was used to treat long-range electrostatic interactions.63 14 Å was used as the real space cutoff in the PME calculations. A switch radius of 12 Å was used for the Lennard-Jones potential. Initial velocity of particles in the system was assigned randomly on the basis of Maxwell distribution. All the data analysis was performed using GROMACS built-in utility. Molecular configurations were visualized using VMD 1.9.1.64 2.2.3. Parallel-Tempering Metadynamics. Well-tempered metadynamics was performed using PLUMED 1.3 patched to GROMACS 4.5.5.65 Two reaction coordinates were used for each peptide, including one collective variable (CV) representing the position of the center-of-mass (COM) of the peptide in the direction perpendicular to the (100) surface (z-direction). To enhance the sampling of binding−unbinding events, which involve slow degrees of freedom orthogonal to the COM separation of the peptide with the surface, we introduced another CV that overcomes the hidden free energy barriers. For VTK, VTK_s and VTK_7E, the second CV is the position of the COM of lysine residue (Lys) in z-direction. For pVTK, Lys is not the primary binding residue; hence, the position of one SerP in z-direction was assigned to the second CV. The introduction of a second distance CV is expected to augment the sampling of peptide− surface separation, while also enabling the exploration of conformational space related to different peptide binding orientations. The width of Gaussian hills was set to 0.5 Å. The height of the hills was set initially at 4 kJ mol−1, with a well-tempered factor of 10. The Gaussian hills were added to the simulation every 500 fs. Parallel-tempering (PT) method was combined with metadynamics to amplify the sampling efficiency of both the binding−unbinding process and the conformational exploration of the peptide (Supporting Information, part 1).66 According to the size of our system, a large number of replicas (∼100) are sufficient to achieve satisfactory exchange probability (>20%). In this case, PT-Meta in Well-tempered

Ensemble (PT-Meta-WTE) was introduced to reduce the computational demand.48,49 The WTE was prepared by 1 ns PT-Meta using system’s potential energy as the sole CV. The biasing potentials accumulated at this stage were recorded and then applied as static bias in the subsequent PT-Meta using the two position CVs (Supporting Information, part 1). The effect of WTE was significant, as demonstrated by increased fluctuation of the potential energy and round-trip speed of PT-Meta (Figure S1). Because of the WTE application, only 24 replicas were required to achieve an exchange probability of 45%, covering the temperature ladder from 300 to 480 K. The exchange between adjacent replicas were attempted every 1 ps. The PT-Meta-WTE simulation was ran to 50 ns for VTK and VTK_7E, 45 ns for VTK_s, and 40 ns for pVTK, as determined by checking the convergence of free energy surface (see below). 2.2.4. Free Energy Calculation. The procedure to recover unbiased free energy surface from metadynamics have been described extensively in the literature67−69 and will only be stated briefly here. The reaction coordinates of 2-D free energy surface were obtained by subtracting the position of the solid−water interface from the position CVs of peptide COM. Based on the unbiased 2-D free energy surface, the 1-D free energy profile as a function of peptide COM-surface separation distance (PSD) was obtained by integrating the 2-D surface over the other reaction coordinate. The position of the surface was calculated as the time average of the topmost Ca2+ coordinates in the z-direction. The boundary between bound and free states of the peptide, zdiv, was set at PSD = 2.85 nm. The binding free energy, ΔGbinding, was then calculated according to the ratio of probability density between the states:

Pbound =

Pfree = 7012

zdiv

1 − z min

1 z max − zdiv

∫z

∫z

zdiv

z max div

exp( − β ΔG(z)) dz

min

exp( − β ΔG(z)) dz

(1)

DOI: 10.1021/acs.langmuir.6b01582 Langmuir 2016, 32, 7009−7022

Article

Langmuir

Figure 3. (a) Equilibrium conformations at the end of 40 ns explicit solvent MD simulation (orange backbone, water molecules not shown for clarity) compared to the lowest energy conformations predicted by RosettaSurface with implicit solvent (blue backbone). Color designation is the same as in Figure 2a. Insets: conformation of peptides using surface representation, in which basic, acidic, neutral polar, and nonpolar residues are colored by blue, red, green, and white, respectively. (b) Distribution of residues belonging to helix, turn, or RC secondary structures for 4000 conformations extracted from 40 ns simulation of each peptide. (c) Evolution of backbone RMSD during 40 ns simulation. (d) Probability distribution of PSD counted over the entire 40 ns.

ΔG binding = −

1 ⎛ Pbound ⎞ ln⎜ ⎟ β ⎝ Pfree ⎠

To minimize the systematic bias of free energy coming from the hysteresis in pulling simulation, the weighted histogram analysis method (WHAM)72 was adapted to compute the binding free energy in conjunction with Jarzynski’s equality (Supporting Information, part 2).

(2)

where β is the Boltzmann temperature factor at 300 K and zmin and zmax are the lower and upper boundaries of the explored phase space of peptide binding/unbinding. To assess the quality of free energy calculation, ΔGbinding was plotted as a function of PT-Meta-WTE simulation time (Figure S2). The result suggests when simulation was extended beyond 40 ns, sampling of bound and free states is sufficient to render a converged binding free energy. The error of ΔGbinding was estimated by the method proposed by Laio et al.70 2.2.5. Steered MD (SMD) Simulation. For each peptide in the adsorbed state, 30 configurations were extracted from the last 20 ns production run of classical MD. These configurations were used to perform 30 repeated SMD simulations independently, by employing the pulling module of GROMACS. During the simulation, the COM of the peptide was forced to move along with a reference point, using a harmonic restraint, ut:

ut = k(zcom(t ) − zcom(0) − vt )2

3. RESULTS 3.1. Structural Analysis Based on RosettaSurface. The lowest energy conformations of the target peptides adsorbed on HAP (100) surface are shown in Figure 2a. Among the four peptides, VTK and VTK_7E exhibit predominantly helical backbone structure, while for VTK_s and pVTK it is random coil (RC). It is clear that alteration of both the sequence and primary amino acid composition has significant effect on the binding structure of the peptides. Statistical analysis on the backbone conformation of the 1000 lowest energy structures extracted from RosettaSurface results corroborates the same observation (Figure 2b). Helical backbone conformation takes 35% of all the 1000 structures for peptides VTK and VTK_7E. In comparison, only 8% and 14% of the structures can be assigned to helix, for VTK_s and pVTK, respectively. For the latter two peptides, RC conformation is dominant and takes 48% and 40% of all the 1000 structures, respectively. The distinction between the backbone conformations of VTK, VTK_s, and pVTK is consistent with the previous results of RosettaSurface simulation.27 We also calculated the root-meansquare-fluctuation (RMSF) of the individual residues in the pool of 1000 lowest energy structures of the four peptides (Figure 2c). VTK_s and pVTK peptides have larger RMSF compared to VTK and VTK_7E. This result is in line with both

(3)

where the spring constant k is assigned to 2000 kJ mol−1 nm−2, zcom(t) is the COM of peptide in z-direction at simulation time t, zcom(0) is the COM of peptide in z-direction at t = 0, and v is the velocity of the reference point. In all the pulling simulations, v was set at 0.2 nm ns−1. At t = 0, the reference point was moved away in the direction normal to the surface. In order to allow complete detachment of the peptide from the surface and to enter the bulk phase of water, a pulling time of 10 ns was required. For the calculation of binding free energy, the values of PSD for all 30 simulations were collected as a function of time. Jarzynski’s equality71 was applied to calculate the free energy difference between bound and unbound states of the peptide during the nonequilibrium pulling process (Supporting Information, part 2). 7013

DOI: 10.1021/acs.langmuir.6b01582 Langmuir 2016, 32, 7009−7022

Article

Langmuir

mean-square deviation (RMSD) of peptide backbone conformation as a function of time (Figure 3c) and the probability distribution of PSD (Figure 3d). For peptide VTK, the backbone orientation of the surfacebound conformation changes significantly compared to the corresponding initial RosettaSurface structure, while still retaining the helical secondary structure (Figure 3a). The rise of backbone RMSD starting at 5 ns (Figure 3c) and the large fluctuation in PSD (Figure 3d) further confirm such a change. Because there is only one charged residue on VTK (Lys) to provide stable binding to the surface, the peptide alternated between the flat (blue) and the inclined (orange) orientations, causing the peptide to “swing” on the surface, with Lys acting like the “hinge”. The observation of this phenomenon suggests that polar side chains, such as Tyr and His, which appeared to bind closely to the surface in RosettaSurface predicted conformations, are weakly binding residues. Although we infer that the interactions of the latter residues with surface ions are challenged by more favorable peptide−solvent interactions, we stress that potential drawback of the fixed point-charge model of the electron-conjugated residues and the surface applies to the MD simulation results as well. Concretely, it has been argued that polarization effect of the charges should be explicitly included in the force field in order to describe intermolecular interactions accurately.74,75 A more comprehensive force field parametrization is thus desirable to address such polarization effect at the interface. For VTK_7E, the blue and the orange backbones overlap with each other (Figure 3c). The backbone RMSD of VTK_7E stays below 1.2 Å, and the PSD does not change significantly in 40 ns (Figure 3c,d). These results suggest that VTK_7E binds to the surface robustly, which may be ascribed to the strong electrostatic interactions between Glu and Lys, near the C- and N-termini of the peptide, respectively, with calcium and phosphate ions. Thus, it appears that the two “surfacerecognition” sites here result in a global energy minimum of peptide binding conformation, rendering the ability of this designed peptide to maintain its secondary structure even in the presence of explicit water adsorption on the surface, while holding high surface-binding affinity. The peptides VTK_s and pVTK, which are predicted by RosettaSurface to have RC backbone conformation, maintain their structural features during MD simulation (Figure 3a,b). Their backbone RMSD values increase to 5 Å, and the mean values at equilibration are around 3 and 3.5 Å respectively for VTK_s and pVTK (Figure 3c). The width of distribution of PSD is also larger than that of VTK_7E (Figure 3d), which is expected typically for a RC peptide. Notably, despite the structural fluctuations of the backbone, the residues that interact with the surface through 40 ns simulation remain the same, namely, Asn and Lys for VTK_s and SerP for pVTK. It is reasonable that the SerP residue bearing a −2 charge becomes a stronger binder than Lys (+1 charge). However, the neutral polar residue, Asn on VTK_s, continues to bind to the surface, in sharp contrast to the situation for VTK. This observation may be attributed to the inherent higher flexibility of RC peptides compared to helical ones. Because of such flexibility, the hydrogen bonds on the backbone only have weak constraint on the conformation dynamics of the surface-bound peptide, which increases the chance of forming multiple bonds to the surface. Generally, the multisite bound state is more easily accessible to RC peptide than to helical peptide.49 Therefore, the binding conformations of RC peptides VTK_s and pVTK

the conformation visualization and statistics (Figure 2a,b), which can be explained by higher ratio of RC conformation that usually relates to more structural fluctuation. Apart from adsorbed backbone structure, we also determined the possible binding amino acid residues on the four peptides, by setting a cutoff value of 8 Å for side chain-to-surface distance. These residues are highlighted by licorice representation in Figure 2a. The features at atomic detail are also comparable to the RosettaSurface low-energy structures reported previously.27 In particular, for VTK, VTK_s, and VTK_7E, the lysine (Lys) residue positions itself within 3 Å to the HAP surface. For VTK_7E, Glu also interacts with the surface (3 Å), which is expected because of its affinity to surface calcium ions. As a result of electrostatic attraction to calcium ions, SerP on peptide pVTK also has a smaller distance (2.5 Å) to the surface. Interestingly, neutral but polar residues, such as asparagine (Asn), histidine (His), and tyrosine (Tyr), tend to bind to the surface as well (from 3 to 5 Å). This indicates that their interactions with surface ions are favored over the interactions with the continuum solvent, which lowers the total potential energies calculated by RosettaSurface. The mechanism for this ability of neutral residues to bind to HAP remains ambiguous, based on previous experimental studies.15,73 We also note that the fixed charge force field used here may be insufficient for treating the interactions between cations and aromatic groups such as the phenol ring in Tyr due to the missing polarization contributions.74,75 Therefore, it is critical to check whether such observation is reproducible in classical MD using explicit solvent. More important, the relative affinity of different amino acid side chains to the surface should be quantified from a thermodynamic perspective by calculating the free energies of adsorption. These two questions merit further investigation in the remaining parts of this paper. 3.2. Binding Conformation Stability in Explicit Solvent. RosettaSurface provides a good starting point to understand the adsorbed conformations of peptides at the surface; however, due to the algorithm design, explicit water molecules and ions cannot be included, such that the simulation may not capture the realistic experimental conditions. To reproduce atomic level interactions at the explicitly solvated interface, additional MD simulation is necessary. Using the lowest energy configurations from RosettaSurface (as shown in Figure 2a) as the initial structures, we performed 40 ns classical MD simulation with explicit water and ions. The extended length of the MD simulation alleviates the often-suffered problem of insufficient sampling in peptide− surface simulations and can capture adequately the conformational dynamics of the surface-bound peptide. Figure 3a shows the comparison of structures obtained from MD (orange backbone) and the corresponding ones from RosettaSurface (blue backbone). The MD structures represent the equilibrium phase of the 40 ns trajectories. The insets in Figure 3a display the molecular surface of the bound peptides rendered from representative MD configurations. The peptides are colored according to the charge type of the residues. All the four peptides maximize the contact of charged or polar residues with the surface. At first glance, the general secondary structure features appear the same, which is also supported by the data in Figure 3b, which quantifies the distribution of secondary structures. However, the surface-bound conformations of all the peptides vary from the RosettaSurface structures, except for VTK_7E. The difference between MD structures with initial RosettaSurface conformation is also characterized by the root7014

DOI: 10.1021/acs.langmuir.6b01582 Langmuir 2016, 32, 7009−7022

Article

Langmuir

Figure 4. Force and conformations of peptides unbinding from surface as simulated by SMD. (a) The instantaneous pulling force (black) and the running average curve over 50 data points (red) plotted as a function of PSD. The individual unbinding events are indicated by the force peaks that were marked by asterisks. (b) Configuration snapshots corresponding to the marked events in force−distance curves. Residues with binding interaction to the surface are highlighted by licorice representation. Water molecules are not shown for clarity. Peptide backbone is shown as an orange ribbon, and the remaining color designations are the same as in Figures 2 and 3.

are more stable than that of helical peptide VTK, as suggested by MD simulation. In summary, for VTK_7E, both the negatively (Glu) and positively (Lys) charged residues are in contact with the surface, suggesting a stable binding conformation. Without the “anchoring” role of additional Glu, VTK peptide forms close contact with only the single Lys residue, indicating again that it is not bound to the surface in a stable conformation. Thus, it is reasonable to postulate that VTK_7E surpasses VTK in binding affinity to HAP (100) surface. For VTK_s and pVTK, multiple peptide−surface contacts are established because of their RC backbones. Therefore, it becomes intriguing, at the thermodynamic level, how multiple contacts can enhance binding. In the next section, the affinities of these peptides to the surface are compared in quantitative details to address this problem. 3.3. Molecular Basis of Binding. We applied SMD as a method to simulate the unbinding process of the peptides from the surface, mimicking AFM pulling experiment (Supporting Information, part 2). The imposed pulling force in the simulation drives the near-equilibrium unbinding of the peptides, from which the interactions of bound residues with the surface can be monitored. Figure 4a illustrates the pulling force as a function of PSD. Critical transition events are indicated by asterisks, the molecular snapshots of which are shown in Figure 4b. These events include (1) the unbinding of Lys on VTK, (2) the unbinding of Asn on VTK_s, (3) the unbinding of Lys on VTK_s, (4) the unbinding of Glu on VTK_7E, (5) the unbinding of Lys on VTK_7E, and (6) the unbinding of SerP on pVTK. The results suggest that Lys manifests as the pivot binding residue of all the peptides, except for pVTK. This observation is in line with the information obtained from RosettaSurface and classical MD. The atomic details of the interactions are shown in Figure 5. As two representative cases, Lys of VTK_s and VTK_7E both form three hydrogen bonds with the phosphate

Figure 5. Atomic details of the binding of (a) VTK_s and (b) VTK_7E to HAP (100) surface. Hydrogen bonds are shown by pink dotted lines. Interatomic interactions between oxygen and calcium ion are shown in green. The molecular surfaces of HAP and peptide are shown transparently in light yellow and white, respectively. The number which appears after the residue name indicates the residue index number in the peptide.

ions on the surface (Figure 5, pink dotted lines). For Asn of VTK_s, the carbonyl oxygen forms electrostatic interaction with surface calcium ion (Figure 5a, green dotted line), while one of the hydrogen atoms of amide forms an additional hydrogen bond with phosphate ion. Additionally, for VTK_7E, oxygen atoms on the Glu side chain form electrostatic interactions with surface calcium ions (Figure 5b, green dotted lines). The SMD pulling simulation at near-equilibrium reveals that the binding of peptides to HAP (100) surface is controlled by direct interactions of individual residues with the surface. Although the strong binding of water to the HAP (100) surface may be expected,57 interfacial water appears not to participate in the local binding configurations between the charged and polar side chains and the surface ions (Figure 5). The displacement of interfacial water may contribute to entropy 7015

DOI: 10.1021/acs.langmuir.6b01582 Langmuir 2016, 32, 7009−7022

Article

Langmuir

Figure 6. 2-D FES contoured based on the two distance CVs. The relative free energy values are indicated by color bar. From the global minima (GM) on FES (dark blue regions), the representing structures are shown as inset snapshots (black arrows). The peptide backbone is colored pink. For HAP crystal surface, only Ca2+ and PO43− ions in the vicinity of binding residues (distance