“Crawling” on the self-assembly system: A molecular simulation of peptide position adjusting over self-assembly block

By combining non-equilibrium molecular dynamics(NEMD), umbrella sampling, and weighted histogram analysis method(WHAM), we calculated the potential of mean force of histidine peptide moving over a self-assembly structure. The reaction coordinate is along the main chain direction of the histidine peptide in the self-assembly structure. It is found that the energy needed for the histidine peptide with 3 and 5 residues while moving along the reaction coordinate is around -2.2 kCal/mol and -7.4 kCal/mol, respectively. And the histidine peptide crawls along the reaction coordinate, performing a snake-like movement. This result could illustrate how histidine peptide adjusts its position during self-assembly process.


Introduction
Self-assembly phenomenon is common in nature and many materials with excellent optical, electrical, catalytic properties such as nonlinear optics, chemical biosensors, information storage materials, and tissue growth scaffolding materials, are produced by self-assembly. Molecular self-assembly is the process by which molecules adopt a defined arrangement without guidance or management from an outside source. Noncovalent interactions (e.g., hydrogen bonding, metal coordination, hydrophobic forces, van der Waals forces, π-π stacking, and/or electrostatic) are crucial to the self-assembly process. [1] Histidine is an α-amino acid that is used in the biosynthesis of proteins. It contains an αamino group, a carboxylic acid group, and an imidazole side chain (which is partially protonated), classifying it as a positively charged amino acid at physiological pH. Histidine was first isolated by German physician Albrecht Kossel and Sven Hedin in 1896. The imidazole ring of histidine is aromatic at all pH values [2] and it can form π stacking interactions [3].
The self-assembly of histidine peptide is by hydrogen bonding between the backbone of the peptide, forming an antiparallel β-sheet second structure. What we are interested in is the self-assembly process. For the initial stage of the self-assembly process is always at random, and the final result is "neat" structure. To find out how histidine peptide adjust itself into neat structure is the target of this article. One guess is that histidine peptide is like Lego peptide [4] which can "slide" on the surface of the assembly only that instead of hydrophobic effect of the residues, there are other mechanisms for histidine peptides. Thus, we perform the following simulation, trying to find out the energy needed for the histidine peptide to adjust its position over the assembly block.

Constructing histidine peptide and histidine self-assembly block
We use full atom model to do our simulation so that a clearer picture can be obtained. First, we obtain the peptides with 15 histidine residues to build an antiparallel β-sheet second structure. Then the structure is equilibrated in water at 300K and 1.0 bar of pressure for 10 ns. Two layers in the middle of the structure are taken out, one as our "reaction" surface, the other being the source of histidine peptide which will be pulled over the reaction surface, showing in Fig 1-A Then 2 types of histidine peptides are "cut" (unnecessary residues removed and only the peptide in the middle of the three chains is kept) from the layer described above, producing 2 different lengths (with 3,5 residues) of histidine peptide. The initial model showing in

Umbrella pulling
To calculate the PMF of the peptide along the reaction coordinate, a non-equilibrium pulling method is used. [5].
When the distance between two groups is changed continuously, work is applied to the system, which means that the system is no longer in equilibrium. However, one can use the Jarzynski relation [6] to obtain the equilibrium free-energy difference ∆G between two distances from many non-equilibrium simulations: where WAB is the work performed to force the system along one path from state A to B, kB is the Boltzmann constant, T is temperature and the angular bracket denotes averaging over a canonical ensemble of the initial state A and β = 1/kBT The pull code applies forces or constraints between the centres of mass of one or more pairs of groups of atoms. In this article, an umbrella pulling is used, in which case a harmonic potential is applied between the centres of mass of two groups. Thus, the force is proportional to the displacement. And the peptide, in our case, is "dragged" by the force along the reaction coordinate. Fig.2 illustrates the umbrella pulling simulation. The alpha-carbon of the first residue is attached to a virtual spring with the force constant being 1500. The pulling is in the direction of the main chain and the other two directions' movements are restrained. The pulling simulation duration is 500 PS of MD, saving snapshots every 1 PS and the distance pulled is 5.0 nm. After the pulling simulation, a series of configurations are generated and selected so that configurations are of 0.2nm spacing. Fig. 2. The pulling simulation. The red arrow represents the force applied by the virtual spring.

Umbrella sampling
The binding energy (ΔGbind) is derived from the potential of mean force (PMF), extracted from a series of umbrella sampling simulations. A series of initial configurations is generated along the reaction coordinate ξ, as described above, each corresponding to a location wherein the ligand is harmonically restrained at increasing centre-of-mass (COM) distance from a reference molecule using an umbrella biasing potential. This restraint allows the ligand to sample the configurational space in a defined region along a reaction coordinate between it and its reference molecule.
In this article, we sample COM distances from 2.0 -5.0 nm (and 1.1 -3.4 nm for the 3residue model) along the x-axis using roughly 0.2-nm spacing. Start by running a brief NPT equilibration in each window. Then each window is sampled 10 ns.

PMF extraction using WHAM
The Weighted Histogram Analysis Method (WHAM) is one of the earliest methods that take into account information from all Intermediate States. The precursor to WHAM and first version of multiple histogram relighting techniques came from Ferrenberg and Swendsen; [7] WHAM was developed later for alchemical simulations [8]. WHAM is included in GROMACS as the wham utility, which this article uses to extract the potential of mean force.

Results of Molecular Dynamics Simulation
The potential of mean force profile is extracted by WHAM. We can thus calculate the energy needed for the histidine peptide to "crawl" on the reaction surface, which is -7.4kCal/mol for peptide with 5 residues of His and -2.2kCal/mol for peptide with 3 residues. From the PMF curve, we can also see that the process is periodic, which is reasonable as the peptide moving along the reaction coordinate, antiparallel β-sheet structure is periodically formed. The bottom of the curve means that such antiparallel β-sheet is formed and the peak means that nearly all of the hydrogen bonds between the main chain of histidine peptide and the assembly block are broken.
From the simulation result we can see that the peptide moves like snake to adjust its positions along the reaction coordinate. The simulation result shows in Fig.3-A,B.   Fig. 3. A: The peptide with 5 residues of His. "crawling" along the reaction coordinate, with labels from "a" to "e" indicating different state on the PMF curve. B: Peptide with 3 residues, from one stable state to another, and PMF curve. C: Backbone molecules' hydrogen bonds amount change in one period, with black representing peptide with 5 His, red 3 His. The Position 1 to 5 are positions shown in A, from position "a" to position "e", respectively.

Conclusions
From the PMF curve obtained above, we can find that the energy needed for the histidine peptide to adjust its position along the reaction coordinate is not much (-7.4kCal/mol and -2.2kCal/mol for peptide with 5 His residues and 3 His residues, respectively). This means the peptide could easily move on the assembly surface and find its stable position by "crawling" and form a neat β-sheet structure. And the difference between the peptides of two different lengths (3 residues and 5 residues), we guess, probably suggests that the energy required to make position adjustment is in associate with the hydrogen bonds amount which is formed between the backbone atoms of histidine peptides. The hydrogen bonds (on backbone atoms) average amount's change in one period of the "crawling" process proves that what we guess is true.
Firstly, the hydrogen bonds of the residue which the virtual spring is attached are broken, causing the hydrogen bonds on the residue next to it broken too. The "free" residue moves to either side of the reaction layer due to the effect of the side chain interaction. From the Fig.3-C we can see the amount change of hydrogen bonds in one moving period, and position 1 to 5 are .
To sum up, the histidine peptide could adjust its position via interacting with other histidine peptide chains in the self-assembly surface, doing a snake-like "crawling". This could indicate that in the self-assembly process, the histidine peptides would form the assembly block in a similar pattern.