5616
Langmuir 2005, 21, 5616-5626
Prediction of the Orientations of Adsorbed Protein Using an Empirical Energy Function with Implicit Solvation Yu Sun,† William J. Welsh,‡ and Robert A. Latour*,† Department of Bioengineering and the Center for Advanced Engineering Fibers & Films (CAEFF), Clemson University, Clemson, South Carolina 29634, and Department of Pharmacology, University of Medicine and Dentistry of New Jersey (UMDNJ), Robert Wood Johnson Medical School and the Informatics Institute of UMDNJ, Piscataway, New Jersey 08854 Received December 14, 2004. In Final Form: March 31, 2005 When simulating protein adsorption behavior, decisions must first be made regarding how the protein should be oriented on the surface. To address this problem, we have developed a molecular simulation program that combines an empirical adsorption free energy function with an efficient configurational search method to calculate orientation-dependent adsorption free energies between proteins and functionalized surfaces. The configuration space is searched systematically using a quaternion rotation technique, and the adsorption free energy is evaluated using an empirical energy function with an efficient grid-based calculational method. In this paper, the developed method is applied to analyze the preferred orientations of a model protein, lysozyme, on various functionalized alkanethiol self-assembled monolayer (SAM) surfaces by the generation of contour graphs that relate adsorption free energy to adsorbed orientation, and the results are compared with experimental observations. As anticipated, the adsorbed orientation of lysozyme is predicted to be dependent on the discrete organization of the functional groups presented by the surface. Lysozyme, which is a positively charged protein, is predicted to adsorb on its ‘side’ on both hydrophobic and negatively charged surfaces. On surfaces with discrete positively charged sites, attractive interaction energies can also be obtained due to the presence of discrete local negative charges present on the lysozyme surface. In this case, ‘end-on’ orientations are preferred. Additionally, SAM surface models with mixed functionality suggest that the interactions between lysozyme and surfaces could be greatly enhanced if individual surface functional groups are able to access the catalytic cleft region of lysozyme, similar to ligand-receptor interactions. The contour graphs generated by this method can be used to identify low-energy orientations that can then be used as starting points for further simulations to investigate conformational changes induced in protein structure following initial adsorption.
Introduction The control of protein adsorption is a key issue in many applications in biomedical engineering and biotechnology. Over the past decade, biomaterials surface design has largely been driven by an effort to avoid nonspecific protein adsorption combined with the covalent bonding of bioactive short peptides, such as arginine-glycine-aspartic acid (RGD), to surfaces to directly control cellular response.1 However, encouraged by the recognition of the role of preadsorbed proteins in the regulation of cell adhesion,2,3 an alternative and potentially more effective strategy of biomaterials surface design is to prevent nonspecific protein adsorption by specifically controlling the adsorption of proteins in a manner to preserve their natural conformation and then utilizing the inherent bioactivity of the proteins to direct cellular response.4 Other applications, such as biosensors,5 biocatalysis,6,7 and bioseparations,8 will also benefit if better control of adsorbed * Author to whom correspondence should be addressed. E-mail:
[email protected]. † Clemson University. ‡ Robert Wood Johnson Medical School and the Informatics Institute of UMDNJ. (1) Hern, D. L.; Hubbell, J. A. J. Biomed. Mater. Res. 1998, 39, 266. (2) Balcells, M.; Edelman, E. R. J. Cell. Physiol. 2002, 191, 155. (3) Goldstein, A. S.; DiMilla, P. A. J. Biomed. Mater. Res. 2002, 59, 665. (4) Latour, R. A.; Hench, L. L. Biomaterials 2002, 23, 4633. (5) Subrahmanyam, S.; Piletsky, S. A.; Turner, A. P. F. Anal. Chem. 2002, 74, 3942. (6) Le Borgne, S.; Quintero, R. Fuel Process. Technol. 2003, 155. (7) Schuler, C.; Caruso, F. Macromol. Rapid Comm. 2000, 21, 750. (8) Sridhar, P. Chem. Eng. Technol. 1996, 19, 398.
protein orientation and conformation can be achieved through surface design. Many experimental techniques, such as circular dichroism,9 attenuated total reflection Fourier transform IR spectroscopy,10,11 infrared-visible sum frequency generation (SFG) vibrational spectroscopy,12 solid-state NMR,13 total internal reflection fluorescence,14 atomic force microscopy,15,16 and scanning tunneling microscopy,17 have been employed to investigate adsorption-induced structural changes in proteins. While much has been learned from these experimental approaches, they are limited in their ability to provide sufficient molecular-level detail to fully understand protein adsorption behavior. Molecular modeling and simulations provide a complementary approach to these experimental techniques that is specifically suited to theoretically investigate the molecular mechanisms involved in complex molecular interactions, such as those that occur during protein adsorption.18-30 Because of the size of systems involved, empirical force (9) Vermeer, A. W. P.; Bremer, M. G. E. G.; Norde, W. Biochim. Biophys. Acta 1998, 1425, 1. (10) Sharp, J. S.; Forrest, J. A.; Jones, R. A. Biochemistry 2002, 41, 15810. (11) Giacomelli, C. E.; Bremer, M. G. E. G.; Norde, W. J. Colloid Interface Sci. 1999, 220, 13. (12) Kim, J.; Somorjai, G. A. J. Am. Chem. Soc. 2003, 125, 3150. (13) Long, J. R.; Shaw, W. J.; Stayton, P. S.; Drobny, G. P. Biochemistry 2001, 40, 15451. (14) Robeson, J. L.; Tilton, R. D. Langmuir 1996, 12, 6104. (15) Kim, D. T.; Blanch, H. W.; Radke, C. J. Langmuir 2002, 18, 5841. (16) Johnson, C. A.; Yuan, Y.; Lenhoff, A. M. J. Colloid Interface Sci. 2000, 223, 261. (17) Haggerty, L.; Lenhoff, A. M. Biophys. J. 1993, 64, 886.
10.1021/la046932o CCC: $30.25 © 2005 American Chemical Society Published on Web 04/30/2005
Prediction of the Orientations of Adsorbed Protein
field-based methods are required for these types of simulations. While molecular simulation methods have been well-developed for the simulation of the behavior of proteins in solution,31 these methods are often not wellsuited for the simulation of protein adsorption behavior, and new methods must be developed for this application. There are three general areas that are typically addressed in the study of protein adsorption: namely, protein orientation, conformation, and bioactivity. In simulations performed to investigate adsorption-induced effects on a protein’s conformation and bioactivity, the first challenge that must be faced is to select initial orientations of the protein on the surface. Previously conducted studies have typically addressed this issue by simply selecting a series of specific starting orientations for the protein simulations.23,24,28 Others have applied a thermodynamic approach based on a specified energy function coupled with a designated search procedure that is applied to sample a limited set of configurations out of the essentially infinite number of possible orientations. Low-energy configurations are then identified from the results for further evaluation.19,27 A similar concept has been used by Zheng et al., who applied Metropolis Monte Carlo (MC) simulations to identify the lowest-energy orientations of lysozyme at different SAM interfaces using a continuum distancedependent dielectric medium.25 As shown in these studies, energy-based methods for evaluating likely adsorbed protein orientations generally require the use of stochastic methods for the calculation of solvation effects because the explicit treatment of solvent molecules in such a simulation contributes such a large amount of overall energy to the system that the effects of adsorption orientation on energy are very difficult to discern, and they involve calculations that are very computationally expensive. Solvation effects, however, are largely responsible for protein adsorption behavior, and thus must be appropriately considered if protein orientation on a surface is to be accurately predicted. Although often employed,19,23-25 simple dielectric constant-based implicit solvation methods borrowed from algorithms developed for the simulation of proteins in solution, such as the distance-dependent dielectric method,32 can be very misleading when applied to protein adsorption applications because they only dampen electrostatic interactions without accounting for the competitive adsorption characteristics of water for charged and polar functional groups of both the protein and the surface. These types of methods (18) Lu, D. R.; Lee, S. J.; Park, K. J. Biomater. Sci.-Polym. Ed. 1991, 3, 127. (19) Noinville, V.; Vidal-Madjar, C.; Sebille, B. J. Phys. Chem. 1995, 99, 1516. (20) Asthagiri, D.; Lenhoff, A. M. Langmuir 1997, 13, 6761. (21) Ravichandran, S.; Madura, J. D.; Talbot, J. J. Phys. Chem. B 2001, 105, 3610. (22) Zhou, J.; Chen, S. F.; Jiang, S. Y. Langmuir 2003, 19, 3472. (23) Raffaini, G.; Ganazzoli, F. Langmuir 2003, 19, 3403. (24) Raffaini, G.; Ganazzoli, F. Langmuir 2004, 20, 3371. (25) Zheng, J.; Li, L.; Chen, S. F.; Jiang, S. Y. Langmuir 2004, 20, 8931. (26) Smith, J. R.; Knight, D.; Kohn, J.; Rasheed, K.; Weber, N.; Kholodovych, V.; Welsh, W. J. J. Chem. Inf. Comput. Sci. 2004, 44, 1008. (27) Wilson, K.; Stuart, S. J.; Garcia, A.; Latour, R. A. J. Biomed. Mater. Res. 2004, 69A, 667. (28) Agashe, M.; Raut, V.; Stuart, S. J.; Latour, R. A. Langmuir 2005, 21, 1103. (29) Kholodovych, V.; Smith, J. R.; Knight, D.; Abramson, S.; Kohn, J.; Welsh, W. J. Polymer 2004, 45, 7367. (30) Smith, J. R.; Kholodovych, V.; Knight, D.; Welsh, W. J.; Kohn, J. QSAR Comb. Sci. 2005, 24, 99. (31) Richards, F. M.; Eisendberg, D. S.; Kuriyan, J. E. Advances in Protein Chemistry. In Protein Simlations; Dagget, V., Ed.; Elsevier: San Diego, CA, 2003; Vol. 66. (32) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227.
Langmuir, Vol. 21, No. 12, 2005 5617
inherently predict trends in adsorption behavior that are very similar to that which would occur in a vacuum environment, which will be distinctly different than the behavior that would occur in aqueous solution. As addressed by Schaefer and co-workers,33 such methods do not properly account for solvation effects and should not be used for simulations in which there are large changes in solvent-accessible surface area. Thus, more realistic methods of efficiently considering solvation effects are needed to calculate the overall energy of protein adsorption and to guide decisions regarding the selection of initial orientations of a protein on a surface for conducting subsequent simulations to study adsorption-induced effects on protein conformation and bioactivity. The objective of this research was therefore to develop an automated energy-based method to predict adsorbed protein orientation on a surface with the capability to efficiently sample selected degrees of freedom of protein orientation with solvation effects implicitly represented in a manner that directly accounts for the free-energy contributions of functional group wetting and dehydration. The developed method is based on the modification of an existing software program (AutoDock34) that was previously developed as a scoring function-based method for ligand-receptor binding and drug design. Scoring function approaches developed for molecular docking, pioneered by Bo¨hm,35 typically involve two primary components: an efficient strategy to search the designated conformational space of the molecular system and the calculation of system energy for each sampled conformation. For ligandreceptor binding applications, the conformational search usually involves the rotation of selected bonds of the ligand, as well as the positioning of the ligand in the receptor pocket,36-42 and energy calculations typically account for both the nonbonded interactions and conformational energy. These conditions, however, are not directly applicable to the study of adsorbed protein orientation, which requires that similar but different strategies be developed to apply scoring function-based methods for this system. In our revised method, the protein and surface structures are fixed and the conformational search centers on the rotation and translation of the protein over the surface. Accordingly, conformation energy is not considered and the empirical energy function only accounts for contributions from the nonbonded interactions between atoms of the protein and the surface in addition to hydration effects. In this paper, the developed program is demonstrated for the adsorption behavior of lysozyme on functionalized alkanethiol self-assembled monolayer (SAM) surfaces. Results are presented in the form of contour graphs of energy versus adsorbed orientation for selected functionalized SAM surfaces, with distinctly different adsorption behavior predicted for each type of surface. The developed (33) Schaefer, M.; Bartels, C.; Karplus, M. Theor. Chem. Acc. 1999, 101, 194. (34) Morris, G. M.; Goodsell, D. S.; Halliday, R. S.; Huey, R.; Hart, W. E.; Belew, R. K.; Olson, A. J. J. Comput. Chem. 1998, 19, 1639. (35) Bo¨hm, H. J. J. Comput. Aid. Mol. Des. 1994, 8, 243. (36) Keenan, S. M.; Welsh, W. J. J. Mol. Graph Model 2003, 22, 241. (37) Tamura, H.; Yoshikawa, H.; Richard, A. M.; Ross, S. M.; DeLisle, R.; Welsh, W. J. Environ. Health Persp. 2003, 111, 1. (38) Ai, N.; DeLisle, R. K.; Yu, S.; Welsh, W. J. Chem. Res. Toxicol. 2003, 16, 1652. (39) Fang, M. Z.; Wang, Y.; Ai, N.; Hou, Z.; Sun, Y.; Lu, H.; Welsh, W. J.; Yang, C. S. Cancer Res. 2003, 63, 7563. (40) Amin, E. A.; Harris, W. R.; Welsh, W. J. Biopolymers 2004, 73, 205. (41) Nagarajan, K.; Zauhar, R.; Welsh, W. J. J. Chem. Inf. Model. 2005, 45, 49. (42) Keenan, S. M.; Geyer, J. A.; Welsh, W. J.; Prigge, S. T.; Waters, N. C. Comb. Chem. High Throughput Screen. 2005, 8, 27.
5618
Langmuir, Vol. 21, No. 12, 2005
Figure 1. Procedure to calculate the minimum adsorption free energy value for designated orientations of a protein on a surface: (A) Specify a position vector from CG of the protein to a selected point on the protein’s surface with respect to the protein’s global coordinate system (pertaining to PDB coordinates) designated by φi and θj; (B) rotate the protein φ ) 180° - φi with θj held constant to make the position vector point toward the surface; (C) rotate the protein about the position vector by ψ and translate the CG of the protein away from the surface by d while calculating the adsorption free energy at each value of ψ and d. For each particular (φi, θj), the minimum interaction energy is then identified as a function of ψ and d. The minimum energy is then plotted as a function of φ and θ to generate a contour graph that characterizes the energy profile of the predicted protein-surface interactions.
methods can be used to calculate overall initial adsorption free energy for any defined protein-surface system. The energy maps can then also be used to identify low-energy orientations as starting points for further simulations to investigate conformational changes following adsorption. Methods Configuration Space Search. To produce contour maps that characterize adsorption free energy of the protein as a function of its orientation over a given surface, an efficient sampling method is needed to consider up to six degrees of freedom (6 DOF: 3 rotations, 3 translations) that define the protein’s orientation over the surface and to calculate adsorption free energy for each designated orientation. If the adsorbent surface models consist of homogeneous functional groups or well-defined adsorption sites at designated positions, then translational searches are not necessary parallel to the surface plane, thus reducing the number of DOFs for the search process from six to four (i.e., three rotations, one translation normal to the surface plane). Configurations for these four DOFs can then be generated by translation and rotation of the protein molecule by predefined small increments in each degree of freedom. To accomplish this, we developed a C++ program to efficiently sample the configuration space by the quaternion rotation technique43 and energy evaluations were then performed at each configuration by a very efficient grid-based method that is described below. While more efficient sampling algorithms could be applied to identify lowenergy orientations, such as the genetic algorithm, a full systematic search enables more detailed orientational energy maps to be generated, visualized, and ensemble averaged adsorption free energies to be calculated. The procedure used to define the orientation of the protein over each surface is shown in Figure 1. The first two DOFs (φi, and θj) determine a point on the protein surface that will face the adsorbent surface by a position vector, which is defined from the center of gravity (CG) of the protein to the designated point on the protein surface. Once φi and θj are designated, the protein molecule is then rotated to orient the position vector perpendicular to the surface. The remaining two DOFs represent the rotation of the oriented protein about the rotated position vector (same as the global Y axis in Figure 1) by ψ and the translation of the protein along the Y axis by distance d. The interaction energy between the protein and surface is then calculated (as described below) as ψ and d are incrementally varied through their respective ranges, and the low-energy position of the protein with respect to ψ and d is determined. In this way, the minimum (43) Watt, A.; Watt, M. Advanced Animation and Rendering Techniques-Theory and Practice; Addison-Wesley: Wokingham, U.K., 1992.
Sun et al. energy value of adsorption free energy is determined for each designated φi and θj orientation, and a three-dimensional (3-D) energy contour map is then obtained by plotting the minimum adsorption free energy of the protein as a function of φ and θ. The contour map can then be used to identify the locations of the local and global energy minima of the system. In our sampling process, φ and θ were incrementally increased by 5° from 0° to 180° and from 0° to 360°, respectively, ψ was incrementally increased by 10° from 0° to 360°, and d was varied by 0.2 Å increments with its range set to vary the distance between the closest atom on the protein’s external surface and the adsorbent surface from 0 to 12 Å (cutoff distance for van der Waals interactions). During the configuration space search, ensemble adsorption free energy (Eads) of the system was obtained by averaging over all adsorbed states of the protein with Boltzmann weighting:44 36 72 36
Eads )
∑∑ ∑
36 72 36
Pijk × E(φi, θj, ψk, deq)/
i)1 j)1 k)1
∑∑∑P
ijk
(1)
i)1 j)1 k)1
with Pijk ) sin φi × exp(-E(φi, θj, ψk, deq)/RT), where E(φi, θj, ψk, deq) is the energy at the adsorbed state (φi, θj, ψk, deq), Pij is the Boltzmann weighting factor, R is the ideal gas constant, and T is the absolute temperature. For a given orientation of (φ, θ, ψ), the protein is considered to be in an adsorbed state when the protein-surface interaction is attractive and the protein-surface separation distance equals to the equilibrium distance, deq, at which the protein-surface interaction is strongest. Empirical Energy Function. The empirical adsorption free energy function was developed by simplifying the functional form of a molecular docking energy function used by AutoDock.34 The AutoDock energy function consists of five terms: a LennardJones (LJ) 12-6 repulsion-dispersion term, a directional 12-10 hydrogen bonding term, a screened Coulombic electrostatic potential term, a conformational energy term, and a desolvation energy term. For our application, we simplified the functional from of this energy calculation to include only two terms. This was accomplished by first incorporating all terms related to secondary bonding interactions (i.e., the LJ, desolvation and hydrogen bonding terms) within the LJ function to form what we will call the combined van der Waals (CvdW) term. This was done by the adjustment of the LJ parameters (re, e) for pairwise interaction energies between atoms of the protein and the surface to represent the combined effect of all of these secondary bondingtype interactions. Second, because we are only interested in calculating the adsorption free energy for the initial contact between a protein and a surface as a function of the protein’s orientation, our method treats the protein and surface as rigid structures, thus eliminating the need for a conformational energy term. Finally, the Coulombic energy term of AutoDock was maintained in our energy function to represent electrostatic interactions. The final functional form of our adsorption free energy function is shown in eq 2:
∑
E ) Csol
i,j
(
ere12 rij12
-
2ere6 rij6
)
+ H + Celec
∑ i,j
332 × qiqj Drij
(2)
where Csol and Celec are the coefficients of the CvdW and electrostatic terms, respectively, with each term being set to unity for these present studies. Csol is unitless and Celec has units of kcal‚Å/mol‚e2. H is a constant that equals 0.0001 kcal/mol, which is introduced to make the CvdW interaction energy between a hydrophilic atom pair not less than zero. The terms re and e are the designated CvdW equilibrium distance and well-depth for a particular atom pair, respectively (see Table 1), and r, q, and D are the interatomic distance (Å), the atomic charge (e), and the distance-dependent dielectric constant, respectively. As in AutoDock,34 a sigmodial distance-dependent dielectric function (44) McQuarrie, D. A. Statistical Thermodynamics; Harper Row: New York, 1976.
Prediction of the Orientations of Adsorbed Protein
Langmuir, Vol. 21, No. 12, 2005 5619
Table 1. Combined van der Waals Equilibrium Distance and Well Depth Parametersa atom pair
distance, re (Å)
well depth, e (kcal/mol)
C-C C-A C-X A-A A-X X-X
4.0 4.0 6.0 4.0 6.0 6.0
0.2250 0.1750 0.0001 0.1450 0.0001 0.0001
a C ) aliphatic carbon; A ) aromatic carbon; X ) all other hydrophilic atoms, including carbons bonded to heteroatom(s) and hydrogens with partial charges.
was adopted for these calculations according to Mehler and Solmajor:45
D)A+
B 1 + k e-λBr
(3)
where B ) 0 - A, 0 ) 78.4, A ) -8.5525, k ) 7.7839, and λ ) 0.003627 Å-1. It is recognized that the free energy value obtained by eq 2 does not include entropic effects for the protein itself; however, this approach to calculate the free energy state of the system represents a standard method used in the protein folding community to represent the free energy state of a protein in a defined conformation in aqueous solution.46 Following the strategy implemented in AutoDock, in order to provide a very efficient means of calculating atom-atom interaction energies, all atom types are classified as being one of three types; designated as C, A, and X. ‘C’ is used to represent hydrophobic aliphatic carbons, ‘A’ is used to represent hydrophobic aromatic carbons, and ‘X’ is used to represent hydrophilic atoms, which includes nitrogens, oxygens, sulfurs, polar hydrogens, and all metals and carbons bonded to non-hydrogen heteroatoms. The CvdW equilibrium distance and well depth parameters for each resulting atom pair type are listed in Table 1. The values of the parameters in Table 1 and the coefficients for each term were adjusted to reproduce the energy vs surface separation distance relationships established in our previous semiempirical quantum mechanical studies of peptide residuesurface interactions, which included both solvation and hydrophobic effects.4,47-49 Figure 2 shows the pairwise interaction energies as a function of interatomic distance for a representative hydrophobic atom pair C-C, a neutral hydrophilic atom pair N-O, and charged N-O pairs under this parameter set. As shown, the desolvation and adsorption of a hydrophobic atom pair is energetically favorable, the desolvation and adsorption of a hydrophilic atom pair is energetically unfavorable, and opposite and same-charged hydrophilic groups exhibit strong attraction and repulsion energy, respectively, due to electrostatic interactions. Grid-Based Energy Calculation. Because of the large number of configurations that must be sampled to map out the orientational space, a very fast method is needed for the calculation of interaction energy for each orientation. To meet this need, a program called “autogrid” in the AutoDock 3.0 software package34 was modified and used to generate precalculated grid-maps that define all possible pairwise atomic interaction energies for a given functionalized surface. By this method, the adsorbent surface is placed at the bottom of a 3-D grid box and a probe atom of each atom type systematically visits each grid point while the pairwise interaction energies of the probe are summed over all atoms of the adsorbent surface within a nonbonded cutoff radius of 12 Å. Separate grids are calculated for the CvdW interactions for each of the three types of atom in the protein and a separate electrostatic potential grid map is generated for a probe atom with +1 e charge. The energy values (45) Mehler, F. L.; Solmajer, T. Protein Eng. 1991, 4, 903. (46) Wu, X. W.; Brooks, B. R. Biophys. J. 2004, 86, 1946. (47) Latour, R. A.; Rini, C. J. J. Biomed. Mater. Res. 2002, 60, 564. (48) Basalyga, D. M.; Latour, R. A. J. Biomed. Mater. Res. 2003, 64A, 120. (49) Gray, J. J. Curr. Opin. Struct. Biol. 2004, 14, 110.
for each atom type are then stored as a function of grid position. Once the grid-map is generated, the energy for a given protein in a given orientation over the surface can be very quickly calculated by summing the interaction energies of individual protein atoms with the surface, which are obtained by trilinear interpolation on the basis of their positions in the corresponding grid-maps, with the electrostatic interaction energies scaled by each atom’s partial charge. Application to a Protein Adsorption System. The methods described above are applicable to characterize the adsorption free energy of any defined protein as a function of its orientation on a homogeneously functionalized surface or a surface with a designated point of adsorption. In this paper, we demonstrate the application of this method to characterize the adsorption behavior of lysozyme on alkanethiol SAM surfaces at 300 K. Lysozyme Model. Atomic coordinates of the protein lysozyme were taken directly from the Protein Data Bank50 (pdb entry: 7LYZ). Partial charges were assigned according to the GROMACS force field51,52 using a united atom model, which makes hydrophobic atoms (C or A) have a partial charge of 0. At pH 7.4, lysozyme has net charge +8 e. To establish a reference to determine the orientation of lysozyme with respect to the surface, the angle between the long axis of lysozyme and the surface plane, referred to as R, is used to determine ‘end-on’ (i.e., R close to 90°) or ‘sideways’ (i.e., R close to 0°) orientations of lysozyme on the surface. To determine the orientation of lysozyme’s bioactive site with respect to the surface, a vector pointing to the active site is defined as extending from the CG of lysozyme to the midpoint of two atoms that span lysozyme’s bioactive cleft. The angle between this vector and the normal vector of the surface plane, referred to as ω, is used to differentiate between ‘face-up’ and ‘face-down’ orientations of lysozyme’s active site with respect to the SAM surface (i.e., if ω < 90°, the bioactive site is oriented toward the solution phase; if ω > 90°, the bioactive site is oriented toward the SAM surface). Adsorbent Surface. SAM surfaces on a gold (111) plane with both homogeneous and mixed surface functional groups were constructed with dimensions 58 Å × 56 Å × 11 Å (thickness) to represent model adsorbent surfaces for this study using previously published methods.4,47,48 For each defined SAM surface, energy minimization and a 10 ps molecular dynamics simulation were first performed in a water box with periodic boundary conditions using GROMACS software with the GROMACS force field51,52 to enable the surface functional groups to equilibrate with respect to one another in a hydrated environment. All SAM surfaces consisted of 143 chains with the formula CH3-(CH2)nX, where X ) CH3, OH, COOH, or NH2, with n ) 7 for OH-terminating group and n ) 8 for CH3, COOH, and NH2terminating groups. United atom groups were used to represent the -CH2- segments of each alkane chain. When modeling COOH- and NH2-SAM surfaces, decisions had to be made regarding the protonation state of these functional groups in the surface plane. The pKa value for a COOH functional group in dilute solution is typically about 4.0, depending on the specific chemistry involved. However, as COOH groups become concentrated, as is the condition for a COOH-SAM surface, deprotonation becomes inhibited and the pKa value shifts upward. On the basis of the work by Creager and Clarke,53 the pKa value of a COOH-SAM surface is estimated to be about 8.7, which then provides a deprotonation/protonation ratio for the COOHterminated SAM chains of about 1/19. This can be represented discretely as 1 out of 20 carboxyl groups being represented as COO- and 19 groups being represented as COOH, or in an average sense with each COO--terminated SAM chain being assigned a charge state of -0.05 e. A similar shift in protonation state, but in the opposite sense, was assumed for the NH2 functional groups on the NH2-SAM surface. Accordingly, to create homogeneously charged surfaces in this work, the COOH- and NH2-SAM (50) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. Nucleic Acids Res. 2000, 28, 235. (51) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Comm. 1995, 91, 43. (52) Lindahl, E.; Hess, B.; van der Spoel, D. J. Mol. Mod. 2001, 7, 306. (53) Creager, S. E.; Clarke, J. Langmuir 1994, 10, 3675.
5620
Langmuir, Vol. 21, No. 12, 2005
Sun et al.
Figure 2. Pairwise interaction energies as a function of interatomic distance for a representative hydrophobic atom pair C-C, a hydrophilic atom pair N-O, and charged N-O pairs. Table 2. Description of SAM Surface Models surface model
description
OH-SAM CH3-SAM COO-SAM 1COO/OH-SAM 1NH3/OH-SAM 3CH3/OH-SAM 1COO/3CH3/OH-SAM 1NH3/3CH3/OH-SAM
143 single-typed OH-terminated chains; Net charge: 0; Partial charge: O(-0.566 e), H(0.303 e), head carbon(0.263 e). 143 single-typed CH3-terminated chains; Net charge: 0 e. 143 single-typed COO-terminated chains; Each chain assigned net charge -0.05 e; Partial charge: O(-0.0396 e), head carbon(0.03395 e), second to head carbon(0.00475 e). 1 COO-terminated chain, net charge: -1; partial charge: O(-0.792 e), head carbon(0.679 e), second to head carbon(0.095 e), placed at the center of otherwise OH-terminated SAM 1 NH3-terminated chain, net charge: +1; partial charge: H(0.326 e), N(-0.073 e), head carbon(0.095 e) placed at the center of otherwise OH-terminated SAM 3 CH3-terminated chains forming a triangle placed at the center of otherwise OH-terminated SAM 1 COO-terminated and 3 CH3-terminated chains forming a diamond and placed at the center of otherwise OH-terminated SAM 1 NH3-terminated and 3 CH3-terminated chains forming a diamond and placed at the center of otherwise OH-terminated SAM
Table 3. Interaction Energies and Orientations of Lysozyme on SAM Surfaces SAM surfaces
Emina (kcal/mol)
OH-SAM CH3-SAM COO-SAM 3CH3/OH-SAM 1COO/OH-SAM 1NH3/OH-SAM 1COO/3CH3/OH-SAM 1NH3/3CH3/OH-SAM
-1.01 -6.56 -8.41 -2.03 -6.30 -1.19 -7.99 -1.84
orientationb φ, θ, ψ (deg) (120, 220, 240) (110, 160, 260) (90, 290, 260) (130, 280, 30) (135, 280, 300) (90, 205, 340) (125, 285, 340) (105, 215, 60)
Rc (deg)
ωd (deg)
Eadse (kcal/mol)
56.7 19.8 27.9 26.3 23.8 67.0 24.9 68.4
2.7 53.2 74.2 51.8 51.3 30.6 56.6 16.2
-0.21 -4.31 -7.32 -0.23 -4.30 -0.60 -5.92 -0.78
a Global minimum free energy. b Orientation corresponding to global minimum free energy. c Angle between the long axis of lysozyme and the surface. d Angle between the surface normal vector and the vector pointing to the bioactive site. e Boltzmann-weighted average adsorption free energy (from eq 1).
surfaces were represented using COO- and NH3-terminated SAM chains, with each function group presenting a net partial charge of -0.05 e and +0.05 e, respectively. However, when single carboxyl or amine functional groups were represented in the mixed SAM surfaces, these groups were assigned charges of 1 e and + 1 e, respectively. The eight SAM surface models used in this study are described in Table 2, with the carbon directly bonded to the X functional group designated as the “head carbon”. In addition to these eight surfaces, lysozyme interactions were also investigated with single chains of each individual type of functionalized alkane.
Results and Discussion Homogeneous SAMs. The global minimum energy orientations and energies of lysozyme on selected SAM
surfaces are summarized in Table 3. The calculated Boltzmann average adsorption free energy for lysozyme on the OH-SAM is small (-0.21 kcal/mol), indicating that adsorption of lysozyme on a neutral hydrophilic surface is very weak. As shown in the contour graph (Figure 3a), the low-energy orientation of lysozyme on an OH-SAM is in the region of 100° e φ e 130° and 200° e θ e 250°. However, as shown in the 3-D plot for lysozyme interaction with the OH-SAM surface (Figure 3a), the energy barriers between the various local minima are so small (
