Introduction

The emergence of the novel coronavirus, SARS-CoV-2, in late 2019 marked the onset of a global health crisis. Initially identified in Wuhan, China, the virus rapidly spread worldwide, causing the coronavirus pandemic 2019 (COVID-19). Characterized by symptoms such as fever, dyspnea, dry cough, and tiredness, SARS-CoV-2 has affected millions of people1, resulting in significant illness and mortality. As a member of the coronavirus family, SARS-CoV-2 is the seventh coronavirus known to infect humans and belongs to the β-coronavirus group2. Unlike SARS-CoV and MERS-CoV, which have higher mortality rates, SARS-CoV-2 poses unique challenges due to its rapid transmission3, resulting in a global pandemic.

Coronaviruses are defined by two groups of proteins: structural and non-structural4. Among the structural proteins, the Spike (S) glycoprotein plays a crucial role in virus attachment and entry into host cells5. Due to its surface exposure and central role in viral pathogenesis, the S glycoprotein is a major target for therapeutic interventions and vaccine development6. Furthermore, recent research indicates that modifying the host cell receptor could enhance its binding ability and assist in preventing the virus from entering cells7.

However, the extensive coverage of S glycoproteins with N- and O-linked glycans presents challenges in targeting specific epitopes for neutralizing antibodies8,9.

The SARS-CoV-2 S glycoprotein is a trimeric protein, with each monomer containing two functional subunits: S1 and S2. The S1 subunit is responsible for binding to the host cell receptor, while the S2 subunit facilitates the fusion of the viral and cellular membranes10. The S-glycoprotein of SARS-CoV-2 directly interacts with the angiotensin-converting enzyme 2 (ACE2) on target cells11. The S1 subunit comprises four domains: the N-terminal ___domain (NTD), receptor-binding ___domain (RBD), and two C-terminal domains (CTD1 and CTD2). The RBD, located in the distal ___domain of the S1 subunit, helps stabilize the prefusion state of the S2 subunit, which contains the fusion machinery12. To interact with ACE2, the RBD of the S1 subunit undergoes hinge-like movements, transitioning between “down” (inaccessible) and “up” (accessible) conformations (see Fig. 1)13,14. These conformational changes are critical for the engagement of host cell receptors14 by the SARS-CoV-2, SARS-CoV, and MERS-CoV viruses12.

Fig. 1
figure 1

SARS-CoV-2 S protein in the down state (left) and up state (right). The RBD ___domain is depicted in yellow, while the rest of the SARS-CoV-2 S trimer is shown in cyan.

Full size image

Each unit of the SARS-CoV-2 S protein contains 19 N-linked and 1 O-linked glycosylation sites, i.e., reacting sites with glycans15,16. Molecular dynamics (MD) simulations have underscored the importance of these glycans in inducing conformational changes in the SARS-CoV-2 S protein9,17,18,19. The conformational masking and glycan shielding of the SARS-CoV-2 S glycoprotein in the down conformation reduce its recognition by the immune response of the infected host14,17.

Generally, two classes of RBD exposure have been reported for the wild type of the S trimer, i.e., all RBDs down, or one up20. A 1:1 ratio of these two classes has been observed using cryo-EM21. During the transition from the down to the up conformation, the fusion-peptide proximal region (FPPR) clashes with the CTD1, which rotates outward along with the RBD. Therefore, a structured FPPR plays an important role in stabilizing the down conformation22, although it occasionally flips out of position, allowing the RBDs to move into the up conformation. Investigating the conformational transition of the SARS-CoV-2 S trimer at the atomic scale helps identify key residues that significantly influence these conformational changes.

It has been proposed that the SARS-CoV-2 virus may trigger a chronic inflammatory response in host cells, thereby contributing to persistent oxidative stress23. A key element in this process is the ACE2 receptor, which plays a critical role in countering oxidative stress. Membrane-bound ACE2 functions by degrading Angiotensin II, a molecule that otherwise promotes the generation of reactive oxygen and nitrogen species (RONS) via membrane-associated NADPH oxidase. The breakdown of Angiotensin II results in the production of Angiotensin 1–7, which serves to inhibit NADPH oxidase activity. However, when the ACE2 receptor binds to the SARS-CoV-2 S-glycoprotein, its capacity to degrade Angiotensin II is compromised. This impairment leads to elevated RONS levels at the cell surface, exacerbating oxidative stress24. Notably, RONS have the ability to oxidize and nitrate amino acids within proteins, which can interfere with their normal functional roles25,26,27,28,29,30. Experiments using an external RONS source (i.e., cold atmospheric plasma) demonstrated that highly reactive amino acids like Met and Cys are predominantly oxidized in aqueous solutions31,32. Therefore, oxidation of these amino acid residues due to oxidative stress by viral infection may interfere with the normal exposure of the SARS-CoV-2 RBD, rendering it easier accessible to the cell receptor and ultimately increasing cell infection33.

In this study, we investigate the conformational changes of the SARS-CoV-2 RBD from the down to the up conformation using targeted molecular dynamics (TMD) simulations, identifying important residues involved in these transitions. By altering specific Cys residues, which are easily oxidized during oxidative stress31,32,34,35, we examine the potential disruption of normal conformation of the SARS-CoV-2 RBD and its accessibility to the receptor. Cysteic acid (CYO) and cystine (formed by a disulfide bond between two cysteine residues) are products of Cys oxidation by hydrogen peroxide (H2O2)36, one of the RONS generated during oxidative stress in response to viral infections24. Generally, intracellular RONS, including H2O2, can oxidize around 5% of Cys residues of proteins to cysteic or Cys sulfenic acid37, and this effect is enhanced under oxidative stress. Furthermore, the most significant amino acid residues and glycans playing an important role in conformational changes in both the native and oxidized states are identified and compared. Finally, using an advanced sampling method known as umbrella sampling (US)38, we perform a comparative analysis of the free energy profiles (FEPs) associated with the translocation between two metastable conformations (i.e., down and up states). These simulations offer crucial insight into how low levels of oxidation influence the conformational transition mechanisms of the SARS-CoV-2 S glycoprotein.

Computational details

The SARS-CoV-2 S trimer is a large molecule, with each monomer (chains A, B, and C) comprising 1273 amino acid residues. Different conformations of its down and up states are available in the Protein Data Bank (PDB). We selected the cryo-EM structures of the SARS-CoV-2 S trimer for the down and up conformations, corresponding to PDB IDs 6VXX14 and 6VSB13, respectively. Since the RBD part has numerous missing residues in both PDB structures, we used the crystal structure of the RBD from PDB ID 6M0J39 to remodel the structures.

To glycosylate the SARS-CoV-2 S trimer, we employed the CHARMM-GUI web server40,41,42. Based on the literature14,15, we selected the most probable glycan structures for both systems. Details of the selected glycans can be found in Table S1 of the Supporting Information (SI).

After preparing the native down and up conformations of the SARS-CoV-2 S trimer, MD simulations were conducted using the Gromacs 2020.2-MODIFIED software43 to equilibrate both structures. The CHARMM36 force field44,45 was employed to describe the interatomic interactions, and the initial structures (model systems) were prepared using the CHARMM-GUI web server40.

The simulation box dimensions were 201 × 192 × 198 Å3 for the down system and 192 × 185 × 204 Å3 for the up system. Periodic boundary conditions were applied in all three spatial directions. The charge of both systems was neutralized by adding chloride and sodium ions to the surrounding water, and explicit water molecules were simulated using the TIP3P model46. Energy minimization was carried out using the steepest descent algorithm for 20,000 steps.

For both systems (i.e., native down and up conformations), three replicas were prepared with different initial atomic velocities using different random seeds. Equilibration was then performed in the NVT ensemble at 310 K for 3 ns with a time step of 1 fs. Subsequently, they were relaxed in the NPT ensemble at 1 atm and 310 K for 200 ns with a 2 fs time step (see the root mean square deviation (RMSD) plots in Figure S1A of the SI). The Nosé-Hoover thermostat47,48,49, with a coupling constant of 1 ps and the isotropic Parrinello-Rahman barostat50, with a coupling constant of 5 ps and a compressibility of 4.5 × 10−5 bar−1, were used to equilibrate the systems.

To maintain the orientation of the SARS-CoV-2 S trimer within the simulation box, the tails of the SARS-CoV-2 S trimer were restrained using a harmonic potential with a force constant of k = 4000 kJ mol−1 nm−2, specifically restraining the heavy atoms in residues 1140–1147. Long-range interactions were computed using the Particle Mesh Ewald (PME) method51,52 with long-range dispersion corrections applied for both pressure and energy. The Verlet list scheme was employed, with a 12 Å cutoff for both electrostatic and van der Waals (VDW) interactions. Considering the cutoff value, interactions occurring within a 12 Å radius between the RBD ___domain of chain A and the surrounding amino acid residues in the down conformation of the SARS-CoV-2 S trimer were assumed to influence structural changes (highlighted in red in Fig. 2A,B).

Fig. 2
figure 2

Structural and solvent accessibility analysis of the SARS-CoV-2 S protein in the down state. (A) The RBD ___domain is highlighted in yellow, with its immediate surrounding region within a 12 Å radius shown in red to emphasize its structural context. All cysteine residues within these domains are marked in purple and represented with VDW spheres. (B) Detailed view of the RBD and its surrounding region within a 12 Å radius, focusing on the areas highlighted in yellow and red. Four cysteine residues demonstrating the highest solvent accessibility, determined by SASA analysis, are again depicted in purple with VDW representation. (C) SASA analysis plot (with corresponding errors in pale color) illustrating the solvent accessibility of all cysteine residues within the yellow and red domains. The four cysteines with the highest SASA values are highlighted with light blue circles (i.e., Cys480-A (Cys480-RBD), Cys166-B, Cys480-C, and Cys488-C), indicating their potential importance and susceptibility to oxidation based on detailed analysis.

Full size image

During oxidation, modifications of amino acid residues within this ___domain (i.e., the region surrounding and including the RBD) can potentially impact conformational changes. Furthermore, the conformational masking and glycan shielding characteristic of the down conformation of the SARS-CoV-2 S trimer expose certain amino acids to the solvent, making them susceptible to oxidation by reactive oxygen species (ROS) generated in the environment. Thus, to identify the amino acid residues most susceptible to oxidation within this ___domain, a solvent accessible surface area (SASA) analysis was conducted using the gmx sasa tool in Gromacs on the last 10 ns of all three replicas of the native down conformations. The SASA analysis was focused on identifying amino acids within the ___domain that are highly accessible to the solvent, as these are more likely to undergo oxidation. However, as there are no Met residues in the selected ___domain, Cys residues, which can be oxidized to CYO, were considered next in susceptibility to oxidation. Based on the SASA results (Fig. 2C), Cys480 in the RBD ___domain, Cys166 in chain B, and Cys480 and Cys488 in chain C were identified as susceptible to be modified to CYO.

The partial charges of CYO residues were identified by determining their protonation states in a physiological environment. The negative partial charge of CYO under physiological conditions is given by the DrugBank database53. Additionally, experimental investigations on Cys oxidation have also confirmed that the oxidation of Cys residues in proteins, induced by ROS, predominantly results in negatively charged CYO residues37.

For negatively charged CYO force-field parameters, a combination of Gaussian 16 software54 and the CHARMM General Force Field (CGenFF)55 was employed, following56. Density Functional Theory (DFT) was applied to optimize and generate the partial charges for CYO using the B3LYP functional with the standard 6-311G* basis set in Gaussian. The resulting partial charges and optimized structure were used to construct topology files for CYO, compatible with the CHARMM36 force field, using the CGenFF program.

Similar to the methods used for the native systems, both the down and up conformations of the oxidized SARS-CoV-2 S trimer were constructed and equilibrated (see the RMSD plots in Figure S1B of the SI). Thus, three replicas were generated for each of the following systems: (i) native down S trimer, (ii) native up S trimer, (iii) oxidized down S trimer, and (iv) oxidized up S trimer. This resulted in twelve systems in total.

Given the substantial energy barriers associated with the RBD transition, conventional MD methods would require impractically long simulation times. Therefore, TMD simulations were employed to accelerate these transitions, while maintaining minimal system perturbation. Using PLUMED-2.6.0 software57,58,59, TMD simulations were performed to transition the SARS-CoV-2 S trimer between its down and up conformations. From the last 100 ns of the trajectories, three structures from each of the three replicas were selected as initial and target conformations for the TMD simulations, resulting in a total of nine structures per system (i.e., native down, native up, oxidized down, and oxidized up S trimers). During the TMD simulations, each of the nine systems in the down state was randomly targeted to one of the nine structures in the up state. Consequently, for each native and oxidized system, nine TMD simulations were conducted for the transition from down to up. Each TMD simulation was run for 5 ns to guide the SARS-CoV-2 S trimer from the down to the up state. This was achieved by applying a spring with a constant free length \(\rho = 0.05\text{ \AA }\) and a linearly increasing elastic constant ranging from 0 to 80,000 kJ mol−1 nm−2. In the TMD simulations, all heavy atoms were subjected to spring forces to guide them from the initial to the final structure.

To investigate the pathway between up and down states, the VMD software60 was employed. Initially, the number of hydrogen bonds (H-bonds) and salt bridges between the RBD and its surroundings in the down conformation of the SARS-CoV-2 S trimer were computed. Subsequently, differences in the pathways of H-bond and salt bridge dynamics during the down-to-up transition between the native and oxidized systems were compared. Quantitative analyses of H-bonds and salt bridges were conducted based on the methodology described by Debiec et al.61, adopting specific criteria: a distance of 3.5 Å for salt bridge formation and 4.5 Å for dissociation, along with a 3 Å cutoff distance for hydrogen bond formation62. The quantities of H-bonds and salt bridges were determined by averaging their occurrence over the last 100 ns of MD simulations, where each interaction was prevalent in more than 10% of all frames.

The US simulation method38 was employed to calculate the FEPs for the translocation of the RBD from the down to the up conformation. For each system (native and oxidized), three trajectories from TMD simulations were selected to extract US windows. The TMD trajectories navigate the complete reaction pathway from the down to up conformation. As a result, the US sampling windows, derived from these trajectories, effectively encompass the entire pathway. According to Fig. 3, the RMSD of the RBD stabilizes after 2 ns for the native system and 5 ns for the oxidized system, demonstrating that the system has attained its final state with minimal fluctuations.

Fig. 3
figure 3

RMSD of the RBD of the SARS-CoV-2 S trimer during the transition from the down to the up conformation for the native (blue) and oxidized (red) systems, averaged over 9 replicas, with errors shown in pale color. The dashed-line ellipsoid illustrates the two distinct proposed regimes.

Full size image

The initial and final states for the TMD simulations were chosen from separate equilibrated replicas: the initial state was derived from three equilibrated down-state replicas, while the target structure was selected from three equilibrated up-state replicas.

We employed almost 30 US windows with a separation of \(<0.1\text{nm}\) for each trajectory of native and oxidized system to build a single FEP. The RBD ___domain was moved from its down to up conformation. To keep the RBD species at the position with respect to the rest of SARS-CoV-2 S protein, a harmonic bias with a force constant of 1000 kJ mol-1 nm-2 was applied to restrict their center of mass (COM) motion along the reaction coordinate. Similar to the equilibration runs mentioned above, US simulations were performed at NPT ensemble for each system, consisting of 5 ns of equilibration and 10 ns of sampling. The umbrella histograms were collected, and an individual FEP was constructed using the weighted histogram analysis method (WHAM)63. The initial and final distances between the COM of two groups (i.e., the RBD and the rest of the protein) for the final FEPs were determined by averaging their COM distances in the equilibrated systems for the down and up conformations where the SARS-CoV-2 S trimer is metastable. In the native system the minimum and maximum distance were equal to 4.734 and 6.645 nm, respectively. In the oxidized system the minimum and maximum distance were equal to 4.729 and 6.626 nm, respectively. The final FEP for each system was obtained by averaging over 3 individual FEPs (using 3 trajectories, see above). Thus, to obtain the final FEPs, 3 × 30 = 90 US windows were used for each native and oxidized system.

Results and discussion

The RMSD values of the RBD during the transition from down to up state revealed that there are different regimes of transitions in the oxidized system (Fig. 3). The averaged RMSD slope distinguishes between rapid and slow transition regimes.

For the first 1.5 ns, the system undergoes rapid changes, while from 1.5 to 4 ns, it moves much slower toward the final state. A detailed investigation of the RMSD is needed to realize the behavior. To further investigate the structural changes observed in the RMSD analysis, we examined specific RBD regions involved in the transition process. As stated in Ref. 3333, the RBD can be indicated in 5 non-overlapping regions: I (Cys336–Cys361), II* (Cys379–Leu390), II&III (Cys391–Cys432), III* (Val433–Pro479 and Tyr489–Cys525), and IV (Cys480–Cys488). Regions III* and IV are particularly crucial as they interact with host cell receptors for cell entry and fusion. The RMSD of these 5 regions during the transition for both the native and oxidized systems show that the two distinct phases observed in the oxidized system are associated with these two domains. Figure 4A,B illustrate the conformational fluctuations of these regions in both the native and oxidized systems, providing insights into the underlying mechanisms driving the transition. The magenta and green curves represent regions III* and IV, respectively, which exhibit significant changes in the oxidized system compared to the native one, highlighting their role in the observed two-regime RMSD behavior in the oxidized system.

Fig. 4
figure 4

RMSD of each region of the native (A) and oxidized (B) RBD of the SARS-CoV-2 S trimer during down-to-up transition, all averaged over nine replicas.

Full size image

Further investigation into the trajectories of TMD simulations revealed that glycans, notably N165B, N343A, and N343B, play a critical role in the slower transition regime observed in the oxidized RBD. Analysis of H-bonds in the down state of both systems showed that in the native system, five glycans, N331A, N343A, N165B, N234B, and N343B, significantly contribute to forming H-bonds (Tables S2, S3, and S4 in the SI), whereas in the oxidized system, only N234B and N343B form H-bonds with the RBD. During the transition, the presence of negatively charged CYO residues plays a key role in attractive interactions with these glycans. In both systems, N234B and N331A, located at the hinge ___domain, remained stable relative to RBD movements. However, N165B was pulled upwards by the RBD in both the native and oxidized systems. A notable difference between the two systems was observed with N165B, which detached from the RBD around 1.6 ns in the TMD simulation of the native system but remained attached throughout the entire TMD simulation in the oxidized system. Figure 5A,B depict the positions of these glycans and the RBD ___domain during the transition from the down to the up state in the native and oxidized systems, respectively. From these figures, it is evident that there were no significant differences between the two systems in terms of N234B and N331A (depicted in white and licorice representations). The distinction primarily lies in the interactions involving oxidized Cys residues (depicted in VDW representation) and glycans N343A and N343B. Specifically, N343A in the RBD was attracted to CYO480 and CYO488 of chain C, hindering easy movement of the RBD. Moreover, CYO480 in the RBD, initially not attracted to any glycans, found N343B during the transition and rotated towards it, preventing complete upward movement. The second phase of RMSD observed after approximately 1.5 ns could be attributed to these interactions, as CYO480 is positioned at the boundary between region III* and region IV. To visualize the dynamic behavior of the system, trajectories from TMD simulations were converted into video format.

Fig. 5
figure 5

The position of RBD decorating glycans in relation to the rest of the SARS-CoV-2 S trimer during down-to-up transition for the native (A) and oxidized (B) systems. Region III*, IV, and the rest of the RBD are illustrated in magenta, green, and yellow colors, respectively. The rest of the SARS-CoV-2 S trimer, is depicted in cyan color. The glycans N331A and N234B are represented in white using licorice representation. The modified Cys residues in the oxidized system (see text) are shown with VDW representation. To enhance clarity on the positions of the RBD, glycans, and (native/oxidized) Cys residues, the top views of both the initial (0 ns) and final (5 ns) states of each system are zoomed in. The position of the host cell receptor is indicated by a dashed circle.

Full size image

These videos, available online at https://github.com/maryamghasemitarei/Transion-of-SARS-CoV-2-Spike-protein, offer insights into the conformational changes and interactions, helping in a deeper understanding of the structural evolution during the transition.

The average total number of H-bonds and salt bridges formed between the RBD ___domain and the remaining SARS-CoV-2 S trimer (with a prevalence exceeding 10% in all frames) is summarized in Tables 1 and 2, respectively. The results indicate that after oxidation, there was no significant change in the number of H-bonds, but there was an increase in the number of salt bridges. This increase suggests that the transition process becomes more challenging in the oxidized protein.

Table 1 Number of H bonds formed between RBD and each chain of SARS-CoV-2 S protein (chains A, B and C), along with the total H bonds. Values are averaged over the last 100 ns of MD simulations, representing bonds present in more than 10% of the frames.
Full size table
Table 2 Number of salt bridges formed between RBD and each chain of SARS-CoV-2 S protein (chain A, B and C), along with the total number of salt bridges. Values are averaged over the last 100 ns of MD simulations, representing bridges present in more than 10% of the frames.
Full size table

In the interface between all pairs forming H-bonds (Tables S2, S3, and S4 in the SI), a total of fourteen pairs were identified within the hinge ___domain. These pairs exhibited minimal movement during the transition from the down to the up conformation, and the distance between their donor and acceptor atoms did not significantly change during this transition (Figure S2 in the SI). During the transition from the down to the up state, the remaining H-bonds between RBD and its surrounding SARS-CoV-2 S trimer were broken and the distance between donor and acceptor atoms in both the native and oxidized systems were increased (Figure S2 in the SI). Similarly, among all pairs forming salt bridges (Table S5 in the SI), four pairs—Arg328-Asp578, Arg319-Asp737, Lys535-Glu583, and Glu471-Lys113—were identified within the hinge ___domain, with minimal change in nitrogen–oxygen distance during the transition (Figure S3 in the SI). However, salt bridges in other pairs—Lys462-Asp198, Lys535-Glu554, and Asp427-Lys986—were disrupted during the transition. The differences in the hydrogen bonding and salt bridge interactions (Tables 1 and 2) align with the observed variations in transition dynamics. Specifically, the increased number of salt bridges in the oxidized system suggests a stronger stabilization of the RBD in its initial conformation, potentially making the transition process more challenging. However, the overall reduction in energy barriers, as seen in the umbrella sampling results, indicates that oxidation ultimately facilitates the transition by altering key stabilizing interactions.

Next, we turn to the results from the US simulations. In Fig. 6, we show the sampled FEPs of RBD translocation. The blue and red curves show the energy profiles for the native and oxidized systems, respectively.

Fig. 6
figure 6

Free energy profiles for the translocation of RBD from down to up state for the native and oxidized structures (oxidation of Cys to CYO). Errors associated with the sampling are presented as shaded regions.

Full size image

As shown in Fig. 6, the transition path along the reaction coordinate (ξ) is relatively complex in both the native and oxidized structures. The total maximum energy barriers for the native and oxidized systems are 62 ± 13 kJ mol⁻1 and 29 ± 11 kJ mol⁻1, respectively. This indicates that the oxidation of cysteine residues to CYO reduces the free energy barrier by approximately 30 kJ mol⁻1, thereby significantly increasing the probability of RBD translocation between the down and up conformations at human body temperature as compared to the native state. Additionally, the oxidized profile (red curve) shows an intermediate local minimum around 6 nm, which could be related to the two regimes observed in the TMD results. The US results (Fig. 6) quantitatively confirm the observations from our TMD simulations. The reduction in free energy barriers in the oxidized system is consistent with the structural alterations observed in the RBD transition pathway. While the increased number of salt bridges initially stabilizes the oxidized RBD, the overall decrease in the energy barrier indicates that oxidation ultimately facilitates the transition. This highlights a dual effect of oxidation: localized stabilization due to increased salt bridges and an overall enhancement of transition probability due to energy barrier reduction. These findings provide a more comprehensive understanding of how oxidation impacts the RBD conformational shift, linking structural perturbations with energetic considerations.

Summary and conclusions

Although RONS are initially produced to combat viral infections, many viruses, including SARS-CoV-2, have evolved mechanisms to exploit this hostile environment. In the case of SARS-CoV-2, oxidative stress is believed to play a central role in the pathogenesis of the disease24. In this study, we have described how oxidation of SARS-CoV-2 S trimer surrounding the RBD ___domain alters its conformational changes between the down and up states and makes the transition energetically more favorable by reducing the relative energy barrier by approximately 30 kJ mol-1. Excessive levels of RONS generated during viral infection can oxidize Cys residues into cysteic acid, an irreversible product of cysteine oxidation64. In our study, the native structure of SARS-CoV-2 S trimer contained Cys in the form of cystine. However, in the oxidized structural complexes, the cystine residues within the 12Å of the RBD were irreversibly converted to cysteic acid. Moreover, when extracellular reactive species are produced during infection, only the oxidation of cystine residues to cysteic acid is possible. The significant reduction of the transition barrier due to oxidation of the thiol group of RBD and its surrounding is thus beneficial for the viral infection.

To the best of our knowledge, this is the first study that explored the impact of cysteine oxidation (i.e., the formation of cysteic acid) on the conformational changes of the SARS-CoV-2 S protein. Additionally, while an excessive level of RONS is required for the oxidation and nitration of aromatic residues, this aspect was not addressed in our study. SARS-CoV-2 exploits host cell machinery to destabilize the redox balance, depleting antioxidants and promoting an environment conducive to viral replication24. Our findings propose an additional mechanism through which oxidative stress induced by SARS-CoV-2 could facilitate viral entry: the oxidation of the S protein trimer may enhance its exposure to cell receptors.

This study provides essential insights into the structural dynamics and conformational changes of the SARS-CoV-2 S protein, focusing on its transitions between the down and up states and the role of cysteine oxidation in these processes. These conformational shifts are critical for viral attachment to host cells, underscoring their importance as potential therapeutic targets.

Using SASA analysis, we identified key residues susceptible to oxidation in the presence of ROS as a result of viral infection. We examined how these modifications affect the SARS-CoV-2 S protein’s conformational transitions through TMD simulations.

The investigation highlights the critical role of glycans in the change of transition dynamics observed in the oxidized RBD during TMD simulations. The main difference lies in oxidized Cys residues interacting with glycans, which hinders RBD movement and prevents complete upward transition. These interactions contribute to the observed changes in RMSD along the TMD simulations.

Additionally, our analysis of H-bonds and salt bridges before and after oxidation shed light on how these changes contribute to stabilizing the SARS-CoV-2 S protein’s structure in each state. While the overall number of H-bonds remained relatively unchanged post-oxidation, we observed an increase in salt bridge formation, further enhancing the stabilization of the protein’s conformation.

The US simulations revealed a reduction in the energy barrier by approximately 30 kJ/mol, thereby significantly enhancing the transition rate in the oxidized system. However, pathway analysis from TMD simulations indicates that the presence of oxidized cysteines and their interactions with other residues and glycans make the transition process more challenging. It is important to note that this challenge stems from the transition pathways observed in TMD simulations and does not conflict with the US simulation results, which highlight the increased transition rate due to the reduced energy barrier. Furthermore, based on the equilibrium up conformation of the oxidized system and analysis of the TMD trajectories, the metastable up conformation of the oxidized system positions the RBD in a partially exposed state as compared to the native system. Therefore, the US simulations focused on the transition rate between the down conformation and this new metastable up conformation. However, the potential interaction of this new up conformation with ACE2 was not explored in this study. Based on the position of the RBD in the oxidized up conformation, it can be assumed that the RBD is sufficiently exposed for potential attachment to ACE2.

This research advances our understanding of the SARS-CoV-2 S protein’s structural dynamics and underscores the potential of targeting post-translational modifications like oxidation for developing antiviral strategies. By elucidating the energy landscape of conformational transitions and identifying key residues affected by oxidation, our findings provide valuable insights for designing therapeutic interventions aimed at combating SARS-CoV-2 infection and transmission.