Water ordering controls the dynamic equilibrium of micelle–fibre formation in self-assembly of peptide amphiphiles

Understanding the role of water in governing the kinetics of the self-assembly processes of amphiphilic peptides remains elusive. Here, we use a multistage atomistic-coarse-grained approach, complemented by circular dichroism/infrared spectroscopy and dynamic light scattering experiments to highlight the dual nature of water in driving the self-assembly of peptide amphiphiles (PAs). We show computationally that water cage formation and breakage near the hydrophobic groups control the fusion dynamics and aggregation of PAs in the micellar stage. Simulations also suggest that enhanced structural ordering of vicinal water near the hydrophilic amino acids shifts the equilibrium towards the fibre phase and stimulates structure and order during the PA assembly into nanofibres. Experiments validate our simulation findings; the measured infrared O–H bond stretching frequency is reminiscent of an ice-like bond which suggests that the solvated water becomes increasingly ordered with time in the assembled peptide network, thus shedding light on the role of water in a self-assembly process.

Thermogravimetric analysis from 25 -200 o C of c16-AHL 3 K 3 -CO 2 H dried film from water (blue) and c16-AHL 3 K 3 -CO 2 H dried film from fibers, 100 mM NH 4 OH (red), leftaxis. Percentage of water trapped in the fibrous assembly determined by plotting the percent difference in mass of the two films (green), right -axis.  The water molecules near the N1, N4, O5, O9, and N11 atoms of Alanine, Histidine, Leucine, and Lysine amino acids, respectively orient themselves to form hydrogen bond with these polar atoms as shown in Supplementary Figure 1. Strong solvation around the hydrophobic region of the PAs hydrophobic tail can also be seen during initial ~150 ps of simulations. This solvation of water molecules near the hydrophobic tail would lead to formation of water-cage-like structure around these hydrophobic groups due to lack of ability of these non-polar hydrophobic groups to form hydrogen bonds with water.

Order Parameters
Supplementary Table 1 and Supplementary Table 2 show the translational and tetrahedral order parameters near hydrophilic and hydrophobic groups of fibers, respectively, during the initial 200 ps of stage 1. We find that the translational and tetrahedral order of water increases and decreases, respectively, with increase in simulation time and approach the values of bulk water. Note, the bulk water has a translational order of ~1.2 and tetrahedral order parameter of ~0.68. 3 This suggests that water molecules near hydrophobic groups of PAs are ordered in very initial stages of simulations (< 150 ps). On the other hand water molecules near hydrophilic groups the translational and tetrahedral order of water increases and decreases, respectively. This suggests increase in ordering of water molecules near hydrophilic groups with increase in simulation time.

Potential of Mean Force of PA Aggregation
The entropic contribution, −T ΔS(r) to the PMF, ΔW(r), at 280 K is obtained through supplementary equation 5. The entropic (−T ΔS(r)) and enthalpic (ΔH(r)) contributions to the PMF at 280 K are shown in Supplementary Figure 4 along with the PMF. The entropic and enthalpic contributions of the free energy act in opposite direction to one another, and the relative proportion of the two contributions depend on the distance between the two PA molecules. Similar observations has been reported in prior studies both at nano and macro-scales. 5,6 As can be seen in Supplementary

Residence Time Probability of Water
Supplementary Figure 6 shows the P res (t) for water molecules in hydrophilic and hydrophobic region of PA molecules. The water molecules near hydrophilic region reside for longer time as compared to water molecules near hydrophobic region. This shows that structure of water molecules near hydrophilic groups is more stable as compared to the water molecules near hydrophobic region. This can be attributed to the presence of water molecules in the PA network that form stronger hydrogen bonds with the amino acids of PA, assisting water molecules to form stable network near hydrophilic region.

Hydrogen Bond Characteristics
In the systems studied, the hydrogen bonds were classified into three types based on donor-acceptor pair: (1)  In the case of type (1) hydrogen bonds, both the donor and acceptor groups are on PA.
The autocorrelation function decays slowly as compare to type (2) and type (3) hydrogen bonds. The slower decay observed in the case of type (1) hydrogen bonds can be attributed to the self-assembly of PAs which leads to strong hydrogen bonding among themselvels. Additionally, β-sheets formed between different PA molecules ( Figure 2 in main manuscript) can also stabilize the PA structure and lead to increase in the hydrogen bond characteristics between PA and PA.Comparison of type (2) and (3) interactions, where the role of both donor and acceptor, respectively, is played by PAs suggests that the hydrogen correlation for type (2) decays slower as compared to hydrogen correlation for type (3). This slower decay dynamic of type (2) hydrogen bonds indicates that the hydrogen bond formed between hydrogen of amide groups of amino acids (Npa(-Hpa)) and oxygen atom of water (Ow) are more stable as compared the hydrogen bond formed between carbonyl oxygen (Opa) and hydrogen of water (Hw).

Order Parameter for Water in Stage 3
Supplementary Table 4 show the translational and tetrahedral order parameters near hydrophilic groups of hexagonally packed fibers during last 10 ns of stage 3.
Supplementary Figure 8 show the snapshots of water molecules in the PA network at the end of stage 3. Both Supplementary Table 4 and Supplementary Figure 8 suggest that water in the hexagonally packed fibers is hydrogen bonded with each other to form ordered network. Such extreme level of ordering of water is typical of that seen for water near solid interfaces.

Vibrational Spectra of Water
Supplementary Figure 10 shows the vibrational spectra of hydrogen of water molecules at the end of stage 1, at the end of stage 3, and bulk water. In the case of pure bulk water, similar to previous studies, the libration, bending, and stretching band shows characteristic peaks at ~390 cm − 1 , ~1820 cm − 1 , and ~3390 cm − 1 , respectively. 8 The libration band for hydrogen of proximal water at the end of stage 3 (~520 cm -1 ) shows a blue shift as compared to the proximal water at the end stage 1 (~460 cm -1 ) and bulk water (~390 cm -1 ) (Supplementary Figure 10 (a)). This suggests that at the end of ~150 ns of stage 3 due to increased ordering of water molecules the hydrogen bond network within water molecules is stronger as compared to the bulk water and proximal water at the end of stage 1 (Supplementary Table 4  Supplementary Figure 11 shows the vibrational spectrum for the oxygen atom of water molecules in bulk water and proximal water at the end of stage 1 and stage 3. The low-frequency (< ~600 cm − 1 ) end of the vibrational spectrum for oxygen atom of water shows two broad bands representing the O···O···O intermolecular bending motions of Hbonded molecules and the O···O intermolecular stretching mode at ~60 cm − 1 and ~ 250-300 cm − 1 , respectively. 8 The peak for the bulk water observed at ~60 cm − 1 shows a blue shift for proximal water present at the end of stage 1 and stage 3 at ~70 cm − 1 and ~90 cm − 1 , respectively. This blue shift suggests the reduction in the tetrahedral order and changes in the distribution of hydrogen bonds in ordered proximal water near PA molecules (Supplementary Table 4). In addition a shoulder for the O···O intermolecular stretching mode was observed at ~260 cm − 1 for bulk water. This peak shows blue shift for the proximal water molecules present at the end of stage 1 and stage 3 at ~270 and ~ 320 cm − 1 . This suggests that a stronger hydrogen-bonding network is present in proximal water near the PA molecules chain as compared to bulk water. 10

Experimental Infra Red Spectra of Dried Sample
In follow the vibrational trends that we observed in our molecular dynamics simulations.

Thermogravimetric Analysis
TGA experiments suggest that ~0.3 mg water is trapped in the fibrous assembly (Supplementary Figure 13). We also dehydrated the sample with excessive vacuum application and found a decrease in the overall signal.

Sources of uncertainty in our model
The following are some of the important factors that one needs to consider while implementing a multi-stage couple all-atom coarse-grained MD approach: Compatibility of all-atom and coarse-grained force-fields: It is very important that both all-atom and coarse-grained force-fields are compatible with each other and will preserve the structural and dynamical properties of the a given system throughout the course of they study. In the past, the combination of both CHARMM and MARTINI has been used successfully to study proteins and lipids. 12 Hence, in the present study we have utilized CHARMM and MARTINI force-fields for all-atom and coarse-grained models, respectively. Indeed, we between two closely packed fibers can be significantly different that the bulk water. Hence, while reinserting water molecules in the simulation cell that is consisting of all-atom self-assembled structure of PA nanofibers, we maintain the higher density of ~1.1 -~1.2 g/cc water between the two fibers. However, the overall water density of kept at 1 g/cc.

Calculations of Order Parameters
We find that the water near PAs form ordered structure. We calculate the translation and orientational order parameters for water molecules. Here are the calculation details:

Translation order
The translational order parameter (t) is given by: where, g OO (r) = the oxygen-oxygen radial distribution function.
s (dimensional variable) =rρ 1/3 , r = radial distance that is scaled by the mean s c = distances which are large enough to have g(sc)≈1. In this work, s c~1 0 Å was chosen.
In an ideal gas, the RDF is equal to 1 and t=0. In a crystal, there is long-range order and g OO (r)=1 over long distances, and order parameter t is large. The physical meaning of t is such that both positive and negative oscillations of g OO (r)=1 and contribute equally to this quantity. In a system with long-range order e.g., a crystal, these oscillations persist over long distances, and hence t is large. The range of the translational order can vary from ~1.2 (for bulk water) to ~3.4 for hexagonal ice.

Orientational order (q)
We calculate the orientational order parameter q to evaluate the differences in the tetrahedrality of the structured water near PAs. The orientational parameter qdescribes the ability of the water to produce tetrahedral arrangements and can be given as: 2 Where, θ jik = angle formed by the lines joining the oxygen atom of a given molecule with that of each of the four nearest neighboring molecules. In this study, the ensemble average is obtained over all the molecules and the timeframes of the simulation. The value of q would effectively range from 0, in a random distribution of molecules to 1, in a perfect tetrahedral arrangement. Note that that the value of this order parameter is not subjected to the choice of any particular local coordinate system in the reference molecule.

Calculations of Potential of Mean Force of PA Aggregation
The PMF between two PAs can be written as a function of the separation between them, W = f(r). The mean force, F(r) between these two atoms is then defined as: The PMF between two PAs can be calculated by using following equation: Where, r = distance between the hydrophobic tails of two different PAs. To calculate the entropy we used the finite difference temperature derivative of the PMF, W(r) at each inter-atomic distance between the group of atoms on hydrophobic tails, r, as:
The enthalpy contribution to the free energy, H(r), can then be obtained from entropy S(r) and the PMF, W(r) at temperature T as: Here, values of T and ΔT were chosen to be 280 and 40, respectively.

Coarse-grained Model of PAs
The MARTINI CG model has been previously used to study the self-assembly of peptides. 7 This model uses a four-to-one mapping scheme, wherein four heavy atoms (e.g. C, O, N) and their associated hydrogen atoms are represented as a single CG bead.
Four water molecules are also mapped to form one water bead and are treated explicitly. In the third stage of self-assembly process we removed all the CG water molecules from the simulation cell. This was followed by back-mapping of the CG model of PA to fully-atomistic model by using VMD's Coarse Grain Builder tool. 9 This all atom structure was further solvated by randomly distributed modified TIP3P water molecules.
Simulations were carried out using NAMD simulation package using NPT ensemble for 150 ns with 1 fs timestep. The parameters for the NPT ensemble are identical to the one described in first stage of the self-assembly simulations.

Calculation of Residence Time Probability P res (t)
Residence probability (P res (t)) of water molecules near the hydrophobic tail and amino acids was calculated at the end of stage 3. Water molecules were divided into two different regions: (1) hydrophilic water molecules (water molecules located within 5.5 Å of hydrophilic groups) and (2) hydrophobic water molecules (water molecules located within 5.5 Å of hydrophobic groups). Note, alanine, histidine, leucine, and lysine amino acids were included in the definition of hydrophilic region of the PA and hydrophobic region is consisting of aliphatic chain. P res (t) was defined as the probability of continuous existence a water molecule in one of the above defined regions for the time interval t 0 and t + t 0 . 9 P res (t) was averaged over several different simulation frames that were collected in intervals of 1 ps over a time frame of last ~10 ns of ~150 ns simulation run.

Hydrogen Bonding Criteria
Hydrogen bonding characteristics were evaluated using following geometric criteria for hydrogen bonding: 10 In supplementary equation 7, R OO and R OH are the distances between oxygen and oxygen and between oxygen and hydrogen of water 1 and 2, respectively, which are involved in hydrogen bond formation. The angle between oxygen of water 1 and oxygen & hydrogen of water 2 (O 1 ⋅⋅⋅⋅⋅O 2 −H 2 ) can be defined by angle φ. In the present study, the cut-off distances, R OO and R OH were defined based on the first minimum in the corresponding radial distribution functions of pure water.
Here, total three types of hydrogen bonds are possible. One hydrogen bond between PA and PA (Type 1) and two types of hydrogen bonds are possible between peptide and water (1) bond between oxygen (Opa) or nitrogen (Npa) atoms of alanine, histidine, leucine, and lysine amino acids and water hydrogen (Hw) (Type 2) and ( A hydrogen bond occupation number (S ij ) can be defined to distinguish between the hydrogen bonded and non-hydrogen bonded atoms.

⎧ ⎨ ⎩
The S ij that describes the existence or nonexistence of bonds between a selected donoracceptor pair ij (supplementary equation 8). A time-dependent autocorrelation function of this state variable S ij can be defined to determine the stability of the hydrogen bonds. We use the continuous correlation functions to define the S ij . In the case of continuous correlation, the hydrogen bond occupation number (S ij ) was allowed only one transition from 1 to 0 when the bond between atom i and j breaks for the first time. After the first breakage of hydrogen bonds between atom i and j, S ij was never allowed to return to 1.

Calculations of Vibrational Spectra of Water
The vibrational spectra for water molecules were calculated by Fourier where, v j (t) represents the velocity at time tfor atom type jand the angular brackets denote an average over the ensemble and over time 0. This power spectrum will also contain anharmonic contributions. A separate MD simulation was conducted on the final configuration generated at the end of the production run of 150 ns (both at the end of stage 1 and at the end of stage 3) to calculate the vibrational spectra. Simulations were carried out with a time-step of 0.5 fs at 300 K for 100 ps with velocities for each atoms being stored at every time-step. Water molecules within 5 Å of PA atoms were selected as proximal water molecules. Similarly, we have also calculated the vibrational spectra for bulk water at 300 K to facilitate comparison with the spectra obtained for proximal water.