Mechanism of the pH-Controlled Self-Assembly of Nanofibers from Peptide Amphiphiles.

Stimuli-responsive, self-assembling nanomaterials hold a great promise to revolutionize medicine and technology. However, current discovery is slow and often serendipitous. Here we report a multiscale modeling study to elucidate the pH-controlled self-assembly of nanofibers from the peptide amphiphiles, palmitoyl-I-A3E4-NH2. The coarse-grained simulations revealed the formation of random-coil based spherical micelles at strong electrostatic repulsion. However, at weak or no electrostatic repulsion, the micelles merge into a nanofiber driven by the β-sheet formation between the peptide segments. The all-atom constant pH molecular dynamics revealed a cooperative transition between random coil and β-sheet in the pH range 6-7, matching the CD data. Interestingly, although the bulk pKa is more than one unit below the transition pH, consistent with the titration data, the highest pKa's coincide with the transition pH, suggesting that the latter may be tuned by modulating the pKa's of a few solvent-buried Glu side chains. Together, these data offer, to our best knowledge, the first multiresolution and quantitative view of the pH-dependent self-assembly of nanofibers. The novel protocols and insights gained are expected to advance the computer-aided design and discovery of pH-responsive nanomaterials.

Introduction

Peptide amphiphiles (PAs) are peptides with a terminal group covalently linked to a long hydrophobic chain. Over the past decade, biocompatible PAs have attracted much attention as building blocks of nanomaterials due to the ability to self-assemble into supramolecular structures at specific solution conditions relevant for biomedical and biotechnological applications (see a recent review[1] and references therein). Using TEM, FTIR, and other techniques, Stupp and co-workers showed that PAs self-assemble into nanofibers[2] upon a change in the external stimuli such as solution pH and ionic concentrations and that the morphology, surface chemistry, and potential bioactivity of the nanomaterials can be controlled by the selection of amino acid sequence or the modification of the alkyl tail.[3] This pioneering work set the stage for designing PA-based nanomaterials with desired properties to target specific applications.[4]

Of particular interest to us are the PAs that can form nanocarriers for delivery of therapeutic or imaging agents in response to endogenous stimuli such as a drop in pH at tumor site.[5] Toward this end, Goldberger and co-workers designed a series of PAs that can self-assemble into nanofibers when the solution pH is decreased from the normal physiological condition 7.4 to 6.6, the mild acidic condition found in the microenvironment of malignant tumor tissues.[6] The PAs used are the charged polypeptides palmitoylated at the N-terminus, CH3-(CH2)14CO-NH-X-Ala3-Glu4-CO-NH2, where X is a hydrophobic residue, such as Ile, Phe, Val, or Tyr, with a varying propensity for β-sheet formation. The Glu residues are present to allow pH response. CD and TEM experiments revealed that, at high pH, PAs at low concentrations are in random-coil states while PAs at high concentrations form micelles. Goldberger and co-workers also conducted titration experiments, which showed that the bulk pKa of the nanofibers remains at 4.7, while the self-assembly transition pH decreases from 6.6 to 6.0, when Ile is substituted by Tyr. Thus, they concluded that the transition pH is not dependent on the pKa of the Glu side chains but rather on the β propensity of the X residue. The transition pH is shifted lower because Tyr has a lower beta propensity than Ile.[7] Below we will refer to the sequence of the PA with X being Ile as PA1.

In order to systematically design PAs that can self-assemble into nanomaterials at specified pH conditions with desired properties, it is desirable to learn the detailed mechanism of the pH-induced self-assembly process and to make quantitative prediction of the transition pH. Using the standard all-atom molecular dynamics (MD) as well as various coarse-grained Monte Carlo and MD simulations, studies led by Stupp, Ratner, and de la Cruz groups have provided valuable insights into the structural and kinetic aspects of the PA self-assembly.[8−11] However, there are two limitations of these studies. First, the secondary structure change that accompanies the nanofiber formation cannot be described by the coarse-grained models used. And importantly, the pH-dependent mechanism of the self-assembly cannot be studied by the standard all-atom MD simulations.

Here we present a multiscale MD study to elucidate the pH-dependent kinetic and thermodynamic aspects of the nanofiber self-assembly from the aforementioned PA synthesized in the Goldberger lab.[6] The kinetic details will be qualitatively mapped out using discrete molecular dynamics with the Protein Intermediate Resolution model (PRIME)[12−14] that was originally developed by the Hall group and recently extended by Nguyen and co-workers (ePRIME).[15] The mechanism of the pH-dependence will be quantitatively uncovered by studying the unfolding of a tetrameric β-sheet model using the recently developed all-atom version[16] of the continuous constant pH molecular dynamics (CpHMD) technique[17,18] with the pH-based replica-exchange protocol.[19] The CG simulations showed that the PAs self-assemble into micelles, which merge into an elongated nanofiber driven by the β-sheet formation at low pH conditions. The all-atom CpHMD simulations revealed a cooperative transition from β-sheet to random coil and were able to determine the bulk pKa as well as individual pKa’s of Glu side chains, both of which are in quantitative agreement with the experiment. Surprisingly, the data showed that the transition pH is determined by the deprotonation of the central Glu residues. The combined coarse-grained and all-atom simulation data offer, to our best knowledge, the first multiresolution and quantitative view of the pH-dependent mechanism of PA self-assembly.

Simulation Methods

Coarse-Grained Discontinuous Molecular Dynamics Simulations

The self-assembly process was simulated using the discontinuous molecular dynamics[20] with our newly developed coarse-grained model ePRIME for the representation of the PA molecule, CH3-(CH2)14CO-Ile-Ala3-Glu4-CO-NH2.[15] ePRIME is an extension of the original PRIME model, in which each amino acid residue is represented by four beads, three of which represent the backbone while a single bead represents the side chain.[12−14] The ePRIME model offers representations for all 20 amino acid side chains by taking into account size, hydrophobicity, and charge using well-established experimental data and modeling parameters.[21,22] The discontinuous molecular dynamics (MD) is orders of magnitude faster than regular molecular dynamics because the forces are modeled by hard-sphere or square-well potentials. Solvent effect is included implicitly as a potential of mean force in the energy function. Hydrogen bonding between the backbone amide and carbonyl groups is represented by a directionally square-well attraction potential of the magnitude εHB, which serves as a reference strength for other types of interactions.

The self-assembly simulations were initiated from a configuration in which 800 PA molecules were randomly placed in a cubic box with periodic boundary dimensions of 250 Å × 250 Å × 250 Å. The system contained 42 400 particles representing a PA concentration of 85 mM. Each simulation was heated at high temperature (reduce temperature of 5.0) to reach an unbiased initial configuration and then quickly cooled to the temperature of interest (reduced temperature of 0.08–0.15), for a constant-temperature production run until equilibration. To mimic the effect of pH, each simulation was conducted at an electrostatic repulsion strength of 0%, 200%, 400%, and 600%, relative to that of the hydrogen bond, εE/εHB. Each simulation was independently repeated ten times to ensure statistical significance. Quantitative results in this paper are the averages of the last 10% of simulation data with error bars taken from the standard deviation. The secondary structures were defined through the implementation of STRIDE.[23] Here we focused on α-helix, β-strand or β-sheet, and random coil that also contains turn. An aggregate is defined if each peptide in a group of PA molecules has at least two interpeptide hydrogen bonds or four hydrophobic interactions with a neighboring PA molecule in the same group.

All-Atom Constant pH Molecular Dynamics Simulations

The initial structure of the tetrameric β-sheet of PA1 was built with PyMOL[24] in which the four PA1 peptides were placed in a well-defined parallel β-sheet conformation with protonated Glu side chains. The β-sheet was immersed in a cubic box of 4296 water molecules, keeping a minimal distance of 9 Å between the solute and each face of the box. The system was then subjected to energy minimization prior to a 100 ns standard molecular dynamics (MD) run with fixed protonation states using the GROMACS package.[25] The CHARMM C36 force field was used to represent the proteins,[26] and the modified TIP3P model was used for water.[27] The MD run was performed with periodic boundary conditions at constant temperature, pressure, and fixed protonation states. The temperature was maintained at 298 K using the Nosé–Hoover thermostat,[28,29] while the pressure was maintained at 1 atm using the Parrinello–Rahman barostat.[30] The LINCS algorithm[31] was used to constrain all bonds containing hydrogen atoms, allowing an integration time step of 2 fs. The coordinates were saved every 2 ps. The particle-mesh-Ewald method[32,33] was used for the calculation of electrostatic interactions, and a cutoff of 12 Å was used for the calculation of van der Waals interactions. If not specified, the default setting was used for all other input parameters.

Final snapshots of the PA1 tetramer were extracted from the above simulations to initiate the all-atom, pH-based replica-exchange (REX)[16] continuous constant pH MD (CpHMD)[17,18] using an in-house modified version of the CHARMM package (c37b1).[34] The same force fields as in the GROMACS simulations were used. The tetramer was resolvated in a cubic box filled with 5013 TIP3P water molecules. Sixteen co-ions were added to level the total charge,[16] while 16 sodium ions were added for neutralizing the system. No additional salt ions were added to avoid the possible convergence issue related to the insufficient sampling of salt ions. In the pH-based REX protocol,[19] 24 replicas were placed in the range of pH 3–11 with an interval of 0.5 pH units. The exchange of pH was attempted every 250 steps or 0.5 ps between replicas adjacent in the pH ladder. Additional replicas were added in order to achieve an exchange ratio of at least 15% at all pH conditions. Each pH replica underwent 12 ns molecular dynamics at constant temperature, pressure, and pH. The temperature was 298 K using the Hoover thermostat,[29] while the pressure was maintained at 1 atm using the Langevin piston pressure coupling algorithm.[35] The SHAKE algorithm was applied to all the hydrogen bonds and angles to allowed a time step of 2 fs, and the coordinates were saved every 0.5 ps. The electrostatic interactions were calculated using the generalized reaction field method with a cutoff of 14 Å. The same cutoff distance was used in the calculation of van der Waals interactions. The REX all-atom CpHMD was repeated once by starting from a different set of initial velocities.

Following our previous work,[18] the protonated and unprotonated forms of a titratable side chain are defined as the states with λ value below 0.1 and above 0.9, respectively. After each replica-exchange attempt, the λ values of all Glu’s were recorded. The unprotonated fractions were calculated for individual Glu’s as well as by taking into account all 16 Glu’s. Fitting the former data at different the pH values to the generalized Henderson–Hasselbalch equation gives the microscopic pKa’s, while fitting the latter data gives the macroscopic or bulk pKa of the β-sheet. The data from the last 2 ns (per replica) was used for analysis and pKa calculations.

Results and Discussion

Kinetic Mechanism of the Self-Assembly Process

Discrete molecular dynamics with the ePRIME model[15] was applied to study the self-assembly process of 800 PA1 molecules (CH3-(CH2)14CO-NH-I-A3-E4-CO-NH2) that were randomly distributed in the simulation box at the beginning. To mimic the varying degree of ionization of the EEEE-stretch under different pH conditions, the strength of the electrostatic repulsion between the peptide segments was varied, relative to that of a hydrogen bond, i.e., εES/εHB. Driven by hydrophobic forces, small spherical micelles were quickly formed, burying the alkyl groups in the interior and exposing the peptide segments on the surface (Figure 1). However, the next sequence of kinetic steps was determined by the strength of the electrostatic repulsion between the EEEE-stretches. When the electrostatic repulsion is strong, i.e., εE/εHB ≥ 200%, representing the high pH conditions, under which some or all of the Glu side chains were deprotonated or charged, the micelles remain in the spherical arrangement with the peptide segments uniformly distributed along the surface to minimize repulsion between neighboring EEEE-stretches (Figure 1A).

Figure 1

Snapshots of the self-assembly process of PA1 molecules as a function of CG simulation time (in reduced time units). (A) A spherical micelle without β-sheet content is formed at a strong electrostatic repulsion (200% εHB). (B) A cylindrical nanofiber with β-sheet content is formed at zero electrostatic repulsion. PA1 molecules are colored as follows: alkyl tail in red, isoleucine in green, alanine in blue, and glutamic acid in pink.

By contrast, in the absence of electrostatic repulsion, representing the low pH conditions, under which the Glu side chains are fully protonated or neutral, hydrogen bonds started to form between the backbone amide and carbonyl groups of the adjacent peptide segments, leading to the formation of β-sheets (Figure 1B). As a result, the peptide segments on the surface of the micelle were no longer uniformly distributed, creating segregation that facilitates the exposure of hydrophobic alkyl groups to the surface. These micelles subsequently merged at the exposed locations through hydrophobic interactions between the alkyl groups. Thus, the early stage of the nanofiber self-assembly was characterized by the sharp increase in the β-sheet content, which accompanies the highly frequent events of micelle merging. Therefore, a kinetic mechanism for the formation of nanofibers is deciphered, involving the following steps: (1) simultaneous formation of multiple spherical micelles driven by hydrophobic interactions between alkyl tails and (2) micelle-merging events prompted by the exposure of the hydrophobic core due to the change in the architecture of the spherical micelles that ultimately leads to the fabrication of one continuous rod-shaped nanostructure. The latter supports the notion that β-sheet promotes fiber formation based on the existing experimental data.[36]

Structural Details of the Self-Assembled Nanofiber in Agreement with Experiment

The in silico observation that PA1 molecules self-assemble into spherical micelles and a rod-shaped nanofiber at conditions of high and low electrostatic repulsion, respectively, is in agreement with the CD and TEM data of Goldberger lab, which showed the formation of spherical micelles at high pH and nanofibers at low pH when the PA concentrations were high (corresponding to the simulation condition).[6] The in silico observation of β-sheet formation in the nanofiber and not in the micelles also agrees with the aforementioned CD data[6] as well as the earlier studies of Stupp and co-workers using FTIR[2] and CD experiments.[4] The rod-shaped nanofiber with exposed peptide segments and buried alkyl tails observed in the simulation is consistent with the model proposed by Stupp and co-workers over a decade ago based on cryo-TEM imaging and staining, which showed PA molecules forming cylindrical nanofiber with the charged amino acids on the surface and alkyl tails buried in the core.[2]

Electrostatic Repulsion Controls the Micelle-to-Nanofiber Transition

Since the intermolecular electrostatic repulsion between the charged Glu side chains favors the formation of micelle, while the backbone hydrogen bonding between peptide segments drives the distortion and subsequent merging of spherical micelles leading to the formation of nanofiber, weakening electrostatic repulsion would shift the equilibrium from micelle to nanofiber. The latter corresponds to a decrease in the degree of ionization of Glu side chains and therefore mimics the condition of lowering pH. In fact, the coarse-grained simulations resulted in a nanofiber when the electrostatic repulsion strength is below 200% relative to the hydrogen bonding (εES/εHB) (Figure 2A). In the absence of electrostatic repulsion (all Glu side chains are neutral), the self-assembled nanofiber contains a β-sheet content of almost 35% (Figure 2A). As the relative electrostatic repulsion strength is at or above 200%, β-sheet content drops below 5% (Figure 2A). The typical size of a β-sheet in the cylindrical nanofiber fluctuates between two and four (Figure 2B). The range of β-sheet content observed in our simulation is in good agreement with the CD data by Stupp and co-workers, who showed the β-sheet population of IKVAV-bearing PAs to be 25 ± 20%.[4] The electrostatic repulsion strength at the micelle-to-fiber transition for PA1 is lower than but consistent with that for the palmitoyl-V3A3E3 molecules (εE/εHB ≤ 300–400%) from our previous work,[15] likely due to the fact that the latter sequence has one fewer charged residue and two extra hydrophobic residues as compared to PA1.

Figure 2

A. β-sheet content at varying strength of electrostatic repulsion in the CG simulation. B. Snapshot of a typical β-sheet in the CG simulation.

Establishing a Tetrameric β-Sheet Model for the All-Atom CpHMD Simulations

The coarse-grained simulations have established the pH-induced β-sheet formation as the hallmark and driving force for the creation of nanofibers. To understand the details of the pH-dependent transition between the random coil (as present in micelle) and β-sheet (nanofiber), we employed the all-atom continuous constant pH molecular dynamics (CpHMD) simulations with the pH replica-exchange sampling protocol. To make the study amenable to the current computational capability, we focused on the pH-dependent stability of a minimal β-sheet model comprising four PA1 peptides (NH2-IAAAEEEE-CO-NH2), consistent with the observed sheet size in the CG simulation (Figure 2). The alkyl tails were omitted as they do not contain ionizable residues and as such do not influence the pH dependence. The tetramer was prebuilt with four parallel β-strands, and all Glu were protonated to represent the least electrostatic repulsion, which is the most favorable condition for the formation of the β-sheet based nanofiber. The tetramer was then subjected to the standard all-atom MD simulations with fixed protonation states to verify stability and establish the starting structure as well as control for the subsequent pH-dependent simulations. Indeed, within 100 ns, the PA1 tetramer remained stable with a well-defined β-sheet along the Glu sequences. A snapshot for the tetramer is shown in Figure 3.

Figure 3

Snapshot of the PA1 tetramer obtained from the 100 ns standard MD simulation with all Glu side chains fully protonated. Two views (related by a 180° rotation) are presented. The Glu side chains are highlighted in blue stick representation. The shade of blue (white to blue) corresponds to the pKa value (0 to 8). The pKa values are extracted from an explicit-solvent CpHMD simulation (see Figure 5).

Figure 5

Titration of Glu residues in the PA1 tetrameric β-sheet. (a) Calculated pKa values of the individual Glu side chains. The black dashed line indicates an average pKa of 5.4. (b) The unprotonated fraction averaged over all the Glu’s of the PA1 tetramer as a function of pH. The curve represents the fitting to the generalized Henderson–Hasselbalch equation. The obtained bulk pKa is 5.4 with a Hill coefficient of 0.6. The error bars correspond to the standard deviation among the unprotonated fractions of all the Glu’s.

Tetrameric β-Sheet that Undergoes a pH-Dependent Unfolding Transition

We next performed the all-atom replica-exchange (REX) CpHMD simulations[16] on the PA1 tetramer starting from the structures extracted from the aforementioned standard MD simulations. Two independent runs of REX CpHMD simulations were performed. Each run of simulations included pH replicas in the pH range 3–10. Simulation of each replica lasted 12 ns per replica. More details are given in Simulation Methods. The discussion below will focus on simulation run 1. Figure 4a shows the residue-based β-sheet content. At pH 3, the β-sheet is most stable and spans residues 2–7 (AAAEEE) in the monomer. As pH is increased to 6.5, the β-sheet is shortened to residues 4–6 (AEEE). As pH is increased further to 8, the β content reduces to below 20% in those regions, indicating that the tetramer is nearly unfolded. It is worthwhile to notice that the β-sheet is less stable in the two edge monomers as compared to the central ones, as a result of solvent exposure.

Figure 4

pH-dependent β-sheet content of the PA1 tetramer. (a) Residue-based β-sheet content calculated as the percentage residence time of a residue in the β-sheet conformation (assigned using the DSSP algorithm[37]). Blue, green, and red represent pH 3, 6.5, and 8, respectively. The four monomers correspond to residues 1–8, 9–16, 17–24, and 25–32, respectively. (b) Total number of backbone hydrogen bonds as a function of pH. Data points are the averages while the error bars represent the standard deviation. Data from the last 2 ns of the all-atom CpHMD simulation run 2 were used.

To quantitatively characterize the pH-dependent unfolding transition, we calculated the total number of backbone hydrogen bonds for each snapshot (Figure 4b). At low pH conditions, pH 3–6, there are on average 11 backbone hydrogen bonds, corresponding to a β-sheet content of 26%, in good agreement with the coarse-grained data (35%, see Figure 2). At high pH conditions, pH 7.5–11, the average number of hydrogen bonds is 4.5, in agreement with the β-content found for individual residues (Figure 4a). Interestingly, a sharp drop in the number of hydrogen bonds occurs in a narrow pH range of 6–7, indicating that the β-sheet starts to unfold at pH 6 and the unfolding is complete at pH 7. The latter is in a quantitative agreement with the CD data, indicating that at 10 μM the PA1 molecules undergo a random-coil to β-sheet transition in the pH range pH 6.8–6.0 with an estimated transition pH of 6.6.[6] To verify convergence, a second set of replica-exchange CpHMD simulations was performed starting from a different velocity seed. The resulting transition pH region remains the same.

pKa Values for the Glu’s in the Central Strands Determine the Transition pH

To further understand the pH-dependent β-sheet formation, we examine the calculated pKa values of individual Glu side chains in the PA1 tetramer and their degree of ionization (also known as unprotonated fractions) at different pH values. The individual pKa’s range from 4.2 to 6.6 with most of them shifted up from the model value of 4.4, giving an average pKa of 5.4 (Figure 5a), which is 0.7 units higher than the measured average pKa of the PA1 solution based on acid-based titration.[6] We note the second set of simulations gave an average pKa of 5.2. The pKa up-shifts are due to the intermolecular electrostatic repulsion between the ionized Glu’s on the adjacent β-strands. The pKa shifts are more pronounced for Glu’s in the central strands (residues 13–16 and 21–24) as compared to those in the edge strands (residues 5–8 and 29–32), because the central Glu’s are shielded from solvent, as evident from the significantly smaller solvent accessible surface area of the carboxylic groups in the protonated form as compared to those of the edge residues (data not shown). Thus, the wide spread in the individual pKa’s is a result of different local environment of Glu side chains, which can be visualized by projecting the pKa values onto the initial structure of the PA1 tetramer (Figure 3). It can be seen that the lowest pKa values are from those Glu side chains in the edge strands pointing out to solvent, while the highest pKa’s are from those Glu side chains in the two central strands that are in a confined configuration (Figure 3). Interestingly, the Glu side chains in the edge strands, which have a relatively higher pKa value (see n = 7 and 31), point inside the tetramer and interact with the Glu’s in the central strands (Figure 3).

Structure analysis revealed that, at low pH conditions, the carboxylic groups of the Glu’s in the central β-strands form hydrogen bonds, which disappear at high pH conditions (data not shown). These hydrogen bonds stabilize the β-sheet as well as the protonated state. Thus, in addition to intermolecular charge–charge repulsion, hydrogen bonding is another contributor to the pKa upshift for the central Glu’s. Interestingly, intermolecular hydrogen bonding between carboxylic side chains and the effect on the pKa shift have also been suggested for fatty acids in the premicellar aggregates based on experimental data[38] and recently confirmed by simulation.[39]

Related to the large variation in the pKa’s is the varying degree of ionization of the individual Glu’s at a specific pH (Figure 5b). The variation is particularly large in the pH region where the unfolding of the β-sheet occurs. For example, at pH 6, the degree of ionization for individual Glu’s can be as low as 30% and as high as 100% (error bars in Figure 5b). Similar effects of interpeptide electrostatic repulsion leading to a wide range of pKa values have been suggested for the β-sheet self-assembly of small peptides containing Glu and Arg residues.[40] It is reasonable to hypothesize that the tetramer β-sheet unfolds when all the Glu side chains become ionized. Indeed, above pH 7, when the unfolding transition is complete (Figure 4), the Glu’s are at least 95% ionized including those in the two central strands (Figure 5b). By contrast, at pH 6, when the β-sheet starts to unfold, the Glu’s in the two edge strands are completed ionized, while those in the two central strands become ionized by at least 30%. The average ionization is about 65%. Thus, these data suggest that, to completely unfold the β-sheet, all Glu’s need to be ionized. Since Glu’s of the central strands ionize at much higher pH than those of the edge strands, their pKa’s determine the transition pH. The average degree of ionization (unprotonated fraction) over all the Glu’s as a function of pH can be fit to the Hill equation, giving an average or bulk pKa of 5.4 (Figure 5b), in agreement with the one obtained by averaging the individual pKa’s (Figure 5a). The Hill coefficient is 0.6 (Figure 5b), indicating a high degree of anticooperativity due to electrostatic repulsion between the charged Glu side chains.

Conclusions

In this study we employed the state-of-the-art coarse-grained (CG) and all-atom molecular dynamics simulations to elucidate the kinetic and thermodynamic mechanisms of the pH-dependent self-assembly of PA1 molecules. The CG simulations revealed that the PA1 molecules initially aggregate to form multiple spherical micelles with the palmitoyl tails buried in the micellar core and the PA1 peptides in the random-coil state exposed to solvent. When the electrostatic repulsion between the peptide segments is weak, corresponding to the low pH conditions when the Glu side chains are (mainly) neutral, interpeptide hydrogen bonds begin to form, which perturbs the spherical shape of the micelle, leading to the exposure of hydrophobic regions. The latter drives the micelles to merge, ultimately leading to the formation of one continuous rod-shaped nanofiber with the PA1 peptides aligned in a β-sheet configuration. The kinetic pathway to nanofiber driven by the β-sheet formation is in agreement with existing experimental work.[2,4,6]

To investigate the pH-dependent thermodynamic aspect, conformational dynamics of a PA1 tetrameric β-sheet model was studied using the all-atom CpHMD simulations. On the basis of the number of backbone hydrogen bonds, the unfolding of the β-sheet occurs between pH 6 and pH 7 with the transition midpoint at about pH 6.5. Remarkably, the individual pKa’s of the Glu’s in the tetramer vary widely, and the calculated average or bulk pKa is 5.3 (5.4 and 5.2 in two independent runs). The Glu’s of the edge strands are generally solvent exposed and have pKa’s close to the model value of 4.4. By contrast, the Glu’s of the central strands, especially those three next to Ala, are shielded from solvent and have pKa’s that are shifted to 6–6.6. Together, these data suggest that ionization of the three Glu’s next to Ala in the central peptides induces unfolding of the β-sheet. Their pKa’s determine the pH condition for the transition from nanofiber to micelle. The cooperative unfolding transition of PA1 β-sheet is reminiscent of the pH-dependent phase behavior of self-assembling peptides containing Glu’s, for which the ionization of single Glu was suggested to be responsible.[40]

Finally, we remark on the caveats of the simulations. In the coarse-grained simulations, the PA concentration (85 mM) was more than 3 orders of magnitude higher than the experimental condition (10 μM). Therefore, alternative kinetic pathways such as the direct formation of nanofiber from single PA molecules can not be studied as it would require orders of magnitude longer simulation time which is currently not feasible. In the all-atom simulations, since only counterions and no additional salt ions were included (to avoid the potential convergence problem due to limited sampling of ions), the ionic strength was lower than that in experiment (150 mM).[6] We suggest that this is a major cause for the calculated average pKa to be about 0.6 units higher than the measured value, as lower ionic strength leads to the overestimation of electrostatic repulsion which increases the pKa’s. Another limitation of the all-atom simulation is that the concentration effect is not accounted for due to the nature of the unfolding simulation. However, despite these limitations, the data show that the combined coarse-grained and all-atom constant pH molecular dynamics simulations present an attractive strategy that may be used to guide the design and discovery of pH-responsive nanomaterials.