Investigations into the role of non-bond interaction on gelation mechanism of silk fibroin hydrogel

: Silk fibroin hydrogel not only has biocompatibility, but also has environmental response ability. It plays an important role in the development of material. The gelation mechanism of silk fibroin hydrogel is very important to textile and medicine fields. The molecular dynamics simulation was used to discuss the structure and non-bond interaction of silk fibroin hydrogel. The results show that the non-bond interactions between silk fibroin molecules and water molecules have certain influence on the formation of silk fibroin hydrogel. According to the hydrogen bond analysis, the hydrogen bonds are mainly formed between random coil peptide fragments at the two ends of silk fibroin molecules and residues 252–254 are the key residues. The electrostatic and polar solvation interactions between silk fibroin molecules plays a major role in cross-linking of the coil segments of two silk fibroin molecules.


Introduction
Silk is a kind of natural protein fiber with high strength. There are two silk fibroins and sericin around them [1,2]. Silk fibroin has excellent biocompatibility and degradability. It has been widely used in textile and biomaterial fields [3][4][5]. It has great potential in developing sustainable and environment-interactive materials [6].
Silk fibroin can be easily processed into various materials, such as nanoparticles and microcapsules [1,[7][8][9][10][11]. Hydrogel is composed of polymer with three-dimensional network structure and medium filled in the gaps of polymer segments. The medium can be gas or liquid, and most of them are liquid in general. Therefore, hydrogels can be regarded as swelling bodies with threedimensional polymer networks containing liquid. Almost all-natural gels and most synthetic gels are hydrogels with water as liquid medium. The silk fibroin hydrogel is a swelling body formed by water contained in the three-dimensional network. At present, silk fibroin hydrogel has been developed to construct functional tissues and has been used to solve many problems in drug release [12][13][14][15][16][17].
Silk fibroin hydrogel not only has biocompatibility, but also has environmental response ability, which plays an important role in the development of smart material. However, the gel time limit its application. It takes several days or even several weeks for the aqueous solution of silk fibroin with the concentration of 0.6-15% to transform into gel at room temperature [1,18]. Therefore, it is necessary to explore the gelation mechanism and promotes the formation of silk fibroin hydrogel.
Due to extensive hydrogen bonding, silk fibroin is a hydrophobic protein, 70% of the total components are alanine and glycine, with high crystallinity [19][20][21]. The silk fibroin molecule consists of light chain, heavy chain and glycoprotein. The ratio of them is 6:6:1 [22][23][24]. It is difficult to detect and quantify the interaction between these molecules by experimental methods. Molecular dynamics can provide atomic-scale details of atomic interaction and structure property, supplementing the details of the relation between them.
A four-stage physical gelation process was proposed by Liu et al [25]. Silk fibroin hydrogels have a multi-domain network with weak domain-domain interactions. By studying the process of silk fibroin gelation, they believed that silk fibroin gelation was divided into four steps. First, silk fibroin molecules aggregate into a nucleus. Then the crystal nucleus grows into β-crystallites. Subsequently, the nanofibers further grow and branch into a single-domain. Finally, these silk fiber networks penetrate each other and form a stable multi-domain. This process shows that the weak interaction plays a very important role in the formation of silk fibroin hydrogel. Therefore, an atomistic silk fibroin hydrogel model comprising of two anti-parallel β-sheet has been built in this work. The radial distribution, secondary structure, hydrogen bond and non-bond interaction have been used to describe the weak interaction and structure of the silk fibroin hydrogel based on the molecular dynamics [26].

Models and methods
The initial model was taken from the protein database (ID: 3UA0, Figure 1a). The molecular model comprises of two crystal structure of the N-terminal domain of silk fibroin with double antiparallel β-sheet [2]. The missing part of this model was supplemented by the homologous modeling method with Chimera ( Figure 1b) [27,28]. Two silk fibroin molecules were stacked in the minimum space ( Figure 1c). One trajectory at temperature 300K had been performed for 0-200ns. The equilibrium of the system was analyzed by energy and RMSD. The rationality of the model was proved, and the gel model was discussed by the hydrogen bonds between silk fibroin molecules. The secondary structure, radial distribution function, hydrogen bond and other parameters were used to characterize the structure and interactions of the gel.
The model was simulated by molecular dynamics at 300 K with Amber software [29]. The process of energy minimization and heating had been performed firstly, and then performed the equilibrium and product simulation in this work (Figure 1d). In the minimization stage, constrained and nonconstrained energy minimization were performed. The steepest descent method [30] was used in the first 2000 steps, and then converted to the 2000-step conjugate gradient method [31,32]. In the minimization of constraining energy, a force of 10.0 kcal/mol was used to constrain the position of the residues. The temperature of system was heated to 300 K from 0 K within 300 ps, and the same force was used for restraint. To prevent the temperature rising too fast, the process was divided into two steps. All bonds involving hydrogen atom were restrained by SHAKE algorithm [33]. The NPT ensemble was used in 5 ns equilibrium simulate. The simulation time of the product simulation is 200 ns with NPT ensemble. The integral time was set to 2 fs and the cutoff distance of non-bonding was set to 10 Å. The sampling interval was 10 ps. Berendsen algorithm [34] was used to control the pressure at 1 bar, and the temperature was controlled with the Langevin method [35] at 300K. Collision frequency for Langevin thermostat was 5.0 ps −1 . Pressure relaxation time was set to 10.0 ps. The force field ff14SB [36] and tip3p water cube box of 10.0 Å were used.

System equilibrium and verification
Root Mean Square Deviation (RMSD) has been designed to measure the coordinate root mean square difference of a specific atom relative to the reference structure, and to characterize the difference between the structure and the initial structure of protein. RMSD value of α-carbon atom of protein molecular chain is an important basis to measure the conformational change of protein in the simulation process [37]. The smaller the value, the smaller the deviation.
RMSD and potential energy of the system were adopted to verify whether the system has reached equilibrium ( Figure 2). The RMSD and potential energy tends to the constant values 12.5 Å and -181000 kcal/mol respectively during 150-200 ns. This indicates that the system reached equilibrium after 150 ns.
The label Total is the amount of every intermolecular hydrogen bonds ( Figure 3). The label M2-M1 indicates the hydrogen bonds which donor are provided by molecule 2 and acceptor are provided by molecule 1. The label M1-M2 is the opposite of the label M2-M1. The intermolecular hydrogen bonds are formed between silk fibroin molecules and this suggests that cross-linking is formed [38]. This model can be used as a silk fibroin gel model.  The radial distribution function was estimated by cpptraj of Amber. The high probability indicates that the particle is likely to agglomerate around the given particle [39,40]. The label M1-M2 represents the radial distribution function between O or N atoms of two silk fibroin molecules ( Figure 4). The label M-W represents the radial distribution function between O or N atoms of silk fibroin molecules and O atoms of water molecules.
The distance 2.6-3 Å represents hydrogen bonds and the distance 3-5 Å represents van der Waals forces. The Figure 5 shown that hydrogen bonds are formed between two silk fibroin molecules and between silk fibroin molecules and water molecules. Weak non-bond interactions such as van der Waals forces also exist. The radial distribution probability between two silk fibroin molecules is lower than that between silk fibroin molecules and water molecules. This indicates that the weak interaction of former is weaker than that of the later.

Conformation analysis
There is a gap between two silk fibroin molecules at the beginning and the random coil gradually fold together and fill the gap (Figure 5a). Some α-helix structures are around random coil during the folding process and β-sheets of the two molecules are relatively stable (Figure 5b). These results denote that cross-linking of silk fibroin molecules mainly occurred in random coil segments.   The distribution probability of 2.6-3 Å increases gradually during simulation (Figure 6a). This indicates that the non-bond interaction between silk fibroin molecules increases gradually with the formation and crosslinking of silk fibroin, while the non-bond interaction between silk fibroin and water molecules doesn't change significantly (Figure 6b). These results suggest that cross-linking occurs between silk fibroins. Because non-bond interaction between silk fibroin increases gradually and is weak in the early stage, the non-bond interaction between silk fibroin is weaker than that between silk fibroin and water molecules (Figure 4).

Secondary structure
It can be seen from the secondary structure probability diagram (Figure 7) that most of them are anti-parallel β-sheet and random coil, with a small amount of helix structure. It is consistent with the conclusion reported in the literature [41]. The β-sheet is a kind of stable secondary structure, which has a positive effect on the stability of silk fibroin structure [42,43]. Most of the interactions between molecules occur in random coils, thus they play a positive role in the gelation formation.
The structure of each residue can be divided into β-sheet and random coil at the stage of 100-200 ns. For the sake of statistics, more than three consecutive residues with anti-parallel β-sheet are defined as β-sheet segments, and others are defined as random coil segments ( Figure 8). The types and number of residues on β-sheet and random coil were counted [44]. Overall, the residues number of random coils is more than that of the β-sheets. There are more nonpolar residues in β-sheet peptide and more polar residues in random coil peptide (Figure 8). Compared with hydrophobic residues, both β-sheet and random coil have more hydrophilic residues. According to the analysis of structure and hydrogen bond, the polar residues of coil and hydrogen bonds of segments promote the aggregation of two silk fibroin molecules and form a cross-linked structure in the gap region.

Hydrogen bond
The total number of intermolecular hydrogen bonds is about 3-5 ( Figure 3). The detail of these hydrogen bonds is necessary to discuss (Table 1 and Figure 9). The residues ALA 50 and MET 253 are hydrophobic [44][45][46][47][48], and the others are hydrophilic. The occupancy of hydrogen bond between residue ASN 272 and TYR 241 is more than 50% in the process of 150-200 ns [49]. This result denotes that it is important for the formation of gel. The hydrogen bond between residue GLN 491 and ALA 50 appears after 170 ns and it is consistent with the literature report [50]. Because of the occupancy of intermolecular hydrogen bonds between the three pairs of residues in M2-M1 are more than 76% in the process of 150-200 ns, they are more important for the formation of gel. With the simulation, the number of intermolecular hydrogen bonds increases gradually. This indicates that silk fibroin molecules aggregate gradually, which promotes further formation of osmotic network.  Figure 9. Residues for stable intermolecular hydrogen bonds (red).
In the initial phase of the simulation, a large number of hydrogen bonds formed between the end of silk fibroin molecules and waters ( Figure 10). With the number of hydrogen bonds between them decreasing gradually, the number of hydrogen bonds between random coils of silk fibroin molecules increased gradually (Figure 3). There is no interaction among coil segments and they have great freedom to move freely in the early stage. The hydrogen bonds between the coil segments and water can pull silk fibroin molecules to aggregate, shorten the average distance of molecules, and create conditions for intermolecular hydrogen bonds forming. The intermolecular hydrogen bonds between silk fibroin molecules increase gradually and form cross-linked structures.
The probability of hydrophobic residues was 38% at the tail-end of molecule 1 and 42.1% at the head-end of molecule 2. It is easier for the tail-end of molecule 1 to form hydrogen bonds with water molecules. Because of the random coil at the head-end of molecule 2 (First2) is longer than that at the tail-end of molecule 1 (End1), there are more hydrophilic residues at End2, which can form hydrogen bonds easily with water molecules. Therefore, the ability to form intermolecular hydrogen bonds with water molecules has little difference ( Figure 10).

Other non-bond interactions
The non-bond interaction of the system was calculated by MMPBSA of AmberTools [51,52]. The contribution of each residue belonging to above hydrogen bonds was analyzed according interaction energies. The negative energy indicates attractive force, while positive indicates repulsive. Because the hydrogen bonding is closely related to electrostatic interaction, the secondary electrostatic interaction model was established by Van der Lubbe [53]. It has been used to predict and explain the hydrogen bond strength of self-assembled systems, and the experimental results are often consistent with the prediction of this model. The interaction energies of these residues were sorted according to electrostatic interactions and shown in Table 2.
The contribution of residues to electrostatic interaction is closely related to interaction of polar solvation. The residues with high electrostatic interaction have high polar solvation energy. The residues SER 252 MET 253 and ARG 254 have the highest electrostatic interaction among these 8 residues. The stronger the electrostatic interaction between residue pairs, the easier to form hydrogen bonds. The residue ALA 50 has the lowest electrostatic interaction, and it is difficult to form hydrogen bonds. Therefore, residues 252-254 can be considered as key residues.
Most of coil are non-sheet peptides, and they are located on gap region in conformation ( Figure  5a). The interaction of electrostatic electricity and polar solvation plays an important role in the aggregation of these residues. With the aggregation of silk fibroin molecules, hydrogen bonds crossand linking among coils are formed gradually and leads to the formation of gel structures. The contribution of electrostatic explains the reason that the silk fibroin solution can be induced to gel rapidly by surfactants [54][55][56].

Diffusion properties
Diffusion coefficient is an important physical parameter of biomolecules, which defines the diffusivity and permeability of the solvent [57]. The larger the molecular diffusion coefficient, the stronger the molecular migration ability [58]. Diffusion coefficient can not only describe the solution viscosity, but also relate to the intermolecular interaction [59]. Diffusion coefficient D can be estimated by Einstein relation: where, the unit of D is cm 2 /s. MSD is mean square displacement from initial position. The dimension degree n is 3 for all MSDs. The slope of MSD is obtained by linear regression. The diffusion coefficient of silk fibroin molecules is about 2.445 × 10 −5 cm 2 /s and that of water is about 2.6624 × 10 −5 cm 2 /s. There is little difference between the diffusion coefficient of water molecules and that of pure water solution (2.275 × 10 −5 cm 2 /s) [60]. The diffusion coefficient of water molecules is slightly higher than that of silk fibroin molecules, but it is still the same order of magnitude, which reflects the good molecular migration ability of silk fibroin hydrogel. The diffusion coefficient of silk fibroin hydrogel in different solvent in is around 0.2-8.1 × 10 −7 cm 2 /s [61][62][63]. The interactions between small molecules and silk fibroin lead to these diffusion coefficients are much lower than that of this work.
These results indicate that the diffusivity of silk fibroin hydrogel formation depends on the interactions between molecules and solvents.

Conclusions
In this paper, the gelation process of silk fibroin and the role of non-bond interaction in the gelation process were analyzed by molecular dynamics. The formation of silk fibroin hydrogel depends on non-bond interactions between silk fibroin molecules and solvent. The hydrogen bonds were mainly formed between the random coil peptide fragments at the ends of different molecules. In the early stage of simulation, hydrogen bonds with water were formed, which leads to silk fibroin molecules assemble and hydrogen bonds form. The residues 252-254 of the random coil were identified as key residues according hydrogen bond and non-bond interactions. In the gelation process of silk fibroin, the electrostatic interaction and polar solvation interaction among the random coil segments play an important role in the cross-linking and effect on diffusion coefficient. This is important for further understanding of the gelation mechanism.