Nanotoxicity of two-dimensional nanomaterials on human skin and the structural evolution of keratin protein

Two-dimensional (2D) materials have been increasingly widely used in biomedical and cosmetical products nowadays, yet their safe usage in human body and environment necessitates a comprehensive understanding of their nanotoxicity. In this work, the effect of pristine graphene and graphene oxide (GO) on the adsorption and conformational changes of skin keratin using molecular dynamics simulations. It is found that skin keratin can be absorbed through various noncovalent driving forces, such as van der Waals (vdW) and electrostatics. In the case of GO, the oxygen-containing groups prevent tighter contact between skin keratin and the graphene basal plane through steric effects and electrostatic repulsion. On the other hand, electrostatic attraction and hydrogen bonding enhance their binding affinity to positively charged residues such as lysine and arginine. The secondary structure of skin keratin is better preserved in GO system, suggesting that GO has good biocompatibility. The charged groups on GO surface perform as the hydrogen bond acceptors, which is like to the natural receptors of keratin in this physiological environment. This work contributes to a better knowledge of the nanotoxicity of cutting-edge 2D materials on human health, thereby advancing their potential biological applications.


Introduction
Skin, the largest organ in the human body, acts as the body's outermost defense barrier against external threats and excessive water loss [1].Skin-related issues, such as skin cancer, epidermolysis bullosa, cuts and burns, etc, can have a Nanotechnology Nanotechnology 35 (2024) 225101 (11pp) https://doi.org/10.1088/1361-6528/ad2c58* Authors to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.detrimental influence on patients' physical and emotional well-being [2].The topical cosmetics and personal-care products industry has evolved into the most inventive and appealing markets, propelling the global economy.The revenue in the beauty and personal-care business is predicted to reach US$571.10 billion in 2023, with a 3.80% annual rise [3].Aside from natural ingredients like polyphenols in cosmetics and drugs, two-dimensional (2D) materials have also been added to enhance the efficacy of these personal-care products due to their exceptional physicochemical properties, such as excellent electronic and thermal conductivity, extensive surface area, strong mechanical strength, and impressive optical transmittance [4,5].Particularly, pristine graphene (PG) and graphene oxide (GO), two prototypes of 2D materials, have long been widely applied in cosmetics, drug delivery, biosensors, optics, and nanocomposites [6][7][8][9], to name a few.In addition, their functionalized composites, also known as graphene-based materials (GBM), have also been explored and applied in the fields above [10].PG consists of a single layer of sp2 uncharged carbon atoms forming 6-membered rings in a honeycomb structure, resulting in its hydrophobicity.GO, however, is a highly oxidized form of chemically modified graphene and can be produced by oxidizing graphite in highly acidic conditions.GO has more excellent water dispersion capabilities than PG because of its randomly placed hydroxyl and epoxy groups on both sides of the carbon atoms basal plane and arbitrarily grafted carboxy on the edge of the graphene basal plane [11].It is reported that these charged groups enable GO to possess an amphiphilic character and superior biocompatibility than PG nanomaterials [12,13].
Despite the fact that the multifarious uses of 2D materials have the potential to revolutionize biomedicine, it is critical to address safety issues before their large-scale applications [14].Once administered in the human body, 2D materials have the potential to interact with various biological macromolecules such as proteins and DNA [15].Li et al [16], revealed the effects of nanomaterial defects in nanotoxicity on YAP65WW protein by molecular dynamics (MD) simulations.Mukhopadhyay et al [17] demonstrated the orientationdependent nanotoxicity of black phosphorene on HIV-1 integrase via MD simulations.Tu et al [18], found that PG and GO present antibacterial performance by inducing the breakdown of both the inner and outer cell membranes of Escherichia coli (E.coli).As reported, 2D materials showed antioxidant properties, which could be used in skin whitening products.For example, Qiu et al [19] demonstrated that graphene oxide and few-layer graphene display substantial antioxidant activity, effectively shielding various biomolecules from oxidation.However, Jayakumar et al [20] reviewed the application of 2D materials in healthcare products and emphasized their cytotoxicity and oxidative stress because they are difficult to remove from the human body.The nanotoxicity of 2D materials on disease-related proteins can also be utilized to treat some diseases, such as Alzheimer's disease (AD) and Parkinson's disease (PD).AD stands out as one of the most prevalent neurodegenerative conditions.A key factor contributing to the onset and progression of AD is the aggregation of amyloid, specifically formed by Aβ 1-42 peptides through self-assembly in the brains of AD patients.Yang et al [21] discovered that GO has the ability to extract Aβ 1-40 monomers from Aβ 1-40 aggregates.This extraction process is driven by hydrophobic interactions and π-π stacking interactions.GO effectively disrupts the mature Aβ aggregation and hinders the accumulation of Aβ 1-40 monomers, thereby alleviating oxidative stress injury caused by Aβ 1-40 aggregates to living neuroendocrine PC12 cells.Exploring the application of 2D materials as therapeutic agents for interacting with disease-related proteins represents a novel and intriguing research frontier.Currently, most inquiries in this field are in the theoretical phase, aiming to understand the nanotoxicity mechanisms of 2D materials on disease proteins.This foundational understanding is crucial for the subsequent advancement of animal studies in this area, paving the way for future positive developments in the field.
All these potential risks restrict the broad implementation of 2D materials in biomedicine.To better understand the interfacial interactions between 2D materials and skin, we have employed MD simulations to investigate its molecular mechanisms.The stratum corneum, the outermost layer of human skin epidermis, is made of dead corneocytes (bricks) and lipid matrix (mortar) in the type of bricks-and-mortar structure.As a result, the transdermal permeation and adsorption of drug molecules could be achieved by the transcellular corneocytes or the tortuous lipid pathway, depending on their physiological features, such as size and hydrophobicity [22,23].Keratin is the most abundant cytoskeletal protein in skin corneocytes, feathers, and wools, due to its unique hierarchical structures, which starts with monomer to heterodimer, tetramer, octamer, and unit filament precursor, and finally constructed in the shape of coiled-coil.Recently, keratin and its complexes have been a research frontier in tissue engineering, drug carriers, and wound healing [24].Therefore, it is inevitable to investigate the nanotoxicity of 2D materials on keratin before their wide applications in biomedicine.
Therefore, we aim to investigate the nanotoxicity of two types of typical 2D materials (PG and GO) on skin keratin at the atomistic scale via MD simulations.Because of the limitations of the cost and operating conditions, experimental approaches make it difficult to investigate the nanoscale mechanisms.MD simulations are now routinely used to assess the nanotoxicity of 2D materials, including PG, GO, MoS 2 , and hexagonal boron nitride [18,[25][26][27][28]. Importantly, this less expensive but more efficient computing technique could provide more detailed data than conventional experiments [29].Note that the secondary structure of a protein refers to local folded structures within a polypeptide caused by interactions between atoms of the backbone.It is closely related to protein's higher structure and physiological properties [30].Previous studies show the adsorption of proteins on 2D material surfaces may perturb the protein's native conformation existing in the intracellular and extracellular compartments [28,31].In our simulation, the changes in skin keratin's secondary structure are considered a critical parameter for assessing nanotoxicity.We also analyzed the driving forces between the keratin and PG or GO and found that vdW, electrostatic interactions, and hydrogen bonding played a vital role in adsorption.The root mean square deviation (RMSD) and the root mean square fluctuation (RMSF) were also discussed to analyze the secondary changes of keratin.Our results revealed that the keratin could be absorbed onto PG in a few nanoseconds rebuilt and partially changed into different motifs due to the strong vdW interactions.It is demonstrated that the keratin could retain part of its original secondary structure upon approaching GO, implying that GO might have superior biocompatibility than PG.As a result, we herein proposed a mechanism of the nanotoxicity of 2D materials, such as PG and GO, by analyzing the conformation changes of skin keratin using all atom molecular dynamics (AAMD) simulations (figure 1).The current work contributes to the understanding of the nanotoxicity of cutting-edge 2D materials on human skin, thereby assisting experimental scientists in exploring the multifaceted applications of 2D materials in biomedicine.

Computational modeling and simulations
AAMD simulations were performed in this work to investigate the interaction mechanisms between skin keratin and PG or GO.To save computer resources, we chose the helical 2B segment of Keratin 1/10, an essential part of the whole skin keratin protein, as a representative model in our simulations.Among its several structure files in the protein data bank (PDB), 6UUI has the highest resolution at 2.07 Å from x-ray crystallography [32].6UUI contains 215 amino acid residues assembled as a left-handed coiled-coil heterodimer.Its secondary structure comprises amounts of α-helix segments, which is an essential sign to evaluate the configuration changes of keratin.Both sizes of PG and GO are 18 × 5 nm 2 .PG was built with the graphene builder plugin of VMD [33] (a visual molecular dynamics software).GO was modified from PG using GOPY [34] according to the molecular formula of C 10 O 1 (COOH) 0.5 (every 100 carbon atoms of PG with 10 epoxy groups, and 5 carboxyl groups), which was consistent with c = n O / n C around 20% (where the c stands for the degree of oxidation of PG, n O and n C stand for the number of oxygen-containing groups and carbons atoms, respectively) and had been used in other peer works [35,36].It is noted that the epoxy groups are arbitrarily located on both sides of the carbon atoms' basal plane.In contrast, the carboxyl groups are grafted randomly on the edge of the graphene basal plane.
The widely adopted force field of optimized potentials for liquid simulations-all atoms (OPLS-AA) [37] in the GROMACS software package (version 2021.4)[38] was adapted in our simulations.The key force field parameters are listed in table S1.All sp 2 carbon atoms in PG and GO were defined as uncharged Lennard-Jones (LJ) spheres with a default cross section (σ cc = 3.55 Å) and depth of the potential well (ε cc = 0.29 kJ mol −1 ) in OPLS-AA force field.Initially, the keratin was placed on the top center of the PG or GO surface parallel to their xz-plane.As shown in figure 2, the initial distance of the center of geometry (COG) between the keratin and graphene basal plane was set at 30 Å.The water molecules were simulated with the standard SPC/E model [39].LINCS algorithm [40] was used to constrain bond lengths.The Verlet was adapted to cutoff-scheme, and the cutoff value for van der Waals (vdW) and electrostatic interactions was set to 12 Å.The V-rescale temperature coupling [41] was used to control the temperature.The particle mesh ewald (PME) method [42] was set to calculate the long-range electrostatic interaction with a 12 Å real space cutoff.Periodic boundary conditions were also conducted in all three directions.Then, all systems were neutralized with counter ions (here were Na + and Cl − ), and the final concentration of NaCl was set at 0.15 M to simulate the physiological condition.All carbon atoms of the basal plane were frozen.The neutralized systems were pre-equilibrated after 5000 steps energy minimization in the NVT ensemble for 400 ps.Leap-frog algorithm [43] was used to integrate the equations of motion with a 2 fs time step.Annealing control was also used to heat the temperature from 0 to 310 K at 100 ps.Eventually, the production simulation continued at 310 K for 200 ns.

Data analysis of MD simulations
Herein, the entire adsorption process of the keratin onto the PG or GO surface was investigated using MD simulations.The non-bonded energy, including vdW and electrostatic interactions, the number of hydrogen bonds between the keratin and GO, the distance of COG, the contact area, and the changes of the secondary structures of keratin during the whole simulation were all considered to investigate such adsorption mechanisms at the atomistic level.The secondary structure of keratin was analyzed by the DSSP (dictionary second structure of protein) method [44], which was adopted in our previous works [26,45].The number of hydrogen bonds was calculated based on angle cutoff at 30°and bond length cutoff at 3.5 Å, which were widely accepted by other peer works [38,46,47].The evolution of the RMSD and the RMSF of the keratin during the whole simulation were calculated for its backbone using the relevant trajectory tool of GROMACS.The visual weak interactions were calculated with Multiwfn [48].The last 50 ns of the trajectory was extracted for calculating the free binding energy between the materials and keratin using gmx_MMPBSA and gmx_MMPBSA_ana [49,50].The vdW, electrostatic interactions, and the contact area between the keratin and PG or GO were all plotted as a function of time by the instinct modules implemented in the GROMACS software package.

Secondary structure changes of the keratin
Since the adsorption of keratin onto 2D materials is timedependent, three sets of 200 ns long trajectories are discussed in detail.Figure 2 shows some representative snapshots of the adsorption process of keratin onto the PG or GO surface during the first 100 ns simulations (the whole 200 ns simulation process is shown in figure S1).
The coiled-coil secondary structures of keratin are upheld by a diverse array of both covalent and non-bonded interactions among residues.Tzeng et al [51] demonstrated that the structure of the protein has a significant influence on regulating its biology activity.In addition, previous studies on nanomaterial toxicity have indicated that alterations in protein conformation, resulting from interactions with nanomaterials, constitute a primary factor contributing to toxicity [16,52].In this section, we calculated the conformational changes of the keratin to access the nanotoxicity of 2D materials.As shown in figures 2(a) and S1, the configuration of the keratin undergoes minor changes in the entire 200 ns simulation process when saturated in 0.15 M NaCl solution without 2D materials.For better comparison, we removed the center of mass rotational and translational velocity around the center of mass of the keratin.The overall configuration of the keratin is as intact as its initial state, implying that the bioactivity and biofunction of keratin could be preserved without 2D materials.However, the adsorption process in PG and GO systems is totally different.As shown in figure 2(b), the keratin is adsorbed onto the graphene basal plane in the first 10 ns, indicating PG possessed strong adsorption on the keratin.After that, these contacted residues never desorb from the graphene until the end of the simulation.It reflects that the vdW interactions between the keratin and graphene play a significant role in adsorption.Subsequently, the remaining domain of the keratin adsorbs onto graphene within 50 ns.At 60 ns, the keratin is completely adsorbed on the surface of PG.From these plots, we can see that the keratin's middle domain (rectangular area) is massively distorted, implying that the biological activities of the keratin might be affected by PG.To understand how the keratin approaches the graphene interface, the closest distance between the graphene basal plane and the keratin residues is calculated (figure 2(b)).In the simulation, we fixed the carbon basal of PG and GO to better analyze the adsorption process, whereas the carboxyl and epoxy groups could freely swing.At 10 ns, the closest distance between the keratin residue and carbon basal declines to 2.6 Å, which is smaller than the diameter of a water molecule (4 Å).This result reflects that the water molecules at the interface are pushed away because of the hydrophobic effect of the graphene basal surface.Du et al [53] proposed that the solid hydrophobic property of PG repulsed water molecules on its surface, creating a dehumidification layer promoting adsorption.In this PG system, vdW is the sole interaction between the keratin and PG because of the uncharged sp 2 carbon atoms.While in the GO system, there are additional driving forces, such as hydrogen bonding and electrostatic interactions to absorb the keratin.As illustrated in figure 2(c), the proper domain of the keratin is adsorbed onto the graphene basal plane at 20 ns.Subsequently, it acts as an anchor to grasp the GO firmly.In contrast, another domain fluctuates as a random swing until it is adsorbed on the surface of GO after 60 ns time of structural adjustment.Unlike PG, the sophisticated interactions between the keratin and GO, including vdW, electrostatic interactions, and hydrogen bonding, cause the tail of the keratin to fluctuate like a random swing for structural adjustment and interface rearrangement.In addition, we observe that the keratin prefers to be adsorbed on the edge of GO due to the charged carboxyl groups (from the top views in figure 2(c)).This finding supports prior research suggesting electrostatic interactions significantly attract charged molecules [54,55].The hydrophobic interaction also contributes to the adsorption, which is incorporated in the vdW interaction of classical nonpolarized MD force field [53].
Furthermore, we analyzed the time evolution of the RMSD and RMSF of the keratin backbone and secondary structure profiles.First, RMSD was calculated to investigate the backbone changes of the keratin during its adsorption process.A larger RMSD reveals the backbone change of keratin is more significant.In figure 3(a), we can observe that the RMSD of the keratin backbone fluctuates in a more extensive range when simulated without PG or GO in the entire simulation.In contrast, the RMSD of the keratin backbone varies at the beginning, but it shows slight variation and converges in PG and GO systems since 60 ns.Our results revealed that 2D materials could immobilize the keratin.These findings are likely to bring fresh insights into the future creation of advanced functional 2D biomaterials by modifying the oxygen-containing groups.Then, the RMSF of each keratin residue was also determined.The RMSF represents the average fluctuation of each residue relative to the reference position in the entire simulation process.Figure 3(b) shows the specific residue changes in amino acid movement and flexibility during the adsorption simulation.These results elucidate that PG and GO strongly influence restricting the structure changes of these residues and might affect the bioactivity and biofunctions of the keratin, particularly in the middle domain of each monomer residue 51 and 162.
Additionally, the detailed secondary structure variations of keratin were calculated using the DSSP (dictionary of protein secondary structure) method [44], which was implanted in the GROMACS software package.As shown in figures 3(c) and (d), the secondary structure of the keratin is retained constant in the absence of 2D materials.Still, a significant loss of α-helix structure in keratin is found with 2D materials.From figure S2, we also observe that the keratin's secondary structures change in the PG system.Our research findings indicate that the secondary structure of keratin remains relatively stable in a pure water system, maintaining a significant portion of its α-helix secondary structure (depicted in blue) and thereby preserving its stable biological activity.However, in the presence of PG systems, there are notable changes in the secondary structure of keratin.Most α-helix structures are converted into different forms, especially 3 10 -helix structures, reflecting that the secondary structure of the keratin is destroyed when absorbed onto PG by strong vdW interactions.Mukhopadhyay et al [17] reported that a portion of 3 10 -helix structures of HIV-1 integrase transformed into turn and random coil structures in phosphorene system, whereas the structure remained nearly unchanged in the pure water system.Du et al [53] found that graphene could cause the turn structures of SARS-CoV-2 spike receptor to become bend structures.Paul et al [56] illustrated that the secondary structure of hen egg white lysozyme remained largely unaffected upon adsorption on graphene, but h-BN significantly disrupted both the α-helix and β-sheet components by interfering with intraprotein hydrogen bonds.Furthermore, several studies have revealed the biological effects of MoS 2 on protein structure and function, for its advancing application in biomedicine.Gu et al [57] demonstrated that the monolayer MoS 2 triggered a conformational shift in polyalanine, resulting in a specific and subsequently stabilized decrease in the α-helix ratio.The interaction between polyalanine and MoS 2 nanosheet was primarily driven by van der Waals forces, exhibiting a rapid binding rate.MoS 2 nanosheet disrupted the intramolecular hydrogen bonds within polyalanine, leading to the breakdown of the initial α-helix structure.This ultimately induced significant conformational alterations, disrupting the secondary structure.All these results revealed that 2D materials might have different impacts on different proteins.From our simulations, we concluded that GO could preserve the secondary structure of keratin in a more extended continuity.It is worth noting that the secondary structure of the keratin changes significantly after 60 ns.It has been reported that the components ratio of the secondary structure is bound to influence the volume, shape, and biological functions of proteins [58], which are needed for the maintenance of the keratin structure integrity, and hinting that their variations are a critical parameter to evaluate the nanotoxicity of PG and GO materials on proteins.Our results, GO is more biocompatible than PG, are consistent with prior research [12,55,59,60].

Driving forces and energy changes during the adsorption process
This section analyzed the non-bonded interaction energies between keratin and 2D materials throughout adsorption.As shown in figure 4(a), the vdW interaction energies in the PG system achieve a stable negative state (around −1894.51kJ mol −1 ) at 125 ns, while the vdW interaction energies in the GO system reach a steady state (around −2090.54kJ mol −1 ) at 150 ns.Owing to the charged functionalized groups on the surface of GO, electrostatic interactions are significant in the GO system.It has been reported  The vdW interaction energies between the keratin and PG (red), the vdW (blue) and electrostatic (green) interaction energies between the keratin and GO, the total non-bonded energy between the keratin and GO (purple) as a function of time.As graphene has no charge and function groups, its total non-bonded energy is only vdW interaction energies.In contrast, the total non-bonded energy in the GO system is the sum of vdW and electrostatic interaction energies.that the intrinsic electron transfer from N to C atoms might cause C 2 N to be less hydrophobic, to some extent, than pure graphene [27].In our simulation, carboxyl and epoxy groups of GO could attract the electrons from C atoms, yielding a less hydrophobic carbon basal than pristine graphene.In addition, some keratin residues might be repulsed by the oppositely charged groups and their steric hindrance effects.As a result, the vdW interaction energies in the GO system are weaker than those in the PG system in the first 150 ns.However, as shown in figure 5(a), the distance of COG in the GO system is closer than that in the PG system, resulting in greater vdW interaction energies between keratin and GO.
Compared with electrostatic interactions, the vdW interactions play the most crucial role in the adsorption process between the keratin and PG or GO.We calculated the average non-bonded energies of the last 50 ns simulation due to their stability.The vdW (total non-bonded) interaction energy between the keratin and PG is −1864.09± 54.41 kJ mol −1 , weaker than the GO system (−2541.66± 91.14 kJ mol −1 ).Meanwhile, the electrostatic energy accounts for 18.3% of the total non-bonded interaction energy in the GO system.The more accurate binding free energy between materials and keratin illustrates that the GO system has a more substantial binding free energy than PG (−1754.64 ± 55.10 kJ mol −1 and −1549.71± 50.75 kJ mol −1 , respectively), identifying this specific GO sheet (C 10 O 1 (COOH) 0.5 ) shows more vital binding ability to the keratin.
To better understand how these weak interactions are distributed at the interface, we visualized the weak interactions between residue GLN162 and its surroundings (figure 6).From figure 3(b), we found that the RMSF of residue GLN162 differs significantly between the two systems.Our results revealed that the vdW interactions (small green icon) contribute to attract the residues, while the steric effects (small red cone) prevent further adsorption in the PG system.In the GO system, some hydrogen bonds (small blue icon) are formed between the -O (of epoxy group) atom and -H (of the side chains) atom.To further elucidate the adsorption process, we calculated the distance of COG and the contact area between the keratin and 2D materials.An important parameter is the distance of COG between keratin and the carbon atoms basal plane projected along the Y-axis (adsorption direction).At the initial state, the distance of COG was set to 30 Å. Two tendencies are fluctuated in the first 20 ns.Following that, the distance of COG between the keratin and PG steadily reduces until the keratin is completely absorbed onto the PG surface.While the GO system exhibits various phenomena, this might be caused by distinct interactions at the interface.The final 50 ns of data are magnified for better observation, and it is found that the distance of COG in the GO system is less than that in the PG system.These results reflect that the keratin could be absorbed onto the surface of GO closer in our simulations.
The contact area is defined as [61]: where SASA protein and SASA 2D material is the solvent accessible surface area (SASA) of the isolated protein and 2D material, respectively, and SASA complex is that of the 2D material-protein complex.It can be seen from figure 5(b), that the average contact area of keratin and PG is 5 Å 2 in 10-35 ns, which is consistent with process in figure 2. However, this contact area rapidly grows to 30 Å 2 by 60 ns.The contact area then stabilizes at roughly 35 Å 2 from 75 ns until the end of the simulation.The average contact area of keratin and GO exhibits a similar growing tendency, albeit slower than the PG system.For example, the initial contact area reflection occurs at 20 ns in the GO system, which is later than the PG system at 10 ns.After that, the contact area fluctuates in 20-80 ns during its approaching period.This might be caused by complicated interactions between the functionalized groups of GO and different charged residues of the side chains of the keratin at the interface.Subsequently, the contact area in the GO system is less than the PG system in the first 125 ns due to the electrostatic repulsion and steric hindrance effects of these functional groups on the GO surface.However, this contact area surpasses the value of the PG system at roughly 140 ns.Due to the aqueous condition of the protein, the hydrophilic residues tend to occur on the outer surface of the keratin, while the hydrophobic residues are inside.In our simulation, the hydrophobic residues proportion increases at the interface, indicating that when the keratin adsorbs onto PG and GO surfaces, the hydrophobic residues might turn to the outer surface, destroying the initial structures of the keratin.

Number of hydrogen bonds.
Hydrogen bonding has been identified as a significant interaction in proteins or protein-ligand complexes, playing a pivotal role in preserving structural stability, flexibility, folding, and binding affinity [62].Therefore, the number of hydrogen bonds was also calculated to assess the nanotoxicity of 2D materials.Therefore, we analyzed the number of intramolecular hydrogen bonds of the keratin as a function of time.As demonstrated in figure 7(a), the number of hydrogen bonds in the keratin in PG or GO systems decreases by 10 and 20, respectively.This is consistent with the secondary structure changes of the keratin in figure S2.The initial proportion of α-helix in the PG and GO system is 97.15%, while after 200 ns, it decreases to 89.19% and 82.73%, respectively.Furthermore, we also quantified the number of hydrogen bonds formed between keratin and GO surface during the entire simulation, as shown in figure 7(b).In the first 20 ns, no hydrogen bonds are found because the keratin is located too far away from these functional groups (the judging criteria of  hydrogen bonding is angle cutoff at 30°and bond length cutoff at 3.5 Å).At 20 ns, the number of hydrogen bonds increases substantially as the contact surface area increases rapidly.After being entirely adsorbed on the GO surface, the amount of hydrogen bonds between the keratin and the functional groups, including epoxy and carboxyl, reaches an equilibrium state at 7.

Conclusions
In summary, we aim to investigate the nanotoxicity of 2D materials on skin keratin at an atomistic scale via classical MD simulations.The adsorption process of the keratin on 2D materials, taking the prototypical 2D PG and GO, for example, was investigated.Our results show that the keratin could be spontaneously absorbed onto the PG or GO surface.However, the driving forces of the adsorption in PG and GO systems are different.Since graphene solely contains uncharged sp 2 carbon atoms, vdW interaction is the major attraction interaction between keratin and PG.In the GO system, a vast number of charged functional groups such as epoxy and carboxyl are randomly distributed on the surface and edge of GO, resulting in diverse interactions, such as vdW and electrostatic interactions, and hydrogen bonding, at the interface of the contact area.The 4.64% and 3.97% α-helix structures formed in keratin are converted to 3 10 -helix structures in PG and GO systems, indicating that the helix diameter is shrunk under the strong adsorption forces.While stronger interaction energy and free binding energy are observed, the conformational change in the GO system is less than that of the PG system, indicating GO has more favorable binding sites and biocompatibility than PG on the keratin structure.Consequently, our MD simulations demonstrate that PG and GO could absorb the keratin through various driving forces, including vdW, electrostatic interactions, and hydrogen bonding.According to our simulation results, GO has better adsorption characteristics and biocompatibility than PG on the keratin structure.Future studies should focus on investigating how functional groups affect the nanotoxicity of GO since it has been reported that the hydrophobicity and biocompatibility of GO are related to its physicochemical properties such as the oxygen content, dispersion and size [14,63].More diverse oxygen-containing groups will be simulated to study further and boost the application of GO and its derivatives in various biological applications.In addition, the issues pertaining to the degradation, excretion, and residual characteristics of 2D materials in vivo parallel challenges encountered in diverse biomedical applications of 2D materials.It is anticipated that an escalating number of scholars will direct their attention to this scientific domain in the future.

Figure 1 .
Figure 1.Schematic illustrating the mechanism of adsorption and structure changes of keratin on PG and GO surface.

Figure 2 .
Figure 2. (a) The movement trajectory of the keratin without 2D materials.(b), (c) The secondary structure of the keratin.The red circle indicates the remarkable protein adsorption area.The rectangular shapes stand for the middle domain of the keratin.The enlarged structure details show the contact area and some top views.Note that water molecules and ions are all deleted for clarity.

Figure 3 .
Figure 3.The (a) RMSD and (b) RMSF of the keratin structures in different systems.The time evolution of the secondary structure such as (c) α-Helix and (d) 3 10 -Helix of the keratin during the adsorption process in different systems.

Figure 4 .
Figure 4. (a), (b)The vdW interaction energies between the keratin and PG (red), the vdW (blue) and electrostatic (green) interaction energies between the keratin and GO, the total non-bonded energy between the keratin and GO (purple) as a function of time.As graphene has no charge and function groups, its total non-bonded energy is only vdW interaction energies.In contrast, the total non-bonded energy in the GO system is the sum of vdW and electrostatic interaction energies.

Figure 5 .
Figure 5.The time evolution of (a) the distance of COG and (b) the contact area between the keratin and PG or GO surface during the adsorption process in different systems.

Figure 6 .
Figure 6.The visual pictures of the weak interactions of the keratin in (a) the PG system and (b) the GO system.Deeper blue stands for stronger attraction, such as Hydrogen and Halogen bonds; green stands for vdW interactions; deeper red stands for stronger repulsion, such as Steric effects.

Figure 7 .
Figure 7.The time evolution of the number of hydrogen bonds (a) in keratin and (b) between keratin and GO surface during the adsorption process.