C60 fullerene against SARS-CoV-2 coronavirus: an in silico insight

Based on WHO reports the new SARS-CoV-2 coronavirus is currently widespread all over the world. So far > 162 million cases have been confirmed, including > 3 million deaths. Because of the pandemic still spreading across the globe the accomplishment of computational methods to find new potential mechanisms of virus inhibitions is necessary. According to the fact that C60 fullerene (a sphere-shaped molecule consisting of carbon) has shown inhibitory activity against various protein targets, here the analysis of the potential binding mechanism between SARS-CoV-2 proteins 3CLpro and RdRp with C60 fullerene was done; it has resulted in one and two possible binding mechanisms, respectively. In the case of 3CLpro, C60 fullerene interacts in the catalytic binding pocket. And for RdRp in the first model C60 fullerene blocks RNA synthesis pore and in the second one it prevents binding with Nsp8 co-factor (without this complex formation, RdRp can’t perform its initial functions). Then the molecular dynamics simulation confirmed the stability of created complexes. The obtained results might be a basis for other computational studies of 3CLPro and RdRp potential inhibition ways as well as the potential usage of C60 fullerene in the fight against COVID-19 disease.

www.nature.com/scientificreports/ symmetric nanostructure, such as C 60 fullerene, has extraordinary geometric affinity with icosahedral viruses 29 . The creation of nanobiomaterials on the same scale and with similar geometry is a fascinating possibility which can be used to foster interactions and build smart nanostructures for inhibition or inactivate virus replication. Although pristine C 60 fullerene has very low solubility in water, it can form a stable colloidal solution 30 , containing both single C 60 molecules (0.72 nm) and their spherical-like aggregates, the size of which is comparable to the average particle size of SARS-CoV-2-about 120 nm in diameter 31 . The mechanism of C 60 fullerene dispersal in aqueous solutions might be explained by a formation of a covalent bond between the hydroxyls and carbons in the C 60 molecule cage, as a result of ultrasound treatment that culminates in a consequent easy C 60 fullerene dissolution 32 . Due to hydrophobicity, C 60 fullerene easily penetrates the biological membrane by passive diffusion or endocytosis 33,34 . C 60 fullerene serves as an enzyme inhibitor, drug delivery vector, contrast agent for MRT and photodynamic therapy 33,[35][36][37] . Due to the presence of double electron-deficient chemical bonds in the structure, C 60 fullerene easily attaches free radicals, i.e. is a powerful antioxidant [38][39][40] , able to effectively exhibit the anti-inflammatory, antibacterial, antitumor, neuro-and radioprotective effects in the in vitro and in vivo systems 35,36,41,42 . Finally, importantly, C 60 fullerenes and their water-soluble derivatives exhibit remarkable antiviral activity 27,28,35,36 . Thus, in this work we first analyzed available 3CLpro and RdRp structures and their possibility to interact with C 60 fullerene employing computational methods. For this purpose, C 60 fullerene was docked into 3CLpro and RdRp according to obtained binding models (in the case of 3CLpro, it is blocking of a catalytic dyad, for RdRp-blocking of RNA synthesis pore and preventing of binding with Nsp8). Then the molecular dynamics (MD) simulation was performed on obtained "C 60 fullerene-3CLpro or RdRp" complexes. Subsequently, received MD trajectories were a subject of MMPBSA (molecular mechanics Poisson-Boltzmann surface area) as well as MMGBSA (molecular mechanics-generalized Born surface area) free-energy analysis. We believe that our computational results shed light on possible ways of 3CLpro and RdRp inhibition by nanobiomaterials, reveal the main flexibility properties of investigated targets and indicate favorable and unfavorable amino acid for interaction in selected binding pockets.

Calculation methods
Construction of "C 60 fullerene-3CLpro or RdRp" systems. According to available structure data 43,44 , the X-ray structures of monomeric 3CLpro protein (Protein Data Bank (PDB) ID 6M2N) 45 and RdRp protein (PDB ID 7BV2) 44 were retrieved from the RCSB PDB 42 . Firstly, all routine water molecules and native ligands were removed from the protein structure. Then prior to identifying binding pockets the protein structure treatment was done by the addition of missing hydrogen, correcting protonation states of amides, repairing side chains, and in the end, energy minimized. To run molecular docking simulation, the possible binding pockets for C 60 fullerene were defined by Caver software cavity computational algorithm 46 and based on literature analysis 14,21,45,[47][48][49][50] . As a result, three and one possible binding pockets were identified for RdRp and 3CLpro, respectively.
To each target molecular docking was carried out utilizing flexible C 60 fullerene molecule and rigid 3CLpro or RdRp molecule. The systematic docking algorithm was used (SDOCK+) 51 , implemented in the QXP package 52 (the method demonstrates all possible conformations of the studied structures with a minimum RMSD (rootmean-square deviation) value 51 ). The maximum number of SDOCK routine steps was set to 300, and the 10 best complexes based on the built-in QXP scoring function 52 were selected for analysis in the next stages of the investigation. The optimal structure of the studied "C 60 fullerene-3CLpro or RdRp" complexes was determined by the following basic criteria: (1) the area of the contact surfaces of the protein and ligand; (2) the distance between the 3CLpro or RdRp and C 60 fullerene; (3) the energy characteristics of the binding in the formed complex. As a result of molecular docking to each study one (best) "C 60 fullerene-3CLpro or RdRp" complex was selected.

MD simulation protocol.
To estimate stability and crucial interactions of obtained complexes after molecular docking, MD simulation was performed. The calculations were done using Gromacs 53 5.1.3 in force field 54 Charmm36. All exploring 3CLpro or RdRp were protonated according to the build-in function in Gromacs 5.1.3. The topology for C 60 fullerene was generated by SwissParam 55 . The complexes obtained after molecular docking were used for MD simulation. Each system was placed into the center of a periodic cubic box which was then filled with TIP3P water molecules. A minimum 0.9 nm distance was maintained between the nearest atom of the complex and the edge of the simulation box so that the complex can fully immerse with water and rotate freely. Then, to neutralize the system and mimic the cellular environment (pH = 7), Na + and Cl − ions were added to bring the ionic concentration to 150 mM. Here, the solvent molecules are replaced with monoatomic ions, randomly. Next, the obtained complex was energy minimized what also relieved any steric clashes. The system was relaxed by applying the steepest descent algorithm (the maximum number of steps was 50,000). Then the equilibration was computed in two stages: NVT was first equalized at 100 ps, with the second NPT equalization of 1 ns. After that, we launched MD simulation within 50 ns.
Note, that MD simulation was performed 3 times to each investigated "C 60 fullerene-3CLpro or RdRp" complex. All calculations were done at the temperature of 300 K and at constant atmospheric pressure.
Binding free energy calculations. The binding energies of each complex were calculated by applying the MMPBSA method of g_mmpbsa toll 56 according to the following equation 57,58 G binding = G complex − G protein + G ligand , www.nature.com/scientificreports/ where G complex : total free energy of "C 60 fullerene-3CLpro or RdRp" complex, G protein/ligand : total free energies of the isolated protein and ligand in a solvent, respectively. And then each parameter was estimated as follows where x-"C 60 fullerene-3CLpro or RdRp" complex or just separate protein or C 60 fullerene, G mm -average molecular mechanics potential energy in a vacuum (comprise energy of bonded (E bonded ) and nonbonded (E nonbonded = E vdW + E elect , where E vdW -van der Waals energy and E elect -electrostatic energy) interactions and calculated based on the molecular mechanics force-field parameters) 59-61 , TΔS-entropic contribution to the system, G solvatation -the solvatation free energy (energy required to transfer a solute from vacuum into the solvent). So, here the Poisson-Boltzmann (PB) equation has solved to estimate polar desolvation energy (G polar ) 62 . In the PB calculations the grid size was 0.5 Å. The solvent dielectric constant and solute dielectric constant value was 80 and 2, respectively. Nonpolar (G nonpolar ) contribution was obtained based on solvent accessible surface area (SASA, Å 2 ) where γ-surface tension of the solvent and b-fitting parameter 63 .
Next MMGBSA technique was used to obtain more precise picture of "C 60 fullerene-3CLpro or RdRp" interaction 64,65 . The essential idea of MMGBSA is the estimation of binding energy between each residue of 3CLpro or RdRp and C 60 fullerene. The algorithm of the MMGBSA method is the same as in MMPBSA and presented above. The whole energy parameters were calculated by utilizing all snapshots from the MD simulation trajectory of length 50 ns.

Result and discussion
Binding pocket determination. According to the surface analysis of 6M2N and available literature data 21,22,[66][67][68] , the only binding pocket capable of interacting with C 60 fullerene was identified (Fig. 1).
It is a catalytic binding pocket. This pocket occupies volume 892 Å 3 and contains amino acids (Met 165, Met 49, Tyr 54, Cys 44, His 41, and Cys 145), which are able to create stacking interactions with C 60 fullerene. Furthermore, taking into account that this binding pocket directly mediate the mutation of different Nsp 21,69,70 , which are crucial for the virus life cycle, it make 3CLpro a perspective target for drug development against SARS-CoV-2 17,21,22 .
As a result of surface and literature investigation of RdRp, three binding pockets were found. It is well known that RdRp plays a key role in SARS-CoV-2 RNA synthesis 14 . Based on our result two binding models between C 60 fullerene and RdRp were indicated: model 1-C 60 fullerene binding in RNA synthesis channel and model 2simultaneous binding of C 60 fullerene with pocket 2 and 3 ( Fig. 2). So, by blocking the RNA synthesis channel the RNA synthesis procedure will be impossible 19,21,71 . That is why the first detected binding pocket locates in the RNA synthesis channel (Fig. 2). On the other hand, without the assistance of Nsp7 and Nsp8 as co-factors RdRp is not able to carry out its initial functions 14,72,73 . Based on these two other pockets were found in the RdRp-Nsp8 binding interface (Fig. 2). All found binding pockets comprise at least two amino acids, which can create any stacking interactions with C 60 fullerene. For example, pocket 1 contains Arg 570, Lys 578 and Tyr 690.   www.nature.com/scientificreports/ RdRp as with 3CLpro, C 60 fullerene filled in selected binding pockets of RdRp and tightly clamped there by different stacking and steric interactions (Fig. 4A,C,E). The inhibition of pocket 1 could cause blocking of RNA synthesis channel (Figs. 2, 4A). Here, C 60 fullerene fits perfectly to the binding pocket and makes π-cation interactions with Arg 570 and Lys 578, T-stacking with Tyr 690, and steric interactions with Asn 497 and Leu 577. Conversely, the C 60 fullerene interaction with pocket 2 or 3 depicts the model which prevents complex formation between RdRp and Nsp8 (Figs. 2, 4C,E). As a result, RdRp is not able to carry out its initial functions. So, in pocket 2 the bottom part of C 60 fullerene stacks between Trp 510, Phe 369 and Leu 372, Leu 515 by stacking and steric interactions, respectively. Additionally, Tyr 516 and Phe 507 are located at the bottom of the binding pocket (Fig. 4C). These two amino acids possibly are capable of holding C 60 fullerene in the current position by stacking interaction. In spite of mentioned above the stability of the obtained complex is questionable because the binding pocket itself isn't deep. Because of such a flat surface geometry of pocket 2, C 60 fullerene could be forced out of this binding pocket. Despite the fact that in pocket 3 almost not presents any aromatic amino acid which is able to create stacking interactions with C 60 fullerene, we think that pocket 3 is promising because of its depth. Here, C 60 fullerene creates steric interactions with Ala 384, Val 331, Val 399, Thr 325, and Leu 271. Also, the binding pocket contains Phe 397 and Tyr 274, which are spatially close to docked C 60 fullerene. So, there is a possibility of stacking interaction with those amino acids.
The molecular docking results suggested that in all selected binding pockets C 60 fullerene is able to create a stable complex with 3Clpro and RdRp targets. The binding with 3CLpro is characterized by one possible binding model. For the RdRp, model 1 (interaction in pocket 1) and model 2 (simultaneous interaction with pockets 2 and 3) were investigated.

MD analysis.
To obtain more accurate results, 50 ns MD simulations were carried out (Figs. 3, 4, 5, 6). The MD results showed that each investigated "C 60 fullerene-3CLpro or RdRp" are stable. The RMSD movement during MD simulation of each complex is in a range 2-3 Å (Fig. 5). Furthermore, in some examples C 60 fullerene is able to form new and more profitable interactions.
3CLpro the simulated complex is shown in Fig. 3B. For 3CLpro, it was observed that C 60 fullerene shifts by 3.2 Å and caused Asn 142 displacement for 4.8 Å toward C 60 fullerene. In contrast, C 60 fullerene forced out Gln 189 for 1.8 Å from the catalytic binding pocket. The binding with amino acids like Met 156, Phe 181, His 164, and others are without any fundamental changes. The more intriguing is that both His 41 and Cys 145 (catalytic dyad) are forced out from their initial position by C 60 fullerene, either. As a result, the integrity of the catalytic dyad is violated and without any doubt, it has a negative impact on 3CLpro functionality. The "C 60 fullerene-3CLpro" complex was then subjected to RMSF (root-mean-squared-fluctuation) investigation. No, any extreme fluctuations compare to the free form of the 3CLpro molecule were detected (Fig. 6A). For example, the flexibility of His 41 in free and bond to C 60 fullerene form are 0.8 and 1.4 Å, respectively. Moreover, in the case of the "C 60 fullerene-3CLpro" complex, the reduction of fluctuation values was determined.
RdRp according to the MD simulation result in a case of pockets 1 and 2 the fluctuations of C 60 fullerene inside both binding pockets (4.1 Å and 3.5 Å, respectively) without significant changes in obtained complexes (Fig. 4B,D) were observed. And, in pocket 3 opposite picture was detected. Here, C 60 fullerene immersed inside the binding pocket by 4.0 Å. Anyway, the key interactions between C 60 fullerene and RdRp in whole models remain (Fig. 4B,D,F). So, almost no changes were indicated during MD simulation for pocket 1. Here, it is possible to say that amino acids Ile 590, Tyr 690, Leu 577, and Gln 574 are rigid and just slightly change their initial position during MD simulation (< 0.75 Å), those amino acids interact with C 60 fullerene by strong steric interactions. Contrarily, Lys 578 and Arg 570 are more flexible and characterized by some displacement (1.5 Å and 1.6 Å, respectively). Interestingly, that the π-cation interaction with Lys 578 isn't stable, but with Arg 570 is stable. Arg 570 is tightly locked to C 60 fullerene. Such difference in the interactions between C 60 fullerene and Lys 578/Arg 570 could be related with the location of those amino acids in the binding pocket. The Lys 578 locates on the edge of the pocket and can freely move (especially, because of Lys amino acid long linker) in contrast Arg 570 locates inside the binding pocket and has less space for movement. In pocket 2 (Fig. 4E) C 60 fullerene stuck in the insight the pocket, creating steric interaction with Leu 372, Leu 515, and stacking with Phe 369 and Trp 510. Notably, during the whole MD simulation trajectory, the displacement of key binding amino acids (e.g. Phe 369, Trp 510, Tyr 516, Leu 315 and Leu 367) is minimal, about 1 Å. Also, we think that in this position C 60 fullerene is held by stacking interaction with Tyr 516, which locates in the bottom of the binding pocket 2. Thus, the distance between C 60 fullerene and Tyr 516 is about 3 Å during all MD simulation. Finally, in pocket 3 inverted picture compare to the above has observed (Fig. 4F). C 60 fullerene shifted inside the binding pocket and tightly stacked among the surrounding hydrophobic amino acids by interacting with RdRp via stacking and steric interactions with Phe 397 and Val 399, Ala 384, Val 331, Thr 325, respectively. Those interactions become possible due to the shift of Phe 397 (4.9 Å), Val 399 (3.5 Å), and Ala 384 (3.2 Å). Other amino acids Val 331 and Thr 235 located in pocket 3 almost not shifted (about 0.5 Å). As for the previous target, here RMSF analysis has been done as well (Fig. 6B). So, the amino acids 324-342 and 364-410 in the case of C 60 fullerene binding in pockets 2 and 3 are more flexible compare to the model, where C 60 fullerene bond to pocket 1 and to unbound RdRp. It is not surprising and could be simply explained by shifting and immersion of C 60 fullerene in the binding pockets 2 and 3. Of course, during the shifting/immersion C 60 fullerene in some level pushes amino acids, and it resulted in their mobility. However, in a model with C 60 fullerene bond pocket 1 mirrored situation is observed. Here, the flexibility of amino acids 840-862 is greater compare to their flexibility in the case of unbound RdRp or bond to pocket 2 and 3. Such a result was obtained due to the fact that in the case of C 60 fullerene binding in pocket 1 RNA molecule isn't present in RNA synthesis pore (there isn't enough space for RNA molecule). www.nature.com/scientificreports/ www.nature.com/scientificreports/ MMPBSA approach. The potential binding affinities within "C 60 fullerene-3CLpro or RdRp" complexes were estimated by the MMPBSA analysis. Table 1 suggests that the G binding in "C 60 fullerene-3CLpro or RdRp pocket 3" complexes are more favorable compare to "C 60 fullerene-RdRp pocket 1/2". Next, for a better understanding of which energy type has a bigger contribution to complex formation and stability, each separate energy component was analyzed. From Table 1 it can be found that E vdW has the largest contribution to the binding energy of the C 60 fullerene with 3CLpro and RdRp. As previously, the E vdW values of "C 60 fullerene-3CLpro or Figure 5. RMSD trajectories of 3CLpro (A) and RdRp (B) complexes with C 60 fullerene: 3CLpro free protein molecule-green, "C 60 fullerene-3CLpro" complex-red, and RdRp free protein molecule-green, "C 60 fullerene-RdRp" pocket 1 (binding model 1)-red, "C 60 fullerene-RdRp" pocket 2/3 (binding model 2)-orange. Figure 6. RMSF of the 3CLpro (A) and RdRp (B) in complex with C 60 fullerene: 3CLpro free protein molecule-green, "C 60 fullerene-3CLpro" complex-red, and RdRp free protein molecule-green, "C 60 fullerene-RdRp" pocket 1 (binding model 1)-red, "C 60 fullerene-RdRp" pocket 2/3 (binding model 2)-orange. www.nature.com/scientificreports/ RdRp pocket 3" complexes are far better compare to "C 60 fullerene-RdRp pocket 1/2". The other calculated energies (E elect , G polar and G nonpolar ) are positive or close to null. That is why their impact on complex formation and especially stability is unfavorable. However, the E elect and G nonpolar contributions are lightly favorably than G polar . Anyway, this slightly bigger effect of E elect and G nonpolar is not much more noticeable in comparison with G polar .

MMGBSA approach.
To get a more precise binding energy characterization of obtained complexes, a perresidue free-energy decomposition study was performed. According to the results presented in Fig. 7 the most favorable contribution in all models is done by amino acids, which are able to create stacking interaction with C 60 fullerene. However, some exceptions are detected (e.g. "C 60 fullerene-RdRp pocket 1"). 3CLpro Fig. 7 A emphasizes the role of the catalytic dyad in the binding process. The His 41 and Cys 145 have a positive effect on binding with C 60 fullerene and the binding energy is -4.63 and -3.88 kJ/mol, respectively. As was mentioned above, due to this, the C 60 fullerene can immerse into the 3CLpro catalytic binding pocket. Nevertheless, the binding contribution of the catalytic dyad isn't most favorable. The amino acids Met 49 (binding energy with C 60 fullerene is − 5.54 kJ/mol), Met 165 (− 8.18 kJ/mol), Leu 50 (-3.99 kJ/mol), Gly 143 (− 2.58 kJ/mol), Asp 187 (− 2.83 kJ/mol) and Gln 189 (− 5.25 kJ/mol) clamp C 60 fullerene in the binding pocket by stacking and steric interactions (Figs. 3B, 7A). Surprisingly, contrary to the above Arg 40 (6.94 kJ/mol) causes the destabilization of the "C 60 fullerene-3CLpro" complex formation. This can be explained by the fact C 60 fullerene prevents the interaction of Arg 40 with water molecules which usually in abundance in any binding pocket.
RdRp the comparison of key amino acids in pocket 1 has shown that steric interactions with C 60 fullerene here are the most favorable (Fig. 7B). Among the whole residue, the interaction between C 60  And finally in this case the C 60 fullerene interaction with Phe 327 is unfavorable. That can be related to the fact that only the peptide backbone of Phe 327 interacts with C 60 fullerene. Furthermore, this part of the backbone by C 60 fullerene isolated from a solvent, what possible is able to cause such tension between those molecular parts/structures. Summarizing, the in silico approach allows simulating the behavior of C 60 fullerene in the binding sites of SARS-CoV-2 coronavirus and thus to predict the therapeutic effect of this unique molecule. A complementary interaction of C 60 fullerene with proteins is a basis of its biomedical effects [74][75][76][77] . So, molecular docking and MD simulation were carried out using Toll-like receptors (TLRs play extremely critical roles in maintaining the immune-homeostasis of the human body) including TLR4 75 . The binding of C 60 fullerene with TLR4 is characterized by complete filling of the hydrophobic pocket of MD-2 domain and the formation of a significant number of stacking interactions (e.g. with Phe 119, Phe 76 and Phe 104). This change of binding site is associated with significant mobility of interacting components: RMSD value for protein is 4.6 Å, and for C 60 fullerene-5.3 Å. The obtained "C 60 fullerene-TLR4" complex was characterized by a high energy: − 50 kJ/mol. In a recent study 78 , spike protein of SARS-CoV-2 was found to interact with the extracellular domain of the cell surface TLRs. Intriguingly, the highest binding affinity and strength were evident in the "spike protein-TLR4" complex. Thus, the usage of C 60 fullerene to inhibit TLR4 as well as 3CLpro and RdRp activation may be an effective strategy to treat COVID-19. www.nature.com/scientificreports/ Figure 7. Per-residue binding energy decomposition of investigated complexes: "C 60 fullerene-3CLpro" (A); "C 60 fullerene-RdRp" pocket 1 (B); "C 60 fullerene-RdRp" pocket 2 (C); "C 60 fullerene-RdRp" pocket 3 (D). www.nature.com/scientificreports/

Conclusion
The computer simulations (docking and molecular dynamics) we presented here suggest that C 60 fullerene is able to block 3CLpro (blocking of catalytic dyad) and RdRp (blocking model 1 and 2) protein targets of SARS-CoV-2 coronavirus with different mechanisms and suppress its functional activity. The simulations revealed that in all investigated complexes C 60 fullerene filled in the binding pocket and stuck there by the stacking and steric interactions. Critically that for 3CLpro C 60 fullerene violated catalytic dyad integrity. All the other changes during simulations there weren't significant. In the case of RdRp pockets 1 and 2, C 60 fullerene just fluctuates inside the binding pockets without fundamental changes of previously obtained complexes. The reverse picture has been observed in RdRp pocket 3, here C 60 fullerene immerses inside the binding pocket and stuck there.
The MMPBSA study has shown that in all cases G binding is more favorable in the case of "C 60 fullerene-3CLpro or RdRp pocket 3" complexes compare to others. And as the main component, E vdW has the biggest contribution to all complexes. Furthermore, the contribution of E elect , G polar , and G nonpolar are questionable and most of all are unfavorable because of that energies proximity to null or in some cases that energies are far bigger than null. And finally based on MMGBSA investigation favorable and unfavorable amino acids for complex formation with C 60 fullerene were detected. The results of the study can provide understanding of 3CLpro and RdRp binding with other nanobiomaterials. Moreover, since the pristine C 60 fullerenes can form a high stable aqueous colloidal solution and exhibit anticoronavirus activity, this expands their use for prophylactic and therapeutic purposes, which requires further in vitro and in vivo testing.