Refined Empirical Force Field to Model Protein–Self-Assembled

Article Cite This: Langmuir 2019, 35, 9622−9633

pubs.acs.org/Langmuir

Refined Empirical Force Field to Model Protein−Self-Assembled Monolayer Interactions Based on AMBER14 and GAFF Pratiti Bhadra and Shirley W. I. Siu* Department of Computer and Information Science, University of Macau, Avenida da Universidade, Taipa, Macau

Downloaded via KEAN UNIV on July 25, 2019 at 07:00:26 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Understanding protein interaction with material surfaces is important for the development of nanotechnological devices. The structures and dynamics of proteins can be studied via molecular dynamics (MD) if the protein−surface interactions can be accurately modeled. To answer this question, we computed the adsorption free energies of peptides (representing eleven different amino acids) on a hydrophobic self-assembled monolayer (CH3−SAM) and compared them to the benchmark experimental data set. Our result revealed that existing biomolecular force fields, GAFF and AMBER ff14sb, cannot reproduce the experimental peptide adsorption free energies by Wei and Latour (Langmuir, 2009, 25, 5637−5646). To obtain the improved force fields, we systematically tuned the Lennard-Jones parameters of selected amino acid sidechains and the functional group of SAM with repeated metadynamics and umbrella sampling simulations. The final parameter set has yielded a significant improvement in the free energy values with R = 0.83 and MSE = 0.65 kcal/mol. We applied the refined force field to predict the initial adsorption orientation of lysozyme on CH3−SAM. Two major orientationsface-down and face-upwere predicted. Our analysis on the protein structure, solvent accessible surface area, and binding of native ligand NAG3 suggested that lysozyme in the face-up orientation can remain active after initial adsorption. However, because of its weaker affinity (ΔΔG = 7.86 kcal/mol) for the ligand, the bioactivity of the protein is expected to reduce. Our work facilitates the use of MD for the study of protein−SAM systems. The refined force field compatible with GROMACS is available at https://cbbio.cis.um.edu. mo/software/SAMFF.

■

clease A,29 bovine serum albumin,30 ubiquitin,31 etc.) have been studied for their adsorption pathway, binding, orientation, conformational changes, and desorption from SAMs using molecular simulations. A prerequisite for meaningful results from MD simulations is the use of an accurate molecular mechanical force field. The development of an accurate empirical force field requires the availability of suitable experimental data to which simulated properties can be compared and parameters of the force field can be optimized. The first benchmark data set characterizing protein−surface interactions was produced by the Latour group using a host−guest peptide model, where the guest amino acid is in the middle of the 9-residue amino acid sequence.32 In their study, the standard state adsorption free ° ) value for each peptide−surface combination energy (ΔGads was determined by surface plasmon resonance spectroscopy. The adsorption behavior of 12 amino acid residues over 9 different functionalized alkanethiol SAMs on gold was studied experimentally. Their results indicated that the peptide adsorption strength on noncharged surfaces correlates positively and linearly to the hydrophobicity of the surface, with specific interactions between functional groups of the peptide and the SAM surfaces playing a secondary role. In total, 108 peptide−surface systems were obtained with their

INTRODUCTION Understanding the protein adsorption process on selfassembled monolayers (SAMs) is fundamental to advance the technologies of biosensors and analytical devices,1−4 drug delivery systems,5,6 and biomaterials.7−9 Protein adsorption is considered to be the first important event when a material is in contact with a biological environment.10−12 Adsorbed proteins can influence properties of the surfaces and interfere with the functionalities of the devices.13 Hence, optimization of the surface structure and chemistry of SAMs to properly control or prevent protein adsorption has been the main goal of rational surface design. While a range of experimental techniques have been developed to study protein adsorption, 14 these techniques cannot provide a full picture at the molecular level. Without a thorough understanding of the protein−SAM surface interactions with atomic detail, surface design can only be approached by experimental trial and error, which is costly and inefficient. Over the last three decades, molecular dynamics (MD) has emerged as an important tool to study biomolecular processes. In particular, it has been used to elucidate the mechanism of protein adsorption on surfaces.15 Since the success of atomistic modeling of SAMs,16,17 the number of studies on simulating SAMs and proteins on SAMs has been increased tremendously. MD simulations can provide pictorial presentation for the formation,18,19 structural arrangement,20 and dynamics of SAMs.21 Numerous proteins (such as lysozyme,22,23 fibrinogen,24,25 cytochrome c,26 lectin,27 neuromedin-B,28 ribonu© 2019 American Chemical Society

Received: May 8, 2019 Revised: June 20, 2019 Published: June 27, 2019 9622

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir experimental ΔG°ads values, making this data set indispensable for validating the adsorption behavior of peptides from simulation methods. Since then, several studies have been conducted to rigorously assess the ability of existing biomolecular force fields to reproduce experimental values of the same peptide− SAM systems in simulations. For example, Vellore et al.33 showed that the combination of CHARMM22 (for non-SAMs molecules) and CGenFF (for SAMs) was able to generate ΔG°ads within 1 kcal/mol of the experimental values for most systems, but for hydrophobic and positively charged surfaces, the deviations were as large as 4 kcal/mol, suggesting that force field modification is desired. Interestingly, despite having experimental data of 12 guest residues, their study only included peptides of 5 guest residues. Though imperfect, CHARMM22 was later found by Collier et al. to be better than OPLS-AA and AMBER94 in simulating the conformational behavior of the model LKβ7 and LKα14 peptides in the adsorbed state and in solution.34 On the other hand, Deighan and Pfaendtner35 showed that when using parallel tempering metadynamics to improve the convergence of simulations, both CHARMM22 and the newer AMBER99SB can produce low free energy conformations of LK-peptides in the adsorbed state that are more consistent with experimental observations than that of OPLS-AA. In that study, AMBER99SB simulations were found to deviate the least from the experiment in terms of the tilt angle values of the leucine side chain on hydrophobic SAMs. It is important to note that due to the lack of experimental data, most of the FF comparison studies for SAMs have been conducted using LK-peptides, which consist of two types of amino acids only. In addition to FF evaluation and comparison, the Latour group also took an initiative to tune interaction parameters of CHARMM22 for protein−SAM simulations. Biswas et al.36 proposed a new approach to present interfacial interactions, named dual-FF, and implemented it in the CHARMM simulation package. This dual-FF method used two separate sets of nonbonded force field parameters within the same simulation: one set for intraphase (liquid−liquid) and another for interphase (liquid−solid) interactions. Additionally, they reparametrized the TIP3P water model for protein−SAM simulations based on the TIP4P-FQ-derived electrostatic properties of water molecules at the water−SAM interface. However, in that study, they determined the nonbonded interaction parameters involving the carbon and hydrogen atoms of CH3 and CH2 only. Following the same approach, Abramyan et al. adjusted nonbonded parameters for the peptide−high-density polyethylene interaction37,38 and Synder et al. for the peptide−silica glass interaction.39 By tuning the energy well depth parameters in the Lennard-Jones (LJ) potentials of 10 amino acids, they obtained the peptide ΔG°ads within a 0.5 kcal/mol agreement with experiments. In addition, Wei and Knotts40 developed a new coarse-grained model of protein−SAM interaction that was parameterized against experimental adsorption free energies of multiple model peptides at different types of SAM surfaces. Here, our goal is to update the parameter set for protein− SAM interactions based on the recent AMBER protein force field ff14sb41 and the generalized force field GAFF42 for simulating proteins on hydrophobic SAM surfaces. In our previous work, we found that GAFF is a better option for modeling hydrophobic SAMs than other biomolecular FFs (CHARMM36, L-OPLS, Slipids, and GROMOS54a7), as it

yields the most agreeable results in the structural parameters and the phase transition behavior of SAMs.43 Therefore, in this work, we have evaluated the performance of AMBER ff14sb and GAFF for protein−SAM simulation using the experimental host−guest peptide system in Wei and Latour’s work.32 As simulations using the standard AMBER and GAFF parameters did not yield satisfactory results in terms of peptide adsorption free energy, we went on further to adjust the pairwise van der Waals interaction parameters between the functional group of selected amino acids and that of the SAM ligand. Furthermore, as a case study and to compare the effect introduced by the modified FF parameters, we have simulated the widely studied hen egg-white lysozyme protein on the CH3−SAM surface.

■

METHODS

Host−Guest Peptides and the SAM Surface. Our simulation setup closely followed the experimental system used by Wei and Latour.32 The SAM surface was made from an assembly of methylterminated alkanethiolates (CH3−(CH2)11−SH) attached on the Au(111) substrate. They were arranged as a 9 × 10 array in ( 3 × 3 )R 30° geometry with ∼5 Å spacing between neighboring sulfur atoms. Following the simulation protocol in our previous study,43 the SAM surface was first equilibrated for 15 ns at 298 K under the NVT condition; then, the peptide was added to the solvent above the equilibrated SAM. The zwitterionic host−guest peptide for free energy study has the sequence TGTG-X-GTGT, where T, G, and X represent threonine, glycine, and the guest residue, respectively. Eleven guest residues were investigated in this study; they can be grouped into three categories: polar (X = T, G, S, N), nonpolar (X = L, F, V, W), and charged (X = R, K, D). The generated peptide was placed on top of the SAM surface with a 5 Å surface separation distance (SSD); it is the distance between the center of mass (COM) of the peptide and the average COM of the SAM’s heavy atoms in the CH3 functional group. Then, TIP3P water44 and NaCl ions were added to reach the physiological condition of 140 mM, and counterions were also added to neutralize the system due to the charged guest residue; these settings yielded a water phase of about 27 Å thickness after equilibration. To mimic the experimental condition of an open system, a thick vacuum slab was included above the water phase. To minimize the periodic boundary effect, the total box length in the Z dimension was made to be three times the height of the water−peptide−SAM system. Then, the complete system was energy minimized in 1000 steps, followed by a 1 ns position restraint of the peptide and a 5 ns unrestrained MD simulation at 298 K in the NVT ensemble. To avoid the problem of water molecules crossing the periodic boundary to interact with the Au substrate, which is likely to happen over long simulations, an additional 5 Å layer of water was placed on top of the water phase. The water molecules were positionally restrained by a force constant of 1000 kJ mol−1 nm2 in all directions. The whole system was equilibrated again for 1 ns using conventional MD simulation, and the final equilibrated structure was used in subsequent free energy simulations. A snapshot of the protein−SAM model is shown in Figure 1. Sampling Method and Peptide Adsorption Free Energy Calculation. Previously, the experimental standard state adsorption free energy (ΔG°ads) values of the eleven host−guest peptides on the methyl-terminated SAM surface were reported in ref 32. To calculate ΔG°ads from molecular simulations, it is required to sample the conformational state of the peptide over the full range of SSD. To achieve this, we used well-tempered metadynamics (WT-MetaD) methods.46 Metadynamics47 is a technique for enhancing sampling in MD simulation. It can sample conformation of the target molecule as a function of a few selected degrees of freedom called the collective variables (CVs). In metadynamics, sampling is accelerated by a history-dependent bias potential. This bias potential can be written as a sum of Gaussians deposited along the system trajectory in the CV space to discourage the system from revisiting configurations that 9623

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir

and 10 ns of the same example peptide clearly showed a trend of convergence (see Supporting Information, Figure S2). Extending the simulation time did not yield different results (data not shown). Therefore, our tests indicate that the selected window interval, force constant, and simulation length for the US simulations are appropriate. The potential of mean force (PMF) of the peptide as a function of SSD was generated using the weighted histogram analysis (WHAM) ° was calculated employing the probability method.49 Specifically, ΔGads ratio method from the sampled distribution of states using the expression explained in ref 33. This method was designed to match closely with the free energy calculation method based on the ratio of concentrations of the peptide on the surface and in the bulk; thus, the computed free energy values can be directly compared to the experiment. The free energy of adsorption is defined as follows ÄÅ É ÅÅ W N ÑÑÑ ÅÅ ÑÑ ° = − RT ÅÅ Ñ ΔGads ÅÅ δP ∑ Pi ÑÑÑ (3) ÅÅÇ b i = 1 ÑÑÖ where Pi and Pb are the probabilities of the peptide being at the interfacial region and in the bulk, respectively. The region where protein−surface interaction becomes negligibly small was considered the bulk, which is typically beyond 20 Å for our systems (see Results and Discussion). Thus, Pb was obtained as the average probability of the peptide occurrences between 21 and 23 Å to the surface. N is the number of bins over the SSD-coordinate space that is used for Pi calculation, where Pi < Pb. δ is the theoretical thickness of the adsorbed layer as given in ref 32 (see Supporting Information, Table S2). W is the bin width used to produce the probability distribution from simulation, and a value of 0.13 Å was selected based on the WHAM results. R is the ideal gas constant, and T is the absolute temperature, which was 298 K. We conducted three independent sets of US simulations, and each replica was started with a different initial structure extracted from the WT-MetaD trajectory. The WHAM analysis was performed using the GROMACS tool gmx wham. The probability distributions of Pi and Pb were computed, and then the PMF profiles were generated using 200 bins. Taken the three independent sets of simulations, we computed the free energy means and the 95% confidence interval (CI) for each host−guest system as follows

Figure 1. Molecular model of the host−guest peptide (TGTG-XGTGT, where X = lysine in this image) at the CH3−SAM surface. The constrained water molecules were drawn as spheres, mobile water as lines, ions as balls, peptide as the surface, SAM as sticks, and the substrate Au atoms as spheres. SSD stands for surface separation distance (see the text). The image was prepared using VMD-1.9.3.45 have already been sampled. The WT-MetaD is a recently improved variant of the MetaD technique, in which the deposition rate of the bias potential decreases over the simulation time. This is achieved by rescaling the Gaussian initial height ω0 according to

W = ω0τG exp(− VG(S , t )/kBΔT )

(1)

where τG is the Gaussian deposition stride, ΔT is the temperature difference, and VG(S,t) is the bias potential accumulated in S over time t, where S denotes the user-defined CV. In this study, 75 ns long WT-MetaD simulations were performed using the SSD as onedimensional CV. The initial height of the Gaussian bias potential is 0.2 kJ/mol (ω0), and the deposition stride is 500 steps (τG), which is 1 ps in the simulation using the 2 fs timestep. The widths of the Gaussian potential are different for different host−guest peptides. For a peptide, the width is equal to half of the natural fluctuation of the SSD at equilibrium (see Supporting Information, Table S1). The bias factor γ = (T + ΔT)/T of the simulation is set to be 6. The simulation temperature (T) is 298 K, whereas ΔT was calculated from the given bias factor during simulation. The umbrella sampling (US) method48 was performed to compute the peptide adsorption free energies. In the US method, a bias potential is applied to drive a system from one thermodynamic state to another along a reaction coordinate (e.g., from an adsorption state to a desorption state). Here, the SSD is used as the reaction coordinate, and a harmonic potential is applied in the windows simulation Vn(S) =

k (S − Sn)2 2

CI = t * ×

σ /n

(4)

where n = 3 and t* = 4.303, σ is the variance. Tuning of the Interfacial Force Field Parameters for Protein−SAM Interaction. The adsorption behavior of protein on surfaces is a physical process that is influenced by the nonbonded interactions between protein and surface atoms. Empirical force fields model these interactions by two separate components: the electrostatic and van der Waals interactions. The former is commonly described by the Coulombic potential and the partial charges of atoms, and the latter is described by the LJ potential with the zeropotential interatomic distance and the depth of the potential well. The atomic partial charges are typically determined from quantum mechanical calculations, while the LJ parameters are determined either from quantum mechanical calculations or adjusted to match experimental data.37,50,51 We note that the LJ terms are particularly subject to errors because the interaction potential of an atom pair in a classical force field is often estimated by applying combination rules (e.g., arithmetic mean for radii and geometric mean for well-depth) on atom-type specific LJ parameters of the two atoms. Because standard LJ parameters are derived for molecules in aqueous solution, they may not be appropriate to represent the interfacial (solid−liquid) interactions between protein molecules and SAMs.52 Furthermore, Abramyan et al. has shown that LJ interaction between the peptide and surface plays a key role in peptide adsorption in terms of energy contribution.37 In AMBER, the LJ interaction is described using the 12-6 LJ potential

(2)

where S is the current SSD, Vn(S) is the bias potential, Sn is the reference SSD distance for the n-th window, and k is the spring constant, which is set to 1000 kJ mol−1 nm−2. The initial structures of the window simulations were extracted from the WT-MetaD trajectory. Using a window spacing of 1 Å, we obtained 24 US windows in total; all of them were simulated for 10 ns. The choices of the US simulation setup and parameters were made to ensure that sampling was sufficient and the generated profiles were converged. As shown in Supporting Information Figure S1, the histogram from window simulations of an example peptide TGTG-F-GTGT using the current parameters showed good overlaps between windows over the entire SSD range. Furthermore, the PMF profiles computed from 6, 8, 9624

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

ÄÅ É ÅÅi y12 i y6ÑÑÑ ÅÅjj σij zz σij zz ÑÑ j j V (rij) = 4ϵijÅÅÅÅjjj zzz − jjj zzz ÑÑÑÑ j rij z ÑÑ ÅÅj rij z ÅÅÇk { k { ÑÖÑ

Article

Langmuir

Within this simulation time, we observed that the protein orientation reached a stable state for initial adsorption, and the simulation was considered converged. Then, by pulling the equilibrated protein away from the SAM surface along Z, we obtained separated starting structures for twenty windows of US simulations. The window interval was 0.1 nm; the pulling was performed using a rate constant of 0.01 nm ps−1 and a force constant of 1000 kJ mol−1 nm2. Each US window was simulated for 5 ns. The mean and confident interval of the protein adsorption free energy was computed from three independent sets of US simulations. All simulations with the lysozyme protein were performed at the experimental condition with pH 7.4. In this condition, the arginine and lysine residues were protonated, but the glutamate and aspartate were deprotonated,53 resulting in a net charge of +8 e. Then, the protein above the SAM surface was fully solvated with water, and eight Cl− counterions were added to neutralize the system. It was then energy minimized, followed by a 2 ns position restraint and a MD simulation until 300 ns. The lysozyme protein was also simulated in bulk water as a control. It was conducted in the isothermal−isobaric ensemble (NPT) at 298 K, and an atmospheric pressure of 1 bar was applied using the Parrinello−Rahman barostat. To study the ligand binding activity of lysozyme in the surfaceadsorbed state, we modeled the ligand-bound protein using the final snapshots from the previous lysozyme MD simulations. This ligand a short sugar chain tri-N-acetyl-D-glucosamine (NAG3)was modeled at the active site of lysozyme by superimposition based on the ligand position given in the crystal structure (PDB 1LZB, 1.5 Å resolution). The ligand-bound protein on the surface was simulated for 100 ns. The free energy of protein−ligand binding was computed by the MM/PBSA method with the g_mmpbsa tool54 using the solute dielectric constant of 2, solvent dielectric constant of 80, and the SASA model for nonpolar solvation energy estimation. All other parameters for the binding energy calculation were set to the default values. The ligand-bound protein in the bulk solution state was simulated directly from its crystal structure (PDB 1LZB). Reported binding free energy values were averaged from three replica simulations; each calculation was performed using 50 snapshots from the last 5 ns of the trajectory. Simulation Conditions. All systems were simulated using periodic boundary conditions in all three dimensions. The particlebased cutoff method (Verlet) was employed for neighbor searching. Long-range electrostatics were addressed by the particle mesh Ewald method, and a cutoff of 1 nm was used as the short-range nonbonded cutoff parameter. Other simulation parameters were the same as in our previous simulation study.43 FF parameters of the SAM ligand were generated using the general AMBER force field (GAFF)42 and the am1-bcc charge model.55 The peptide and protein were simulated using the standard AMBER ff14sb parameters41 and the water model using TIP3P.44 In the lysozyme study, modified interfacial force field parameters were used. Packmol,56 ProtPOS,53 and in-house developed scripts were used for preparing initial structures and conducting analysis. The GROMACS version 5.0.7 MD engine57 was used to perform classical and US simulations, while PLUMED version 2.2.458 was used to perform metadynamics simulations. All simulations were performed at 298 K with a time step of 2 fs for equilibration and production, in which data were collected at 1 ps intervals. Protein Orientation. We define the protein orientation based on two angles: α is the angle between the long axis of the protein and the surface plane where the long axis of the protein was taken as the vector of Asn46-Cα to Ala9-Cα. ω is the angle between the surface normal and the vector from the protein center to the active site, where the COM of two catalytic residues, Glu35 and Asp52, was used to define the location of the active site. When α is close to 0°, the protein is referred as side-on; when it is close to 90°, it is referred as end-on. When ω is less than 90°, the protein is referred as face-up; otherwise it is face-down.

(5)

where ϵij is the well depth, a measure of how strongly two atoms i and j attract each other. σij is the distance at which the intermolecular potential between the two atoms is zero, referred as the van der Waals radius. AMBER uses the standard Lorentz−Berthelot combination rules

σij =

σi + σj

ϵij =

2

ϵiϵj

(6)

where i and j represent two interacting atoms. To improve the accuracy of the protein−SAM interaction, we tuned the potential well depths (ϵ) of selected problematic interactions. Because the van der Waals radius (σ) of an atom depends on the atomic volume, it remained unchanged. The peptide adsorption free energy is expected to be influenced not only by the peptide−SAM interaction but also by the water−SAM interaction. In our previous simulation study,43 we have shown that contact angle values obtained from water droplet simulations using GAFF on CH3−SAM are in good agreement with experiments. The result suggested that TIP3P-GAFF nonbonded interaction parameters are appropriate. Therefore, in this study, we focused on the interaction parameters of the peptide and SAM. The experimental free energies of adsorption were used as the benchmark data. First, we obtained the adsorption free energy values using the standard parameters of the GAFF + AMBER ff14sb force field. Then, we compared them to the experimental values to identify the interaction to be corrected. The parameters of the guest residue which had large square error (SE) were subjected to further parameterization. Adjusted parameters were applied, and the affected peptide−SAM simulations were rerun. This process was repeated iteratively until the computed peptide adsorption free energies reached close to the experimental values. A modification was performed by scaling the interaction parameter of atom pairs that were involved in a problematic interaction by an optimization factor (Of) ϵij = Of ×

ϵiϵj

(7)

The optimization factor modifies the result from the Lorentz− Berthelot combination rule calculation for specific atom pair interactions and provides the adjusted value for the interaction parameter. Repeated simulations and adsorption free energy calculations were executed to find the set of optimal Of values. We have introduced new atom types for modified atoms to avoid the influence of modification in the interaction of other atom pairs. Lysozyme Protein System. As a case study and a comparison of the effect introduced by the modified FF parameters, we have simulated the widely studied lysozyme protein on the CH3−SAM surface. The initial orientation of the lysozyme protein (PDB 1AKI) on CH3−SAM was predicted by ProtPOS.53 ProtPOS is a selfcontained, lightweight, and easy-to-use software package for predicting the preferred orientation of proteins on a given surface based on the user-selected FF. It searches quickly for the low-energy protein conformation in all six degrees of freedom of the protein using particle swarm optimization. Each successful run returns the lowest energy conformation of the protein on the surface in PDB format. The following PSO parameters were used in each ProtPOS run: n = 200, c1 = c2 = 1.193, w = 0.721, r = 10. In total, ten ProtPOS runs were conducted. Then, the predicted conformations were subjected to clustering analysis based on pairwise root-mean-square deviation (rmsd) values of the residue minimum distance profile. Finally, the cluster-representative protein orientations were taken as the initial protein structures for the MD and free energy simulations were carried out subsequently. The free energy of protein adsorption was predicted using US simulations. To obtain equilibrated structures, we simulated the ProtPOS predicted protein orientations for 300 ns with conventional MD. 9625

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir

Table 1. Experimental and Simulated Adsorption Free Energies (kcal/mol) of the 11 Host−Guest Peptides (TGTG-X-GTGT, Where X is the Guest Residue) on the Methyl Terminated Alkathiol SAM (CH3−(CH2)11−SH)a standard

modified

guest residue

−X−b

exptc

parameterd

SEe

parameterd

SEe

charged

−R− −K− −D− −L− −V− −W− −F− −T− −G− −S− −N−

−4.15±0.55 −3.34±0.39 −3.54±0.60 −3.87±0.69 −4.4±0.31 −3.89±0.34 −4.16±0.16 −2.76±0.28 −3.4±0.39 −2.75±0.23 −4.33±0.62

−1.336±0.9 −1.574±1.2 −1.786±0.7 −3.121±1.2 −3.368±1.5 −3.839±1.4 −3.985±0.8 −2.005±1.4 −2.193±0.5 −2.363±0.8 −1.676±1.1 Rf 0.317

7.919 3.119 3.078 0.561 1.065 0.003 0.031 0.755 1.207 0.387 7.045 MSEg 1.507

−4.537±0.9 −2.300±0.4 −2.948±1.4 −3.121±1.2 −3.368±1.5 −3.839±1.4 −3.985±0.8 −2.005±1.4 −2.193±0.5 −2.363±0.8 −3.658±1.4 Rf 0.831

0.150 1.082 0.350 0.561 1.065 0.003 0.031 0.755 1.207 0.387 0.452 MSEg 0.648

nonpolar

polar

Simulations were performed using the standard GAFF + AMBER ff14sb parameters and the tuned interaction parameters (see Table 2). bOne ° from ref 32. dSimulated letter animo acid code used (R:Arg, K:Lys, D:Asp; L:Leu, V:Val; W:Trp, F:Phe; T:Thr, G:Gly, S:Ser). cExperimental ΔGads ΔG°ads from US simulations in this study; N = 3, the mean ± 95% CI. eSquare error (SE). fPearson’s correlation coefficient (R). gMean square error (MSE). a

■

RESULTS AND DISCUSSION Adsorption Free Energy of Peptide and Force Field Parameter Modification. Using the standard parameters of GAFF + AMBER ff14sb, we performed US simulations for each of the 11 host−guest peptides in 24 different separation distances from the CH 3 −SAM surface. The window trajectories were then subjected to WHAM analysis to obtain the sampled distribution of states and the PMF profile. Finally, ΔGads ° of peptide adsorption were estimated using the probability ratio method,33 and therefore, the computed values can be directly compared to the experimental ΔGads ° values reported by Wei and Latour.32 As shown in Table 1, adsorption free energies of host−guest peptides obtained from standard parameter simulations are poorly correlated to the experiment with a Pearson correlation coefficient (R) of 0.317 and a mean square error (MSE) of 1.507 kcal/mol. Among the different amino acid classes, charged residues (R, K, and D) exhibit the largest deviation from the experiment with an MSE of 4.7 kcal/mol, whereas nonpolar residues (L, V, W, and P) show the best agreement with only an MSE of 0.42 kcal/mol. For polar residues, except Asn, other residues are in good agreement with the experiment. Thus, a large gap of errors is clearly seen between the charged residues (including the polar Asn) and the rest. The evaluation of the result with standard parameters suggested that the adsorption free energy of the host−guest peptides with X = R, K, D, and N needed to be substantially increased. R and N deviated the greatest from the experiment (SE ≈ 7 kcal/mol) compared to K and D (SE ≈ 3 kcal/mol). We observed that side chains of K, R, and N consist of one or more amine groups (NH, NH2, NH3), whereas side chains of D and N contain a carbonyl group (CO). As proteins on the surface mainly interact with the functional group of SAM, we enhanced the interactions of these residues with the surface by modifying the LJ potential of each atomic pair VLJ(rij), where i is from the CH3 functional group of SAM and j from the amine or carbonyl functional groups of the four amino acids (see Figure 2). The parameter change was made consistently in the group of atoms with similar physicochemical properties by introducing the same scaling factor (Of) to the well depth

Figure 2. Functional groups of selected amino acids (Lys, Asp, Asn, and Arg) and the CH3−SAM ligand for which reparameterization was conducted to fit the experimental ΔG°ads values. Atom labelings correspond to the atom names in Table 2.

calculation of the atomic pair. By systematically adjusting this factor and testing the effect of the new parameters with repeated US simulations, we successfully brought the computed adsorption free energy into closer agreement with experiments. The final set of modified parameters is presented in Table 2. In total, we modified LJ parameters of 14 interaction pairs involving 7 different atom types, in which 19 atoms were from protein residues and 4 atoms from the SAM ligand. To prevent affecting the original protein−water and protein−protein interactions, 7 new atom types for protein residues were introduced. The changes in the PMF profiles after parameter modification are shown in Figure 3. The largest change in the PMF profile was the peptide with guest residue R. It has a much deeper and wider free energy well at the SAM interface. This change was expected because the R side chain contains three amine groups. Energy wells of peptides with guest residues K, D, and N were also lowered to some extent. Experimental data are not available for the other nine natural amino acids residues (X = A, M, Q, E, P, I, C, Y, and H). Among them, the functional group atoms of glutamic acid (E) and glutamine (Q) are similar to aspartic acid (D) and 9626

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir

Table 2. Summary of Modified Well-Depth Parameters of LJ Potential (ϵij) between SAMs and Amino Acid Molecules SAM atom name (atom typea) C12 (c3)

H12A, H12B, H12C (hc)

−X−

amino acid, IUPACb name, (atom typec→ new atom typed)

standardc ϵij

modifiedd (Of)e ϵij

−R− −K− −N−

NE, NH1, NH2 (N2 → NF2) NZ (N3 → NF3) ND2 (N → NF)

0.5706

0.8559 (1.5)

−R− −K− −N−

HE, HH11, HH12, HH21, HH22 (H → HF) HZ1, HZ2, HZ3 (H → HF) HD21, HD22 (H → HF)

0.1734

0.2601 (1.5)

−D− −N−

OD1, OD2 (O2 → O2F) OD1 (O → OF)

0.6342

0.9513 (1.5)

−D− −N−

CG (C → CF) CG (C → CF)

0.4058

0.9131 (2.25)

−R− −K− −N−

NE, NH1, NH2 (N2 → NF2) NZ (N3 → NF3) ND2 (N → NF2)

0.2162

0.8646 (4)

−R− −K− −N−

HE, HH11, HH12, HH21, HH22 (H → HF) HZ1, HZ2, HZ3 (H → HF) HD21, HD22 (H → HF)

0.0657

0.2628 (4)

−D− −N−

OD1, OD2 (O2 → O2F) OD1 (O → OF)

0.2402

0.9610 (4)

−D− −N−

CG (C → CF) CG (C → CF)

0.1537

0.9924 (6.5)

GAFF atom type. bIUPAC: International Union of Pure and Applied Chemistry. cAMBER ff14SB atom type and interaction parameter. dNew atom type with modified interaction parameter. eOf: optimization factor; see eq 7.

a

Figure 3. PMF profiles extracted from one set of US simulation using standard (solid line) and modified (dashed line) force field parameters. Left, middle, and right panel are PMFs of charged, nonpolar, and polar residues, respectively. All profiles reached equilibrium beyond 20 Å of SSD, justifying the choice of the system setup.

asparagine (N), respectively; thus, their parameters were modified for consistency (see Supporting Information, Table S3). Tyrosine (Y) is an aromatic amino acid that is classified into the group of phenylalanine (F) and tryptophan (W) by their similarity in structures and physicochemical properties. The adsorption free energy of F and W were reproducible by standard parameters; therefore, it can be assumed that standard parameters are also able to represent proper adsorption behavior of Y. Similarly, modification of alanine (A) and isoleucine (I) is not required; they are classified into the nonpolar aliphatic group of valine (V) and leucine (L). It is difficult to comment on the protein−SAM interaction of sulfur-containing amino acids cysteine (C), methionine (M), and the cyclic amino acid proline (P) without new experimental data. Histidine (H) is a polar amino acid with an amine group, but the amine group is in a ring conformation, which is different from the linear structure of the lysine (K)

and arginine (R) side chains. Hence, we did not make any correction for H due to lack of experimental data. As shown in Figure 4, the Pearson correlation coefficient (R) between the full set of simulated and experimental free energy values increased from 0.317 to 0.831 with the modified parameters. Adsorption free energies obtained from modified parameters were within approximately 1 kcal/mol of experimental values for most of the residues. The MSE was improved from 1.507 to 0.648 kcal/mol. Therefore, the modified parameters showed considerably better performance than the standard parameters in representing protein−SAM interactions. Lysozyme on CH3−SAM: A Case Study. To explore the effect of the modified parameters on protein−SAM interactions, we simulated the adsorbed state of a lysozyme molecule on the CH3−SAM surface. Adsorbed orientations of the protein were predicted by ProtPOS.53 ProtPOS is an 9627

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir

The ProtPOS scores suggested that the probability of the protein being adsorbed in the face-up orientation was slightly more preferred than the face-down orientation. To check whether the predicted orientations were stable, we performed MD simulations of 300 ns for both adsorbed orientations. The initial and final snapshots of the face-up and face-down simulations are displayed in Figure 5, where the two catalytic

Figure 4. Comparison of the experimental ΔG°ads and predicted values from US simulations. The region within 1.0 kcal/mol difference between the observed and computed results is indicated by two dashed lines.

efficient conformational sampling method to predict the lowenergy orientations of protein at SAM surfaces. Taking the predicted structures further in simulation studies reduces the risk of sampling unsuccessful adsorption events (e.g., the protein bounce-off from the surface). The lysozyme−SAM predicted structures were subjected to 300 ns simulations; they were subsequently analyzed and compared to previous computational studies. It should be mentioned that proteins on hydrophobic surfaces can undergo large conformational changes that may induce protein unfolding.40 However, as the energetic barriers between the different adsorbed states are typically large,59 simulation of the entire adsorption process is very challenging. Here, our relatively short MD simulations capture only the initial adsorbed states of single protein on the surface. Our purpose is to compare the orientations and respective bioactivities of the protein upon initial attachment to the surface. Predicted Protein Adsorption Orientations. We performed ProtPOS predictions using the standard and modified parameters. With the standard parameters ProtPOS predicted two clusters of protein orientations; both of them are side-on and face-down orientations of the protein (see Supporting Information, Table S4). Interestingly, with the modified parameters, ProtPOS predicted two clusters of different protein orientations as shown in Table 3. Both orientations were side-on, meaning that the long axis of the protein is parallel to the SAM surface. One predicted protein orientation (cluster 1) is face-up, where the active site is opened to the bulk solution, whereas the other one (cluster 2) is face-down, where the active site is fully blocked by the SAM surface.

Figure 5. Snapshots of protein orientations on the CH3−SAM surface: (A) face-up and (B) face-down. Initial structures were predicted by ProtPOS using the modified parameters. Proteins are shown in cartoon style with secondary structure coloring. Important residues at the active site are shown using the stick model and colored: Glu 35 (red) and Asp 52 (green); Trp 62 (blue) and Trp 63 (blue).

residuesGlu 35 and Asp 52are highlighted to indicate the location of the active site. During the course of 300 ns simulation, the protein was slightly rotated, but it remained attached to the SAM surface as side-on. This slight adjustment is expected as ProtPOS predicts based on the rigid protein and rigid surface (see Figure 6). Adjusted protein orientations in both simulations resulted in an increase of protein−surface contacts. Using a 5 Å distance cutoff, we found that 15 residues are in close contact with the SAM surface when the protein is in face-up orientation, whereas 6 residues are in close contact when it is in the face-down orientation (see Table 3). The two protein orientations predicted using the standard parameters were also simulated for 300 ns. Their initial and final snapshots of the simulations are shown in the Supporting Information, Figure S3 and the SSD over simulation time in Figure S4. Unlike the cases with the modified parameters, lysozyme underwent many orientational changes. In one run, starting from the side-on orientation, it rolled over to adopt an end-on orientation for almost 200 ns. Afterward, it rolled back

Table 3. Lysozyme initial Adsorption Orientations Predicted by ProtPOS53 Using Modified Parameters and Surface Contact Residues Obtained from Equilibrated Structurea ProtPOS score (kcal/mol) cluster 1 (57.2%) cluster 2 (42.8%)

−137.2 −95.0

predicted initial adsorption orientationb α = −30.07°, ω = 30.78° side-on, face-up α = 22.74°, ω = 118.56° side-on, face-down

adsorptiond free energy (kcal/mol)

major contact residuec Arg14, Asp66, Gly67, Arg68, Arg73, Asn74, Asn77, Ile78, Pro79, Ser81, Ala82, Ser85, Ser86, Asp87, Thr89 Asn37, Asn44, Arg45, Thr47, Asn113, Asn114

−28.8±0.07 −25.1±1.26

a Surface contact residues are protein residues within 5 Å of the surface. bSee the Methods section for definitions. cResidue that is within 5 Å distance from the surface. Average distance of each residue was calculated from the last 50 ns of the 300 ns long simulation. Contact residues that were also observed in Wei et al.23 are indicated by underline. dThe mean (±95% CI) from three independent sets of US simulations.

9628

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir

Figure 6. SSD profiles of the two protein orientations on the CH3−SAM surface using the modified parameters.

Figure 7. Comparison of the lysozyme adsorption behavior simulated using the standard and refined force fields.

to a side-on orientation, similar to the one predicted by ProtPOS (see Supporting Information, Figure S5). In the second run, the protein also went to adopt an end-on orientation. Then, it desorbed from the surface at about 40 ns, rotated, and readsorbed onto the surface with another endon orientation. It further rotated and assumed a different sideon orientation (see Supporting Information, Figure S6). Figure 7 compares the SSD and orientation of the protein (in terms of the angles α and ω) simulated by the standard and refined force fields. It is clear that the protein simulated using the standard parameters exhibits high level of instability and orientational fluctuation on the surface. In contrast, using the modified parameters, the protein shows progressive attachment behavior to the hydrophobic surface in relatively stable orientations. The drastically different adsorption behavior of the protein with the modified parameters is due to the enhanced protein−surface interactions of residues, including Arg, Asp, and Asn, that involved in surface contacts. Although it is unclear from experiments whether or not lysozyme experiences such high instability on hydrophobic surfaces, it is unambiguous that the preferred and stable initial orientations of the protein on CH3−SAM (under the most diluted condition) are side-on orientations.60 With the modified parameters, such orientations can be obtained in shorter simulation time. However, more experimental data, such as free energy of protein adsorption, are needed to confirm about the observed adsorption behavior of the protein on the CH3− SAM surface. Energetics of Adsorption Orientations. To further quantify the energetic differences between the two orientations predicted by the modified parameters, we performed US simulations to obtain the PMF profiles of the adsorption process. The result showed that the adsorption free energies of the two protein orientations were similar: −28.8±0.07 kcal/mol for face-up and −25.1±1.26 kcal/mol for face-down. However, this difference of about 4 kcal/mol amounts to an approximately 500-fold increase in the probability of the protein being adsorbed in the face-up orientation than in the face-down orientation. In addition, because the transition from one to the other state involves at least a rotation of the protein

by about 60° with breaking and reforming a number of residue-surface contacts, the energetic barrier separating these two initial adsorbed states is likely to be high.59 Therefore, the protein will tend to remain close to its initial adsorbed state. Unfortunately, which of the two initial orientations are more preferred in reality is not an easy question to answer. Previously, there were reports of different opinions regarding which surface of the lysozyme molecule is adsorbed on the hydrophobic SAM surface. By an examination with a chromatographic retention technique, Fausnaugh and Regnier reported that only amino acid substitutions that occurred on the protein surface opposite the active cleft affected lysozyme retention on a hydrophobic interaction column.61 They found that residues 41−102 and 75−89 are located at contact areas between the protein and the stationary phase surface. Interestingly, this second range of residues coincides with our observed contact residues in the face-up orientation. In contrast, using side-chain-specific chemical modification and bioactivity assays, Fears et al.62 found that four out of six Trp residues were labeled after adsorption, suggesting that they were still solvent accessible; however, the two Trp residues (Trp 28 and Trp 108) were identified as unlabeled in the experiment, suggesting that they were likely to be in contact with the SAM surface. Our predicted face-down protein orientation also indicated that Trp 108 is in close proximity to the SAM surface with the active site being sterically blocked by the surface. It is worth noting that although both orientations seem possible according to different experiments, previous computational studies using conventional atomistic23,63,64 or coarsegrained65 MD simulations on the same or similar hydrophobic surfaces reported more often the face-down orientation of the protein as its initial adsorption orientations. This may be due to the use of general biomolecular force fields (e.g., CVFF,63 OPLS-AA,23,64 and MARTINI65) which has underestimated charged nonpolar interactions at the SAM interface. Effect of Adsorption Orientation on Protein Bioactivity. Fears at el. revealed that lysozyme reduced its specific hydrolytic activity against bacterial (Staphylococcus aureus) peptidoglycan at the adsorbed state on the CH3−SAM surface; 9629

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir

Figure 8. Comparison of lysozyme (apo) in the face-up adsorbed state and in the bulk solution state in 300 ns simulations. (A) Backbone rmsd with respect to the final solution state structure and (B) superposition of the final structures at 300 ns. (C) Average SASA of active cleft residues over the simulation time. (D) Average SASA of each active cleft residue over the last 50 ns of the simulations.

Figure 9. Snapshots of the protein−NAG3 complex at 10, 50, and 100 ns in simulations of the bulk solution state (upper panel) and adsorbed state (lower panel). NAG3 is represented by the stick model with yellow color; the protein is in the cartoon model. The residues involved in hydrogen bonding with NAG3 in the crystal structure are shown in the sphere model: Asn59 (green), Trp62 (pink), Trp63 (gray), Asp101 (red), Asn103 (tan), and Ala107 (blue).

the amount of reduction was ≈85% to that in the bulk solution state.62 In simulations, we noticed that two catalytic residues Glu35 and Asp52 were faced toward the surface in the case of the face-down orientation; therefore, adsorbed lysozyme in this orientation was completely inactive due to inaccessibility of the catalytic site. On the contrary, in face-up orientation, catalytic residues situated opposite to the surface and were exposed to solvent. To see if the face-up protein is still functionally active, we compared the structural properties of lysozyme in both ligandfree and ligand-bound forms in the bulk solution state and faceup adsorbed state. Lysozyme binds N-acetyl-β-glucosamine (NAG)n oligomers; it hydrolyzes the β-(1,4)-glycosidic linkage and produces tri-N-acetyl-D-glucosamine (NAG3). We modeled NAG3 into the active site of the equilibrated lysozyme structure in the face-up adsorbed state taken the crystal ligandbinding position (PDB ID: 1LZB) as a reference. The ligand-

bound complex on the SAM surface was simulated for 100 ns. For comparison, the ligand-bound complex in the bulk solution was also simulated for 100 ns. Figure 8a shows that in the ligand-free form, the protein in the adsorbed state had a slight increase in rmsd compared to that in the solution state (0.22 and 0.14 nm, respectively). The deviation occurred at the loop region, where residues Asp66− Asn74 landed on the SAM surface during adsorption. A small deviation of the backbone was observed near residues Asn103−Ala107. Apart from the backbone deviation, we noticed that the active site was more exposed in the adsorbed state compared to its solution state. The active cleft residues, including Glu35, Asp52, Leu56, Gln57, Ile58, Trp62, Trp63, Ile98, Asp101, Ala107, and Trp108 (all within 0.5 nm from NAG3 in the crystal structure), were found to have increased solvent accessible surface areas (SASAs). This gave a total increase of 0.2 nm2 in average SASA of all active cleft residues 9630

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir in the face-up adsorbed state with respect to its bulk solution state. To investigate whether the exposed active site has an effect on the capability of lysozyme to bind native ligands, we simulated the NAG3-bound proteins for 100 ns and monitored the distance between the ligand and protein active site. The simulations were repeated three times with different starting velocities. Snapshots of the simulations are shown in Figure 9. The adsorbed complex was slightly tilted but it maintained contacts with the surface; the binding site was more closed than its initial conformation, which resembles its conformation in the solution state. The free energy of binding was evaluated using the MM/PBSA method from the last 5 ns of trajectory. In the crystal structure, the hydrolyzed ligand NAG3 was found to interact with lysozyme via stacking interactions to Trp62 and hydrogen bonding to Trp59, Trp62, Trp63, Asp101, Asn103, and Ala107.66,67 From the free energy decomposition analysis (see Figure 10), we observed that the binding

that lysozyme can reduce or lose its biological function, with high probability, upon initial adsorption.

■

CONCLUSIONS In this work, we assessed the popular AMBER force field against the experimental adsorption free energies of eleven host−guest peptide on CH3−SAM. The result of the comparison revealed that a correction of the force field was required to reproduce protein or peptide activities on SAM surfaces. To this end, we systematically adjust the LJ parameters of the amine and carbonyl groups of selected protein side chains (including Lys, Arg, Asp, and Asn), and the CH3 functional group of the SAM ligand. Modified parameters were tested via repeated simulations and free energy calculations. The final parameter set achieved close agreement in the peptide adsorption free energy with a MSE of 0.65 kcal/ mol compared to 1.5 kcal/mol using standard parameters. As a case study and an application of the refined parameter set, we predicted that the initial lysozyme adsorption poses on the CH3−SAM surface by ProtPOS and subsequently simulated them to assess the bioactivity of the protein. Analysis from the ProtPOS prediction and MD simulations revealed that the protein can adsorb onto the surface in both face-down and face-up orientations. In the face-down orientation, as the active site of the protein is sterically blocked by the surface, its bioactivity is likely to be abolished completely. On the contrary, in the face-up orientation, the active site of the protein is opened to the bulk solution. Nevertheless, the results of protein−ligand binding free energy calculation by MM/ PBSA indicate that the adsorbed protein will bind the NAG3 ligand weaker by ≈8 kcal/mol compared to its free form in the bulk solution. Therefore, the protein activity after initial adsorption can be largely reduced. Force field is undoubtedly one of the most crucial factors to accurately simulate molecular processes. As our first attempt to improve AMBER and GAFF for protein−SAM simulations, we decided to focus on the parameterization of the CH3−SAM molecule. This molecule with the methyl functional group is an important component of many biosensor devices and drug delivery systems. Alkanethiolates are also the basis of other popular SAM molecules and the silane-based molecules widely used to coat silicon wafers, quartz slides, and glass. Thus, the refined force field developed in this work can be applied on simulating SAM molecules with similar moieties or further extended to include new moieties. With improved representation of molecular interactions, new insights of the protein adsorption mechanism and pathway can be obtained, which may contribute to the rational design of nanosurfaces. In our lab, validation of force field for other SAM molecules such as the hydrophilic (OH−) and charged (NH2− and COOH− SAM) moieties are planned; simulation studies of important model proteins for biosensors68 such as protein A and bovine serum albumin are currently underway.

Figure 10. Free energy contribution of each active cleft residue of lysozyme to the binding of NAG3. Starred residues are those involved in hydrogen bonding with the ligand in the crystal structure. Free energy values were averaged from three replica simulations, using 50 snapshots each taken from the 95−100 ns trajectory.

interactions of the adsorbed lysozyme with these residues were slightly weaker except Asp101. The total binding free energy obtained for the adsorbed state was −27.38±2.75 kcal/ mol, which indicates an overall weaker binding than that in the solution state (−35.24±1.62 kcal/mol) and thus a reduced bioactivity of the adsorbed protein, though in face-up orientation, is still expected. To sum up, our results of ligand-free and ligand-bound simulations suggest that regardless of the adsorption orientation, the protein will reduce or lose its ability to function. In the face-down orientation, the active cleft of lysozyme is blocked and not accessible by the ligand, while in the face-up orientation, the accessible but enlarged active cleft has led to a reduction of the protein’s capability to bind the ligand compared to that in the native state. Hence, an overall reduction, but not completely abolished, in the function of the lysozyme protein upon initial adsorption onto the hydrophobic SAM can be expected. Though our single molecule study cannot give a full picture to the activity reduction of lysozyme on the hydrophobic SAM for direct comparison with the bioactivity measurement experiment reported by Fears et al.,62 as that may involve other complicated factors such as blockage of the active site by protein−protein interaction and protein denaturation over long incubation time, our results suggest

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.langmuir.9b01367. Additional data for the simulation parameters and peptide adsorption free energy calculation; additional modified force field parameters; ProtPOS predictions 9631

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir

■

using the standard parameters; additional figures for testing of simulation parameters; and simulation snapshots and analysis of the standard force field simulations (PDF)

(13) Banerjee, I.; Pangule, R. C.; Kane, R. S. Antifouling Coatings: Recent Developments in the Design of Surfaces that Prevent Fouling by Proteins, Bacteria, and Marine Organisms. Adv. Mater. 2010, 23, 690−718. (14) Schreiber, F. Structure and Growth of Self-Assembling Monolayers. Prog. Surf. Sci. 2000, 65, 151−257. (15) Ozboyaci, M.; Kokh, D. B.; Corni, S.; Wade, R. C. Modeling and Simulation of Protein−Surface Interactions: Achievements and Challenges. Q. Rev. Biophys. 2016, 49, No. E4. (16) Bhatia, R.; Garrison, B. J. Phase Transitions in a MethylTerminated Monolayer Self-Assembled on Au{111}. Langmuir 1997, 13, 765−769. (17) Rai, B.; P, S.; Malhotra, C. P.; Pradip; Ayappa, K. G. Molecular Dynamic Simulations of Self-Assembled Alkylthiolate Monolayers on an Au(111) Surface. Langmuir 2004, 20, 3138−3144. (18) Kim, H.; Saha, J. K.; Zhang, Z.; Jang, J.; Matin, M. A.; Jang, J. Molecular Dynamics Study on the Self-Assembled Monolayer Grown from a Droplet of Alkanethiol. J. Phys. Chem. C 2014, 118, 11149− 11157. (19) Ahn, Y.; Saha, J. K.; Schatz, G. C.; Jang, J. Molecular Dynamics Study of the Formation of a Self-Assembled Monolayer on Gold. J. Phys. Chem. C 2011, 115, 10668−10674. (20) Yan, C.; Yuan, R.; Pfalzgraff, W. C.; Nishida, J.; Wang, L.; Markland, T. E.; Fayer, M. D. Unraveling the Dynamics and Structure of Functionalized Self-Assembled Monolayers on Gold using 2D IR Spectroscopy and MD Simulations. Proc. Natl. Acad. Sci. U.S.A. 2016, 113, 4929−4934. (21) Vemparala, S.; Karki, B. B.; Kalia, R. K.; Nakano, A.; Vashishta, P. Large-Scale Molecular Dynamics Simulations of Alkanethiol SelfAssembled Monolayers. J. Chem. Phys. 2004, 121, 4323−4330. (22) He, Y.; Hower, J.; Chen, S.; Bernards, M. T.; Chang, Y.; Jiang, S. Molecular Simulation Studies of Protein Interactions with Zwitterionic Phosphorylcholine Self-Assembled Monolayers in the Presence of Water. Langmuir 2008, 24, 10358−10364. (23) Wei, T.; Carignano, M. A.; Szleifer, I. Lysozyme Adsorption on Polyethylene Surfaces: Why are Long Simulations Needed? Langmuir 2011, 27, 12074−12081. (24) Agashe, M.; Raut, V.; Stuart, S. J.; Latour, R. A. Molecular Simulation To Characterize the Adsorption Behavior of a Fibrinogen γ-Chain Fragment. Langmuir 2005, 21, 1103−1117. (25) Liamas, E.; Kubiak-Ossowska, K.; Black, R.; Thomas, O.; Zhang, Z.; Mulheran, P. Adsorption of Fibronectin Fragment on Surfaces Using Fully Atomistic Molecular Dynamics Simulations. Int. J. Mol. Sci. 2018, 19, 3321. (26) Zhou, J.; Zheng, J.; Jiang, S. Molecular Simulation Studies of the Orientation and Conformation of Cytochrome C Adsorbed on Self-Assembled Monolayers. J. Phys. Chem. B 2004, 108, 17418− 17424. (27) He, C.; Zhang, H.; Lin, C.; Wang, L.; Yuan, S. A Molecular Dynamics Study on the Adsorption of a Mussel Protein on two Different Films: Polymer Film and a SAM. Chem. Phys. Lett. 2017, 676, 144−149. (28) Xie, Y.; Liu, M.; Zhou, J. Molecular Dynamics Simulations of Peptide Adsorption on Self-Assembled Monolayers. Appl. Surf. Sci. 2012, 258, 8153−8159. (29) Liu, J.; Yu, G.; Zhou, J. Ribonuclease A Adsorption onto Charged Self-Assembled Monolayers: A Multiscale Simulation Study. Chem. Eng. Sci. 2015, 121, 331−339. (30) Mücksch, C.; Urbassek, H. M. Forced Desorption of Bovine Serum Albumin and Lysozyme from Graphite: Insights from Molecular Dynamics Simulation. J. Phys. Chem. B 2016, 120, 7889− 7895. (31) Srinivasan, K.; Banerjee, S.; Parimal, S.; Sejergaard, L.; Berkovich, R.; Barquera, B.; Garde, S.; Cramer, S. M. Single Molecule Force Spectroscopy and Molecular Dynamics Simulations as a Combined Platform for Probing Protein Face-Specific Binding. Langmuir 2017, 33, 10851−10860.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +853 8822-4452. Fax: +853 8822-2426, +853 3974285. ORCID

Shirley W. I. Siu: 0000-0002-3695-7758 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Kristyna Pluhackova for sharing the GROMACS files for AMBER14 force field. This project is supported by the Science and Technology Development Fund from Macao S.A.R. (grant no. 066/2016/A) and the University of Macau (grant nos. MYRG2014-00104-FST and MYRG2015-00212-FST). Simulations were performed using the High Performance Computing Cluster (HPCC) facilities provided by the Faculty of Science and Technology (FST) and the Information and Communication Technology Office (ICTO) of the University of Macau.

■

REFERENCES

(1) Flink, S.; Veggel, F. C. J. M. v.; Reinhoudt, D. N. Sensor Functionalities in Self-Assembled Monolayers. Adv. Mater. 2000, 12, 1315−1328. (2) Chaki, N. K.; Vijayamohanan, K. Self-Assembled Monolayers as a Tunable Platform for Biosensor Applications. Biosens. Bioelectron. 2002, 17, 1−12. (3) Samanta, D.; Sarkar, A. Immobilization of Bio-Macromolecules on Self-Assembled Monolayers: Methods and Sensor Applications. Chem. Soc. Rev. 2011, 40, 2567−2592. (4) Bhakta, S. A.; Evans, E.; Benavidez, T. E.; Garcia, C. D. Protein Adsorption onto Nanomaterials for the Development of Biosensors and Analytical Devices: A Review. Anal. Chim. Acta 2015, 872, 7−25. (5) Mani, G.; Johnson, D. M.; Marton, D.; Feldman, M. D.; Patel, D.; Ayon, A. A.; Agrawal, C. M. Drug Delivery from Gold and Titanium Surfaces using Self-Assembled Monolayers. Biomaterials 2008, 29, 4561−4573. (6) Pinholt, C.; Hartvig, R. A.; Medlicott, N. J.; Jorgensen, L. The Importance of Interfaces in Protein Drug Delivery - Why is Protein Adsorption of Interest in Pharmaceutical Formulations? Expert Opin. Drug Deliv. 2011, 8, 949−964. (7) Roach, P.; Farrar, D.; Perry, C. C. Interpretation of Protein Adsorption: Surface-Induced Conformational Changes. J. Am. Chem. Soc. 2005, 127, 8168−8173. (8) Latour, R. Biomaterials: Protein-Surface Interactions. Encyclopedia of Biomaterials and Biomedical Engineering; CRC Press, 2005; Vol. 1, pp 270−278. (9) Senaratne, W.; Andruzzi, L.; Ober, C. K. Self-Assembled Monolayers and Polymer Brushes in Biotechnology: Current Applications and Future Perspectives. Biomacromolecules 2005, 6, 2427−2448. (10) Anderson, J. M. Biological Responses to Materials. Annu. Rev. Mater. Res. 2001, 31, 81−110. (11) Wang, K.; Zhou, C.; Hong, Y.; Zhang, X. A Review of Protein Adsorption on Bioceramics. Interface focus 2012, 2, 259−277. (12) Li, D.; Chen, H. Regulation of Protein/Surface Interactions by Surface Chemical Modification and Topographic Design. ACS Symposium Series; American Chemical Society, 2012; Vol. 1120, pp 301−319. 9632

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633

Article

Langmuir

(53) Ngai, J. C. F.; Mak, P.-I.; Siu, S. W. I. ProtPOS: A Python Package for the Prediction of Protein Preferred Orientation on a Surface. Bioinformatics 2016, 32, 2537−2538. (54) Kumari, R.; Kumar, R.; Lynn, A. g_mmpbsa-A GROMACS Tool for High-Throughput MM-PBSA Calculations. J. Chem. Inf. Model. 2014, 54, 1951−1962. (55) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC model: II. Parameterization and Validation. J. Comput. Chem. 2002, 23, 1623−1641. (56) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157−2164. (57) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (58) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (59) Latour, R. A. Perspectives on the simulation of protein-surface interactions using empirical force field methods. Colloids Surf., B 2014, 124, 25−37. (60) Wertz, C. F.; Santore, M. M. Adsorption and Reorientation Kinetics of Lysozyme on Hydrophobic Surfaces. Langmuir 2002, 18, 1190−1199. (61) Fausnaugh, J. L.; Regnier, F. E. Solute and Mobile Phase Contributions to Retention in Hydrophobic Interaction Chromatography of Proteins. J. Chromatogr. A 1986, 359, 131−146. (62) Fears, K. P.; Sivaraman, B.; Powell, G. L.; Wu, Y.; Latour, R. A. Probing the Conformation and Orientation of Adsorbed Enzymes using Side-Chain Modification. Langmuir 2009, 25, 9319−9327. (63) Raffaini, G.; Ganazzoli, F. Protein Adsorption on a Hydrophobic Surface: a Molecular Dynamics Study of Lysozyme on Graphite. Langmuir 2010, 26, 5679−5689. (64) Wei, T.; Carignano, M. A.; Szleifer, I. Molecular Dynamics Simulation of Lysozyme Adsorption/Desorption on Hydrophobic Surfaces. J. Phys. Chem. B 2012, 116, 10189−10194. (65) Yu, G.; Liu, J.; Zhou, J. Mesoscopic Coarse-Grained Simulations of Lysozyme Adsorption. J. Phys. Chem. B 2014, 118, 4451−4460. (66) Maenaka, K.; Matsushima, M.; Song, H.; Sunada, F.; Watanabe, K.; Kumagai, I. Dissection of protein-carbohydrate interactions in mutant hen egg-white lysozyme complexes and their hydrolytic activity. J. Mol. Biol. 1995, 247, 281−293. (67) Zhong, Y.; Patel, S. Binding structures of tri-N-acetyl-βglucosamine in hen egg white lysozyme using molecular dynamics with a polarizable force field. J. Comput. Chem. 2013, 34, 163−174. (68) Briand, E.; Salmain, M.; Herry, J.-M.; Perrot, H.; Compère, C.; Pradier, C.-M. Building of an Immunosensor: How can the Composition and Structure of the Thiol Attachment Layer Affect the Immunosensor Efficiency? Biosens. Bioelectron. 2006, 22, 440− 448.

(32) Wei, Y.; Latour, R. A. Benchmark Experimental Data Set and Assessment of Adsorption Free Energy for Peptide−Surface Interactions. Langmuir 2009, 25, 5637−5646. (33) Vellore, N. A.; Yancey, J. A.; Collier, G.; Latour, R. A.; Stuart, S. J. Assessment of the Transferability of a Protein Force Field for the Simulation of Peptide-Surface Interactions. Langmuir 2010, 26, 7396−7404. (34) Collier, G.; Vellore, N. A.; Yancey, J. A.; Stuart, S. J.; Latour, R. A. Comparison between Empirical Protein Force Fields for the Simulation of the Adsorption Behavior of Structured LK Peptides on Functionalized Surfaces. Biointerphases 2012, 7, 24. (35) Deighan, M.; Pfaendtner, J. Exhaustively Sampling Peptide Adsorption with Metadynamics. Langmuir 2013, 29, 7999−8009. (36) Biswas, P. K.; Vellore, N. A.; Yancey, J. A.; Kucukkal, T. G.; Collier, G.; Brooks, B. R.; Stuart, S. J.; Latour, R. A. Simulation of Multiphase Systems Utilizing Independent Force Fields to Control Intraphase and Interphase Behavior. J. Comput. Chem. 2012, 33, 1458−1466. (37) Abramyan, T. M.; Snyder, J. A.; Yancey, J. A.; Thyparambil, A. A.; Wei, Y.; Stuart, S. J.; Latour, R. A. Parameterization of an Interfacial Force Field for Accurate Representation of Peptide Adsorption Free Energy on High-Density Polyethylene. Biointerphases 2015, 10, 021002. (38) Abramyan, T. M.; Hyde-Volpe, D. L.; Stuart, S. J.; Latour, R. A. Application of Advanced Sampling and Analysis Methods to Predict the Structure of Adsorbed Protein on a Material Surface. Biointerphases 2017, 12, 02D409. (39) Snyder, J. A.; Abramyan, T.; Yancey, J. A.; Thyparambil, A. A.; Wei, Y.; Stuart, S. J.; Latour, R. A. Development of a Tuned Interfacial Force Field Parameter Set for the Simulation of Protein Adsorption to Silica Glass. Biointerphases 2012, 7, 56. (40) Wei, S.; Knotts, T. A. A Coarse Grain Model for ProteinSurface Interactions. J. Chem. Phys. 2013, 139, 095102. (41) Maier, J. A.; Martinez, C.; Kasavajhala, K.; Wickstrom, L.; Hauser, K. E.; Simmerling, C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11, 3696−3713. (42) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (43) Bhadra, P.; Siu, S. W. I. Comparison of Biomolecular Force Fields for Alkanethiol Self-Assembled Monolayer Simulations. J. Phys. Chem. C 2017, 121, 26340−26349. (44) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (45) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33−38. (46) Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: a Smoothly Converging and Tunable Free-Energy Method. Phys. Rev. Lett. 2008, 100, 020603. (47) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 826−843. (48) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 1977, 23, 187−199. (49) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (50) Jämbeck, J. P. M.; Lyubartsev, A. P. Update to the General Amber Force Field for Small Solutes with an Emphasis on Free Energies of Hydration. J. Phys. Chem. B 2014, 118, 3793−3804. (51) Bastien, L. A. J.; Price, P. N.; Brown, N. J. Intermolecular Potential Parameters and Combining Rules Determined from Viscosity Data. Int. J. Chem. Kinet. 2010, 42, 713−723. (52) Nikitin, A.; Milchevskiy, Y.; Lyubartsev, A. AMBER-II: New Combining Rules and Force Field for Perfluoroalkanes. J. Phys. Chem. B 2015, 119, 14563−14573. 9633

DOI: 10.1021/acs.langmuir.9b01367 Langmuir 2019, 35, 9622−9633