Computational Investigation of Structural Dynamics of SARS-CoV-2 Methyltransferase-Stimulatory Factor Heterodimer nsp16/nsp10 Bound to the Cofactor SAM

Introduction

Coronaviruses (CoVs) are considered as an etiological agent for causing severe acute respiratory syndrome (SARS) in humans. In the past two decades, SARS-CoV and MERS-CoV (middle east respiratory syndrome coronavirus) were responsible for the epidemic in 2003 and 2012, respectively. Recently, in December 2019, the novel coronavirus (COVID-19) pandemic, caused by SARS-CoV-2, first broke out in Wuhan city of China and has been spreading worldwide (Bogoch et al., 2020; Guan et al., 2020; Lin X. et al., 2020). So far, ~42 million people worldwide are infected by SARS-CoV-2, including ~1.1 million deaths. In India, the COVID-19 tally has already crossed 7.7 million, including ~0.1 million deaths. However, neither prophylactic vaccines nor any direct antiviral drugs are available to effectively treat the human and animal coronavirus disease (Lu, 2020; Pillaiyar et al., 2020; Sheahan et al., 2020).

CoVs are single-stranded positive-sense RNA viruses (Eckerle et al., 2010; Fehr and Perlman, 2015) belonging to the family Coronaviridae and possess the largest genome (26.4–31.7 kb) (Woo et al., 2010). They are classified into four genera, namely Alphacoronaviruses, Betacoronaviruses, Gammacoronaviruses, and Deltacoronaviruses (King et al., 2011). CoVs can infect both humans and animals (Lun and Qu, 2004; Coleman and Frieman, 2014), and can cause diseases like Hepatitis and Pneumonitis in mouse (Weiss and Leibowitz, 2011) and neurologic & respiratory diseases in humans (Arabi et al., 2015; Huang et al., 2020). So far, seven distinctive strains of human coronaviruses (HCoV) have been disclosed that includes 229E and NL63 (Alphacoronaviruses), and OC43, HKU1, SARS-CoV (2002), MERS-CoV (2012), and SARS-CoV-2 (2019) (Betacoronaviruses) (Weiss and Navas-Martin, 2005; Zeng et al., 2018; Singh et al., 2020). SARS-CoV-2 shares a strong correlation with the bat CoV RaTG13 having a 96.2% genome sequence identity, suggesting its possible evolution from bat CoVs (Wu et al., 2020; Zhou P. et al., 2020). SARS-CoV-2 displays a sequence similarity of 80% with SARS-CoV, whereas only 50% with MERS-CoV (Lu et al., 2020; Wu et al., 2020; Zhou Y. et al., 2020). Currently, most of the therapeutic options that are available for controlling COVID-19 is based on the previous knowledge and information gathered from SARS-CoV and MERS-CoV.

Most viruses or eukaryotic cellular mRNA possess the 5′-end capping mechanism that plays a vital role in mRNA splicing, translation initiation, stability, and intracellular RNA transport (Furuichi and Shatkin, 2000). The capping of the 5′-end occurs through a sequential enzymatic process. This involves three enzymes, such as RNA guanylyltransferase (GTase), RNA triphosphatase (TPase), and RNA guanine-N7-methyltransferase (N7-MTase). This generates a cap-0 structure (m7GpppN). It is further methylated at the 2′-O position of mRNA by 2′-O- methyltransferase (2′-O-MTase) and generates the cap-1 (m7GpppNm) and cap-2 (m7GpppNmNm) structures. This mimicking of the eukaryotic mRNA capping mechanism ultimately helps the virus evade the host innate immune system (Furuichi and Shatkin, 2000; Wang et al., 2015). Both the MTase uses S-adenosyl-L-methionine (SAM or AdoMet) as a donor of methyl and gives a by-product, S-adenosyl-L-homocysteine (SAH or AdoHcy).

The genome of SARS-CoV-2 is comprised of 10 ORFs (Open Reading Frame). ORF1ab encodes the replicase polyprotein 1ab (PP1ab), which gets cleaved by the two viral proteases, PL^pro (papain-like protease) and 3CL^pro (3-C like protease) at the N-terminus and C-terminus, respectively, to form all the 16 non-structural proteins (nsp1, nsp2, nsp3 by PL^pro, and nsp4-nsp16 by 3CL^pro). The remaining ORFs are associated with encoding structural proteins, such as spike (S), envelope (E), membrane (M), and nucleocapsid (N) proteins, as well as other accessory proteins (Harcourt et al., 2004; Yang et al., 2005; Chen et al., 2020; Wu et al., 2020). Previous structural and biochemical characterization studies on SARS-CoV (2002) showed that in the nsp10-nsp16/SAM complex, nsp10 enhances the methylase activity of nsp16 by increasing the stability of the SAM-binding pocket (Chen et al., 2011; Wang et al., 2015).

The present study involves analyzing protein-protein interactions in the heterodimeric structure formed by nsp10 and nsp16 of SARS-CoV-2. The recently solved crystal structure of SARS-CoV-2 nsp16/nsp10 heterodimer bound to SAM (PDB: 6W4H) (Minasov et al., 2020) is used in the current study. The three-dimensional structure of the complex is shown in Figure 1. Here, we have used the standard molecular dynamics (MD) simulations in conjunction with the molecular-mechanics Poisson-Boltzmann surface area (MM-PBSA) method to identify the critical residues involved in the complex formation. Further, we conducted MD simulations of nsp16/SAM, elucidating the importance of nsp10 in the binding of SAM to nsp16. The detailed structural analysis and inter-molecular interactions reveal the interface of nsp16/ns10, which provides a distinctive feature for coronaviruses, can be exploited as an attractive target in developing specific antiviral drugs for controlling COVID-19. MD simulations were conducted to elucidate the role of nsp10 in the interaction of SAM to nsp16. Also, the conformational changes in nsp10 and nsp16 due to their association were investigated.

FIGURE 1

Figure 1. Crystal structure of SARS-CoV-2 nsp10-nsp16 complexed with cofactor S-Adenosyl Methionine (SAM) (PDB: 6W4H). The nsp10, interacting part of nsp10, interacting part of nsp16, nsp16, and the cofactor SAM in a complex is shown in color teal, gold, dark sky blue, light green, and yellow, respectively. The surface represents an electrostatic surface.

Materials and Methods

We retrieved the experimental coordinate of the SARS-CoV-2 nsp16/nsp10 complex from Protein Data Bank (PDB), which was crystallized at 1.8 Å (PDB ID 6W4H) (Minasov et al., 2020). The complex structure includes S-adenosyl-L-methionine (SAM). The protonation states were determined using PROPKA 3.1 (Olsson et al., 2011), and the corresponding residues were modified accordingly. The following three systems were simulated for the current study: complex (nsp16/nsp10/SAM), nsp16_SAM (nsp16/SAM), and nsp10_apo (apo nsp10). The starting configurations of nsp16_SAM and nsp10_apo were constructed manually from nsp16/nsp10/SAM (PDB: 6W4H).

Missing hydrogens in the crystal structure were built using the Leap module of AmberTools19 (Case et al., 2018), and a proper amount of Na⁺ was added to neutralize the system. All the systems were solvated in a periodic octahedron TIP3P water box (Price and Brooks, 2004) with a 10 Å buffering distance from all directions. The proteins (nsp10 and nsp16) were assigned the Amber ff14SB force field parameters (Maier et al., 2015) while the cofactor SAM was modeled using the updated Generalized Amber Force Field (GAFF2) (Wang et al., 2004) and AM1-bcc (Jakalian et al., 2002) charges. The bond lengths having hydrogen atoms were kept fixed using the SHAKE algorithm (Kräutler et al., 2001), which allowed 2 fs time-step for MD simulations. The Particle-Mesh Ewald (PME) summation (Darden et al., 1993) scheme was employed to compute the long-range interactions. The non-bonded cut-off was set to 10 Å.

Firstly, the solvated systems were subjected to an energy minimization using 500 steps of steepest descent, followed by another 500 steps of the conjugant gradient algorithm. During the minimization, the solute atoms were kept fixed with a restraint force of 2.0 kcal mol⁻¹Å⁻². The second stage of minimization was carried out without any restraint force. After the minimization, each system was gradually heated from 0 to 300 K in the NVT ensemble, and the solute atoms were restrained with force constant of 2.0 kcal mol⁻¹Å⁻². The temperature was maintained using a Langevin thermostat (Pastor et al., 1988; Loncharich et al., 1992). Subsequently, each system was simulated for 50 ps at 300 K at a constant pressure of 1 atm using the Berendsen barostat (Berendsen et al., 1984). In this stage also, the same restraint force was applied to the solute atoms. Before the production run, we equilibrated each system by conducting 1 ns MD simulation in the NPT ensemble without restraining the system. Finally, the production simulation was carried out for 1 μs using the pmemd.cuda module of AMBER18. Overall, 100,000 snapshots were generated and used for the analysis using the Cpptraj module (Roe and Cheatham, 2013) of AMBER18 (Case et al., 2018).

To study the effect of salt concentration on the binding, we simulated the complex system at 300 K for 1 μs at two different salt concentrations (0.15 and 0.25 M). The desired salt concentration was achieved by adding an appropriate number of monovalent ions (Na⁺ and Cl⁻) obtained with the protocol suggested by Machado and Pantano (2020). Similarly, we have also investigated the effect of temperature on the binding by conducting MD simulation of the complex at a relatively elevated temperature of 310 K.

Trajectory Analysis

All analyses, including root-mean-square deviations (RMSDs), root-mean-square fluctuations (RMSFs), radius of gyration (R_g), solvent accessible surface area (SASA), etc. were performed using the Cpptraj module (Roe and Cheatham, 2013) of AMBER18. A distance of ≤ 3.5 Å and an angle cut-off of ≥120° were used for hydrogen bond calculations. The same criterion was used in earlier studies (Jeffrey, 1997; Hu et al., 2016; Shi et al., 2018; Sanachai et al., 2020).

Principal Component Analysis

One of the widely used unsupervised data reduction schemes is principal component analysis (PCA) (Ichiye and Karplus, 1991). First, a covariance matrix (C) is constructed from the atomic fluctuations of C_α-atoms of each residue (Amadei et al., 1996). The elements C_ij of the covariance matrix C are defines as

$\begin{array}{ll}{C}_{ij}=<((x_i-<x_i>)(x_j-<x_j>))>& (1)\end{array}$

where x_i and x_j denotes the instant coordinates of the ith or jth atoms, while < x_i > and < x_j > denote the mean of ith or jth atoms over the ensemble. The diagonalization of the covariance matrix yields orthogonal eigenvectors and the corresponding eigenvalues. The resulting eigenvectors are termed as principal components (PCs) and indicate the direction of the movement, while the corresponding eigenvalue describes the amplitude of motions. We adopted the same methodology for PCA, as discussed in our earlier studies (Jonniya et al., 2019; Sk et al., 2020c).

The free energy landscape (FEL) was generated based on the following Equation (Frauenfelder et al., 1991):

$\begin{array}{ll}{G}_i=-k_{B}T\text{ ln}\left(\frac{N_i}{N_m}\right)& (2)\end{array}$

where k_B represents the Boltzmann constant, and T is the absolute temperature. N_i is the population of the ith bin, and N_m is the population of the most populated bin. The 2-dimensional FEL was constructed using PC1 and PC2 as the reaction coordinate.

Residual Network Analysis of Protein

The residual network analysis approach is widely used to explore the viral fitness and resistance development of protein structure. The network analysis of protein structures (NAPS) (Chakrabarty et al., 2019) server (http://bioinf.iiit.ac.in/NAPS/) was used to identify key residue interactions in the residual network and the network-based hydrophobic contacts from the simulation trajectories. Here, we used the protein-protein complex option, followed by the C_α network type and unweighted edge weight. An edge represented the distance between a pair of C_α atoms within the lower and upper thresholds (default upper threshold = 7 Å; lower threshold = 0 Å). We considered the default parameters for the long-range interaction networks and minimum residue separation of 1. The hydrophobicity indices of 20 amino acids range from −4.5 to 4.5 (Kyte and Doolittle, 1982). These values were colored in gradients of red to show the hydrophobic residues, while the gradients of blue were used to display the hydrophilic residues.

Binding Free Energy and Alanine Scanning

The interaction energy between nsp10 and nsp16, and between the cofactor SAM and nsp16 in its both monomer and the complex were computed using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) methodology (Kollman et al., 2000; Wang et al., 2006; Kar et al., 2007a,b, 2011, 2013; Hou et al., 2011). MMPBSA.py script available in the AmberTools19 was used for the analysis. Details of the MM-PBSA protocol were provided in our previous studies (Kar and Knecht, 2012a,b,c,d; Kar et al., 2013; Jonniya et al., 2019; Jonniya and Kar, 2020; Roy et al., 2020; Sk et al., 2020a,d), and the same protocol was adopted here.

The binding free energy is estimated by using the following equations of the MM-PBSA scheme:

$\begin{array}{ll}\Delta G_{\text{bind}}=\Delta H-T\Delta S\approx \Delta E_{\text{internal}}+\Delta G_{\text{solv}}-T\Delta S& (3)\end{array}$

$\begin{array}{ll}\Delta E_{\text{internal}}=\Delta E_{\text{covalent}}+\Delta E_{\text{elec}}+\Delta E_{\text{vdW}}& (4)\end{array}$

$\begin{array}{ll}\Delta G_{\text{solv}}=\Delta G_{\text{pol}}+\Delta G_{\text{np}}& (5)\end{array}$

where ΔE_internal, ΔG_solv, and TΔS represent the total internal energy, desolvation free energy, and conformational entropy, respectively. Further, the internal energy is composed of ΔE_covalent (bond, dihedral, and angle), ΔE_elec (electrostatic) and ΔE_vdW (van der Waals) and the desolvation free energy is composed of polar (ΔG_pol) and non-polar solvation energies (ΔG_np). Here, ΔG_pol was estimated from the Poisson-Boltzmann (PB) equation. The dielectric constant of the solute and solvent was set to 1.0 and 80.0, respectively. ΔG_np was estimated from the following Equation (6),

$\begin{array}{ll}\Delta G_{np}=\gamma (SASA)+b& (6)\end{array}$

Here, SASA represented the solvent-accessible surface area and was estimated using the LCPO (linear combination of pairwise overlap) algorithm (Weiser et al., 1999) with a probe radius of 1.4 Å. The surface tension coefficient, γ and offset (b) value were set to 0.00542 kcal.mol⁻¹.Å⁻² and 0.92 kcal.mol⁻¹, respectively (Gohlke et al., 2003). For the MM-PBSA calculation, 2000 snapshots selected from the last 300 ns trajectory with a frequency of 15 ps were employed.

The configurational entropy was calculated using the normal mode analysis (NMA) method (Karplus and Kushick, 1981; Rempe and Jónsson, 1998; Xu et al., 2011), and ~300 snapshots were used in the calculation due to high computational cost. Each configuration was energy minimized in an implicit solvent (nmode_igb = 1) using a maximum of 50,000 steps, and a target root-mean-square gradient of 10⁻⁴ kcal mol⁻¹Å⁻¹ via mmpbsa_py_nabnmode (Hawkins et al., 1995, 1996; Kar and Knecht, 2012a,b,c; Kar et al., 2013). NMA was applied to calculate the vibrational entropy (S_vib) using the following equation (Carlsson and Åqvist, 2006);

$\begin{array}{ll}{S}_{vib}=R[\frac{x}{(e^x-1)}-ln(1-e^{-x})];x=\frac{h\upsilon }{kT}& (7)\end{array}$

Where R denotes the universal gas constant, k is the Boltzmann constant, and T is the absolute temperature. The Planck constant and the vibrational frequency are denoted as h and ν, respectively.

Further, the decomposition of the total binding free energy at the residue level was conducted by using the Molecular Mechanics Generalized Born Surface Area (MM-GBSA) scheme (Kar and Knecht, 2012a,b,c; Kar et al., 2013; Jonniya and Kar, 2020) as per the following equation.

$\begin{array}{ll}\Delta G_{\text{residue}}=\Delta E_{\text{vdW}}+\Delta E_{\text{ele}}+\Delta v{G}_{\text{pol}}+\Delta G_{\text{np}}& (8)\end{array}$

The total contribution from each residue can also be defined as the sum of van der Waals (ΔE_vdW), electrostatic (ΔE_ele), polar (ΔG_pol), and non-polar (ΔG_np) terms. This method was proposed by Gohlke et al. (2003).

Finally, the Computational Alanine Scanning (CAS) was performed for some essential residues. This method yields the energy difference between the wild type and mutant (alanine) variants.

$\begin{array}{ll}\Delta \Delta G_{\text{bind}}=\Delta G_{\text{mutant}}-\Delta G_{\text{wild}}& (9)\end{array}$

The basic principle involves in AS (alanine scanning) is the substitution of a residue with alanine that has an impact on the side-chain beyond C_β and not in the main-chain. In vitro, AS has been proven an advantageous mutagenesis method in finding hotspot residues in protein-protein interfaces. CAS is an excellent alternative approach to the in vitro experimental alanine scanning (Massova and Kollman, 1999; Moreira et al., 2007). Therefore, in this study, CAS was applied to illustrate further the importance of specific residues except alanine mentioned in the binding decomposition free energy of nsp16 of COVID-19.

Results and Discussion

To explore the mechanism underlying the dimerization of nsp16 and nsp10 as well as the preferential binding of the cofactor SAM to the nsp16/nsp10 heterodimer compared to nsp16 (nsp16_SAM), a conformational free energy landscape (FEL) and binding free energy calculations were performed using the MD/MM-PBSA scheme.

Overall Structural Dynamic Features

To explore the thermodynamic stability of each system and to ensure the rationality of the sampling method, we monitored the structural and energetic properties during the entire 1 μs production simulation. The time evolution of the root-mean-square deviations (RMSDs) of the protein backbone atoms, which reflects the stability of the system, were calculated relative to the initial configuration and shown in Figure 2A. It is worth mentioning here that a stable RMSD does not always provide stable energy profiles; hence we have also verified the potential and total energy of each system from the respective MD trajectory (data not shown).

FIGURE 2

Figure 2. (A) The root-mean-square deviations (RMSDs) of the backbone atoms relative to their initial coordinates as a function of simulation time, (B) the root-mean-square fluctuations (RMSFs) of C_α atoms for each residue in the complex and monomers of nsp16_SAM and nsp10_apo.

The RMSD plots indicated that all the studied systems reached equilibrium after ~100 ns. The average RMSD value varies between 2.8 and 4.0 Å for all systems (see Table 1). The highest and lowest deviations were obtained for nsp10_apo and complex, respectively. An intermediate RMSD value of 3.1 Å was obtained for nsp16_SAM. Overall, the nsp16/nsp10/SAM complex showed more stability throughout the simulations than nsp16_SAM or nsp10_apo. This suggests that the dimerization leads to the overall stabilization of the complex.

TABLE 1

Table 1. The average value of protein backbone RMSD and solvent-accessible surface area (SASA) obtained from MD simulations.

Moreover, we also calculated the temporal RMSDs of heavy atoms of the cofactor SAM (see Supplementary Figure 1A) and the backbone atoms of residues within 5 Å around SAM in the binding pocket (see Supplementary Figure 1B). It is evident from Supplementary Figure 1B that the RMSD values of SAM and the binding pocket for the complex were stable up to 850 ns of the simulation. After that, in the last 150 ns, a sudden increase in RMSD was observed. Overall, the average RMSD value was 0.7 Å and 0.8 Å for SAM and the binding pocket, respectively (see Supplementary Table 1). It indicates that SAM binds strongly onto the cofactor binding cavity of the heterodimer nsp16/nsp10. However, in the case of nsp16_SAM, the RMSD values of SAM and the binding cavity were relatively stable up to 100 ns. After that, the RMSD value increased by ~3 times compared to the complex system (see Supplementary Table 1 and Supplementary Figure 1). These high values of RMSDs suggest that SAM does not bind to nsp16 monomer (nsp16_SAM), which is in agreement with the experimental study by Chen and coworkers for SARS-CoV (Chen et al., 2011).

Additionally, three loops, namely 71–79, 100–108, and 130–148, which comprises the SAM binding pocket as shown in the crystal structure, were also analyzed and found to be stable for both complex and nsp16_SAM (see Supplementary Figures 2A–C). However, the cap-binding groove of nsp16 in SARS-CoV, as well as SARS-CoV-2, is mainly composed of two flexible loops, comprised of residues 26–38 and 130–148, while in the case of flavivirus, the NS5 MTase was relatively stable with α-helices (A1, A2, and half of αD) along the cap-binding groove (Egloff et al., 2002). Among the two loops involved in the cap-binding, only the loop 26–38 exhibited differences in the complex and nsp16_SAM system, as revealed from our simulations. The 26–38 loop is located near the interface of nsp16 and nsp10. On the other hand, the 130–148 loop is found near the SAM binding pocket, away from the nsp16/nsp10 interface. As shown in Supplementary Table 1, the average RMSD value of this loop was higher in nsp16_SAM (3.1 Å) compared to the complex (1.4 Å). The time evolution of RMSD of the 26–38 loop and the corresponding potential mean force (PMF) were displayed in Supplementary Figures 3A,B. From Supplementary Figure 3A, it was observed that in the case of nsp16_SAM, the RMSD value of the loop 26–38 increased during the initial 300 ns and remained stable thereafter. On the other hand, in the complex case, the RMSD of the loop was found to increase during the initial 300 ns and gradually decreased up to 600 ns of the simulation and finally reached equilibrium. We computed the potential of mean force (PMF), taking the loop's RMSD as a reaction coordinate and shown in Supplementary Figure 3B. The primary low energy structure of the loop was observed at ~3.2 Å for nsp16_SAM. However, for the complex system, the global energy minimum was found at a relatively lower RMSD (1.4 Å). Besides, we detected two secondary minima (~1.9 Å, and ~2.8 Å) for the 26–38 loop in the complex system, and the energy barriers of the adjacent states were ~1.1 and 1.2 kcal/mol, respectively. Overall, our results suggest that the dimerization resulted in the loop rearrangement near the interface.

Next, we measured the root-mean-square fluctuations (RMSFs) of Cα atoms for all systems (see Figure 2B) to elucidate the effect of dimerization on the residual flexibility. It is evident from Figure 2B that all three systems displayed a similar trend in RMSF patterns. However, we found some dynamic fluctuations in different loop regions, including the N- and C-terminals. A relatively large fluctuation was observed for various loop regions located at residues 26–38, 74–80, 104–115, and 130–148. The loop regions 26–38 and 104–115, which contribute to the interface of the heterodimer, showed higher fluctuations in complex compared to nsp16SAM, which suggests that the dimerization leads to the loop rearrangements. In contrast, the 74–80 loop region, located near SAM binding pocket, stabilized after the heterodimerization. However, the cap-binding groove of the 130–148 loop region exhibited no differences in both complex and apo forms, which agrees with the RMSD analysis (Supplementary Figure 2C). Overall, it depicts the binding of nsp10 to nsp16 favors and stabilizes the binding pocket of SAM, leading to its high activity.

Finally, we also measured the solvent-accessible surface area (SASA) of all systems and reported in Table 1. The average SASA values were found to be 18928.1, 13467.3, and 6663.8 Ų, for the complex, nsp16_SAM, and nsp10_apo, respectively. This suggested that the dimerization of nsp16 and nsp10 covered ~1,203 Ų surface area in total, indicating a very stable interaction. Our simulation results also favor an experimental study that showed the value of the solvent-exposed surface area for the heterodimer complex of nsp16/nsp10 as 19710 Ų (Rosas-Lemus et al., 2020).

Free Energy Landscape (FEL) of SAM Binding

The RMSD profile of SAM and the binding pocket (Supplementary Figure 1) for the complex and nsp16_SAM systems suggested that SAM binds more strongly to the heterodimer nsp16/nsp10 than the monomer nsp16. The 2-dimensional free energy landscape (FEL) was generated for the complex and nsp16_SAM systems to explore the strong and weakly bound states of the SAM molecule. The RMSD of heavy atoms of SAM and the center of mass (CoM) distance between SAM and nsp16, which characterizes the displacement of SAM from the substrate-binding pocket of nsp16, were used as reaction coordinates for the construction of FEL. The value of a CoM distance <14.5 Å denotes the strong binding of SAM toward nsp16. The two-dimensional FEL of the complex (nsp16/nsp10/SAM) and monomer (nsp16_SAM) systems were shown in Figures 3A,B, and the corresponding positions of SAM in the binding cavity were depicted in Figures 3C,D. FEL of the complex and monomer systems suggested that the conformational state of SAM in two systems prefer disparate configurations. In the case of the heterodimeric complex system, the underlying FEL was characterized by a single global minimum (see Figure 3A), which corresponds to the strongly bound state of SAM. On the other hand, in the case of nsp16_SAM, the global free energy minimum was obtained at a CoM distance of ~17.5 Å, which corresponds to a weakly or unbound state of SAM. The secondary minimum was detected at a CoM distance of ~13.5 Å, which corresponds to a strongly bound state of SAM to nsp16_SAM. Overall, our results suggested that although nsp16 monomer (nsp16_SAM) favored two binding states, the weakly bound or unbound state is energetically more favorable than the strongly bound state. These results illustrated a stronger binding of SAM with the nsp16/nsp10 heterodimer than nsp16 alone (nsp16_SAM), which agrees with the experimental results obtained for SARS-CoV (Chen et al., 2011).

FIGURE 3

Figure 3. Free energy landscape (FEL) and representative SAM position in the complex (nsp16/nsp10/SAM) and nsp16_SAM (nsp16/SAM). (A) FEL of SAM bound with nsp16 in the complex, (B) FEL of SAM bound with monomer nsp16_SAM. Representative structures of (C) SAM with nsp16 in the complex and (D) SAM with monomer nsp16_SAM.

Principal Component Analysis of nsp16/nsp10

The principal component analysis (PCA) was carried out for both nsp16 and nsp10 when they are in complex and monomeric states. FELs of sub-units nsp16 (nsp16_complex) and nsp10 (nsp10_complex) in the complex were compared with the corresponding monomer simulations (nsp16_SAM and nsp10_apo). Each eigenvector and eigenvalue was plotted in the decreasing order, as shown in Supplementary Figure 4. For all cases, the first few eigenvectors describe the collective motion of local fluctuations. Comparing the four systems indicated that the first few PCs that describe the properties of movements were not the same. The first two eigenvectors encapsulated for 42, 67, 46, and 70% of overall movements in nsp16_complex (nsp16 in complex), nsp10_complex (nsp10 in complex), nsp16_SAM, and nsp10_apo, respectively. Similarly, the first ten eigenvectors accounted for 80–90% of the total motions in all four systems.

Next, we constructed the 2-dimensional FEL at 300 K for each system using the first two principal components (PC1 and PC2) as reaction coordinates and shown in Figures 4A–D. The sampling space of four systems was different, as evident from Figure 4. The conformational sampling space of nsp16_complex (nsp16 in nsp16/nsp10/SAM), depicted in Figure 4A, was found to be restricted compared to the other three systems, which sampled more expansive conformational space. From Figure 4A, we observed a global minimum of 62.7% occupancy and a secondary minimum of 37.3% occupancy, which suggested that the nsp16_complex system is more stabilized in the presence of nsp10. In the monomeric form of nsp16_SAM, a wider conformational basin was sampled (see Figure 4C), and the energy barriers between adjacent minima were ~1.5 kcal/mol. On the other hand, nsp10 sampled more phase space than nsp16 in the complex as well in its monomer apo form, as shown in Figures 4B,D. However, nsp10_apo has three distinct conformations with an energy barrier of 4.0 kcal/mol suggesting structural disparity in both nsp10_complex and nsp10_apo as compared to nsp16.

FIGURE 4

Figure 4. Two-dimensional FEL generated by projecting the first two principal components, PC1 and PC2 for (A) nsp16_complex (nsp16 in nsp16/nsp10/SAM), (B) nsp10_complex (nsp10 in nsp16/nsp10/SAM), (C) nsp16_SAM (nsp16/SAM), and (D) nsp10_apo. The representative structures are shown on the left panel.

Residual Network Analysis

The Network Analysis of Protein Structure (NAPS) server provided a visual examination of sub-network based on the physicochemical properties of the protein residues to find more details about the interaction between nsp16 and nsp10 (see Figure 5). Herein, we explored the 3D interaction network of hydrophobic residues of nsp16 and nsp10. The results showed that the number of hydrophobic interaction networks between nsp16 and nsp10 was initially very high, and after 100 ns, the number of networks reduced to ~2–3. The hydrophobic interaction networks, such as (V104, A71) and (P80, V42), were found as strong and stable contacts throughout the simulation time (see Supplementary Table 2). We also calculated the total number of interactions and noticed that the number of interaction networks was more or less conserved throughout the first 500 ns (see Supplementary Figure 5). In general, it can be concluded from our analyses that nsp16 and nsp10 interact with a very high affinity.

FIGURE 5

Figure 5. Sub-network representation of network 3D view of (A,B) initial structure and (C,D) final structure of the complex. Each residual C_α atom is represented by a sphere where the red sphere indicates hydrophobic residues. Blue and green color spheres correspond to residues of nsp16 and nsp10, respectively. The yellow line edges represent the contact network of hydrophobic residues between nsp16 and nsp10.

Binding Free Energy of SAM Bound to nsp16/nsp10 Heterodimer and nsp16 Alone

The differences in binding affinity of the SAM molecule toward nsp16 alone (nsp16_SAM) and nsp16/nsp10 heterodimer (complex) were calculated by utilizing the MM-PBSA scheme. The MM-PBSA scheme provides various components contributing to the total binding energy (ΔG_bind), such as van der Waals interactions (ΔE_vdW), electrostatic interactions (ΔE_ele), polar solvation energy (ΔG_pol), non-polar solvation free energy (ΔG_np), and configurational entropy (TΔS). All these components of the binding free energy of SAM to nsp16_SAM or nsp16/nsp10 (complex) were shown in Table 2. The binding free energy of SAM to the heterodimer nsp16/nsp10 was estimated as −6.8 (± 0.9) kcal/mol, which agrees with the experimental result (−7.4 ± 0.1 kcal/mol) (Lin S. et al., 2020). On the other hand, SAM was found to bind very weakly (ΔG_bind = −0.6 ± 0.8 kcal/mol) to nsp16 monomer (nsp16_SAM), as evident from Table 2. This is consistent with Figures 3A,B. A similar mode of SAM binding was observed for SARS-CoV nsp16/nsp10 or nsp16 alone. For both cases, the intermolecular van der Waals interactions (ΔE_vdW), electrostatic interactions (ΔE_ele), and non-polar solvation free energy (ΔG_np) favored the binding of SAM. In contrast, polar solvation energy (ΔG_pol) and entropy (TΔS) disfavored the complexation. The electrostatic interaction energy was more favorable than the van der Waals interactions for both cases. In nsp16_SAM, although ΔE_ele was much more favorable (−76.8 kcal/mol) than ΔE_vdW (−28.6 kcal/mol). However, the disfavorable polar solvation energy (ΔG_pol = 84.7 kcal/mol) overcompensated the intermolecular electrostatic interaction energy. Hence, the overall polar contribution (ΔE_ele + ΔG_pol) was found to be unfavorable (7.9 kcal/mol) to the SAM binding. In contrast, the overall non-polar contribution (ΔE_vdW + ΔG_np) was found to be −32.1 kcal/mol. It implies that hydrophobic interactions drive the binding. In the case of the complex system (nsp16/nsp10/SAM), the binding affinity of SAM to nsp16 increases, as evident from Table 2. Here, the electrostatic interactions energy (ΔE_ele = −181.1 kcal/mol) was more favorable than the van der Waals energy (ΔE_vdW = −40.4 kcal/mol). The total polar contributions (ΔE_ele + ΔG_pol) value was found to be favorable (−2.1 kcal/mol) to the binding of SAM, which is in contrast to what had been observed for nsp16_SAM. Further, the overall polar energy was less favorable than the net non-polar contribution (ΔE_vdW + ΔG_np = −44.5 kcal/mol). Hence, in the complex, the main force behind the binding of SAM to nsp16 is hydrophobic interactions. Overall, the binding free energy analysis showed that the binding affinity of the SAM molecule to nsp16 increases with the presence of nsp10. This implies that nsp10 acts as a stimulator to bind SAM to nsp16/nsp10 of SARS-CoV-2, which agrees with the experiment (Viswanathan et al., 2020). Wang and coworkers obtained a similar result for SARS-CoV and MERS-CoV (Wang et al., 2015).

TABLE 2

Table 2. Energetic components of the total binding free energy of SAM bound to the heterodimer nsp16/nsp10 (complex), and nsp16 alone (nsp16_SAM).

Further, to explore the significant residues of nsp16 in the binding of SAM and to evaluate the differences in the complex (nsp16/nsp10) and monomer (nsp16_SAM), the per-residue decomposition of the binding free energy was performed using the MM-GBSA approach. All the residues having > 1.5 kcal/mol of energetic contributions were considered important and shown in Figure 6. As seen in Table 3, the number of amino acids contributing to binding with SAM is high in the complex as compared to nsp16_SAM. Residues such as Leu100, Asp99, Cys115, Met131, Phe149, and Asp114 are common in cases of complex and nsp16_SAM. Apart from these residues, additional residues, such as Asn43, Tyr47, Gly71, Tyr132, and Ser74, also played a significant role in the binding of SAM to the complex heterodimer. These results are consistent with the higher binding affinity of SAM in the complex than the nsp16 monomer. All these residues were also considered as SAM-engaging in an experimental study by Lin S. et al. (2020). It is worth noting here that all these SAM-interacting residues are conserved both in SARS-CoV and MERS-CoV, suggesting a conserved binding mode of SAM in these viruses. The conserved SAM-binding mechanism to nsp16/nsp10 also opens the possibility of developing a pan-CoV inhibitor by targeting this SAM-binding pocket.

FIGURE 6

Figure 6. Per-residue decomposition free energy of nsp16 for the binding of SAM (in the presence of nsp10 and without nsp10) and respective binding pocket drawn from MD snapshots. (A,B) nsp16 with SAM in the presence of nsp10, (C,D) nsp16 with SAM without nsp10.

TABLE 3

Table 3. Per-residue decomposition of the total binding free energy (kcal/mol) of SAM to the complex (nsp16/nsp10) and monomer (nsp16_SAM) systems.

Hydrogen Bond Interaction Between nsp16 and SAM

Next, we investigated the hydrogen bond (H-bond) interactions between SAM and nsp16 in the complex and apo forms (see Table 4). The % occupancy calculated from the respective MD trajectories reflected the stability of H-bonds. As seen in Table 4, residues with > 50% H-bond occupancy in the MD simulations were recognized for the complex (nsp16/nsp10) compared to monomer nsp16_SAM. The residues, including Asp130, Tyr47, Gly71, and Asp99, showed more H-bond stability in the complex. Overall, a stable H-bond is formed between nitrogen and oxygen atoms of SAM with Asp130, Tyr47, Gly71, and Asp99 residues of nsp16, respectively. Hence, the identified significant residues like Asp130, Tyr47, Gly71, and Asp99 of nsp16 may aid in the development of SAM competitive inhibitors.

TABLE 4

Table 4. Hydrogen bonds formed between SAM and receptor (nsp16/nsp10 or nsp16).

The changes in the distance of the atoms forming H-bond between SAM and nsp16 were also determined and shown in Supplementary Figure 6. In the complex, the distance between the oxygen atom of Asp130, Tyr47, and Gly71 of nsp16 and the nitrogen atom of SAM illustrated the average distance of 3.01, 2.89, and 3.01 Å, respectively, suggesting strong H-bond interactions. In the case of nsp16_SAM, the average distances for these atoms were 12.48, 14.64, and 11.06 Å, respectively, indicating a lack of formation of H-bonds in the nsp16 monomer. However, for Asp99, the distance between its oxygen atom (OD1 and OD2) and the oxygen atom (O3) of SAM illustrated an average distance of 2.9 Å and 3.2 Å for the complex and monomer, respectively. In contrast, the average distance between Asp99 (OD2) and SAM (O2) was higher for monomer (4.2 Å) than the complex (3.1Å). All these results were consistent with the occupancy analysis. Further, they emphasized that Asp130, Tyr47, and Gly71 of nsp16 were key residues in the binding of SAM and for the 2′-O-MTase activity of nsp16.

Energetics of nsp16/nsp10 Complexation

To evaluate further the binding free energy of the heterodimer (nsp16/nsp10), ΔG_bind was estimated using the MD/MM-PBSA approach and reported in Table 5. From Table 5, the binding free energy (ΔG_bind) between nsp16 and nsp10 was found to be −47.4 kcal/mol. The various components of the total binding free energy indicate that the intermolecular van der Waals interactions (ΔE_vdW), electrostatic interactions (ΔE_ele), and non-polar solvation free energy (ΔG_np) favors the binding of nsp10 and nsp16. While polar solvation free energy (ΔG_pol) disfavor the binding. The electrostatic interactions (ΔE_ele) energy is higher (−429.4 kcal/mol) than the van der Waal interactions (ΔE_vdW) (−90.4 kcal/mol) (see Supplementary Figure 7). However, the disfavouring components of polar solvation energy (ΔG_pol) compensate for the ΔE_ele being a value of 481.7 kcal/mol. Hence, the total polar interactions (ΔE_ele + ΔG_pol) disfavor the binding between nsp10 and nsp16 with a value of 52.3 kcal/mol. In contrast, the total non-polar contributions (ΔE_vdW + ΔG_np) favor the binding between nsp16 and nsp10 (−100.0 kcal/mol). Therefore, the binding between nsp10 and nsp16 is mainly driven by hydrophobic interactions. Overall, the binding free energy analysis depicts, the binding affinity of the SAM in the complex nsp16/nsp10 (ΔG_bind = −46.6 kcal/mol) is similar to that between nsp10 and nsp16 (ΔG_bind = −47.4 kcal/mol) of SARS-CoV-2. These observations were in agreement with its similar homology virus, SARS-CoV (2002). A similar binding affinity was seen between nsp10 and nsp16, as well as between SAM and the nsp16/nsp10 complex of SARS-CoV (Chen et al., 2011).

TABLE 5

Table 5. Energetic components of the dimerization energy (G_bind) between nsp16 and nsp10.

To further investigate the critical residues involved in the binding between nsp10 and nsp16, the binding free energy was decomposed at the individual residue level using the MM-GBSA approach and shown in Figure 7. Here, we considered only residues having the energy value ≥ 1.0 kcal/mol and listed in Supplementary Table 3. The key residues from nsp10 include Leu45, Ala71, Val42, Met44, Tyr96, Gly69, Thr47, Arg78, Gly70, Gly94, and Pro59. Similarly, the critical residues from nsp10 include Ile40, Val104, Ala83, Val78, Met247, Val44, Gln87, Arg86, Lys76, Lys38, Val84, and Met41. It indicates that hydrophobic interactions significantly control the binding between nsp10 and nsp16. Our findings agree with a recent experimental study (Lin S. et al., 2020) where it was reported that a total of 31 residues (Asn40, Val42, Lys43, Met44, Leu45, Cys46, Thr47, Pro59, Gly70, Ala71, Cys77, Arg78, Lys93, Gly94, and Tyr96 in nsp10 and Lys38, Gly39, Ile40, Met41, Val44, Lys76, Val78, Pro80, Ala83, Arg86, Gln87, Val104, Ser105, Asp106, Leu244, and Met247 in nsp16) mediates the intimate nsp16-nsp10 interaction. The per-residue decomposition of the binding free energy was further validated by other computational approaches, such as FTmap (Kozakov et al., 2015). The final complex structure obtained from the MD trajectory was analyzed using FTmap, and Gln87, Arg86, Val104, Met247 of nsp16, and Leu45, Met44, Pro59, Arg78, Tyr96 of nsp10 were identified as the hotspot residues, which again agrees with the MM-GBSA analysis.

FIGURE 7

Figure 7. (A) Decomposition of ΔG on a per-residue basis for the complex, (B) corresponding residual position of nsp16 in complex, and (C) residual position of nsp10 in the complex.

Interaction Analysis Between nsp16 and nsp10

In the above finding of binding free energy analysis, we have seen that nsp10 acts as a stimulator in the binding of the SAM to nsp16. Hence, to further explore the binding interactions between nsp16 and nsp10, H-bond and hydrophobicity analysis was calculated. As seen in Table 6, the H-bond occupancy reflects the stability of H-bond formation between protein-protein in the MD simulations. The strong H-bonds were formed between residues Ala83 (nsp16) and Tyr96 (nsp10) and between Gln87 (nsp16) and Leu45 (nsp10). Other residues, including Asp106 (nsp16), formed two H-bonds with Ala71 and Gly94 of nsp10, Lys38 (nsp16) to Lys43 (nsp10), and Ser105 (nsp16) to Lys93 (nsp10). Hydrophobic interactions also played an important role in protein-protein or protein-ligand interactions (Jonniya et al., 2020; Roy et al., 2020; Sk et al., 2020d). Different H-bonds and hydrophobic interactions from the stable structure of the nsp16-nsp10 obtained from the MD simulations were plotted via Ligplot (Wallace et al., 1995) and shown in Figure 8.

TABLE 6

Table 6. The hydrogen bonds formed between nsp16 and nsp10^* in the complex and the corresponding average distance and percent determined using the production trajectories in the MD simulations.

FIGURE 8

Figure 8. Protein-protein interaction diagrams for the nsp16/nsp10 dimerization. The upper panel corresponds to nsp10, and the lower panel corresponds to nsp16 in the complex. The plot was generated with the dimplot module of LigPlot+. Hydrogen bonds are shown as a lime green dotted line, and hydrophobic bonds are represented in red.

The percentage occupancy of the residual contacts with the cut-off 3.9 Å for the protein-protein is also listed in Supplementary Table 4. Overall, these results highlighted the significant residues forming strong interactions between nsp16 and nsp10, which may aid in the development and design of inhibitors that could block these protein-protein interactions by inhibiting the 2′-O-MTase activity.

Computational Alanine Scanning

The changes in the binding free energy were computed as in Equation (4) (see Method section) after replacing the residue of WT to alanine. The impact of a single mutation in the protein-protein interaction is mainly reflected by the more positive value of the ΔΔG_bind. In our study, we have conducted CAS for nsp16 by considering those residues with the decomposition free energy > 1.5 kcal/mol, as given in Supplementary Table 3. The binding free energy components calculated from the CAS mutagenesis for residues I40A, V104A, R86A, V78A, V44A, M247A, and Q87A of nsp16 are listed in Table 5 and compared with the WT. As seen in Supplementary Table 3, although the ΔG_bind of WT is higher than mutants, different energy components of mutants follow the same trend as WT. The ΔE_ele is higher than ΔE_vdW, but disfavouring polar solvation energy ΔG_pol compensates ΔE_ele. As shown in Figure 9A, for all the systems, the binding free energy is mainly coming from the hydrophobic interactions, where the total non-polar energy (ΔE_vdW + ΔG_np) is higher than the total polar energy (ΔE_ele + ΔG_pol). Figure 9B reflects significant residues in the protein interface calculated from CAS mutagenesis. It shows that ΔΔG_bind for mutants, I40A, V104A, and R86A are comparatively high, suggesting that primarily these residues play a significant role in the heterodimer formation, which is in accordance with the decomposition of energy. The CAS mutagenesis results depict that in addition to hydrophobic residues I40 and V104, the hydrophilic residue R86 also plays a vital role in the binding of nsp10-nsp16.

FIGURE 9

Figure 9. (A) Binding free energy components of the wild type and seven mutations (alanine scanning) in the complex, (B) alanine scanning mutagenesis analysis of complex.

Temperature and Salt Concentration Effect on the Binding Affinity

Next, we investigated the effect of temperature and salt concentration in the binding of nsp16/nsp10 and nsp16/nsp10/SAM using the MM-PBSA method. This method presents itself to be a perfect model to study the effect of conformational dynamics of the same protein evolved in different environmental conditions to the fate of the ligand, ultimately decided by the binding affinity or dissociation.

Supplementary Tables 5, 6 displayed the binding free energy components for nsp16-nsp10 and nsp16/nsp10/SAM, respectively, at two different temperatures. As shown in Supplementary Table 5, the predicted binding free energies between nsp10 and nsp16 (ΔG_bind) was found to be −47.7 kcal/mol, and −52.2 kcal/mol at 300 and 310 K, respectively. It suggests that nsp16 binds more strongly to nsp10 at 310 K than 300 K. It can also be noted from Supplementary Table 5 that at the elevated temperature (310 K), the magnitude of total polar interactions (ΔE_ele + ΔG_pol) decreases (52.3 to 49.6 kcal/mol) and oppose less unfavorably to the protein-protein association. At a higher temperature, the non-polar interactions (ΔE_vdW + ΔG_np) decreased from −100 to −102 kcal/mol and contributed more favorably to the association of nsp16-nsp10. Similarly, Supplementary Table 6 suggested that the estimated binding free energies (without entropy) of the cofactor SAM with nsp16/nsp10 (ΔG_bind) were −46.6 and −13.2 kcal/mol at 300 and 310 K, respectively. It was further evident from Supplementary Table 6 that after increasing the temperature, the binding affinity of SAM toward nsp16/nsp10 decreased. At room temperature (300 K) both total polar (ΔE_ele + ΔG_pol = −2.1 kcal/mol) and non-polar (ΔE_vdW + ΔG_np = −44.5 kcal/mol) interactions favored the binding of SAM with nsp16. After increasing the temperature to 310 K, the total polar interaction (ΔE_ele + ΔG_pol) became unfavorable to the binding of SAM to nsp16/nsp10. However, the total non-polar contribution (ΔE_vdW + ΔG_np) still favored the binding. As compared to room temperature, the strength of the non-polar (ΔE_vdW + ΔG_np) interactions decreases significantly (−44.5 to −23.5 kcal/mol). At 310 K, we have noticed that up to 400 ns, SAM binds strongly to nsp16. After that, SAM got detached from the binding cavity (see Supplementary Figure 8).

The salt concentration also influenced the protein-protein interactions as well as the binding of SAM to nsp16/nsp10, which are evident from Supplementary Tables 5, 6. As shown in Supplementary Table 5, the predicted binding free energies of nsp10 with nsp16 (ΔG_bind) were found to be −44.8 and −37.7 kcal/mol, for the salt concentration of 0.15 and 0.25 M, respectively. After increasing the salt concentration, it was observed that the total polar (ΔE_ele + ΔG_pol) and non-polar interactions (ΔE_vdW + ΔG_np) increased from 54.1 to 55.3 kcal/mol and −98.9 to −93 kcal/mol, respectively. The rise in unfavorable interactions caused the destabilization of the protein-protein interactions. Similarly, the binding of SAM to the heterodimer was found to be affected by the increased salt concentration (see Supplementary Table 6). The estimated binding free energy of SAM (ΔG_bind) was −23.7 and −22.1 kcal/mol for the salt concentration of 0.15 and 0.25 M, respectively. These values are significantly higher than what was obtained for the neutral salt concentration (−46.6 kcal/mol). ΔE_vdW and ΔE_elec were found to contribute less favorably to the binding of SAM to nsp16/nsp10 in comparison to neutral simulations. A similar trend was observed when the temperature was increased from 300 to 310 K. The time evolution of CoM distance between nsp16 and SAM was displayed in Supplementary Figure 9. It is evident from Supplementary Figure 9 that with the increased salt concentration, the CoM distance between nsp16 and SAM increased and remained stable around 17 Å and 18 Å for 0.15 and 0.25 M, respectively. This implies that the cofactor SAM moves away from the binding cavity with the increased salt concentration.

Conclusions

Herein, we have employed extensive MD simulations of 1 μs along with the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) method to study the binding mechanism of the cofactor SAM to the nsp16/nsp10 heterodimer and sub-unit nsp16 alone (nsp16_SAM) of SARS-CoV-2. Our MD/MMPBSA calculations suggested that SAM binds strongly to the nsp16/nsp10 heterodimer and fails to bind to the nsp16 monomer. It emphasized that nsp10 helps in the strong interaction between SAM and nsp16 to execute the 2′-O-MTase activity. This finding agrees well with the experimental study, and a similar observation was made for other coronaviruses, including SARS-CoV and MERS-CoV. We have also investigated the interactions between nsp16 and nsp10. Our study revealed that the binding of nsp10 stabilizes the complex (nsp16/nsp10) structure. Further, our study showed that hydrophobic interactions are critical for the heterodimer association. Thus, apart from the active site of nsp16, the interface of nsp16/nsp10 can also be considered as a potential target site in the design of antiviral drugs such as peptide inhibitors. It includes Ile40, Val104, Ala83, Val78, Met247, Val44, Gln87, Arg86, Lys76, Lys38, Val84, and Met41 from nsp16, and Leu45, Ala71, Val42, Met44, Tyr96, Gly69, Thr47, Arg78, Gly70, Gly94, and Pro59 from nsp10. Besides, the stable hydrogen bond between Ala83 (nsp16) and Tyr96 (nsp10), and between Gln87 (nsp16) and Leu45 (nsp10) were important in the nsp16-nsp10 interface. The CAS study reveals that residues I40A, V104A, R86A, V78A, V44A, M247A, and Q87A of nsp16 were considered as hot spot residues for the association of nsp16-nsp10. ΔΔG_bind for mutants, I40A, V104A, and R86A are comparatively high, suggesting the significance of these residues. Finally, we investigated the effect of temperature and salt concentration on the binding of SAM to the heterodimer. Our study revealed that the SAM binding was impaired due to an increase in temperature or salt concentration. Finally, our study provides a comprehensive understanding of the dynamic and thermodynamic process of binding nsp16 and nsp10 that will contribute to the rational design of inhibitors targeting nsp16/nsp10.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

Project conceived and supervised by PK. Simulations and analyses were performed by MS, NJ, RR, and SP. Manuscript was written by MS, NJ, and PK. All authors contributed to the article and approved the submitted version.

Funding

This work was partially supported by the Department of Biotechnology, Govt. of India (Grant number BT/RLF/RE-entry/40/2014, DBT Ramalingaswami Re-entry Fellowship), and Department of Science and Technology (DST), Govt. of India (Grant number ECR/2017/000010). MS would like to thank DST, Govt. of India, for providing fellowship under the INSPIRE Fellowship Scheme (DST/INSPIRE Fellowship/2017/IF170145).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This manuscript has been released as a pre-print at ChemRxiv (Sk et al., 2020b).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2020.590165/full#supplementary-material

References

Amadei, A., Linssen, A., de Groot, B., van Aalten, D., and Berendsen, H. (1996). An efficient method for sampling the essential subspace of proteins. J. Biomol. Struct. Dyn. 13, 615–625. doi: 10.1080/07391102.1996.10508874

PubMed Abstract | CrossRef Full Text | Google Scholar

Arabi, Y., Harthi, A., Hussein, J., Bouchama, A., Johani, S., Hajeer, A., et al. (2015). Severe neurologic syndrome associated with Middle East respiratory syndrome corona virus (MERS-CoV). Infection 43, 495–501. doi: 10.1007/s15010-015-0720-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Berendsen, H. J., Postma, J. V., van Gunsteren, W. F., Dinola, A., and Haak, J. (1984). Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690. doi: 10.1063/1.448118

CrossRef Full Text | Google Scholar

Bogoch, I. I., Watts, A., Thomas-Bachli, A., Huber, C., Kraemer, M. U., and Khan, K. (2020). Potential for global spread of a novel coronavirus from China. J Travel Med. 27:taaa011. doi: 10.1093/jtm/taaa011

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlsson, J., and Åqvist, J. (2006). Calculations of solute and solvent entropies from molecular dynamics simulations. Phys Chem. Chem. Phys. 8, 5385–5395. doi: 10.1039/B608486A

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Huang, Y., Ben-Shalom, S. R., Brozell, D. S., Cerutti, T. E., Cheatham, I., et al. (2018). AMBER 2018. San Francisco, CA: University of California.

PubMed Abstract | Google Scholar

Chakrabarty, B., Naganathan, V., Garg, K., Agarwal, Y., and Parekh, N. (2019). NAPS update: network analysis of molecular dynamics data and protein–nucleic acid complexes. Nucleic Acids Res. 47, W462–W470. doi: 10.1093/nar/gkz399

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Su, C., Ke, M., Jin, X., Xu, L., Zhang, Z., et al. (2011). Biochemical and structural insights into the mechanisms of SARS coronavirus RNA ribose 2′-O-methylation by nsp16/nsp10 protein complex. PLoS Pathog. 7:e1002294. doi: 10.1371/journal.ppat.1002294

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y. W., Yiu, C.-P. B., and Wong, K.-Y. (2020). Prediction of the SARS-CoV-2 (2019-nCoV) 3C-like protease (3CL pro) structure: virtual screening reveals velpatasvir, ledipasvir, and other drug repurposing candidates. F1000Research 9:129. doi: 10.12688/f1000research.22457.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Coleman, C. M., and Frieman, M. B. (2014). Coronaviruses: important emerging human pathogens. J. Virol. 88, 5209–5212. doi: 10.1128/JVI.03488-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Darden, T., York, D., and Pedersen, L. (1993). Particle mesh Ewald: An N· log (N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. doi: 10.1063/1.464397

CrossRef Full Text | Google Scholar

Eckerle, L. D., Becker, M. M., Halpin, R. A., Li, K., Venter, E., Lu, X., et al. (2010). Infidelity of SARS-CoV Nsp14-exonuclease mutant virus replication is revealed by complete genome sequencing. PLoS Pathog. 6:e1000896. doi: 10.1371/journal.ppat.1000896

PubMed Abstract | CrossRef Full Text | Google Scholar

Egloff, M. P., Benarroch, D., Selisko, B., Romette, J. L., and Canard, B. (2002). An RNA cap (nucleoside-2′-O-)-methyltransferase in the flavivirus RNA polymerase NS5: crystal structure and functional characterization. EMBO J. 21, 2757–2768. doi: 10.1093/emboj/21.11.2757

PubMed Abstract | CrossRef Full Text | Google Scholar

Fehr, A. R., and Perlman, S. (2015). “Coronaviruses: an overview of their replication and pathogenesis,” in Coronaviruses. Methods in Molecular Biology, Vol. 1282, eds H. Maier, E. Bickerton, and P. Britton (New York, NY: Humana Press), 1–23. doi: 10.1007/978-1-4939-2438-7_1

PubMed Abstract | CrossRef Full Text | Google Scholar

Frauenfelder, H., Sligar, S. G., and Wolynes, P. G. (1991). The energy landscapes and motions of proteins. Science 254, 1598–1603. doi: 10.1126/science.1749933

PubMed Abstract | CrossRef Full Text | Google Scholar

Furuichi, Y., and Shatkin, A. J. (2000). Viral and cellular mRNA capping: past and prospects. Adv. Virus Res. 55:135. doi: 10.1016/S0065-3527(00)55003-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Gohlke, H., Kiel, C., and Case, D. A. (2003). Insights into protein–protein binding by binding free energy calculation and free energy decomposition for the Ras–Raf and Ras–RalGDS complexes. J. Mol. Biol. 330, 891–913. doi: 10.1016/S0022-2836(03)00610-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Guan, W. J., Ni, Z. Y., Hu, Y., Liang, W. H., Ou, C. Q., He, J. X., et al. (2020). Clinical Characteristics of Coronavirus Disease 2019 in China. N Engl J Med. 382, 1708–1720. doi: 10.1056/NEJMoa2002032

CrossRef Full Text | Google Scholar

Harcourt, B. H., Jukneliene, D., Kanjanahaluethai, A., Bechill, J., Severson, K. M., Smith, C. M., et al. (2004). Identification of severe acute respiratory syndrome coronavirus replicase products and characterization of papain-like protease activity. J. Virol. 78, 13600–13612. doi: 10.1128/JVI.78.24.13600-13612.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hawkins, G. D., Cramer, C. J., and Truhlar, D. G. (1995). Pairwise solute descreening of solute charges from a dielectric medium. Chem. Phys. Lett. 246, 122–129. doi: 10.1016/0009-2614(95)01082-K

CrossRef Full Text | Google Scholar

Hawkins, G. D., Cramer, C. J., and Truhlar, D. G. (1996). Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium. J. Phys. Chem. 100, 19824–19839. doi: 10.1021/jp961710n

CrossRef Full Text | Google Scholar

Hou, T., Wang, J., Li, Y., and Wang, W. (2011). Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations. J. Chem. Inf. Model. 51, 69–82. doi: 10.1021/ci100275a

CrossRef Full Text | Google Scholar

Hu, G., Ma, A., Dou, X., Zhao, L., and Wang, J. J. (2016). Computational studies of a mechanism for binding and drug resistance in the wild type and four mutations of HIV-1 protease with a GRL-0519 inhibitor. Int. J. Mol. Sci. 17:819. doi: 10.3390/ijms17060819

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Wang, Y., Li, X., Ren, L., Zhao, J., Hu, Y., et al. (2020). Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. Lancet 395, 497–506. doi: 10.1016/S0140-6736(20)30183-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ichiye, T., and Karplus, M. (1991). Collective motions in proteins: a covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations. Proteins 11, 205–217. doi: 10.1002/prot.340110305

PubMed Abstract | CrossRef Full Text | Google Scholar

Jakalian, A., Jack, D. B., and Bayly, C. I. (2002). Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J. Comput. Chem. 23, 1623–1641. doi: 10.1002/jcc.10128

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeffrey, G. A. (1997). An Introduction to Hydrogen Bonding. New York, NY: Oxford University Press.

PubMed Abstract | Google Scholar

Jonniya, N. A., and Kar, P. (2020). Investigating specificity of the anti-hypertensive inhibitor WNK463 against With-No-Lysine kinase family isoforms via multiscale simulations. J. Biomol. Struct. Dyn. 38, 1306–1321. doi: 10.1080/07391102.2019.1602079

PubMed Abstract | CrossRef Full Text | Google Scholar

Jonniya, N. A., Sk, M. F., and Kar, P. (2019). Investigating phosphorylation-induced conformational changes in WNK1 kinase by molecular dynamics simulations. ACS Omega 4, 17404–17416. doi: 10.1021/acsomega.9b02187

PubMed Abstract | CrossRef Full Text | Google Scholar

Jonniya, N. A., Sk, M. F., and Kar, P. (2020). A comparative study of structural and conformational properties of WNK kinase isoforms bound to an inhibitor: insights from molecular dynamic simulations. J. Biomol. Struct. Dyn. doi: 10.1080/07391102.2020.1827035. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, P., and Knecht, V. (2012a). Energetic basis for drug resistance of HIV-1 protease mutants against amprenavir. J. Comput. Aided Mol. Des. 26, 215–232. doi: 10.1007/s10822-012-9550-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, P., and Knecht, V. (2012b). Energetics of mutation-induced changes in potency of lersivirine against HIV-1 reverse transcriptase. J. Phys. Chem. B 116, 6269–6278. doi: 10.1021/jp300818c

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, P., and Knecht, V. (2012c). Origin of decrease in potency of darunavir and two related antiviral inhibitors against HIV-2 compared to HIV-1 protease. J. Phys. Chem. B 116, 2605–2614. doi: 10.1021/jp211768n

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, P., and Knecht, V. (2012d). Mutation-induced loop opening and energetics for binding of tamiflu to influenza N8 neuraminidase. J. Phys. Chem. B 116, 6137–6149. doi: 10.1021/jp3022612

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, P., Lipowsky, R., and Knecht, V. (2011). Importance of polar solvation for cross-reactivity of antibody and its variants with steroids. J. Phys. Chem. B 115, 7661–7669. doi: 10.1021/jp201538t

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, P., Lipowsky, R., and Knecht, V. (2013). Importance of polar solvation and configurational entropy for design of antiretroviral drugs targeting HIV-1 protease. J. Phys. Chem. B 117, 5793–5805. doi: 10.1021/jp3085292

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, P., Seel, M., Hansmann, U. H., and Höfinger, S. (2007a). Dispersion terms and analysis of size-and charge dependence in an enhanced poisson– Boltzmann approach. J. Phys. Chem. B 111, 8910–8918. doi: 10.1021/jp072302u

PubMed Abstract | CrossRef Full Text | Google Scholar

Kar, P., Wei, Y., Hansmann, U. H., and Höfinger, S. (2007b). Systematic study of the boundary composition in Poisson Boltzmann calculations. J. Comput. Chem. 28, 2538–2544. doi: 10.1002/jcc.20698

PubMed Abstract | CrossRef Full Text | Google Scholar

Karplus, M., and Kushick, J. N. (1981). Method for estimating the configurational entropy of macromolecules. Macromolecules 14, 325–332. doi: 10.1021/ma50003a019

CrossRef Full Text | Google Scholar

King, A. M. Q., Adams, M. J., Carstens, E. B., and Lefkowitz, E. J. (2011). Virus Taxonomy: 9th Report of the International Committee on Taxonomy of Viruses. San Diego, CA: Academic Press

Google Scholar

Kollman, P. A., Massova, I., Reyes, C., Kuhn, B., Huo, S., Chong, L., et al. (2000). Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889–897. doi: 10.1021/ar000033j

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozakov, D., Grove, L. E., Hall, D. R., Bohnuud, T., Mottarella, S. E., Luo, L., et al. (2015). The FTMap family of web servers for determining and characterizing ligand-binding hot spots of proteins. Nat. Protoc. 10, 733–755. doi: 10.1038/nprot.2015.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Kräutler, V., van Gunsteren, W. F., and Hünenberger, P. H. (2001). A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 22, 501–508. doi: 10.1002/1096-987X(20010415)22:5<501::AID-JCC1021>3.0.CO;2-V

CrossRef Full Text | Google Scholar

Kyte, J., and Doolittle, R. F. (1982). A simple method for displaying the hydropathic character of a protein. J. Mol. Biol. 157, 105–132. doi: 10.1016/0022-2836(82)90515-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, S., Chen, H., Ye, F., Chen, Z., Yang, F., Zheng, Y., et al. (2020). Crystal structure of SARS-CoV-2 nsp10/nsp16 2′-O-methylase and its implication on antiviral drug design. Signal Transduct Target Ther. 5:131. doi: 10.1038/s41392-020-00241-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, X., Gong, Z., Xiao, Z., Xiong, J., Fan, B., and Liu, J. (2020). Novel coronavirus pneumonia outbreak in 2019: computed tomographic findings in two cases. Korean J. Radiol. 21, 365–368. doi: 10.3348/kjr.2020.0078

PubMed Abstract | CrossRef Full Text | Google Scholar

Loncharich, R. J., Brooks, B. R., and Pastor, R. W. (1992). Langevin dynamics of peptides: the frictional dependence of isomerization rates of N-acetylalanyl-N′-methylamide. Biopolymers 32, 523–535. doi: 10.1002/bip.360320508

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, H. (2020). Drug treatment options for the 2019-new coronavirus (2019-nCoV). Biosci. Trends 14, 69–71. doi: 10.5582/bst.2020.01020

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, R., Zhao, X., Li, J., Niu, P., Yang, B., Wu, H., et al. (2020). Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding. Lancet 395, 565–574. doi: 10.1016/S0140-6736(20)30251-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lun, Z.-R., and Qu, L.-H. (2004). Animal-to-human SARS-associated coronavirus transmission? Emerg Infect Dis. 10:959. doi: 10.3201/eid1005.040022

PubMed Abstract | CrossRef Full Text | Google Scholar

Machado, M. R., and Pantano, S. (2020). Split the charge difference in two! a rule of thumb for adding proper amounts of ions in MD simulations. J. Chem. Theory Comput. 16, 1367–1372. doi: 10.1021/acs.jctc.9b00953

PubMed Abstract | CrossRef Full Text | Google Scholar

Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713. doi: 10.1021/acs.jctc.5b00255

PubMed Abstract | CrossRef Full Text | Google Scholar

Massova, I., and Kollman, P. A. (1999). Computational alanine scanning to probe protein– protein interactions: a novel approach to evaluate binding free energies. J. Am. Chem. Soc. 121, 8133–8143. doi: 10.1021/ja990935j

CrossRef Full Text | Google Scholar

Minasov, G., Shuvalova, L., Rosas-Lemus, M., Kiryukhina, O., Wiersum, G., Godzik, A., et al. (2020). 1.80 Angstrom Resolution Crystal Structure of NSP16 - NSP10 Complex from SARS-CoV-2. Available online at: http://www.rcsb.org/structure/6W4H (accessed March 18, 2020).

Moreira, I. S., Fernandes, P. A., and Ramos, M. J. (2007). Computational alanine scanning mutagenesis—an improved methodological approach. J. Comput. Chem 28, 644–654. doi: 10.1002/jcc.20566

PubMed Abstract | CrossRef Full Text | Google Scholar

Olsson, M. H., Søndergaard, C. R., Rostkowski, M., and Jensen, J. H. (2011). PROPKA3: consistent treatment of internal and surface residues in empirical p K a predictions. J. Chem. Theory Comput. 7, 525–537. doi: 10.1021/ct100578z

PubMed Abstract | CrossRef Full Text | Google Scholar

Pastor, R. W., Brooks, B. R., and Szabo, A. (1988). An analysis of the accuracy of Langevin and molecular dynamics algorithms. Mol. Phys. 65, 1409–1419. doi: 10.1080/00268978800101881

CrossRef Full Text | Google Scholar

Pillaiyar, T., Meenakshisundaram, S., and Manickam, M. (2020). Recent Discovery and Development of Inhibitors Targeting Coronaviruses. Drug Discov. Today. 25, 668–688. doi: 10.1016/j.drudis.2020.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, D. J., and Brooks, C. L. III. (2004). A modified TIP3P water potential for simulation with Ewald summation. J. Chem. Phys. 121, 10096–10103. doi: 10.1063/1.1808117

PubMed Abstract | CrossRef Full Text | Google Scholar

Rempe, S. B., and Jónsson, H. (1998). A computational exercise illustrating molecular vibrations and normal modes. Chem. Educ. 3, 1–17. doi: 10.1007/s00897980231a

CrossRef Full Text | Google Scholar

Roe, D. R., and Cheatham, T. E. III. (2013). PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data. J. Chem. Theory Comput. 9, 3084–3095. doi: 10.1021/ct400341p

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosas-Lemus, M., Minasov, G., Shuvalova, L., Inniss, N. L., Kiryukhina, O., Wiersum, G., et al. (2020). The crystal structure of nsp10-nsp16 heterodimer from SARS-CoV-2 in complex with S-adenosylmethionine. bioRxiv. doi: 10.1101/2020.04.17.047498

PubMed Abstract | CrossRef Full Text | Google Scholar

Roy, R., Ghosh, B., and Kar, P. (2020). Investigating conformational dynamics of lewis Y oligosaccharides and elucidating blood group dependency of cholera using molecular dynamics. ACS Omega 5, 3932–3942. doi: 10.1021/acsomega.9b03398

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanachai, K., Mahalapbutr, P., Choowongkomon, K., Poo-Arporn, R. P., Wolschann, P., and Rungrotmongkol, T. (2020). Insights into the binding recognition and susceptibility of tofacitinib toward janus kinases. ACS Omega 5, 369–377. doi: 10.1021/acsomega.9b02800

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheahan, T., Sims, A., Leist, S., Schäfer, A., Won, J., Brown, A., et al. (2020). Comparative therapeutic efficacy of remdesivir and combination lopinavir, ritonavir, and interferon beta against MERS-CoV. Nat. Commun. 11:222. doi: 10.1038/s41467-019-13940-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, S., Zhang, S., and Zhang, Q. (2018). Insight into binding mechanisms of inhibitors MKP56, MKP73, MKP86, and MKP97 to HIV-1 protease by using molecular dynamics simulation. J. Biomol. Struct. Dyn. 36, 981–992. doi: 10.1080/07391102.2017.1305296

PubMed Abstract | CrossRef Full Text | Google Scholar

Singh, S., Sk, M. F., Sonawane, A., Kar, P., and Sadhukhan, S. (2020). Plant-derived natural polyphenols as potential antiviral drugs against SARS-CoV-2 via RNA-dependent RNA polymerase (RdRp) inhibition: an in-silico analysis. J. Biomol. Struct. Dyn. doi: 10.1080/07391102.2020.1796810. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Sk, M. F., Jonniya, N. A., and Kar, P. (2020a). Exploring the energetic basis of binding of currently used drugs against HIV-1 subtype CRF01_AE protease via molecular dynamics simulations. J. Biomol. Struct. Dyn. doi: 10.1080/07391102.2020.1794965. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Sk, M. F., Jonniya, N. A., Roy, R., Poddar, S., and Kar, P. (2020b). Computational investigation of structural dynamics of SARS-CoV-2 methyltransferase-stimulatory factor heterodimer nsp16/nsp10 bound to the cofactor SAM. ChemRxiv [Preprint]. doi: 10.26434/chemrxiv.12608795.v1

CrossRef Full Text | Google Scholar

Sk, M. F., Roy, R., Jonniya, N. A., Poddar, S., and Kar, P. (2020c). Elucidating biophysical basis of binding of inhibitors to SARS-CoV-2 main protease by using molecular dynamics simulations and free energy calculations. J. Biomol. Struct. Dyn. doi: 10.26434/chemrxiv.12084207. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Sk, M. F., Roy, R., and Kar, P. (2020d). Exploring the potency of currently used drugs against HIV-1 protease of subtype D variant by using multiscale simulations. J. Biomol. Struct. Dyn. doi: 10.1080/07391102.2020.1724196. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Viswanathan, T., Arya, S., Chan, S. H., Qi, S., Dai, N., Misra, A., et al. (2020). Structural basis of RNA cap modification by SARS-CoV-2. Nat. Commun. 11:3718. doi: 10.1038/s41467-020-17496-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Wallace, A. C., Laskowski, R. A., and Thornton, J. M. (1995). LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions. Protein Eng. Design Select. 8, 127–134. doi: 10.1093/protein/8.2.127

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Wang, W., Kollman, P. A., and Case, D. A. (2006). Automatic atom type and bond type perception in molecular mechanical calculations. J. Mol. Graphics Model. 25, 247–260. doi: 10.1016/j.jmgm.2005.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., and Case, D. A. (2004). Development and testing of a general amber force field. J. Comput. Chem. 25, 1157–1174. doi: 10.1002/jcc.20035

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Sun, Y., Wu, A., Xu, S., Pan, R., Zeng, C., et al. (2015). Coronavirus nsp10/nsp16 methyltransferase can be targeted by nsp10-derived peptide in vitro and in vivo to reduce replication and pathogenesis. J. Virol. 89, 8416–8427. doi: 10.1128/JVI.00948-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiser, J., Shenkin, P. S., and Still, W. C. (1999). Fast, approximate algorithm for detection of solvent-inaccessible atoms. J. Comput. Chem. 20, 586–596.

Google Scholar

Weiss, S. R., and Leibowitz, J. L. (2011). “Coronavirus pathogenesis,” in Advances in Virus Research, eds K. Maramorosch, A. J. Shatkin, and F. A. Murphy (Elsevier) 81, 85–164. doi: 10.1016/B978-0-12-385885-6.00009-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss, S. R., and Navas-Martin, S. (2005). Coronavirus pathogenesis and the emerging pathogen severe acute respiratory syndrome coronavirus. Microbiol. Mol. Biol. Rev. 69, 635–664. doi: 10.1128/MMBR.69.4.635-664.2005

CrossRef Full Text | Google Scholar

Woo, P. C., Huang, Y., Lau, S. K., and Yuen, K.-Y. (2010). Coronavirus genomics and bioinformatics analysis. Viruses 2, 1804–1820. doi: 10.3390/v2081803

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, C., Liu, Y., Yang, Y., Zhang, P., Zhong, W., Wang, Y., et al. (2020). Analysis of therapeutic targets for SARS-CoV-2 and discovery of potential drugs by computational methods. Acta Pharm. Sinica B 10, 766–788. doi: 10.1016/j.apsb.2020.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, B., Shen, H., Zhu, X., and Li, G. (2011). Fast and accurate computation schemes for evaluating vibrational entropy of proteins. J. Comput. Chem. 32, 3188–3193. doi: 10.1002/jcc.21900

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, H., Xie, W., Xue, X., Yang, K., Ma, J., Liang, W., et al. (2005). Design of wide-spectrum inhibitors targeting coronavirus main proteases. PLoS Biol 3:e324. doi: 10.1371/journal.pbio.0030428

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, Z. Q., Chen, D. H., Tan, W. P., Qiu, S. Y., Xu, D., Liang, H. X., et al. (2018). Epidemiology and clinical characteristics of human coronaviruses OC43, 229E, NL63, and HKU1: a study of hospitalized children with acute respiratory tract infection in Guangzhou, China. Eur. J. Clin. Microbiol. Infect. Dis. 37, 363–369. doi: 10.1007/s10096-017-3144-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, P., Yang, X. L., Wang, X. G., Hu, B., Zhang, L., Zhang, W., et al. (2020). A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature 579, 270–273. doi: 10.1038/s41586-020-2012-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Hou, Y., Shen, J., Huang, Y., Martin, W., and Cheng, F. (2020). Network-based drug repurposing for novel coronavirus 2019-nCoV/SARS-CoV-2. Cell Discov. 6:14. doi: 10.1038/s41421-020-0153-3

PubMed Abstract | CrossRef Full Text | Google Scholar