Biocompatibility of 2D Silicon Nitride: Interaction at the Nano-Bio interface

Determining potential abilities of nanostructures to induce toxicity to biological molecules is still a convoluted challenge in the realm of nanomedicine. Based on the unprecedented achievements of twodimensional nanomaterials in nearly all areas of applied sciences particularly medicine, we carried out all-atom molecular dynamics simulations to assess the biologically-important, yet-unmapped issue of the biocompatibility of 2D, hexagonal \b{eta}-Si3N4 nanosheet via investigating its possible cross interactions with both human serum albumin (HSA) and p53 tumor suppressor. Examining the conventional MD indicators in the presence and absence of the monolayer revealed that hexagonal Si3N4 nanosheet weakly binds to these two proteins without inducing any important, dramatic change to their secondary structures, revealing accordingly the biological compatibility of the monolayer in case it is released as therapeutics or carriers in vivo. This finding was also broadly supported by the related, time-dependent behaviors of the protein-monolayer as well as the protein-water interaction energies.


Introduction
The stupendous advent of nanoscience and nanotechnology has enabled human being to precisely manipulate matter at one of the deepest levels of reality, the nanoscale, with a pervasive impact on nearly all branches of physical sciences, from materials science, to electronic devices, to biotechnology, leading to the emergence of multifarious interdisciplinary fields at the interface of physics, chemistry, and biology. Among such uncharted territories is nanobiotechnology [1], which has made a categorically great contribution to modern medicine and therefore to the dawn of the nanomedicine epoch via adding functionalities to nanomaterials and then by interfacing them with biological structures, aiming at diagnosing, preventing, and treating diseases on molecular scales as well. These nanomaterials, whether as therapeutics or carriers, inevitably interact with biological entities within human body, and one ultimate goal of such so-called nanobiosystems is then their in vivo applications [2]. To this end, they must accordingly conquer several intractable challenges, which an important of them is the biocompatibility prerequisite, in the sense that nanobiomaterials should not exhibit any toxic or injurious effect on biological systems (cells, proteins, tissues, etc.). Therefore, in vivo toxicological evaluation of nanomaterials is cardinal for exploiting them in nanomedicine.
As yet, sp 2 carbon nanomaterials, particularly fullerenes [3], carbon nanotubes [4], and graphene [5], have excited the ardor of scientists in the area of nanomedicine because they are best suited for drug delivery [6][7][8], sensing biological targets [9,10], biomedical imaging [11], and cancer treatment [12]. Novel 2D layered nanomaterials such as MoS 2 [13], boron nitride [14], WS 2 [15], and graphite-carbon nitride [16] have also captivated many researchers in biosensing and nanomedicine because of their structural similarities to graphene as a paragon of layered nanomaterial. As a result, except silica [17,18], minimal attention has so far been paid to silicon-based nanomaterials in terms of their potential medicinal applications. Indeed, the macroscopic state of Si 3 N 4 has been proved to be biocompatible and stable in vivo, and such properties, when combined with its phenomenal mechanical features [19][20][21], make Si 3 N 4 an intriguing ceramic implant material, being truly useful in some healthcare applications, particularly in orthopedic surgery [22].
A recent investigation on Si 3 N 4 nanostructures carried out in 2020 [23] has unveiled a new, 2D member of this family, showing the yet-continuing promisingness of silicon nitride materials in the post-graphene age. In the same work, it has been proved that β-Si 3 N 4 nanosheets exhibit a semiconducting behavior with a band gap of about 2 eV. Taking into account the findings that semiconducting MoS 2 [24] or graphene [25] monolayers indeed destabilize amyloid beta fibrils [26], we therefore hypothesized that 2D Si 3 N 4 may exhibit the same effect based on the structural and electronic similarities. As a result, the present work was devoted to evaluating the biocompatibility of β-Si 3 N 4 monolayers as a first step in testing our hypothesis in case they are released in vivo. We carried out the present investigation via examining possible cross interactions of the monolayer with two specific homo sapiens proteins, namely human serum albumin (HSA)-as the most abundant protein in human blood plasma, constituting about half of the serum protein-and the p53 tumor suppressor protein described as the guardian of the genome due to its role in conserving genome stability by preventing mutations [27]. We applied all-atom molecular dynamics (MD) simulations, and calculated and examined a number of important conventional MD indicators as described in Sec. 2.

Computational details
Initial atomic positions of HSA and p53 proteins were taken from the RCSB Protein Data Bank (PDB) [28] with the entry codes 3B9M and 1TUP, respectively. All the classical MD simulations were carried out by NAMD (version 2.14b2) [29], in parallel, on Debian-style [30] Linux [31] systems using Open MPI v.3.1.6 [32], with the July 2018 update of CHARMM36 [33,34] force fields. The VMD program (version 1.9.4a9) [35] was also used for post-processing. The two nanobiosystems including HSA-monolayer and p53-monolayer (abbreviated as HSA-ML-W and p53-ML-W; ML for monolayer and W for water) were then solvated in two boxes with TIP3P [36] water with dimensions of 17.8 × 15.5 × 8.7, and 17.8 × 15.5 × 7.5 nm 3 respectively, under periodic boundary conditions with a unit-cell padding of about 2 nm to decouple the periodic interactions, as illustrated in Fig. 1. . The α-helices, β-strands, and random coils turns have been shown in pink, yellow, and cyan white, respectively.
We also solvated the two proteins in the same water boxes (abbreviated as HSA-W and p53-W), this time without Si 3 N 4 monolayer, and accordingly referred to them as the reference (control) trajectories to which the results of the HSA-ML-W and p53-ML-W simulations were compared. The Na + and Cl − ions were randomly distributed (replaced by the same number of water molecules) in order for the entire systems to be electrostatically neutral, as tabulated in Table 1.
The switching and cutoff distances of 1.0 and 1.2 nm were also used for truncating non-bonded van der Waals interactions, respectively. Particle-mesh Ewald (PME) [38] with grid dimensions of 180 × 160 × 88 and 180 × 160 × 80 were applied to systems containing HSA and p53 for long-range interactions, respectively. Four minimization simulations for HSA-ML-W, p53-ML-W, HSA-W, and p53-W were carried out, each for 400000 conjugate gradient steps (0.4 ns). The next four MD simulations were carried out, each for 50 ns, with the integration time step of 1 ns within the NPT ensemble at 310 K and 1.01325 bar using Langevin forces with the Langevin damping constant of 2.5 ps along with the Nosé-Hoover Langevin piston pressure control. For systems without the silicon nitride monolayer (HSA-W and p53-W), the dielectric constant was set to 1. The hydrogen donor-acceptor distance and the angle cutoff have been respectively set to 3 Å and 20 • to identify hydrogen bonds. The solvent-accessible surface area (SASA) [39] was calculated for each system using the rolling-ball algorithm [40] with the radius of 1.4 Å for the probe sphere. The interaction energy (E int = electrostatic+van der Waals energies) between protein-monolayer as well as protein-water in the presence and absence of the monolayer were also calculated as functions of time for each protein in order to check whether stable nano-bio complexes are made.

Force field of β-Si N 4 nanosheet
The initial force field describing the bulk phase of β-Si 3 N 4 was constructed using the values provided by Wendel et al. [41]. The force field associated to 2D Si 3 N 4 was then obtained via calibrating that of the bulk phase in a way that the (molecular-mechanical) dielectric constant of the monolayer obtained by NAMD was fitted to its (quantummechanical) analogue (namely, κ 2.37) calculated using QUANTUM ESPRESSO [42]. To this end, we adopted a self-consistent [43], plane-wave, pseudopotential approach [44] at the PBE-GGA [45] level of density-functional theory (DFT) [46]. Scalar-relativistic ultrasoft pseudopotentials [47,48] (Si.pbe-n-rrkjus_psl.1.0.0.UPF and N.pben-rrkjus_psl.1.0.0.UPF [49]) generated by Rappe-Rabe-Kaxiras-Joannopoulos (RRKJ) [50] pseudization method with nonlinear core correction [51] were used to model the core electrons. The valence shells of the N and Si atomic species were also described by 2s2p and 3s3p orbitals, respectively. The pseudo-wavefunctions and charge density were expanded in a plane-wave basis set with kinetic energy cutoff values of 90 and 360 Rydberg (Ry) respectively, for which energy convergence was optimally achieved. We started from the β-Si 3 N 4 primitive cell of the bulk phase containing 6 silicon and 8 nitrogen atoms in a hexagonal Bravais lattice under periodic boundary conditions with experimental lattice constant of a 0 = 7.82 Å, and then applied a vacuum space of > 15 Å along z to produce the 2D structure, as shown in Fig. 2. The equilibrium (zero-pressure) lattice constant of the nanosheet was also calculated using Murnaghan's isothermal equation of state [52,53], leading to the theoretical, PBE-GGA value of a = 8.262 Å, which is about 5.65% larger than that of the bulk phase (a 0 ).  [54]. The purple parallelogram is the associated unit cell, containing six Si and eight N atoms in a hexagonal Bravais lattice.
We carried out four density-functional molecular dynamics (DFMD) simulations (i.e., two minimization and two finite-temperature free dynamics simulations) within the Car-Parrinello (CP) approach at 310 K, one with applying a zero electric field (E = 0), and the other with E = 16 kcal (mol.Å.e) along +z to change the dipole moment (p) of the monolayer. These CPMD simulations were preceded by electronic minimization processes for each value of E in order to bring the electronic wavefunctions on their ground states relative to the starting atomic configurations. The two setups (E = 0 and E = 16) were minimized after 100 damped-dynamics steps (0.012 ps) with the time-step of ∆t = 0.12 fs and with the electron damping value (= damping frequency times ∆t) of 0.1. The fictitious electron mass in the CP Lagrangian [55] was set to 1000 a.u. (1000 times the rest mass of electron) to guarantee the validity of the adiabatic approximation [56]; the mass cutoff of 2.5 Ry was also chosen for the Fourier acceleration effective mass to keep the quality of simulations from being adversely affected, as well as to minimize the electron drag effect. The four DFMD simulations were started from the same minimized structure in terms of ionic degrees of freedom, in that the value of the total force exerted on each atom eventually became < 10 −4 eV Å. The propagation time was also chosen about 0.36 ps (3000 Verlet steps). Electronic equations of motion were also accelerated using a preconditioning scheme [57]. We ignored at least the first 0.12 ps (1000 steps) of the simulations for thermalization and reliable statistical averaging. Both electronic and ionic contributions were taken into account in estimating the average value of the total dipole moment. The dielectric constant was also calculated using where ∆p = |p (E=0) − p (E=16) |, 0 = 2.398 × 10 −4 mol.e 2 (kcal.Å) is the vacuum permittivity, and Ω is volume of the monolayer.  Fig. 3(a), it is seen that the HSA structure has remained stable during the simulations with all-atom RMSDs of < 0.5 and 0.6 nm from the reference crystal structure (3B9M) on interaction with the monolayer (HSA-ML-W) and in aqueous lonely (HSA-W), respectively. The RMSD of HSA-ML-W also takes smaller values compared to the other over the last 35 ns, demonstrating that interaction with Si 3 N 4 nanosheet considerably reduces thermal fluctuations of HSA due to binding to the nanosheet, and accordingly decreases the protein's conformational change more than that of HSA-W.

Results and discussion
Time dependence of the radius of gyration illustrated in Fig. 3(b) also reveals the fact that the overall structural compactness of HSA on interaction with the monolayer is slightly (by ∼ 1.0 Å) smaller than those of HSA-W or the initial crystal structure. The per-residue RMSF [ Fig. 3(c)] averaged over all frames further indicates that interaction with the nanosheet considerably decreases the flexibility of all regions of HSA, and therefore makes its secondary structure resistant to any change raised by thermal fluctuations in aqueous. Consistently, the number of hydrogen bonds within the protein in both HSA-W and HSA-ML-W [ Fig. 3(d)] exhibit no remarkable change over the entire trajectory compared to each other. More precisely, the number of hydrogen bonds averaged over the last 45 ns is about 147 for HSA-W and 139 for HSA-ML-W, indicating a decrease as negligible as 5.5% on interaction with the monolayer. The calculated SASA [ Fig. 3(e)] over the last 40 ns of the two trajectories also show convergence in a way that HSA-ML-W takes larger values compared to HSA-W on average, in agreement with the gyration radius [ Fig. 3(b)]. As a result, HSA strongly binds onto the surface of Si 3 N 4 monolayer, forming a stable complex.
The secondary structure of HSA as a function of time has also been provided in Fig. 4 in the presence and absence of the silicon nitride monolayer.  E, B, H, G, I, and C also stand for turn, extended configuration, isolated bridge, α-helix, 3 10 -helix, π-helix, and coil, respectively. The α-helix-rich structure (the predominant, violet areas) of HSA is unmistakably clear.
Consistent with the previous analyses, no dramatic change is accordingly observed, and the α-helix-rich structure of HSA [ Fig. 1(a)] is clearly preserved on interaction with the monolayer.
Examining the corresponding Ramachandran plots also confirms the preceding observations as illustrated in Fig. 5. That the RMSD of p53-ML-W takes smaller values over the entire time-span reveals that interaction with 2D Si 3 N 4 dramatically decreases fluctuations of the protein caused by thermal energy, therefore, abates the tendency toward any conformational change.
Examining the time dependence of radius of gyration [ Fig. 6(b)] further indicates that the overall structural compactness of p53 in the presence and absence of the monolayer are nearly the same over the entire trajectory. The RMSF curves, illustrated in Fig. 6(c), show a considerable decrease in the fluctuating pattern values by ∼ 0.7 nm, and therefore in the overall flexibility of the protein on interaction with the monolayer. The peak in the middle (at residue 209) of the p53-ML-W curve is for the arginine residue, which exhibits a high degree of flexibility at this point based on the fact that it is located on a turn secondary structure with a distance of about 2.8 nm from the monolayer, indicating no effective binding between them as well.
Interaction with Si 3 N 4 nanosheet has also no (significant) impact on the number of hydrogen bonds within p53, and consequently on the secondary structure of p53 as seen in Fig. 6 . Indeed, such a discrepancy is not so important as to considerably change the secondary structure of the protein.
The secondary structure of p53 has been illustrated in Fig. 7 in the presence and absence of the silicon nitride monolayer over the entire trajectory. Nearly the same, β-rich patterns is observed in the two subfigures consistent with the tertiary structure of p53 [ Fig. 1(b)], demonstrating the minimal, insignificant change in the related secondary structure on interaction with the monolayer.
The associated Ramachandran plots illustrated in Fig. 8 are also in agreement with the previous findings. In the same way as HSA, but to a higher extent, the distributions of the dihedral angles (dominantly on the extended region) over different areas are very close to each other comparing p53-W and p53-ML-W at t = 0 and 50 ns, and comparing p53-W and p53-ML-W at t = 50 ns. As a result, no considerable change in the related secondary structure could be observed in case it is in contact with the monolayer.
We finally estimated the time dependence of the interaction energy (E int ) between protein-monolayer on one hand, and protein-water in the presence and absence of the monolayer on the other, as illustrated in Fig. 9. The E int values for p53-and HSA-monolayer are respectively about −151 and −510 kcal mol, which are dramatically larger (therefore weaker bindings) than those of the protein-water interactions, showing that these biological proteins do not make stable complexes on interaction with the monolayer. That HSA has the lower value is also an indication of the fact that it is a considerably larger protein compared to p53, which accordingly led to larger SASA values comparing Figs. 3(e) and 6(e).
Interaction with the Si 3 N 4 nanosheet also decreases the p53-water binding by about 16.7% from −3574 (in p53-W) to −2977.7 kcal mol (in p53-ML-W). In contrast, the presence of the monolayer increases the HSA-water binding by 1.85% from −18201.4 (in HSA-W) to −18537.7 (in HSA-ML-W). However, non of these percent values is significant so as to dramatically affect the protein-water bindings.

Conclusions
All-atom molecular dynamics simulations were applied to investigate the biocompatibility of 2D, hexagonal β-Si 3 N 4 monolayer via examining its possible impacts on both HSA (human serum albumin) and p53 antitumor protein. We accordingly calculated and examined a number of important MD indicators including RMSD, radius of gyration, per-residue RMSF, number of hydrogen bonds, solvent-accessible surface area, and secondary structure of each protein in a contrasting fashion in the presence and absence of the monolayer. Results verified that the secondary structures of these proteins remain nearly intact on interaction with Si 3 N 4 nanosheet. Examining the associated protein-monolayer and protein-water interaction energies in the presence and absence of the monolayer further revealed that these biological proteins do not make stable complexes with the monolayer. The presence of Si 3 N 4 also affected both HSA-and p53water bindings very marginally. It was accordingly inferred that hexagonal β-Si 3 N 4 nanosheet is indeed a biocompatible material and could then be used as a therapeutic or carrier for in vivo applications.

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.

Data Availability Statement
The data supporting the findings of the present investigation are available on reasonable request from the corresponding author.