Structural Basis of the Potential Binding Mechanism of Remdesivir to SARS-CoV-2 RNA-Dependent RNA Polymerase.

Starting from late 2019, the coronavirus disease 2019 (COVID-19) has emerged as a once-in-a-century pandemic with deadly consequences, which urgently calls for new treatments, cures, and supporting apparatuses. Recently, because of its positive results in clinical trials, remdesivir was approved by the Food and Drug Administration to treat COVID-19 through Emergency Use Authorization. Here, we used molecular dynamics simulations and free energy perturbation methods to study the inhibition mechanism of remdesivir to its target SARS-CoV-2 virus RNA-dependent RNA polymerase (RdRp). We first constructed the homology model of this polymerase based on a previously available structure of SARS-CoV NSP12 RdRp (with a sequence identity of 95.8%). We then built a putative preinsertion binding structure by aligning the remdesivir + RdRp complex to the ATP bound poliovirus RdRp without the RNA template. The putative binding structure was further optimized with molecular dynamics simulations. The resulting stable preinsertion state of remdesivir appeared to form hydrogen bonds with the RNA template when aligned with the newly solved cryo-EM structure of SARS-CoV-2 RdRp. The relative binding free energy between remdesivir and ATP was calculated to be -2.80 ± 0.84 kcal/mol, where remdesivir bound much stronger to SARS-CoV-2 RdRp than the natural substrate ATP. The ∼100-fold improvement in the Kd from remdesivir over ATP indicates an effective replacement of ATP in blocking of the RdRp preinsertion site. Key residues D618, S549, and R555 are found to be the contributors to the binding affinity of remdesivir. These findings suggest that remdesivir can potentially act as a SARS-CoV-2 RNA-chain terminator, effectively stopping its RNA replication, with key residues also identified for future lead optimization and/or drug resistance studies.

Introduction

Coronavirus disease 2019 (COVID-19) is a disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)—a novel coronavirus that has been causing a once-in-a-century pandemic.[1] SARS-CoV-2[2] belongs to the family Coronaviridae, which includes RNA viruses such as severe acute respiratory syndrome coronavirus (SARS-CoV,[3] which caused a pandemic in 2003) and Middle East respiratory syndrome-related coronavirus (MERS-CoV,[4] which has caused a continuing epidemic since 2012). The mortality rate of COVID-19 (in the range of ∼1–6%) is believed to be less deadly than SARS (∼10%)[5] or MERS (∼40%);[5] however, its reproductive number (R0) has been estimated to be 2.0–6.5,[6] higher than SARS and MERS. COVID-19 has been spreading to all continents with multiple epicenters. While certain physical treatment has been shown to assist patients to fight this disease with their own immune systems,[7] no proven remedies exist so far, causing high mortality rates especially in senior groups.[8] This raises high and urgent demand to screen for potential drugs through either drug-repurposing or novel drug development.[9−13]

Remdesivir is a nucleotide analogue that mimics the structure of adenosine. It was originally developed by Gilead Sciences, Inc., to treat Ebola virus disease.[14] Even though it has not passed the phase 3 clinical trial of Ebola treatment, it showed promising improvement over the mortality rate of this deadly disease.[15] In the case of Ebola virus, remdesivir was found to act as an RNA-dependent RNA polymerase (RdRp) binding substrate that replaces ATP in the polymerization process before terminating it.[14] Such an agent is also known as a “chain terminator”. Upon entering the body, remdesivir is hydrolyzed and phosphorylated through metabolism, using the core of the molecule as a nucleoside (GS-441524).[14] We term the hydrolyzed and phosphorylated remdesivir as “RemTP”. Like other nucleotide analogues, remdesivir could potentially be utilized as a broad-spectrum antiviral drug[16] due to the structural similarities of RdRp’s of various viruses.[17] For example, it was clinically tested against MERS-CoV and showed significant efficacy.[18] Currently, phase 3 trials of remdesivir against COVID-19 are under progress in the U.S.[19] It has been shown that remdesivir shortens the time of recovery among COVID-19 patients,[19] leading to the recent approval of it for Emergency Use Authorization (EUA) by Food and Drug Administration (FDA). However, the detailed inhibition mechanism of remdesivir remains unknown. Therefore, we carried out a physics-based molecular modeling study on the binding mechanism between remdesivir and SARS-CoV-2 RdRp. Note that the RdRp complex of coronaviruses has multiple nonstructural protein (NSP) units, such as NSP12, NSP8, and NSP7, among which NSP12 possesses some minimal polymerase activity on its own.[20] A previous study revealed that remdesivir mainly substituted ATP as substrate to NSP12 of SARS-CoV RdRp.[21] The NSP8 and NSP7 cofactors, although found to significantly enhance the polymerase activity of SARS-CoV RdRp,[22] are structurally distant from the nucleoside triphosphate (NTP) binding site.[17] To build a minimal model, we narrowed down our search of the RemTP binding site to NSP12 of SARS-CoV-2 RdRp. We employed homology modeling to first construct the tertiary structure of SARS-CoV-2 NSP12. The initial binding mode of ATP to SARS-CoV-2 NSP12 was subsequently determined by structural alignment to the ATP bound poliovirus RdRp due to the structural resemblance of viral RdRp’s. However, we must emphasize that our proposed binding structure should be considered as a checkpoint structure before the RNA template is inserted (which is therefore considered as a “preinsertion binding mode”), as was the case for the poliovirus RdRp structure.[23] We then performed molecular dynamics (MD) simulations to validate the identified preinsertion binding mode. Upon locating the stable binding mode, we further carried out free energy perturbation (FEP) calculations to identify the key residues in the binding process and estimate the binding affinity of RemTP and ATP to SARS-CoV-2 NSP12.

Method

Homology Modeling

The structure of SARS-CoV NSP12 RdRp was obtained from the Protein Data Bank (PDB ID: 6NUR).[17] The sequence of SARS-CoV-2 NSP12 RdRp (SARS-COV-2 NSP12) was obtained from entry YP_009725307.1 at NCBI.[24] Sequence alignment and homology modeling were performed with MODELER 9.23,[25] with unresolved structures on the N-terminus and C-terminus truncated (gray residues on Figure A), which should not affect our current binding affinity calculations due to their structural distances from the binding pocket.

Figure 1

(A) Sequence alignment between 6NUR (PDB ID) and SARS-CoV-2 NSP12 RdRp (SARS-COV-2 NSP12). Conserved residues are colored in black and labeled with asterisks (*). Mutated residues are colored as follows: brown as hydrophobic, green as hydrophilic, red as positively charged, and blue as negatively charged. The N-terminus and C-terminus are truncated in SARS-COV-2 NSP12 due to the lack of structural information in 6NUR. Residues (r.) 896–905 are filled with a random coil. (B) The homology model of SARS-COV-2 NSP12 (magenta) is overlaid with 6NUR (silver). (C) Based on the study on SARS-CoV NSP12 (6NUR), the domains of SARS-COV-2 NSP12 are colored as follows: green for N-terminal extension (r. 1–397, green), yellow for fingers domain (r. 398–581, 628–687), blue for the palm domain (r. 582–627, 688–815), and red for the thumb domain (r. 816–919). (D) The RMSD from the simulation using the apo form of SARS-COV-2 NSP12.

Molecular Dynamics (MD) Simulation

Based on the cocrystal structure of poliovirus RdRp with ATP[23] (PDB ID: 2ILY), we prepared the initial structure of SARS-CoV-2 NSP12-ATP complex by aligning the “fingers” domain of the two proteins based on the alignment of the sequence (i.e., one from poliovirus and one from COVID-19). Note that this is the initial structure, which will be optimized by MD simulations. The complex structure of RemTP was then aligned with that of ATP for the corresponding simulations. Chelating Mg2+ ions are often needed for these RdRp’s, which are positioned in our simulations based on previous studies on class I RNA polymerase.[26] A total of nine different simulations were performed: one for the apo form of SARS-COV-2 NSP12, two (independent runs) for SARS-COV-2 NSP12-ATP (no Mg2+), two for SARS-COV-2 NSP12-ATP (with Mg2+), two for SARS-COV-2 NSP12-RemTP (no Mg2+), and two for SARS-COV-2 NSP12-RemTP (with Mg2+). The parameters of RemTP were generated by CHARMM CGenFF[27] and listed in Table S2. All other parameters were taken from CHARMM36.[28] Simulations were performed with GROMACS 5.1.2.[29] van der Waals interactions were treated with a switching distance of 10 Å and a smooth cutoff distance of 12 Å. Electrostatic interactions were treated with particle mesh Ewald with a grid size of 1 Å. All simulations lasted about 100 ns before the free energy perturbation calculations.

Free Energy Perturbation (FEP) Calculation

The thermodynamic cycles of our FEP calculations were illustrated in Figure S1. A softcore potential was applied in FEP, and more details can be found in our previous studies.[30] The λ windows were set as (0.00, 0.00001, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 0.999, 0.9999, 0.99999, 1.00) to avoid the well-known FEP end-point catastrophe. NAMD 2.13[31] was used to perform the FEP calculations. The starting FEP complex structures of SARS-COV-2 NSP12 binding with NTP (ATP or RemTP) were selected using the clustering algorithm based on the RMSD of NTP while SARS-COV-2 NSP12 was aligned. The largest cluster of each simulation was taken to carry out the FEP calculations. All other simulation parameters were the same as MD simulations. Each simulation was repeated 10 times with different random seeds to obtain the standard error of the free energy changes.

Results and Discussion

Homology Model of SARS-COV-2 NSP12

We constructed the homology model of SARS-COV-2 NSP12 by first aligning the sequences between SARS-COV-2 NSP12 (Wu et al.,[24] NCBI: YP_009725307.1) and the recently resolved SARS-CoV NSP12 RdRp structure[17] (PDB ID: 6NUR). The sequence identity was 95.8%, with 131 unresolved terminal residues, 10 unresolved hinge residues, and 24 mutated residues (Figure A). The 3D structure of SARS-COV-2 NSP12 was then created with the MODELER package without the 131 unresolved terminal residues. The 10 unresolved hinge residues (residues 896–905) were estimated with the default algorithm, thus resulting in a random coil configuration. Because of the high sequence conservation, the backbone RMSD between 6NUR NSP12 and SARS-COV-2 NSP12 was a mere 0.25 Å (Figure B). The similarity between two proteins has been verified by the newly solved cryo-EM structure (PDB ID: 6M71) of SARS-COV-2 NSP12 (the RMSD is 0.55 Å between PDB IDs 6NUR and 6M71).[32] The newly added random coil and 24 mutated residues have been highlighted in Figure S2. It was previously discovered that RdRp’s of various viruses share the same morphology with common building blocks:[17] the fingers domain (amino acids 398–581, 628–687), thumb domain (amino acids 816–919), and palm domain (amino acids 582–627, 688–815, Figure C). The “grip” (a hole formed between these three domains) served as the binding site for the RNA template and nucleoside triphosphates (NTPs). Additionally, the SARS-CoV NSP12 also featured a unique N-terminal extension, which was mostly conserved in SARS-COV-2 NSP12 (colored by green in Figure C). Comparing Figure C and Figure S2, we notice that the majority of the mutations occur at the N-terminal extension and the palm domain. Meanwhile, the fingers domain and the thumb domain remain highly conserved. Consequentially, the RNA/NTP binding grip is highly conserved. As a simple test on the structural stability, we ran a 100 ns simulation on the apo form of SARS-COV-2 NSP12. The RMSD (compared to the initial structure, Figure D) remains below 3 Å, indicating that the constructed protein is stable.

Binding Mode of the Substrates

No crystal structures exist yet for the ATP-bound SARS-CoV-2 NSP12. Therefore, we referred to another ATP-bound RdRp structure from poliovirus (PDB ID: 2ILY) utilizing the structural conservation of viral RdRp’s. The binding structure of ATP to poliovirus RdRp is shown in Figure A, with residues in close contact with ATP labeled. The triphosphate part of ATP mostly interacts with positively charged residues (we define residues to interact with the substrate when the closest heavy atom distance is below 5 Å, same below), such as K159, R174, R163, K167, K172, and K359. The negatively charged residue D323 in the vicinity implies the existence of a chelated ion (Mg2+), suggested by other studies.[26] The nucleoside part of ATP (adenosine) interacts with more diverse residues, such as K61, I176, E177, D238, and S288, which are mostly hydrophilic with only one hydrophobic residue. To approximate the binding mode of ATP in SARS-COV-2 NSP12, we aligned the fingers domain of poliovirus RdRp with that of SARS-COV-2 NSP12 because this was the main constituent of the grip. The aligned SARS-COV-2 NSP12 is shown in Figure B with both ATP and RemTP (referred to as NTP for either; note that the residue numbers are drastically different because of the size of the protein: 932 residues in SARS-COV-2 NSP12 versus 461 residues in poliovirus RdRp). As expected, the alignment ensured that triphosphates of NTPs are in close contact with positively charged residues K545, R551, R553, and R555. Nucleosides of NTPs are surrounded by residues like T556, V557, A558, C622, D623, and S682, which are still mostly hydrophilic with only some hydrophobic. Our alignment of the preinsertion site here indicates once again the conservation of the binding grip from viral RdRp’s.

Figure 2

In part A, the binding mode of ATP to poliovirus RdRp (PDB ID: 2ILY) is illustrated. The ATP molecule is highlighted in stick models with yellow (carbon), blue (nitrogen), red (oxygen), and tan (phosphorus). Residues that directly contact ATP (within 5 Å of the molecule) are illustrated in thin stick models colored by gray (carbon), blue (nitrogen), red (oxygen), and yellow (sulfur). The rest of the RdRp is illustrated with a new cartoon model with VMD. The fingers domain, thumb domain, and palm domain are labeled to guide the eyes. By aligning the fingers domain, the binding mode of ATP and RemTP to SARS-COV-2 NSP12 is constructed as the initial structure for our simulations (illustrated in part B). Note that in the figure the constructed ATP and RemTP completely overlap. All other illustrations are the same as those in part A.

Binding Mode Optimization

MD simulations were then carried out from the initial binding complex structures prepared above. We ran two sets of simulations for both ATP and RemTP—one set without Mg2+ and one set with Mg2+—to examine the structural influence of Mg2+. The initial position of Mg2+ was constructed based on previous studies on class I RNA polymerase of poliovirus.[26] We must emphasize that at least two Mg2+ ions are needed for the catalytic process in RNA replication, with one structurally bound to NTP.[26,33] However, the exact location of the second Mg2+ remains in debate.[33−35] The inclusion of a second Mg2+ also tended to destabilize the NTP binding (to hepatitis C virus RdRp).[33] Therefore, in our current work, we included at most one Mg2+ ion. More studies are needed in terms of understanding the functions of Mg2+ in SARS-CoV-2 RdRp.

We performed clustering analysis based on the RMSD of NTPs during the simulations, with SARS-COV-2 NSP12 aligned. We noticed that, by adding the Mg2+, the largest cluster size was slightly larger than that without Mg2+ (Figure S1A vs S1B). This indicates that Mg2+ might have stabilized the binding between ATP and SARS-COV-2 NSP12. As a result, the determined “bound state” (defined as the largest binding cluster observed in our simulations) was likely from the first set of simulations of SARS-COV-2 NSP12-ATP with Mg2+ (Figure S1B, upper panel). A closer examination (Figure A) showed that the “bound state” did not deviate much from the initial structure (which was aligned to 2ILY). The triphosphate of ATP was found to interact with S549, K551, R553, R555 (from the fingers domain), R836 (from the thumb domain), and Mg2+. Particularly, the first solvation shell of Mg2+ consisted of one to two phosphate oxygens and four to five water oxygens. Mg2+ was situated between the triphosphate group and the palm domain of NSP12, neutralizing the slightly negatively charged protein surface, featuring D618, D761, and E811. The adenosine of ATP was found to interact with M542, T556, V557, A558, and S682 (from the fingers domain), similar to its initial structure. This indicates a remarkable consistency between the preinsertion binding modes of ATP to RdRp’s from two species: poliovirus and SARS-CoV-2. Additionally, we plotted the solvent accessible surface of SARS-COV-2 NSP12 upon binding to ATP (Figure B). ATP was found to reside well inside a local pocket within the grip. However, the “extra space” as seen on the right side of the adenosine group might suggest a possible druggable target. To further test this assumption, we performed molecular docking with several adenosine/guanosine analogues (see more details in the Supporting Information). Interestingly, several molecules, including RemTP, with an enlarged nucleoside, occupied the top of the list ordered by docking scores (see Table S1 for details).

Figure 4

In part A, we illustrate the most stable binding mode of SARS-COV-2 NSP12-ATP (with Mg2+) based on the clustering result (Figure S1B, upper panel). All representations are the same as those in Figure except that Mg2+ is drawn with a green sphere and that the initial structure of ATP is drawn in transparent pink. In part B, the representation of SARS-COV-2 NSP12 is switched to the solvent accessible surface, highlighting the local pocket.

Likewise, we performed the MD simulations for the putative binding structures of RemTP in SARS-COV-2 NSP12. Similar to that of ATP, the largest cluster size of RemTP with Mg2+ is significantly larger than that without Mg2+ (Figure A vs Figure B). The “bound state” was therefore chosen as the largest cluster of the first set of simulations with Mg2+ (Figure B, upper panel). A closer inspection of this structure showed that it deviated from the ATP binding site somewhat (Figure A) while staying as the most stable binding state in our simulations. This was mainly due to the binding of the additional nitrile group (—C≡N) in RemTP to the shallow pocket formed by K545, Y546, and A547 (from the fingers domain). The triphosphate part of the molecule mainly interacted with positively charged residues K551, R553, and R555 (from the fingers domain), K621 and K798 (from the palm domain), and R836 (from the thumb domain) along with Mg2+ and S549 (from the fingers domain). Similar to the ATP binding state, Mg2+ was situated between the triphosphate group of RemTP and the palm domain of NSP12, with a solvating layer of water separating the Mg2+ and the slightly negatively charged protein surface. The base part of the RemTP mainly interacted with K545, A547, S549, R555, and V557 (from the fingers domain). Interestingly, D618 (from the palm domain), a previously identified residue that was crucial for the SARS-CoV RdRp activity,[21] was found to be directly affected by the binding of RemTP (by RemTP forming a hydrogen bond with K798, which could have originally formed a hydrogen bond with D618). Another interesting finding was that the configuration of RemTP in the preinsertion binding state seemed to “block” the grip of SARS-CoV-2 NSP12 (where the RNA template inserts, Figure B). Therefore, it might potentially slow down the activity of RdRp, on top of its supposed function that RemTP could act as a terminating nucleotide on the RNA.[14] Our optimized binding structure agrees with the RemTP binding structure hypothesized by Gao et al.[32] which showed that R555, V557, and D618 were among the key binding residues.

Figure 3

We identify the clusters of the ligand binding state with the RMSD calculations on the ligand alone (in this case, RemTP), while SARS-COV-2 NSP12 is aligned. A cutoff of 2 Å is chosen. The clustering results from the two simulations with the SARS-COV-2 NSP12-RemTP complex without Mg2+ are plotted in part A. The clustering results from the two simulations with the SARS-COV-2 NSP12-RemTP complex and Mg2+ are plotted in part B. The clusters are plotted using their “persistence time”, defined as the simulation time during which the cluster persists. Larger clusters are plotted with big crosses and red colors. Smaller clusters are plotted with small crosses and blue colors. Upon obtaining the largest clusters of each simulation, we perform free energy perturbation (FEP) calculations using the thermodynamic cycle in part C.

Figure 5

In part A, we illustrate the most stable binding mode of SARS-COV-2 NSP12-RemTP (with Mg2+) based on the clustering result (Figure B, upper panel). All representations are the same as those in Figure except that Mg2+ is drawn with a green sphere, the most stable binding structure of ATP (from Figure A) is drawn in transparent pink, and D618 is emphasized in red with an orange arrow. In part B, the representation of SARS-COV-2 NSP12 is switched to the solvent accessible surface, highlighting the local pocket.

Kd of RemTP Is ∼100 Times Stronger than ATP in the Preinsertion Site

We followed the thermodynamic cycle illustrated in Figure S1C and Figure C to calculate the relative binding free energy of ATP and RemTP (with the zero-point energy reference as the MTP molecule: methanol triphosphate). The binding free energy of the most stable binding structures was then calculated and averaged over 10 runs of FEP calculations. The relative binding free energy (ΔΔG) of ATP was found to be −4.88 ± 0.62 kcal/mol, while that of RemTP was −7.68 ± 0.57 kcal/mol. The difference between the two was ∼2.80 ± 0.85 kcal/mol, equating to ∼94 times of the difference in Kd values under 310 K. It should be noted that FEP can converge slowly particularly on systems where the environments of the target alchemical modification undergo slow response fluctuations; for this purpose, various advanced free energy computing strategies were developed, such as funnel metadynamics[36] and orthogonal space random walk (OSRW).[37] In this study, we adopted a simple approach following our previous studies[38−41] with multiple runs starting from the most probable binding cluster. It would also be desirable in the future to study the ligand–RdRp binding process with more advanced sampling techniques, such as replica exchange,[42] replica exchange with solute tempering,[43] and metadynamics[44] to help reduce the enormous computational resources required (it takes hours if not days in experiment to measure Kd values typically), and better understand the underlying mechanism of the therapeutic effects of remdesivir.

Nevertheless, our calculations strongly suggest that RemTP would almost completely replace the “native ligand” ATP when both presenting at the binding grip of SARS-COV-2 NSP12. This finding supplements the previous study where RemTP was found to bind with the Ebola RdRp.[45] It also suggests that, if administrated correctly and highly metabolized, remdesivir might indeed be a promising drug to significantly lower the reproduction of SARS-CoV-2. This finding agrees with the study by Yin et al.[46] where RemTP completely blocks the activity of SARS-CoV-2 RdRp with 1/10 concentration of ATP. Additionally, the backbone RMSD of our modeled RdRp is only 3.3 Å from the solved cryo-EM structure with the inserted RNA template,[46] indicating a good overall agreement with the experimental structure. More importantly, a structure overlay of our predicted preinsertion state of RemTP and the postinsertion cryo-EM structure (Figure S3) reveals that RemTP is situated perfectly to form hydrogen bonds with the primer RNA, ready to be incorporated in the polymerization process.

Conclusion

In this study, we constructed the homology model of SARS-CoV-2 NSP12 RdRp with high sequence identity (95.8%). All key residues at the RNA/NTP binding site were conserved between SARS-CoV and SARS-CoV-2. Although there has been no cocrystal structure yet for ATP bound with either SARS-CoV or SARS-CoV-2 RdRp’s, we successfully constructed a model for ATP/remdesivir binding with SARS-CoV-2 NSP12 based on a previous cocrystal structure of poliovirus RdRp. The relative binding free energy of ATP (the baseline is MTP, same below) was subsequently calculated to be −4.88 ± 0.62 kcal/mol with the presence of Mg2+. The active metabolite of remdesivir (RemTP) was found to have a relative binding free energy of −7.68 ± 0.57 kcal/mol, which is significantly stronger than ATP. The ∼100-fold difference in the Kd value might decisively block ATP out of the preinsertion site when RemTP is in the vicinity. Subsequently, RemTP could act as an effective SARS-CoV-2 RNA-chain terminator, stopping its RNA replication. More studies are still needed to explore the details of this process, especially on how NSP7 and NSP8 cofactors mediate the entry and incorporation of RemTP. Finally, the previously identified crucial residue D618 was affected by the binding of RemTP, indicating a possible secondary effect of this drug that would slow down the activity of RdRp, thus also helping cure COVID-19.