Effect of Nanohydroxyapatite on Silk Fibroin-Chitosan Interactions. Molecular Dynamics Study

1. Introduction

It is widely recognized that the interactions between proteins and polysaccharides are pivotal in both living organisms [1] and the development of novel biomaterials [2]. Among the polysaccharide analogues, proteoglycans deserve special attention. These biomolecules are composed of a protein and glycosaminoglycan (GAG) moieties bound to the protein by the serine and threonine fragments [3]. They are typically found on the cell surface and in the extracellular matrix (ECM). Proteoglycans are essential for the strength and flexibility of cartilage and tendons. The absence of proteoglycans is associated with heart and respiratory failure, defects in skeletal development, and increased tumor metastasis [4].

Silk threads are considered to be one of the strongest natural fibers. They exhibit a high tensile strength of 4.8 GPa, very good flexibility [5,6,7] and high thermal stability [8,9,10]. Silk is more resistant to loads than Kevlar (the most potent manufactured fiber to date) [7,11,12]. Spiderweb silk exhibits exceptional mechanical properties [10,13]. However, silk produced on an industrial scale is primarily sourced from domesticated mulberry silkworms (B. mori) due to their ease of cultivation and higher thread yield from their cocoons compared to spider webs [10,13]. Natural silk thread is composed of two tightly packed bundles of fibroin fibers glued with a sericin sheath, soluble in hot water. Fibroin protein is a dominant component of raw silk is, accounting for 70% to 80% of its weight [14,15,16]. The other components are sericin, wax, fats and dyes [12,17]. Fibroin, has a high content of amino acids, with serine comprising approximately 12% [18], whereas collagen contains four times less of it [19]. The repetitive sequence of amino acids (alanine and glycine) already in the primary structure makes silk fibroin a highly ordered protein. An essential repeating motif of the fibroin amino acid sequence is –(GLY-SER-GLY-ALA-GLY-ALA)_n. These chains, connected by hydrogen bonds, form the antiparallel β-sheet structure. Due to the periodic primary structure, there are glycine residues in the secondary structure on one side of the surface of the β-sheet, interacting with the glycine residues of the adjacent chain. Similarly, on the other side of the surface of the β-sheet, there are serine and alanine residues interacting with similar residues of another neighboring chain (Figure 1). This quasi-crystalline structure is interrupted by stretches of sequence containing other amino acids (e.g. TYR, ASP, ARG and VAL).

Silk fibers with a high proportion of regular glycine-alanine fragments tend to be more flexible due to weak interactions between chains in the β-sheet structure, but less stretchable as the polypeptide chains are under significant tension in this secondary structure. Fiber properties such as gloss and softness depend on differences in the size and characteristics of amino acid residues within irregular structural fragments. Crystalline regions within fibers, marked as A in Figure 1, consist of tightly packed motifs of the antiparallel β-sheet, aligned along the fiber axis. Dispersion angle measurements revealed the presence of crystallites with dimensions of 2 x 2 x 2 nm [20]. However, these dimensions increase to 5-7 nm in length at slower reeling speeds [23]. Additionally, the analysis identified non-crystalline regions with two distinct morphologies: (i) oriented regions with an α-helix structure (B , Figure 1) and (ii) non-oriented amorphous areas (C, Figure 1). A closer insight into the less-ordered regions aligned with the fiber axis reveals the presence of the 31-helix motif (which grants the fibers elasticity) and a greater abundance of β-bends/turns and potentially β-spirals. These β-structures, prevalent in spider silk, might be responsible for the rapid shrinking of the thread in water (super contraction effect) [24,25,26].

Silk fibroin consists of two main chains: a heavier H-chain (around 300 kDa in B. mori) with hydrophobic nature and a lighter L-chain (25 kDa) with hydrophilic nature. These chains are linked by a disulfide bridge. The structure also includes a chaperone protein, P25, which is a fibrohexamer glycoprotein. The structure of the silk thread is supported by strong hydrogen bonds, especially in the crystalline areas, and by Van der Waals interactions. According to Kaplan models, the primary structure's chain is predominantly hydrophobic, interspersed with small repetitive sections of a hydrophilic nature [27]. These alternating hydrophilic and hydrophobic interactions lead to the formation of micelles, where water is bound within the hydrophilic domains (approximately 8%). These micelles, measuring up to 200 nm in size, ultimately form fibroin globulins. When subjected to various physical or chemical stimuli such as changes in pH, temperature fluctuations, variations in metal ion concentration, increased fibroin concentration, ultrasound exposure, or mechanical shearing, disordered and β-helical regions transform into ordered β-sheet structures. Concentrated globulins transition into a liquid crystal phase, with particles aligning parallel to one another as they flow, solidifying upon exposure to air to form a thread.

Due to its hydrophobic nature, and the presence of inter- and intramolecular strong hydrogen bonds in many crystalline regions, natural silk fibers dissolve only in a few solvents including concentrated acids: H₃PO₄, HCOOH, H₃SO₄, and HCl and aqueous solutions of concentrated salts: LiSCN, LiBr, CaCl₂, Ca(SCN)₂, ZnCl₂, NH₄SCN, CuSO4, Ca(NO₃)₂ [17]. Despite the solubility limitation, it is possible to obtain many different forms of silk material, e.g., films, foams, hydrogels, powders, nanofibers, and microcapsules [18]. Silk is used not only in the clothing industry for producing fabrics, anti-allergic bedding, bulletproof vests, ropes, and parachutes but also in medical applications. Silk fibers (catgut) have been used as surgical threads for years. There is a lot of interest in an aqueous solution of silk fibroin, which, thanks to its biocompatibility, biodegradability, and similarity to collagen and elastin, is successfully used in tissue engineering. Silk fibroin is a promising scaffold material for cell growth and regeneration in various damaged tissues, including connective, nervous, cartilage, and bone tissues [7]. Furthermore, it has a wide range of other biomedical and pharmaceutical applications including drug carriers [28,29,30] and hydrogels for transparent wound dressings [31,32]. Noteworthy, additional beneficial properties, such as antibacterial properties, increased strength, resistance to degradation, viscosity suitable for bioprinting, and others, can be achieved by adding other polymers such as chitosan, hyaluronic acid, cellulose, sodium alginate, starch, carrageenan, collagen, gelatin, PVA, PVP, poly(lactic acid) [33,34].

Among polysaccharide analogues, chitosan demonstrates exceptional biocompatibility and biodegradability [35,36]. It is derived from chitin, an abundant and renewable resource, and is widely acknowledged for its safety profile [36]. Its versatility and potential or actual application has been widely demonstrated in numerous fields including agriculture [37,38], pharmacy and medicine [30,36,37,39,40,41,42,43], food industry [36,37,44], textile industry [36,37,45,46], biotechnology [47], cosmetics [37,48,49,50,51], paper production [37,52] and various green and sustainable technologies [53,54,55]. Chitosan emerges as a highly desirable substance, owing to its non-toxicity, cationic nature, biocompatibility and antibacterial properties, alongside its unique chemical composition [35].

Our previous research employed molecular dynamics simulations to explore complexes of chitosan/chitin with collagen [56,57,58] and hyaluronan [59]. Understanding these interactions is crucial for developing strategies for cartilage regeneration [60]. Likewise, chitosan-fibroin composites are relevant form the biomedical and pharmaceutical perspectives [34,61,62,63,64,65]. Notably, chitosan/fibroin composites containing hydroxyapatite exhibit exceptional tissue regeneration potential due to their unique properties [61,63]. This study focuses on unraveling the intermolecular interactions between fibroin and chitosan in the presence of HPO₄²;⁻ and Ca²;⁺ ions using the molecular dynamics simulations. The main goal is to elucidate the factors that contribute to the remarkable stability of these composites.

2. Methods

All molecular dynamics computations were carried out using an YASARA software [66,67], a comprehensive and versatile tool for modeling biopolymer complexes including polysaccharide/peptide systems [56,57,58,59,68,69,70]. The structure of the chitosan was taken from PubChem and modified to obtain polymer of total mass of 200 kDa (1050 monomers). Structure of silk fibroin heavy chain was obtained from https://alphafold.ebi.ac.uk/entry/B0FRJ4 (Fibroin heavy chain Fib-h). In order to optimize structure further it was equilibrated for 5 ns. Additionally, a single chitosan chain was included in the simulation box, and five separate simulations, each lasting 40 ns, were conducted. Finally, 8 fibroin chains were added to chitosan for each system. Two cases have been studied: 1) with the presence of Ca²⁺ and HPO₄^2- ions (HPO₄²;⁻/Ca²;⁺) and 2) with the presence of acetic acid. Point charges for silk fibroin atoms were assigned using the AMBER14 force field [71], while charges for chitosan atoms were assigned using the GLYCAM06 force field [72].

In the subsequent step all molecular systems were used for the molecular dynamics computations. The applied approach included the water solvation effects and acid-base characteristics. For this reason, the procedure encompassed the adjustment of hydrogen bonds structural features to achieve appropriate protonation microstates of the peptide at pH=7.0 [73]. After equilibration and energy minimization, molecular dynamics (MD) simulations were run for 1 nanosecond. The simulations employed appropriate force fields for all atoms in the studied molecular systems: AMBER14 for fibroin, GLYCAM06 for chitosan, and TIP3P for water. Standard AMBER settings were used, including a 1.0 nm van der Waals cutoff and the Particle Mesh Ewald algorithm [74] for calculating electrostatic interactions. The equations of motion were integrated using time steps of 1.25 femtoseconds and 2.5 femtoseconds for bonded and non-bonded interactions, respectively, under the conditions of temperature (T) set at 310 K and pressure (P) at 101325 Pa within the NPT ensemble [74]. The simulation duration time was set to 200 ns to ensure the sufficient information for subsequent analysis. Achieving the optimal structure of the complex was confirmed based on several structural parameters including RMSD off all atoms evolution, radius gyration and numbers of intermolecular interactions (supplementary Figures S1-S2). In order to determine binding energy and other critical characteristics, 1500 data points were meticulously chosen from the 50-200 ns timeframe subsequent to root mean square deviation/displacement (RMSD) analysis.

3. Results and Discussion

This study examined two chitosan-silk fibroin complexes. One complex incorporated Ca²⁺ ions and HPO₄²⁻ to mimic nanohydroxyapatite, while the other contained acetic acid. The selection of acetic acid environment, as a reference simulation conditions was guided by two considerations: its well-documented role as a stabilizer for fibroin [62] and its ability to dissolve both fibroin and chitosan [62,75]. The exemplary optimized structures of chitosan-silk fibroin complexes are presented in Figure 2. The detailed structural and geometric parameters, including hydrogen bond lengths, interaction distances, and binding energies are provided in the supplementary materials. As anticipated, the basic amino groups in chitosan undergo protonation in the presence of acetic acid, as well as in the presence of HPO₄^2- ions. The role of charged ions in stabilizing these large biomolecular systems presents an intriguing aspect for investigation. Molecular dynamics simulations reveal that the formation of complexes with both CH₃COOH and HPO₄²;⁻/Ca²;⁺ is energetically favorable. However, the total binding energy, encompassing all stabilizing interactions, is slightly higher in the former system (7452 ± 397 kJ/mol) compared to the latter system (6483 ± 561 kJ/mol). Despite these differences, both values are exceptionally high. Notably, the binding energies calculated using similar procedures for chitosan-collagen complexes, reported in our previous studies [56,58], did not exceed 850 kJ/mol. This suggests a significantly lower mutual affinity between chitosan and collagen compared to chitosan-fibroin systems.

Similarly to chitosan, fibroin undergoes protonation in acidic conditions [76,77], as evidenced by computational results. While both biomolecules become positively charged through protonation, the formation of stable complexes is possible due to the presence of oppositely charged ions that neutralize the positive charge. Our previous studies [58] have demonstrated a similar behavior in case of collagen-chitosan systems. Importantly, both chitosan and proteins such as collagen and fibroin exhibit polyelectrolytic properties [78,79,80], leading to charge migration within the polymer chain and fluctuations in charge density within the biomolecule. These properties can be manifested in the complex molecular systems. Nevertheless, a similar calculations performed in the previous studies [58] have shown that even oligomeric chitosan chains exhibit polyelectrolytic features. Since polyelectrolytes contain many ionic and polar groups, they are capable of forming various intermolecular contacts. Remarkably, the complexes investigated in this study reveal a dense network of interactions between chitosan and fibroin. The YASARA program employed in this study was capable of identifying various primary types of interactions in the considered fibroin-chitosan complexes, such as ionic contacts, hydrogen bonds, and hydrophobic interactions. Of particular significance are the strong ionic contacts These interactions are known to be crucial for stabilizing the fibroin-chitosan complexes [81]. As shown in Figure 3, in case of both complexes, containing acetic acid and HPO₄²;⁻/Ca²;⁺ ions various ionic contacts between chitosan and fibroin appears. However, a significantly higher number of these interactions are observed under the influence of HPO₄²;⁻/Ca²;⁺ ions compared to the acetic acid solution. Similarly, Figure 4 reveals a greater frequency of hydrogen bonds between fibroin and chitosan in the presence of HPO₄²;⁻ and Ca²;⁺ ions. Noteworthy, hydrogen bonds and electrostatic interactions in fibroin blends containing different polysaccharides including chitosan were analyzed using experimental tools [81]. As indicated by these studies, interactions between the protein and polysaccharides can indeed modify the protein structure itself. From the perspective of bone regeneration, it seems to be beneficial to utilize materials structurally similar to natural bone tissue. Fibroin emerges as a promising candidate in this context. It is well established that the crystal phase formed during bone mineralization exhibits a characteristic inter-crystalline center spacing of 0.362 nm [82]. Fortunately, this spacing aligns with the β-sheet conformation observed in fibroin [17]. According to the results of infrared spectroscopy measurements, the complexation of fibroin with chitosan favors the β-sheet conformation of the peptide [81,83,84]. The β-sheet order was also found in the case of the complexes optimized in this study as visualized by the red band (Figure 2). As it can be inferred form Figure 2 and 3, the contribution of the ionic contacts and hydrogen bond between chitosan and β-sheet fragments is quite significant comparing to the contacts formed by other fragments of the peptide (coil, turn). Furthermore, these interactions involve diverse positions in the primary structure of peptide. The effective arrangement of the chitosan chains with respect to the fibroin fragments characterized by the β-sheet order is shown in Figure 2.

The visualization of all main types of interactions formed in the considered systems (ionic contacts, hydrophobic contacts and hydrogen bonds) was presented on Figure 5 and Figure 6. The optimized structures of fibroin-chitosan complexes are stabilized by strong hydrogen bonds and ionic interactions, which are characterized by short interatomic distances ranging from 1.85 to 2.11 Å for hydrogen bonds and from 3.83 to 3.92 Å for ionic contacts. As it was established, when HPO₄^2- ions are present in the system, the ionic bridges can be formed. These interactions involve protonated amino groups in chitosan and LYS fibroin residues (Figure 5). Interestingly, the formation of analogical systems involving CH₃COO^- ions was not observed.

The analysis of molecular contacts involving Ca²⁺ seems to offer a valuable insight into the stabilization of potential biomaterials used for bone regeneration. It is well known that, the addition of calcium salts is advantageous for the solubilization of fibroin [17], indicating Ca²⁺/peptide interactions. Indeed, several types of ionic bridges (ASP-Ca²⁺-ASP, ASP-Ca²⁺-GLU, ASP-Ca²⁺-GLY, GLU-Ca²⁺-GLU, GLU-Ca²⁺-GLY) were observed in the case of optimized molecular structures of fibroin-chitosan complexes. An exemplary contacts of this type were presented on Figure 6c and 6d. Noteworthy, an experimental study by Koeppel et al. [85] confirms that the carboxylic groups in ASP and GLU in silk peptides form ionic contacts with Ca²⁺ and K⁺, which has a significant impact on their rheological properties.

To identify hydrogen bonds, default YASARA protocol based on energy calculation with a specified threshold of −6.25 kJ/mol was applied [86]. The molecular dynamics simulations revealed that ALA, ARG, ASN, ASP, GLN, GLU, GLY, LEU, PRO, SER, THR, TYR, and VAL form the strongest hydrogen bonds, whether the complexes are stabilized by CH₃COOH or Ca²⁺/HPO₄²⁻. The highest binding energy values were observed in the case of ARG, reaching 22.30±3.45 kJ/mol for acetic acid-stabilized complexes and 19.49±5.86 kJ/mol for systems stabilized by Ca²⁺ and HPO₄²⁻ ions. Among these hydrogen bonds, interactions formed by ASP were found to be the most significant contributors, constituting approximately 26.3% of all hydrogen bonding interactions. A considerable contribution to complex stabilization was also noted for SER residues. In the case of complexes stabilized by acetic acid, hydrogen bonds involving SER exhibited a binding energy of 18.16±5.45 kJ/mol, accounting for 10.5% of all hydrogen bonds, while those stabilized by Ca²⁺ and HPO₄^2- ions accounted for 10% (18.36±5.56 kJ/mol). Additionally, THR residues were observed to form relatively strong hydrogen bonds, with binding energies of 17.72±5.91 kJ/mol (CH₃COOH) and 18.80±5.60 kJ/mol (Ca²⁺/HPO₄^2-). Noteworthy, both SER and THR play an important structural role in silk fibroin. Since SER and THR residues contain active hydroxyl groups, they serve as attachment sites for carbohydrates. Furthermore, SER serves as the phosphorylation site [21]. In hydrogen bonding, SER and THR can act as both hydrogen donor and acceptor. Moreover, SER is often found in tight bends or loops, faces outward, forming hydrogen bonds with molecules of the medium (water) [21]. THR residue, which is more hydrophobic than SER, often occurs in an antiparallel β-sheet, facing the inside of the molecule [21].

Unlike ionic interactions and hydrogen bonds, hydrophobic contacts exhibit less significant differences in distribution between complexes stabilized by HPO₄²⁻/Ca²⁺ and acetic acid (Figure 7). Nevertheless, their contribution to the complex stabilization is much less important comparing to previously discussed short-distance interactions. According to the applied default criteria [86] used for the hydrophobic contacts detection in the YASARA software, this type of weak contacts involve carbon atoms surrounded by three or two hydrogen or carbon atoms, or C-H moieties in carbon rings [87]. As it was expected, the majority of hydrophobic contacts are formed by aromatic residues in fibroin chain (PHE and TYR, supplementary Table S4). The interatomic distances determined for these interactions range from 4.33 to 4.54 Å (CH₃COOH) and from 4.25 to 4.66 Å (Ca²⁺/HPO₄²⁻).

It is worth noting that, in addition to the aforementioned main types of interactions, a relatively weak but identifiable cation-π contacts can be also distinguished. The distances determined for these contacts are relatively large and ranged from 3.80-4.19 Å (in the case of the complexes containing Ca²⁺ and HPO₄²⁻) and from 3.91-4.23 Å (in the case of the complexes containing CH₃COOH). These interactions involve the positively charged ammonium group in chitosan and the moieties in fibroin, containing π electrons Our analysis reveals that HPO₄²;⁻/Ca²;⁺-influenced complexes involve cation-π contacts formed by HIS (11.0%), PHE (39.3%), and TYR (49.7%). Interestingly, complexes formed with CH₃COOH exhibit interactions with a wider range of residues, including ASP (7.0%), GLU (4.1%), GLY (0.2%), HIS (4.3%), PHE (27.3%) and TYR (57.1%).

4. Conclusions

The study focused on the intermolecular interactions in fibroin-chitosan complexes in the presence of Ca²⁺ and HPO₄²⁻ ions, which simulate nanohydroxyapatite and in the presence of acetic acid. As we found, molecular dynamics simulations performed using the YASARA program revealed a chemically intuitive description of the structural features and interactions between the studied biomolecules. Our analysis showed that the formation of fibroin-chitosan complexes is highly energetically favorable. Furthermore these systems are much more stable than the previously studied chitosan-collagen complexes. In order to explain the exceptionally high affinity of fibroin for chitosan, the distributions and structural features of the key interaction types were analyzed. All molecular systems considered in this study were found to have a dense network of intermolecular contacts. Additionally, the optimized structures analysis indicated that the complexation of fibroin with chitosan maintains the β-sheet conformation. These computational insights are consistent with empirical observations regarding the favorable mechanical properties and high stability of fibroin-chitosan blends.

Molecular dynamics simulations of biomolecule pairs offer valuable insights into the properties of new materials. The promising results obtained for fibroin-chitosan complexes in the presence of Ca²⁺/HPO₄²⁻ and acetic acid suggest the need for further research on the effects of other ions or molecules. These studies could lead to the development of novel fibroin composites with tailored properties for a variety of applications.