Distribution and Structure Analysis of Fibril-Forming Peptides Focusing on Concentration Dependency

We focus on the concentration dependency of fibril-forming peptides, which have the potential of aggregation by themselves. In this study, we performed replica-exchange molecular dynamics simulations of Lys-Phe-Phe-Glu (KFFE) fragments, which are known to form fibrils in experiments under different concentration environments. The analysis by static structure factors suggested that the density fluctuation of the KFFE fragments becomes large as the concentration increases. It was also found that the number of β-structures and oligomers also increases under a high concentration environment. Hence, a high concentration environment of fibril-forming peptides is likely to cause protein aggregation.


INTRODUCTION
Protein aggregation has become one of the most important research subjects in recent years because it is known that the phenomenon leads to serious diseases such as Alzheimer's disease, Huntington's disease, Parkinson's disease, and type II diabetes, etc., which are referred to as amyloidoses. Generally, the structure of aggregated proteins, which is known as amyloid fibril, is a cross-β-sheet structure with the β-strands perpendicular to the fibril axis. 1−3 The fibrillization mechanism has been elucidated by many researchers. 4−7 According to a kinetic analysis, the process of amyloid formation follows the reaction time of the sigmoidal form. 8 The behavior of the reaction has a lag time, which is observed before a rapid growth phase and exhibits a feature of nucleated polymerization. 9 As with general crystallization phenomena, the addition of seeds shortens the lag time. In addition, it is suggested that amyloid fibril formation is a property determined by the concentration of peptides or proteins relative to the thermodynamic solubility. 8,10 If the concentration of dissolved proteins is larger than the thermodynamic solubility, the protein becomes more stable in the amyloid state than in its native or dissolved state. Namely, this also implies the presence of the critical concentration. Moreover, if the nucleation process in the amyloid formation as the initial stage is similar to the general crystallization, the effect of the density fluctuation is also important. The presence of large density fluctuation at the critical temperature affects the route to the nucleation by decreasing the free-energy barrier. 11−15 The analysis of the fibrillization mechanism has also been done using molecular simulations. 16−22 However, there are differences between simulations and experiments in regard to protein aggregations. One is the concentration. 22 Generally, concentrations of fibrils such as amyloid β, which is involved in Alzheimer's disease, and α-synuclein, which is involved in Parkinson's disease in cerebrospinal fluid, are 10 −11 −10 −9 M. In the case of in vitro studies, the order of concentrations is 10 −6 −10 −4 M. That of in silico studies is 10 −3 −10 −1 M. Another difference is the time scale. Indeed, a lag time of the amyloid formation in experiments is over 10 h. 8 Thus, some thermodynamic physical quantities obtained from simulations may not be quantitatively in agreement with those of the actual in vivo or in vitro environments. Despite the difficulty of comparing the concentration of protein aggregation, however, molecular simulations are very helpful as a way of analyzing the inter-and intramolecular structures in detail. We have performed a qualitative comparison of molecular simulations between different sizes of the simulation boxes, which include the same number of fragments. In our previous results, the concentration dependency of an amyloid-β peptide fragment (Aβ 25−35 ) by a molecular simulation has been studied. 23 The simulation under the high concentration of eight Aβ 25−35 resulted in an increase of β-structures. In this study, we examine the concentration dependency of fibril-forming tetrapeptides using a molecular dynamics (MD) simulation. In addition to atomic-level analysis, moreover, we focus on density fluctuation at the initial stage of the aggregation process. The peptide is the KFFE fragment, which consists of only four amino-acid residues, which has been known as a minimum-size fibril-forming peptide. 24 The dimerization mechanism for some kinds of tetrapeptides including the KFFE fragment has been studied by molecular simulations. 25−28 The simulation method that we employed is replica-exchange molecular dynamics (REMD), 29 which is one of the efficient conformational sampling techniques. In order to perform simulations at different concentrations, we simply prepared simulation boxes of several different sizes and the same number of fragments (= 30), in other words, the different number density of fragments (see Figure 1). By comparing the simulation results with different concentrations, we examined the differences of the distribution and structure characteristics of the fragments.
This article is organized as follows. Section 2 describes the details of the methodology, and Section 3 gives the results and discussion. Finally, Section 4 presents our conclusions.

METHODS
Static Structure Factor. In order to examine the distribution characteristics of KFFE fragments during the simulations, we calculated the static structure factor, which is a useful tool for scattering studies such as X-ray, electron, and neutron diffraction experiments. In this study, we investigated the density fluctuation of the KFFE fragments by using the static structure factor S(q) defined by which is related to the fluctuation in a quantity ρ̃and the number of fragments N. ⟨···⟩ is the ensemble average. See the Appendix for the derivation of this formula. We can separate the number density ρ into the average ρ 0 (= N/V) and the fluctuation δρ around ρ 0 where ⟨δρ(r)⟩ ≡ 0 by definition. In the case of q → 0, we obtain the static structure factor by eq 2 where the delta function δ(q) in the first term in the second line can be treated as zero because the wavenumber q does not have a value of zero at the actual simulation system. Here, δN is the number fluctuation (δN = N − ⟨N⟩). Namely, in the case of a low-q limit, if S(q) is large, the density and number fluctuations are also large. The number fluctuation in the general thermodynamic relation 30 is Here, k B , T, and κ T are, respectively, Boltzmann's constant, temperature, and isothermal compressibility: From eqs 3 and 4, S(q) also has a relationship with κ T . Therefore, we may consider that the high number density of the fragments increases κ T . Radial Distribution Function. As with the static structure factor analysis, the radial distribution function is also helpful in characterizing the properties of dense fluids and materials.
Here, Δn i (r) is the number of other peptides except reference peptide i in the volume (= 4πr 2 Δr) of a spherical shell of radius r and thickness Δr. V is the volume of the system. Orientational Order Parametar. To investigate the orientational characteristics of fragments, we calculated an order parameter Q for a pair of KFFE fragments by using the second Legendre polynomial P 2 as follows where θ ij is an angle formed by vectors i and j, which represent the orientation of each KFFE fragment, and r is a distance between two centers of mass of each KFFE fragment (see Figure 7(b) and its caption for a detailed definition). Structural Analysis of the β-Structure. We employ the DSSP algorithm 31 in order to examine the β-structures formed from the KFFE fragments during the simulations. The secondary structures such as helix structure and β-structure are held in shape by hydrogen bonds in the polypeptide backbone. In the DSSP algorithm, the hydrogen bonds in the backbone atoms are defined by using the distances between the CO group and N−H group as follows where r ON , r CH , r OH , and r CN are the distances of O−N, C−H, O−H, and C−N, respectively. The hydrogen bond is identified if E is less than −0.5 kcal/mol. In a β-structure, segments of a polypeptide chain line up next to each other, forming a sheetlike structure held together by hydrogen bonds. The βstructure is parallel, pointing in the same direction, or antiparallel, pointing in the opposite direction of polypeptide chains. Moreover, we consider a cluster of oligomers, which consists of the KFFE fragments forming β-structures. Here, the fragments are consecutively connected to each other by the hydrogen bonds in the β-structures. Namely, an oligomer n(nmer) indicates that n KFFE fragments have β-structures and form a network by the hydrogen bonds.
Simulation Details. In order to study the concentration dependence of the KFFE (Ace-Lys-Phe-Phe-Glu-NH 2 ) fragments, we have performed REMD simulations of the system of 30 KFFE peptides with four different volumes (50 3 , 60 3 , 80 3 , and 100 3 Å 3 ) with periodic boundary conditions. For our simulations, the NAMD 2.1 program package was used. 32 The force field for the KFFE fragment was AMBER ff14SBonlysc, 33 and the solvent model was the GB/SA (Generalized Born/ Solvent-Accessible Surface Area) solvent model with the Generalized Born of OBCII parameter set. 34 For all bonds involving the hydrogen atoms in the fragments, we imposed the constraints by using the SETTLE algorithm 35

RESULTS AND DISCUSSION
In this study, we performed four REMD simulations of different number densities, 240, 139, 59, and 30 (10 −6 /Å 3 ). The acceptance ratio for each replica lied between 0.2−0.4 in all cases (details are listed in Table S1 of the Supporting Information). There were no significant changes in the acceptance ratio for different number densities, and all simulations were performed properly (see also Figure S1 and Table S2).
Static Structure Factor. In Figure 2, we show the q dependence of the static structure factors S(q) of the KFFE fragments in eq 16 in the Appendix from the simulations with four number densities. When the wavenumber is less than about 0.3 Å −1 , S(q) values are more than 1.0, and the lower the wavenumber is, the larger S(q) tends to become. In addition, the values of S(q) at the high concentration of the number density become larger in magnitude.
From these results, we consider that the simulation system with high number density of KFFE fragments has high number fluctuations in Figure 2. In addition, we calculated S(q) for several temperatures from 300 to 450 K at number density 240 (10 −6 /Å 3 ) in Figure 3. It seems that S(q) becomes low at the low-q limit, while the fluctuations also become low as the temperature increases. In order to investigate the S(q) values in more detail, we plotted the S(q) values for each direction (k x , k y , and k z ) of wave vector q in Figure S2 of the Supporting Information. The results show that there is no dependence on the direction of wave vector q at any number densities. Since these S(q) values are isotropic, we consider that the results in Figures 2 and 3 are sufficient to indicate the characteristics of the S(q) values.
In order to examine the finite-size effects on the density fluctuations, we also performed an addtional REMD simulation, which had almost the same number density as one of the calculated systems but had a different number of KFFE fragments (20 KFFE peptides) and a different box size (about 70 Å 3 ). In Figure S3 of the Supporting Information, we compared between the two systems. The results show that these factors are in good overall agreement.
Radial Distribution Function. In Figure 4, the g(r) of KFFE fragments at 300 K is illustrated. All the simulation systems of the number density from 240 to 30 (10 −6 /Å 3 ) have a comparatively large peak between 5 and 10 Å. Moreover, there are two small peaks in the large peak for the case of the highest number-density system (see Figure 4(a)). These two peaks are at r = 5 and 10 Å. We also calculated g(r) of the truncated structures, which have only dimers and only trimers for the highest density system, ρ 0 = 240 (10 −6 /Å 3 ), in Figure 5. We used the DSSP program 31 for the criteria for secondary structure formations. In the case of only dimers of KFFE fragments, there is clearly one peak at r = 5 Å. On the other hand, there is another peak at r = 10 Å as well as the peak at r = 5 Å in the case of only trimers. In other words, as a dimer is a subset of a trimer, trimers have not only the peak at r = 10 Å but also the peak at r = 5 Å from dimers. Namely, we see that the characteristic distances of dimers and trimers are r = 5 and 10 Å, respectively. It seems that the high density system has dimers and trimers of KFFE fragments frequently. We compare g(r) at temperatures from 300 to 400 K in Figure 6. As the temperature increases, the peaks of dimers and trimers decrease. Especially, the peak of dimers clearly disappears. Due to the temperature rise, the peaks have a constant value from around 10 Å. This indicates that the peptides are randomly distributed.
Orientational Order Parameters. In Figure 7, we show the orientational order parameter Q of KFFE fragments at 300 K. In Figure 7(a), we see that the simulation system with high number density of KFFE fragments has high order parameter Q. The tendency is particularly pronounced when r is short (less than 6.0 Å). On the other hand, Q values of all the cases are close to zero when r becomes larger than 10 Å. The results show that the orientations of peptides are aligned at short distances between the peptides, and the tendency is greater for higher number density environments.
Structural Analysis of the β-Structure. In Figure 8, we show the number of β-structures during the simulations for four different volumes as a function of number density. The structures included all the β-bridge and β-ladder structures, which were identified by the DSSP program. 31 We see that the frequency of the β-structure increases as the number density increases. The number of oligomers is also listed in Table 1. Most numbers of oligomers increase as the number density increases. Moreover, the number of n (= n-mers) also increases. Here, n-mer indicates the cluster, in which n fragments are connected by the hydrogen bonds. As references, we illustrate some snapshots including the oligomers of the KFFE fragments in the case of number densities 240 (10 −6 / Å 3 ) in Figure 9 and 59 (10 −6 /Å 3 ) in Figure 10. For the obtained β-structures, we separate for more differences of parallel or antiparallel (see Table 2) and bridge or ladder (see Table 3). Usually, it is considered that the β-structure of KFFE fragments is stable in the form of antiparallel conformation because of the electrostatic interactions of the side chains of the terminal residues K (lysine) and E (glutamic acid). Although the obtained β-structures are in agreement with the general view, the number of parallel conformations slightly increases in the case of the high number density. The high concentration of the KFFE fragments may result in an increase of the other interactions such as hydrophobicity. In the case of the parallel β conformation, the electrostatic forces are repulsive. However, under the high-density conditions, KFFE fragments are difficult to move freely due to their own Lennard-Jones repulsive term. Generally, the Lennard-Jones repulsive force is much stronger than the electrostatic force. Thus, some fragments form the parallel β conformations by a hydrophobic effect, even though they are electrostatically unstable. This tendency is represented in Table 2.
In the past experimental results of KFFE, KFFK, and EFFE fragments, 24 fibril formation of KFFE was more pronounced, indicating that charge attractions are important for fibril formation. The peptides KFFK and EFFE do not form fibrils when incubated individually. However, coincubation of equimolar amounts of KFFK and EFFE produces fibrils as detected by EM (electron microscopy). Therefore, the authors have considered that β-strand structure in solution and attractive electrostatic interactions are required for fibrillo-genesis. Our results are consistent with these experimental results, with the antiparallel structure predominant in Table 2. Moreover, the model of a β-structure predicted by the experimental results 24 is stabilized not only by electrostatic interactions between the side chains of K (lysine) and E (glutamic acid) but also by van der Waals interactions between F (phenylalanine) side chains. We see β-structure conformations of the oligomers of the KFFE fragments in our simulations, e.g., Figure S4 of the Supporting Information. The side chains of the phenylalanine are close to each other. The conformations of the oligomers obtained from our simulations as well as the predicted model by the experiment suggest the importance of the hydrophobic interaction by the side chains of phenylalanine.
The formation of KFFE produces sedimentable amyloid fibrils (1.2−1.6 nm in width) detected by EM. 24 We estimated the average width (i.e., length of the end-to-end distance of  KFFE fragments) in the heptamer conformation (see Figure  S4), which was 1.351 nm. Thus, the fibril conformations of KFFE seem to have a stack of one continuous fragment.
In Table 3, we list the fragment of the bridge or ladder in the β-structures. From the table, we see that the β-ladder structures increase in contradiction to decreasing the β-bridge structures as the number density becomes high. We consider that the stabilities of not only the oligomerizations but also βstructures between fragments themselves increase. In Figure  11, we show the free energy of the oligomer formation which is computed directly from the distribution of oligomers of KFFE fragments, following ref 36. We chose to compute the series of equilibrium constants corresponding to the addition of one fragment (F 1 ) to an oligomer (F n ): Equilibrium constants for the above series of reactions are given by which gives where k B is the Boltzmann's constant and T is the temperature. In Figure 11, we see that the stability of forming oligomers increases as the number density increases. Moreover, as the size of the oligomer increases, the KFFE fragments tend to be added more easily.
Increased β-Structure in a High-Density Environment. In the simulation of the case of high-number density 240 (10 −6 /Å 3 ), some snapshots during the simulation have high S(q) values (≥9.0) unlike the case of low-number density 30 (10 −6 /Å 3 ). In Figure 12, we show the probabilities at high S(q) (≥9.0) and low S(q) (<9.0) for each number of βstructures. We see that in the case of low S(q) there are various numbers of β-structures with a peak of 4. On the other hand, in the case of high S(q), with 11 as a large peak, there are βstructures in the range of 7 to 18. Namely, when S(q) is low, the KFFE fragments have various conformations including fewer or more β-structures. However, when S(q) is high, the     conformations are biased toward the states of many βstructures. Although we could not confirm the fibrillization of the fragments in these simulations fully, we consider that a state that contains a lot of β-structures is a precursor of the fibril formation and that a high S(q), that is, a high density fluctuation, leads to the fibrillization that precedes the increased β-structures.

CONCLUSIONS
In this study, we performed the simulations focusing on the concentration dependence of a fibril-forming fragment, KFFE. By using statistical mechanical analysis, we showed the difference of the distribution of the fragments with four simulations, which have different fragment concentrations. The static structure factor analysis showed that the density fluctuations of the fragments become large when the fragment concentration becomes high. In addition, the isothermal compressibility may also increase. By the radial distribution function analysis, the distances between fragments of dimers and trimers are clearly characterized in the system of the high number density. On the other hand, the numbers of dimers and trimers of the fragments are obviously larger than the other oligomers (n-mers, n > 2) for all values of the number density by the structural analysis of the β-structure. Namely, it is suggested that the number of fragments, which have distances between fragments of dimers and/or trimers and do not have the β-structure of dimers and/or trimers, increases relatively when the concentration of fragments becomes low. For the structural analysis of parallel or antiparallel β-structure, the number of parallel conformations increases in the simulation of high number density. One of the reasons is that the hydrophobic interaction increases. Thus, we expect that the free-energy barrier of the nucleation decreases at high concentration by the presence of the large density fluctuations in crystal nucleation, and then the fibrillization proceeds smoothly. In a future work, we will proceed to the next simulations of Amyloid β peptides using an explicit solvent model in order to analyze the distribution and the structural characteristics of the peptides and the solvent. The spatial Fourier transform of the number density ρ, ρ, is expressed by the integral in volume V and given by The number density ρ is defined by Here, the position vector r i represents that of the center of mass of fragment i, q = (2π/L) (q x , q y , q z ) is the wave vector, where L is the box length, q x , q y , q z are integers, q k * = q −k , and δ(···) is the Dirac delta function. From eq 13, eq 12 becomes The ensemble average formally gives ρ ′ = ⟨ ′ ⟩ n r r r r ( , ) In the case where the system is both homogeneous and isotropic, the dependence on the coordinate vectors, r and r′, Figure 11. Change in free energy on adding a single fragment to an oligomer in the case of the simulation at 300 K with number densities of 240, 139, 59, and 30 (10 −6 /Å 3 ), which correspond to volumes of 50 3 , 60 3 , 80 3 , and 100 3 (Å 3 ), respectively. Figure 12. Probablitiy of the static structure factor S(q) of KFFE fragments for each number of β-structures at 240 (10 −6 /Å 3 ). Red and black bars stand for S(q) larger than 9.0 and less than 9.0, respectively. is given by the relative distance, r = |r − r′|. We then obtain the radial distribution function g(r) as follows: ρ = ′ = ′ g r g n r r r r ( ) ( , ) 1 ( , ) 0 2 (2) (19) For an isotropic system, eq 19 can be simplified by recognizing that the radial distribution function depends only on the particle separation and should also be normalized by the area of the sphere of radius r drawn around a reference particle: Here, Δn i (r) is the number of other peptides except reference peptide i in the volume (= 4π r 2 Δr) of a spherical shell of radius r and thickness Δr.
■ ASSOCIATED CONTENT