Tracking the interaction between single-wall carbon nanotube and SARS-Cov-2 spike glycoprotein: A molecular dynamics simulations study.

COVID-19, a newly discovered type of coronavirus, is the cause of the pandemic infection that was first reported in Wuhan, China, in December 2019. One of the most critical problems in this regard is to identify innovative drugs that may reduce or manage this global health concern. Nanoparticles have shown a pivotal role in drug delivery systems in recent decades. The surface of nanoparticles could be covered by a layer composed of different biomolecules (e.g., proteins and macromolecules) following the incubation with a biological fluid. This protein-rich layer is called "Protein Corona." In this study, an all-atom molecular dynamics simulation was used for investigating the monomeric B domain of the spike glycoprotein due to its role in the accessibility of the spike glycoprotein to single-wall carbon nanotubes (SWCNTs). The interaction energy values between the carbon nanotube and B domain of the viral spike glycoprotein were evaluated. The obtained results, based on Lennard-Jones potentials, demonstrated that SWCNTs had an affinity to the B domain of the S1 subunit in the spike glycoprotein. The adsorption of SWCNTs on the B domain surface led to a significant change in solvent-accessible surface, internal hydrogen bonds, and finally in the tertiary structure, which could provide a reasonable method to impede the interaction between the angiotensin-converting enzyme II and SARS-CoV-2 spike glycoprotein. A decrease in the mean square displacement of the B domain was shown after the adsorption of SWCNTs as a result of increasing the hydrophobic-hydrophilic properties of the B domain. The arrangement of SWCNTs on the B domain surface and their interaction using the 2-acetamido-2-deoxy-β-d-glucopyranose group (988, 991, and 992) demonstrated that a change in the affinity of the S1 subunit could be used as a barrier to viral replication. The analysis of the SWCNT-B domain complex indicated that the presence of SWCNTs is able to cause alterations in the S1 subunit of the spike protein, and these nanotubes could be employed for further in-vitro and in-vivo antiviral studies. Also, SWCNTs are able to be utilized in drug delivery systems.

Introduction

The cylindrical structure of carbon nanotubes (CNTs) has a nanoscale diameter. CNTs were discovered accidently by Oberlin and Sumio Iijima in 1976 [1]. Biomedical materials are widely used in medicine, such as tumor therapy, and among these agents, nanomaterials have attracted much attention in clinical applications [[2], [3], [4]]. Nanomaterials have a high surface-to-volume ratio in comparison with the bulk materials and are chemically modifiable, making them an ideal medium for the attachment of biological molecules [5]. CNTs are one of the members of carbon-based nanomaterials [6,7]. These materials are broadly applied in nanomedicine because of their potential to interact directly with biomolecules, including peptides, lipids, nucleic acids, and proteins for site-specific drug delivery [8,9]. The dynamic interaction between CNTs and the biological system built to query is an important yet frequently underestimated problem for these nanoscale instruments [10]. A group of nanoparticles in biological environments are masked by a layer of proteins called “protein corona,” altering nanoparticle bio-distribution and toxicity [11]. The formation of protein corona is considered a pitfall in nanomedicine and clinical translation [12]. Protein binding to nanoparticles may have a detrimental effect on the affinity, structure, and function of proteins, masking and redefining the identification of nanoparticles [13,14]. In this study, the formation of the protein corona complex was performed to analyze the interaction between the SARS-CoV-2 spike receptor-binding domain and single-wall carbon nanotubes (SWCNTs) as well as structural changes and uptake of the protein. Three zoonotic coronaviruses have breached the species boundary to cause serious pneumonia in humans in the last two decades as follows:

SARS-CoV (severe acute respiratory syndrome coronavirus) was linked to an epidemic in 2002–2003 and a couple more cases in 2004 [15].

MERS-CoV (Middle-East respiratory syndrome coronavirus) is a coronavirus currently circulating in the Arabian Peninsula [16].

The COVID-19 pandemic is caused by SARS-CoV-2 coronavirus-19 [17].

The COVID-19 epidemic, which has caused an acute health crisis, has infected millions of people worldwide since December 2019 [18]. SARS-CoV-2 consists of four structural proteins named nucleocapsid (N), matrix (M), envelope (E), and spike (S) to protect its single-positive strand RNA (ribonucleic acid) genome [19]. According to recent studies, SARS-CoV-2 uses the angiotensin-converting enzyme II (ACE2) cellular receptor for cellular binding and internalization [20,21]. According to recent research, SARS-CoV-2 employs the ACE2 cellular receptor for cellular binding. This interaction between the spike protein and ACE2 occurs at the first point of viral entry to initiate viral fusion to the host cell [22]. The transmembrane spike glycoprotein, which forms homotrimers protruding from the viral surface, is responsible for Coronavirus entry into host cells [23]. The S protein is a viral fusion class I protein with a trimeric and crown-like structure that protrudes from the virus envelope and targets diverse host cell receptors in various organisms [24]. S1 and S2 are two subunits of the S protein that, during the first stage of viral infection, are cleaved. The S1 subunit is in charge of the S protein's interaction with the host receptor (ACE2), while the S2 subunit facilitates the viral fusion [25]. In the S protein, the S1 subunit plays a crucial role in the formation of two forms of spike proteins: inaccessible (closed) and accessible (open) proteins [26]. An asymmetric reconstruction of the trimmer with a single subunit B domain establishes an accessible version of the S protein. The S protein trimers in the accessible form are found in severe infectious diseases caused by coronaviruses; on the other hand, in the common cold, the inaccessible conformation is mostly found [27]. Regarding the significant role of the B domain in the accessible type of SARS CoV-2 spike glycoprotein, the complex of protein corona consisting of single-wall carbon nanotubes and the B domain of the S subunit was established, and their interactions were analyzed using the computational methods (Fig. 1 ).

Fig. 1

A schematic representation of interactions of spike glycoprotein with the ACE2 in the present and absent the single wall nanotube. (a) a trimeric structure of spike glycoprotein accessible type in interaction with ACE2, (b) B monomer of the spike glycoprotein in interaction with SWCNT.

Models and methodology

Design of system model

The accessible type of the S protein was obtained from the Protein Data Bank (PDB code: 6VYB) [28]. The Modeler software version 9.25 was used to model the B domain to fix the missing residues [29]. Regarding the fact that the spike protein is a glycoprotein, the number of NAGs for each domain was calculated in the crystalized structure. Since the function of NAGs in protein behavior is not entirely apparent, NAGs were considered for further analyses. The CHARMM-GUI software was utilized to add 19 NAG (2-acetamido-2-deoxy-beta-d-glucopyranose) groups to the B domain structure [30]. The carbon atoms on SWCNTs were used as uncharged in metallic conductivity, with a length and diameter of 10 and 4 nm, respectively (Table 1 ) [31]. As shown in the initial configuration of the B domain- SWCNT complex, the average distance of SWCNTs and the protein was set to be 2 nm (Fig. 2 ). In order to determine how the wall of SWCNTs affects the B domain, two system models were prepared (Table 1).

Table 1

Systems studied and the corresponding details in this study.

No. of simulation	Model systems	Length (ns)	No. of water	Box size	Box type
MD1	B domain	100	266463	20.61563*20.61563*20.61563	Cubic
MD2	B domain -SWCNT	100	361213	22.17581 *22.17581*22.17581	Cubic

Fig. 2

(a) B chain of the spike glycoprotein before simulation in the absence of the SWCNT, (b) Merged structure of SWCNT and the B chain of the spike glycoprotein before simulation with average distance 2 nm.

Simulation methods for B domain spike glycoprotein-SWCNT interactions

All molecular dynamics simulations were performed by the GROMACS software [32]. The CHARMM36 force-field software was employed to parameterize the protein [33]. The interatomic potential for carbon atoms in SWCNTs was calculated using the harmonic cosine of the bend angle, double twist potential, vdW interactions, and LJ parameters [34]. Also, 1 bar pressure and 310 K temperature were applied to all simulation systems. A Nose-Hoover thermostat with a coupling time of 0.5ps was used to control the temperature [35]. The pressure control was carried out by coupling the simulation cell to a Parrinello-Rahman barostat, with a coupling time constant of 5 ns. Isotropic pressure coupling was performed in all directions [36]. Appropriate numbers of sodium and chloride ions were added to each simulation box to neutralize the systems. In all simulation systems, periodic boundary conditions were used along all simulations box axes, and the transferable intermolecular potential 3 points (TIP3P) water model [37] was applied to solvate the systems. All-atom bond lengths were linked using the LINCS algorithm [38]. The Lennard-Jones potential was used to calculate the van der Waals (vdW) interaction with a cut-off radius of 1.2 nm, as recommended [39]. At first, while the positions of the protein and SWCNTs were restrained, the equilibration was conducted in the NVT ensemble for 5 ns, and then the equilibration was continued in the NPT ensemble for 10 ns in both MD1 and MD2 system. At the end of NVT and NPT ensembles, the restrained positions were removed. After equilibration runs, both simulation systems were run for 100 ns The characteristics of each simulation system are depicted in Table 1(Fig. 2).

Results

Protein stability

The Root-Mean-Square Divergence (RMSD) and the Mean Square Displacement (MSD) parameters were used to investigate the B domain stability (Fig. 3 ). The variations in the average of backbone RMSD of the B domain during conformational equilibration are shown on the outer surface of CNTs and in the water (Fig. 3a). The B domain conformation effectively reached equilibrium within 80 ns, as displayed in Fig. 3a, and the simulation time scale of 100 ns was enough to create stable interactions between the protein and CNTs. The fluctuations of the B domain in the time scale of 25–60ns in the absence of SWCNTs are demonstrated by the black graph (Fig. 3a). However, the presence of SWCNTs at the vicinity of the B domain renders substantial stability to the B domain structure (red graph). After 80 ns, a relative equilibrium was established in both systems, denoting that the B domain fluctuation is still evident in the B domain in the aqueous medium after 80 ns compared with the B domain-SWCNT complex. The movement of the B domain in both systems is shown in Fig. 3b. The B domain movement in both systems was demonstrated in the first 80 ns, and the increase of MSD was apparent. It is noticeable that in the first 20 ns, the MSD value of MD2 was slightly higher than MD1, but after 20 ns, the B domain in the water shows a further increase.

Fig. 3

Evolution of (a) the backbone RMSD that in the presence of SWCNT the stability of the protein is apparent (red line), (b) The MSD of protein in both MD1 and MD2 systems that in the absent SWCNT increasing of the protein movement is shown as it is confirmed in the RMSD (black line).

Interaction between SWCNTs and proteins

The Lennard-Jones (LJ) potential and the electrostatic interactions between the B domain and SWCNTs were measured to evaluate the penetration mechanism of SWCNTs. Fig. 4 b illustrates the LJ potential between the B domain and SWCNTs. The attraction force between the B domain and SWCNTs is verified by the negative values of LJ potential after 38 ns. This high potential of interaction could interpret the B domain-SWCNT interaction. The electrostatic potential between the B domain and SWCNTs is not notable. After 38ns of simulation, the results showed that the average interaction energy between SWCNTs and the B domain was −395.56 ± 22kJ/mol-1 showing that the B domain has significant Van-der-Waals interactions with SWCNTs. This indicates that B domain-SWCNT interaction is formed after 38ns due to decreased distance between the B domain and SWCNTs (Fig. 4c and d). Upon the approach of SWCNTs to the B domain, the interaction was increased and remained stable until the end of the simulation. The distance is decreased less than 2 nm in the minimum distance within 18–20ns; however, since the complex was not yet stable, fluctuations within 38ns can be observed (Fig. 4c). The adsorption of the B domain on the surface of SWCNTs can also be measured using the “number of contacts” parameter. The following expression was used to calculate the number of contacts created between the B domain and SWCNTs [40]:where NSWCNT and N B-domain are the total numbers of atoms in SWCNTs and the B domain, respectively. The rj is the jth distance of atoms in the B domain from the ith atom in SWCNTs. The number of contacts between the B domain and the surface of SWCNTs is calculated by a simulation time function as shown in Fig. 4d. This graph reveals that the number of contacts in the B domain-SWCNT complex is significantly increased after 38 ns and then remained constant until the end of the simulation. The results showed a strong correlation between the interaction energy and the number of contacts between the B domain and SWCNTs (Fig. S2).

Fig. 4

(a) Comparison of interaction between SWCNT and B domain of the accessible spike glycoprotein, (b) Lennard-Jones (LJ) potential energy profiles between the B domain of the accessible spike glycoprotein and the SWCNT, (c) The minimum distance of B domain of the accessible spike glycoprotein and SWCNT, and (d) Number of close contacts of the B domain of the accessible spike glycoprotein atoms within 0.5 nm from the SWCNT surface.

Protein structural changes

After the adsorption of SWCNTs on the B domain, it is essential to analyze how the secondary structure of the peptide is influenced. The secondary structure contents in MD1 and MD2 are plotted against the simulation time (Fig. 5 a). In the presence of SWCNTs, the coil and β-bridge proportions of the peptide were reduced compared with the absence of SWCNTs (MD1). In contrast, the bend and turn structures were slightly increased in the presence of SWCNTs (MD2). The adsorption of SWCNT did not induce significant structural changes in the B domain when the difference in the secondary structure was assessed. However, in the presence of SWCNTs, the hydrophilic-hydrophobicity properties of the peptide exhibit a slight increase (Fig. 5b). Due to the structural alterations in the B domain, hydrogen bonds were also investigated (Fig. 5c). The internal hydrogen bonds and their interactions in the B domain with water molecules were analyzed in MD1 and MD2. There were no hydrogen bonds observed between the B domain and carbon atoms of SWCNTs because SWCNTs cannot form hydrogen bonds. Our observations revealed that the average number of hydrogen bonds between the B domain and water molecules in MD2 was higher than that of the average number of hydrogen bonds between the B domain and water molecules in MD1 (2228.06 ± 35 vs. 2182 ± 44, respectively). Subsequently, the average number of hydrogen bonds in MD2 was decreased compared with MD1 (674.11 ± 17 vs. 689.43 ± 21, respectively). The Lj interactions between the B domain and SWCNTs corroborated this fact in MD2.

Fig. 5

(a) Comparison of secondary structure contents in the presence and absence of SWCNT. (b) Time evolution of the B domain's hydrophilic-hydrophobicity level. (c) The number of hydrogen bonds between protein and B domain-water in the different systems, that the average of hydrogen bond for each group mentioned in graphs. (d) solvent-accessible surface of the B domain-no SWCNT and protein-SWCNT.

Moreover, the water molecules form higher numbers of hydrogen bonds with the B domain in MD2 in comparison with MD1. According to our findings, the presence of SWCNTs plays a significant role in the internal hydrogen bonds in the B domain and the interaction of the B domain with water molecules. Furthermore, we compared the solvent-accessible surface of the peptide in two positions: the protein in the absence of SWCNTs and the protein in the presence of SWCNTs. The configuration changes in the B domain were also evaluated when SWCNT is present (Fig. 5d). In MD1, the black graph shows the B domain in the absence of SWCNTs, in which a 20-nm decrease was found in the solvent-accessible surface of the peptide after 38 ns with significant fluctuations. In contrast, in the presence of SWCNTs, the red graph clearly shows higher stable fluctuations in the solvent-accessible surface of the B domain after 38 ns, according to the previous sections. Such stability is due to the interaction between the B domain and SWCNTs.

Tracking the SWCNT on the B domain

According to the results, the location of SWCNTs on the B domain was shown to explain the spatial changes in the B domain. As previously mentioned, each domain of the spike glycoprotein trimmer is made of two parts, namely the S1 and S2 subunits (Fig. 1). The presence of SWCNTs at the vicinity of the B domain induces the Lj interactions with the S1 subunit, according to a schematic comparison between the MD1 and MD2 structures (Fig. 6, Fig. 7, Fig. 8 a-c). The S1 subunit is a crucial agent in the formation of interactions between the spike glycoprotein and ACE2. The interaction of SWCNTs with the B domain causes a conformational change in the S1 subunit and alters the configuration of the active site in the Spike-ACE2 complex. In contrast, in MD1, some shifts are found in the S2 subunit in the absence of the SWCNTs during simulation. Such shifting does not play a direct role in the binding of the B domain to the ACE2 enzyme (Fig. 7 a-c). The amino acid positions of the binding site in the spike glycoprotein-ACE2 complex are represented as spheres (Fig. 6, Fig. 7, Fig. 8). In recent studies, the binding site contains Thr 453, Leu 455, Phe 456, Ala 475, Cys 480, Gly 482, Val 483, Glu 484, Gly 485, Phe 486, Asn 487, Cys 488, Thr 489, Phe 490, Gln 493, Ser 494, Thr 495, Gly 496, Gln 498, Thr 500, Asn 501, Gly 502, and Thr 505 [26,27]. It demonstrates how the presence of SWCNTs that are adjacent to the simulated B domain changes the structure of the binding site (Fig. 6a–c). The disparity in the tertiary structure of the B domain in the absence of SWCNTs is apparent (Fig. 8). The stability of the S1 subunit and the compatibility of the B domain-ACE2 binding site (purple and red circles) indicate no difference in this binding site due to the lack of changes in the S2 subunit in MD1.

Fig. 6

a-c shows the change in the form of the S1 subunit during the MD2 simulation, Figure a and c and from the side view and in b from the top, the green chain is B chain of Spike glycoprotein simulated in SWCNT presence, and the blue-red is the Spike glycoprotein (PDBID:6VYB).

Fig. 7

A-c shows the S1 and S2 subunit form changes during the MD1 simulation with circular and rectangular, respectively. Figure a and c, and from the side view and in b from the top, the purple chain is the B chain of Spike glycoprotein simulated in water, and the blue-red is the Spike glycoprotein (PDBID:6VYB).

Fig. 8

A-c shows the S1 and S2 subunit form changes during the MD1 simulation with circular and rectangular, respectively. Figure a and c, and from the side view and in b from the top, the purple and green chain is the B chain of Spike glycoprotein simulated in the water and presence SWCNT, respectively. The blue-red is the Spike glycoprotein (PDBID:6VYB).

The radial distribution functions (RDF) between the NAG groups and SWCNTs were calculated to investigate SWCNT distribution around NAG atoms. In MD2, NAG 987, NAG 988, NAG 991, NAG 992, and NAG1005, their peaks were approximately 1.2, 0.4, 1.2, 1.2, and 1 nm, respectively (Fig. 9 a). In the MD2, NAG987, NAG988, NAG991, NAG992, and NAG1005, their peaks were approximately 1.2, 0.4, 1.2, 1.2, and 1 nm, respectively. This phenomenon indicates the Lj interactions between the NAGs and SWCNT atoms. Subsequently, to further investigate the Lj interactions between SWCNTs and the NAG groups, the hydration of the NAGs groups was also determined in both systems (MD1 and MD2). As shown, the number of contacts in NAG 988, 991, and 992 with water molecules is decreased at the vicinity of SWCNTs (Fig. 9b–d). Hydration is reduced throughout the simulation time for the NAG988, as depicted in the red graph (Fig. 9b). In the NAG991 group, the hydration difference in the first 40 ns was not considerable; after this time, a significant decrease was evident in the interaction between NAG991 and water molecules (Fig. 9c). An apparent fluctuation was also observed for the contact between water molecules and NAG992, indicating the instability of NAG 992 in SWCNTs during the simulation time in both systems. Nevertheless, in MD2, after 40 ns, the reduction of the contacts between water molecules and NAG 992 was detected (Fig. 9d). The contacts of NAG987 and NAG1005 with water molecules are not noticeable, even though these two NAGs show two graphs in the RDF graph (Fig S3).

Fig. 9

(a) Radial distribution functions plots between the NAG 987, 988, 991, 992, and 1005 groups and the SWCNT. (b, c, and d the number of contacts between water and NAG 988, 991, and 992 respectively in both MD1 and MD2 systems.

Discussion

Over the last decade, carbon nanotubes as a carrier in drug delivery systems have emerged as a fundamental tool in drug delivery, which the combination of SWCNTs and proteins results in the formation of structures with unparalleled efficiency [41,42]. COVID-19 rapidly spreads all over the world from late 2019. The accessible spike glycoprotein of SARS-CoV-2 has become an essential protein in viral infection [28]. The present study investigated a single-wall carbon nanotube for the interaction with the accessible type of spike glycoprotein involved in the COVID-19 infection. As shown in Fig. 1, the S glycoprotein is composed of homotrimer and monomer peptide forms. An asymmetric reconstruction of the trimmer with a single subunit B domain is the accessible type of SARS-CoV-2 spike glycoprotein [28].

Given that the monomeric B domain of accessible type has a significant effect on the severity of pathogenicity and the importance of protein corona complex, in our study, two systems were designed to assess the interaction between SWCNTs and the B domain of SARS-CoV-2 spike glycoprotein. The stability of the B domain at the vicinity of SWCNTs is displayed in Fig. 3a, while in the absence of SWCNTs, the B domain shows an increase in the range of 28–60 ns in the root mean square deviation. Regarding the changes in the movement of the B domain (Fig. 3b), it was shown that the MSD value in MD1 was higher than MD2, while the B domain structure in MD2 was more open than MD1 (Fig. 5d), indicating an increase in the radius of gyration in MD2 (Fig. S1). As a result, the MSD value was decreased in MD2. An increase in the number of contacts after 38 ns confirms the reduction in the minimum distance between the B domain and SWCNTs (Fig. 4c and d). This phenomenon was validated by Lennard-Jones potential energy interaction between the B domain and SWCNTs. The Lennard-Jones interaction between the B domain and SWCNTs is altered slightly in random coil, beta-sheet, beta-bridge, and turn domains (Fig. 5a). These changes led to a decrease in the internal hydrogen bond of the B domain in MD2 but an increase in the interaction of water molecules with the B domain (external hydrogen bonds). An increase in the hydrophobic-hydrophilic properties of the B domain as a result of these alterations is another explanation for the Rg B domain (Fig. 5b). The adsorption of the B domain on SWCNTs can cause changes in the tertiary structure and subsequently in the binding site of the spike glycoprotein-ACE2 complex (Fig. 6, Fig. 7, Fig. 8). The effect of the presence of the NAG groups on SWCNTs increased the stability of their interactions so that the NAG 988 was in a stable interaction with SWCNTs after reducing the distance (Fig. 9a). A reduction in the interaction of the NAG group with water molecules is observed as a result of the presence of SWCNTs (Fig. 9b–d).

Conclusion

The high prevalence of SARS-CoV-2 necessitates the development of an effective drug with appropriate drug delivery as well as integration into care protocols in the world. The binding affinity between a single-wall carbon nanotube and the B domain of an accessible spike glycoprotein was evaluated in this study using all-atom molecular dynamics as an effective and efficient method. In the present study, two systems were employed, namely the simulation of the B domain of spike glycoprotein in the absence of SWCNTs (MD1) and the simulation of the B domain in the presence of SWCNTs (MD2). Simulations revealed that functionalizing the SWCNT can affect the stability of the spike protein. The Lennard-Jones potential energy between the B domain and SWCNTs was also identified. The energy exchange between the B domain and SWCNTs was discovered to improve protein stability without causing significant changes in the secondary structure of the protein. Our findings have shown that SWCNTs are adsorbed on the B domain to decrease internal hydrogen bonds, while they cause significant improvements in the solvent-accessible surface of the protein. Therefore, the results indicate that the tertiary structure of the B domain is altered in the interaction of the B domain with SWCNTs, leading to changes in the interaction of the Spike glycoprotein B domain and ACE2. All of these findings were validated using the GROMACS and VMD molecular dynamics simulations. These results can be valuable in the future for innovative drug delivery systems.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement

Data will be available on request.

Declaration of competing interest

The authors declare no conflict of interest.