Computationally Efficient Methodology for Atomic ... - ACS Publications

Article pubs.acs.org/JPCB

Computationally Efficient Methodology for Atomic-Level Characterization of Dendrimer−Drug Complexes: A Comparison of Amine- and Acetyl-Terminated PAMAM Ariela Vergara-Jaque,† Jeffrey Comer,‡,§ Luis Monsalve,† Fernando D. González-Nilo,§,∥,‡ and Claudia Sandoval*,§ †

Center for Bioinformatics and Molecular Simulation, Universidad de Talca, 2 norte 685, Talca-Chile Fraunhofer Chile Research, Mariano Sánchez Fontecilla 310 piso 14, Las Condes, Santiago-Chile § Universidad Andres Bello, Center for Bioinformatics and Integrative Biology, Facultad de Ciencias Biológicas, Av. Republica 239, Santiago-Chile ∥ Centro Interdisciplinario de Neurociencia de Valparaíso, Universidad de Valparaíso, Gran Bretaña 287, Valparaíso-Chile ‡

S Supporting Information *

ABSTRACT: PAMAM dendrimers have been widely studied as a novel means for controlled drug delivery; however, computational study of dendrimer−drug complexation is made difficult by the conformational flexibility of dendrimers and the nonspecific nature of the dendrimer−drug interactions. Conventional protocols for studying drug binding have been designed primarily for protein substrates, and, therefore, there is a need to establish new protocols to deal with the unique aspects of dendrimers. In this work, we generate cavities in generation-5 polyamidoamine (PAMAM) dendrimers at selected distances from the center of mass of the dendrimer for the insertion of the model drug: dexamethasone 21phosphate or Dp21. The complexes are then allowed to equilibrate with distance between centers of mass of the drug and dendrimers confined to selected ranges; the free energy of complexation is estimated by the MM-GBSA (MM, molecular mechanics; GB, generalized Born; SA, surface area) method. For both amine- and modified acetyl-terminated PAMAM at both low and neutral pH, the most favorable free energy of complexation is associated with Dp21 at distance of 15−20 Å from the center of mass of the dendrimer and that smaller or larger distances yield considerably weaker affinity. In agreement with experimental results, we find acetyl-terminated PAMAM at neutral pH to form the least stable complex with Dp21. The greatest affinity is seen in the case of acetyl-terminated PAMAM at low pH, which appears to be due a complex balance of different contributions, which cannot be attributed to electrostatics, van der Waals interactions, hydrogen bonds, or charge−charge interactions alone.

■

INTRODUCTION

branched structures with a globular architecture. They may have hydrophobic or hydrophilic internal cavities, and their terminal groups can be modified with different molecular moieties in order to decrease cytotoxicity, improve bioavailability, and optimize transport of particular drugs.12−14 Dendrimers may encapsulate a drug through intermolecular interactions or, alternatively, form covalent bonds with a drug, and in this way enhance the uptake, cellular trafficking,15 and controlled release of the drugs. Polyamidoamine (PAMAM)10,11 is one of the most studied classes of dendrimers,16,17 and has widely been investigated for delivery of drugs targeting a diverse set of diseases.18 Studies have

A recent trend in the pharmaceutical industry has been the focus on improving the properties of drugs that are currently used in various treatments, rather than focus on the creation of new drugs. Side effects in most drugs are considerable, in some cases more severe than the malady the drug is chosen to treat. To provide a solution to this problem, diverse nanocarriers have been proposed, such as micelles,1 polymersomes,2−4 carbon nanotubes,5,6 nanospheres,7,8 and many others. These nanocarriers can reduce the side effects while improving targeting and controlled release all without changing the chemical structure of the drug.9 Dendritic nanoparticles are actually one of the most studied for use as nanocarriers,10,11 because these can enhance the transport and release of drugs due to their particular multifunctional structure. Dendrimers are synthetic hyper© 2013 American Chemical Society

Received: January 2, 2013 Revised: April 10, 2013 Published: May 3, 2013 6801

dx.doi.org/10.1021/jp4000363 | J. Phys. Chem. B 2013, 117, 6801−6813

The Journal of Physical Chemistry B

Article

needs to compare only the complexed state and the free state to obtain the binding free energy and sampling performed over the quasi-continuum of variously defined intermediate states is not needed. Furthermore, MM-GBSA has the advantage of permitting easy separation of various contributions to the free energy, such as the van der Waals, electrostatic, and solvation contributions. However, MM-GBSA is somewhat ad hoc due to the fact that while the configurational sampling is usually performed in explicit solvent MD, the free energy calculation relies on the approximate implicit-solvent Generalized Born method38 and a phenomenological treatment of hydrophobic effects.39 Similar statements hold when the Generalized Born implicit solvent is replaced with the formally more accurate implicit solvent of the Poisson−Boltzmann equation. In practice, the so-called MM-PBSA method, which differs from MM-GBSA in that numerical solution of the Poisson− Boltzmann equation replaces the Generalized Born formalism, does not necessarily provide better qualitative rankings of macromolecule−ligand systems.40 Thus, with the aim of describing dendrimer−drug complexation at the atomic level, in this work we develop a computational protocol to characterize the structural and energetic properties that facilitate the encapsulation of drugs, such as Dp21, by unmodified and acetylated PAMAM-G5 dendrimers (Figure 1) at different protonation states. We use

shown that functionalization of the terminal groups of PAMAM by acetyl moieties reduces toxicity,19,20 while conferring an improvement in the rate of intracellular uptake,21 because it diminishes the electrostatic interactions between the dendrimer and the phospholipid headgroups in the cell membrane.22−25 Dendrimer−drug noncovalent complexes are established mainly by electrostatic, hydrogen bonding, and hydrophobic interactions.26 The great versatility of PAMAM dendrimers to encapsulate and interact through electrostatic interactions makes these nanoparticles interesting as drug delivery systems.27,28 Several studies have been carried out on the transport and delivery of different model drugs with dendrimers, but molecular level descriptions on the nature of specific interactions or encapsulation mechanism using computational methodologies are few. Tanis and Karatasos26 carried out atomistic molecular dynamics (MD) simulations of generation-3 (G3) PAMAM dendrimers and ibuprofen. This study revealed that dendrimer−drug complex formation depends on protonation state of both PAMAM dendrimer and the drug. Considering the previous examples, protonation states and functionalization of the dendrimer’s terminal groups are important parameters to take into account in drug delivery applications. Yang et al.22 experimentally evaluated encapsulating Dexamethasone 21-phosphate (Dp21) in generation-5 fully acetylated and amine-terminated PAMAM, finding that the amine-terminated dendrimer can host Dp21 at both acidic and neutral conditions, while the acetylated dendrimer only encapsulates it under acidic conditions. Moreover, the orientation of Dp21 in the dendrimer cavities depends on the degree of protonation of the amine groups of the dendrimers and the phosphate group of the drug. Atomistic MD simulation affords the ability to estimate free energies while simultaneously allowing one to determine the atomic scale interactions that contribute to these energies. For example, Maingi et al.29 used umbrella sampling to calculate potentials of mean force for several drugs as a function of distance between the center of mass of the dendrimer and the center of mass of the drug. These calculations showed the importance of the position of the drug within the dendrimer, as large differences in the binding energy (tens of kcal/mol) with respect to different drug positions within the dendrimer were observed. However, presumably due to computational expense that the additional umbrella sampling windows would have required, the potentials of mean force were not continuously extended to large enough distances that the dendrimer and drug were effectively in their free states; thus, the binding free energies could not be clearly determined.30 Similarly, most rigorous free energy methods, including umbrella sampling followed by the weighted histogram analysis method,31,32 adaptive force biasing,33,34 metadynamics,35 and free-energy perturbation,36 require quasi-continuous sampling between all energetically relevant complexed states and the free state for the binding free energy to be determined. For example, to determine the binding free energy between a dendrimer and drug by umbrella sampling, one must restrain the distance between the two molecules to many different values to obtain a continuous set of well-sampled overlapping distributions from separations of a few angstroms to separations where the two are far enough apart to be effectively free. The MM-GBSA (MM, molecular mechanics; GB, generalized Born; SA, surface area) method,37 on the other hand, can be much cheaper than more rigorous methods in that one

Figure 1. Model of (A) amine and (B) acetyl-terminated PAMAM dendrimer in a box containing explicit water molecules and 0.15 M NaCl. The dendrimers structures are shown by a space-filling representation with carbon, nitrogen, and oxygen atoms shown in green, blue, and red, respectively. Acetyl groups are shown in purple (right panel). The Na+ and Cl− ions are depicted as yellow and cyan spheres. The chemical structure of Dexamethasone 21-phosphate disodium (Dp21) is shown in detail in the center of the figure. The encapsulation of Dp21 was simulated in both systems at low and neutral pH.

the MM-GBSA method to obtain semiquantitative comparisons of the free energies of dendrimer−drug complexes at three selected ranges of the dendrimer−drug distance, in order to capture some of the position-dependence of the drug affinity within the dendrimer. Our results are consistent with existing experimental data.

■

COMPUTATIONAL METHODS For this study, the interaction of the drug Dp21 with amineand acetyl-terminated PAMAM dendrimers of the fifth generation at low and neutral pH was characterized through MD simulations. As a first step, we describe the construction of the atomistic models of the dendrimers, together with the protocol used to equilibrate each structure. Next, construction of cavities in the dendrimers and the insertion of the drug by docking molecular are detailed. Finally, we present the 6802



Article

Table 1. Chemical and Physical Properties of the Simulated Complexes complexes PAMAM-G5/Dp21 (Low pH)

PAMAM-G5/Dp21 (Neutral pH)

PAMAM-G5-Ace/Dp21 (Low pH)

PAMAM-G5-Ace/Dp21 (Neutral pH)

interaction point core middle surface core middle surface core middle surface core middle surface

acetyl groups

dendrimer charge (e)

0 0 0 0 0 0 128 128 128 128 128 128

+254 +254 +254 +128 +128 +128 +126 +126 +126 0 0 0

Dp21 charge (e)

no. water molecules

−1 −1 −1 −2 −2 −2 −1 −1 −1 −2 −2 −2

36872 39650 37855 28156 25479 26428 31037 25308 25835 22305 23575 27211

no. sodium ions 107 115 110 82 74 77 90 74 75 65 68 79

no. chloride ions 360 368 363 208 200 203 215 199 200 63 66 77

piston method.48 Short-range electrostatic and van der Waals interactions were truncated smoothly from 10 to 12 Å, while long-range electrostatic forces were computed using the particle-mesh Ewald (PME) summation method.49 All MD simulations were performed using the program NAMD2.850 with the standard ions and water of the CHARMM2742 force field, and analyzed with VMD51 software. The potential energy of each system was monitored throughout the simulation in order to determine the time necessary to reach a stable state. Data collected along the trajectories were used to calculate molecular properties such as: radius of gyration, asphericity, solvent-accessible surface area (SASA), and the radial distribution function (RDF) of ions. Drug Insertion into the Cavities. For the purpose of generating binding sites to encapsulate Dp21 within the dendrimers, we built empty cavities within the amine and acetyl-terminated PAMAM dendrimers (at low and neutral pH) beginning with their equilibrated conformations. MD simulations, incorporating a tclForces script designed to apply forces to a small number of preselected atoms, were carried out to generate cavities of the appropriate size (12 Å) to encapsulate the Dp21 in different regions of the dendrimers. The cavities were formed at three distinct distances from the center of the dendrimers: (i) around the core, (ii) on middle region, and (iii) on the outer surface of the dendrimer, yielding 12 separate systems. The angular positions of the centers of the cavities were chosen at random. Then, we performed docking simulations using the program AutoDock 4.052 to insert the drug into the cavity of each of the twelve systems and to find the preferred binding conformations of Dp21. Dp21 structure was parametrized considering pKa1 = 1.9 and pKa2 = 6.41.22 The parameters were obtained from the ParamChem online server, using the CGenFF force field.53 Therefore, the charges of the drug and the dendrimers were obtained from ParamChem and Paratool,43 respectively. In preparation for docking, a grid box was built in each cavity large enough to accommodate free motion of the drug (60 × 60 × 60 Å3). The grid parameters were generated using AutoGrid 4.052 and the Lamarckian genetic algorithm was used to perform a search of the conformational space of the drug. Ten conformations of the drug were generated for each of the 12 systems, which were finally analyzed with the program AutoDockTools54 to select the best scoring docked conformation. Molecular Dynamics Simulation of Dendrimer−Drug Complexes. The lowest energy conformation of each complex obtained from the docking calculations was used as starting

methodology used to simulate the dendrimer−drug complexes, including the MM-GBSA method used to evaluate the binding free energy. Description of the Dendrimer Models. The initial structure for each dendrimer used for this study was built with ICM software,41 employing an in silico protocol similar to the divergent synthesis of dendrimers (i.e., the dendrimer is assembled from a multifunctional core). Atomistic models of ethyldiamine-core and NH2-terminated G5 PAMAM dendrimers were built in two distinct protonation states: all amines protonated (low pH) and only primary amines protonated (neutral pH). Acetylated-terminated (100%) G5 PAMAM dendrimers were generated under the same conditions. It was necessary to parametrize each structure because the dendritic structures are not included in conventional force fields. As strategy of parametrization, the dendrimer was divided into three segments: the core, intermediate dendrons, and the outer surface. Each segment was parametrized separately and then the structures were assembled. The parameters were constructed based on the CHARMM42 force field using the Paratool plugin.43 The molecular geometry of each segment was optimized using DFT theory with the Hartree−Fock functional and 6-31+G* basis set.44 The force constants for bond length, bond angles, dihedral angles, and improper angles were determined from single point calculations and by comparison to atoms with similar chemical environments in the CHARMM2742 force field. The restrained electrostatic potential (RESP) method45 was used to fit the electrostatic potential and to get atomic partial charges. All quantum chemical calculations were carried out using the Gaussian 03 software package46 and then imported into Paratool.43 Molecular Simulation of PAMAM Dendrimers. The dendritic structures parametrized above (PAMAM-NH2 and acetylated PAMAM at low and neutral pH) were minimized and equilibrated through MD simulation. As starting point, the dendrimers structures were solvated in a periodic box (104 × 110 × 100 Å3) of pre-equilibrated TIP3P47 water molecules. Sodium and chloride ions (0.15 M NaCl) were added to the aqueous phase to mimic physiological conditions and to guarantee electric neutrality. The initial configuration of the systems (Figure 1) was optimized by means of an energy minimization algorithm, followed by an equilibration and relaxation in a 15 ns MD simulation at 310 K in the isobaric− isothermal ensemble. Constant temperature was enforced using a Langevin thermostat with a damping coefficient of 1 ps−1. Constant pressure (1 atm) was enforced using the Langevin 6803



Article

Figure 2. (A) Radius of gyration and (B) SASA as a function of simulated time for four systemsamine-terminated PAMAM at low (Low_pH) and neutral pH (Neutral_pH), and acetylated PAMAM at low (ACE_Low_pH) and neutral pH (ACE_Neutra_pH).

equilibrium of each simulated complex, was calculated the radial distribution of amine groups around the drug (see Supporting Information). MM-GBSA Calculations of Dendrimer−Drug Complexes. To estimate the relative binding free energy of the dendrimer−drug complexes we employed the MM-GBSA (MM, molecular mechanics; GB, generalized Born; SA, surface area) method.37 In preparation for the MM-GBSA calculations, snapshots of the system configuration were extracted from MD trajectories at 10 ps intervals, and the explicit water molecules and ions were removed. The MM-GBSA analysis was performed on three subsets of each system: the dendrimer alone, the drug alone, and the complex. For each of these subsets, the free energy was calculated as follows:

point to evaluate the specific interactions between dendrimers and the drug in MD simulations. Each dendrimer−Dp21 complex was solvated in a cubic box of TIP3P water molecules and NaCl (0.15 M) to mimic physiological conditions. A total of twelve simulations were carried out, with three different initial positions of the Dp21 within amine- or acetyl-terminated PAMAM having the protonation states typical of low or neutral pH. A summary of the properties of these systems is given in Table 1. The simulations protocols were the same as those described in the Molecular Simulation of PAMAM Dendrimers section. It is not possible obtain adequate sampling of Dp21 configurations throughout the entire PAMAM-G5 dendrimer in simulations shorter than several hundred nanoseconds due to the relatively low diffusivity of Dp21 within the dendrimer and time scale of conformational changes for PAMAM-G5. Therefore, for subsequent calculations of binding free energy obtaining converged free energy values, we chose to restrict the sampling to certain ranges of distances between the centers of mass of the drug and dendrimer. Restricting the range of the reaction coordinate is a common strategy to accelerate convergence in free energy calculations.36 Here, we applied commonly used boundary potentials to restrict the sampling to three different regions: 0 to 10 Å, 15 to 20 Å, and 25 to 30 Å, henceforth called the “core”, “middle”, and “surface”, respectively. Specifically, the following boundary potential was applied as a function of the distance, d, between the center of mass of the dendrimer and the center of mass of Dp21 using the Colvars module34 of NAMD ⎧1 2 ⎪ K (d − d0) , d < d0 2 ⎪ ⎪ V (d ) = ⎨ 1 2 ⎪ 2 K (d − d1) , d > d1 ⎪ ⎪ otherwise ⎩ 0,

Gtot = HMM + Gsolv − pol + Gsolv − np − T ΔSconf

(2)

where HMM, Gsolv−pol, Gsolv−np, and Sconf corresponded to the sum of bonded and Lennard-Jones energy terms, the polar contribution of solvation energy, the nonpolar contribution to the solvation energy, and the conformational entropy, respectively. Both HMM and Gsolv−pol were calculated by using NAMD 2.850 with the generalized Born implicit solvent model55 and the parameters described above. The dielectric constant of the solvent was set 78.5. Gsolv−np was calculated as a linear function of the solvent-accessible surface area (SASA), which was calculated with a probe radius of 1.4 Å.39 The binding free energy of amine- and acetyl-terminated G5 PAMAM dendrimers complexed with Dp21 was calculated by the difference ΔG bind = Gtot(dendrimer − drug) − Gtot(dendrimer) − Gtot(drug)

(1)

(3)

where Gtot values were the averages over the simulation. The initial portions of the ∼48 ns simulations were excluded to permit ΔGbind to approach a steady value. The change in the conformational entropy should, in principle, be included; however, it was not computed, due to the large computational cost. Furthermore, due to the fact that we employed the singletrajectory approach,37 the bonded contribution to ΔGbind was identically zero. Finally, we compared the relative binding affinities to experimental results and performed structural analyses for each system in study.

2

where K = 20 kcal/(mol Å ) and d0 and d1 were chosen in a manner consistent with the definitions of the regions above. The quoted free energy values can only be seen to apply to intervals that were sampled and are not necessarily representative of the binding free energy for the complex as a whole (see Supporting Information). It is possible that important regions of the reaction coordinate, including perhaps the value with the lowest free energy, has not been sampled. Finally, to appraise conformational stability and thermodynamic 6804



Article

Table 2. Calculated Average Values of the Principal Moments of Inertia, Aspect Ratios, and Asphericity of Amine- and AcetylTerminated PAMAM Dendrimers at Low and Neutral pHa principal moments of inertia Iz (kDa nm2)

dendrimers Amine-terminated PAMAM Acetyl-terminated PAMAM

low pH neutral pH low pH neutral pH

161 134 190 210

± ± ± ±

12 8 10 5

Iy (kDa nm2) 129 107 137 124

± ± ± ±

10 4 3 3

aspect ratios

Ix (kDa nm2) 108 93 116 113

± ± ± ±

8 6 4 2

Iz/Iy 1.3 1.3 1.4 1.7

± ± ± ±

0.2 0.1 0.1 0.1

asphericity

Iz/Ix

δ

± ± ± ±

0.010 0.016 0.022 0.041

1.5 1.4 1.6 1.9

0.2 0.2 0.1 0.1

δb 0.01 0.03

a

The values were averaged from the last 8 ns of MD simulation. The associated uncertainties are the standard deviations over the same 8 ns interval. The asphericity values are expressed with two significant figures after the decimal point to show a clear difference. bValues published by Lee et al.56

Figure 3. Two views of each dendrimer structure at the end of the simulations (t = 15 ns) for amine-terminated and acetylated PAMAM at low and neutral pH. The view at left is taken from the direction at which the dendrimer has the greatest apparent area; the view at right is taken from a direction rotated by 90° from the former. Ions and water molecules are not shown.

■

RESULTS AND DISCUSSION To characterize the models of amine- and acetyl-terminated PAMAM dendrimers at the molecular level, we performed a series of MD simulations that were subsequently subjected to structural and energetic analyses. First, we present the structural properties of the dendrimers in the absence of the drug at low and neutral pH. Next, we describe the construction of cavities for encapsulating drug within the dendrimers and subsequent docking of the drug into the cavities. Finally, molecular dynamics simulations of the dendrimer−drug complexes were performed; the resulting trajectories were analyzed via MMGBSA, and the specific dendrimer−drug interactions were characterized. Structural Characterization of the Dendrimers. Radius of Gyration. The radii of gyration (Rg) were calculated as a function of simulated time, as shown in Figure 2A. We furthermore calculated the mean radius of gyration (Rg), including only the last 8 ns of the 15 ns equilbration simulations, where Rg appeared to take on a steady value. For amine-terminated PAMAM, we obtained the values 25.2 ± 0.2 Å and 23.7 ± 0.2 Å at low and neutral pH, respectively, while, for acetyl-terminated PAMAM, we obtained 24.4 ± 0.3 Å at low pH and 25.3 ± 0.2 Å at neutral pH. Previous experiments56 at

neutral pH suggested a slightly smaller radius of gyration for 90% acetylated G5 PAMAM (23.5 Å) than we obtained here (25.3 ± 0.2 Å) for 100% acetylated PAMAM. This discrepancy might be explained by the difference in the number of acetyl groups or imply that the dendrimer structure is not yet fully equilibrated in our simulations. Solvent-Accessible Surface Area (SASA). Although we did not observe large differences in the radius of gyration between the systems ( amine-terminated PAMAM at low pH > acetylterminated PAMAM at neutral pH. This ranking was in agreement with the more accurate MM-GBSA method, which is described in the next section, as well as with the experimental results of Yang et al.,22 which showed a lower affinity of Dp21 for acetylated PAMAM at neutral pH than in the other three 6807



Article

of favorable ionic interactions. In the other three systems, ionic interactions were present because the phosphate group of Dp21 carries a negative charge at both pH values considered here, while amine-terminated PAMAM has at least some positively charged amines at both low and neutral pH, and the acetylterminated PAMAM has positively charged tertiary amines at low pH. Although the electrostatics made a significant contribution to the affinity of the dendrimer−drug complexes, the magnitude of van der Waals contribution (ΔHvdW MM ) appeared to be greater. The amine-terminated PAMAM at neutral pH and acetylterminated PAMAM at low and neutral pH show significantly vdW values than the amine-terminated more negative ΔHMM PAMAM at low pH. However, the differences in values of ΔGbind between the four systems depended on a balance between the van der Waals and electrostatic contributions: neither contribution alone was sufficient to rank the affinities among the systems. On the other hand, the contribution of solvation (ΔGsolv), which is the sum of a polar component, calculated by the Generalized Born implicit solvent method, and a nonpolar component, calculated from the solvent accessible surface area, did not show dramatic differences among the systems and did not affect the final ranking of affinities of the complexes. In the MM-GBSA method, the change in the conformational entropy is usually calculated by normal-mode analysis, which is often the most computationally expensive portion of the method. However, in many cases, the inclusion of the conformational entropy term does not improve predicted differences in ΔGbind between similar complexes.40 For this reason, the conformational entropy is often neglected when only comparisons of ΔGbind between similar complexes are desired (see, for example, ref 39). In this work, we were most interested in the differences in ΔGbind among the systems considered, and, therefore, did not calculate the conformational entropy due to its computational expense. Note that without the inclusion of the conformational entropy, the MM-GBSA method systematically gives unrealistically low values for the absolute ΔGbind, overestimating the stability of complexes40 relative to the uncomplexed form. Thus, the absolute ΔGbind values shown here are probably unrealistically low, and the affinities of Dp21 for the dendrimers are likely uniformly weaker than implied by these values. Characterization of Complexes. To better understand the origin of the binding free energies, we performed an atomicscale analysis of interactions between the dendrimers and Dp21. In particular, we identified all hydrogen bonds formed between the dendrimers and Dp21 as well as contacts between charged amine groups of the dendrimers and charged phosphate group of the drug. Hydrogen bonding analyses were carried out for each complex using a cutoff of 3.5 Å for the donor−acceptor distance and a cutoff of 120° for the donor−hydrogen− acceptor angle. The mean numbers of hydrogen bonds between the dendrimers and Dp21 for each simulation are shown in Figure 4. Most notably, we found comparatively few hydrogen bonds for acetyl-terminated PAMAM at neutral pH. Furthermore, more hydrogen bonds were identified when the Dp21 was in the middle region than the others. In these respects, larger numbers of hydrogen bonds were correlated with low values of ΔGbind calculated by MM-GBSA; however, while acetyl-terminated PAMAM at low pH showed the lowest

systems. Furthermore, AutoDock values suggested a weaker affinity for Dp21 at the surface than in the middle region of the dendrimers, which is also in agreement with the MM-GBSA calculations. Free Energy of Binding Estimated by MM-GBSA. We sought to obtain reliable estimates of the free energy of binding for each dendrimer−drug complex. There exist many computational methods of varying efficiency and accuracy to determine free energies of binding;36 however, the most robust methods are computationally expensive, especially for large systems such as the PAMAM-G5 dendrimers. MM-GBSA provides moderately accurate results at comparatively low computational cost.40 Much of the advantage of MM-GBSA is that it is an end-point methodonly the free and complexed states need to be considered, unlike more robust methods which usually require consideration of variously defined intermediate states. Furthermore, MM-GBSA has the advantage of permitting easy separation of various contributions to the free energy, such as the van der Waals, electrostatic, and solvation contributions. Table 3 shows details of the MM-GBSA calculations, indicating the contribution of each term to the binding free energy (ΔGbind) of the dendrimer−drug complexes. The ΔGbind values show that in all cases (both dendrimers at both pH values), the affinity of Dp21 is higher in the middle region of the dendrimers than at the surface or near the core. Moreover, in agreement with the experiments of Yang et al.,22 we find that the weakest affinity is in the case of acetyl-terminated PAMAM at neutral pH (−19.0 ± 0.4 kcal/mol). Interestingly, MM-GBSA calculation predicts the highest affinity in the case of acetylated PAMAM at low pH, where ionic interactions might appear to be weaker than for amineterminated PAMAM. At large separations between the Dp21 and dendrimer, the interaction should be describable by the Debye−Hückel model for two point charges Velec =

qdrug qden 4πε0εrR drug,den

exp( −κR drug,den)

(4)

where qdrug and qden are the respective charges of the drug and dendrimer, ε0 is the electric constant in vacuum, εr is the dielectric constant of the solvent, Rdrug,den is the distance between the drug and dendrimer, and κ is the inverse Debye length. The latter may be expressed as κ=

2NAe 2I ε0εrkBT

(5)

where I is the ionic strength. This model would predict the strongest affinity in the systems with the most negative qdrugqden. Thus, we would predict that amine-terminated PAMAM at low or neutral pH would have the highest affinity for Dp21, with qdrugqden = −254e and −256e, respectively, while we would expect the weakest affinity for acetyl-terminated PAMAM at neutral pH where qden = 0 (see Table 1). However, the interaction is much more complicated when the Dp21 is embedded in the dendrimer. Indeed, the highest affinity is found for acetylated PAMAM at low pH, where qdrugqden = −126e. The Δelec MM values (see Table 3) show favorable electrostatic contributions when the Dp21 is in the middle region for the cationic dendrimers, but a large positive value for the neutral dendrimer. The relatively low affinity of acetyl-terminated PAMAM for Dp21 at neutral pH can be attributed to the lack 6808



Article

Figure 6. Intermolecular hydrogen bonds for Dp21 interacting with amine-terminated PAMAM at (A) low and (B) neutral pH. The snapshots show the most common hydrogen bonds formed between the dendrons and the drug in the core, middle region, and surface of each dendrimer. Below each snapshot is shown the number of hydrogen bonds as functions of time.

ΔGbind, amine-terminated PAMAM at neutral pH had significantly more hydrogen bonds than the other systems. Figure 4 shows the percentage of hydrogen bonds between the dendrimers and Dp21 involving specific kinds of atoms. With the exception of acetyl-terminated PAMAM at neutral pH, Dp21 participated in hydrogen bonding mainly through its phosphate oxygens. At low pH, some of these oxygens could serve as both donors and acceptors, while they could only serve as acceptors at neutral pH, where they were unprotonated. In

contrast, few hydrogen bonds involving the phosphate oxygen atoms were observed for acetyl-terminated PAMAM at neutral pH. Instead, a majority of the contributions came from the two hydroxyl oxygen atoms (not associated with the phosphate group), the two carbonyl oxygen atoms, and the fluorine atom. Some of these hydrogen bonds had significantly larger absolute prevalences in acetyl-terminated PAMAM at neutral pH than for the other systems. For example, the mean number of hydrogen bonds involving fluorine was 42% for acetyl6809



Article

Figure 7. Intermolecular hydrogen bonds for Dp21 interacting with acetyl-terminated PAMAM at (A) low and (B) neutral pH. The snapshots show the most common hydrogen bonds formed between the dendrons and the drug in the core, middle region and surface of each dendrimer. Below each snapshot is shown the number of hydrogen bonds as functions of time.

terminated PAMAM at neutral pH in the core region, while it was 0% for the other systems in the same region. Exemplary images of hydrogen bonding between the components of the complexes are shown in Figure 6 for amine-terminated PAMAM and in Figure 7 for acetylterminated PAMAM. The chemical structures of the Dp21 and dendrimers depended on pH and acetylation; therefore, possible hydrogen bonds also depended on these factors. Acetylation removes the primary amines of PAMAM; thus, the

acetyl-terminated PAMAM trivially showed no hydrogen bonds involving primary amines. At low pH, the tertiary amines were important hydrogen bond donors, while at neutral pH, where they are unprotonated and could no longer serve as donors, they made negligible contributions to hydrogen bonding. The primary amines made negligible contributions at low pH, except in the surface region, owing to the fact that the primary amines were distributed mostly near the surface at this pH. At neutral pH, the primary amines were more uniformly 6810



Article

Table 4. Hydrogen Bonds between the Dendrimers and Dp21 in MD Simulationsa % PAMAM atom complexes PAMAM-G5/Dp21 (LowpH)

PAMAM-G5/Dp21 (Neutral pH)

PAMAM-G5-Ace/Dp21 (LowpH)

PAMAM-G5-Ace/Dp21 (Neutral pH)


mean h-bond 1.28 2.90 1.11 3.66 5.67 3.30 1.72 3.51 2.62 0.57 0.38 0.24

NAP

(±0.04) (±0.09) (±0.17) (±0.22) (±0.06) (±0.30) (±0.09) (±0.13) (±0.16) (±0.08) (±0.03) (±0.04)

5 1 36 87 48 70 0 0 0 0 0 0

NAT

NAm

75 67 19 0 0 0 73 49 44 3 30 27

%Dp2l atom OAm

9 17 36 9 51 27 17 32 38 89 17 48

11 14 9 4 1 3 9 19 18 8 53 25

OP

OH

85 87 85 79 85 85 92 88 66 1 2 0

9 9 7 12 11 10 5 8 26 9 65 36

OC 6 4 7 9 2 5 3 3 3 47 29 50

F 0 0 1 0 2 0 0 1 5 42 4 14

a Mean values were calculated after discarding the first 10 ns of the simulations. Statistical errors in the means (given in parentheses) were calculated according to method presented by Flyvbjerg.59 The percentages of the hydrogen bonds involving different kinds of atoms are also given. For PAMAM, we distinguish hydrogen bonds involving the nitrogen atoms of the primary amines (NAP), tertiary amines (NAT), and amide groups (NAm), as well as the oxygen of the amide groups (OAm). For Dp21, we distinguish hydrogen bonds involving one of the four oxygens of the phosphate group (OP), one of the two hydroxyl groups not attached to the phophorous atom (OH), one of the two carbonyl oxygen atoms (OC), or the fluorine atom (F).

distributed and, therefore, made significant contributions to the hydrogen bonding in all regions. The largest number of hydrogen bonds was observed for Dp21 in the middle region of amine-terminated PAMAM at neutral pH. In this case, the hydrogen bonds are split almost evenly between the protonated primary amines and the nitrogen of the amide groups. Clearly, interactions between the negatively charged phosphate group of Dp21 and positively charged amines of the dendrimers made an important contribution to the free energy of binding. Figure 5 shows the mean number of nitrogen atoms in a charged amine group within 4.5 Å of the phosphorus atom of Dp21, which are referred to here as P−N contacts. The change of pH leads to changes in the charges of the dendrimers and the drug (Table 1). Although the number of charged amine groups of amine-terminated PAMAM was almost a factor of 2 smaller at neutral pH than at low pH, a larger number of P−N contacts were observed at neutral pH, which was likely due to the greater charge of the Dp21 phosphate and greater flexibility of the dendrimer at neutral pH. Moreover, the system with the greatest affinity as suggested by MM-GBSA, acetyl-terminated PAMAM at low pH, did not show a greater number of P−N contacts than the other systems. Our results show that ionic interactions and hydrogen bonding, while somewhat correlated with the stability of the complexes, cannot be completely responsible for the calculated ΔGbind values. Other factors more difficult to characterize, such as hydrophobic interactions, could be important for the calculated ranking of affinities.

Table 5. Mean Numbers of Contacts between the Phosphate of Dp21 and Charged Amine Groups of the Dendrimersa complexes PAMAM-G5/Dp2l (Low pH)

PAMAM-G5/Dp2l (Neutral pH)

PAMAM-G5-Ace/Dp2l (Low pH)

PAMAM-G5-Ace/Dp2l (Neutral pH)


P−N contacts 0.96 1.92 0.56 2.30 2.15 1.76 1.46 1.68 0.93

(±0.03) (±0.04) (±0.09) (±0.15) (±0.05) (±0.12) (±0.09) (±0.10) (±0.05) 0 0 0

a

Because acetyl-terminated PAMAM at neutral pH has no charged amine groups, the value of zero is shown for the corresponding simulations. Statistical errors were calculated according to method presented by Flyvbjerg.59

to estimate the free energy of binding. An important aspect of our protocol was the positioning and confinement of the drug within a certain range of distances from the center of mass of the dendrimers. We found that the free energy of binding determined by MM-GBSA as well as structural properties, such as the prevalence of hydrogen bonds, varied significantly with the range of distance chosen. The most favorable free energy of binding was found at intermediate distances (15−20 Å between the centers of mass of the drug and dendrimer). Several measures, including ΔGbind calculated by MM-GBSA, the number of hydrogen bonds, and the number of charged amine−phosphate contacts, suggested particularly weak affinity between Dp21 and the acetyl-terminated dendrimer at neutral pH, consistent with experimental results.22 We found the dendrimer−drug affinity to be a complex balance of several factors and could not be characterized only by the generalized Born solvation energy, electrostatic and van der Waals energy, hydrogen bonds, or contact between charged groups. We anticipate that the protocol presented here will prove useful as an efficient means to estimate affinity between drugs

■

CONCLUSIONS In the present work, we characterized the structure of amineand acetyl-terminated PAMAM dendrimers at low and neutral pH, showing decidedly more open structures at low pH which allowed greater contact with the solvent. We observed a particularly compact disc-like shape for acetyl-terminated PAMAM at neutral pH. The observed structural differences were not apparent from the Rg values. Furthrmore, we demonstrated a novel and efficient protocol for constructing dendrimer−drug complexes. The MM-GBSA method was used 6811



Article

(15) Goldberg, D.; Ghandehari, H.; Swaan, P. Pharm. Res. 2010, 27, 1547−1557. (16) Esfand, R.; Tomalia, D. Drug Discovery Today 2001, 6, 427−436. (17) Menjoge, A.; Kannan, R.; Tomalia, D. Drug Discovery Today 2010, 15, 171−185. (18) Lim, J.; Lo, S.; Hill, S.; Pavan, G.; Sun, X.; Simanek, E. Mol. Pharmaceutics 2012, 9, 404−412. (19) Waite, C.; Sparks, S.; Uhrich, K.; Roth, C. BMC Biotech. 2009, 9, 38−47. (20) Fant, K.; Esbjorner, E.; Jenkins, A.; Grossel, M.; Lincoln, P.; Nordén, B. Mol. Pharmaceutics 2010, 7, 1734−1746. (21) Kolhatkar, R.; Kitchens, K.; Swaan, P.; Ghandehari, H. Bioconjugate Chem. 2007, 18, 2054−2060. (22) Yang, K.; Weng, L.; Cheng, Y.; Zhang, H.; Zhang, J.; Wu, Q.; Xu, T. J. Phys. Chem. B 2011, 115, 2185−2195. (23) Majoros, I.; Keszler, B.; Woehler, S.; Bull, T.; Baker, J., Jr. Macromolecules 2003, 36, 5526−5529. (24) Stasko, N.; Johnson, C.; Schoenfisch, M.; Johnson, T.; Holmuhamedov, E. Biomacromolecules 2007, 8, 3853−3859. (25) Lee, H.; Larson, R. J. Phys. Chem. B 2006, 110, 18204−18211. (26) Tanis, I.; Karatasos, K. J. Phys. Chem. B 2009, 113, 10984. (27) Cheng, Y.; Wu, Q.; Li, Y.; Xu, T. J. Phys. Chem. B 2008, 112, 8884−8890. (28) Hu, J.; Fang, M.; Cheng, Y.; Zhang, J.; Wu, Q.; Xu, T. J. Phys. Chem. B 2010, 114, 7148−7157. (29) Maingi, V.; Kumar, M. V. S.; Maiti, P. K. J. Phys. Chem. B 2012, 116, 4370−4376. (30) Jorgensen, W.; Severance, D. J. Am. Chem. Soc. 1990, 112, 4768−4774. (31) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992, 13, 1011−1021. (32) Roux, B. Comput. Phys. Commun. 1995, 91, 275−282. (33) Darve, E.; Rodrguez-Gómez, D.; Pohorille, A. J. Chem. Phys. 2008, 128, 144120. (34) Hénin, J.; Fiorin, G.; Chipot, C.; Klein, M. J. Chem. Theory Comput. 2009, 6, 35−47. (35) Bussi, G.; Laio, A.; Parrinello, M. Phys. Rev. Lett. 2006, 96, 90601. (36) Chipot, C.; Pohorille, A. Free Energy Calculations; Springer: Berlin, 2007. (37) Homeyer, N.; Gohlke, H. Mol. Informatics 2012, 31, 114−122. (38) Onufriev, A.; Bashford, D.; Case, D. A. Proteins: Struct., Funct., Bioinf. 2004, 55, 383−394. (39) Abroshan, H.; Akbarzadeh, H.; Parsafar, G. J. Phys. Org. Chem. 2010, 23, 866−877. (40) Hou, T.; Wang, J.; Li, Y.; Wang, W. J. Chem. Inf. Model. 2011, 51, 69. (41) Abagyan, R.; Totrov, M.; Kuznetsov, D. J. Comput. Chem. 1994, 15, 488−506. (42) MacKerell, A. D., Jr.; et al. J. Phys. Chem. B 1998, 102, 3586− 3616. (43) Saam, J.; Ivanov, I.; Walther, M.; Holzhütter, H.; Kuhn, H. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 13319. (44) Becke, A. Phys. Rev. A 1988, 38, 3098. (45) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10268−10280. (46) Frisch, M. J. et al. Gaussian Inc.: Wallingford CT, 2009. (47) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926−935. (48) Feller, S.; Zhang, Y.; Pastor, R.; Brooks, B. J. Chem. Phys. 1995, 103, 4613−4621. (49) Essmann, U.; Perera, L.; Berkowitz, M.; Darden, T.; Lee, H.; Pedersen, L. J. Chem. Phys. 1995, 103, 8577. (50) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781−1802. (51) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33−38.

and dendrimers, which may be essential for rational design of dendrimer-based nanocarriers. In applications, it is often desirable for the dendrimer to carry many molecules of the drug; thus, in future work, we will apply our protocol to systems containing many copies of the drug molecules. Further work also needs to be done to evaluate the accuracy of MMGBSA in dendrimer systems, which may be particularly challenging for the method due to the flexibility of dendrimers. A useful strategy would be to make comparisons of MM-GBSA calculations and calculations performed by more accurate methods.36

■

ASSOCIATED CONTENT

S Supporting Information *

Conformational stability of dendrimer−drug complexes and boundary forces in the MD trajectories used for MM-GBSA. The topology and parameter files for amine-terminated and acetylated PAMAM dendrimers. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.V.-J. thanks the Doctoral Program of Applied Sciences at the Universidad de Talca, as well as CONICYT-Chile for a doctoral fellowship. The authors acknowledge support from Proyecto ́ ACT1107 and the Centro Interdisciplinario de Anillo Cientifico ́ a Millennium Institute supported Neurociencia de Valparaiso, by the Millennium Scientific Initiative of the Ministerio de Economia,́ Fomento y Turismo.

■

REFERENCES

(1) Zhao, L.; Xu, F.; Chen, H.; Tang, G.; Hu, X. J. Appl. Polym. Sci. 2011, 123, 1509−1517. (2) Discher, D.; Ortiz, V.; Srinivas, G.; Klein, M.; Kim, Y.; Christian, D.; Cai, S.; Photos, P.; Ahmed, F. Prog. Polym. Sci. 2007, 32, 838−857. (3) Loverde, S.; Pantano, D.; Christian, D.; Mahmud, A.; Klein, M.; Discher, D. Curr. Opin. Solid State Mater. Sci. 2011, 15, 277−284. (4) Pang, Z.; Gao, H.; Yu, Y.; Chen, J.; Guo, L.; Ren, J.; Wen, Z.; Su, J.; Jiang, X. Int. J. Pharm. 2011, 415, 284−292. (5) Ghalandari, B.; Monajjemi, M.; Mollaamin, F. J. Comput. Theor. Nanosci. 2011, 8, 1212−1219. (6) Hadad, A.; Azevedo, D.; Caetano, E.; Freire, V.; Mendonca, G.; Neto, P.; Albuquerque, E.; Margis, R.; Gottfried, C. J. Phys. Chem. C 2011, 115, 24501−24511. (7) Bromley, S.; Flikkema, E. Comput. Mater. Sci. 2006, 35, 382−386. (8) Costache, A.; Sheihet, L.; Zaveri, K.; Knight, D.; Kohn, J. Mol. Pharmaceutics 2009, 6, 1620. (9) Chen, A. Z.; Chen, M. Y.; Wang, S. B.; Huang, X. N.; Liu, Y. G.; Chen, Z. X. J. Appl. Polym. Sci. 2012, 124, 3728−3736. (10) Naylor, A.; Goddard, W., III; Kiefer, G.; Tomalia, D. J. Am. Chem. Soc. 1989, 111, 2339−2341. (11) Tomalia, D.; Naylor, A.; Goddard, W., III Angew. Chem., Int. Ed. Engl. 1990, 29, 138−175. (12) Mintzer, M.; Dane, E.; O’Toole, G.; Grinstaff, M. Mol. Pharmaceutics 2011, 9, 342−354. (13) Giovino, C.; Ayensu, I.; Tetteh, J.; Boateng, J. Int. J. Pharm. 2012, 428, 143−151. (14) Manchun, S.; Dass, C.; Sriamornsak, P. Life Sci. 2012, 90, 381− 387. 6812



Article

(52) Morris, G.; Goodsell, D.; Halliday, R.; Huey, R.; Hart, W.; Belew, R.; Olson, A. J. Comput. Chem. 1998, 19, 1639−1662. (53) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; MacKerell, A. D. J. Comput. Chem. 2010, 31, 671−690. (54) Morris, G.; Huey, R.; Lindstrom, W.; Sanner, M.; Belew, R.; Goodsell, D.; Olson, A. J. Comput. Chem. 2009, 30, 2785−2791. (55) Tanner, D.; Chan, K.; Phillips, J.; Schulten, K. J. Chem. Theory Comput. 2011, 7, 3635−3642. (56) Lee, H.; Baker, J., Jr.; Larson, R. J. Phys. Chem. B 2006, 110, 4014−4019. (57) Maiti, P.; Cagin, T.; Wang, G.; Goddard, W., III Macromolecules 2004, 37, 6236−6254. (58) Ballauff, M.; Likos, C. N. Angew. Chem., Int. Ed. 2004, 43, 2998− 3020. (59) Flyvbjerg, H. Adv. Comput. Simul. 1998, 88−103.

6813


Computationally Efficient Methodology for Atomic ... - ACS Publications

Recommend Documents