J. Phys. Chem. B 2005, 109, 15145-15149
15145
Accurate Determination of Pyridine-Poly(amidoamine) Dendrimer Absolute Binding Constants with the OPLS-AA Force Field and Direct Integration of Radial Distribution Functions Yong Peng and George A. Kaminski* Department of Chemistry, Central Michigan UniVersity, Mount Pleasant, Michigan 48859 ReceiVed: March 8, 2005; In Final Form: June 3, 2005
OPLS-AA force field and direct integration of intermolecular radial distribution functions (RDF) were employed to calculate absolute binding constants of pyridine molecules to amino group (NH2) and amide group hydrogen atoms in and first generation poly(amidoamine) dendrimers in chloroform. The average errors in the absolute and relative association constants, as predicted with the calculations, are 14.1% and 10.8%, respectively, which translate into ca. 0.08 and 0.06 kcal/mol errors in the absolute and relative binding free energies. We believe that this level of accuracy proves the applicability of the OPLS-AA, force field, in combination with the direct RDF integration, to reproducing and predicting absolute intermolecular association constants of low magnitudes (ca. 0.2-2.0 range).
I. Introduction The OPLS-AA force filed has been successfully used in a wide variety of molecular simulations. This includes calculations of binding constants for systems ranging from small molecules to protein-ligand complexes.1 Usually, the statistical perturbation theory is used to compute changes in free energy as one ligand or guest molecule is slowly changed into another one.2 The resulting free energy change is then translated into the relative association constant for the pair of ligands involved. Although this method is capable of yielding very accurate results, successful application of this technique requires that the association constants for both the initial and the final systems are high enough for the complexes to stay together all the time during the course of the simulations. Therefore, the method in its original form is not applicable to calculation of binding constants of magnitudes close to one, when the complexes can stay dissociated for a significant fraction of the time. One the other hand, there are chemically significant cases when predicting association constants of low magnitudes is highly desirable. Complexes involving dendrimers represent one of such cases. They are distinguished from other macromolecules by their unique and well-defined structure derived from the dendritic topology rather then the traditionally recognized linear, branched, and cross-linked architecture. This highly branched structure features three architecture components. First a core is set at the center of the dendrimer. Radially attached to the central core are a number of monomeric repeating units, their length defining the dendrimer generation. The exterior part is at the outermost position along the chain and bears functional groups, e.g., NH2.3 Given the variety of possible terminal groups, numbers of generations, and core types, the dendrimers present a wide range of conformations and physical properties. This enables applications in the fields of photo- and electroactive materials,4 bioactive agents,5 catalysis,6 encapsulation and delivery,7 sensors/diagnoistics,8 and membrane chemistry.9 A significant * To whom correspondence should be addressed. E-mail address:
[email protected].
Figure 1. Generation 0 PAMAM dendrimer.
biological interest in the encapsulation of molecules by dendrimers is caused by their great potential in the area of membrane transport and drug delivery.10-12 Studies show that guest molecules can bind to a dendrimer by either hydrophobic or hydrophilic interactions.13 On the other hand, in the case of dendrimers with highly dense-packed terminal groups, encapsulation can occur as the guest molecules is trapped inside a “dendritic box” (or the space between the core and the surface of a dendrimer).14 Therefore, obtaining detailed microscopic information about the binding mechanism and relative association constants for different binding sites of dendrimers is important. The dendrimers studied in the presented work are of the poly(amidoamine) or PAMAM Starburst type, generations 0 (G0) and 1 (G1), as shown in Figures 1 and 2. The core that characters a PAMAM dendrimer is ethylenediamine (EDA). The branches are formed from N-(2-aminoethyl) acrylamide units. The PAMAM has been selected for this study because of the high affinity to external agents with complementary hydrogen bonding moieties. Presented below are results of simulations of the PAMAM dendrimers in complexes with pyridine, as a model for biologically significant interactions. The rest of the article is organized as follows. Described in section II is the methodology used in the calculations. In section III, we present results of the calculation and discuss their significance. The conclusions are given in section IV. II. Computational Methods A. OPLS-AA Force Field. We used the standard OPLSAA force field.1a,15,16 This force field has been shown to allow
10.1021/jp0511956 CCC: $30.25 © 2005 American Chemical Society Published on Web 07/16/2005
15146 J. Phys. Chem. B, Vol. 109, No. 31, 2005
Peng and Kaminski
Figure 2. Generation 1 PAMAM dendrimer.
a high accuracy in a wide variety of calculations, including small organic molecules (and explicitly pyridine) as well as peptides and proteins,1 which have significant structural similarity to the PAMAM dendrimers. Since the formalism of the OPLS-AA force field has been presented in detail elsewhere,1a we will limit ourselves to a brief description. The total energy expression is composed of the bond stretching and angle bending components Ebond and Eangle, the nonbonded term Enb and the torsional energy Etorsion. Both bond stretching and angle bending terms have the following harmonic form:
Ebond )
∑i Kb,i(ri - r0,i)2
(1)
Eangle )
∑i Kθ,i(θi - θ0,i)2
(2)
Here ri and θi, and r0,i and θ0,i denote the actual and equilibrium values for bond i and valent angle i, respectively. Kbi, and Kθi stand for the force constants. The nonbonded term includes pairwise Coulomb and Lennard-Jones intra- and intermolecular interactions between atoms i and j as shown in eq 3:
Enb )
∑i ∑j {qiqje2/rij + 4ij[(σij/rij)12 + (σij/rij)6]}fij
(3)
Geometric combining rules for the Lennard-Jones parameters were employed, with σij ) (σiiσjj)1/2 and ij ) (iijj)1/2. Scaling factor fij equals to zero when the interaction sites belong to the same valence bond or angle (1,2- and 1,3-interactions). When the atoms are separated by three bonds and belong to the same dihedral angle (1,4-interaction), fij ) 0.5. For all the other cases fij equals 1.0. Finally, the torsion term is given by a Fourier series using the dihedral angles φi and coefficients V1, V2, and V3:
Etorsion )
V2i V1i [1 + cos(φi)] + [1 - cos(2φi)] + 2 2 V3i [1 + cos(3φi)] (4) 2
∑
Once again, all the OPLS-AA force field parameters were adopted without modification. B. Liquid Simulations. Monte Carlo simulation of generations 0 and 1 PAMAM dendrimers were carried out with BOSS Version 4.3 program16 on 1.85 GHz Xeon workstations in our laboratory. Geometries for the solutes were taken from gas-
phase energy minimizations of the systems, also conducted with the BOSS software. Condensed state simulations were run in a chloroform cluster with one PAMAM dendrimer molecule as a solute and a number of pyridine molecules to create a ca. 1.0 M concentration of pyridine (close to the experimental conditions17). Spherical chloroform clusters were used to represent the solvation. The radii of the clusters were 31 and 35 Å for the G0 and G1 cases, respectively. While the main dendrimer-chloroform-pyridine simulations were performed using spherical clusters, the number of pyridine molecules to be included in order to achieve the desired pyridine solution concentration (close to the experimental conditions) was determined in a different way. To find the correct number of the pyridine molecules, we performed a simulation of pyridine solution in chloroform using a cubic box with periodic boundary conditions. The only components of this system were chloroform and pyridine molecules, with no dendrimer or other solute molecules present. The concentration of the pyridine molecules was determined by the ratio of the number of pyridine molecules to that of the chloroform. This ratio corresponds to the 1.0 M pyridine concentration, and it was then used for both the 31 and 35 Å clusters. A ca. 32 × 32 × 32 Å3 solvent box with periodic boundary conditions imposed was employed and the simulations were run at constant temperature of 298.15 K and pressure of 1 atm. (NPT ensemble). We found that in order to achieve the desired concentration of the pyridine molecules, we needed to use 245 molecules of chloroform and 22 pyridine molecules in a box with a total of 267 molecules. This Monte Carlo simulation consisted of 2 × 106 configurations of equilibration followed by 2 × 106 configurations of averaging used to determine the average volume and thus to ensure the correct concentration of pyridine. The target Monte Carlo acceptance ratio was 50%. Monte Carlo simulations of the PAMAM dendrimer and pyridine in the chloroform solution were carried out in the NPT ensemble at a pressure of 1 atm and a temperature of 298.15 K. The pyridine molecules were fully flexible. The ranges of the sampling over the rotational and translational movement of the solute and solvent molecules were adjusted to produce acceptance ratio of approximately 40% for the Monte Carlo configurations. Intermolecular nonbonded interactions were truncated at 11 Å with quadratic smoothing over the last 0.5 Å for the electrostatics. A standard correction was made to the energy for not including Lennard-Jones interactions beyond the cutoff radius.18 Each simulation consisted of at least 1.6 × 107 Monte Carlo configurations of equilibration followed by 4.2 × 107 configurations of averaging with Metropolis sampling for the generation 0 system and 6.0 × 107 configurations of
Simulating Pyridine Binding to PAMAM Dendrimers
J. Phys. Chem. B, Vol. 109, No. 31, 2005 15147
averaging for the generation 1. This ensured convergence of the radial distribution functions, and thus the association constants, as shown in the Results section below. C. Direct RDF Integration. Since the pyridine-dendrimer association constants are low, a significant number of the pyridine molecules is unbound during the simulations. Therefore, we propose using a method for computing the binding constants, which is different from the widely employed but not usable in this case statistical perturbation theory. The utilized technique is based upon integration of the intermolecular radial distribution functions (RDF). We are interested in computing radial distribution functions between the pyridine nitrogen atom and the hydrogen atoms at the terminal amino and more buried amide groups. In the canonical ensemble the radial distribution function for atoms 1 and 2 g(r1,r2) is defined as
g(r1,r2) )
N(N - 1) F2ZNVT
∫dr3 dr4 ... drN exp(- βE(r1,r2,...rN)) (5)
This function represents the probability density of finding the particle 1 at r1 and particle 2 at r2, at the same time. N and F stand for the total number of particles in the system and their number density, respectively. β ) 1/RT, where T is the temperature in degrees K. ZNVT is the partition function. Normally, this probability density depends on relative position of the atoms 1 and 2 b r )b r2 - b r1 and can be averaged to depend on the interatomic distance r ) |r b1 - b r2| only. Let us denote the RDF between the pyridine nitrogen atom and the -NH2 hydrogen as g1(r) and the RDF between the same nitrogen and the dendrimer amide hydrogen atom as g2(r). Then the total numbers of the corresponding N‚‚‚H pairs within a distance R1 from the hydrogens can be expressed as shown in eqs 6 and 7:
NNH2...pyridine ) NCONH...pyridine )
∫0
R1
∫0R
1
2
g1(r)F4πr dr
(6)
g2(r)F4π r2 dr
(7)
Here F is the density of the pyridine. Furthermore, one can obtain the total number of the pyridine molecules from radial distribution function by integrating to infinity.
Npyridine1 )
∫0∞ g1(r)F4π r2 dr
(8)
Npyridine2 )
∫0∞ g2(r)F4π r2 dr
(9)
Let us now introduce ratios X as XNH2 ) NNH2...pyridine/Npyridine and XCONH2 ) NCONH2...pyridine/Npyridine. These ratios can be computed from integration of the radial distribution functions as given in eqs 6-9. Then the total numbers of bonds per hydrogen atom for amino and amide groups will be given by ZNH2 ) XNH2N and ZCONH ) XCONHN, respectively, where N is the total number of the pyridine molecules (known from the computational setup). Therefore, the association constant can be estimated from the ratios of bound/unbound dendrimer amino or amide hydrogens, ZNH2/(1 - ZNH2) and ZCONH/(1 - ZCONH), if we divide them by the concentration of the unbound pyridine, which is computed as the difference between the total number of the pyridine molecules and the number of the bound ones,
Figure 3. Gas-phase optimized structure of the G0 PAMAM dendrimer with the longest distance between heavy atoms shown.
resulting in approximately 1 M in this work:
Kb(NH2) ) [ZNH2/(1 - ZNH2)]/1M
(10)
Kb(CONH) ) [ZCONH/(1 - ZCONH)]/1M
(11)
and for the free energies of association
∆Gb(NH2) ) -RT ln[Kb(NH2)]
(12)
∆Gb(CONH) ) - RTln[Kb(CONH)]
(13)
III. Results and Discussion Shown in Figures 3 and 4 are the lowest energy conformers obtained from the conformational search in gas phase for the and the first generation (G0 and G1) of the PAMAM dendrimers, respectively. It can be seen that these low-generation PAMAM structures are not as spherical as higher generation dendrimers. The largest atom-to-atom distances (for heavy atoms) are ca. 18 and 26 Å for the zeroth- and the first-order molecules. Therefore, once can expect that the pyridine molecules have access not only to the surface amino groups but also to the somewhat more buried amides. Since our method for computing the dendrimer-pyridine association constants depends on the accuracy of the computed radial distribution functions, it was crucial to ensure that the RDF’s were completely converged. Data in Figures 5-8 demonstrate that the convergence has, in fact, been achieved. Several conclusions can be drawn by looking at these radial distribution function graphs. First, the RDFs clearly converge completely by the end of the simulations (4.2 × 107 configurations of averaging for the G0 case and 6 × 107 configurations of averaging for the G1 one). Moreover, some curves are not even clearly visible on the graphs as they largely coincide with the final RDFs. The radial distribution functions for aminogroup hydrogen and the pyridine nitrogen seem to converge slower than those for the amide hydrogens, but even the former ones are indistinguishable from the final results after 2.8 × 107 and 4 × 107 configurations of averaging for the zeroth- and first-order denerations, respectively. Second, while the RDFs for the amide group contain just one first peak for the binding to the pyridine nitrogen, the NH2 groups clearly produce two
15148 J. Phys. Chem. B, Vol. 109, No. 31, 2005
Peng and Kaminski
Figure 6. H-N radial distribution functions for the G0 dendrimer NH2 hydrogens and pyridine nitrogen atoms. The legend shows the corresponding numbers of the Monte Carlo configurations of averaging.
Figure 4. Gas-phase optimized structure of the G1 PAMAM dendrimer with the longest distance between heavy atoms shown.
Figure 7. H-N radial distribution functions for the G1 dendrimer amide hydrogens and pyridine nitrogen atoms. The legend shows the corresponding numbers of the Monte Carlo configurations of averaging.
Figure 5. H-N radial distribution functions for the G0 dendrimer amide hydrogens and pyridine nitrogen atoms. The legend shows the corresponding numbers of the Monte Carlo configurations of averaging.
peaks. This is natural, as the second peak signifies binding to the neighboring hydrogen. However, the end of the first peak is at only about 2.75-3.00 Å distance, the peaks are well separated, and the first one can be used easily to determine the thermodynamic average of the number of the H‚‚‚N bonds with the pyridine molecules and the association constant and free energy using eqs 6-13 above. There is still a certain ambiguity as to at which distance should the R1 in the eqs 6 and 7 be placed (the cutoff distance for the bound hydrogen-nitrogen pair). We will discuss this issue below. Finally, the thickest lines on the graphs correspond to the longest simulations, and these are the radial distribution functions we used in the actual calculations of the association constants. The resulting absolute association constants and free energies are shown in Table 1. It can be seen that the computed results
are very close to the experimental ones. The absolute errors for the association constants for the zeroth-order PAMAM dendrimer are 14.5% and 8.3% for the pyridine binding to the amide and NH2 group hydrogens, respectively. This corresponds to the free energy errors of 0.08 and 0.05 kcal/mol, which means the presented technique is very accurate, and, which is also very important, that the OPLS-AA force field is very adequate in describing such dendrimer-guest molecule interactions. The errors for the G1 dendrimer are 8.5% and 25.0% for the association constants and 0.09 and 0.08 kcal/mol, which is also a very respectable result. If we consider the relative association constants (the ratio of the association constants of pyridine with the amide and amino hydrogens), the target experimental numbers for the zeroth and first generations are 3.64 and 5.12, respectively, and the computed ones are 3.85 and 4.42. The errors are thus 5.7% and 10.8%, or 0.03 and 0.06 kcal/mol. Therefore, the presented method permits to reproduce both the absolute and the relative association constants in a very good agreement with the experimental data. There is one more issue that should be discussed here. The cutoff distances R1 to define the bound hydrogen-nitrogen
Simulating Pyridine Binding to PAMAM Dendrimers
J. Phys. Chem. B, Vol. 109, No. 31, 2005 15149 reasonable values, the association energy is still in a good agreement with the experimental data. IV. Conclusions
Figure 8. H-N radial distribution functions for the G1 dendrimer NH2 hydrogens and pyridine nitrogen atoms. The legend shows the corresponding numbers of the Monte Carlo configurations of averaging.
TABLE 1: Comparison of Experimental and Computed Association Constants and Free Energies (in kcal/mol) for Binding between the PAMAM Dendrimer Hydrogens and Pyridine Moleculesa experiment
this work
Kb
∆G
Kb
∆G
R1, Åb
G0-amide G0-NH2 G1-amide G1-NH2
1.31 0.36 0.82 0.16
-0.160 0.606 0.118 1.087
1.504 0.392 0.893 0.202
-0.242 0.556 0.067 0.949
3.27 2.73 3.75 2.97
Experimental data from ref 17. b Integration limit in eqs 6 and 7.
TABLE 2: Computed Association Constants for Binding between the PAMAM Dendrimer Hydrogens and Pyridine Molecules
a
Acknowledgment. This work was supported by the U.S. Army Research Laboratory, Grant DAAD19-03-2-0012. References and Notes
pair with N
a
We have computed absolute association constants and free energies of association with pyridine for the and first generation poly(amidoamine) dendrimers in chloroform. OPLS-AA force field has been used. Since the association constants are rather low and have the order of magnitude of 1.0, we have employed the direct integration of intermolecular radial distribution functions (RDF) to obtain the values of the association constants. The results are in an excellent agreement with the available experimental data, with the accuracy in the association free energies being better than 0.1 kcal/mol. This demonstrates that the OPLS-AA force field is adequate and robust in such dendrimer-heterocycle simulations and that the direct RDF integration is viable technique in reproducing and predicting association constant values in the ca. 0.2-2.0 range. We hope that this work will contribute to the field of dendrimer research and will be useful in development of applications of dendrimer host-guest complexes.
H bound to pyridine N
R1, Åa
Kb
G0-amide G0-amide G0-amide G1-amide G1-amide G1-amide
3.03 3.27 3.39 3.15 3.39 3.75
1.312 1.504 1.614 0.644 0.756 0.893
Integration limit in eqs 6 and 7.
couples (the values employed are given in the last column of Table 1) were chosen to be at the minima of the radial distribution functions in Figures 5-8. Nevertheless, these values are somewhat arbitrary, and to ensure the robustness of the direct RDF integration technique, we have also studied the stability of the computed association constant values with respect to the choice of the cutoff distance. The results of these studies are given in Table 2. The minima in the RDFs for the amino group cases are more than clear. This is why we only used one R1 value for each calculation of the NH2 binding constant and do not present a table similar to Table 2. However, the amide group cutoff distances are chosen more arbitrarily. As can be seen from Table 2, the results are sufficiently stable with respect to the choice of R1. The ratios of the maximum and minimum association constants values are 1.23 and 1.39 for the G0 and G1 dendrimers, respectively. This translates into free energy differences of only 0.12 and 0.19 kcal/mol. Therefore, even if we push the limits and set the integration cutoff to some less
(1) See, for example: (a) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. J. Am. Chem. Soc. 1996, 118, 11225; (b) Tirado-Rives, J.; Orozco, M.; Jorgensen, W. L. Biochemistry 1997, 36, 7313; (c) Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. J. Phys. Chem. B 2001, 105, 6474. (d) Tominaga, Y.; Jorgensen, W. L. J. Med. Chem. 2004, 47, 2534. (2) (a) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420; (b) Price, D. J.; Jorgensen, W. L. J. Comput. Aid Mol. Des. 2001, 15, 681. (3) (a) Tomalia, D. A.; Naylor, A. M.; Goddard, W. A., III. Angew. Chem. 1990, 102, 119; (b) Tomalia, D. A.; Hedstrand, D. M.; Ferritto, M. S. Macromol. 1991, 24, 1435. (4) (a) Kim, E.; Kim, K.; Yang, H.; Kim, Y. T.; Kwak, J. Anal. Chem. 2003, 75, 5665. (b) Huang, B.; Tomalia, D. A. J. Lumin. 2005, 111, 215,. (5) (a) Khopade, A. J.; Caruso, F.; Tripathi, P.; Nagaich, S.; Jain, N. K. Int. J. Pharm. 2002, 232, 157; (b) Tomalia, D. A.; Huang, B.; Swanson, D. R.; Brothers, H. M.; II; Klimash, J. W. Tetrahedron 2003, 59, 3799. (6) Astruc, D.; Chardac, F. Chem. Rew. 2001, 101, 2991. (7) Morgan, M. T.; Carnahan, M. A.; Immoos, C. E.; Ribeiro, A. A.; Finkelstein, S.; Lee, S. J.; Grinstaff, M. W. J. Am. Chem. Soc. 2003, 125, 15485. (8) Gong, L. Z.; Hu, Q. S.; Pu, L. J. Org. Chem. 2001, 66, 2358. (9) Kovvali, A. S.; Sirkar, K. K. Ind. Eng. Chem. Res. 2002, 41, 2287. (10) Kukowska-Latallo, J. F.; Bielinska, A. U.; Johnson, J.; Spindler, R.; Tomalia, D. A.; Baker, J., Jr. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 4897. (11) Turro, C.; Niu, S.; Bossman, S. H.; Tomalia, D. A.; Turro, N. J. J. Phys. Chem. 1995, 99, 5512. (12) Jansen, J. F. G. A.; deBrabander, E. M. M.; Meijer, E. W. Science 1994, 266, 1226. (13) (a) Newkome, G. R.; Moorefield, C. N.; Baker, G. R.; Saunders, M. J.; Grossman, S. H. Angew. Chem. Int. Edit. Engl. 1991, 30, 1178,; (b) Hawker, C. J.; Wooley, K. L.; Frechet, J. M. J. J. Chem. Soc., Perkin Trans. 1993, 1, 1287; (c) Stevelmans, S.; van Hest, J. C. M.; Jansen, J. F. G. A.; van Boxtel, D. A. F. J.; de Brabander-van den Berg, E. M. M.; Meijer, E. W. J. Am. Chem. Soc. 1996, 118, 7398. (d) Watkins, D. M.; Sayed-Sweet, Y.; Klimash, J. W.; Turro, N. J.; Tomalia, D. A. Langmuir 1997, 13, 3136. (14) (a) Johan, F. G. A.; Jansen, E. W.; Meijer, E. W.; de Brabandervan den Berg, E. M. M. J. Am. Chem. Soc. 1995, 117, 4417. (b) Jansen, J. F. G. A.; de Brabander-van den Berg, E. M. M.; Meijer, E. W. Science 1994, 266, 1226. (15) Jorgensen, W. L.; McDonald, N. A. THEOCHEM 1998, 424,. (16) Jorgensen, W. L. BOSS Version 4.3. Yale University: New Haven, CT, 2001. (17) Santo, M.; Fox, M. A. J. Phys. Org. Chem. 1999, 12, 293. (18) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J. Am. Chem. Soc. 1984, 106, 6638.
