Influence of repeating sequence on structural and thermal stability of crystalline domain of bombyx mori silk fibroin

Thermophysical properties of polar-antiparallel β-sheet crystalline domains of B.mori silk fibroin has been investigated via molecular dynamics (MD) simulations between 300 K and 700 K. In general, the type of interactions determines the character of the thermal expansion in corresponding directions except that the crystalline domains with serine residue contracts along the chain direction resulting with a negative thermal expansion (NTE) coefficient of α y ≈ − 10 − 4 K − 1 at 300 K. The heat capacity at constant pressure C P increases sublinearly up to about 650 K with an almost continuous decrease in the rate of change of entropy d S d T P . Energetic behavior of the β-sheet crystalline units is determined by the main chain below 400 K then thermally activated motion of side chain (HG1) of serine residue becomes effective up to the degradation point. These are particularly important in designing thermally controlled composite materials.

This repetitive primary sequence undergoes secondary structure of antiparallel β-sheets (figure 1(d)) that give rise to the extraordinary physicochemical properties. The crystalline segments, which are embedded into the amorphous region, are aligned after the natural spinning process (figure 1(d)) of silk larvae ( figure 1(a)). This leads to formation of fibers ( figure 1(b)) with a high tensile strength and modulus, while the amorphous flexible chains generally contribute to the elasticity of the silk fibers [28][29][30][31]. Therefore, it is of great importance to understand the silk structure before spinning, which is called Silk I (figure 1(c)), and after spinning Silk II (figure 1(d)), contribution of the repeating elements, and the role of individual amino acid to its unique structure and remarkable physicochemical properties.
Beside the great mechanical properties, it is know that the fiber is structurally stable up to 423 K under thermal treatment and undergoes a glassy state at the temperature of T g =428 K [32,33]. The thermal degradation temperature of the silk fibers occurs at about 523 K [32]. Moreover, in case of silk fibroin film, the glass transition of film form was observed for water cast [34] and dried samples [35] at 450 K and 451 K respectively. The crystallization of β-sheets was found to be at around 485-493 K [35][36][37]. In addition to the onset of degradation in the heavy chain (∼350 K) and light chain (∼400 K), β-sheet crystallinity remained stable up to 500 K [38,39]. It was reported that both strain-to-failure and the tensile strength of silk fiber decrease above ∼400 K while the modulus increases up to ∼500 K. Differential scanning calorimetry experiments show that B. mori fibroin has distinctive water evaporation peak just below the 373 K and thermal degradation peak at about 610 K with an indistinctive endothermic peak at about 500 K [40]. The glass transition of pure noncrystalline silk was observed at 451 K and a lower glass transition for the silk-water system [36]. Fast scanning chip calorimetry experiments have shown that the silk fibroin film displays a glass transition at ∼474 K and the onset of melting of native degummed fibers at around 610 K [41].
There are limited numbers of experimental and theoretical researches on the temperature evaluation of physical properties of B.mori silk fibers and the crystalline units of silk II. XRD study of thermal dependence of structure and mechanical properties revealed that the position of the B. mori [002] reflection is practically constant, and the peptide chain expansion coefficient along the c-axis (fiber axis) is significantly below the [210] and [010] expansion coefficients [38]. On the other hand, the structural and thermal properties of dragline silk such as spider silk and some wild silkworm silk (Samia Cyntia ricini) were mainly studied by the method of molecular dynamic simulation (MD) [42][43][44][45][46]. It is well known that there are two experimental model structure proposed for the crystalline unit of B.mori silk fiber as polar-antiparallel and antipolar-antiparallel β-sheets [24,27,[47][48][49]. Although these models are mainly composed of highly repetitive [GAGA(GS)] sequence, the polar-antiparallel form of (Gly-Ala) dipeptides is taken as the only repeating element to construct the model structure for the crystalline segments in examining the fiber formation B.mori silk by MD simulation technique [50,51]. The unit cell of crystalline domains was first described by Marsh et al [47] as regular antiparallel β-sheets with the rectangular unit cell parameters of a 0 =9.40 Å, b 0 =6.97 Å (chain direction), and c 0 =9.20 Å and the space group P . 2 1 However, the crystallites of silk II is better and more precisely characterized as less regular antiparallel β-sheets composing the crystallites by stacking of two antipolar-antiparalel sheets with the space group and parameters, essentially, the same as those of March et al [48,49].
Inclusion of serine residue in repeating dipeptide (GX) as in the ratio of 3G:2A:1S, [(GAS) 3:2:1 ], which is the most dominant sequence among the other repeating segments, clearly, will have direct influence on the studies by MD technique. Investigation of thermal properties such as thermal expansion, heat capacity, etc are of particular interest in order to further increase our understanding on the distinctive structural arrangement of simple amino acids as is known that the heat induces the formation of anti-parallel β-sheet domains. Its intriguing mechanical and biochemical properties have enhanced experimental and theoretical interest in understanding the structural conformation and the mechanics of the formation and the adaptation of silk fibroin in variety of application area, such as tissue engineering. In that respect, Asakura et al [27] pointed out in their perspective that there are important questions to be answered: What is the role and contribution of individual amino acids, repeating sequences and crystalline domains on the structural formation and functionality of silk fiber? How does Silk I evolves to Silk II? Lack of intensive theoretical research on the temperature evaluation of structural and physical properties in literature provide an opportunity to investigate crystalline segment of silk fibroin. This study aims to investigate the thermophysical properties, such as specific heat, rate of entropy change, enthalpic strain, and the coefficient of thermal expansion at atomic level via molecular dynamic simulation technique. The role of serine residue and different repeating elements were assessed by constructing a repetitive unit cell of polar-anti-parallel β-sheet crystalline regions made up by two different units; {[(GA) 3:3 ] 2×1×2 } 3×2×3 and {[(GAS) 3:2:1 ] 2×1×2 } 3×2×3 (figure 1(g)).

Construction of unit cell and simulation box
In order to obtain information on the role and contribution of individual amino acids and different module on the structural formation, functionality and physical properties under thermal treatment, two different crystalline domains were constructed from the repeating unit cells of The unit cell (a 0 , b 0 , c 0 ) is orthogonal with the monoclinic symmetry and P 2 1 space group, and consists of 8 amino acid residues of 4 polypeptide chains placed in a polar-antiparallel form. One more amino acid has been added to one end of each chain in order to construct the repetitive unit cell and an infinite sample in all directions for the purpose of this work (figure 2(a)). Thus, new unit cell contains 12 amino acid and 4 polypeptide chains of repeating dipeptides (GX) ( figure 2(a)). For the first model, which has been commonly  [50]. The simulation box including serine residues from different views. The fiber axis was denoted as y (green) and the other two axes, x (red) and z (blue) indicate the intra-sheet and inter-sheet directions, respectively. used in literature for MD simulation studies, X is, simply, taken as alanine. The second model structure, which includes serine residue, was made up by adding an oxygen atom to the carbon atom of functional group which gives the polypeptide chain of hexapeptide in the ratio of 3G:2A:1S ( figure 2(b)). Then hydrogen atoms were added to the side chains and charge neutrality was satisfied ( figure 2(c)). The unit cell in figure 2(c) contains 12 amino acid and 4 polypeptide chains which consist of 102 and 104 atoms for [(GA) 3:3 ] 2×1×2 and [(GAS) 3:2:1 ] 2×1×2 models, respectively. This procedure has been achieved by the help of VESTA (Visualization System for Electronic and Structural Analysis) [52] and AVAGADRO [53] software.
The simulation box was constructed by multiplying the unit cell six times along chain direction (fiber axis) b 6 , 0 and three times along the intra-sheet (hydrogen-bond) direction a 3 0 and intra-sheet (sheet stacking) direction c 3 . 0 Thus, each of the simulation boxes consists of 36 polypeptide chains, 432 amino acid residues. The periodic boundary conditions are applied to replicate the simulation box throughout space to form an infinite lattice in order to minimize the surface effects. The both ends of each polypeptide chain has to be connected along the chain direction. Therefore, we have defined the bonding between Carbon atom of carbonyl group (C=O) at one end and Nitrogen atom of amino (-NH 2 ) groups at other end. Through this procedure lattice samples of {[(GA) 3:3 ] 2×1×2 } 3×2×3 and {[(GAS) 3:2:1 ] 2×1×2 } 3×2×3 infinite in size were successfully created.

Simulations
Following the construction of the repetitive unit cell and simulation box for the models single crystals of polarantiparallel β-sheet crystallites {[(GA) 3:3 ] 2×1×2 } 3×2×3 and {[(GAS) 3:2:1 ] 2×1×2 } 3×2×3 , MD simulations were performed using GROMACS 5.1.4 package [54] at temperatures varying from 300 K up to 700 K. Prior to running the constant pressure-temperature (NPT) simulation, the constant temperature-volume ensemble (NVT equilibration simulation) was employed to relax the system at the temperature of 300 K for 50 ps period. The NPT simulation was then run at temperatures varying from 300 K up to 700 K with 20 K interval using final configuration of crystal structure as an input data for the simulation. Long range electrostatics interactions were handled with the Particle Mesh Ewald (PME) method. The NPT simulation was conducted for 50 ns for each temperature with a time step of 1 fs. The periodic boundary conditions were applied to replicate the simulation box of´á    Figures 3 and 4, respectively, contain comparative summary of thermal expansion and relative thermal expansion of two different model crystalline domains between 300 K-700 K. This denotes that not only directional thermal expansion has been inspected but also the role of serine residue and different repeating element on the structure has been assessed under thermal strain.
In general, the type of interaction determines the character of thermal expansion in the corresponding directions. Maximum expansion is less than 1% along the h-bond and chain directions. Most of the volumetric expansion comes from the sheet stacking direction as a result of weak van der Walls interaction. Simultaneous operation of multitude of hydrogen bonding and the covalent bonding give rise a cooperative effect resulting  with high thermal strength and stability along the intra-sheet direction and the chain direction for both model, respectively.
A nonconventional behavior is observed along the chain direction of crystalline domain with serine that it is contracts about 0.5% between 300 K and 420 K then it stays stable. Note that inclusion of serine residue in the model crystalline domains creates the hydrogen bonding constructed by the serine amino acids between the neighboring sheets of (GAS) 3:2:1 hexapeptides. Thermally activated torsional motion of rich side chain (HG1) of serine residue increases the dimensionality of the hydrogen bonding between the sheets, that causes the contraction along the chain direction. The types of h-bond between β-sheets of (GAS) 3 and are shown in figure 6. The estimated linear TECs are compared the experimental data of x-ray diffraction of along the h-bond and sheet stacking directions for β-sheet crystal unit of spider silk (BSS)) [56], α form of silk fibroin (ASF), β form of poly(γ-methyl D-glutamate (PMDG) and poly(L-alanine) (PLA) [57] which are also given in figure 6 as inset for the corresponding directions as well as in table 2. Table 2    This shows that the inclusion of a rich side chain like serine will give rubbery behavior to β-sheets crystallites as a result of thermally activated torsional motion of rich side chain (HG1) of serine residue by increasing the dimensionality of the hydrogen bonding between the sheets. Hence, this causes the contraction of β-sheets crystallites along the chain direction, which will also contribute to the contraction of the silk fibers together with amorphous regions. The negative thermal expansion materials are of particular importance in designing controlled thermal expansion composite materials.

Heat capacity C P and energetics of the crystallite
In order to investigate thermal stability thoroughly, we have also look at the behavior of the heat capacity C , P rate of change of entropy Energetically, serine residue has less than 2 kJ mol −1 contribution per serine up to about 450 K suggesting that the main chain in the β-sheet crystal unit is the determining factor and the richness of the side chain of the  serine has almost no effect on energetic of the system as the extension of the serine residue has rather small magnitude of thermal vibrations and torsions. Then, it starts increasing with temperature indicating that the side chains become the determining factor. Talking hypothetically, although it is not possible to observe denaturation process as we don't use the reactive potential, the main chain becomes the determining factor after about 600 K which corresponds to the experimental denaturation temperature of the crystal structure [38][39][40].

Conclusions
We have assessed the role of the serine residue and two different repeating elements in thermophysical properties of polar-antiparallel β-sheet crystalline domains of B.mori silk fibroin via molecular dynamics simulations between 300 K and 700 K. As expected, thermal expansion is anisotropic, and the type of interactions determines the equilibrated size of the domains and the character of thermal expansion in the corresponding directions. Inclusion of serine amino acid results with the contraction in the crystalline domains along the chain direction between 300 K and 420 K as a result of excess transverse hydrogen bonding between the sheets. That gives rise to a nonconventional 'negative' thermal expansion coefficient a » ---10 K y 4 1 at 300 K. This is of particular importance in designing controlled thermal expansion composite materials.
The heat capacity increases sublinearly up to about 650 K leading to almost continuous decrease in the rate of change of entropy. Although pronounced thermal vibrations and torsions causes comparable large variation in C P at about 450 K, there is no any indication of structural transformation. At low temperature region below 450 K, the main chain in the β-sheet crystal unit is the determining factor; the richness of the side chains has no effect on energetic of the system. Then the side chain difference becomes evident with a linear increase in enthalpy change per excess residue up to 600 K. Talking hypothetically, the energetic of excess side chain is dissolved during the denaturation of the crystallites and, again, the main chain is the determining factor.