Mesoscopic Simulation of Adsorption of Peptides in a Hydrophobic

Jan 28, 2005 - on nine protected peptides covering a wide range of structures. There is a ... on adsorption, especially in hydrophobic systems. Some e...
1 downloads 0 Views 290KB Size
Download PDF
Anal. Chem. 2005, 77, 1243-1252

Mesoscopic Simulation of Adsorption of Peptides in a Hydrophobic Chromatography System Kosta Makrodimitris,†,‡ Erik J. Fernandez,† Thomas B. Woolf,‡ and John P. O’Connell*,†

Department of Chemical Engineering, University of Virginia, Charlottesville, Virginia 22904, and Departments of Physiology, Biophysics and Biophysical Chemistry, and Biomedical Engineering, Johns Hopkins University, Baltimore, Maryland 21205

Mesoscopic simulations using Langevin dipoles on a lattice for the solvent and calculated partial charges for the solute have been used to estimate free energies of adsorption from data on reversed-phase chromatography on nine protected peptides covering a wide range of structures. There is a single parameter, the effective solvent dipole moment, that is fit to data for one peptide and used to predict properties of the other eight peptides. Good agreement of adsorption chemical potentials, including order of chromatographic retention times, is found for calculations that are Boltzmann-averaged over a set of orientations. In addition, the results suggest that there are preferential orientations for each peptide at the model hydrophobic chromatographic surface. Estimation methods for adsorption based on molecular descriptors and hydrophobicity scales are shown to be unreliable for these systems. With refinements and extensions, this simulation method should be applicable to solvents containing salt, such as in hydrophobic interaction chromatography, and to larger solutes including proteins. The manipulation of protein adsorption to hydrophobic surfaces is of growing importance for novel biomaterials1 and biomembranes and for purification of biomolecules at scales from microanalytical to preparative.2-4 The ability to predict or rapidly identify optimum separation conditions for new and difficult chromatographic separations depends on knowing the relationship of molecular structure to the thermodynamics of adsorption and conformation. This can be challenging especially for hydrophobic chromatography used for high-resolution peptide and protein purification, because it is a complex function of many variables, including salt and cosolvent types and concentrations, as well as pH and temperature.5 The behavior of such systems is further * To whom correspondence should be addressed. E-mail: [email protected]. † University of Virginia. ‡ Johns Hopkins University. (1) Ratner, B. D.; Hoffman, F. J.; Schoen, J. E.; Lemon, F. Biomaterials Science. An Introduction to Materials in Medicine; Academic Press: New York, 1996. (2) Queiroz, J. A.; Tomaz, C. T.; Cabral, J. M. J. Biotechnol. 2001, 87, 143159. (3) Norde, W. In Biopolymers at Interfaces, 2nd ed.; Malmsten, M., Ed.; Marcel Dekker: New York, 2003; Vol. 75, pp 21-45. (4) Hearn, M. T. W. In HPLC of Biological Macromolecules: Methods and Applications, 2bd ed.; Gooding, K. M., Regnier, F. E., Eds.; Marcel Dekker Inc.: New York, 2001; pp 99-245. (5) Gagnon, P.; Grund, E. Biopharm-Technol. Business Biopharm. 1996, 9, 54. 10.1021/ac048812r CCC: $30.25 Published on Web 01/28/2005

© 2005 American Chemical Society

complicated by the fact that, in reversed-phase chromatography (RPC), partial or complete unfolding of proteins can occur.6-10 Though hydrophobic interaction chromatography (HIC)11,12 involves surfaces of lower hydrophobicity and, thus, more often retains stable protein conformations,6 this effect can still be observed. Unfortunately, very limited data are available to document protein behavior as a function of many of these variables on either chromatographic13 or even on well-defined surfaces.14 A variety of theoretical and computational approaches have been taken to treat such systems with species ranging from small molecules to peptides to proteins. For peptides and stable proteins that retain their native structure upon adsorption, the solvophobic theory, initially developed for interphase partitioning of small species, has been a starting point for describing the effects of salt type and concentration on adsorption.15,16 Preferential interaction theory,17 which was used initially to study the interactions of solutes and solvent with proteins, has also been extended to describe protein retention in HIC.18,19 Simplified statistical mechanical models based on lattice structures of ligands have been used to treat their responses to environmental conditions.20 However, these theoretical approaches have not generally accounted for the effects of solute or surface molecular structure on adsorption, especially in hydrophobic systems. Some empirical models, based on hydrophobicity scales, have been successful in predicting retention for unstructured peptides for a single set of conditions.21,22 Asenjo and co-workers23,24 have recently used (6) Fausnaugh, J. L.; Kennedy, L. A.; Regnier, F. E. J. Chromatogr. 1984, 317, 141-155. (7) McNay, J. L. M.; Fernandez, E. J. Biotechnol. Bioeng. 2001, 76, 233-240. (8) Oroszlan, P.; Blanco, R.; Lu, X. M.; Yarmush, D.; Karger, B. L. J. Chromatogr. 1990, 500, 481-502. (9) Tibbs, T. L.; Fernandez, E. J. J. Colloid Interface Sci. 2002, 259, 27-35. (10) White, S. H.; Wimley, W. C. Annu. Rev. Biophys. Biomol. Struct. 1999, 28, 319-365. (11) Shaltiel, S.; Er-El, Z. Proc. Natl. Acad. Sci. U.S.A. 1973, 70, 778-781. (12) Hjerten, S. J. Chromatogr. 1973, 87, 325-331. (13) Machold, C.; Deinhofer, K.; Hahn, R.; Jungbauer, A. J. Chromatogr., A 2002, 972, 3-19. (14) Ostuni, E.; Grzybowski, B. A.; Mrksich, M.; Roberts, C. S.; Whitesides, G. M. Langmuir 2003, 19, 1861-1872. (15) Vailaya, A.; Horva´th, C. J. Chromatogr., A 1998, 829, 1-27. (16) Melander, W.; Horvath, C. Arch. Biochem. Biophys. 1977, 183, 200-215. (17) Arakawa, T.; Timasheff, S. N. Biochemistry 1982, 21, 6545-6552. (18) Roettger, B. F.; Myers, J. A.; Ladisch, M. R.; Regnier, F. E. Biotechnol. Prog. 1989, 5, 79-88. (19) Perkins, T. W.; Mak, D. S.; Root, T. W.; Lightfoot, E. N. J. Chromatogr., A 1997, 766, 1-14. (20) Dill, K. A. J. Phys. Chem. 1987, 91, 1980-1988. (21) Meek, J. L. Proc. Natl. Acad. Sci. U.S.A. 1980, 77, 1632-1636. (22) Biswas, K.; DeVido, D.; Dorsey, J. J. Chromatogr., A 2003, 1000, 637-655.

Analytical Chemistry, Vol. 77, No. 5, March 1, 2005 1243

protein structure information to calculate an average surfaceaccessible hydrophobicity, which was correlated with adsorption on hydrophobic surfaces. Quantitative structure-activity relationship methods have been used to identify the molecular properties that seem to play dominant roles in adsorption of known proteins and predict adsorption for unknown or “test” proteins bound to charged surfaces,25 though none seem to have been applied to protein adsorption on hydrophobic surfaces. Molecular simulation is an attractive alternative approach because it may provide both thermodynamic and detailed information about protein-surface interactions. The results should be able to test hypotheses about the most significant molecular interactions involved, as well as to guide development of new adsorbing (or nonadsorbing) surfaces. However, current computational capabilities restrict atomistic simulations to only limited structural elements of the proteins of interest. For the case of RPC, certain effects of alkyl ligand length and degree of substitution on ligand conformation and mobility have been investigated.26 The free energy of adsorption of an idealized small spherical solute to an RPC surface has also been calculated.27,28 Finally, the radial distribution functions of a peptide near an alkylsilica RPC sorbent have been calculated using a four-step docking atomistic procedure.29 However, none of these approaches seem likely to be extended to full-scale protein adsorption on hydrophobic surfaces. As a result, approximate simulation methods based on mesoscale properties or with fields have been used to increase the solute sizes that can be studied. While some molecular details may be sacrificed, certain important aspects can be elucidated. For example, an effective dielectric solvent medium has been able to capture the structural aspects of unfolding of some albumin subdomains on a model graphite surface.30 A continuum electrostatic approximation plus a Hamaker treatment of van der Waals interactions has been applied to the quantitative treatment of protein adsorption to charged, ion-exchange chromatography surfaces.31,32 Simplified treatments of electrostatics, along with explicit van der Waals forces from quantum chemistry calculations, have demonstrated preferred binding orientations through Monte Carlo computations.33,34 However, the focuses of these studies are generally structural, and there seem not to have been studies attempting to predict adsorption strength that can be validated by experimental adsorption or chromatographic data.32 Hydrophobic adsorption has many common issues with other phenomena in biophysical systems, especially those involving interfaces and membranes. Various theoretical and computational (23) Lienqueo, M. E.; Mahn, A.; Asenjo, J. A. J. Chromatogr., A 2002, 978, 7179. (24) Lienqueo, M. E.; Mahn, A.; Vasquez, L.; Asenjo, J. A. J. Chromatogr., A 2003, 1009, 189-196. (25) Mazza, C. B.; Sukumar, N.; Breneman, C. M.; Cramer, S. M. Anal. Chem. 2001, 73, 5457-5461. (26) Yarovsky, I.; Aguilar, M. L.; Hearn, M. T. W. Anal. Chem. 1995, 67, 21452153. (27) Klatte, S. J.; Beck, T. L. J. Phys. Chem. 1996, 100, 5931-5934. (28) Slusher, J. T.; Mountain, R. D. J. Phys. Chem. B 1999, 103, 1354-1362. (29) Yarovsky, I.; Hearn, M. T. W.; Aguilar, M. I. J. Phys. Chem. B 1997, 101, 10962-10970. (30) Raffaini, G.; Ganazzoli, F. Langmuir 2003, 19, 3403-3412. (31) Roth, C. M.; Lenhoff, A. M. Langmuir 1993, 9, 962-972. (32) Roth, C. M.; Lenhoff, A. M. Langmuir 1995, 11, 3500-3509. (33) Hallgren, E.; Kalman, F.; Farnan, D.; Horvath, C.; Stahlberg, J. J. Chromatogr., A 2000, 877, 13-24. (34) Zhou, J.; Chen, S. F.; Jiang, S. Y. Langmuir 2003, 19, 3472-3478.

1244

Analytical Chemistry, Vol. 77, No. 5, March 1, 2005

treatments have been used for such systems, with species ranging from small molecules to peptides to proteins.35-37 Similar progress and limitations to adsorption studies have arisen in these systems as well. In this paper, we describe a first step toward quantitatively estimating protein adsorption to hydrophobic surfaces. Identified as the adsorptive Langevin lattice dipole (ALLD) method, our mesoscale simulation incorporates atomistic description of solutes with a lattice-dipole description of both the solvent and uncharged surface regions.38 The concept is similar to that of Grossfield et al.39 for peptides in cells with membranes. Here, we have computed the thermodynamics of adsorption of protected peptides onto a model RPC surface and directly tested the results with new experimental adsorption data for nine peptides covering a wide range of structure, retention, and hydrophobicity. Aspects of preferred orientations in adsorption as well as comparisons with empirical descriptor and hydrophobicity scales are discussed.

MATERIALS AND METHODS Peptides. Fifty different peptides were evaluated for comparing experiment with simulation. Because the current stage of the simulation method is not intended to treat full charges on the solute or the surface, only N- and C-terminally protected peptides were considered. To be efficient in measurement and calculation expense and effort, nine protected peptides covering a wide range of hydrophobic characteristics and practical retention volumes were chosen, purchased from Sigma-Aldrich (St. Louis, MO), and used without further purification. The properties of these peptides, mostly calculated by the MOE software40 are shown in Table 1, ordered by retention volume (Vr) with indexes from one to nine. Included are the Meek hydrophobicity scale21 (HS; 23-88), radius of gyration (RG; 0.45-0.68 nm), molecular volume (MV; 0.360.73 nm3), accessible hydrophobic surface area (HSA; 6.8-10.6 nm2), dipole moment (DM; 1.9-5.7 D), and number of atoms (#A; 49-103). One indication of the range of structure-property variations is that the least adsorbing peptides 1 and 2, as well as the more strongly adsorbing peptide 7, all have larger hydrophobic surface areas, volumes, and sizes than the others. Chromatography. Reversed-phase chromatography was performed using C18-alkyl bonded silica particles and columns that were a gift of Novo Nordisk (Gentofte, Denmark). The silica particles were 15 µm in diameter with average pore size of 20 nm. Isocratic chromatography was performed with a mobile phase of 70 v/v% H2O, 30 v/v% acetonitrile, and 0.01% trifluoroacetic acid, pH 2.2. The temperature was controlled at 298 ( 0.5 K. Replicate measurements gave a precision in retention volumes of ∼(2.5%. In the chromatography experiments, the equilibrium adsorption constant, K, is determined from peptide and inert tracer retention volumes, Vr and Vo, respectively, and the column mobilephase volume, VM

K ) (Vr - Vo)/VMφ

(1)

(35) Dill, K. A.; Flory, P. J. Proc. Natl. Acad. Sci. U.S.A. 1980, 77, 3115-3119. (36) Wilson, M. A.; Pohorille, A. J. Am. Chem. Soc. 1994, 116, 1490-1501. (37) Zangi, R.; de Vocht, M.; Robillard, G.; Mark, A. Biophys. J. 2002, 83, 112124. (38) Warshel, A.; Russell, S. T. Q. Rev. Biophys. 1984, 17, 283-422.

Table 1. Protected Peptide Formulas and Experimental Results K (eq 1)

∆G° (eq 3)

HS

RG (nm)

MV (nm3)

HSA (nm2)

DM (D)

#A

9.3

7.5

-5.0

23.52

0.63

0.56

9.23

1.87

80

2

11.3

9.6

-5.6

31.57

0.66

0.57

9.68

3.52

81

3 4 5 6 7

11.7 19.6 31.4 35.7 36.3

10.0 18.4 30.7 35.3 35.9

-5.7 -7.2 -8.5 -8.8 -8.9

49.44 27.97 41.53 58.47 59.18

0.44 0.47 0.56 0.51 0.67

0.41 0.36 0.43 0.48 0.73

6.86 6.98 7.69 8.10 10.63

3.37 5.02 3.27 3.40 1.91

60 49 59 68 103

8 9

96.7 151.4

99.5 157.0

-11.4 -12.5

25.42 87.57

0.47 0.49

0.46 0.53

7.34 8.20

5.68 5.55

65 74

protected peptide

index

N-methoxysuccinyl-Ala-Ala-ProVal-p-nitroanilide (AAPV) N-methoxysuccinyl-Ala-Ala-ProMet-p-nitroanilide (AAPM) Z-Pro-Leu-Gly amide (3-PL) Z-Met-Gly-ethyl ester (MG) Z-Val-Tyr-methyl ester (VY) pGlu-Phe-Leu-p-nitroanilide (FL) Boc-b-Ala-Trp-Met-Asp-Pheamide (AWMDF) Z-Gly-Gly-Leu-p-nitroanilide (GGL) N-formyl-Met-Leu-Phe-benzyl ester (MLF)

1

Vr (mL)

where the phase ratio, φ, is the ratio of stationary to mobile-phase volumes, VS and VM, respectively.

φ ) VS/VM ) (V - VM)/VM

(2)

In these equations, V is the total column volume calculated from the dimensions and VM ) Vo - Vexcol, where Vexcol is the volume of the chromatographic system that is external to the column. Sodium nitrate was used to determine Vo and Vexcol; the value of φ was 0.434. The free energy (or chemical potential) difference between the mobile and stationary phases was calculated from the equilibrium constant as

∆G0 ) -RT ln K

(3)

Table 1 shows retention volumes (Vr) measured for each of the peptides under isocratic elution conditions and resulting adsorption equilibrium properties. The retention volume precision yields an uncertainty in free energy difference of 0.05-0.10 kJ mol-1. In the figures that follow, the nine peptides are indexed from 1 (least retained) to 9 (most retained) in order of the retention data in Table 1. While there may be uncertainty about the appropriate phase volume ratio to use in solute partitioning between solvent and surface, since comparisons are made in ∆G h0 any change in the phase volume is an additive constant and is thus canceled out. SIMULATION METHODOLOGY We have implemented a mesoscale simulation method in order to ultimately treat larger solutes such as proteins with some level of structural detail without requiring excessive computational time. Thus, the region from the mobile phase to the surface of the stationary-phase substrate has been represented as a lattice of dipoles. The peptide solutes have been represented atomistically by a combination of partial charges obtained from computational chemistry software. The thermodynamics of adsorption were then determined from the difference between the chemical potential value found from replacing lattice dipoles with the solute in only (39) Grossfield, A.; Sachs, J.; Woolf, T. B. Proteins: Struct. Funct. Genet. 2000, 41, 211-223. (40) In Molecular Operating Environment; Chemical Computing Group Inc.: Montreal, Canada, 2003.

solvent and that found in the interfacial region. Details of the method are given in this section. Solvent Representation. Lattice-dipole models use orientationally variable dipoles placed at specified grid points in order to represent the environment that changes when a solute is inserted into the system. In the present case, a rectangular simulation box was set up longer in the direction normal to the surface (Z) than along the surface (XY) with periodic 2D boundary conditions for interactions in xy planes. The dimensions were LX ) LY ) 3.2 nm and LZ ) 6.4 nm in accordance with other simulations of chromatographic systems27 The mean lattice spacing was chosen as Sd ) 0.4 nm, giving 1024 dipole grid points. Larger box sizes (LX ) LY ) 4.8 nm, LZ ) 6.4 nm) and (LX ) LY ) LZ ) 6.4 nm) gave almost identical results for our peptides but demanded much more CPU time. However, this is still much less than would be needed for a molecular dynamics simulation where ∼2600 water molecules would be required.41 Similar to the method of Grossfield et al.,39 our dipole moments on the lattice varied in the Z direction as

µi )

[ (

) ]

Zi - A ∆µ tanh +1 2 B

(4)

where Zi is the position of dipole i and ∆µ is the difference of dipole moments between the solvent and surface. The parameters A and B set the location and steepness of the dipole profile; here we set A ) 1.6 nm and B ) 0.4 nm, giving the variation shown in Figure 1. Also, we have set Z ) -0.4 nm as the center of mass of the surface silicon atoms of the substrate. As the figure suggests, the function is chosen so that the increase of dipole moment is in the region near the end of the C18 ligands attached to the surface. For perspective, a structure of the VY peptide is also shown; its radius of gyration is 0.5 nm. As with Grossfield et al.,39 a shell of diameter of Dd ) 0.2 nm was put around each dipole. All dipoles for which a portion of the solute overlapped this shell were removed from the lattice to make the cavity. Both our work and that of others suggest that the final results are relatively insensitive to the specifics of the lattice structure. (41) Rossky, P. J.; Karplus, M. J. J. Am. Chem. Soc. 1979, 101, 1913-1937.

Analytical Chemistry, Vol. 77, No. 5, March 1, 2005

1245

The energy of a discrete dipole, with (vector) dipole moment µ, in a field, E, is a function of the angle φ between the dipole and the field:

U (φ) ) - µ‚E ) - µE cos φ

(5)

In a fluid of multiple dipoles, to obtain an appropriate average dipole that accounts for the influence of thermal randomization of dipole orientations, the permanent ideal dipole moment is represented by a Langevin function

〈µ〉 ) µE ˆ [L(yL)]

(6)

L(yL) ) coth(yL) - 1/yL

(7)

where L(yL) is

with the Langevin argument Figure 1. Variation of solvent lattice dipoles. Relative sizes of an alkyl C18 chain and the protected peptide VY compared with lattice spacing.

yL ) µE/kBT

(8)

with E is the field from the other solvent dipoles plus the external field from the solute and E ˆ is the unit vector in the field direction.44,45 This calculation is done for all of the solvent dipoles remaining after the solute is introduced. Since each solvent dipole is influenced by the other solvent dipoles, it is necessary to do an iterative calculation to arrive at a self-consistent field and orientation of the dipoles on the lattice. On lattice point j at iteration n, the magnitude and the orientation of the electric field Enj and the dipole moment µnj is solved selfconsistently via the following equations38

Solute Representation. The CHARMM force field42 was employed for the atomistic representation of the amino acid structure of the peptide solutes. Our solutes also had protected groups such as NO2 that were not amino acids and, thus, not represented in the CHARMM program. For these, we used the partial bond charge increment method of the Merck molecular force field.43 It is based on formal charges, atom types, electronegativities, bond lengths, and bonded neighbors of the force field in order to assign partial charges in groups and atoms that CHARMM, AMBER, or other force fields cannot handle. Every peptide was treated exhaustively with three methods to effect an energy minimization: steepest descent, conjugate gradient, and truncated Newton as included with the MOE utility.40 We kept the configurations of the atoms in the peptides fixed for all calculations at positions and orientations relative to the surface. Other research groups38 have successfully used X-ray and minimized structures with lattice Langevin models in solutions to model proteins and enzymes. In addition, our previous work7 with proteins indicates that solute conformational effects cause broadened and asymmetric chromatography peaks. With our peptides, the peaks were sharp and symmetric. Self-Consistent Lattice-Dipole Model for Electrostatic Interactions. The response of aqueous solvent to insertion of a solute is mostly determined by polarization effects: ionic, atomic, and orientational. The last is of greatest interest here, though atomic polarization may be of significance when the solvent has large concentrations of organic cosolvents, and ionic polarization could arise if the silica substrate is solvent-exposed. Orientational polarization occurs when the field from an inserted solute induces alignment of the permanent dipoles of the solvent that rotate freely and randomly in solute-free solvent. The result is both energetic and entropic.

where rij and rjk are separation distances of charges from dipoles and of dipoles from dipoles, qi is the partial charge i in the solute according to the CHARMM-MERCK force fields, µnj is the lattice dipole j of the solvent at iteration n, µ0j is the permanent dipole moment of lattice point j calculated from eq 4, and the iteration control parameter, R, was chosen to be 0.1. Chemical Potential Energy Terms. We have used four contributions to the chemical potential change from solute substitution into the lattice: electrostatic, entropic, cavity, and repulsion. The first is obtained from direct electrostatic effects due to orientational changes of the solvent dipoles, the second

(42) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. J. Comput. Chem. 1983, 4, 187-217. (43) Halgren, T. A. J. Comput. Chem. 1996, 196, 490-615.

(44) Bo¨ttcher, C. J. F. Theory of Electric Polarization, 2nd ed.; Elsevier: Amsterdam, 1973. (45) Debye, P. Polar Molecules; The Chemical Catalog Co.: New York, 1929.

1246

Analytical Chemistry, Vol. 77, No. 5, March 1, 2005

peptide

qirij

i

rij3

E0j )

∑

(9)

µ1j ) µ0j E ˆ 0j [L(yL)]0 Enj

n

) -∇φ(rj) )

E0j

+

∑ j*k

(10)

3(µnj ‚rjk)rjk

-

r5jk

µn+1 ) Rµnj + (1 - R)µ0j E ˆ nj [L(yL)]n j

µnj r3jk

(11)

(12)

from the consequent entropic effect of dipole alignments with the solute field, the third from removal of the solvent dipoles at the solute lattice points, and the fourth from solute exclusion from the stationary phase. It is possible to consider dipolar effects from induced dipoles and from induced dipole-induced dipole interactions, but the forms of these are uncertain. and their influences are likely to be small due to cancellations in the replacement process. Electrostatic Term. The electrostatic chemical potential difference is a difference in energy from the permanent lattice dipoles from changes in the electrostatic field by insertion of the solute’s charges and dipoles is

∆Gelectro ) -

∑µ ‚E j

0

(13)

j

j

The dipole moments µj and fields E0j are those obtained after the n iteration of eq 11 and 12 converged to

|∆Gelectron - ∆Gelectron-1| ∆Gelectron

< 0.5%

(14)

Entropy Term. In the NVT (canonical) ensemble, the partition function of a dipole in an electrostatic field is46

Q(N,V,β) )

∫e

-U(E,µ)/kBT

dµ )

∫e

-yL

dy )

4π sinh(yL) yL (15)

with the resulting entropy being

|

∂ ln Q(N,V,T) ∂T

S ) kBT

+ kB ln Q(N,V,T)

(16)

∑∑ i

j

3kBT(4π0)2rij6

(18)

where 0 is the dielectric permittivity of vacuum. Repulsive Term. The repulsion energy due to overlap forces on a solute near the surface has been estimated by47 charges

∑

∆Uwall )

i

πFiσi12

45(zi + Sd)9

(19)

where Fi ) 21/2/σi3, σi is the van der Waals radius of solute atom i, and i is the energy coefficient. The values for σi were taken from Bondi49 while values for i were taken from Daragan and Mayo.50 (For the atoms C, H, N, O, and S, the values of i are 0.280, 0.176, 0.264, 0.289, and 0.482 kJ mol-1, respectively). Orientations and Averaging for Comparison with Experiment. To establish solute orientation, its coordinate axes are defined so that the z axis is from the center of geometry to the atom furthest from the center and the x and y axes are arbitrarily defined normal to z. Figure 2 shows an orientation in which z is tilted relative to the lattice coordinate Z. Different alignments may be examined to evaluate possible preferential orientations for solute adsorption. To obtain average values at the given Z location and different z-Z alignments, calculations were performed at different rotations about the z axis, since the solute partial charges and dipoles are not symmetrically positioned and thus are at different distances to the lattice dipoles at different rotational positions. The orientational averaging was done via

∑∆G

N,V

) kB[-yL coth(yL) - ln(yL) + ln sinh(yL) + ln(4π) +1]

2µ0,i2µ0,j2

Excl Incl

∆Gcav )

0 〈∆Gi,Z 〉)

0 i,Z

0 exp(- ∆Gi,Z /kT)

i

∑

(20) 0 exp(-∆Gi,Z /kT)

i

When the solute replaces the solvent, y changes from zero to the final value so the dipole entropy change is

∆Sd ) S - S0 ) kB{-yL coth(yL) - ln(yL) + ln[(sinh(yL)] + 1} (17) Cavity Term. The chemical potential change from replacement of solvent by solute is rigorously a complex function of size and interaction effects, especially in aqueous solvents.47 It is common to use van der Waals sizes and surface areas to estimate this term.48 However, this introduces additional parameters and is not always successful. Following Grossfield et al.,39 we have simplified this term to be a Keesom energy for interactions of permanent dipoles that are replaced by the solute (Excl) and those that remain (Incl). (46) McQuarrie, D. A. Statistical thermodynamics, 2nd ed.; University Science Books: Mill Valley, CA, 1973. (47) Israelachvili, J. Intermolecular and Surface Forces, 2nd ed.; Academic Press: San Diego, CA, 1991. (48) Levy, R. M.; Zhang, L. Y.; Gallicchio, E.; Felts, A. K. J. Am. Chem. Soc. 2003, 125, 9523-0530.

0 where ∆Gi,Z ) [∆Gcav + ∆Gelectro - T∆S0 + ∆Uwall]i,Z. 0 calculations were done at 14 unique In the present case, ∆Gi,Z orthogonal alignments of each of the solute Cartesian coordinate vectors ex, ey, ez at angles of nπ/4 (n ) 0, 1, 2) radians relative to the Z-lattice axis. This was intended to provide measures for all orientations of the solute, as indicated in Figure 2. At each of these 14 orientations, calculations were made at 8 unique rotational positions about the z-axis at every π/16 radians, some of which are shown in Figure 2. The 4-fold symmetry of the cubic lattice representation can be used to obtain average values for all other rotations. Thus, a total of 112 orientations were examined. Tests were done for several peptides with the number of z-Z alignments ranging from 1 to 10 in order to see where the simulated adsorption equilibrium constants reach a plateau and orientationally averaged results did not change significantly with additional alignments. The number of eight rotational positions about the z-axis was chosen as the minimum number needed and more orientations would not change the results.

(49) Bondi, A. J. Phys. Chem. 1964, 68, 441-450. (50) Daragan, V. A. M., K. H. Prog. NMR Spectrosc. 1997, 31, 63-105.

Analytical Chemistry, Vol. 77, No. 5, March 1, 2005

1247

Figure 2. Example alignment of peptide solute x, y, and z internal axes relative to solvent dipole lattice X, Y, Z axes.

Figure 4. Favorable (a) and unfavorable (b) orientations of VY peptide toward a hydrophobic surface at left.

Figure 3. Chemical potential profiles for the 14 unique alignments of the VY protected peptide with the lattice Z axis.

As an example of orientation dependence, Figure 3 shows Z variations of the chemical potential in the 14 orientations for the peptide VY. Several orientations have a minimum with negative 0 ∆Gi,Z as a function of Z, while others do not. Orientations 2 and 13 show extremes of unfavorable and favorable behavior, respectively. Figure 4 shows molecular structures for these orientations, (a) for 13 and (b) for 2. The difference is the nonpolar surface and polar solvent environments of the amino acids tyrosine (polar) and valine (nonpolar). 1248

Analytical Chemistry, Vol. 77, No. 5, March 1, 2005

The adsorption chemical potential at a given location Z, ∆∆G0Z, is the difference in the orientation averages of eq 20 at a position Z and at Z ) ∞ (in only solvent)

∆∆G0Z ) 〈∆Gi,Z0〉 - 〈∆Gi,∞0〉

(21)

0 This averaging means that the most negative values of ∆Gi,Z are weighted the greatest while positive values are insignificant. Examination of detailed results suggests strategies that might be implemented to save considerable computer time by ignoring negligible orientations.

The standard state of µ ) 0 for the solute environment exists at both the substrate surface and isolated in the ideal gas state. Thus, the quantity ∆∆G0Z is the result of a process of taking the orientationally averaged solute from the solvent to the ideal gas and then averaging the results of inserting it in different orientations of the coordinate Z. This process is also equivalent to adsorption of a peptide that can assume all orientations in both the mobile and stationary phases.51 There are differences in the chemical potential values in solvent only that are associated with rotation about the z axis, though there is little variation with z-Z alignment. Since neither rotation nor orientation should affect the values, these differences are apparently artifacts of the lattice system. However, with a sufficient number of rotational values, the average should be reliable, as we found. The function of eq 21 shows minimums at some position Zmin; this value can be compared to the experimental numbers obtained via eq 3. Comparisons were made with ∆∆G0Z at its minimum for different lattice dipole differences, ∆µ, between the solvent and surface regions. The value of ∆µ ) 1.055 gave agreement with the experimental value for adsorption of VY (index 5). Once this value was fit, the solvent and surface environments were kept the same in all calculations; only the peptide representation via the CHARMM force field varied in order to assess the predictive value of the ALLD method. Technical Details of Simulation Calculations. The algorithm for the above simulation method was developed in C++ with a programming style that was object-oriented. The development of the code with such a class hierarchy has important advantages in maintenance, inheritance, improvement, and extension. The debugger used was Microsoft Visual Studio.NET C++ and the Microsoft Developer Network library http:// msdn.microsoft.com/. The simulation runs were done with Microsoft C++.NET 2003 (Windows XP OS), Intel C++ (Red Hat Linux OS), and The Portland Group C++ (Red Hat Linux OS) compilers. The hardware architecture was an Intel Xeon dual CPU at 2.8 GHz and a Linux cluster of 1.53 GHz AMD Athlon K7 MP CPUs. The runs for a specific case took a few hours depending on the system and method. Simulation errors from the finite-sized lattice with limited number of the dipoles and the convergence scheme is estimated to be less than 1%. RESULTS Chemical potential profiles and averaged adsorption free energy changes have been obtained with the methodology described above for all nine peptides. These results are compared with experiment. Predicted chemical potentials for adsorption. Figure 5 0 shows adsorption chemical potential profiles, 〈∆Gi,Z 〉, along with 0 〉, for the values at separations representing only solvent, 〈∆Gi,∞ nine peptides. All show a minimum between 2.0 and 2.4 nm from the solid surface, the location depending on the solute size and force field partial charges. The minimum is away from the wall because the polar portions of the peptides prefer regions of larger dipoles, while the nonpolar regions prefer to be closer to the surface. Note that the ordering of chemical potentials here is not (51) Chen, S. F.; Liu, L. Y.; Zhou, J.; Jiang, S. Y. Langmuir 2003, 19, 28592864.

Figure 5. Profiles of angle-averaged chemical potentials of nine peptides from the surface to the mobile phase (values at 4.8 nm calculated without surface).

the same as the adsorption order because it is the difference 0 , which is compared to the between the two values, ∆∆GZmin 0 0 experimental value, ∆G . Figure 6 shows ∆∆GZmin results for the nine peptides along with comparisons to experiment. Values are given in Table 2. The errors are not much larger than the experimental uncertainty. The retention order is generally preserved except for peptides 1 and 2, the least strongly adsorbed species, which are predicted to be too strongly adsorbed relative to peptide 3. This may be due to insufficient configurations for appropriate averaging, by underestimating the cavity effect with eq 18, or by ignoring other effects such as van der Waals forces. Temperature and Solvent Dipole Dependence of Chemical Potential Differences for VY. To examine the influence of temperature and solvent dipole on the results, we have calculated 0 ∆∆GZmin for VY at 293, 298, 313, and 328 K with solvent dipole of 0 1.055. We find ∆∆GZmin is linear in temperature. From the slope and the intercept of this plot, the energy change of adsorption is constant at ∼-24 kJ mol-1 and the entropy change is ∼-52 J mol-1 K-1. The thermodynamic interpretation of this is that the retention is exothermic and enthalpy driven while there is Analytical Chemistry, Vol. 77, No. 5, March 1, 2005

1249

Table 2. Comparisons of Predicted ∆G0 Values from Descriptors and Hydrophobicity Scales molecular descriptors ∆G0exp

ALLD

-5.0 (1)* -5.6 (2) -5.7 (3) -7.2 (4) -8.5 (5) -8.8 (6) -8.9 (7) -11.4 (8) -12.5 (9)

-1.3a (2)* -1.4 (3) -0.2 (1) -0.4 (4) 0.0 (5) -0.2 (6) -0.7 (7) 0.3 (8) 1.0 (9)

log

P56,b

log

P57,b

3.1a (2)* -1.0a (4)* 4.1 (1) 0.3 (3) 2.1 (3) 1.8 (1) 1.9 (5) 2.7 (2) 0.0 (8) 0.0 (7) 2.6 (6) -0.6 (8) 5.0 (4) 1.6 (6) 4.2 (7) 5.0 (5) 3.8 (9) 2.6 (9)

hydrophobicity scales

dipolec

HASA/fd

HASA/Ne

ASA58,f

CE59,g

0.1a (1)* -3.5 (6) -3.0 (4) -5.8 (7) 0.0 (3) 0.0 (5) 3.9 (2) -3.3 (9) -1.9 (8)

-3.4a (4)* -2.9 (8) -2.0 (2) -1.1 (3) 0.0 (7) 0.4 (6) 1.9 (1) 3.0 (5) 3.7 (9)

-2.4a (4)* -2.2 (7) -1.0 (2) -1.9 (9) 0.0 (8) 1.1 (6) 3.3 (1) 4.1 (3) 5.0 (5)

-1.6a (2)* -1.1 (3) -0.7 (1) -0.1 (4) 0.0 (6) -1.4 (8) 0.4 (5) -2.2 (9) 2.5 (7)

-5.0a (4)* -4.4 (5) -3.7 (3) 2.4 (1) 0.0 (2) -4.8 (6) -9.8 (9) -2.2 (7) -5.8 (8)

P55,h

RPLC60,i

RPLC21,i

-0.8a (1)* -0.2a (2)* 0.2a (1)* -0.9 (3) 0.3 (3) -0.8 (3) -0.2 (2) 0.2 (4) -4.4 (5) -1.2 (4) 3.0 (1) 1.5 (2) 0.0 (5) 0.0 (6) 0.0 (4) -4.1 (7) -2.1 (9) -11.9 (9) 0.0 (6) 1.6 (5) -3.7 (7) -2.0 (9) 1.3 (7) -0.8 (6) -0.5 (8) 2.4 (8) -5.3 (8)

* (Order of adsorption strength). a (∆G0pred - ∆G0exp), kJ mol-1, uncertainties ( 0.05-0.10 kJ mol-1. b log P (octanol/water partition coefficient, including implicit hydrogens).56,57 c Dipole (peptide dipole moment). d HASA/f, (Hydrophobic accessible surface area fraction). e HASA/N (hydrophobic accessible surface area divided by the number of atoms). f Accessible surface area (ASA).58 g Contact energy (CE).59 h Partitioning (P).55 i Chromatographic techniques (RPLC).21,22,60

Figure 6. Comparison of simulation (filled symbols) and experimental results (open symbols) of adsorption chemical potential differences. Experimental uncertainties within the size of the symbols.

Figure 7. Profiles of angle-averaged chemical potential differences of VY peptide as a function of solvent dipole moment.

negligible heat capacity change.52 Both the electrostatic and the cavity terms include contributions from enthalpy and entropy as evident from the temperature dependence in their analytical equations, eqs 13 and 18. 0 The influence of solvent dipole on ∆∆GZmin is indicated in Figure 7 for VY. When the dipole is lowered, the adsorption is 0 less strong, with high sensitivity of ∆∆GZmin to ∆µ. This suggests that the model can also predict desorption as a less polar solvent responds less attractively to the solute field, as would occur with added cosolvents or added salt, as used for elution in hydrophobic chromatography. This aspect is being explored. DISCUSSION Contributions to Chemical Potential Changes and Orientational Averaging. Figure 8 shows Z profiles of the four contributions to the chemical potential (angle-averaged) in eq 23 for the VY peptide. All contributions but the repulsion are asymptotic to zero at the substrate surface and to fixed values at large separations (only solvent environment). While repulsion alone might determine the closest approach to the substrate, all (52) Dill, K. A.; Bromberg, S. Molecular Driving Forces; Garland Science, Taylor & Francis: New York, 2003.

1250 Analytical Chemistry, Vol. 77, No. 5, March 1, 2005

Figure 8. Angle-averaged contributions to computed chemical potential difference for the VY peptide. (0) Electrostatic; (4) entropic; (]) cavity; (×) repulsion; (O) total.

of the other terms actually establish the minimum location and value. The magnitude of the electrostatic term is largest, while the range of variation of the cavity term is the longest. The relative contributions for the entropy, cavity, and electrostatic terms at the calculated minimum positions of all nine peptides are shown in Figure 9. While the entropy terms make

Figure 9. Fractional contributions of absolute values of entropy (>0, eq 17), cavity (>0, eq 18), and electrostatic, (