Phenomenological models of Bombyx mori silk fibroin and their mechanical behavior using molecular dynamics simulations

Materials Science & Engineering C 108 (2020) 110414

Contents lists available at ScienceDirect

Materials Science & Engineering C journal homepage: www.elsevier.com/locate/msec

Phenomenological models of Bombyx mori silk fibroin and their mechanical behavior using molecular dynamics simulations

T

⁎

Mrinal Patel, Devendra K. Dubey , Satinder Paul Singh Department of Mechanical Engineering, Indian Institute of Technology Delhi, New Delhi 110016, India

A R T I C LE I N FO

A B S T R A C T

Keywords: Bombyx mori silk fibroin Biomaterial Molecular dynamics Computational materials science Nanoscale mechanics β-Sheet crystallite

Bombyx mori silk fibroin (B. mori SF) is a promising biopolymer for use in biomedical applications such as tissue engineered grafts and as a load bearing biopolymer with biocompatible and bioresorbable properties. B. mori SF is a hierarchical bio- macro-molecule made up of amino acid (residue) chains consisting of a crystalline phase and an amorphous phase arranged in a specific order. Understanding about the mechanical behavior of B. mori SF at multiple length scales is of importance when developing tissue grafts, which requires a deeper understanding of the mechanics of its nanostructure. Four phenomenological models of B. mori SF nanostructures were developed, based on crystalline and amorphous phase connectivity. Tensile loading based mechanical behavior analysis of these models were performed using molecular dynamics (MD) simulations and compared with existing results from literature for selection of best performing model. Elastic modulus of ~7.4GPa and tensile strength of ~340 MPa were obtained for this model. Analysis of results reveals that deformation mechanisms in B. mori SF at nanoscale are a combination of tensile and shear deformations, wherein, the tensile deformation of amorphous region results into excessive extension of B. mori SF, whereas, shear deformation of crystalline region results into a high tensile strength. Overall, this work is instrumental in development of a right computational nanoscale model of SF nanostructure and provides deeper insights into the mechanistic interactions and mechanisms between amorphous and crystalline regions of B. mori SF, which would be useful for further studies of silk based biomaterials.

1. Introduction Silk fiber is a natural protein polymer, finding extensive application in the textile industry. Apart from this, silk has been used as biomaterial for surgical sutures for centuries [1]. The application of silk as suture material has been possible due to its biocompatible and biodegradable nature [2]. Along with properties like biocompatibility and controllable biodegradability silk also has unique mechanical properties with a combination of strength and extensibility suitable for scaffolds/grafts in biomedical applications [3]. This collection of properties has motivated the study of silk as biomaterial for tissue engineered grafts. Silks are spun into fiber by Arthropods and Lepidoptera larvae such as silkworms (e.g. Bombyx mori), spiders, scorpions, mites, fleas [4,5]. The composition (amino acid residue sequence) and molecular structure of silk varies with the source leading to variation in mechanical properties as well [2]. For example, spider dragline silk produced by spiders (Nephila clavipes) is one of the strongest biomaterial with high strength, extensibility and toughness, even higher than some of the engineering material like steel [6]. B. mori silk also has high tensile

⁎

strength (300–800 MPa), elastic modulus (7–18 GPa), toughness (6 × 104 J/kg), and extensibility (15–19%) [7,8]. Along with having suitable mechanical properties for biomedical application and biocompatibility, it is also an easily available material as it is produced by domesticated silkworm. 1.1. Hierarchical structure of silk The mechanical and other physical properties of materials are governed by its structure at atomic level, i.e., the nanostructure of the material. The knowledge about the nanostructure of materials is important to understand the source and mechanism of their mechanical and other functional properties [9]. In proteins, the primary structure or the amino acid sequence governs the molecular structure and arrangement of the peptide chains, i.e., the secondary structure. Further, the secondary structure of the proteins governs its mechanical and other physical properties [10,11]. At macroscale both B. mori silk fiber and spider dragline silk fiber have similar structure, i.e., both have core-shell structure. In B. mori silk fiber two parallel fibroin (fibrous

Corresponding author. E-mail address: [email protected] (D.K. Dubey).

https://doi.org/10.1016/j.msec.2019.110414 Received 2 March 2019; Received in revised form 31 October 2019; Accepted 7 November 2019 Available online 09 November 2019 0928-4931/ © 2019 Elsevier B.V. All rights reserved.

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Fig. 1. Hierarchical structure of B. mori Silk Fiber. (a) B. mori silk worm cocoon. (b) B. mori silk fiber consists of two fibroin fibers forming the core coated with sericin protein holding the core fibers together. (c) The fibroin fibers consist of numerous nano-fibrils. (d) β-Sheet crystallite of size around 20X20X20 Å and (e) Amorphous chain.

conformation forming the parallel and anti- parallel β-sheet crystallites, respectively [16–18]. In B. mori silks the fractional ratio of antiparallel to parallel β-sheet crystallites is ~2:1 [19]. As illustrated in Fig. 1(d), the repetitive domains from different molecular chains align/arrange together to form the β-sheet crystallites (intermolecular) and the nonrepetitive domains form the amorphous region. This combination of various molecular chains in a nano-structural region gives rise to an intertwined molecular network. The (Alanine)n repeats in spider dragline silk is suggested to have higher propensity for interlocking the adjacent chain via hydrogen bonds and form β-sheet crystallites which in turn is one of the factors responsible for difference in mechanical behavior of B. mori silk and spider dragline silk [20–22]. Due to this unique fundamental units (primary structure) of B. mori silk it essential to study it from the building block level. The atomistic model of spider dragline silk nanostructure involving crystalline and semi-amorphous regions has been developed [9] and investigated using computational techniques to understand the molecular mechanism of deformation under mechanical loading. Also for B. mori SF, a representative atomistic model (1SLK and 2SLK) for β-sheet crystallites has been reported [23]. But studying only the crystalline region will not provide an exact interpretation of the mechanical behavior of B. mori SF due to the unknown factor of amorphous region and connectivity between crystalline and amorphous regions. To study and establish an understanding of the molecular mechanism of deformation under mechanical loading and the structure-property correlations in B. mori SF, there is a need to develop an atomistic model of B. mori SF nanostructure which includes both crystalline and amorphous regions. The current work focuses on developing a computational model for B. mori SF and understands the molecular mechanics in B. mori SF nanostructure under tensile deformation.

protein) fibers form the core, coated with sericin (glue-like protein) forming the shell, as shown in Fig. 1(b), whereas, the spider dragline silk fiber consists of protein monofilament of spidroin (containing Major Ampullate Spidroin 1 and 2) forming the core, coated with MAspidroin like protein, glycoprotein and lipids forming different layers of the shell [12–14]. These spidroin and fibroin cores are composed of several nano-fibrils with diameter of nano-meter length scale. Fig. 1b shows the schematic representation of B. mori silk fiber with two fibroin fibers coated with sericin protein where the fibroin fiber is composed of numerous nano-fibrils, respectively. The primary and secondary structures of silk such as amino acid residue sequence, size of crystallites, packing of peptide chains in crystal and alignment of peptide chains in the molecule have been reported in earlier experimental studies describing structural details of silk at multiple hierarchical levels [15,16]. For example, the secondary structures of B. mori SF and spider spidroins have been found to be similar in terms that both the silks have large repetitive domains flanked by non-repetitive domains. These repetitive domains form the crystalline region and the non-repetitive domains form the amorphous or semi-amorphous region in both B. mori SF and spider spidroins, respectively [4,12,13]. Fig. 1c shows the schematic representation of secondary structure of the molecules in the B. mori SF nano-fibril, where β-sheet crystallites, shown in Fig. 1(d), forming the crystalline region is embedded in a matrix formed by amorphous chains, shown in Fig. 1(e). Along with the various similarities in the secondary structure of B. mori SF and spider spidroins, there is a major fundamental dissimilarity in their primary structure. The amino acid residue sequence of the B. mori SF and spider spidroins has been found to be very different from each other. For example, in B. mori SF, the repetitive domain is majorly composed of (Glycine-Alanine)n repeats whereas, in spider spidroins the repetitive domain mainly consists of (Alanine)n repeats [12]. These repetitive domains are responsible for forming a well-defined crystalline region. Both the silkworm silk and the spider dragline silk contain the repetitive domain arranged in parallel and anti- parallel 2

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

and leading to great extensibility, whereas, the crystalline region governs the mechanical behavior at large deformation and the tensile strength of spider silk [6]. The force distribution analysis performed in an earlier study at multiple scales using computational model of spider silk at nano-scale and continuum scale (using FE) showed that the mechanical behavior of spider silk is highly dependent on the specific size and relative arrangement of the crystalline and amorphous subunit. [49] Various studies have been performed on native B. mori SF and regenerated SF fibers, their structure–property relationship, enhancement of mechanical properties (in fiber, film and hydrogel form), behavior of regenerated SF at air–water interface, the glass transitions of naturally spun silk fibers using dynamic mechanical thermal analysis [27,29,38,42–45]. From these studies, it was seen that the fracture stress and strain of regenerated silk fiber depend on the drawing conditions due to the structural changes induced in the regenerated fibers, such as increase in β-sheet content, crystallinity and molecular orientation. It was suggested that the evolution of structure during the spinning process determines the mechanical performance of silk fibers [27]. The elastic modulus of the pre-stretched fiber was found to be increasing in the direction of stretching due to increase in molecular chain orientation with increased draw ratio [43] and also the conformation transition from random coil to β-sheet in alcohol environment produced the regenerated SF with small-sized and uniformly distributed crystallites throughout the sample having excellent strength and extensibility [38]. The chemical environment and the spinning parameters such as spinning speed, draw-down ratio during the spinning process of regenerated SF was shown to be governing the mechanical behavior of silk fibers by governing the nano-structure of fiber. Investigations focusing on the structure – mechanical property relationship of B. mori SF, the dynamic mechanical properties, effect of macroscopic deformation on molecular structure of B. mori SF and the conformational changes induced in regenerated B. mori SF molecular structure upon treatment with ethanol has also been performed [24,28,37]. The results from these investigations suggested that the initial deformation originates from the chain stretching and sliding in the amorphous regions and as the chains are extended further, the crystalline regions also start to take up the load and eventually when the load is too high to be accommodated, the fiber fails [37]. Under mechanical stress silks tend to dissipate strain energy by transforming the secondary structure [28]. The yield strain of around 3%, tensile strength of 400 MPa and failure strain of 22% were also obtained. Forced reeling of silk from immobilized silkworms under steady and controlled conditions was performed and it was observed that fibers produced are superior than naturally spun ones in terms of mechanical behavior. At faster spinning speeds, stronger (but more brittle) fibers were obtained, whereas at slower speeds weaker, more extensible fibers were produced. It was suggested that silkworm silks with mechanical properties comparable to spider silks can be produced by adjusting the harvesting parameters [8]. Role of different molecular segments (crystalline and amorphous) in governing the overall mechanical properties and behavior of B. mori SF was also investigated by adapting a linear viscoelastic model to the semi-crystalline morphology of silk. It was suggested that the amorphous segment plays an important role in mechanical behavior of B. mori SF. The high extensibility of the silk was linked majorly with the disordered phase with little contribution from the crystallites also [31,36]. Studies focused on modulating the mechanical properties of B. mori SF have also been performed by controlling the extent of crystallinity (e.g. using methanol), by varying the spinning process of the naturally spun silk (e.g. changing the environmental condition like temperature and humidity or force-spinning at a controlled drawing rate), by controlling the morphology (fibroin sponges, hydrogel, films, mats and blended/cross-linked fibroin) of regenerated silk [3,26,32,33,39]. A considerable amount of experimental studies have been performed to understand the mechanical behavior of B. mori SF along with

1.2. Studies on mechanical behavior of B. mori SF In past few decades the mechanical properties and behavior of B. mori SF have been extensively studied using experimental and computational techniques [8,24–45]. The computational studies focused on analyzing the mechanical behavior of β-sheet crystallites of B. mori SF using MD simulations. The results from these studies showed that the βsheet crystallites fail by shearing of different sheets forming the crystallite, wherein the β-strands are pulled out of the crystallite in stickslip motion. It was also suggested that both intra-sheet H-bonding and inter-sheet side-chain interactions contribute equally to the strength of β-sheet crystallite. These investigations showed that the antiparallel βsheet crystallite models were found to be stronger and stiffer than the parallel one [40]. The rupture force of the β-sheet crystallites obtained in different studies is in the range of 2000–3600 pN. It was also observed that the rupture force of β-sheet crystallites decreases with loading rate [41]. In these studies, only the crystalline region has been investigated but to analyse the mechanical behavior of B. mori SF, it is necessary to consider the crystalline and amorphous regions both. The role of H-bonds in deformation behavior of β-sheet crystallites under shear loading and the size effect of β-sheet crystallites on the overall mechanical properties of B. mori SF has been investigated to understand the deformation mechanism and the contribution of weak H-bonds in the strength of the silk fibers. The results showed that the size of the β-sheet nanocrystals must be confined to a few nanometers to ensure the maximum cooperative deformation of H-bonds and it was illustrated that the β-sheet crystallites of 2 to 4 nm size range exhibit highest stiffness, toughness and strength [30,34]. Also a study to develop a method for designing the silk fiber with desired mechanical properties through a proper amino acid sequencing and processing parameters using multiscale simulations and genetic synthesis of protein sequence has been performed. The protein network was analyzed using computational methods and its effect on mechanical properties of fiber. The results were validated experimentally using the synthetically fabricated protein with genetic design. The observation from this study suggested a 1 to 1 ratio of hydrophobic to hydrophilic blocks for a connected polymer network where longer sequences were found to be stronger and having higher Young's modulus before shear flow, whereas, after the shear flow, maximum stress and tensile strain were found to be much higher for shorter sequences [35]. The mechanical behavior of spider silk has also been thoroughly investigated at multiple length scales to understand the nano-structural and molecular mechanics of deformation of spider silk and its sequencestructure-properties correlation. Computational model of spider silk at nano-scale using REMD (Replica Exchange MD) was developed and at continuum scale using FE was also developed and different loading scenarios of silk crystallites were investigated. The size effect of constituent building blocks on the nanostructure and mechanical behavior of Silk fiber such as effect of poly-Alanine strand length on the behavior of β-sheet crystallites was also studied [6,9,14,20,22,46–51]. The results from these studies showed that the poly-Alanine domains form βsheet crystallites, and Glycine-rich domains form β-turns and 31 helix type structures and the size-effect study showed that the poly-Alanine strand of minimum 6–8 residues has a well-defined and clearly visible β-sheet crystallites and β-sheet crystallites with poly-Alanine strand above 8 residues shows no significant improvement in strength and structure. Also the crystallite with cross-sectional are of around 1 nm2 shows highest fiber toughness [20,46,47,50]. A coarse-grained mesoscale modelling of spider silk nanostructure was also performed to study the molecular and nano-structural mechanism of deformation under tensile loading. It was showed that βsheet crystals geometrically confined to few nano-meters (3 nm) combined with semi-amorphous region (highly extensible), is important to achieve unique mechanical properties of spider silk and the semiamorphous region plays a significant role in defining the mechanical behavior of spider silk at small deformations by uncoiling under load 3

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

The repetitive domain is majorly made up of Glycine-X repeats (94% of repetitive domain), where X is Alanine in 64%, Serine in 22%, Tyrosine in 10%, Valine in 3%, and Threonine in 1.3% of the repeats [11] where, Glycine-X repeats act as building block for the crystalline region of B. mori SF. Two hexapeptides GAGAGS and GAGAGY make 70% of the crystalline region with 433 and 120 copies, respectively. Also the GAGAGS hexapeptide forms around 56% of the crystal forming domain [11]. Many of the previous computational studies with a focus on analyzing the crystalline region of B. mori SF studied a β-sheet crystallite comprising of β-strands formed by hexapeptide GAGAGS [25,30,40,41,58]. The β-sheet crystallites present in B. mori SF can be either intermolecular or intramolecular, but since the crystal forming domains are predominantly long (~300–400 residues), the likelihood of forming intermolecular β-sheet crystallites is higher than the intramolecular β-sheet crystallites (18%) [3,59]. X-ray diffraction and Solid-state NMR studies have shown the presence of both parallel or anti-parallel β-sheet crystallites in a ratio of ~1:2 in B. mori SF [3,16,19]. Also, the computational studies have shown that (also mentioned in Section 1.2) the antiparallel β-sheet crystallite models are stronger and stiffer than the parallel ones [40]. At molecular level, the amorphous chains and the peptide chain forming the β-sheet crystallites are connected through covalent bonds [60]. Further, in an earlier study, it has been observed that β-sheet crystals of few nanometers size achieve high stiffness, strength and toughness as compared to larger nanocrystals and the optimum size of nanocrystals was found to be 1–2 nm along the fiber axis direction and 2–4 nm in direction perpendicular to fiber axis and along the sheet [30]. The overall mechanical properties of silk are governed by repetitive domains that play a key role in providing stiff domains embedded in soft and extensible amorphous region in B. mori SF.

some computational studies to understand the nanoscale mechanics and mechanical behavior of crystalline region of B. mori SF. However, there is lack of a detailed atomistic computational model of B. mori SF as a semi-crystalline composite (i.e. with both crystalline and amorphous phases). Also, B. mori SF has similar secondary structure as the spider spidroins (in terms that both the silks have large repetitive domains and non-repetitive domains), some insights can be drawn from the computational studies focused on understanding the structure-mechanical behavior correlation in spider silk. However, the correlation between the secondary structure, such as, the relative alignment of the crystalline and amorphous regions in the B. mori SF fiber at nanoscale along with the interactions/connectivity present between the crystalline and amorphous regions and the role of two fundamental components in mechanical behavior of B. mori SF, still remains obscure for B. mori SF. This understanding about the mechanical behavior of B. mori SF at multiple length scales is vital to aid the ongoing research in the field of developing scaffolds and tissue grafts for biomedical application. The applications of B. mori SF in biomedical field may require its use in different forms (e.g. as a matrix material in bio-nano-composites) or morphologies (e.g. hydrogel, film, sponge, 3D scaffolds) and certain functionalities need to be introduced in it at the same time conserving its intrinsic properties. The imminent accessibility of computational model of B. mori SF nanostructure and insights into the mechanisms of deformation can be instrumental for tailor-made material designing with mechanical properties conforming with those of the respective tissues (both hard and soft tissues). To facilitate this, investigations need to be performed at nanometer length-scale. Molecular dynamics method is one of the promising tools to study such phenomenon at nanoscale. For this purpose, an atomistic model comprising of both crystalline and amorphous regions needs to be developed. In this study, four computational atomistic models of B. mori SF have been proposed. The models consist of both crystalline and amorphous regions. Different configurations/connectivity between amorphous chains and crystalline units based on density of B. mori SF and crystalline to amorphous ratio have been considered. The mechanical properties and behavior under tension have been determined and compared with the previously reported results in literature.

2.2. Atomistic models of B. mori SF The complete atomistic model of B. mori SF, although not fully known yet, would be computationally very expensive and infeasible to study. Therefore, a representative model of β-sheet crystallite and amorphous chains reported in earlier studies is considered here [11,15,23,61].

2. Method and framework 2.2.1. Modelling the crystalline and amorphous regions The crystalline region of B. mori SF is majorly made up of the hexapeptide GAGAGS and is shown to be the building block of the crystalline region in earlier studies [25,30,40,41,58]. The coordinates of β-sheet crystallite used to build the models have been obtained from Protein Databank (PDB Id - 1slk) [23,62]. The last residue of each strand has been mutated from alanine to serine using Visual Molecular Dynamics software (VMD) [63]. This is done to obtain the GAGAGS sequence. Also extra two layers are added to get a 5 layer β-sheet crystallite, as shown in Fig. 3(a), to achieve the nanocrystal size in the range of 2 to 4 nm as described earlier in Section 2.1. Hence each βsheet crystallite consists of 5 β-sheets arranged one top of another at 5.6 Å distance and each β-sheet, as shown in Fig. 3(b), consists of 5 antiparallel β-strands at 4.8 Å distance, i.e., a total of 25 β-strands. Each βstrand, as shown in Fig. 3(c), is a hexapeptide made of GAGAGS sequence, length 23.6 Å. To model the amorphous domain, a segment of 42 residues (linker 6 from earlier study [11]), GAGAGAGAGAGTGSSGFGPYVAHGGYSGYEYAWSSESDFGTGS has been considered. Although, linker 6 (from earlier study [11]) has been chosen to model the amorphous chains in the current study, the residue sequence of all the eleven linkers [11] are mostly same with hardly a difference of a few residues. The amorphous chain has been shown in Fig. 3(d). The amorphous chain with the above mentioned sequence has been modeled using PEPFOLD [64–66]. Free modelling feature of PEPFOLD has been used to model the amorphous chain in neutral pH (pH 7) environment. PEPFOLD uses the structural alphabet concept, a

In order to develop an atomistic level mechanical behavior understanding of B. mori SF, molecular dynamics method is used as it is a potential tool to study nanoscale material behavior [52–54]. The computational methodology followed here has three steps, (i) atomistic modelling of the two components (crystalline and amorphous regions) of B. mori SF, (ii) creating four models of B. mori SF based on connectivity between crystalline and amorphous regions and (iii) MD simulations for tensile deformation test of the resulting structures followed by analysis of results. 2.1. Structure of B. mori SF The B. mori SF macro-molecule comprises of three segments called as Heavy Chain (~350 kDa), Light Chain (26 kDa) and P25 gene (25 kDa) in ratio of 6:6:1 [11,15]. The Heavy Chain and the Light Chain are linked by a single disulfide link while the P25 gene has non-covalent interactions with the Heavy Chain and the Light Chain [55]. Heavy chain is made up of 5263 residues where Glycine (G) is present in 45.9%, Alanine (A) in 30.3%, Serine (S) in 12.1%, Tyrosine (T) in 5.3%, Valine (V) in 1.8%, and 4.7% of the other amino acids [11]. The heavy chain has twelve repetitive domains which are Glycine rich (~45%) forming the crystalline region separated by eleven short linker domains (42–43 residues). These short linkers are non-repetitive in nature and form amorphous region [11,56]. Fig. 2 shows a schematic representation for organization of amino acid residues in the primary structure of B. mori SF heavy chain. 4

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Fig. 2. Schematic representation of fine organization of the primary structure of B. mori SF heavy chain, showing the 12 repetitive and 11 amorphous regions, respectively. An approximate amino acid sequence of the R10 region is illustrated by modules i, ii, and iii. Module iv depicts a typical amino acid sequence in the amorphous region. (Reprinted with permission from ref. [57]. Copyright 2005 American Chemical Society).

states of prototype fragments using the sOPEP coarse grained representation (one bead per side chain). The peptide chain is built residue by residue. Once the peptide chain is complete, 30,000 steps of Monte-Carlo simulation are used for optimization [64].

probabilistic framework for generalizing the concept of secondary structure. While modelling in PEPFOLD, peptides are divided into fragments of four amino acids. For each such fragment, geometric descriptors are calculated and a series of states are identified that best describe the conformation of peptide chain, using standard Hidden Markov Models algorithms such as the Viterbi algorithm or the Forward Backward algorithm. One conformation is generated from each series of

2.2.2. Computational model of B. mori SF nanostructure To develop a computational model of B. mori SF nanostructure, the

19.2Å

Glycine (G) Serine (S) 22.4Å

Serine (S) Glycine (G)

Amorphous chain (d) β-sheet crystallite (a)

Threonine (T) Glycine (G) Carbon Nitrogen Alanine (A) Oxygen Hydrogen

Single β-sheet (b)

β-strand (c)

Serine (S)

Fig. 3. Building blocks of SF molecule. (a) β-Sheet crystallite consisting of 5 β-sheets. (b) Single β-sheet consisting of five β – strands (hexapeptides). (c) All atom structure of one hexapeptide GAGAGS (forming the β – strand). The residues displayed in blue, green, yellow represent Glycine, Alanine and Serine residues, respectively. (d) Amorphous chain used for creating SF model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 5

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

interactions only, representing the nano-structural region where the neighboring hexapeptide GAGAGS and the amorphous chains are from different molecules. Fig. 4(b) shows model II-CA having Amorphous-L chains covalently connected with β-strands in crystallites and no covalent connectivity between the three Amorphous-P chains. This models the nano-structural region where the neighboring hexapeptide GAGAGS and the Amorphous-L chains are from the same molecules which is then surrounded by a matrix of amorphous chains (Amorphous-P chains) from different molecules. In model II-CA, the interface between AmorphousL chains and β-sheet crystallites have both bonded and non-bonded interactions. Fig. 4(c) shows model III-AA, it has no covalent connectivity between Amorphous-L chains and β-strands in crystallites. The Amorphous-P chains are covalently bonded together thus forming one continuous chain. Hence, in model III-AA, the interface between Amorphous-L chains and β-sheet crystallites have only non-bonded interactions. Fig. 4(d) shows model IV-AACA, having Amorphous-L chains covalently bonded with β-strands in crystallites and the Amorphous-P chains are also covalently bonded together forming one continuous chain to model the nano-structural region where the neighboring hexapeptide GAGAGS and the Amorphous-L chains are from the same molecules. As a result of this connectivity in model IV-AACA, the interface between Amorphous-L chains and β-sheet crystallites have both bonded and non-bonded interactions between them. The covalent bonds between Amorphous-L chains and the β-strands in crystallite and in-between the three Amorphous-P chains, in model IV-AACA, III-AA, II-CA have been created using BIOVIA Materials Studio.

required combination of crystalline and amorphous regions is to be unified taking care of certain parameters, such as ratio of crystalline and amorphous regions, density of B. mori SF and the mass and size of the two regions, along with the residue sequence, size and composition of crystalline and amorphous regions. The density of B. mori SF has been reported to be around 1.3 g/cc, with the amount of crystalline to amorphous ratio to be around 1:1 [67]. Based on these information, it is calculated that for each β-sheet crystallite made of 25 GAGAGS strands, there should be around 4 amorphous chains with the mentioned 42 residues to maintain the ratio of crystalline and amorphous chains to 1:1 and 1.3 g/cc density. The four idealized 3-D atomistic computational models of B. mori SF at nanoscale based on connectivity between crystalline and amorphous regions have been conceptualized. The extent of connectivity between the crystalline and amorphous domains at any nano-structural region is not known. The connectivity between the crystallites and amorphous chains has been devised to capture the inter-atomic interactions between different neighboring regions, i.e., the hexapeptide GAGAGS (in β-sheet crystallite) and amorphous chains in the nano-structural region are part of the same molecule or different molecules. If the hexapeptide GAGAGS and the amorphous chain are from the same molecule, both bonded and non-bonded interaction will exist between them, whereas, if the neighboring amorphous chain is not the part of same molecule as the hexapeptide GAGAGS, only non-bonded interactions will exist between them. The crystallites are aligned along the fiber axis direction separated with amorphous chains in all the three directions to replicate a system where β-sheet crystallites are arranged in amorphous region. Such a model is necessary to capture the effect of both crystalline and amorphous regions and interactions between them, on deformation mechanism of B. mori SF. As shown in Fig. 4, the amorphous chains represented in red (denoted as ‘Amorphous-L’), provides separation and connectivity between two crystallites (represented in blue) along the fiber axis direction and the amorphous chains represented in green (denoted as ‘Amorphous-P’), provide separation between the adjacent crystallites in other two directions, i.e., perpendicular to fiber axis to ascertain the crystallite size of 23.6 Å × 22.4 Å × 19.2 Å. All the models consist of 8 crystallites separated by amorphous regions. The interface between the two regions has the covalent connectivity in two of the models along with the non-bonded interaction (in all the models) between the peptide chains in crystalline and amorphous regions. The models differ in connectivity between different crystallites and amorphous chains and are named as follows based on connectivity:

2.3. Simulation details The MD simulations in this work have been performed using LAMMPS [68], an open source MD code by Sandia National Laboratories. CHARMM27 force field [69–71] is used for simulation as it has been used extensively in earlier studies on proteins like collagen and silk [20,22,25,30,41,46,72,73] and shown to correctly model their inter-atomic interactions. The CHARMM27 potential energy function comprises of intramolecular energy (Uintra) and intermolecular energy (Uinter), given by:

U = Uintra + Uinter The intramolecular interaction potential is expressed as:

I. No interconnection between crystalline and/or amorphous chains (I-NC) II. Crystalline-Amorphous interconnection (II-CA) III. Amorphous-Amorphous interconnection (III-AA) IV. Both Amorphous-Amorphous and Crystalline-Amorphous interconnections (IV-AACA)

Uintra =

∑

Kb (bij − b0 )2 +

bonds

(1 + cos(nχ − δ )) +

∑

K θ (θijk − θ0 )2 +

angle

∑ impropers

∑

Kχ

dihedrals

Kimp (φ − φ0 )2 +

∑ KUB (Sik − S0 )2 UB

where the first term in the energy function represents bond stretching energy, second term accounts for bond angle interactions, third term is for dihedrals (or torsion angles) and fourth term is for improper, i.e., out of plane bending. The last term is Urey-Bradley component, that is, cross-term accounting for angle bending using 1, 3 non-bonded interactions. Here, Kb, Kθ, Kχ, Kimp and KUB are force constants, b0, θ0, φ0 and S0 are equilibrium bond length, equilibrium bond angle, equilibrium improper torsion angle and equilibrium distance between atoms separated by two covalent bonds (1,3 distance) respectively. bij is the distance between atoms i and j, θijk is bond angle between sites i and j, k, and χ and φ are the dihedral angle and improper torsion angle and Sik is the distance between atoms i and k, respectively [70]. The Intermolecular potential involves Van der Waals interactions and electrostatic interactions, and is expressed as:

Fig. 4a–d shows the 2D schematic representation of the four hypothesized phenomenological atomistic computational models of B. mori SF, where blue squares represent β-sheet crystallites, the red circles represent the Amorphous-L chains providing connectivity between the crystallites along the fiber axis direction and the green circles represent Amorphous-P chains separating the adjacent crystallites in direction perpendicular to the fiber axis. In all the four models the interface between Amorphous-P chains and β-sheet crystallites have nonbonded interactions only. Fig. 4(a) shows model I-NC. It has no covalent connectivity between Amorphous-L chains the β-strands in the crystallites and also in between Amorphous-P chains. As a result, the interface between Amorphous-L chains and β-sheet crystallites and the interface between Amorphous-P chains and β-sheet crystallites have non-bonded 6

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Fig. 4. 2D schematic representation the four hypothesized phenomenological models of B. mori SF nanostructure (a) I-NC; (b) II-CA; (c) III-AA; (d) IV-AACA. 12

⎡ Rminij ⎞ ⎜ ⎟ r ⎣ ⎝ ij ⎠

Uinter =  εij ⎢ ⎛

6

Rminij ⎞ ⎤ − ⎜⎛ ⎟ ⎥ + ⎝ rij ⎠ ⎦

∑

direction (z-axis) as shown in Fig. 6. The tensile deformation here is carried-out by imposing a uniform strain on the simulation cell at any point of time. Also, the strain rates used in earlier studies focusing on understanding the molecular mechanism of polymer systems (both amorphous and crystalline), are of comparable order to the strain rates used in current work. [74–77]. The strain rates used here are typical for MD simulations in terms of, first, the length-scale of the system under study, and second, the localization of the simulation cell to a particular section of the macro-scale body under study. The MD simulations for tensile deformation are performed in NVT ensemble at the 300 K until failure. The stress calculated here is Virial Stress, which is established from the Clausius Virial Theorem along with Maxwell velocity distribution [78] and is given by:

qi qj erij

where the Van der Waals interaction is modeled using a Lennard-Jones potential and the electrostatic interaction with Columbic potential. qi and qj are the partial atomic charges, εij is the Lennard–Jones (LJ) well depth, Rminij is minimum interaction radius [70]. Each B. mori SF simulation cell (atomistic model) is simulated under periodic boundary conditions (PBCs) in the direction of fiber axis (zaxis) and non-periodic for the other two directions (x and y axis). The initial atomic velocities are assigned from Gaussian distributions corresponding to the desired temperature, i.e. 300 K. The simulation cells are first allowed to undergo energy minimization using Conjugate Gradient Method for 10 ps. Thereafter dynamics is performed and Velocity-Verlet scheme is used to integrate Newton's equations of motion with time step of 0.5 fs. The cutoff distance for non-bonded interaction is set to 13 Å. Before deformation, the models are first equilibrated in NPT ensemble at 300 K and 1 bar pressure for 5 ns. Constant temperature is maintained by a Nose−Hoover thermostat and constant pressure is maintained using a Nose−Hoover barostat. The stability of the models is verified from the observation of total energy and volume of the system. It is observed that total energy and volume of the models have been stabilized, the same has been shown in Fig. 5, in this time with minimal pressure fluctuations. Further, the models are equilibrated in NVT ensemble at 300 K until the total energy and pressure of the system have been stabilized. Tensile deformation tests are then performed at three different strain rates, i.e., 5 × 10−4, 5 × 10−5, 5 × 10−6 per fs in fiber axis

σij =

1 V

⎡1

N

⎤

α α α ∑ ⎢ 2 ∑ (rβi − r αi )F αβ j + m vi v j ⎥ α

⎣

β=1

⎦

where σij calculates stress experienced by atom α due to neighboring atoms β (β = 1 to N), i and j represent x, y and z direction. V is volume of the system, riα and riβ are position of atom α and β in direction i. Fjαβ is force experienced by atom α due to atom β in direction j, mα is mass of atom α, viα and vjα are velocities of atom α in direction i and j [79]. 3. Results and analysis Earlier molecular dynamic simulation studies have focused on understanding the mechanical behavior of B. mori SF using the crystalline region only [25,30,40,41]. This gives information only about the 7

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Y

~71 Å

~71 Å

Fig. 5. (a) Energy vs. Time plot for equilibration simulation of model IV-AACA. During the 5 ns NPT equilibration simulation, the total energy of the system is observed to decrease and then become constant showing a stable configuration of the system. (b and c) The volume/size of the simulation cell is also observed to decrease from 69 × 69 × 98 Å to 71 × 71 × 75 Å and thereafter becomes constant.

≡

X Z

~75 Å

~75 Å

Fiber Axis Direction – Direction of applied tensile deformation Fig. 6. Three dimensional schematic of the model IV-AACA along with the atomic model. Tensile deformation is applied along the fiber axis direction which, in this case, is aligned along the Z-axis.

crystalline and amorphous chains.

crystalline part and not the B. mori SF as a whole. Present study is focused on understanding the behavior of B. mori SF nanostructure under tensile deformation by creating a nanoscale computational atomistic model of B. mori SF including both crystalline and amorphous regions along with exploring the extent of connectivity between different

3.1. Stress-strain behavior of B. mori SF models MD simulation for tensile deformation using atomistic model 8

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Fig. 7. Stress-strain plots of the four models I-NC, II-CA, III-AA and IV-AACA at different strain rates, (a) 5 × 10−4 per fs, (b) 5 × 10−5 per fs, (c) 5 × 10−6 per fs, (d) Normal, shear and equivalent stress for model IV-AACA at 5 × 10−6 per fs strain rate. (e) Schematic diagram for stress-strain curve showing four distinct regimes.

and normal stress in z-direction is not very significant. The same has been shown in Fig. 7(d). The stress-strain curve (Fig. 7) for models IV-AACA, III-AA and II-CA show four distinct regimes named as initial elastic regime (ending with intrinsic yield point), strain softening regime, strain hardening regimes and fracture regime, whereas, the stress-strain plot for model I-NC shows initial elastic regime followed by strain softening regime and ending in fracture regime with no strain hardening regime at all the three strain rates. Yield point here in the four model range from [0.13–0.25], [0.21–0.25] and [0.6–0.7] strain values for strain rates 5 × 10−6, 5 × 10−5 and 5 × 10−4 per fs respectively. Fig. 8 shows the

including both crystalline and amorphous regions has been employed to determine the modulus, tensile strength and yield strain of B. mori SF as a whole and provide an insight about the molecular mechanics involved in deformation. Molecular dynamic simulations for tensile test are performed on all the four models at different strain rates of 5 × 10−6, 5 × 10−5 and 5 × 10−4 per fs. Fig. 7 shows the tensile stress vs. strain plot for four models at the different strain rates. Also the stress values reported in this work are component of virial stress in z-direction (fiber axis direction) for all the models. Based on the calculated stress values it is seen that the shear stresses obtained during tensile test are negligible due to which the difference between the equivalent Von-Mises Stress 9

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Undeformed Structures

At yield point

At the end of strain softening regime

At the start of fracture regime

IV -AACA

0 strain

0.13 strain

1.3 strain

0.5 strain

III-AA

0 strain

0.18 strain

0.6 strain

0.20 strain

0.75 strain

1 strain

II-CA

0 strain

0.8 strain

I-NC

0 strain

0.25 strain

0.5 strain

1 strain

Fig. 8. Snapshots of the models IV-AACA, III-AA, II-CA and I-NC at different strain levels during deformation at strain rate of 5 × 10−6 per fs.

rates (5 × 10−6, 5 × 10−5 and 5 × 10−4 per fs) followed by yield point. Analyzing the deformation simulation trajectory along with the stress-strain plot, it is observed that up to the yield point, the deformation is majorly concentrated in amorphous region only by uncoiling and sliding of chains, as highlighted in Fig. 9(b). Also upon analysis based on the number of H-bonds, a steep drop in number of Hbonds has been observed until yield strain. Fig. 10(a) shows the variation in number of H-bonds with deformation (strain) for model IVAACA during the tensile deformation at different strain rates (5 × 10−6 per fs, 5 × 10−5 per fs and 5 × 10−4 per fs). It is seen that, the number of H-bonds decreases rapidly in the range of strain corresponding to the elastic regime and thereafter remains fairly constant with minimal fluctuations. This rapid drop in number of H-bonds upto the yield strain shows that the non-bonded interaction plays a vital role in defining the initial elastic modulus of B. mori SF. Fig. 10(b) shows the evolution in potential energy and its different components (non-bonded energy, bond stretching energy, bond angle energy and dihedral energy) for model IV-AACA during tensile deformation at 5 × 10−6 per fs strain rate. The non-bonded energy here is

snapshots of all the four models in undeformed state, at yield point, at the end of strain softening regime and at the start of fracture regime. For strain rates 5 × 10−5 and 5 × 10−6 per fs, following the yield point a clear strain softening regime is observed before the strain hardening regime for models II-CA, III-AA and IV-AACA whereas, in case of 5 × 10−4 per fs strain rate, models IV-AACA and II-CA show initial linear elastic regime followed immediately strain hardening regime without any intermediate s train softening regime, and, model IIIAA shows a flat plateau like regime between the yield point and the strain hardening regime. Model I-NC shows only the initial elastic regime followed by yield point, strain softening regime and thereafter ending with fracture regime with no strain hardening regime. The four regimes observed in the stress-strain curves has been discussed further and described in association with the deformation trajectory of the four models in the following subsections. 3.2. Elastic regime All the four models show the initial elastic regime at the three strain 10

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Fig. 9. Snapshots from deformation trajectory of model IV-AACA at (a) 0 strain and (b) 0.13 strain (yield point) at 5 × 10−6 per fs strain rate, uncoiling/stretching of amorphous chains is observed as highlighted, with no noticeable deformation in crystalline domain.

amorphous chains in elastic regime. From all the above observations, the elastic regime here can be correlated with deformation in amorphous region and the initial elastic modulus of B. mori SF is a result of non-bonded interactions between the chains in amorphous region and with chains in crystalline region. After the elastic regime also the potential energy has a major contribution from non-bonded energy and the increase in non-bonded energy at a lower rate can be correlated with the observation that there is almost no change in number of H-bonds in model IV-AACA after the yield point, as shown in Fig. 10(a), i.e., after elastic regime the nonbonded energy has almost no contribution from H-bonds. Also, after the

a combination of van der Waal's interaction and electrostatic interaction. The potential energy of the model during tensile deformation is increasing sharply in the initial elastic regime. The potential is increasing further in strain softening and strain hardening regimes also but at a lower rate and becomes almost constant after around 1.5 strain, i.e., in the fracture regime. In the elastic regime, the bond-stretching, bond angle and dihedral energies are fairly constant and the sharp increase in the potential energy is due to the increase in the non-bonded energy. This sharp increase in non-bonded energy together with the steep drop in number of H-bonds up to the yield strain (as shown in Fig. 10a), can be correlated to the chain slippage and uncoiling of the

Energy vs. strain for model IV-AACA at 5x10-6 per fs strain rates

Number of H-bonds with increasing strain for model IV-AACA at different strain rates No. of H-bonds corresponding to yield strain

5X10-6 per fs

Energy at yield strain

5X10-5 per fs

5X10-4 per fs

(a)

(b)

Fig. 10. (a) Variation in number of H-bonds with increasing strain for model IV-AACA during the tensile deformation at different strain rates, 5 × 10−6 per fs, 5 × 10−5 per fs and 5 × 10−4 per fs. (b) Variation of energy with strain for model IV-AACA during the tensile deformation at 5 × 10−6 per fs strain rates. 11

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

(a) IV-AACA: Start of strain hardening regime (at 0.5 strain)

(b) III-AA: Start of strain hardening regime (at 0.75 strain)

Start of strain hardening regime (c) II-CA: Start of strain hardening regime (at 0.6 strain) (d) Stress Vs Strain: Strain rate 5X10-6 per fs Fig. 11. Snapshots of the models at the start of strain hardening regime (a) IV-AACA, (b) III-AA, (c) II-CA at strain rate of 5 × 10−6 per fs, respectively. (d) Stressstrain plot of the four models at 5 × 10−6 per fs strain rate demarcating the start of strain hardening regime.

rate of 5 × 10−6 per fs. It is observed that, in model IV-AACA, both the types of amorphous regions (Amorphous-L and Amorphous-P) undergo equivalent amount of deformations, leading to the inference that load is carried by both the types of amorphous regions equally. In models IIIAA, it is seen that Amorphous-P chains have deformed slightly more than the Amorphous-L chains whereas, in model II-CA, Amorphous-L chains are seen to be deforming and virtually no deformation in Amorphous-P chains is seen. From these observations, it can be said that in model III-AA and II-CA, one of the amorphous regions is contributing more in carrying load and hence providing more resistance to deformation than the other. Also it is seen in Figs. 7c and 11d that for model IV-AACA and III-AA, the stress corresponding to yield strain is same, this can be attributed to the covalent connectivity between the Amorphous-P chains and the non-bonded interaction between the Amorphous-L chains and the crystallites. In model III-AA also, the deformation in both the types of amorphous regions is almost equivalent. The Amorphous-P region is observed to deform slightly more than the Amorphous-L region. Here Amorphous-L region is deforming in the strain softening regime as the non-bonded interaction between the crystallite and the Amorphous-L region is strong enough to pull the chains apart resulting into stretching

elastic regime, the bond-stretching and bond angle energies are also contributing slightly in the increase of the potential energy. The increase in bond-stretching and bond angle energy is indicating that the bond length and bond angle are moving away from equilibrium resulting into stretching of polymer chain. Also the dihedral energy is steadily decreasing, which is indicative of energy dissipation due to rotation about the backbone chain towards the lower energy conformations (trans). The potential energy is increasing in all the three regimes (the elastic regime, strain softening regime and strain hardening regime) with major contribution from non-bonded energy which means, the non-bonded interactions are playing a dominating role during the tensile deformation. 3.3. Strain softening regime Upon visual analysis of the deformation trajectory, it is seen that, in all the models, the amorphous region is uncoiling and stretching significantly providing considerable permanent deformation in strain softening regime, whereas no visible deformation is observed in the crystalline region. Fig. 11 shows the snapshots of the models at the end of strain softening regime during the tensile deformation at the strain 12

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

and β-sheets being pulled-out of the crystallite resulting into a shear type deformation mechanism. From this, it can be said that, strain hardening regime represents that the load is now majorly carried by the crystalline region. The β-sheet crystallites, being inherently stiffer than the amorphous region, need larger load to deform resulting into the strain hardening. From these observations, the strain hardening regime can be correlated to the distortion of β-sheet crystallites. The distortion of the β-sheet crystallites in strain hardening regime can also be established from energetics point of view. Fig. 14(a) and (b) showing the variation in non-bonded energy of crystalline region, Amorphous-L chains and Amorphous-P chains. It is observed from Fig. 14(a) that in strain hardening regime, initially (0.5–0.75 strain) the non-bonded energy of the Amorphous-L chains is increases and stagnates thereafter, whereas, the non-bonded energy of the crystalline region keeps on increasing at a higher rate than the amorphous chains. Also in Fig. 14(b) it is seen that the interaction energy of crystalline region and the Amorphous-L chains keeps on increasing steeply until around 1.2 strain (until the start of fracture regime). From this steep increase in energy, it can be said that the interaction between crystalline and Amorphous-L regions plays a dominating role in load transmission in the system. The interaction energy of crystalline region and the Amorphous-P chains plays a little role in deformation stagnates after around 0.5 strain. In β-sheet crystallites, the β-strands are held together by H-bonds which are responsible for its high rupture strength. During distortion of the crystallites, the β-strands are seen to be continuously sliding out of the crystallite. A few of the β-strands (mostly in the inner domains of the crystallite) are also seen to be pulled out of the crystallite in stickslip type motion. Earlier computational studies investigating the failure mechanism and mechanical behavior of the crystalline region of B. mori SF or the spider dragline silk has shown that the crystallites exhibit stick-slip type mechanism [30,40,41] or second catch mechanism [80] during deformation, wherein the β-strands slide out of the crystallites as a result of applied force and reform the hydrogen bonds with new neighbors. This stick-slip motion or second catch, leads to energy dissipation and force peaks (corresponding to each step of stick-slip motion) in the force–displacement curve. Such force/stress peaks associated with the stick-slip of the β-strands are not observed in stressstrain plots (Fig. 7) in the current work, this can be attributed to two factors, (i) in the current work, the mechanism of force application is different from that of the above mentioned literatures. The tensile deformation here is carried-out by imposing a uniform strain on the simulation cell at any point of time (as mentioned in Section 2.3) whereas, in the above mentioned literatures, the crystallites are deformed by pulling a spring attached to the crystallite, with a constant velocity (also known as SMD method) [30,40,41] or electrostatically charged pulling force [80] and (ii) the stress values reported here are the overall stress values for the complete system, which is a combination of crystalline and amorphous region, whereas the force-displacement plot reported in the above mentioned literature are corresponding

and uncoiling. However, for II-CA configuration, it is seen that the one of the amorphous regions undergo more deformations as compared to the other, meaning that a major portion of load borne by only one of the amorphous regions. Here Amorphous-L region is observed to deform more than the Amorphous-P region unlike model IV-AACA and III-AA. This can be attributed to the lack of covalent connectivity between Amorphous-P chains in model II-CA unlike model IV-AACA and III-AA. For model I-NC, it is seen that the both the types of amorphous regions do not undergo any significant deformations as compared to the other models, meaning that in I-NC, very little portion of load is borne by both the amorphous regions, leading to a low strength. Also in model INC no strain hardening regime has been observed, after the yield point the stress keeps on decreasing with increase in strain. In model III-AA, it is also seen that, the strain softening regime ends at 0.75 strain whereas for model IV-AACA and II-CA, the strain softening regime ends around 0.5 strain. This implies that in model III-AA the load transfer from amorphous region to crystalline region is slower as compared to IV-AACA and II-CA. This can be attributed to the difference in connectivity of the Amorphous-L region and the β-strands in crystallites. In model III-AA, Amorphous-L chains are not covalently connected to the crystallites unlike IV-AACA and II-CA. From the above observations it can be said that the non-bonded interactions (Van der Waal's forces, H-bonds and electrostatic interactions) in amorphous region play an important role in providing flexibility to B. mori SF until yield point and in extension/stretching during the strain softening regime. It can be also said that, both Amorphous-L and Amorphous-P chains play a key role in energy absorption during deformation. Fig. 12 shows the crystalline region of model IV-AACA at 0.5 strain (50% strain) during the tensile deformation at the strain rate of 5 × 10−6 per fs, i.e., at the end of strain softening regime or start of strain hardening regime. It is seen that, until this point, the crystalline region still largely retains its structural integrity. From this point onwards the β-strands are seen to be beginning to slip out of the crystallites. 3.4. Strain hardening regime The amorphous chains separating the crystallites continue to uncoil and stretch in strain hardening regime also. Upon correlating the simulation trajectory with the stress-strain behavior of the models, it is observed that, in models IV-AACA, III-AA and II-CA, the crystallites are also distorting along with the uncoiling of the amorphous chains in the strain hardening regime. It is seen that, the segments of Amorphous-L chains in the vicinity of crystallites are nearly straightened (specifically, in model IV-AACA, the ends of the Amorphous-L chains which are covalently connected to the crystallite), as shown in Fig. 13. This localized straightening of the amorphous chains causes the backbone chains to take up the load which is further transferred to the crystallites. The distortion of crystallites is seen to be occurring by various β-strands

β –strand not connected with the amorphous chain (blue)

β –strand connected to amorphous chains (yellow) starts to slip out of the crystallite

β –strand covalently connected with the amorphous chain (yellow) Fig. 12. Snapshot of the crystalline region of model IV-AACA at 0.5 strain (50%) at 5 × 10−6 per fs strain rate, at this point β-strands start to slip out of the crystallites. Until this point the crystalline region still retains its structural integrity but the slip-out of β-strands mark the initiation of distortion in crystallites. 13

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Locally uncoiled Amorphous-L chains near the covalent connection to βstrands in crystallites Upon further deformation the load is transferred to the β-sheet crystallites

Fig. 13. Snapshots of deformation trajectory of model IV-AACA at 0.5 strain (50%).

Non-bonded interaction energy: IV-AACA

Non-bonded energy: IV-AACA

(b)

(a)

Fig. 14. (a) Variation in in non-bonded energy with strain for crystalline region, Amorphous-L chains and Amorphous-P chains for model IV-AACA during the tensile deformation at strain rate 5 × 10−6 per fs. (b) Variation in non-bonded interaction energy between crystalline and Amorphous-L chains and crystalline and Amorphous-P chains for model IV-AACA during the tensile deformation at strain rate 5 × 10−6 per fs.

Fig. 15. Snapshots of β-sheet crystallite of model IV-AACA at 1.3 strain. The crystallites are observed to be significantly distorted.

14

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

(a) IV-AACA: Start of fracture regime (at 1.3 strain)

(b) III-AA: Start of fracture regime (at 1 strain)

(c) II-CA: Start of fracture regime (at 0.8 strain)

Start of fracture regime (d) I-NC: Start of fracture regime (at 0.4 strain)

(e) Stress Vs Strain: Strain rate 5X10-6 per fs Fig. 16. Snapshots of the models at start of fracture regime (a) IV-AACA, (b) III-AA, (c) II-CA, (d) I-NC at strain rate of 5 × 10−6 per fs and (e) stress-strain plot of the four models at 5 × 10−6 per fs demarcating the start of fracture regime.

pronounced and distinct, where as in III-AA and II-CA, the strain softening regime is followed by a smaller strain hardening regime with a plateau and hump like regimes, respectively. Also, in model IV-AACA, the fracture regime is starting at a higher strain and stress value whereas in model III-AA and II-CA the fracture regime is starting at a lower stress values. This can be attributed to the collective covalent connectivity between the crystallite and Amorphous-L chains and between Amorphous-P chains in IV-AACA producing simultaneous distortion of crystalline region and the stretching of both the types of

to the crystalline region only. Fig. 15 shows the deformed crystalline region of model IV-AACA at the end of strain hardening regime where the β-strands are shown to be pulled out of the crystallites. Further analysis shows that, in model IV-AACA, III-AA and II-CA, the β-sheet crystallites have reached significant distortion, by the end of strain hardening regime. Fig. 16(a)–(d) shows the snapshots of the models under the tensile deformation at the strain rate of 5 × 10−6 per fs at the start of fracture regime (end of the strain hardening regime). It is observed that the strain hardening regime in IV-AACA is more 15

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Fig. 17. (a) Elastic modulus, (b) tensile strength and (c) yield strain of the four models at different strain rates.

along with the distortion of crystallites only one of the amorphous regions (i.e., Amorphous-L chains) is deforming resulting into even lower tensile strength than model IV-AACA and III-AA. The distortion of the crystallites and uncoiling of the amorphous chains both, as mentioned above, take part in dissipation of the strain energy whereas, the pull-out of β-strands and β-sheets causing distortion of crystallites plays a major role in achieving a high strength. The failure mechanism for the crystallites observed here has also been reported in earlier studies investigating the crystalline region of silkworm silk or spider dragline silk [6,25,30]. Some of the computational studies investigating the intramolecular β-sheet crystallites have also reported the unfolding of the intramolecular bonded β-sheets into a single strand as source of high strength and strain-hardening observed in silks [21,81]. Such an unfolding of intramolecular β-sheets is not observed in intermolecular β-sheet crystallites (used in current study) due to lack of any bonded connection (β-turns) between the adjacent strands. In intermolecular β-sheet crystallites, pull-out of β-strands (or β-sheets) is observed as the deformation mechanism. In case of 5 × 10−4 per fs strain rate, models IV-AACA and II-CA show initial linear elastic regime followed immediately strain hardening regime without any intermediate strain softening regime, and, model III-AA shows a flat plateau like regime between the yield point and the strain hardening regime. The reason for this to happen could be the participation of both bonded and non-bonded interaction in load transfer from Amorphous-L chains to crystalline region in model IVAACA and II-CA, whereas in III-AA only non-bonded interaction is responsible for load transfer which at higher strain rates is causing a drawing like occurrence in III-AA. In model I-NC, no deformation is observed in amorphous regions during the strain hardening regime as compared to IV-AACA, III-AA and II-CA and the β-sheet crystallites also distort at lower stress and strain

amorphous regions during strain hardening regime, whereas in III-AA, the absence of covalent connectivity between the crystallite and Amorphous-L chains, leaves only the non-bonded interaction between crystallite and Amorphous-L chains responsible for the load transfer between amorphous and crystalline region, resulting into a gradual and weaker dissemination of load. In model II-CA, along with the crystallites only one of the amorphous regions (Amorphous-L chains) is contributing in carrying the load. The Amorphous-P chains here are seen to be undergoing almost no deformation/uncoiling as the inter-chain nonbonded interaction between the adjacent chains is not strong enough to overcome the intra-chain non-bonded interaction and hence virtually providing no resistance to deformation. From these observations, it can be said that model IV-AACA has more resistance to deformation than III-AA and II-CA. In model IV-AACA, the collective covalent connectivity, as mentioned above, is bringing into larger distortion of crystalline region. Various β-strands are pulled out from the crystallites as a result of combined bonded and non-bonded interaction between Amorphous-L chains and crystallites and the non-bonded interaction at the interface of Amorphous-P chains and the crystallites is contributing in shearing of β-sheets during deformation. This distortion of crystallites along with the stretching of both the types of amorphous regions during strain hardening regime is resulting into a higher tensile strength for model IV-AACA. In model III-AA also both the amorphous regions are seen to be uncoiling, but the lack of covalent connection between Amorphous-L chains and crystallites is causing lesser distortion in crystallites, as shown in Fig. 16(a) and (b). That is, the non-bonded interaction between the amorphous chains and the crystallites is not strong enough for extensive load transfer from amorphous chains to crystalline region hence causing lesser distortion of crystallites and resulting into lower tensile strength in III-AA than model IV-AACA. While in model II-CA, 16

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Table 1 shows the mechanical properties of the four B. mori SF nanostructure models obtained in present work from MD simulations of tensile deformation at three different strain rates, and its comparison with the values for silk fiber reported in earlier studies. For model IIIAA, it is seen that the elastic modulus values are comparable with the values obtained for model IV-AACA. Also, the elastic modulus values show an increasing trend starting from a non-bonded configuration (model I-NC) to a completely bonded one (IV-AACA) indicating that the highest elastic modulus is achieved when the Amorphous-P chains are bonded, whereas, the bond between Amorphous-L chains and crystalline region is more effective in imparting a high tensile strength and yield strain. Also the difference between the values of elastic modulus, tensile strength and yield strain are more prominent at higher strain rate. The elastic moduli obtained for model IV-AACA from simulations are 15.5 ± 0.5, 10.4 ± 1.3 and 7.4 ± 0.3 GPa corresponding to strain rates of 5 × 10−4, 5 × 10−5 and 5 × 10−6 per fs, respectively, which are found to be well within the range of 7 to 18 GPa elastic moduli values of B. mori SF reported in earlier experimental studies [7,39,59,67,84]. This range is reported for different type of B. mori. SF, i.e. either native silk or artificially spun silk fiber, forced reeling or regenerated silk, different spinning conditions, at both macro and micro lengthscales. The tensile properties for B. mori SF model IV-AACA at strain rate of 5 × 10−6 per fs, i.e., elastic modulus value of 7.4 ± 0.3 GPa, tensile strength value of ~340 MPa and yield strain value of ~13%, are in close association with the tensile properties of B. mori SF reported elsewhere [4,7,39,59,67,84]. Furthermore, the comparison of tensile properties for B. mori SF model IV-AACA at strain rate of 5 × 10−6/fs with the tensile properties of spider dragline silk and flagelliform silk (N. clavipes, L. Hesperus, A. diadematus and N. edulis), reported in earlier studies, shows that the elastic modulus value of B. mori SF model IVAACA (7.4 ± 0.3 GPa) is close to the range of elastic modulus values reported for the spider dragline silk (8–17 GPa) [82,83] and much higher than the values reported for flagelliform silk (~0.003 GPa) [10,12,85] whereas the tensile strength values for B. mori SF model IVAACA (~340 MPa) is much lower than the reported tensile strength values for the spider dragline silk (800–1700 MPa) [82,83] and comparable to the tensile strength of the flagelliform silk (500–700 MPa) [10,12,85]. Fig. 19 shows the comparison of the stress-strain plot for B. mori SF single filament of ~10 μm in diameter [8], B. mori silk fiber (degummed cocoon silk) of ~12 μm in diameter [86], and stress-strain plot obtained for the B. mori SF model IV-AACA at 5 × 10−6 per fs strain rate obtained from MD simulation in current study. It is observed that the behavior of stress-strain plot obtained from MD simulation is in close correlation with stress-strain behavior of B. mori SF reported elsewhere [8]. The stress-strain plot obtained from the MD simulation is strictly following the plot reported in literature upto around 0.03 strain value (3% strain), also in this range, the stress-strain plots are typically linear. After 0.03 strain value, there is visible deviation in the stress-strain plot from MD simulation from that of the one reported in literature. This difference in the behavior of the stress-strain plots can be attributed to the difference between stress calculations in macro-scale tensile test and MD Simulations. In the experimental result, the stresses obtained are average values involving combined contributions from inter-fibril sliding, inter-molecular sliding and molecular deformations. However, in MD simulations, only a nano-scale sample is considered and the stresses obtained include only the inter-molecular sliding and molecular deformations. Though the reported tensile properties and behavior for B. mori SF are for various lengthscales, the close correlation of the tensile properties and behavior obtained from MD simulations for a nanoscale size SF material model suggests that this nanoscale model is capturing the building block B. mori SF material properties.

values as major load is borne by crystalline region only and the load is enough to overcome the H-bonds in β-sheet crystallites resulting into a much lower tensile strength. In the configuration I-NC, little deformation is observed in amorphous region during the strain hardening regime as compared to IVAACA and III-AA and β-sheet crystallites distort at lower stress and strain values since major load is borne by crystalline region and the load is enough to overcome the H-bonds in β-sheet crystallites. From all these observations it can be said that, the H-bonds in βsheet crystallites and the covalent connectivity between the Amorphous-L the crystalline regions and the covalent connectivity between the Amorphous-P chains plays a vital role in deformation mechanism of B. mori SF in strain hardening regime. Fig. 17(a–c) shows the elastic modulus, ultimate tensile strength and yield strain of the four models and at different strain rates. From Fig. 17(a), it is observed that model IV-AACA shows the highest elastic modulus in all the cases, whereas, model I-NC has the lowest elastic modulus. It can also be seen that model IV-AACA shows highest tensile strength and yield strain. The elastic modulus, strength and yield strain values are seen to be increasing with increase in strain rate. At lower strain rates the deformation is majorly dominated by the non-bonded interaction whereas at higher strain rates, along with the uncoiling of turns (overcoming the non-bonded interaction), the backbone chains are also seen to be stretching significantly all through the deformation and hence substantially participating in load carrying. The same can be substantiated from the variation of bonded and non-bonded energy of the model IV-AACA during deformation at two strain rates (5 × 10−5/ fs and 5 × 10−6/fs), as shown in Fig. 18. It is observed that the at higher strain rate both bonded and non-bonded interaction are increasing substantially with deformation whereas at lower strain rates, only the non-bonded energy is seen to be increasing considerably with deformation. This augmented stretching of the backbone chains during deformation at higher strain rates is resulting into higher elastic modulus and tensile strength of the models. Also the stretching of backbone of the molecular chains require larger loads due to the presence of stronger covalent bond (compared to electrostatic or van der Waal's interaction) resulting into higher tensile strength at higher strain rates. The increase in yield strain observed with increase in strain rates can be attributed to the relaxation time at different strain rates. During deformation at lower strain rates the system has more relaxation time (time to uncoil or rearrange to new stable configuration) than at higher strain rate. That is, at lower strain rate, the viscous behavior is dominating resulting into permanent deformation at lower strain values than at higher strain rates.

Fig. 18. Variation in bonded and non-bonded energy with strain for model IVAACA during the tensile deformation at two strain rates 5 × 10−5 per fs and 5 × 10−6 per fs. E_b and E_nb represent bonded energy and non-bonded energy, respectively. 17

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

Table 1 Comparison of tensile properties of the four B. mori SF nanostructure models at different strain rates obtained in present work, and the values for silk fiber reported elsewhere.

Present work

Strain rates

Parameter/model

Modulus (GPa)

Tensile strength (MPa)

Yield strain (%)

5 × 10−4 per fs

I-NC II-CA III-AA IV-AACA I-NC II-CA III-AA IV-AACA I-NC II-CA III-AA IV-AACA Water extracted Native degummed – – – – – –

12.2 ± 0.4 13.4 ± 0.4 14.6 ± 0.5 15.5 ± 0.5 7.8 ± 0.9 8 ± 0.9 9.4 ± 1.3 10.4 ± 1.3 6.4 ± 1.1 6.8 ± 0.9 7.1 ± 1.3 7.4 ± 0.3 8a 13.6 7b – 7 10–17 8–17 0.003

3300 4700 5300 6500 810 820 850 1100 270 270 330 340 – – 300 500 600 300–740 800–1700 500–700

70 60 70 60 25 25 21 21 25 20 18 13 – – 15 15 18 – – –

5 × 10−5 per fs

5 × 10−6per fs

Reported in literature

Kaplan et al. [39] Fritz et al. [67] Shao et al. [8] X. Liu et al. [12] L.D. Koh et al. [3] Dragline silk [83,84] Flagelliform silk [10]

a b

Electrospun silk fiber. Non-transgenic silk fiber.

properties by modulating the extent of crystallinity, spinning process and the morphology (for regenerated silk) of B. mori SF have also been reported in previous studies, however, the values reported in present work correspond to naturally spun B. mori SF [8,26,33,37,67,87]. It is also observed that the strain rate has a significant effect on the mechanical behavior of SF at nanoscale. For different strain rates, the stress-strain plots obtained from MD simulations are also different. The differences in the plot beyond the initial linear regime can be understood by analyzing the deformation from a polymer standpoint. In a polymer system, during tensile deformation lower strain rates favours the viscous flow behavior of the polymer while the higher strain rates favours the elastic load bearing behavior of the polymer. Beyond the yielding of material, in case of deformation at lower strain rates, the uncoiling of the amorphous chains play a dominating role as compared to the distortion of the β-sheet crystallites while in case of higher strain rates, the shear distortion of the β-sheet crystallites play a governing role. The study also reveals that the model IV-AACA has higher elastic modulus and tensile strength values than other three models at any particular strain rate. This implies that the connectivity, that is, covalent or non-bonded connectivity, between the crystalline and amorphous region plays an important role in governing the mechanical behavior of B. mori SF. The results obtained in this study should be considered along with the computational assumptions made in MD simulations. First, the nanoscale model developed in this study consists of 8 β-sheet crystallites to model the crystalline region and linker 6 from the full length residue sequence to model the amorphous regions and not the whole 5263 residues of the B. mori SF heavy chain as it is computationally very expensive to model and simulate its large size. Second, the MD simulations were performed with the models in completely anhydrous conditions. Also the simulations performed on the developed models are only tensile in nature. In order to obtain better insights of nanoscale mechanical properties and behavior of SF and mechanisms involved in deformation further qualitative investigations are needed. Primarily, the deformation for all the four models is observed to be a combination of tensile and shear deformation in the simulations. Initially the amorphous region undergoes tensile deformation due to uncoiling, sliding and stretching of the chains resulting into the shear deformation of the crystalline region. This behavior is in accordance with the mechanisms proposed elsewhere [3,30,37]. Also despite the differences in primary structure (amino acid residue sequence) of spider

Fig. 19. Comparison of the stress-strain plot for B. mori silk fibroin single filament of ~10 μm in diameter [8], B. mori silk fiber (degummed cocoon silk) of ~12 μm in diameter [86] and stress-strain plot obtained for the model IV-AACA at 5 × 10−6 per fs strain rate from MD simulation in present study.

4. Discussion and conclusions In this study, the deformation mechanisms of B. mori SF have been investigated using a set of MD simulations. For this purpose, four different computational atomistic models composed of both crystalline and amorphous regions are developed and tested for mechanical properties and deformation under tensile loading at three different strain rates of 5 × 10−6, 5 × 10−5 and 5 × 10−4 per fs. The elastic modulus values obtained from MD simulations are in the range of 6 to 15 GPa for the models at different strain rates. The elastic modulus values for the four models are in the following order: I-NC < IICA < III-AA < IV-AACA at all the three strain rates. The obtained tensile properties (elastic modulus of ~7.4 GPa and tensile strength of ~340 MPa) and the stress-strain behavior for B. mori SF model IV-AACA (at 5 × 10−6 per fs strain rate) has been found to be in close correlation with the range of tensile properties and behavior reported elsewhere [4,12,36,39,59,67]. Enhanced mechanical 18

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

silk from B. mori SF, the computational studies for spider silk also show similar deformation mechanism at nanoscale. During tensile deformation, both crystalline and amorphous regions in B. mori SF nanostructure play a defining role at different strain levels. The elastic modulus of B. mori SF is governed by the non-bonded interactions between the chains in amorphous region and with chains in crystalline region, whereas, the tensile strength of B. mori SF is governed by the crystalline region. The strain softening in B. mori SF is a result of uncoiling and stretching of the amorphous chains. At higher levels of strain, β-sheet crystallites start deforming in shear distortion modes resulting into strain hardening. Insights gained here are important in advancing the understanding about multi-scale mechanics of B. mori SF and for studies focused on developing silk based biomaterials with tailored properties. Conclusions drawn in this work complement the findings from earlier experimental studies conforming to the reported mechanical properties and predicted mechanics of deformation. The computational framework developed here can be further used for carrying out investigations on silk or silk like materials to develop a broader understanding of structure-property relationships, and for studies on silk based nanocomposites. Overall, this work further contributes towards advancing the understanding of B. mori SF mechanics at nanoscale while modelling both amorphous and crystalline phases together under tensile deformations by developing atomistic models and a MD computational framework.

[7] J. Pérez-Rigueiro, C. Viney, J. Llorca, M. Elices, Mechanical properties of single-brin silkworm silk, J. Appl. Polym. Sci. 75 (2000) 1270–1277, https://doi.org/10.1002/ (SICI)1097-4628(20000307)75:10<1270::AID-APP8>3.0.CO;2-C. [8] Z. Shao, F. Vollrath, Surprising strength of silkworm silk, Nature 418 (2002) 741, https://doi.org/10.1038/418741a. [9] S. Keten, M.J. Buehler, Atomistic model of the spider silk nanostructure, Appl. Phys. Lett. 96 (2010) 25–28, https://doi.org/10.1063/1.3385388. [10] C.Y. Hayashi, N.H. Shipley, R.V. Lewis, Hypotheses that correlate the sequence, structure, and mechanical properties of spider silk proteins, Int. J. Biol. Macromol. 24 (1999) 271–275, https://doi.org/10.1016/S0141-8130(98)00089-0. [11] C.-Z. Zhou, F. Confalonieri, M. Jacquet, R. Perasso, Z.-G. Li, J. Janin, Silk fibroin: structural implications of a remarkable amino acid sequence, Proteins Struct. Funct. Genet. 44 (2001) 119–122, https://doi.org/10.1002/prot.1078. [12] X. Liu, K. Zhang, Silk fiber — molecular formation mechanism , structure- property relationship and advanced applications, Oligomerization Chem. Biol. Compd. (2014) 69–102, https://doi.org/10.5772/57611. [13] T.R. Scheibel, J.G. Hardy, L.M. Ro, Polymeric materials based on silk proteins, Poly 49 (2008) 4309–4327, https://doi.org/10.1016/j.polymer.2008.08.006. [14] O. Tokareva, M. Jacobsen, M. Buehler, J. Wong, D.L. Kaplan, Structure-functionproperty-design interplay in biopolymers: spider silk, Acta Biomater. 10 (2014) 1612–1626, https://doi.org/10.1016/j.actbio.2013.08.020. [15] T. Asakura, K. Okushita, M.P. Williamson, Analysis of the structure of Bombyx mori silk fibroin by NMR, Macromolecules 48 (2015) 2345–2357, https://doi.org/10. 1021/acs.macromol.5b00160. [16] Y. Takahashi, M. Gehoh, K. Yuzuriha, Structure refinement and diffuse streak scattering of silk (Bombyx mori), Int. J. Biol. Macromol. 24 (1999) 127–138, https://doi.org/10.1016/S0141-8130(98)00080-4. [17] E. Sintya, P. Alam, Self-assembled semi-crystallinity at parallel β-sheet nanocrystal interfaces in clustered MaSp1 (spider silk) proteins, Mater. Sci. Eng. C 58 (2016) 366–371, https://doi.org/10.1016/j.msec.2015.08.038. [18] E. Sintya, P. Alam, Localised semicrystalline phases of MaSp1 proteins show high sensitivity to overshearing in beta-sheet nanocrystals, Int. J. Biol. Macromol. 92 (2016) 1006–1011, https://doi.org/10.1016/j.ijbiomac.2016.07.081. [19] T. Asakura, M. Okonogi, Y. Nakazawa, K. Yamauchi, Structural Analysis of Alanine Tripeptide With Antiparallel and Parallel-sheet Structures in Relation to the Analysis of Mixed-sheet Structures in Samia cynthia ricini Silk Protein Fiber Using Solid-state NMR Spectroscopy, (2006), pp. 6231–6238, https://doi.org/10.1021/ ja060251t. [20] G. Bratzel, M.J. Buehler, Sequence-structure correlations in silk: poly-Ala repeat of N. clavipes MaSp1 is naturally optimized at a critical length scale, J. Mech. Behav. Biomed. Mater. 7 (2012) 30–40, https://doi.org/10.1016/j.jmbbm.2011.07.012. [21] N. Du, Z. Yang, X.Y. Liu, Y. Li, H.Y. Xu, Structural origin of the strain-hardening of spider silk, Adv. Funct. Mater. 21 (2011) 772–778, https://doi.org/10.1002/adfm. 201001397. [22] S. Keten, M.J. Buehler, Nanostructure and molecular mechanics of spider dragline silk protein assemblies, J. R. Soc. Interface 7 (2010) 1709–1721, https://doi.org/ 10.1098/rsif.2010.0149. [23] S.A. Fossey, G. Nfmethy, K.D. Gibson, H.A. Scheraga, Conformational energy studies of beta-sheets of model silk fibroin peptides . 1 . Sheets of poly ( Ala-Cly ) chains, Bioplymers 31 (1991) 1529–1541, https://doi.org/10.1002/bip. 360311309. [24] X. Chen, Z. Shao, D.P. Knight, F. Vollrath, Conformation transition kinetics of Bombyx mori silk protein, PROTEINS Struct. Funct. Bioinforma. 68 (2007) 223–231, https://doi.org/10.1002/prot.21414. [25] Y. Cheng, L.D. Koh, D. Li, B. Ji, M.Y. Han, Y.W. Zhang, On the strength of beta-sheet crystallites of Bombyx mori silk fibroin, J. R. Soc. Interface 11 (2014) 20140305, , https://doi.org/10.1098/rsif.2014.0305. [26] P. Dubey, S. Murab, S. Karmakar, P.K. Chowdhury, S. Ghosh, Modulation of selfassembly process of fibroin: an insight for regulating the conformation of silk biomaterials, Biomacromolecules 16 (2015) 3936–3944, https://doi.org/10.1021/acs. biomac.5b01258. [27] G. Fang, Y. Huang, Y. Tang, Z. Qi, J. Yao, Z. Shao, Insights into silk formation process: correlation of mechanical properties and structural evolution during artificial spinning of silk fibers, ACS Biomater. Sci. Eng. 2 (2016) 1992–2000, https:// doi.org/10.1021/acsbiomaterials.6b00392. [28] J. Guan, D. Porter, F. Vollrath, Silks cope with stress by tuning their mechanical properties under load, Polymer (Guildf) 53 (2012) 2717–2726, https://doi.org/10. 1016/j.polymer.2012.04.017. [29] J. Guan, Y. Wang, B. Mortimer, C. Holland, Z. Shao, D. Porter, F. Vollrath, Glass transitions in native silk fibres studied by dynamic mechanical thermal analysis, Soft Matter 12 (2016) 5926–5936, https://doi.org/10.1039/c6sm00019c. [30] S. Keten, Z. Xu, B. Ihle, M.J. Buehler, Nanoconfinement controls stiffness, strength and mechanical toughness of Β-sheet crystals in silk, Nat. Mater. 9 (2010) 359–367, https://doi.org/10.1038/nmat2704. [31] I. Krasnov, I. Diddens, N. Hauptmann, G. Helms, M. Ogurreck, T. Seydel, S.S. Funari, M. Müller, Mechanical properties of silk: interplay of deformation on macroscopic and molecular length scales, Phys. Rev. Lett. 100 (2008) 2–5, https:// doi.org/10.1103/PhysRevLett.100.048104. [32] B.D. Lawrence, S. Wharram, J.A. Kluge, G.G. Leisk, F.G. Omenetto, M.I. Rosenblatt, D.L. Kaplan, Effect of hydration on silk film material properties, Macromol. Biosci. 10 (2010) 393–403, https://doi.org/10.1002/mabi.200900294. [33] R. Nazarov, H.J. Jin, D.L. Kaplan, Porous 3-D scaffolds from regenerated silk fibroin, Biomacromolecules 5 (2004) 718–726, https://doi.org/10.1021/ bm034327e. [34] Z. Qin, M.J. Buehler, Cooperative deformation of hydrogen bonds in beta-strands and beta-sheet nanocrystals, Phys. Rev. E 82 (2010) 1–9, https://doi.org/10.1103/

Data availability All information/data required to reproduce the results reported in this article are mentioned/cited in the manuscript. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors gratefully acknowledge High Performance Computing (HPC) facility of IIT Delhi, using which part of this work was performed. Authors also thank Prof. Sourabh Ghosh, Textile Engineering Department, IIT Delhi for fruitful discussions, which helped improve the quality of work. Appendix A. Supplementary data Video of tensile deformation of model IV-AACA at 5 × 10−6 per fs strain rate. Supplementary data to this article can be found online at https://doi.org/10.1016/j.msec.2019.110414. References [1] C. Vepari, D.L. Kaplan, Silk as a biomaterial, Prog. Polym. Sci. 32 (2007) 991–1007, https://doi.org/10.1016/j.progpolymsci.2007.05.013. [2] J. Melke, S. Midha, S. Ghosh, K. Ito, S. Hofmann, Silk fibroin as biomaterial for bone tissue engineering, Acta Biomater. 31 (2015) 1–16, https://doi.org/10.1016/j. actbio.2015.09.005. [3] L.D. Koh, Y. Cheng, C.P. Teng, Y.W. Khin, X.J. Loh, S.Y. Tee, M. Low, E. Ye, H.D. Yu, Y.W. Zhang, M.Y. Han, Structures, mechanical properties and applications of silk fibroin materials, Prog. Polym. Sci. 46 (2015) 86–110, https://doi.org/10.1016/j. progpolymsci.2015.02.001. [4] G.H. Altman, F. Diaz, C. Jakuba, T. Calabro, R.L. Horan, J. Chen, H. Lu, J. Richmond, D.L. Kaplan, Silk-based biomaterials, Biomaterials 24 (2003) 401–416, https://doi.org/10.1016/S0142-9612(02)00353-8. [5] O. Hakimi, D.P. Knight, F. Vollrath, P. Vadgama, Spider and mulberry silkworm silks as compatible biomaterials, Compos. Part B Eng. 38 (2007) 324–337, https:// doi.org/10.1016/j.compositesb.2006.06.012. [6] A. Nova, S. Keten, N.M. Pugno, A. Redaelli, M.J. Buehler, Molecular and nanostructural mechanisms of deformation, strength and toughness of spider silk fibrils, Nano Lett. 10 (2010) 2626–2634, https://doi.org/10.1021/nl101341w.

19

Materials Science & Engineering C 108 (2020) 110414

M. Patel, et al.

[60] D.Z. Bermudez, R.F.P. Pereira, M.M. Silva, Bombyx mori silk fibers: an outstanding family of materials, Macromol. Mater. Eng. 300 (2015) 1171–1198, https://doi. org/10.1002/mame.201400276. [61] Y. Cheng, L.D. Koh, D. Li, B. Ji, Y. Zhang, J. Yeo, G. Guan, M.Y. Han, Y.W. Zhang, Peptide-graphene interactions enhance the mechanical properties of silk fibroin, ACS Appl. Mater. Interfaces 7 (2015) 21787–21796, https://doi.org/10.1021/ acsami.5b05615. [62] RCSB Protein Data Bank, https://www.rcsb.org/. [63] A.D. William Humphrey, K.S., VMD:visual molecular dynamics, J. Mol. Graph. 14 (1996) 33–38, https://doi.org/10.1016/0263-7855(96)00018-5. [64] A. Lamiable, P. Thévenet, J. Rey, M. Vavrusa, P. Derreumaux, P. Tufféry, PEPFOLD3: faster de novo structure prediction for linear peptides in solution and in complex, Nucleic Acids Res. 44 (2016) W449–W454, https://doi.org/10.1093/nar/ gkw329. [65] Y. Shen, J. Maupetit, P. Derreumaux, P. Tufféry, Improved PEP-FOLD approach for peptide and miniprotein structure prediction, J. Chem. Theory Comput. 10 (2014) 4745–4758, https://doi.org/10.1021/ct500592m. [66] P. Thévenet, Y. Shen, J. Maupetit, F. Guyon, P. Derreumaux, P. Tufféry, PEP-FOLD: an updated de novo structure prediction server for both linear and disulfide bonded cyclic peptides, Nucleic Acids Res. 40 (2012) 288–293, https://doi.org/10.1093/ nar/gks419. [67] C. Fu, Z. Shao, V. Fritz, Animal silks: their structures, properties and artificial production, Chem. Commun. 43 (2009) 6515–6529, https://doi.org/10.1039/ b911049f. [68] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. (1995), https://doi.org/10.1006/jcph.1995.1039. [69] CHARMM Force Field Files, http://mackerell.umaryland.edu/charmm_ff.shtml. [70] N. Foloppe, A.D. MacKerell, All-atom empirical force field for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data, J. Comput. Chem. 21 (2000) 86–104, https://doi.org/10. 1002/(SICI)1096-987X(20000130)21:2<86::AID-JCC2>3.0.CO;2-G. [71] A.D. MacKerell, N. Banavali, N. Foloppe, Development and current status of the CHARMM force field for nucleic acids, Biopolymers 56 (2000) 257–265, https:// doi.org/10.1002/1097-0282(2000)56:4<257::AID-BIP10029>3.0.CO;2-W. [72] D.K. Dubey, V. Tomar, Role of the nanoscale interfacial arrangement in mechanical strength of tropocollagen-hydroxyapatite-based hard biomaterials, Acta Biomater. 5 (2009) 2704–2716, https://doi.org/10.1016/j.actbio.2009.02.035. [73] A.K. Nair, A. Gautieri, S.-W. Chang, M.J. Buehler, Molecular mechanics of mineralized collagen fibrils in bone, Nat. Commun. 4 (2013) 1724–1729, https://doi.org/ 10.1038/ncomms2720. [74] S. Deng, Tensile deformation of semi - crystalline polymers by molecular dynamics simulation, Iran. Polym. J. 26 (2017) 903–911, https://doi.org/10.1007/s13726017-0577-2. [75] D. Hossain, M.A. Tschopp, D.K. Ward, J.L. Bouvard, P. Wang, M.F. Horstemeyer, Molecular dynamics simulations of deformation mechanisms of amorphous polyethylene, Polymer (Guildf) 51 (2010) 6071–6083, https://doi.org/10.1016/j. polymer.2010.10.009. [76] A.V. Lyulin, N.K. Balabaev, M.A. Mazo, M.A.J. Michels, Molecular dynamics simulation of uniaxial deformation of glassy amorphous atactic polystyrene, Macromolecules 37 (2004) 8785–8793, https://doi.org/10.1021/ma049737p. [77] K. Yashiro, T. Ito, Y. Tomita, Molecular dynamics simulation of deformation behavior in amorphous polymer : nucleation of chain entanglements and network structure under uniaxial tension, Int. J. Mech. Sci. 45 (2003) 1863–1876, https:// doi.org/10.1016/j.ijmecsci.2003.11.001. [78] D.H. Tsai, The virial theorem and stress calculation in molecular dynamics, J. Chem. Phys. 70 (1979) 1375–1382, https://doi.org/10.1063/1.437577. [79] A.K. Subramaniyan, C.T. Sun, Continuum interpretation of virial stress in molecular simulations, Int. J. Solids Struct. 45 (2008) 4340–4346, https://doi.org/10.1016/j. ijsolstr.2008.03.016. [80] M. Pahlevan, M. Toivakka, P. Alam, Electrostatic charges instigate ‘concertina-like’ mechanisms of molecular toughening in MaSp1 (spider silk) proteins, Mater. Sci. Eng. C 41 (2014) 329–334. [81] P. Alam, Protein unfolding versus β-sheet separation in spider silk nanocrystals, Adv. Nat. Sci. Nanosci. Nanotechnol. 5 (2014), https://doi.org/10.1088/20436262/5/1/015015. [82] P.M. Cunniff, S.A. Fossey, M.A. Auerbach, J.W. Song, D.L. Kaplan, W.W. Adams, R.K. Eby, D. Mahoney, D.L. Vezie, Mechanical and thermal properties of dragline silk from the spider Nephila clavipes, Polym. Adv. Technol. 5 (1994) 401–410. [83] Y. Liu, A. Sponner, D. Porter, F. Vollrath, Proline and Processing of Spider Silks, (2008), pp. 116–121. [84] D. Ebrahimi, O. Tokareva, N.G. Rim, J.Y. Wong, D.L. Kaplan, M.J. Buehler, Silk-its mysteries, how it is made, and how it is used, ACS Biomater. Sci. Eng. 1 (2015) 864–876, https://doi.org/10.1021/acsbiomaterials.5b00152. [85] X. Hu, K. Vasanthavada, K. Kohler, S. Mcnary, A.M.F. Moore, C.A. Vierra, Molecular mechanisms of spider silk, 63 (2006) 1986–1999, https://doi.org/10.1007/s00018006-6090-y. [86] C. Guo, J. Zhang, X. Wang, A.T. Nguyen, X.Y. Liu, D.L. Kaplan, Comparative study of strain-dependent structural changes of silkworm silks : insight into the structural origin of, Small 1702266 (2017) 1–8, https://doi.org/10.1002/smll.201702266. [87] L. Wang, R. Nemoto, M. Senna, Changes in microstructure and physico-chemical properties of hydroxyapatite-silk fibroin nanocomposite with varying silk fibroin content, J. Eur. Ceram. Soc. 24 (2004) 2707–2715, https://doi.org/10.1016/j. jeurceramsoc.2003.09.006.

PhysRevE.82.061906. [35] N. Rim, E.G. Roberts, D. Ebrahimi, N. Dinjaski, M.M. Jacobsen, Z. Mart, M.J. Buehler, D.L. Kaplan, J.Y. Wong, Predicting silk fiber mechanical properties through multiscale simulation and protein design, ACS Biomater. Sci. Eng. 3 (2017) 1542–1556, https://doi.org/10.1021/acsbiomaterials.7b00292. [36] T. Seydel, K. Kölln, I. Krasnov, I. Diddens, N. Hauptmann, G. Helms, M. Ogurreck, S.G. Kang, M.M. Koza, M. Müller, Silkworm silk under tensile strain investigated by synchrotron X-ray diffraction and neutron spectroscopy, Macromolecules 40 (2007) 1035–1042, https://doi.org/10.1021/ma0624189. [37] J. Sirichaisit, V.L. Brookes, R.J. Young, F. Vollrath, Analysis of structure/property relationships in silkworm (Bombyx mori) and spider dragline (Nephila edulis) silks using raman spectroscopy, Biomacromolecules 4 (2003) 387–394, https://doi.org/ 10.1021/bm0256956. [38] D. Su, M. Yao, J. Liu, Y. Zhong, X. Chen, Z. Shao, Enhancing mechanical properties of silk fibroin hydrogel through restricting the growth of β - sheet domains, ACS Appl. Mater. Interfaces 9 (2017) 17489–17498, https://doi.org/10.1021/acsami. 7b04623. [39] M. Wang, H.-J. Jin, D.L. Kaplan, G.C. Rutledge, Mechanical properties of electrospun silk fibers, Macromolecules 37 (2004) 6856–6864, https://doi.org/10.1021/ MA048988V. [40] S. Xiao, W. Stacklies, M. Cetinkaya, B. Markert, F. Gräter, Mechanical response of silk crystalline units from force-distribution analysis, Biophys. J. 96 (2009) 3997–4005, https://doi.org/10.1016/j.bpj.2009.02.052. [41] C. Xu, D. Li, Y. Cheng, M. Liu, Y. Zhang, B. Ji, Pulling out a peptide chain from βsheet crystallite: propagation of instability of H-bonds under shear force, Acta Mech. Sinica 31 (2015) 416–424, https://doi.org/10.1007/s10409-015-0404-y. [42] Y. Yang, C. Dicko, C.D. Bain, Z. Gong, R.M.J. Jacobs, Z. Shao, A.E. Terry, F. Vollrath, Behavior of silk protein at the air – water interface, Soft Matter 8 (2012) 9705–9712, https://doi.org/10.1039/c2sm26054a. [43] J. Yin, E. Chen, D. Porter, Z. Shao, Enhancing the toughness of regenerated silk fibroin film through uniaxial extension, Biomacromolecules 11 (2010) 2890–2895, https://doi.org/10.1021/bm100643q. [44] Q. Yuan, J. Yao, X. Chen, L. Huang, Z. Shao, The preparation of high performance silk fiber/fibroin composite, Polymer (Guildf) 51 (2010) 4843–4849, https://doi. org/10.1016/j.polymer.2010.08.042. [45] Q. Yuan, J. Yao, L. Huang, X. Chen, Z. Shao, Correlation between structural and dynamic mechanical transitions of regenerated silk fi broin, Polymer (Guildf) 51 (2010) 6278–6283, https://doi.org/10.1016/j.polymer.2010.10.046. [46] G. Bratzel, M.J. Buehler, Molecular mechanics of silk nanostructures under varied mechanical loading, Biopolymers 97 (2011) 408–417, https://doi.org/10.1002/bip. 21729. [47] M. Cetinkaya, S. Xiao, G. F, Effects of crystalline subunit size on silk fiber mechanics, Soft Matter 7 (2011) 8142–8148, https://doi.org/10.1039/c1sm05470h. [48] M. Cetinkaya, S. Xiao, F. Grater, Bottom-up Computational Modeling of SemiCrystalline Fibers : From Atomistic to Continuum Scale, (2011), pp. 10426–10429, https://doi.org/10.1039/c0cp02936j. [49] M. Cetinkaya, S. Xiao, B. Markert, W. Stacklies, F. Gra, Silk fiber mechanics from multiscale force distribution analysis, 100 (2011) 1298–1305, https://doi.org/10. 1016/j.bpj.2010.12.3712. [50] T. Giesa, M. Arslan, N.M. Pugno, M.J. Buehler, Nanoconfinement of spider silk fibrils begets superior strength, extensibility, and toughness, Nano Lett. 11 (2011) 5038–5046, https://doi.org/10.1021/nl203108t. [51] I. Su, M.J. Buehler, Nanomechanics of silk: the fundamentals of a strong, tough and versatile material, Nanotechnology 27 (2016) 1–15, https://doi.org/10.1088/09574484/27/30/302001. [52] S.A. Adcock, J.A. McCammon, Molecular dynamics: survey of methods for simulating the activity of proteins, Chem. Rev. 106 (2006) 1589–1615, https://doi.org/ 10.1021/cr040426m. [53] J.L. Klepeis, K. Lindorff-Larsen, R.O. Dror, D.E. Shaw, Long-timescale molecular dynamics simulations of protein structure and function, Curr. Opin. Struct. Biol. 19 (2009) 120–127, https://doi.org/10.1016/j.sbi.2009.03.004. [54] W.F. Van Gunsteren, H.J.C. Berendsen, Computer simulation of molecular dynamics: methodology, application and perspectives in chemistry, Angew. Chem. Int. Ed. Engl. 29 (1990) 992–1023, https://doi.org/10.1002/anie.199009921. [55] S. Inoue, K. Tanaka, F. Arisaka, S. Kimura, K. Ohtomo, S. Mizuno, Silk fibroin of Bombyx mori is secreted , assembling a high molecular mass elementary unit consisting of H-chain , L-chain , and P25 , with a 6 : 6 : 1 molar ratio *, J. Biol. Chem. 275 (2000) 40517–40528, https://doi.org/10.1074/jbc.M006897200. [56] C.Z. Zhou, F. Confalonieri, N. Medina, Y. Zivanovic, C. Esnault, T. Yang, M. Jacquet, J. Janin, M. Duguet, R. Perasso, Z.G. Li, Fine organization of Bombyx mori fibroin heavy chain gene, Nucleic Acids Res. 28 (2000) 2413–2419, https://doi.org/10. 1093/nar/28.12.2413. [57] T. Asakura, K. Ohgo, T. Ishida, P. Taddei, P. Monti, R. Kishore, Possible implications of serine and tyrosine residues and intermolecular interactions on the appearance of silk I structure of Bombyx mori silk fibroin-derived synthetic peptides: high-resolution 13C cross-polarization/magic-angle spinning NMR study, Biomacromolecules 6 (2005) 468–474, https://doi.org/10.1021/bm049487k. [58] M. Lee, J. Kwon, S. Na, Mechanical behavior comparison of spider and silkworm silks using molecular dynamics at atomic scale, R. Soc. Chem. 18 (2016) 4814–4821, https://doi.org/10.1039/C5CP06809F. [59] J.M. Gosline, P.A. Guerette, C.S. Ortlepp, K.N. Savage, The mechanical design of spider silks: from fibroin sequence to mechanical function, J. Exp. Biol. 202 (1999) 3295–3303 (doi:JEB2216).

20