# Effect of Nanohydroxyapatite on Silk Fibroin-Chitosan Interactions--Molecular Dynamics Study
## 2 Methods  
All molecular dynamics computations were carried out using an YASARA software [78; 79], a comprehensive and versatile tool for modeling biopolymer complexes, including polysaccharide/peptide systems [61; 62; 63; 64; 80; 81; 82]. The structure of the chitosan was taken from PubChem (PubChem Chem CID: 129662530) and modified to obtain polymer of total mass of 200 kDa (1250 monomer units). The deacetylation degree of the chitosan applied in this study was 100%. Notably, chitosan with a high degree of deacetylation is commonly used in fibroin composites [83; 84; 85; 86; 87]. A high degree of deacetylation results in a higher number of protonated amino groups in aqueous solutions, thus enhancing solubility and biocompatibility [88]. The structure of the silk fibroin heavy chain was obtained from [[https://alphafold.ebi.ac.uk/entry/B0](https://alphafold.ebi.ac.uk/entry/B0) FRJ4]([https://alphafold.ebi.ac.uk/entry/B0](https://alphafold.ebi.ac.uk/entry/B0) FRJ4) accessed on 6 May 2024 (fibroin heavy chain Fib-h). In order to optimize the 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 Ca2+ and HPO42- ions (HPO42-/Ca2+) and (2) with the presence of acetic acid. Point charges for silk fibroin atoms were assigned using the AMBER14 force field [89], while charges for chitosan atoms were assigned using the GLYCAM06 force field [90]. Notably, studies employing the AMBER14 force field in peptide simulations frequently yielded results aligned with experimental observations [91; 92; 93]. This force field was modified compared to previous versions to better reproduce structural data obtained from NMR measurements [94]. The GLYCAM06 force field, utilized in this study alongside AMBER14, has shown a good quality of prediction of structural features in polysaccharides [95; 96]. 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 the structural features of hydrogen bonds to achieve appropriate protonation microstates of the peptide at pH = 7.0 [97]. The neutral pH was set for all systems considered in the study. The charge of the generated structures used for molecular dynamics simulations was set to zero. 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 [98] 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 101,325 Pa within the NPT ensemble [98]. The [PERSON] barostat and thermostat [99], featuring a relaxation time of 1 fs, were employed to ensure constant pressure and temperature. To mitigate the impact of rescaling, YASARA avoids utilizing the highly variable instantaneous microscopic temperature and pressure for velocity rescaling in each simulation step (a characteristic of the classic [PERSON] thermostat). Instead, it computes the scaling factor based on [PERSON]'s formula using the time-averaged temperature. The simulation duration time was set to 200 ns to ensure sufficient information for subsequent analysis.