Free Energy Profile and Mechanism of Self-Assembly of Peptide

Article pubs.acs.org/JPCB

Free Energy Profile and Mechanism of Self-Assembly of Peptide Amphiphiles Based on a Collective Assembly Coordinate Tao Yu and George C. Schatz* Department of Chemistry, Northwestern University, 2145 Sheridan Road, Evanston, Illinois 60208-3113, United States S Supporting Information *

ABSTRACT: By combining targeted molecular dynamics (TMD) simulations, umbrella sampling, and the weighted histogram analysis method (WHAM), we have calculated the potential of mean force (PMF) for the transition between the bound and free states of 90 peptide amphiphiles (PAs) in aqueous solution, with the bound state corresponding to a cylindrical micelle fiber. We specifically consider a collective reaction coordinate, the radius of gyration of the PAs, to describe assembly in this work. It is found that the free energy, enthalpy, and entropy differences between the free and bound states are −126 kcal/mol, −185 kcal/mol, and −190 cal/(mol K), respectively, for the self-assembly process. This indicates that the driving force to form the micelle structure is enthalpic. The enthalpic driving forces originate from several factors, including the conformational energy of PAs and the electrostatic and van der Waals interaction energy between solvent molecules and between solvent and PAs. Among these interactions, the solvent electrostatic interaction is the dominating one, contributing 54% of the total driving force. The PMF profile can be recognized as involving two stages of assembly: (1) PAs initially approach each other in mostly random configurations and loosely aggregate, resulting in significant desolvation and initiation of head−tail conformational reorganization; (2) PAs undergo a conformational disorder-toorder transition, including forming secondary structures that include more β-sheets and fewer random coils, along with tail−head core−shell alignment and condensation that leads to total exclusion of water from the core. The PMF decreases slowly in the first stage, but rapidly in the second. This study demonstrates a hierarchy of assembly steps in which PA structural changes, solvation, and redistribution of solvent molecules play significant roles in the PA self-assembly process.

■

INTRODUCTION Peptide amphiphiles (PAs) have been being widely used to construct various functional supramolecular structures.1−6 One of these promising supramolecular structures, cylindrical micelle nanofibers, has attracted lots of attention for scientific research, as the residue sequence can be programmed for biological activity that leads to applications in bone healing, wound repair, and nerve growth. The cylindrical micelle structures are built up through self-assembly processes (Figure 1a) in which the peptide amphiphiles spontaneously aggregate in aqueous solution under appropriate pH values or ion concentration conditions, as reported by Stupp and coworkers,3,4,7 and others.5,8−11 This self-assembly process occurs because the PAs are amphiphilic, consisting of a hydrocarbon tail segment and a peptide head segment. In the cylindrical fiber structure, the hydrocarbon tail segments aggregate together to form a hydrophobic core, while the peptide head segments form a jacket shell around the core. The peptide segments are partially solvated by water, and they also form different secondary structures, such as β-sheet secondary structures that are important to fiber stability. Many theoretical studies have been performed to understand (and predict) PA self-assembly;12−19 however, the conditions of assembly are complex, with sensitivity to choice and number of residues in addition to alkane tail length, pH, temperature, and salt concentration. In addition, there are no solid conclusions as © 2013 American Chemical Society

to what type of interactions (hydrophobic effect, hydrogen bonding leading to β-sheets, and electrostatic interactions) play a more important role in the self-assembly process. Recently, Lee et al. performed a fully atomistic-level simulation of a PA fiber based on a seeded molecular dynamics technique,20,21 and the work reveals atomistic structural details in the description of micelles based on the peptide sequence: SLSLAAAEIKVAV. Furthermore, a coarse-grained model22 was developed to simulate the self-assembly process starting from a homogeneous mixture of PAs in aqueous solution. This study suggested that spherical micelle structures are initially formed, and then these micelles merge into cylindrical micelles and finally aggregate to fibers. Although these simulation studies provide structural information for the PA nanofiber as well as possible mechanisms of the assembly process, they do not address the following fundamental but crucial questions: (1) What is the free energy landscape in self-assembly process? (2) Are the driving forces for assembly enthalpic or entropic? (3) What type of interactions, electrostatic, van der Waals, or others, dominate assembly from a homogeneous aqueous solution? Addressing these questions is important for understanding the Received: May 16, 2013 Revised: July 2, 2013 Published: July 3, 2013 9004

dx.doi.org/10.1021/jp404835q | J. Phys. Chem. B 2013, 117, 9004−9013

The Journal of Physical Chemistry B

Article

radius of gyration of the PAs as the prior defined RC, labeled as Rg in the present work. It is found that the free energy change for the self-assembly process is 126 kcal/mol for 90 PAs, demonstrating that the nanofiber structure is more stable, but only by a few kT per PA. In addition, we carry out a thermodynamic analysis which shows that the changes in free energy, enthalpy, and entropy are −126 kcal/mol, −185 kcal/ mol, and −190 cal/(mol K), respectively. Thus PA selfassembly is an enthalpy-driven process, and we find that the driving forces include the conformational energy of the PAs and electrostatic and van der Waals interaction energy between solvent molecules and between solvent and PAs. Among these interactions, the solvent electrostatic interaction is the dominating one, contributing 54% out of the total driving forces.

■

COMPUTATIONAL DETAILS 1. Constructing the Nanofiber Structures. The peptide sequence is chosen as Ser-Leu-Ser-Leu-Ala-Ala-Ala-Glu-Ile-LysVal-Ala-Val; see Figure 1b for a detailed structure. This is the same as was used in earlier work, and the fiber micelle structure was constructed using the same seeded procedure as in Lee et al.21 wherein circular disks of PAs are constructed to initiate formation of the micelle structure. Each circular disk is defined with nine PAs radially placed on a plane (orange 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 (gold 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 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 molecular dynamics calculations. The 90 PA molecules in the seed structure were solvated in a water box using the SOLVATE21 module31 implemented in VMD.32 Periodic boundary conditions were used in the simulations, using a box of dimensions of 190, 190, and 48 Å in x, y, and z directions, respectively. The box was filled with 47 130 water molecules based on the modified TIP3P potential.33,34 Note (Figure 1b) that there are three charged residues close to the C-terminal in the head segment of each PA; they are Glu (residue 8), Lys (residue 10), and Val (residue 13, as C-terminal). The net charge of the PA is therefore −1. To neutralize the system, 90 Na+ ions are added. In addition to these sodium ions, another 133 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 high-energy contacts. Molecular dynamics simulations were carried out using NAMD2.9.35 A 2 ns molecular dynamics simulation at 310 K with an 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 dynamics36 at a temperature of 310 K with a damping coefficient γ = 5 ps−1. Pressure was maintained at 1 atm using the Langevin piston method36,37 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

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) The initial structure of the PAs for MD simulations. Nine PAs are placed radially in the first layer (orange) where their tails are pointing inward. The second layer also has nine PAs (gold) and is rotated by 20° relative to the first layer. The angle between neighboring PAs is 40°, and the distance between layers is 5 Å.

mechanism of self-assembly, and especially in the further design of experiments to generate more complex nanostructures that are generated from mixtures of different PAs or of PAs with other molecules. Very recently, Yu et al. studied the potential of mean force (PMF) or free energy profile for the assembly of a single PA into a cylindrical micelle using steered molecular dynamics simulations (SMD).23 A pairwise interaction analysis in that work demonstrated 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. This work provides insight to understand the energy landscape and driving forces for the PA assembly process; however, the conclusions are limited due to the relatively simple model used in which the PMF depends on the coordinates of a single PA rather than describing the simultaneous assembly of many PAs. The seeded molecular dynamics studies of Lee and coworkers20,21 started with a structure that was already partly assembled, but within this constraint all PAs in the calculation underwent collective assembly to give the final micelle structure. The coarse-grained calculations22 mentioned previously showed collective assembly from homogeneous solution, although initially involving spherical structures and then cylindrical structures. We infer from this that the assembly process is likely collective, so using a collective coordinate is likely to give the most realistic description of self-assembly. However, even with a collective coordinate, one can imagine different assembly mechanisms, such as (1) a one-step model, where PAs that are homogeneously distributed in solution (free state) aggregate collectively to form a fiber structure (bound state); (2) a two-step model, where the PAs form intermediate states, for instance, spherical micelles, and then these aggregate to the final fiber structure. In the present work, we employed both targeted MD simulations24−26 and umbrella sampling27−29 to calculate the PMF based on a collective reaction coordinate (RC). We are inspired by previous studies of protein folding;30 we choose the 9005



Article

method with a 1 Å grid width.38 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,39 allowing a 2 fs time step. Atomic coordinates were saved every 10 ps for the trajectory analysis. 2. Constructing the Disperse Structure of PAs. We built up a model disperse structure of PAs to mimic a homogeneous configuration of fully solvated PAs. First, we equilibrated a single PA in a water box of dimensions of 90, 60, and 60 Å in x, y, and z directions for 100 ps, under the same conditions as for the assembled micelle structure, with 0.15 mol/L NaCl at 310 K. We then duplicated this equilibrated PA structure to generate 90 PAs that were randomly distributed in a box with dimensions 150, 150, and 48 Å in the x, y, and z directions. Then these randomly distributed PAs were transferred into a water box of dimensions 190, 190, and 48 Å. (The larger box size in the x−y directions was used to eliminate interaction between PAs in the simulation box and in its periodic images, and to avoid the PAs diffusing out of the simulating box during the simulation. The radius of gyration analysis described below is simplified with these assumptions.) This box was filled with 42 982 water molecules and 211 Na+ and 121 Cl− to neutralize the system and make the NaCl concentration be 0.15 mol/L. A 2 ns molecular dynamics simulation at 310 K with an NPT ensemble was performed to equilibrate the system and the equilibrated PA structure was used as the reference structure in the targeted molecular dynamics simulation discussed in the next section. 3. Umbrella Sampling. For umbrella sampling, the radius of gyration was selected as the reaction coordinate, to characterize the transition between the ordered nanofiber structure and disordered disperse PA structure. In particular, the radius of gyration Rg was defined based on selecting all the CA carbons in the head part and all the carbons in the tail segment of the PAs. On the basis of the previous steps, the fiber and disperse structures of PAs were used as the initial and final structures for the umbrella sampling, with radii of gyration of 27.2 and 47.9 Å, respectively. Note that the correlation between PAs in the fiber and disperse structures is chosen randomly. To generate intermediate structures with Rg that were between 27.2 and 47.9 Å, we employed a targeted molecular dynamics (TMD) simulation25 in which the transition between targeted structures is constrained by the tuning the root-mean-square deviation (rmsd). In TMD, a harmonic potential is added to the system as defined by U (X , t ) =

1 Nk[RMSD(X ) − RMSD(t )]2 2

In the umbrella sampling, we divided the Rg simulation range between 27 and 48 Å into 105 windows, such that the separation between centers of the neighboring windows was 0.2 Å. From trajectory files obtained in the TMD simulation, we selected configurations whose Rg values were 27.1, 27.3, 27.5, ..., 47.5, 47.7, and 47.9 Å, respectively. The constraint on Rg at the center of each window in the umbrella sampling was achieved by adding harmonic potentials to all the CA carbons in the head part and all the carbons in the tail segment of the PAs. The harmonic potential is defined as follows Ui =

1 k[R g(t ) − R gi ]2 2

(2)

Rig

where k, Rg(t), and are the force constant, radius of gyration at time t, and the radius of gyration at the center of ith window, and i varies from 1 to 105. In the simulation, k was set to 500 kcal/mol, and Rig was equal to 27.1, 27.3, 27.5, ..., 47.5, 47.7, and 47.9 Å, respectively. The umbrella sampling was carried out in two stages. First, a 2 ns NPT ensemble simulation for each window was performed to equilibrate the window at 310 K and 1 atm. Next, in the production period another NPT ensemble 2 ns simulation was performed to collect the Rg distribution in each window at the same temperature and pressure, with the Rg value saved every 40 fs. Once again, all the simulation parameters and details are the same as in the previous description of the nanofiber simulations. The total simulation time in the umbrella sampling step is 420 ns. 4. Free Energy Calculation Using WHAM. The weighted histogram analysis method40 was employed to analyze the Rg distribution data as obtained by the umbrella sampling at 310 K. In the WHAM calculation using the program written by Grossfield,41 the upper and lower boundaries were set to 27 and 48 Å. The number of bins was equal to 840, and the convergence tolerance was set to 0.00001.

■

RESULTS AND DISCUSSION Figures 2 and 3 display snapshots of side and top views of the micelle, showing the movement of PAs from the fiber structure,

(1)

where N, k, RMSD(X), and RMSD(t) are the number of atoms constrained, force constant for the constraint on each atom, the rmsd value of the configuration X with respect to the reference structure, and the targeted rmsd that decreases to zero linearly with time. In this study, all the CA carbons in the head part and all the carbons in the tail segment of PAs were selected for applying the constraint so that N is equal to 2520 in the simulation. The force constant k was set as 100 kcal/mol Å2. The fiber and disperse structures obtained in the previous steps were used as the starting point and reference structure, respectively. The TMD simulation was performed to drive the PAs from fiber structure to disperse structure in 10 ns.

Figure 2. Representative configurations taken from MD simulations of PAs transferring from the bound state to the free state with different radii of gyration equal to (a) 27.2, (b) 30.9, (c) 34.9, (d) 38.9, (e) 42.9, and (f) 47.9 Å. Note that three unit cells are pictured, and the red and blue represent the tail and head segment, respectively. 9006



Article

These configurations reveal important details of the disassembly process. For example, in comparing the configurations in Figure 3, parts a and b, it is clear that the initial stages of disassembly involve asymmetrical deformation of the fiber core into a more flattened structure in the xz plane. This structural change is found to be associated with an asymmetric PA secondary structure distribution in the bound state as shown in Supporting Information Figure S2. This figure shows that the β-sheets are radially directed and angularly asymmetric in the bound state. The highlighted range (red in Supporting Information Figure S2) along the x-axis shows PAs that are random-coil structures. These are more flexible and thus are easier to be deformed in the transition. After the initial flattening of the fiber in going from configuration a to b of Figure 3, the PAs gradually separate (mostly) in going to c to d to e. Only in Figure 3f have the tail regions become completely separated, as will be further explored in discussions below. Another trend that is clear in Figure 3c−e is that a few PAs dissociate completely from the micelle at an early stage of disassembly, basically following the isolated PA dissociation pathway that was studied in our earlier work.23 This indicates that there is a distribution of disassembly pathways possible for a given radius of gyration, with a few PAs undergoing complete dissociation from the micelle while most undergo gradual dissociation through various levels of clustering. This is a more detailed picture than was found in the coarse-grained modeling,22 where only the small cluster component of the disassembly process was observed. One important structural change in going from the bound state to the free state involves the relative positions of the tail (red) and head (blue) segments of the PAs in Figures 2 and 3. In the bound state, the tail and head segments of the PAs form an ideal “core−shell” alignment, where all the tail segments are packed together in the hydrophobic core of the fiber, with no water present, and all the peptide head segments surround the core and form a solvated hydrophilic shell. By contrast, in the free state, this well-organized tail−head pattern is broken down, and both the tail and head segments are randomly distributed and solvated. To understand the structural history for the tail− head order-to-disorder transition process, we investigated the correlations between the radius of gyration of the head segments of the PAs, labeled as Rhg, and the radius of gyration of the tail segments of the PAs, labeled as Rtg. Figure 4a plots these two variables as a function of overall radius of gyration Rg. This shows that in the bound state Rhg is larger than Rg, while Rg is larger than Rtg; and in the free state Rhg, Rtg, and Rg are roughly equal to each other. These results reflect the ordered head−tail core−shell architecture in the bound state and the disordered head−tail random distribution in the free state, respectively. Parts a and b of Figure 4 indicate that the boundto-free transition process can be roughly viewed as involving two stages. In the first stage, where Rg is between 27 and 40 Å, Rhg is greater than Rtg. Both Rhg and Rtg linearly increase with Rg, but Rtg increases faster than Rhg, as shown in Figure 4a. This demonstrates that the tail segment displays a larger structural change compared to the head segment in this stage. In the second stage between 41 and 48 Å, the trends of Rhg and Rtg reverse, with Rhg is a little smaller than Rtg. Also, Rtg increases more slowly than Rhg. Figure 4b is generated from Figure 4a by renormalizing Rhg and Rtg using the definition (Rhg − min[ Rhg ])/(max[ Rhg ] − min[ Rhg ]) and (Rtg − min[ Rtg ])/ (max[ Rtg ] − min[ Rtg ]), respectively. This plot shows quantitative structural changes for head and tail segments, and

Figure 3. Representative configurations taken from MD simulations of PAs transferring from the bound state to the free state with different radii of gyration equal to (a) 27.2, (b) 30.9, (c) 34.9, (d) 38.9, (e) 42.9, and (f) 47.9 Å. Here one unit cell is pictured along the fiber axis, and the red and blue represent the tail and head segment, respectively.

as shown in Figures 2a and 3a and denoted as the “bound state”, into the disperse structure, as shown in Figures 2f and 3f and denoted as the “free state”. Here and throughout the rest of the paper we will describe the process in terms of disassembly (from bound to free states) rather than assembly. Figures 2 and 3b−e show intermediate configurations between the bound and free states for various radii of gyration. Note that these configurations were obtained first by the TMD simulation, followed by equilibrating 2 ns with harmonic forces to constrain the corresponding Rg values. The TMD results are based on changes in the RMSD with respect to the free state; hence, they are highly dependent on the which configuration was chosen for the free state. Since the free state involves a random distribution of the 90 PAs, there are tremendous number of configurations which could be used as the reference frame for the TMD calculation, but only one of them was used in the present work. Consequently, the intermediate structures generated by TMD based on the one RMSD reference configuration and a relatively short integration, 10 ns, might not sample a large enough RMSD configurational space. However, this drawback is reduced by the subsequent equilibration step. By equilibrating the TMD generated intermediate states using constraints for the radius of gyration, the PAs can visit a large RMSD configurational space, as shown in Supporting Information Figure S1 (and this plot also shows that the RMSD has largely stabilized after the 2 ns equilibration for each Rg value). This demonstrates that the selected intermediate states can sample further RMSD changes subject to the constrained Rg values, thus providing flexibility in states associated with the transition between bound state and free states. 9007



Article

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

Figure 4. (a) Radius of gyration of the head segment and tail segments of the PAs as a function of the overall PA radius of gyration. (b) Renormalized results from panel a using the definition given in the text.

we see that the transition of head and tail segments to their free state is 41% and 80% complete, respectively, in the first stage. Note that the bound-to-free disassembly not only makes an order-to-disorder transition for the head−tail organization but also leads to changes in the secondary structure of the head segment in the PAs. Figure 5 shows the averaged secondary structure distribution for the bound and free states using 25 frames for each state. This demonstrates that there is a significant secondary structure change between the two states, with both the turn and β-sheet structures decreasing a lot from the bound to the free state, but not completely disappearing in the free state. By contrast, the random-coil structure increases a lot from the bound to the free state. In addition, we monitored the dynamical history of the secondary structural change during the bound to free transition. Figure 6 shows that both the increasing of the disordered structure (random coil) and decreasing of the ordered structure (β-sheet strand and turn) along the Rg can be thought of as occurring in two stages. In stage one, where Rg is between 27 and 35 Å, the secondary structures change dramatically; while in stage two (Rg is between 35 and 48 Å), the secondary structures show small changes with random fluctuations. This indicates that the ordered secondary structures are stable only when the PAs are close to each other. In addition, comparing Figures 4 and 6, we see a possible reason why Rhg increases slowly in the small Rg range. In this range, the ordered secondary structures lock the interaction between PAs and make the head segment of the PAs more rigid. Note that this dynamical secondary-structurechange feature can only be observed using an atomistic-level simulation; the previous coarse-grained model simulation22

Figure 6. Percent of ordered (β-sheet and turn) and disordered (random coil) secondary structure as a function of PA radius of gyration.

cannot reveal this point because it does not allow relaxation of the secondary structure due to a restriction of the Martini force field.42,43 Associated with the large structural reorganization of PAs, the solvent molecules also undergo a large redistribution in order to adapt to the PA structural changes in the bound-tofree state transition process. Figure 7 shows the radial and number distribution (calculated by integrating the radial distribution) of solvent molecules, including water and ions, with respect to the center of the water box for different boundto-free transition stages that are labeled according to their Rg values. Note that in the initial bound state, all the PAs are condensed in the center of the simulation water box, so water is only found in the peptide and outside the micelle. Once the PAs start to transition to the free state by moving toward the edge of the water box, the center of the fiber starts to be 9008



Article

accessible to solvent molecules. Thus as expected, the distribution shifts left with increasing Rg. Note that the radial distribution displays big jumps left in the 5 and 35 Å radius range corresponding to the change Rg = 42.9 Å to Rg = 47.9 Å. This result agrees with the behavior in Figure 3, going from configuration e to f, where the PAs are diluted by solvent and somewhat reorganized. This solvation of the PAs enhances PA−solvent interactions but reduces solvent−solvent interactions as we discuss further below. Another way to study the correlations between PAs and solvent is provided by examining the solvent-accessible surface area (SASA). Supporting Information Figure S3 shows the SASA for the same group of selected PAs in the bound and free state, calculated using Shrake−Rupley algorithm.44 It is found that the both the head and tail segments are much better solvated in the free state compared with the bound state. The dynamical history of the solvation effect along the reaction coordinate is displayed in Supporting Information Figure S4. It is found that the SASA of both head and tail segments continuously increases as Rg increases, although there is a slight decrease for the tail segment near the free state configuration. This behavior indicates that the PAs are gradually exposed to solvent as the disassembly occurs, except that at the very end the PAs are able to refold to give a slight decrease. The PMF profile as a function of Rg is displayed in Figure 8a. This demonstrates that in this one-step transition model the bound state is more stable than the free state, with a free energy difference of about 126 kcal/mol, corresponding to 1.4 kcal/ mol per PA. The PMF profile can be viewed as involving two stages. In the first stage, for Rg between 27 and 37 Å, the PMF monotonically increases and reaches a local maximum point at Rg = 36.7 Å where the PMF is 114 kcal/mol. Thus, 93% of the

Figure 7. Radial distribution (a) and number distribution (b) of solvent molecules (water, Na+, and Cl−) with respect to the center of the simulation water box.

Figure 8. Potential of mean force (PMF) profile, enthalpy change (total interaction energy change), and the entropy change along the reaction coordinate Rg of the PAs. The PMF was calculated using umbrella sampling method. 9009



Article

Figure 9. PA−PA interaction (a), solvent molecule interaction (b), and PA−solvent interaction (c) along the reaction of coordinate, Rg.

the current NPT ensemble simulation, T is constant and the change in the volume is negligible; therefore, the terms ΔEkinetic(Rg) and PΔV(Rg) can be considered to be zero in the calculations. Equation 3 can then be reduced to following:

total free energy change for disassembly is expended to get to this state. Note that, in this range, the PAs display a large amount of structural changes, including both head−tail segment reorganization and secondary structure relaxation, as shown in Figures 3 and 4. Therefore, PA structural changes are one of the reasons for the free energy cost in this stage. In the second stage, where Rg is in the range between 37 and 48 Å, the PMF profile displays a relatively flat profile increasing by around 11 kcal/mol between 37 and 47.9 Å corresponding to 8% of the total free energy change. Although the PMF is relatively flat in the second stage, other properties show stronger variation. This is clear in Figures 6 and 7, and we will see more examples below. The PMF (free energy) profile is determined by both enthalpy and entropy changes along the reaction coordinate. Considering here the self-assembly process, i.e., free-to-bound state transition, the relationship between the thermodynamics quantities (G, H, and S) in the NPT ensemble is given by

ΔG(R g) = ΔU (R g) − T ΔS(R g)

(4)

where ΔU(Rg) can be obtained directly from the average potential energy, and then ΔS(Rg) is determined from eq 4 and the PMF. On the basis of these calculations, the enthalpy and entropy components are plotted versus Rg in Figure 8, parts b and c. These results show different variation with Rg than we saw for the PMF, with one stage between 27 and 43 Å, and the other between 43 and 48 Å. In the first stage, the interaction energy and entropy components continuously increase and reach maxima near Rg = 43 Å. In the second stage, they both decrease dramatically between 43 and 45 Å and fluctuate as the free state is approached between 45 and 48 Å. Between these two stages, there exists an apparent transition state near Rg = 43 Å, which creates both enthalpy and entropy barriers for the bound-to-free state transition process. Comparing the free energy profile with that for enthalpy and entropy (multiplied by temperature), it is found that ΔU(Rg) and TΔS(Rg) are similar in magnitude to each other but considerably larger (factor of 20) than ΔG(Rg). This result demonstrates that there exists a significant enthalpy−entropy compensation effect45 for the selfassembly process, by which the enthalpy and entropy changes balance each other, leading to a free energy change that is smooth near Rg = 43 Å. In addition, the enthalpy and entropy changes for self-assembly can be calculated, and the resulting ΔH and ΔS are −185 kcal/mol and −190 cal/(mol K), respectively. As expected, the enthalpy is larger in the bound

ΔG(R g) = ΔH(R g) − T ΔS(R g) = ΔE(R g) + P ΔV (R g) − T ΔS(R g) = ΔE kinetic(R g) + ΔU (R g) + P ΔV (R g) − T ΔS(R g) (3)

where ΔG(Rg), ΔH(Rg), ΔEkinetic(Rg), ΔU(Rg), and ΔS(Rg) are the differences in free energy, enthalpy, kinetic energy, interaction energy, and entropy between bound and free states as a function of Rg. T, P, ΔV(Rg) are the temperature, pressure, and change of volume of the system in the simulation. Under 9010



Article

Table 1. Interaction Energies for the Free and Bound Statesa free state error

21413 5800 −4389 0 −499968 49698 −47539 −4097 −479081

21 42 14 0 129 59 84 17 108

average PAs

solvent

PA−solvent

conformational electrostatic van der Waals conformational electrostatic van der Waals electrostatic van der Waals

total

differenced

bound state

b

c

b

error

21064 −51 −6214 0 −510306 49180 −31603 −1336 −479266

23 11 21 0 113 56 60 12 91

average

c

−349 −5851 −1825 −10338 −519 15936 2761 −185

a

Units are kcal/mol. bThe average values were calculated based on 25 trajectories. cThe error value was the standard deviation of the mean, calculated based on 25 trajectories. dThe difference is defined in terms of bound state minus free state.

we have shown that the PAs become much more dilute (as shown in Figure 3, parts e and f), which means that more solvent molecules interact strongly with the head groups (PA− solvent decreases) and less strongly with each other (solvent− solvent interaction increases). This is why parts b and c of Figure 9 change in opposite directions. The increase in PA−PA and solvent−solvent interactions as Rg increases is responsible for stabilization of the bound state, so we can label them as “bound-state-philic” forces. By contrast, the opposite trend of the PA−solvent interaction along the reaction coordinate indicates that the free state is more favorable, so we call it a “free-state-philic” force. It is the competition (and nearly complete cancellation) between the two types of forces that determines the total interaction pattern along the reaction coordinate. Now consider decomposing UPA and USol into three terms, namely, the contribution from conformational, electrostatic, and van der Waals interaction energies for the PAs and solvent, respectively. Detailed values for each type of energy term are tabulated in Table 1. Note that the solvent includes Na+, Cl−, and water molecules, and the water molecules are represented by the TIP3P water force field33,34 where there is no conformational relaxation. Therefore, the conformational energy of the solvent is zero for both free and bound states. Focusing on the energy difference in the last column of Table 1, the terms with negative signs represent driving forces to form the bound state, and those with positive signs are driving forces to form the free state. The conformational portion of the total difference indicates that there is a great structural reorganization transferring from free to bound states, including both head−tail disorder-to-order alignment and secondary structure formation. Such a structural reorganization makes the electrostatic interaction of the PAs change from repulsion in the free state to weak attraction in the bound state. In addition, the PA atoms are closer to each other in the bound state; thus, the van der Waals energy of the PAs is lower compared with the free state. The positive sign of the PA−solvent portion of the total difference is attributed into the fact that PAs are much better solvated and thus the PA−solvent contact surface is larger in the free state. We have renormalized the bound-state driving forces by dividing each term by the summation of all the terms, and from this we can extract information about the relative magnitude that each type of interaction contributes to the selfassembly process. As shown in Figure 10a, the solvent−solvent electrostatic interaction is assigned as the dominating source of the bound-state driving force. Also, the total electrostatic

state, while the entropy is larger in the free state. This demonstrates that the overall transition from free to bound states of the PAs is an enthalpy-driven process. Furthermore, a large enthalpic barrier for the assembly process is found in Figure 8. Although this is exactly canceled by entropic effects such that at 300 K the free energy profile is flat, this situation could change considerably at lower temperature. Consider the enthalpy and entropy at Rg = 44.1 Å, where the resulting values of ΔH⧧ and ΔS⧧ are 703.0 kcal/mol and 2.4 kcal/K mol, respectively. If ΔH⧧ and ΔS⧧ are temperature-independent, then if the temperature decreases to 284.6 K, then there will be a 20 kcal/mol free energy barrier to form a fiber. In order to understand the enthalpic driving force we decompose the total interaction energy into three components as follows: U = UPA + USol + UPA−Sol

(5)

where U, UPA, USol, and UPA−Sol are total interaction, interaction between PAs, interaction between solvent molecules, and interaction between PAs and solvent molecules, respectively. Figure 9a displays the PA−PA interaction energy as a function of Rg. The result shows a three-stage variation with Rg, with an initial rapid increase up to Rg = 35 Å, then a slower increase to Rg = 44 Å, and then a more rapid increase to the fully dissociated state. The overall increase is 8030 kcal/mol, which is 89 kcal/mol per PA. This PA−PA interaction originates from the conformational (bond, angle, and dihedral angle) energy, the electrostatic energy, and the van der Waals energy, all of which depend directly on the structure of the PAs, but based on our earlier discussion it is clear that the first stage results from dramatic changes in these interactions as the PAs are reorganized, while the second and third stages arise from separation of the PAs and further reorganization of the head− tail segments. Parts b and c of Figure 9 show the interaction between solvent molecules and the interaction between PAs and solvent as Rg is increased. These figures also display three stages, where in stages one and three, these two types of interactions change rapidly, while in the stage two, their change is small. It is found that the solvent−solvent interaction increases a lot (by 10 900 or 121 kcal/mol per PA) during the PAs’ bound-to-free transition process, while the PA−solvent interaction decreases a lot (by 18 700 or 207 kcal/mol per PA) for the PA’s bound-tofree transition. These effects can be attributed to the fact that the PAs are much better solvated when they are in the last stage of disassembly, as has been demonstrated in Figures 3 and 7. Note that, in the last stage, where Rg is between 43 and 48 Å, 9011



Article

ing more secondary structures and tail−head core−shell alignment associated with a rapidly decreasing PMF. In the first stage, the PMF decreases smoothly (and slowly), but we found that the solvent−solvent and PA−solvent interaction energies show additional structure as a function of Rg in which PA−solvent and solvent−solvent interactions change more rapidly at long range, >45 Å, than at shorter range. This longrange variation arises from desolvation of the PA head groups as the PAs come together. An interesting aspect of this behavior is that desolvation occurs at long range, followed by head−tail reorganization, and finally secondary structure formation and PA condensation to remove water from the tail region. The generality of this hierarchy of assembly processes will be important to explore in future work. We conclude by noting that the picture provided by the PMF refers to self-assembly that is locally in equilibrium at each point of assembly; however, the energies released, especially in the second stage, are such that nonequilibrium effects and dissipation are expected to occur. However, the overall assembly to give micelle structures does not appear to end up in kinetic structures given that the final state of micelle assembly is found to be stable in the simulation, so apparently any disequilibrium structures are rapidly healed out.

■

ASSOCIATED CONTENT

* Supporting Information S

Figure S1, the RMSD associated with the 2 ns equilibration process after the TMD calculation; Figure S2, additional information about the spatial secondary structure of the assembling fiber; Figures S3 and S4, properties of the solvent-accessible surface area. This material is available free of charge via the Internet at http://pubs.acs.org.

Figure 10. Percent of each interaction term contributing to the total driving forces for the bound state.

interaction, including both solvent and PAs, contributes to about 86% of the total driving force, as displayed in Figure 10b, while the total van der Waals interaction makes a contribution about 12%.

■

■

CONCLUSION In the present work, we proposed the radius of gyration Rg for describing the reaction coordinate for PA self-assembly, in which the PAs make a transition between free and bound states. By combining TMD simulations, the umbrella sampling method, and the weighted histogram analysis method, the potential of mean force for the transition between the bound and free states is obtained as a function of the reaction coordinate, Rg. By analyzing the interaction energy change between the bound and free states, the enthalpy change for selfassembly was also estimated. It is found that the free energy, enthalpy, and entropy differences between the free and bound states of 90 PAs are −126 kcal/mol, −185 kcal/mol, and −190 cal/(mol K), respectively, for the free-to-bound self-assembly process. This indicates that the driving force to form the fiber structure is enthalpic. The enthalpic driving force originates from the conformational energy of the PAs and from electrostatic and van der Waals interaction energies between solvent molecules and between solvent and PAs. Among these interactions, the solvent electrostatic interaction is the dominating one, contributing to 54% of the total driving force. Overall the mechanistic picture for assembly is that the PAs transfer from free state to bound state occurs in two stages based on the PMF profile: (1) PAs approach each other while aggregation, desolvation, and head−tail conformational reorganization takes place (this stage can be viewed as being dominated by PA movement toward the center of the micelle, while solvent moves away); (2) PAs organize together to make a conformational disorder-to-order transition, including form-

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the National Science Foundation (Grant CHE-1147335) (for methods development) and by the DOE NERC EFRC (DE-SC0000989) (for applications). We thank One-Sun Lee for useful comments.

■

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−4254. (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−3623. (3) Hartgerink, J. D.; Beniash, E.; Stupp, S. I. Self-Assembly and Mineralization of Peptide-Amphiphile Nanofibers. Science 2001, 294, 1684−1688. (4) 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−5138. (5) Shimizu, T.; Masuda, M.; Minamikawa, H. Supramolecular Nanotube Architectures Based on Amphiphilic Molecules. Chem. Rev. 2005, 105, 1401−1443. (6) Vauthey, S.; Santoso, S.; Gong, H. Y.; Watson, N.; Zhang, S. G. Molecular Self-Assembly of Surfactant-Like Peptides to Form 9012



Article

by Umbrella Integration with Harmonic Fourier Beads. J. Am. Chem. Soc. 2009, 131, 1706−1716. (28) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187−199. (29) Wilhelm, M.; Mukherjee, A.; Bouvier, B.; Zakrzewska, K.; Hynes, J. T.; Lavery, R. Multistep Drug Intercalation: Molecular Dynamics and Free Energy Studies of the Binding of Daunomycin to DNA. J. Am. Chem. Soc. 2012, 134, 8588−8596. (30) Whitford, P. C.; Sanbonmatsu, K. Y.; Onuchic, J. N. Biomolecular Dynamics: Order−Disorder Transitions and Energy Landscapes. Rep. Prog. Phys. 2012, 75, 076601. (31) Grubmuller , H. Solvate, 1.2 ed.; Theoretical Biophysics Group, Institute for Medical Optics, Ludwig-Maximilian University: Munich, Germany, 1996. (32) Humphrey, W.; Dalke, A.; Schulten, K. Vmd: Visual Molecular Dynamics. J. Mol. Graphics Modell. 1996, 14, 33−38. (33) 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−935. (34) 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−11975. (35) 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. (36) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant-Pressure Molecular-Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177− 4189. (37) 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−4621. (38) 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−10092. (39) Andersen, H. C. Rattle: A “Velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. J. Comput. Phys. 1983, 52, 24−34. (40) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules 0.1. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (41) Grossfield, A. WHAM, version 2.0.6; Unviersity of Rochester: Rochester, NY, 2012. (42) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The Martini Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812−7824. (43) Monticelli, L.; Kandasamy, S. K.; Periole, X.; Larson, R. G.; Tieleman, D. P.; Marrink, S. J. The Martini Coarse-Grained Force Field: Extension to Proteins. J. Chem. Theory Comput. 2008, 4, 819− 834. (44) Shrake, A.; Rupley, J. A. Environment and Exposure to Solvent of Protein AtomsLysozyme and Insulin. J. Mol. Biol. 1973, 79, 351− 371. (45) Cornish-Bowden, A. Enthalpy−Entropy Compensation: A Phantom Phenomenon. J. Biosci. 2002, 27, 121−126.

Nanotubes and Nanovesicles. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 5355−5360. (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−1200. (8) Buerkle, L. E.; Rowan, S. J. Supramolecular Gels Formed from Multi-Component Low Molecular Weight Species. Chem. Soc. Rev. 2012, 41, 6089−6102. (9) Palmer, L. C.; Stupp, S. I. Molecular Self-Assembly into OneDimensional Nanostructures. Acc. Chem. Res. 2008, 41, 1674−1684. (10) Ulijn, R. V.; Smith, A. M. Designing Peptide Based Nanomaterials. Chem. Soc. Rev. 2008, 37, 664−675. (11) Woolfson, D. N.; Mahmoud, Z. N. More Than Just Bare Scaffolds: Towards Multi-Component and Decorated Fibrous Biomaterials. Chem. Soc. Rev. 2010, 39, 3464−3479. (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− 10398. (13) 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, 054905. (14) 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−447. (15) Tsonchev, S.; Schatz, G. C.; Ratner, M. A. HydrophobicallyDriven Self-Assembly: A Geometric Packing Analysis. Nano Lett. 2003, 3, 623−626. (16) Tsonchev, S.; Schatz, G. C.; Ratner, M. A. ElectrostaticallyDirected Self-Assembly of Cylindrical Peptide Amphiphile Nanostructures. J. Phys. Chem. B 2004, 108, 8817−8822. (17) 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−431. (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−2334. (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.; Liu, Y. M.; Schatz, G. C. Molecular Dynamics Simulation of Beta-Sheet Formation in Self-Assembled Peptide Amphiphile Fibers. J. Nanopart. Res. 2012, 14, 936. (21) 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−3683. (22) 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−4913. (23) Yu, T.; Lee, O. S.; Schatz, G. C. Steered Molecular Dynamics Studies of the Potential of Mean Force for Peptide Amphiphile SelfAssembly into Cylindrical Nanofibers. J. Phys. Chem. A [Online early access]. DOI: 10.1021/jp401508w. Published Online: Mar 19, 2013. http://pubs.acs.org/doi/full/10.1021/jp401508w. (24) Ferrara, P.; Apostolakis, J.; Caflisch, A. Targeted Molecular Dynamics Simulations of Protein Unfolding. J. Phys. Chem. B 2000, 104, 4511−4518. (25) Schlitter, J.; Engels, M.; Kruger, P.; Jacoby, E.; Wollmer, A. Targeted Molecular-Dynamics Simulation of Conformational ChangeApplication to the T ↔ R Transition in Insulin. Mol. Simul. 1993, 10, 291 ff. (26) Ferrara, P.; Apostolakis, J.; Caflisch, A. Computer Simulations of Protein Folding by Targeted Molecular Dynamics. Proteins 2000, 39, 252−260. (27) Khavrutskii, I. V.; Gorfe, A. A.; Lu, B. Z.; McCammon, J. A. Free Energy for the Permeation of Na+ and Cl− Ions and Their Ion-Pair through a Zwitterionic Dimyristoyl Phosphatidylcholine Lipid Bilayer 9013


Free Energy Profile and Mechanism of Self-Assembly of Peptide

Recommend Documents