Steered Molecular Dynamics Studies of the Potential of Mean Force

Mar 19, 2013 - The Journal of Physical Chemistry A · Advanced Search ..... Tell Tuttle .... There and back again: The tale of 2 asteroid sample-return...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCA

Steered Molecular Dynamics Studies of the Potential of Mean Force for Peptide Amphiphile Self-Assembly into Cylindrical Nanofibers Tao Yu, One-Sun Lee, and George C. Schatz* Department of Chemistry, Northwestern University, 2145 Sheridan Road, Evanston, Illinois 60208-3113, United States ABSTRACT: Steered molecular dynamics (SMD) simulations were applied to determine the potential of mean force for the self-assembly of peptide amphiphile (PA) nanofibers, specifically considering a single PA adding to a growing cylindrical nanofiber at 310 K. It is found that the free energy, enthalpy, and entropy differences for this assembly process are −67 kcal/mol, −71.5 kcal/ml, and −14.5 cal/(mol K), respectively, and therefore that enthalpy provides the driving force for self-assembly to form a fiber. A pairwise interaction analysis shows that both electrostatic and van der Waals interactions play important roles in the self-assembly process, with the van der Waals interaction being the larger effect. The mechanistic picture that emerges from this work is that as the PA is pulled from the fiber, the interaction evolves through three stages: (1) initially electrostatic interactions between the charged head of the pulled PA and other PAs, and between the pulled PA and solvent are dominant, (2) after the charged head emerges, the rest of the peptide comes out, with both PA-solvent electrostatic interactions and van der Waals interactions being significant, and (3) in the last step, the alkane tail emerges, dominated by van der Waals interactions with either peptide or solvent.

■

INTRODUCTION The self-assembly of peptide amphiphiles (PAs) plays a significant role to design functional supramolecular structures.1−6 A very interesting structure, cylindrical nanofibers are formed by the self-assembly of PAs in aqueous solution as reported by Stupp and co-workers,4,5,7 and others, as has been reviewed.6,8−11 By appropriately choosing the residues in the peptide, these nanoscale fibers can be programmed for biological activity, and as a result, they have been used for many medicinal purposes, such as bone healing, wound repair, and nerve growth. The PAs consist of two segments, a hydrophobic tail segment and a peptide head segment, with the latter containing ionizable residues that lead to aggregation as pH is lowered. When the nanofiber structures form, the hydrophobic tail segments are packed together in the center of the fiber, while portions of the peptide head segments form beta-sheet secondary structures, and other portions allow electrostatic interactions that stabilize the cylindrical structures. The use of theory to predict the self-assembly of PA nanofibers has been of interest for a long time,12−19 as the choice and number of residues in the peptide that leads both to cylindrical micelle assembly and to biological activity has been hard to predict without modeling. There has also been uncertainty as to what interactions (hydrophobic effect, hydrogen bonding leading to beta sheets, and electrostatic interactions) are more important in driving self-assembly. Unfortunately, it has been difficult to develop atomistic theories capable of assembling a complete micelle.16 Recently, a seeded molecular dynamics (MD) simulation was developed by Lee et al.,20−22 which enables atomistic structural detail in the description of cylindrical nanofibers based on the peptide © 2013 American Chemical Society

sequence: SLSLAAAEIKVAV. They also investigated the influence of amino acid sequence on the secondary structures of the PAs in the nanofiber structures. Furthermore, a coarse grained model was developed which can simulate the selfassembly process of this and other PAs into cylindrical micelle fibers starting from a homogeneous mixture of PAs in aqueous solution. The results indicate that spherical micelle structures are initially formed, and then these micelles combine to give cylindrical micelles, followed by further aggregation to form fibers. None of the research so far has studied the driving forces, interaction energies, and thermodynamics of the PA selfassembly process. However, this type of research is significant to understand the mechanism of self-assembly, including the further design of experiments to control the topology of the self-assembled PAs, and how this can be modified by variations in temperature, solvent, templating, and nonequilibrium effects. In the present work, we study the potential of mean force (PMF) or free energy profile for the self-assembly of PAs into cylindrical micelles using steered molecular dynamics simulations (SMD)23−25 based on a simple model that considers the aggregation of a single PA and fiber composed of the same PAs to give a fiber with one additional PA. The peptide sequence is chosen to be the same as in the earlier work of Lee et al.,20 so that we can build on the earlier results. On the basis of the Special Issue: Joel M. Bowman Festschrift Received: February 11, 2013 Revised: March 16, 2013 Published: March 19, 2013 7453

dx.doi.org/10.1021/jp401508w | J. Phys. Chem. A 2013, 117, 7453−7460

The Journal of Physical Chemistry A

Article

earlier coarse-grained modeling,20 addition of a single PA to a micelle is likely to be an important process in the earliest stages of self-assembly where small micelles grow by adding additional PA’s. In addition, this model is important for understanding micelle stability, as dissociation of a PA from a micelle is a possible mechanism for micelle destruction. Micelle destruction is a known process that arises when the pH is raised such that charged residues are neutralized. SMD is an important method for obtaining the PMF or free energy profile for various biological processes, and has been widely used for investigating protein folding,26−28 protein− protein interactions,29,30 and self-assembly issues.31−33 The SMD simulation is based on coupling external mechanical forces to a biomolecular system and applying Jarzynski’s equality34,35 to calculate the PMF (an equilibrium property) using a nonequilibrium simulation. In this work, we set up the simulation by pulling out a single PA from the preformed micelle nanofiber structure. We investigate free energy changes based on the Jarzynski equality formalism and using a cumulant expansion approximation up to second order. The PMF demonstrates that, although the self-assembly of PAs is a physical process, with no covalent bond breaking or forming, it still costs a large free energy, around 67 kcal/mol, to dissociate a single PA from a stable fiber structure. On the basis of a further thermodynamic analysis, the self-assembly process is found to be primarily driven by enthalpic effects, which dominate over entropic changes that destabilize the micelle. The pairwise interactions between the pulled PA and the remaining PAs in the fiber or between the pulled PA and the surrounding solvent, including water, Na+, and Cl− ions, are also analyzed along the reaction coordinate. The results demonstrate that the pulled PA in the fiber (bound state) and out of the fiber (free state) display different interactions with its environment. We also show that (1) both electrostatic and van der Waals interactions contribute to the driving forces to form the stabilized nanofiber structure, with the van der Waals interaction energy being larger than the electrostatic one; (2) the electrostatic interaction of the pulled PA with its environment in the bound state is dominated by the head part of the PA, and the van der Waals interaction is attributed to both head and tail segments of the PA, with the tail’s contribution larger than that from the head. The overall PMF is then composed of segments that reflect changes in these interactions as the PA is pulled out. This provides a semiquantitative analysis of a process that has not previously been available.

Figure 1. (a) Schematic representation of PA self-assembly into cylindrical micelles. (b) PA used in MD simulations, where the PA sequence is SLSLAAAEIKVAV. (c) Initial structure of the PAs used in the MD simulations, including views from the top (left) and side (right) of the micelle. Nine PAs are placed radially in the first layer (purple) where their tails are pointing inward. The second layer also has 9 PAs (green) and is rotated by 20° relative to the first layer. The angle between neighboring PAs is 40°, and the distance between layers is 5 Å.

structure. The procedure leading to these choices of seed parameters was described by Lee et al, where it was demonstrated that this leads to the formation of a stable micelle structure using MD calculations. The 90 PA molecules in the seed structure were solvated in a water box using the SOLVATE21 module36 implemented in VMD.37 Periodic boundary conditions were used in the simulations, using a box of dimensions of 154 Å, 154 Å, and 48 Å in the x, y, and z directions, respectively. This box was filled with 25729 water molecules based on the modified TIP3P potential.38,39 Note (Figure 1b) that there are three charged residues close to the C-terminal in the head segment of each PA: Glu (residue 8), Lys (residue 10), and Val (residue 13, as C-terminal). The net charge of the PA is minus one. To neutralize the system, 90 Na+ ions are added. In addition to these sodium ions, another 73 Na+ and Cl− ions are added to make the concentration of Na+ ion be 0.15 M (to match experimental conditions). A 1000 step energy minimization was applied to the solvated system to remove the high-energy contacts. MD simulations were carried out using NAMD2.9.40 A 2 ns MD simulation at 310 K with a NPT ensemble was performed to equilibrate the system. In the production period, the system was simulated for 40 ns using the NPT ensemble and Langevin dynamics41 at a temperature of 310 K with a damping coefficient γ = 5 ps−1. Pressure was maintained at 1 atm using the Langevin piston method41,42 with a piston period of 100 fs, a damping time constant of 50 fs, and a piston temperature of 310 K. No atomic coordinates were restrained during the production period. Full electrostatics was employed using the particle-mesh Ewald method with a 1 Å grid width.43 Nonbonded interactions were calculated using a group-based cutoff with a switching function and updated every 10 time steps. Covalent bonds involving hydrogen were held rigid using the SHAKE algorithm,44 allowing a 2 fs time step. Atomic coordinates were saved every 10 ps for the trajectory analysis. 2. SMD Simulations. The equilibrated nanofiber structure obtained by the simulation in the previous section was resolvated in another water box using the SOLVATE2136 module implemented in VMD.37 Again periodic boundary conditions were used in the simulations, using a box of

■

COMPUTATIONAL DETAILS 1. Construct Nanofiber Structures. As noted above, the peptide sequence chosen for each PA is Ser-Leu-Ser-Leu-AlaAla-Ala-Glu-Ile-Lys-Val-Ala-Val (see Figure 1b for a detailed structure). To define the micelle, we used the same procedure as in Lee et al.20 in which circular disks of PAs are initially constructed to “seed” the micelle structure. To define one circular disk, nine PAs are radially placed on a plane (purple PAs in Figure 1c) with tails pointing inward. The starting structure of the backbone of each PA is assumed to be an extended conformation. The angle between adjacent PAs is 40°. The second disk (green PAs in Figure 1c) is taken to be identical to the first but rotated by 20° relative to the first and with the distance between layers taken to be 5 Å. A total of 10 disks that alternate between the first and second structures are placed along the fiber axis to define the complete starting 7454

dx.doi.org/10.1021/jp401508w | J. Phys. Chem. A 2013, 117, 7453−7460

The Journal of Physical Chemistry A

Article

each initial r value, five separate independent simulation runs were carried out, so the total number of trajectories for all pulling distances is 5 × 9 = 45. In each run, the pulling speed is 2.5 Å/ns. Such a slow speed was applied in order to reduce the nonequilibrium effect in the simulation. The force constant for the pulling process was 200 kcal/mol/Å. Such a large value of force constant is used to satisfy the stiff spring approximation and guarantee small deviation of the reaction coordinate from the constrained position. All these simulations were performed using NAMD version 2.9,40 and visualizations were made using VMD.31 3. PMF Calculation Based on Jarzynski’s Equality. The PMF is a potential as a function of a predefined reaction coordinate, the gradient of which represents the negative of the effective force on the molecule averaged over all the configurations at a given point parametrized by the reaction coordinate. In the SMD simulation, under the stiff spring approximation, the PMF is identical to the change of free energy along the reaction coordinate. In the SMD simulation, the external nonequilibrium work done by the pulling force can be obtained using the following:

dimensions of 178 Å, 124 Å, and 48 Å in the x, y, and z directions, respectively. In this module, the water box contained 25132 water molecules based on the modified TIP3P potential.32,33 To neutralize the system, 90 Na+ ions are added, and 71 Na+ and Cl− ions are added to make the concentration of Na+ ions be 0.15 M. A 1000 step energy minimization was applied to the solvated system to remove high-energy contacts. MD simulations were carried out using NAMD2.9.40 A 2 ns MD simulation at 310 K with a NPT ensemble was performed to equilibrate the system. Figure 2 shows a representative configuration of a micelle with highlighting used for the PA that we chose to pull (further

W = −kv

∫0

t

dt ′[x(t ′) − x0 − vt ′]

(1)

Figure 2. Representative configurations taken from MD simulations of a PA being pulled out from a cylindrical micelle nanofiber: (a) initial state before pulling force added to the PA where the reaction coordinate r is 29 Å; (b,c) intermediate states when the PA is partially out of the fiber; (d) final state where the PA is completely out of the fiber.

where W, k, and v are the work, force constant for pulling, and pulling velocity. x(t′) and x0 are the reaction coordinate at t′ in the simulation and the relative initial position of the center of mass of the pulled PA, respectively. The associated force f to accomplish this work was calculated using

details of the choice are given later). Using this PA to define the x-direction of a coordinate system, the SMD involved pulling the center of mass of this PA along the x-direction while the other 89 PAs were constrained by applying a harmonic force along the x-direction at the last carbon of the hydrophobic tail of each PA with the spring constant 20 kcal/mol/Å2. With this constraint, the center of mass of the rest of the fiber structure is fixed along the x-coordinate with only small fluctuations in the SMD simulation. External harmonic forces were employed for all the atoms in the pulled PA to constrain the respective center of mass of the PA. The reaction coordinate, denoted r, is defined as the relative distance of the center of mass of the pulled PA with respect to the center of mass of the rest of fiber. In the present work, the pulling range was found to vary from r = 29.2 Å when the PA is equilibrated with the micelle to r = 74.4 Å when fully dissociated. This range was divided into nine windows, and in each window the pulling distance is 5.5 Å, with 0.5 Å out of the 5.5 Å being used to overlap with the next window. The beginning position of each SMD pulling window therefore corresponds to r = 29.2, 33.3, 38.5, 43.2, 48.5, 54.0, 58.7, 63,5, and 68.9 Å . These values were obtained in the following way. First, we employed SMD to pull out the selected PA from r = 29.2 to r = 74.4 Å using a fast pulling speed, 2500 Å/ns, to generate a series of intermediate configurations along the reaction coordinate. Then we select the intermediate configurations at r = 29.2 (the initial bound state), 34, 39, 44, 49, 55, 59, 63, and 68 Å to run MD simulations by adding constraints to hold the r close to the selected values. After equilibrating 2 ns, we picked the last frame in each trajectory to define the beginning position for the SMD calculation. The r values associated with these frames had the values given above, 29.2, 33.3, 38.5, 43.2, 48.5, 54.0, 58.7, 63.5, and 68.9 Å. For

The equality between work and change in free energy was established by Jarzynski22,23 in terms of the formula

f = −k[x(t ′) − x0 − vt ′]

ΔG = −β −1 ln⟨exp( −βW )⟩

(2)

(3)

where β is equal to kT, and k is Boltzmann factor. ΔG is the free energy difference. In practice, the second order cumulant expansion of eq 3 is widely used to calculate the free energy difference using a limited number of SMD simulation trajectories. This converts eq 3 to the following: ΔG = ⟨W ⟩ −

β [⟨W 2⟩ − ⟨W ⟩2 ] 2

(4)

In the present work, we employed eq 4 to calculate the PMF profile for each window based on five SMD simulation trajectories. The final PMF profile along the whole reaction coordinate was reconstructed from the PMFs in the nine windows by using a matching condition: P′i + 1(x) = Pi + 1(x) ×

Pi(xa) Pi + 1(xb)

(5)

where xa ≈ xb and Pi+1(x) and P′i+1(x) are probability distributions in the (i + 1)th window before and after matching using the ith window. Pi(x) is proportional to exp(−βΔGi).

■

RESULTS AND DISCUSSION Figure 2 shows snapshots of the simulations as viewed from the top of the micelle, showing the single PA at different pulling stages as defined by the reaction coordinate r. Each configuration was obtained by equilibrating the system for 2 ns with a constraint added to the system to fix the center of mass of the PA and fix the last carbon atom in the tails of the 7455

dx.doi.org/10.1021/jp401508w | J. Phys. Chem. A 2013, 117, 7453−7460

The Journal of Physical Chemistry A

Article

remaining 89 PAs in the fiber along the x-direction only. Figure 2a displays the fully self-assembled state of the pulled PA, denoted the “bound” state. In the bound state, the pulled PA experiences a surrounding environment that depends strongly on its location in the micelle. The head segment of the pulled PA is surrounded by both solvent and the other eight PAs as nearest neighbors, and the tail segment of the pulled PA is packed in the hydrophobic core of the fiber and only surrounded by the tails belonging to the other PAs. On the other hand, Figure 2d displays the configuration where the pulled PA is completely out of the fiber, and the pulled PA is surrounded by solvent molecules only. In addition, the pulled PA has different conformations that are adapted into different environments at each pulling state. Figure 2 shows that the head part (the hydrophilic peptide) has mostly a random coil secondary structure in all pulling stages, while the hydrophobic tail part has different orientations in different configurations, being extended in Figure 2a and more globular in Figure 2d. It is found that during the pulling process, the rest of the PAs still form a stable fiber structure, and the external forces used in the SMD simulation are just a small local perturbation to the PAs. Thus the forces do not cause any damage to the fiber structure, and in Figure 2d the hole left by the leaving PA has completely healed over. The radii of the fiber in the bound and free states are both around 44 Å. We employed the stride module in VMD31 to assign the secondary structure for the fiber in the free and bound states. This module is based on an algorithm that combines the use of both hydrogen bond energies and statistically derived backbone torsional angle information. Figure 3 shows the secondary structure distribution for the bound and free state of the fiber. This demonstrates that the dissociation of a single PA, which has a random coil structure, does not cause any significant change in the secondary structure. In particular, the beta-sheet strand structure is stable enough and is well retained under the pulling condition. This property indicates that the perturbative forces applied in the SMD simulation introduce only local effects and do not interrupt the internal properties of the fiber on a larger scale. To further study this, a local structure change analysis between the bound and free states was performed. Figure 4a,b displays the configurations of the eight PAs which are the nearest neighbors of the pulled PA in the two different states. This shows that the conformations of both the tail and head segments of the nearest PA have significantly rearranged after pulling the PA out of the fiber, and all eight PAs are more packed in the free state; however, most of the PA’s have only moved by small amounts. For the PMF calculation, note that the self-assembled fiber structure is not an ideal cylinder, showing a very large inhomogeneous secondary structure distribution. Therefore, pulling the PAs with different secondary structures will give different PMF profiles. It has been demonstrated by both experimental and simulation studies that the dominating secondary structure in the cylinder fiber is random coil. For this reason, we chose the PA with most random coil structure for the SMD simulation. The PMF profile as a function of the reaction coordinate r is displayed in Figure 5. This shows that the PMF continuously increases as the pulled PA transits from bound to free state, although some bumpy features are seen, especially near r = 39 Å. The bumpy features might be related with the fact that the three charged residues (Glu, Lys, and Val as charged Cterminal) are completely pulled out of the fiber when r is

Figure 3. Distribution of secondary structure of the PAs in the fiber along the peptide sequence in (a) the bound state and (b) the free state.

Figure 4. Representative configuration taken from MD simulations showing the eight PAs’ (a) head (backbone atoms) and (b) tail (carbon atoms) segment, which are the nearest neighbors of the pulled PA in the bound state (green) and free state (blue).

around 39 Å, as shown in Figure 6. The total free energy difference between the bound and free state is around 67 kcal/ mol. This value demonstrates that the bound state is considerably more stable than the free state. Note that the range where r is between 30 and 40 Å corresponds to the regime where the PMF increases by about half of the total energy change from bound to free (about 35 kcal/mol) and it also increases the most rapidly. According to Figure 5, this is where the charged residues on the PA are being pulled out of the micelle and where the tail segment is pulled out of the hydrophobic core of the micelle. The region from 40 to 55 Å corresponds to emergence of the rest of the peptide, and we see that the free energy increase is nearly 30 kcal/mol for this portion of the PMF. Finally, the PMF becomes completely flat 7456

dx.doi.org/10.1021/jp401508w | J. Phys. Chem. A 2013, 117, 7453−7460

The Journal of Physical Chemistry A

Article

the average of the electrostatic and van der Waals interactions. Figure 8a shows the electrostatic interaction energy between

Figure 5. PMF profile along the reaction coordinate r. The PMF is calculated by using the second-order cumulant approximation of Jarzynski’s equality.

Figure 6. Representative configuration taken from MD simulations showing relative positions of charged resides, Glu, Lys, and Val in the pulled PA when r = 29 Å (orange) and r = 39 Å (red).

at 70 Å, which means that emergence of the alkane tail takes place between 55 Å and 70 Å, requiring 2−3 kcal/mol. Figure 7 displays the force spectrum as a function of r. This shows that at the very beginning of the pulling process, in the

Figure 8. The (a) electrostatic energy and (b) van der Waals energy between the pulled PA and the rest of the fiber (red), and between the pulled PA and the solvent (blue) as a function of r.

the pulled PA and other PAs in the fiber, and also the electrostatic interaction energy between the pulled PA and solvent as a function of r. Note that even in the bound state, the hydrophilic peptide segment of the pulled PA is still partially surrounded by solvent. Indeed, the figure shows that for r = 30 Å, the electrostatic interaction arises about equally from the fiber and solvent, labeled as E-PA-fiber and E-PA-solvent, respectively. Thus the total electrostatic interaction energy at that point is −655.3 kcal/mol in the bound state, of which the E-PA-fiber component is −342.8 kcal/mol and E-PA-solvent component is −312.5 kcal/mol, respectively. Figure 8a also shows that E-PA-fiber monotonically decreases as r increases, and then flattens near r = 55 Å where the peptide portion has been pulled out of the fiber. Note that E-PA-fiber profile displays a large jump between 30 and 39 Å, which corresponds to the range that the charged residues in the head of the PA have separated from the fiber. By contrast, although some fluctuations exist, the magnitude of E-PA-solvent increases more gradually with increasing r until it begins to flatten at about r = 50 Å and becomes completely flat at 70 Å. Figures 9a and 10a show the decomposition of the electrostatic energies in Figure 8a into peptide (head) and alkane (tail) interactions. These figures confirm that the electrostatic energies are dominated by peptide interactions with the other PA’s and with the solvent. Note in Figure 8a that in the free state, E-PAfiber is close to zero, and E-PA-solvent is −629.6 kcal/mol. Comparing the bound and free states, the overall electrostatic

Figure 7. Pulling force along the reaction coordinate r, defined as the relative distance between the center of mass of the pulled PA and the rest of the fiber. The force constant is 200 kcal/mol Å2, and the pulling speed is 2.5 Å/ns.

range where r is between 30 and 35 Å, the force used to pull out the PA is very large (∼600 pN). By contrast, near the end of the pulling process at r = 70 Å, the pulling force is fluctuating around zero. The interactions between the pulled PA and fiber can be attributed to electrostatic and van der Waals nonbonded interactions. We used five trajectories with r values around 29.2, 33.3, 38.5, 43.2, 48.5, 54.0, 58.7, 63.5, and 68.9 Å to calculate 7457

dx.doi.org/10.1021/jp401508w | J. Phys. Chem. A 2013, 117, 7453−7460

The Journal of Physical Chemistry A

Article

Figure 9. The (a) electrostatic energy and (b) van der Waals energy between the head of the PA (the peptide part of the PA) and the rest of the fiber (red), and between the head of the PA and the solvent (blue) as a function of r.

Figure 10. The (a) electrostatic energy and (b) van der Waals energy between the tail of the PA (the hydrocarbon part of the PA) and the rest of the fiber (red), and between the tail of the PA and the solvent (blue) as a function of r.

energy change for the fiber + PA → fiber assembly process is −25.4 kcal/mol. As should be apparent from Figure 2, the hydrophobic tail segment of the pulled PA experiences the largest environmental changes between the bound and free states. In the bound state, the tail is buried in the hydrophobic core of the fiber, where all the tails of the PAs are packed together and form a water exclusive nonpolar environment. In the free state, the tail is now exposed to water and surrounded by a polar environment. Since the tail atoms only have small charges on them, the electrostatic interaction of the tail with its environment in both the bound and free states and in states in-between is very small, and the electrostatic interaction for the pulled PA is dominated by the hydrophilic peptide segment of the pulled PA, as shown in Figure 9a and Figure 10a. Considering that the number of atoms in the peptide segment is 191, the electrostatic energy is −3.4 kcal/mol and −3.3 kcal/mol per atom in the bound and free states, respectively. Figure 8b shows the van der Waals interaction energy between the pulled PA and the fiber, and the van der Waals interaction energy between the pulled PA and solvent as a function of r. This shows that the van der Waals energy components are generally smaller than the corresponding electrostatic energy components. We denote the van der Waals interaction energy with fiber and solvent as VdW-PA-fiber and VdW-PA-solvent, respectively. In the bound state, the total van der Waals interaction energy is −125.4 kcal/mol, of which the VdW-PA-fiber component accounts for 88% (−110.0 kcal/ mol) and the VdW-PA-solvent component is 12% (−15.4 kcal/

mol). This result means that the VdW-PA-fiber is dominant for the total van der Waals interaction energy as the PA is pulled from the bound state. Meanwhile, Figure 8b shows that the magnitude of the VdW-PA-fiber interaction monotonically decreases to zero as r increases, while the magnitude of the VdW-PA-solvent increases. These results are further confirmed in Figures 9b and 10b, which show the decomposition of the van der Waals interaction into head and tail components. We see that the head component is generally larger than the tail component, and both components have similar dependence on r, with the fiber component decreasing with r and the solvent component increasing in magnitude. Overall, in the free state, the VdW-PA-fiber energy is zero, and the VdW-PA-solvent is −79.3 kcal/mol. Comparing the bound and free state, the van der Waals energy change (for the fiber + PA → fiber assembly process) is −46.1 kcal/mol. Unlike the electrostatic interactions, the hydrophobic tail segment of the pulled PA plays an important role in the van der Waals interaction for the pulled PA in both the bound and free state, as displayed in Figure 10b. In the bound state, all the tails of the PAs are packed together due to the hydrophobic effect. As a result, the van der Waals energy in this state is about −36.0 kcal/mol, about half the head-fiber interaction in magnitude. Since the tail contains 46 carbon and hydrogen atoms, the individual van der Waals interaction energy per atom is 0.8 kcal/mol. In the free state, the tail atoms touch with solvent directly, and the van der Waals energy is 24.1 kcal/mol for the whole tail and 0.5 kcal/mol per atom. The stronger magnitude 7458

dx.doi.org/10.1021/jp401508w | J. Phys. Chem. A 2013, 117, 7453−7460

The Journal of Physical Chemistry A

Article

of the tail van der Waals interaction with the fiber demonstrates that the tail segment prefers the hydrophobic environment. For the head segment of the pulled PA, Figure 9b shows that in the bound state, the van der Waals interaction energy is −89.4 kcal/mol (74.1 kcal/mol from the head interacting with fiber and 15.3 kcal/mol from interacting with the solvent) corresponding to −0.5 kcal/mol per atom. In the free state, the van der Waals interaction is reduced to −55.3 kcal/mol in total and 0.3 kcal/mol per atom. Recall that the difference in the electrostatic and van der Waals energies between the bound and free states is 25.4 and 46.1 kcal/mol, respectively. Since there is no significant conformational change for the fiber and the pulled PA between the free and bound states, the self-energy, including the bond, angle, and dihedral angle interactions should not change much, either. Therefore the dominant terms determining the binding enthalpy come from the electrostatic and van der Waals interactions between the fiber and pulled PA. Thus the total binding enthalpy for assembly can be estimated as −71.5 kcal/ mol. Combining this information with the free energy result (−67 kcal/mol for the assembly process) given earlier, we can extract the entropic contribution to PA assembly using the thermodynamic relation: ΔG = ΔH − T ΔS

comparable interaction energies, each being about half of the total free energy change, while energy associated with the third step is quite small, only a few kcal/mol. The ratio of the first two energies is close to unity in this application, but one could imagine that this could vary considerably for other choices of peptides, as determined by the relative importance of charged head groups and hydrogen-bonded interactions that lead to beta-sheet and other peptide structures. Also, the weak interaction associated with step 3 will play a crucial role in PA aggregation, basically determining whether aggregation can occur at all. Indeed other choices of alkane tail lengths, and of environmental parameters such as temperature and salt concentration are likely to influence this interaction strongly. The present calculations can thus provide an accurate approach for studying PA self-assembly that can be used to assess the importance of electrostatic and van der Waals interactions that govern self-assembly.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

(6)

ACKNOWLEDGMENTS The authors appreciate useful discussions with Dr. Xiaohu Li about the PMF calculation. This research was supported by the National Science Foundation (Grant CHE-1147335, support for T.Y.), and by the DOE NERC EFRC (DE-SC0000989, support for O.-S.L).

where ΔG, ΔH, and ΔS are the differences in free energy, enthalpy, and entropy between the bound and free states. The resulting ΔS is found to be −14.5 cal/(mol K). Also, the ratio between ΔH and TΔS is 15.9. As expected, the entropy is largest in the free state, which means than the transition from free state to bound states of the single PA is an enthalpy-driven process. Both electrostatic and van der Waals play important roles for the binding process, with the van der Waals interaction energy providing the larger contribution for the binding enthalpy.

■

REFERENCES

(1) Murakami, Y.; Nakano, A.; Fukuya, K. Stable Single-Compartment Vesicles with Zwitterionic Amphiphile Involving an Amino-Acid Residue. J. Am. Chem. Soc. 1980, 102, 4253−54. (2) Murakami, Y.; Nakano, A.; Yoshimatsu, A.; Uchitomi, K.; Matsuda, Y. Characterization of Molecular Aggregates of Peptide Amphiphiles and Kinetics of Dynamic Processes Performed by SingleWalled Vesicles. J. Am. Chem. Soc. 1984, 106, 3613−23. (3) Vauthey, S.; Santoso, S.; Gong, H. Y.; Watson, N.; Zhang, S. G. Molecular Self-Assembly of Surfactant-Like Peptides to Form Nanotubes and Nanovesicles. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 5355−60. (4) Hartgerink, J. D.; Beniash, E.; Stupp, S. I. Self-Assembly and Mineralization of Peptide-Amphiphile Nanofibers. Science 2001, 294, 1684−88. (5) Hartgerink, J. D.; Beniash, E.; Stupp, S. I. Peptide-Amphiphile Nanofibers: A Versatile Scaffold for the Preparation of Self-Assembling Materials. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 5133−38. (6) Shimizu, T.; Masuda, M.; Minamikawa, H. Supramolecular Nanotube Architectures Based on Amphiphilic Molecules. Chem. Rev. 2005, 105, 1401−43. (7) Behanna, H. A.; Donners, J. J. J. M.; Gordon, A. C.; Stupp, S. I. Coassembly of Amphiphiles with Opposite Peptide Polarities into Nanofibers. J. Am. Chem. Soc. 2005, 127, 1193−200. (8) Buerkle, L. E.; Rowan, S. J. Supramolecular Gels Formed from Multi-Component Low Molecular Weight Species. Chem. Soc. Rev. 2012, 41, 6089−102. (9) Woolfson, D. N.; Mahmoud, Z. N. More Than Just Bare Scaffolds: Towards Multi-Component and Decorated Fibrous Biomaterials. Chem. Soc. Rev. 2010, 39, 3464−79. (10) Palmer, L. C.; Stupp, S. I. Molecular Self-Assembly into OneDimensional Nanostructures. Acc. Chem. Res. 2008, 41, 1674−84. (11) Ulijn, R. V.; Smith, A. M. Designing Peptide Based Nanomaterials. Chem. Soc. Rev. 2008, 37, 664−75.

■

CONCLUSION In this work, we employed SMD molecular dynamics to investigate the PMF for the transition between the bound and free states of a PA that is pulled from a cylindrical micelle fiber structure that is constructed from the same PAs. By doing a thermodynamics analysis, it is found that the free energy, enthalpy, and entropy differences between the bound and free states are −67 kcal/mol, −71.5 kcal/ml, and −14.5 cal/(mol K), respectively, for the fiber + PA → fiber assembly process. This indicates that the driving force for self-assembly is enthalpic, with both electrostatic and van der Waals interactions playing important roles. The electrostatic interaction is dominated by the charged head segment of the PA interacting with other PAs and with solvent. However, for the van der Waals interaction both the head and tail segments of the PA are important, with interactions with the fiber decreasing as the PA emerges and interactions with the solvent increasing. Overall, the mechanistic picture that emerges from this work is that as the PA is pulled from the fiber, the interaction evolves through three stages: (1) initially electrostatic interactions between the charged head of the pulled PA and PAs in the fiber, and between the PA and solvent are dominant, (2) after the charged head emerges, the rest of the peptide comes out, with both PA−solvent electrostatic interactions and van der Waals interactions being significant, and (3) in the last step, the alkane tail emerges, dominated by van der Waals interactions with either peptide or alkane. The first two steps involve 7459

dx.doi.org/10.1021/jp401508w | J. Phys. Chem. A 2013, 117, 7453−7460

The Journal of Physical Chemistry A

Article

(12) McCullagh, M.; Prytkova, T.; Tonzani, S.; Winter, N. D.; Schatz, G. C. Modeling Self-Assembly Processes Driven by Nonbonded Interactions in Soft Materials. J. Phys. Chem. B 2008, 112, 10388−98. (13) Tsonchev, S.; Niece, K. L.; Schatz, G. C.; Ratner, M. A.; Stupp, S. I. Phase Diagram for Assembly of Biologically-Active Peptide Amphiphiles. J. Phys. Chem. B 2007, 112, 441−47. (14) Tsonchev, S.; Schatz, G. C.; Ratner, M. A. HydrophobicallyDriven Self-Assembly: A Geometric Packing Analysis. Nano Lett 2003, 3, 623−26. (15) Tsonchev, S.; Schatz, G. C.; Ratner, M. A. ElectrostaticallyDirected Self-Assembly of Cylindrical Peptide Amphiphile Nanostructures. J. Phys. Chem. B 2004, 108, 8817−22. (16) Tsonchev, S.; Troisi, A.; Schatz, G. C.; Ratner, M. A. On the Structure and Stability of Self-Assembled Zwitterionic Peptide Amphiphiles: A Theoretical Study. Nano Lett. 2004, 4, 427−31. (17) Solis, F. J.; Stupp, S. I.; de la Cruz, M. O. Charge Induced Pattern Formation on Surfaces: Segregation in Cylindrical Micelles of Cationic−Anionic Peptide-Amphiphiles, J. Chem. Phys. 2005, 122. (18) Velichko, Y. S.; Stupp, S. I.; de la Cruz, M. O. Molecular Simulation Study of Peptide Amphiphile Self-Assembly. J Phys Chem B 2008, 112, 2326−34. (19) Zhang, S. M.; Greenfield, M. A.; Mata, A.; Palmer, L. C.; Bitton, R.; Mantei, J. R.; Aparicio, C.; de la Cruz, M. O.; Stupp, S. I. A SelfAssembly Pathway to Aligned Monodomain Gels. Nat. Mater. 2010, 9, 594−601. (20) Lee, O. S.; Stupp, S. I.; Schatz, G. C. Atomistic Molecular Dynamics Simulations of Peptide Amphiphile Self-Assembly into Cylindrical Nanofibers. J. Am. Chem. Soc. 2011, 133, 3677−83. (21) Lee, O. S.; Cho, V.; Schatz, G. C. Modeling the Self-Assembly of Peptide Amphiphiles into Fibers Using Coarse-Grained Molecular Dynamics. Nano Lett. 2012, 12, 4907−13. (22) Lee, O. S.; Liu, Y. M.; Schatz, G. C. Molecular Dynamics Simulation of Beta-Sheet Formation in Self-Assembled Peptide Amphiphile Fibers, J. Nanopart. Res. 2012, 14. (23) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free Energy Calculation from Steered Molecular Dynamics Simulations Using Jarzynski’s Equality. J. Chem. Phys. 2003, 119, 3559−66. (24) Park, S.; Schulten, K. Calculating Potentials of Mean Force from Steered Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120, 5946−61. (25) Calderon, C. P. On the Use of Local Diffusion Models for Path Ensemble Averaging in Potential of Mean Force Computations, J. Chem. Phys. 2007, 126. (26) He, C. Z.; Genchev, G. Z.; Lu, H.; Li, H. B. Mechanically Untying a Protein Slipknot: Multiple Pathways Revealed by Force Spectroscopy and Steered Molecular Dynamics Simulations. J. Am. Chem. Soc. 2012, 134, 10428−35. (27) Lu, H.; Schulten, K. Steered Molecular Dynamics Simulations of Force-Induced Protein Domain Unfolding. Proteins 1999, 35, 453−63. (28) Sharma, D.; Perisic, O.; Peng, Q.; Cao, Y.; Lam, C.; Lu, H.; Li, H. B. Single-Molecule Force Spectroscopy Reveals a Mechanically Stable Protein Fold and the Rational Tuning of Its Mechanical Stability. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 9278−83. (29) Guzman, D. L.; Randall, A.; Baldi, P.; Guan, Z. B. Computational and Single-Molecule Force Studies of a Macro Domain Protein Reveal a Key Molecular Determinant for Mechanical Stability. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 1989−94. (30) Shen, M. L.; Guan, J.; Xu, L. A.; Yu, Y. Y.; He, J. W.; Jones, G. W.; Song, Y. T. Steered Molecular Dynamics Simulations on the Binding of the Appendant Structure and Helix-Beta 2 in DomainSwapped Human Cystatin C Dimer. J. Biomol. Struct. Dyn. 2012, 30, 652−61. (31) Hwang, H.; Schatz, G. C.; Ratner, M. A. Steered Molecular Dynamics Studies of the Potential of Mean Force of a Na+ or K+ Ion in a Cyclic Peptide Nanotube. J. Phys. Chem. B 2006, 110, 26448−60. (32) Khurana, E.; Nielsen, S. O.; Ensing, B.; Klein, M. L. SelfAssembling Cyclic Peptides: Molecular Dynamics Studies of Dimers in Polar and Nonpolar Solvents. J. Phys. Chem. B 2006, 110, 18965−72.

(33) Cheng, C. L.; Zhao, G. J. Steered Molecular Dynamics Simulation Study on Dynamic Self-Assembly of Single-Stranded DNA with Double-Walled Carbon Nanotube and Graphene. Nanoscale 2012, 4, 2301−05. (34) Jarzynski, C. Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett. 1997, 78, 2690−93. (35) Jarzynski, C. Equilibrium Free-Energy Differences from Nonequilibrium Measurements: A Master-Equation Approach. Phys. Rev. E 1997, 56, 5018−35. (36) Grubmuller, H. Solvate, 1.2 ed.; Theoretical Biophysics Group, Institute for Medical Optics, Ludwig-Maximilian University: Munich, 1996. (37) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics Modell. 1996, 14, 33−38. (38) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−35. (39) Mackerell, A. D.; Wiorkiewiczkuczera, J.; Karplus, M. An AllAtom Empirical Energy Function for the Simulation of Nucleic-Acids. J. Am. Chem. Soc. 1995, 117, 11946−75. (40) Kale, L.; Skeel, R.; Bhandarkar, M.; Brunner, R.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. NAMD2: Greater Scalability for Parallel Molecular Dynamics. J. Comput. Phys. 1999, 151, 283−312. (41) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant-Pressure Molecular-Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177−89. (42) Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. Constant-Pressure Molecular-Dynamics Simulation - The Langevin Piston Method. J. Chem. Phys. 1995, 103, 4613−21. (43) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald - An N·log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−92. (44) Andersen, H. C. Rattle: A “Velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. J. Comput. Phys. 1983, 52, 24−34.

7460

dx.doi.org/10.1021/jp401508w | J. Phys. Chem. A 2013, 117, 7453−7460