Molecular Simulation Study of Peptide Amphiphile Self-Assembly

Feb 6, 2008 - Hydrogen bonding occurs when two partially negatively charged atoms share a partially positively charged hydrogen atom.17,18 The strongl...
0 downloads 8 Views 511KB Size
Download PDF
2326

J. Phys. Chem. B 2008, 112, 2326-2334

Molecular Simulation Study of Peptide Amphiphile Self-Assembly Yuri S. Velichko, Samuel I. Stupp, and Monica Olvera de la Cruz* Department of Materials Science and Engineering, Northwestern UniVersity, EVanston, Illinois 60208 ReceiVed: June 7, 2007; In Final Form: NoVember 11, 2007

We study the self-assembly of peptide amphiphile (PA) molecules, which is governed by hydrophobic interactions between alkyl tails and a network of hydrogen bonds between peptide blocks. We demonstrate that the interplay between these two interactions results in the formation of assemblies of different morphology, in particular, single β-sheets connected laterally by hydrogen bonds, stacks of parallel β-sheets, spherical micelles, micelles with β-sheets in the corona, and long cylindrical fibers. We characterize the size distribution of the aggregates as a function of the molecular interactions. Our results suggest that the formation of nanofibers of peptide amphiphiles obeys an open association model, which resembles living polymerization.

I. Introduction Self-assembly of complex molecules is an attractive approach for constructing functional supramolecular structures.1 Selfassembled protein fibers, for example, are involved in many intra- and extracellular functions, including cell deformation, division, and movement. They also contribute to the cell’s structural integrity. Creating three-dimensional cell scaffolds with the biological, chemical, and mechanical properties of native tissue, such as collagen, the major structural component of the extracellular matrix, is necessary to achieve the full potential of regenerative medicine.2,3 A new class of peptide amphiphile molecules (PA) that self-assemble into nanofiber networks has been developed recently, 4-10 which have been shown to promote cell growth. Several other biomedical targets have been also approached by peptide amphiphile self-assembly.8,9 Besides their importance as biomimetic materials, peptide amphiphiles are gaining interest in other material science areas, including electronic materials. The formation of fibers through the self-assembly of peptide amphiphiles is thus a field of modern technological and scientific interest. Peptide amphiphile molecules consist of a hydrophobic alkyl tail, an amino acid hydrophobic peptide block, and a bioactive polar epitope head group. When the peptide sequence includes amino acids with high β-sheet formation propensity, long cylindrical nanofibers are observed. Peptide amphiphile nanofibers usually have a diameter of about 6-8 nm. Their typical length is a few microns. The core of such fibers is composed by the hydrophobic block of the PA molecules. A peptide block, which favors a β-sheet formation,8 stabilizes the cylindrical fiber through formation of a network of hydrogen bonds. The self-assembly of peptide amphiphiles is remarkably tolerant to the changes of the peptide block when a robust β-sheet-forming sequence is included. A number of bioactive head groups have been incorporated without compromising the cylindrical nanostructures.7 The self-assembly of PA molecules is driven mostly by the interplay of hydrogen bonding and hydrophobic forces over a wide range of length scales. The selfassembling property of peptide amphiphiles shows great promise in the development of novel materials. Though hydrogen bonds * To whom correspondence should be addressed. E-mail: m-olvera@ northwestern.edu.

clearly play a key role in the formation of high-aspect-ratio fibers, the self-assembly process is not fully understood. In this paper, we focus on understanding and exploiting the role of molecular interactions in peptide amphiphile selfassembly. The rest of this paper is organized as follows. In section II, we discuss the coarse-grained model and simulation technique we use to study the self-assembly of peptide amphiphile molecules. In section III, we show results and compare our model with some basic models of molecular association. We conclude by highlighting the main results of this study. II. Model and Simulation Technique A. Model of Peptide Amphiphiles. We construct a coarsegrained model for PA molecules using the united atom model,11 which takes into account only the basic chemical structure of the molecules (Figure 1). Nearby similar segments along the chain of the molecule are combined into an effective monomer unit. To connect the coarse-grained model with the atomistic model, we introduce three kinds of monomer units, hydrophobic (H), peptide (P), and epitope (E). The system we examine in our study consists of Np ) 1080 polymer chains, each with Nm ) NH + NP + NE monomer units. For our study, we set NH ) 4, NP ) 5, and NE ) 1, which is close to the composition of most peptide amphiphiles. Each monomer unit is characterized by its diameter σ ) 1 and mass m ) 1. The system is placed into a box of volume L3 with L ) 50σ; periodic boundary conditions are applied in all directions. Nonbonded short-range interactions between hydrophobic (H) monomer units are described by the Morse potential11-13,15

[ (

)

(

)]

δrij δrij - 2exp -R σ σ

UM(rij) ) HH exp -2R

(1)

where δrij ) rij - σ and rij is the distance between the ith and jth monomer units. At rmin /σ ) 1 + ln(2)/R, the potential reaches a minimum that is equal to U(rmin) ) -HH. Other shortrange interactions between monomer pairs (P-P, E-E, P-E, H-P, etc.) include only the excluded volume of magnitude ij ) xiijj approximated using the Lorentz-Berthelot mixing rule11-14, where ii determines the strength of short-range interactions between similar monomers in kBT units. The

10.1021/jp074420n CCC: $40.75 © 2008 American Chemical Society Published on Web 02/06/2008

Peptide Amphiphile Self-Assembly

J. Phys. Chem. B, Vol. 112, No. 8, 2008 2327

Figure 1. (a, b) The chemical structure of peptide amphiphile molecules with oppositely charged head groups. (c) The united atom model of a peptide amphiphile molecule. The hydrophobic tail is colored gray, the β-sheet-forming peptide sequence is in yellow, and the terminal epitope head group is red.

constant R ) 24 determines the shape of the potential and how fast it goes to zero. For the chemically connected monomer units, to maintain the distance between monomers and to prevent chains from crossing each other, a finite extendable nonlinear elastic (FENE) potential 11-13,15 is used

UF(li) )

[ ( )]

kFl2o li ln 1 2 lo

2

(2)

where li ) |r bi - b ri-1| is the bond length and kF ) 50ij/σ2 is the spring constant. The FENE potential is harmonic at its minimum, but the bonds cannot be stretched beyond lo ) 1.5σ. The intrinsic backbone stiffness of the PA molecules11 is incorporated by a three-body interaction potential

[

r i) ) gi 1 + US(b

]

(b r i-1 - b r i)(b r i+1 - b r i) li-1li

2

(3)

where gi is the force constant. For the block of hydrophobic units, gH ) 1, while the peptide block is more rigid and gP ) 100. B. Hydrogen Bonding. The important feature of the interactions between peptide amphiphile molecules is the ability to form a network of hydrogen bonds between peptide units. Single-component peptide amphiphiles assemble only at high salt concentrations or when the epitope head groups are electroneutral (that is, acidic (-) peptide amphiphiles assemble at high pH values and basic (+) peptide amphiphiles at low pH conditions) to overcome electrostatic repulsions. For the sake of simplicity, we consider electroneutral PA molecules, taking into account only hydrophobic interactions and hydrogen bonding. Hydrogen bonding can be understood with rigorous quantum mechanical calculations. Hydrogen bonding occurs when two partially negatively charged atoms share a partially positively charged hydrogen atom.17,18 The strongly electronegative atom attracts the electron cloud from the hydrogen and leaves the atom with a positive partial charge. Hydrogen bonding is frequently considered as an intermolecular force that exists between two partial electric charges of opposite polarity. However, hydrogen bonding has a partial covalent nature. When the hydrogen bond is formed, the redistribution of the electron

Figure 2. The schematic representation of the β-sheet. Three different molecules, i, j, and k, are connected by hydrogen bonds. Covalent bonds are shown as solid lines, while hydrogen bonds are shown by dashed lines. The orientation of the i and j molecules determines the angle ϑij between vectors ∆r bi ) b ri+1 - b ri-1 and ∆r bj ) b rj+1 - b rj-1. Two hydrogen bonds connect molecule i with molecules j and k to form the φi angle.

charge density restricts the possibility of forming another bond. Another important feature of the hydrogen bond is that the formation of the bond depends on the relative orientation of atoms in space.19 From a computational point of view, hydrogen bonding remains a challenge18,20-23. Hydrogen bond models can be classified into two general groups. The first one is geometric in nature. Two molecules are considered to be bound if they are close to each other by a given distance.21,22,24 This condition is frequently coupled with an angular criterion. The second group is of energetic nature. In this case, two molecules are bonded if their interaction energy is stronger than a critical value,18,22,23,25 chosen on the basis of experimental evaluation. Our coarse-grained model of hydrogen bonding includes features of both of these models. Additionally, it requires an activation of the monomer before hydrogen bond formation.18,20,26-28 The geometric criterion limits the intermolecular distance between two P-monomer units to be less than a maximum bond length lo. The energetics of the hydrogen bond formation includes the activation energy Va for initiation of the bonding and intramolecular energetics Vi of the new hydrogen bond. The activation energy depends on the spacial orientation of the donor and acceptor groups and, for simplicity of the model, is defined as Va ) β cos(ϑij), where β is the hydrogen-bonding energy and ϑij is the angle between adjacent monomer segments, that is, cos(ϑij) ) (∆r bi‚∆r bj)/(∆ri∆rj), where ∆r bi ) b ri+1 - b ri-1 (Figure 2). The angular dependence of the activation energy results in the formation of the hydrogen bond only when molecules are aligned in parallel and are oriented in the tail to head group direction. The hydrogen-bonding energy β determines the forming and breaking rate constant of the bond. The strength of the hydrogen bonds between different atoms can vary from a few kJ/mol to a few tenths of a kJ/mol.18 By changing the amino acid sequence of the peptide block, it is also possible to control the strength of the hydrogen bonding between PA molecules.29 Each amino acid includes a hydrogen donor N-H and acceptor CdO group. Each N-H group in the backbone of one amino acid establishes a hydrogen bond with the CdO group in the backbone of another. For the simplicity of our model, we assume that each P-monomer is capable of forming two hydrogen bonds. Therefore, the intramolecular energetics includes a stretching energy of the new bond as well as a bending

2328 J. Phys. Chem. B, Vol. 112, No. 8, 2008

Velichko et al.

energy of the β-sheet if one hydrogen bond is already formed (Figure 2). The hydrogen bond is described as a regular chemical bond when two monomer units are bonded and is modeled using a modified FENE potential (eq 2)

UBond(l) ) UF(l) - UF(σ)

(4)

The assembly of PA molecules into planar structures can be easily introduced into the model by defining a proper value for dihedral angles (Figure 2) of the backbone.15,16 The PA molecules lie in the plane of the β-sheet, and every hydrogen bond connects a donor group of one molecule with an acceptor on the next molecule. The planarity of the β-sheets is enforced by an angular potential applied (eq 3) to the angle φi formed by the two hydrogen bonds (Figure 2) that connect two different molecules with the same monomer unit of the third PA molecule. This angle is determined from cos(φi) ) (∆r bki‚∆r bij)/ (∆rki∆rij), where ∆r bij ) b ri - b rj (Figure 2). This potential determines a bending energy of the β-sheet. The energy of the new bond and the bending energy of the β-sheet constitute the intramolecular energetics Ve of the hydrogen bond. C. Computer Simulation Technique. Many processes in molecular systems take place on different time scales. For example, the bond length fluctuations and torsional angle changes occur on femto- to picosecond time scales.11 On the other hand, the characteristic time of molecular association is determined by a collective energy barrier and varies over a great range.27 Appearance of different time scales significantly complicates molecular simulation, especially for reaching equilibrium in finite computational time scales. Simulation of systems that include interactions of different natures, for example, short-range forces and hydrogen bonding, faces similar challenges. A Monte Carlo method is the simplest way to approach such problems24 due to the possibility to carry out different physically unlikely movements. On the other hand, it excludes collective molecular motion, making it hard to equilibrate large systems. The mixed Monte Carlo-stochastic dynamics method samples local energy minima and enables the system to cross large energy barriers.30,31 A simulation is carried out within the framework of the mixed Monte Carlo-stochastic dynamics method that includes a stochastic dynamics step and a Monte Carlo step in an alternating manner.30,31 The stochastic dynamics step is performed in Cartesian space by solving Langevin’s equations using the Langevin dynamics-velocity Verlet method 11-13,15,31

mi

d2b ri dt

2

)-

db ri ∂U(b,t) r - γ mi + B Γi(t) ∂b ri dt

(5)

where mi is the mass and b ri is the position vector of the ith monomer unit or colloidal particle, U(r b,t) is the interaction potential affecting the monomer unit or the particle at a time t, γ ) 1 is the friction coefficient (damping constant in reciprocal time units), and B Γi(t) is a stochastic force satisfying the conditions of Gaussian white noise 〈Γ Bi(t)Γ Bi(t′)〉 ) 2miγikBTδ(t - t′), where kB is Boltzmann’s constant, T the temperature, and δ(...) the delta-function. Integration is performed with the time step ∆t ) 0.0025 γσ2m/T. In Appendix A, we justify the choice of the simulation parameters. In our simulation, at each time step, an additional Monte Carlo step13,30-32 is performed to create or break a hydrogen bond for each randomly chosen monomer unit. Formation of the bond

is accepted according to the Metropolis criterion,11 with a probability of

Π+{Ri,βj} )

g{Ri,βj} 1 + g{Ri,βj}

>R

(6)

where g{Ri,βj} ) exp(-u{Ri,βj}/kBT), u{Ri,βj} is the bond energy, Ri is a reactive site of the ith monomer unit, and 0 < R < 1 is a random number. The probability of breaking the hydrogen bond is taken to be

Π-{Ri,βj} )

1 1 + g{Ri,βj}

(7)

Each system is equilibrated during τq ) 5 × 104 γσ2m/T time steps. The equilibration process is accompanied by a gradual decrease of temperature (temperature quenching) from Tmax ) 5 to Tmin ) 1.11-13 After reducing the temperature, another τq time steps are allowed to elapse before statistics are collected. Finally, τq time steps are used to perform measurements. III. Results and Discussion A. The Closed and Open Associations Limits. The selfassembly of peptide amphiphile molecules embraces two classical models of aggregation.26 The first model is driven by hydrogen bonding and describes the aggregation of molecules into the one-dimensional β-sheets that can be described using an open association model. It involves step-by-step aggregate growth, A1 + Am f Am+1, where Am stands for a multimer of m unimers A1 that results in a wide distribution of aggregates. The number concentration of multimers obeys the mass action law, Cm+1 ) kcm-1 , where k is the equilibrium constant of the 1 reaction.26-28,33,34 The second model is driven by the attraction between hydrophobic tails that results in formation of finite size, stable, and soluble micellar aggregates.35-37 Formation of monodisperse micellar aggregates can be described in terms of closed association models.37,38 The importance of micelle polydispersity was studied during the same periods.36 In the very low concentration limit, when the entropy dominates over the interaction energy, the amphiphilic molecules interact weakly. Thus, the solution consists of mostly unassociated molecules (unimers). The transition from the unimer phase to the micellar phase occurs at the critical micelle concentration (CMC), Cmc ) k*Cm1 c where k* is an equilibrium constant of the association reaction, mc is a micelle’s aggregation number, and C1 is a concentration of unimers.26 The amphiphilic molecules form micelles with a densely packed hydrophobic core screened from the solvent molecules by the peptide blocks forming the corona of the micelle. The critical micelle concentration depends on both the concentration of molecules, C1, and the strength of the interaction between hydrophobic units, /HH ≡ ln(k*/kBT). A similar transition happens when micellization is investigated as a function of temperature. The temperature at which the formation of micelles occurs is defined as the critical micelle temperature (CMT) at a constant mass concentration.39 The self-assembly of peptide amphiphile molecules is driven mainly by hydrogen bonding and hydrophobic interactions. Competition between these forces shows mixed features peculiar to both association models. Figure 3 illustrates the aggregates size distribution for different values of hydrophobic HH and hydrogen bonding β energies. The calculated aggregate size distribution, Cm, for weakly hydrophobic molecules (HH ) 0.5, Figure 3a) obeys the power low dependence in analogy with

Peptide Amphiphile Self-Assembly

J. Phys. Chem. B, Vol. 112, No. 8, 2008 2329

Figure 3. The aggregate size distribution, Cm, for (a) HH ) 0.5, (b) HH ) 2, and (c) HH ) 3 and different values of hydrogen-bonding energy β. Open association (a) yields a wide distribution of aggregates, in contrast to monodisperse aggregates formed by the (c) closed association.

the open association model. With an increase in the degree of hydrophobicity (HH ) 2, Figure 3b), the distribution becomes wider and deviates from the mass action law. Even at high values of hydrogen-bonding energy, the aggregate size distribution does not reveal any sign of peak formation. With a further increase in the degree of hydrophobicity (Figure 3c, HH ) 3), a clear peak in the aggregate size distribution is observed, which corresponds to the formation of finite size aggregates. Contrary to previous systems, a change in the hydrogen-bonding energy does not substantially change the peak position or the shape of the aggregates distribution for β e 2. On the other hand, for β > 2, the peak becomes wider and shifts to larger m values. Experimental measurements, for example, light scattering, small-angle neutron scattering, and X-ray scattering,35 can only determine the average size of aggregates. The number-average aggregation number, mn, and the weight-average aggregation number, mw, can be calculated using the aggregate size distribution function

mn ≡

mcm ∑ m cm ∑ m

, mw ≡

m2cm ∑ m mcm ∑ m

(8)

Since the simulation is performed at constant concentration and the CMT is a function of the degree of hydrophobicity, HH, the inverse pair interaction energy, 1/HH, measures the relative “distance” from the CMT. The standard critical micelle temperature can be found in the limit of zero hydrogen-bonding energy, β ) 0. Figure 4a shows the weight-average aggregation number as a function of hydrogen-bonding energy β. The transition to the micellar phase (when the distribution Cm versus m shows a second peak as in Figure 3c) occurs in the range 2.5 < /HH < 3, where the weight-average aggregation number mw displays a sharp transition associated with crossing the critical micelle temperature (CMT). The weight-average aggregation

Figure 4. (a) The weight-average aggregation number and (b) the aggregate polydispersity index as a function of hydrogen-bonding energy β for different values of the hydrophobic energy HH. The dashed line corresponds to mw/mn ) 2.

number shows almost no dependence on the hydrogen-bonding energy for β < 1, and therefore, the critical micelle temperature also shows no dependence. With a further increase in the hydrogen-bonding energy, β > 1, the β-sheet size distribution widens (Figure 3b), and the number-average aggregation number, mn, increases. Surprisingly, the mw indicates the strongest dependence on hydrogen-bonding energy near the critical micelle temperature (Figure 4a). Nucleation is the dominant mechanism of micellization. The β-sheets act as seeds of nucleation due to the lower entropy in comparison with single molecules. Simultaneous formation of β-sheets of different sizes results in the appearance of a nontrivial dependence of the aggregates size on both the β and HH energies. For HH e /HH, the weight-average aggregation number mw shows the development of aggregates from the solution of unimers, dimers, and trimers to the aggregates composed of a few hundreds of molecules. On the other hand, above the CMT, HH g /HH, an increase in the hydrogenbonding energy results in steady monotonic aggregate growth. The differences in the nature of the curves give rise to important differences in the nature of the aggregation process for systems with different combinations of hydrogen bonding and hydrophobic interactions. B. Polydispersity of Aggregates. While different associations yield different distributions of aggregates and different average values, the ratio mw/mn (known as the polydispersity index or PDI) is a convenient parameter that can be used to distinguish systems with different mechanisms of aggregation.33,40 The PDI is varying in the range from 1 when all aggregates have the same size to much higher values for polydisperse systems. For the open association scenario, Carother’s approximation33,40 limits PDI to values of mw/mn e 2 and below, while for the closed association scenario, it is characterized by mw/mn = 1.

2330 J. Phys. Chem. B, Vol. 112, No. 8, 2008

Velichko et al.

Figure 5. The aggregate polydispersity index as a function of hydrophobic energy HH. The dashed line corresponds to mw/mn ) 2.

Figures 4b and 5 show the dependence of the polydispersity index on the hydrogen-bonding energy β and hydrophobic energy HH. For small values of hydrogen-bonding energy β < 1, the ratio of mw/mn increases with HH e /HH from 1 to about 2, and then, for HH g/HH, it decreases back to unity again. This suggests the closed-association scenario of micellization and formation of monodisperse micellar aggregates (Figure 6a) above CMT. The peak of PDI around HH = /HH can be explained by the formation of various seeds of nucleation, that is, small nonequilibrium aggregates. On the other hand, for larger values of β > 1, the system consists of β-sheets of different size, and increases in the hydrophobic energy HH result in more complex aggregation. Whereas in the limit of β ) 0 the micelles differ only in the size that is determined by the number of molecules, for the case β > 0, the composition of β-sheets in the corona adds an additional degree of distinction (Figure 6b). Thus, two micelles composed from the exactly the same number of molecules can be different by the composition of β-sheets in the micellar corona. The hydrogen bonding also modifies the fraction of unimers. With an increase in β, free molecules aggregate due to the hydrogen bonding and form a polydisperse solution of β-sheets. The number-average aggregation number mn is considerably influenced by the fraction of low-molecular-weight aggregates and does not clearly reflect the effect of hydrogen bonding on the micellar fraction. Therefore, the polydispersity index varies over a wide range. Figure 6b shows a snapshot of the system for HH ) 2.5 and β ) 2.84 that corresponds to the peak of the polydispersity index. One can see a coexistence of micelles with β-sheets in the corona together with the polydisperse fraction of freely floating molecules and β-sheets. While hydrogen bonding leads to formation of long β-sheets (Figure 6c), the presence of hydrophobic tails on one end of the β-sheets significantly changes the local net hydrophobicity of the aggregates. This results in the formation of various amorphous aggregates from β-sheets to possibly jamming structures41 for relatively small HH e /HH. Especially interesting is the limit when HH g /HH. While strong hydrophobic interaction results in the formation of monodisperse micellar aggregates (Figure 6a), the hydrogen bonding leads to formation of β-sheets in the corona (Figure 6b). For small and intermediate values of the hydrogen-bonding energy, 1 < β < 3, the solution of such micelles is also monodisperse; the PDI is close to unity. On the other hand, with an increase in the hydrogen-bonding energy, β > 3, the polydispersity index approaches 2, suggest-

Figure 6. Snapshots from molecular simulations of peptide amphiphiles. (a) The spherical micelle, (b) the micelle with β-sheets on the outside forming the corona, (c) the β-sheets, and (d) the fiber aggregate. (e) Snapshot of the system that corresponds to the position of the peak of the polydispersity index (HH ) 2.5 and β ) 2.84). Monomer units are not shown full size for ease of viewing.

ing the change in the association mechanism from closed association to open association. We would like to note that PDI as a function of the hydrogen-bonding energy shows almost the same dependence for both HH ) 1 and 3. While in the first case (HH ) 1) we observe formation of β-sheets (Figure 6c), a number of long one-dimensional cylindrical fibers (Figure 6d) are found in the second case (HH ) 3). C. Internal Energy and Heat Capacity. Amphiphilic molecules segregate into micelles when the hydrophobic interaction energy is greater than the one corresponding to the critical micelle temperature. Micelle formation includes the dehydration of the hydrophobic tails (melting of the ordered shell of water molecules around the molecule) to form segregated domains which constitute the core of the micelle.42,43 In our simulation, we do not implicitly take into account the solvent molecules, and thus, only the intermolecular interaction between PA molecules contributes to the internal energy of the system. Figure 7a shows the internal energy normalized by HHNP as a function of hydrophobic interaction energy HH. In the limit of weak hydrogen bonding, β f 0, the normalized internal energy exhibits the jump from the zero value that characterizes the solution of free noninteracting molecules to a nonzero value. This nonzero value is defined by the equilibrium structure of the micelle and is equal to the average number of contacts in the core of the aggregate. The hydrophobic effect associated

Peptide Amphiphile Self-Assembly

J. Phys. Chem. B, Vol. 112, No. 8, 2008 2331 Formation of β-sheets in the micellar corona also leads to opening of the hydrophobic core due to the compact and dense structure of the β-sheets. When the linear size of the β-sheet overcomes the size of the micelle, the micelle becomes unstable. The growth of the heat capacity with an increase of the hydrogen-bonding energy, β, illustrates this effect; the peak in the heat capacity (Figure 7 b) determines the instability condition. For large values of the hydrogen-bonding energy, the micellar aggregate changes shape from spherical micelles to cylindrical fibers in which all β-sheets are aligned parallel to the cylinder axis. D. Static Structure Factor. In order to analyze the selfassembly of PA molecules at different length scales, we calculate the static structure factor S(q). The static structure factor describes the intensity of the incident beam scattered from the aggregates of size λ as a function of the scattering wave vector q ) 2π/λ. The static structure factor represents the resultant amplitude of scattering from the density heterogeneities. We calculate the static structure factor for all hydrophobic units

SH(q) ) Figure 7. Internal energy and heat capacity normalized by NP as a function of the pair attraction energy HH for different values of the hydrogen-bonding energy β.

with the micellization results in the release of latent heat.44 Figure 7b shows a broad peak in the heat capacity at /HH corresponding to the critical micelle temperature. In the limit of weak hydrogen bonding, β f 0, the peak position /HH is independent of the value of the hydrogen-bonding energy, and the nucleation is the dominant mechanism of aggregation. On the other hand, with further increase in the hydrogen-bonding energy, β, the peak position shifts in the direction of smaller /HH energies. At the same time, the normalized internal energy is decreasing with an increase in the hydrogen-bonding energy β due to the formation of hydrogen bonds and β-sheets. The thermodynamics of peptide amphiphile self-assembly is especially interesting due to the ability of these molecules to form a network of intermolecular hydrogen bonds. The hydrogen bonds can absorb heat into their internal degrees of freedom and thus can also store the energy. Figure 7b shows that the heat capacity is significantly increasing with the energy of hydrogen bonding β. The internal energy is clearly divided between two groups corresponding to the systems below and above the CMT. Whereas, in the limit of small hydrogenbonding energies, one can see almost no dependence on β, the internal energy is monotonically decreasing with further increases in the hydrogen-bonding energy. It is noticeable, that below the CMT, HH < /HH, the internal energy is decreasing with β much faster. These results indicate a resistance to hydrogen bonding when the micelles (hydrophobic aggregates) are formed. The reason for such behavior is that the hydrogen bonding requires a parallel arrangement of molecules, whereas the hydrophobic interaction results in the formation of spherical micelles with radial orientation of the molecules. Thus, to form a β-sheet among the peptide blocks in the corona, an additional force to rearrange the micelle and align at least the two nearest molecules must be applied. Formation of the rigid planar β-sheets on the micellar surface (Figure 6b) in analogy with the tiling of the surface of the sphere with finite size rods applies an additional geometric constraint.45 Increasing the hydrogen-bonding energy increases the average size of the β-sheet and the role of this geometric constraint.

1 N2pN2H

〈

|∑m,n ∑i,j eıqb‚∆br NP N H

mn ij

|〉 2

(9)

where ∆b r mn r mi - b r nj , b r mi ) (xmi ,ymi ,zmi ) is a Cartesian vector ij ) b that determines the position of the ith monomer unit of the mth amphiphilic molecule, and SH(q) is averaged among different directions. The wave vector is limited by two major length scales, the size of the simulation box, qmin ) 2π/L, and the size of the monomer unit, qmax ) 2π/σ. For many homogeneous systems, it is sufficient to know that for small wave vectors, q 2, the scattering mostly comes from the bond length that results in the peak around q = qmax (not shown on the Figure 8). In the narrow region for q > 1, we observe a typical Porod scattering46 when S(q) ∝ q-4, which is associated particularly with the scattering from the interface or the surface of aggregates. The scattering at intermediate values of the wave numbers, 0.1 < q < 1, when qR < 1, where R is the radius of gyration of the aggregates, is most informative. Figure 8a shows the static structure factor for the system above the CMT (HH ) 3) for different values of the hydrogen-bonding energy β. One can see the formation of a peak at finite q*. The peak comes from the interaggregates scattering, and the peak position corresponds to the inverse size of the aggregate. It is smeared out due to the large number of aggregates in the system. The wave vector, q*, is decreasing with β, indicating the growth of micellar aggregates in size. Figure 8b shows the dependence of the wave vector q* as a function of the hydrogen-bonding energy. In contrast to the very weak dependence on the hydrophobic energy, HH, the wave vector q* decreases with increasing hydrogen-bonding energy β and diminishes for some intermediate values β > /β. This corresponds to the crossover transition of finite size micellar aggregates to infinite aggregates. As mentioned earlier, in this region, β > /β and HH g /HH, we find long cylindrical fibers with β-sheets parallel to the fiber. IV. Conclusions Self-assembly is a reversible formation of large molecular aggregates. The equilibrium structure of peptide amphiphile aggregates is determined by the molecular structure and the

2332 J. Phys. Chem. B, Vol. 112, No. 8, 2008

Velichko et al.

Figure 9. Schematic phase diagram (β,HH) implied by results of simulation. The phase diagram includes regions with (a) free molecules, (b) spherical micelles, (c) micelles with β-sheets on the outside forming the corona, (d) long cylindrical fibers, (e) stacks of parallel β-sheets, (f) single β-sheets, and (g) the amorphous aggregate phase. Figure 8. (a) Double logarithmic plot of the spherically averaged static structure factor for different values of hydrogen-bonding energy β and HH ) 3. (b) The peak position, q*, as a function of hydrogen-bonding energy β.

balance of intermolecular forces. We study self-assembly of peptide amphiphile molecules by molecular simulation techniques with a model that incorporates molecular structure and inter- and intramolecular interactions. In particular, the model takes into consideration hydrophobic interactions and hydrogen bonding. We show that our model, despite its simplicity, captures various aspects of the peptide amphiphile’s selfassembly and allows the simulation to be done in a reasonable time frame. We show that “pure” hydrophobic interaction (the limit of β ) 0) results in the formation of monodisperse finite-size micelles (Figure 9a,b). Micellization follows a closed association scenario that involves a critical micelle temperature. The nucleation is the dominant mechanism of aggregation in this limit. On the other hand, “pure” hydrogen bonding (the limit of HH ) 0) leads to step-by-step aggregation of molecules into one-dimensional β-sheets. It follows an open association scenario that is characterized by a wide aggregate size distribution. We show that our model accurately describes both of these limits. The ability of the peptide blocks to form a network of hydrogen bonds significantly changes the phase behavior of peptide amphiphiles. At the same time the critical micelle temperature remains the most important parameters of the system. Below the CMT, hydrogen-bonding results in the formation of β-sheets with a broad aggregate size distribution (Figure 9f). The connection of amphiphilic molecules into β-sheets increases the number of hydrophobic tails and, as a result, increases the local concentration of hydrophobic units. This leads to the strengthening of the hydrophobic attraction between β-sheets and results in formation of rolls or stacks of β-sheets. Similar behavior was found in many different systems.47,48 With a further increase in the hydrogen-bonding

energy, long β-sheet aggregates form, resulting in various amorphous aggregates (Figure 9g). Above the CMT, hydrogen bonding changes the micellar structure by forming the β-sheets in the micellar corona. The steric interaction between β-sheets results in liquid crystalline ordering49 of β-sheets on the surface of the micelle. Formation of long, rigid, planar β-sheets in the spherical layer of the micellar corona applies additional geometric constraints. The curvature of the micellar surface induces geometric frustrations in the ordering of β-sheets.45 With further increase in the hydrogen-bonding energy, this results in instability and in breaking of the spherical symmetry of the micelle. This instability results in the formation of one-dimensional long cylindrical fibers where β-sheets are aligned in parallel along the fiber. Figure 9 summarizes our results into a compact schematic phase diagram, plotted in terms of the hydrogen-bonding energy, β, and the energy of interaction between hydrophobic units, HH. The phase diagram consists of regions of a homogeneous solution of free noninteracting molecules, single β-sheets connected laterally by hydrogen bonds, stacks of parallel β-sheets, spherical micelles, and micelles with β-sheets on the outside forming a corona, as well as long cylindrical fibers. A dimensionless packing parameter, S ) V/aolc, where V is a chain volume, ao is a head group area, and lc is a chain length, is convenient to predict the morphology of the aggregates.50 For many PAs, this parameter varies in the range 1/2 < S < 1, which corresponds to bilayers,51 whereas mostly all PAs form cylindrical nanofibers. This is not very surprising due to a complex dependence of the packing parameter on molecular length, as well as some thermodynamic parameters such as temperature and concentration.50,52,53 A variety of different aggregates found in our study suggest the possibility of achieving other morphologies. For example, a stack of β-sheets may grow along three orthogonal directions. In certain limits of concentrations and interaction energies, we

Peptide Amphiphile Self-Assembly

J. Phys. Chem. B, Vol. 112, No. 8, 2008 2333 shows the dependence of the normalized internal energy, Eint/ Np, and the weight-average aggregation number, mw, versus the quenching time τq for the case of β ) 4 and HH ) 4. One can see that both functions approach plateau for τq > 4 × 104, corresponding to structural relaxation of the system. On the other hand, due to the complexity of the system, we expect a strong finite size effect for large values of the interaction energies (HH > 3 and β > 3) when self-assembly results in the formation of long aggregates. Thus, we do not provide detailed statistical analysis of aggregates and limit our study to suggesting only trends in aggregate development. On the other hand, these results are very interesting and shed light on the self-assembly of peptide amphiphiles.

Figure 10. The weight-average aggregation number and the internal energy normalized by the number of molecules, E{int}/Np, as a function of quenching time, τq.

may find a bilayer structure in analogy with morphological transformations in block copolymers in selective solvent.54 We can identify a few sets of problems in need of further investigation. Under certain ionic conditions, the peptide amphiphile has a charged head group exposed to the surface of the aggregates. It was shown experimentally7 that the assembly of peptide amphiphiles into fibers in aqueous solution can be triggered by a change of ionic conditions (by pH or ionic strength changes). In our study, we do not consider electrostatic interactions between epitope head groups or with salt ions, which may induce different effects in fiber 55-59 and micelle formation.60 Second, the structure of the β-sheets is controlled by two dihedral angles. In our simulation, we choose their values to obtain planar β-sheets. A small change in the equilibrium value of one of these angles may induce the β-sheet helicity. The selfassembly of β-sheets with a helical twist may result in the formation of aggregates with a more complex morphology, in particular, nanofibers with a helical pitch.38,61 Third, the absence of solvent does not allow us to account for hydrogen bonding between peptide units and solvent molecules.28,43,62 However , we hope to extend our investigation with a systematic study of electrostatic effects in self-assembling peptide amphiphile systems; the analysis of solvent effects and PA-solvent hydrogen bonding requires significant computational resources. Acknowledgment. This work was supported by the MRSEC programs of the National Science Foundation (DMR-0520513 and DMR-0605427) at the Materials Research Center of Northwestern University, the National Institutes of Health (NCI 5-U54-CA119341, NIDCR 5-R01-DE015920 and NIBIB 5-R01EB003806), and the Department of Energy (DE-FG0200ER45810). Y.S.V. thanks Dr. Krista L. Niece, Megan Greenfield, and Dr. Liam C. Palmer for helpful and stimulating discussions. Appendix A: System Equilibration In preparing simulations to study the self-assembly of peptide amphiphiles, care must be taken to ensure that the lengths of the runs and the time step are properly chosen. The equilibration process is accompanied by a gradual decrease of temperature from Tmax ) 5 to 1. To ensure equilibrium conditions at lower temperatures, the proper quenching time τq must be found. To calculate the quenching time, we performed the set of simulation runs with quenching times varying in the range from 103 to 105 γσ2m/T. After reducing the temperature, another τq time steps were allowed to elapse before statistics were collected and perform measurements were once again conducted. Figure 10

References and Notes (1) Lehn, J. M. Science 2002, 295, 2400. (2) Ohgushi, H.; Miyake, J.; T. Tateishi, Tissue Engineering of Cartilage and Bone; Bock, G., Goode, J., Eds.; John Wiley and Sons: Chichester, U.K., 2003. (3) Cyranoski, D. Nature 2003, 424, 870. (4) Zubarev, E. R.; Pralle, M. U.; Li, L.; Stupp, S. I. Science 1999, 283, 523. (5) Zubarev, E. R.; Pralle, M. U.; Sone, E. D.; Stupp, S. I. J. Am. Chem. Soc. 2001, 123, 4105. (6) Hartgerink, J. D.; Beniash, E.; Stupp, S. I. Science 2001, 294, 1684. (7) Niece, K. L.; Hartgerink, J. D.; Donners, J. J. J. M.; Stupp, S. I. J. Am. Chem. Soc. 2003, 125, 7146. (8) Behanna, H. A.; Donners, J. J. J. M.; Gordon, A. C.; Stupp, S. I. J. Am. Chem. Soc. 2005, 127, 1193. (9) Rajangam, K.; Behanna, H. A.; Hui, M. J.; Han, X.; Hulvat, J. F.; Lomasney, J. W.; Stupp, S. I. Nano Lett. 2006, 6, 2086. (10) Palmer, L. C.; Velichko, Y. S. M. de la Cruz, M. O.; Stupp, S. I. Philos. Trans. R. Soc. London, Ser. B 2007, 365, 1417. (11) Binder, K. Monte Carlo and Molecular Dynamics Simulations in Polymer Science; Oxford University Press, Inc.: New York, 1995. (12) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford Science Publications: Oxford, U.K., 1987. (13) Frenkel, D.; Smit, B. Understanding Molecular Simulation. From Algorithms to Applications; Academic Press: London, 1996. (14) Marrink, S. J.; de Vries, A. H.; Mark, A. E. J. Phys. Chem. B 2004, 108, 750. (15) Leach, A. R. Molecular Modeling. Principles and Applications; Prentice Hall: New York, 2001. (16) Branden, C.; Tooze, J. Introduction to Protein Structure, 2nd ed.; Garland Publishing: New York, 1999. (17) Bratos, S.; Ratajcczak, H.; Viot, P. Hydrogen-Bonded Liquids; Dore, J. C., Texeira, J., Eds.; Kluwer Academic Publishers: The Netherlands, 1991. (18) Ladanyi, B. M.; Skaf, M. S. Annu. ReV. Phys. Chem. 1993, 44, 335. (19) Taylor, R.; Kennard, O. J. Am. Chem. Soc. 1982, 104, 5063. (20) Borgis, D.; Tarjus, G.; Azzouz, H. J. Chem. Phys. 1992, 97, 1390. (21) Luzar, A.; Chandler, D. Nature 1996, 379, 55. (22) Jedlovszky, P.; Turi, L. J. Phys. Chem. B 1997, 101, 5429. (23) Jedlovszky, P.; Vallauri, R. Mol. Phys. 1997, 92, 331. (24) Chertovich, A. V.; Ivanov, V. A.; Khokhlov, A. R.; Bohr, J. J. Phys.: Condens. Matter 2003, 15, 3013. (25) Nguyen, H. D.; Hall, C. K. Proc. Natl. Acad. Sci., U.S.A. 2004, 101, 16180. (26) Nyrkova, I. A.; Semenov, A. N. Eur. Phys. J. E 2005, 17, 327. (27) Nyrkova, I. A.; Semenov, A. N. Macromol. Theory Simul. 2005, 14, 569. (28) Cook, J. L.; Hunter, C. A.; Low, C. M.; Perez-Velasco, R. A.; Vinter, J. G. Angew. Chem. 2007, 46, 3706. (29) Pashuck, E.; Stupp, S. I. Private communication. (30) Chandler, D. J. Chem. Phys. 1978, 68, 2959. (31) Guarnieri, F. J. Math. Chem. 1995, 18, 25. (32) Nedea, S. V.; Frijns, A. J. H.; van Steenhoven, A. A.; Markvoort, A. J.; Hilbers, P. A. J. Phys. ReV. E 2005, 72, 016705. (33) Sˇ olc, K.; Elias, H.-G. J. Polym. Sci. 1973, 11, 137. (34) Huyskens, P.; VanBeylen, M.; Verheyden, H. Pure Appl. Chem. 1996, 68, 1535. (35) Leibler, L.; Orland, H.; Wheeler, J. C. J. Chem. Phys. 1983, 79, 3550. (36) Kao, C. H. R.; de la Cruz, M. O. J. Chem. Phys. 1990, 93, 8284. (37) Shusharina, N. P.; Saphonov, M. V.; Nyrkova, I. A.; Khalatur, P. G.; Khokhlov, A. R. Phys. Chem. Chem. Phys. 1996, 100, 857.

2334 J. Phys. Chem. B, Vol. 112, No. 8, 2008 (38) Berl, V.; Schmutz, M.; Krische, M. J.; Khoury, R. G.; Lehn, J.-M. Chem.sEur. J. 2002, 8, 1227. (39) Astafieva, I.; Zhong, X. F.; Eisenberg, A. Macromolecules 1993, 26, 7339. (40) Odian, G. Principles of Polymerization, 2nd ed.; Wiley Intersciences: New York, 1981. (41) Loscar, E. S.; Borzi, R. A.; Albano, E. V. Phys. ReV. E 2006, 74, 051601. (42) Lumry, R.; Rajender, S. Biopolymers 1970, 9, 1125. (43) Bukk, A.; Høye, J. S.; Hansen, A. Biophys. J. 2001, 81, 710. (44) Chen, L.-J.; Sheu, Y. H.; Li, P. J. J. Phys. Chem. B 2004, 108, 19096. (45) (a) Nelson, D. R. Nano Lett. 2002, 2, 1125. (b) Vitelli, V.; Nelson, D. R. Phys. ReV. E 2006, 74, 021711. (46) Bale, H. D.; Schmidt, P. W. Phys. ReV. Lett. 1984, 53, 596. (47) Nyrkova, I. A.; Semenov, A. N.; Aggeli, A.; Boden, N. Eur. Phys. J. B 2000, 17, 481. (48) Aggeli, A.; Nyrkova, I. A. M.; Bell, R.; Harding, L.; Carrick, T.; McLeish, C. B.; Semenov, A. N.; Boden, N. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 11857. (49) de Gennes, P. G.; Prost, J. The Physics of Liquid Crystals; Clarendon: Oxford, U.K., 1993. (50) Israelachvili, J. N. Intermolecular and Surface Forces; Academic Press: London, 1992.

Velichko et al. (51) Yaliraki, S. N.; Longo, G.; Gale, E.; Szleifer, I.; Ratner, M. A. J. Chem. Phys. 2006, 125, 074708-1 (52) Szleifer, I.; Kramer, D.; Benshaul, A.; Gelbart, W. M.; Safran, S. A. J. Chem. Phys. 1990, 92, 6800. (53) Kumar, V. V. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 444. (54) Lodge, T. P.; Pudil, B. Hanley, K. J. Macromolecules 2002, 35, 707. (55) Tsonchev, S.; Schatz, G. C.; Ratner, M. A. Nano Lett. 2003, 3, 623. (56) Solis, F. J.; Stupp, S. I.; de la Cruz, M. O. J. Chem. Phys. 2005, 122, 054905. (57) Velichko, Y. S.; de la Cruz, M. O. Phys. ReV. E 2005, 72, 041920. (58) Velichko, Y. S.; de la Cruz, M. O. J. Chem. Phys. 2006, 124, 214705. (59) Kohlstedt, K. L.; Solis, F. J.; Vernizzi, G.; de la Cruz, M. O. Phys. ReV. Lett. 2007, 99, 030602. (60) (a) Huang, C.; de la Cruz, M. O.; Delsanti, M.; Guenoun, P. Macromolecules 1997, 30, 8019. (b) Huang, C.; de la Cruz, M. O. Macromolecules 2002, 35, 976. (61) Li, L.-S.; Jiang, H.; Messmore, B. W.; Bull, S. R.; Stupp, S. I. Private communication. (62) Sum, A. K.; Faller, R.; de Pablo, J. J. Biophys. J. 2003, 85, 2830.