Remdesivir Strongly Binds to Both RNA-Dependent RNA Polymerase and Main Protease of SARS-CoV-2: Evidence from Molecular Simulations.

The outbreak of a new coronavirus SARS-CoV-2 (severe acute respiratory syndrome-coronavirus 2) has caused a global COVID-19 (coronavirus disease 2019) pandemic, resulting in millions of infections and thousands of deaths around the world. There is currently no drug or vaccine for COVID-19, but it has been revealed that some commercially available drugs are promising, at least for treating symptoms. Among them, remdesivir, which can block the activity of RNA-dependent RNA polymerase (RdRp) in old SARS-CoV and MERS-CoV viruses, has been prescribed to COVID-19 patients in many countries. A recent experiment showed that remdesivir binds to SARS-CoV-2 with an inhibition constant of μM, but the exact target has not been reported. In this work, combining molecular docking, steered molecular dynamics, and umbrella sampling, we examined its binding affinity to two targets including the main protease (Mpro), also known as 3C-like protease, and RdRp. We showed that remdesivir binds to Mpro slightly weaker than to RdRp, and the corresponding inhibition constants, consistent with the experiment, fall to the μM range. The binding mechanisms of remdesivir to two targets differ in that the electrostatic interaction is the main force in stabilizing the RdRp-remdesivir complex, while the van der Waals interaction dominates in the Mpro-remdesivir case. Our result indicates that remdesivir can target not only RdRp but also Mpro, which can be invoked to explain why this drug is effective in treating COVID-19. We have identified residues of the target protein that make the most important contribution to binding affinity, and this information is useful for drug development for this disease.

Introduction

An outbreak of a new coronavirus appeared in Wuhan, China, at the end of 2019 and is spreading rapidly in many countries,[1,2] resulting in a pandemic announced by WHO in March 2020.[3,4] 2019 coronavirus disease (COVID-19) causes severe acute respiratory syndrome (SARS) with pathological symptoms such as coughing, fever, shortness of breath, and pneumonia,[5] and critically ill patients may develop a cytokine storm syndrome.[6−8] Compared to the 2002 SARS epidemic caused by SARS coronavirus (SARS-CoV), the COVID-19 mortality rate is lower,[9] but the number of infected cases and deaths is much higher.[2] As of 5 August 2020, more than 18.3 million cases of infection and about 700 thousand deaths were recorded, and thousands more are struggling for their lives in hospitals across the globe. In addition, the reproduction rate of COVID-19 is higher than SARS, which increases the risk of this disease.[10]

Experiments have shown that the sequence similarity between SARS-CoV and the new coronavirus called SARS-CoV-2[11] is about 79%, and both viruses belong to the beta genus of the coronavirus family.[11,12] The virion has a sphere-like shape comprising a single positive strand of RNA, four structural proteins, spike (S) protein, nucleocapsid (N) protein, membrane (M) protein, envelope (E) protein, and non-structural proteins (nsp)[13,14] (Figure ). RNA genome enveloped by the N protein plays a crucial role in virial replication and transcription.

Figure 1

(Upper panel) Schematic representation of SARS-CoV-2 RNA sequence; 5′ untranslated region (5′ UTR); 3′ UTR; open reading frames (ORF): 1a, 1b, 3a, 3b, 6, 7a, 7b, 8a, 8b, and 9b; S (spike); E (envelope); M (membrane); N (nucleocapsid); and non-structural proteins (NSP) from 1 to 16. (Middle panel) PDB structures of the main protease (Mpro) and RdRp. (Bottom panel) Schematic representation of the SARS-CoV-2 virion structure.

Like SARS-CoV, the entry of SARS-CoV-2 into the host cell begins by attaching the S protein on the virial surface to the angiotensin-converting enzyme 2 (ACE2) of the host cell.[15,16] Therefore, S proteins and ACE2 are considered as one of the major drug targets.[17] Once entered, virus replication starts using host cell resources.

The replication and transcription are facilitated by the assembly of non-structural proteins (nsp), which are produced as a result of the cleavage of viral polyproteins encoded by open reading frame 1a (ORF1a) and ORF1b.[18,19] Canonical RNA-dependent RNA polymerase (RdRp or also known as nsp12) (Figure ) plays a crucial role in the replication and transcription of the SARS-CoV-2 virus because it catalyzes the synthesis of viral RNA. Therefore, RdRp became an important target for drug development to combat coronavirus infections.[17,20] Note that the function of nsp12 is supported by nsp7 and nsp8.[21,22] Using cryo-electron microscopy, Gao et al.[23] and Yin et al.[24] resolved the structure of RdRp in complex with nsp7 and nsp8. Its active site consists of seven A-G motifs (Figure ), where nsp12 performs its function.

Regarding the mechanism of infection and pathogenicity of SARS-CoV-2, proteases play an important role in viral structure assembly and replication.[25,26] In coronaviruses, ORF1a encodes a main protease (Mpro) (Figure ), which is also called a chymotrypsin-like cysteine protease (3CLPro).[27,28] Mpro has a mass of about 33.8 kDa and is embedded in the nsp5 region, which is encoded by the SARS-CoV-2 RNA sequence (Figure , upper part). When RNA is translated into protein, Mpro itself is cleaved from the entire protein sequence via autocleavage.[29−31]

Mpro is composed of a homodimer, which is divided into two protomers that have three different domains.[32] Domains I and II have an anti-parallel structure of β-sheets, while the third domain included α-helices that are connected in parallel with domain II from one side to the other with a loop region. The substrate binding site of Mpro is situated between domains I and II, in which residues His41 and Cys145 are dominant in catalytic activity.[33−36] Mpro plays a key role in coordinating viral replication and transcription of the virus life cycle. It cleaves the major part of polyproteins and releases proteins that have replicative function such as RdRp and RNA-processing domains.[37] Therefore, Mpro becomes a prime target for drugs for SARS-CoV-2.[38,39]

In general, S protein, ACE2, TMPRSS2 (transmembrane protease serine 2), 3CLpro, RdRp, and PLpro (papain-like protease) are widely considered as main targets for antiviral drugs against SARS including SARS-CoV-2 and other coronavirus infections.[17,39] There is currently no new drug developed to treat COVID-19, but some old medications have been shown to be effective like dexamethasone[40] (https://www.nature.com/articles/d41586-020-01824-5), Avifavir, and remdesivir (https://www.nature.com/articles/d41586-020-01295-8). Remdesivir is an adenosine analogue that inhibits viral RNA polymerase with RdRp as its target. It has antiviral activity against multiple variants of the Ebola virus in both cell experiments and monkey models.[41,42]In vitro experiments indicated that remdesivir inhibits SARS-CoV and MERS-CoV viruses by interfering the polymerase function of RdRp.[43−45] Recent evidence suggests that remdesivir improves the status of severe COVID-19 patients,[46] which has forced the European Medicines Agency (EMA) to approve it for treatment for people over 12 years old (https://www.ema.europa.eu/en/human-regulatory/overview/public-health-threats/coronavirus-disease-COVID-19/treatments-vaccines-COVID-19#remdesivir-section). Therefore, understanding the molecular mechanism of the interaction between remdesivir and RdRp and other possible targets is important for the development of COVID-19 therapy.

A recent experiment[44] has shown that remdesivir effectively inhibits the activity of SARS-CoV-2 in vitro, but its target has not been identified. An interesting question that emerged is what is the binding affinity of remdesivir to RdRp and can it strongly bind to Mpro, which is one of the most important targets for drug design of COVID-19. To answer this question, we performed molecular docking, steered molecular dynamics (SMD) simulations, and umbrella sampling. Our SMD results revealed that remdesivir strongly binds to both targets. The SMD method is good for obtaining relative binding affinities,[47] but it is impractical to use to access the equilibrium building free energy ΔGbind since a huge number of trajectories are required. Therefore, we used the umbrella sampling to estimate ΔGbind of Mpro and RdRp, which is in good agreement with the experimental data reported by Wang et al.[44] Importantly, we showed that remdesivir can strongly bind not only to RdRp but also to Mpro, which partly explains why remdesivir is effective in COVID-19 treatment.

Materials and Methods

Structures of Remdesivir and Targets

The structure of remdesivir was taken from PubChem data bank with CID 121304016, and the corresponding 2D and 3D presentations are shown in Figure . It contains 77 atoms, and their indices are given in Figure S1 in the Supporting Information.

Figure 2

(Top panel) 2D structure of Remdesivir, and red sticks divide it into six blocks. Names of these blocks are shown in Figure S1. (Bottom panel) PDB structures of Mpro (6LU7) and RdRp (7BTF). The seven motifs of the RdRp active site are shown in a colored box.

The Mpro and RdRp structures were retrieved from the Protein Databank (PDB) with PDB ID 6LU7[27] and 7BTF[23] (Figure ), respectively. Mpro has one chain with three domains, while RdRp contains three chains corresponding to nsp12 (chain A), nsp7 (chain B), and nsp8 (chain C).

The RdRp PDB file lacks some residues at the termini and some residues that are not at the ends (Table S1 and Figure S2). We have considered two models. For the first model, we added only nonterminal missing residues, while for the second model, all missing residues were added using the MODELLER program package.[48,49] A reason for considering two separate models is that terminal missing residues are more flexible than nonterminal ones. Since the results obtained for these two models are essentially the same, for the sake of clarity, we mainly discuss the first model unless otherwise stated. The effect of missing terminal residues will be briefly discussed for comparison.

Docking Simulation

PDBQT files prepared by AutoDock Tool 1.5.4[50] were used to dock remdesivir to the Mpro (6LU7)[27] and RdRp (7BTF)[23] binding site. AutoDock Vina version 1.1[51] was utilized for docking simulation. For a global search, the exhaustiveness was set to 600, which was sufficient to achieve reliable results, and the dynamics of receptor atoms was neglected.

Molecular Dynamics Simulation

Molecular dynamics (MD) simulation was performed using GROMACS 2020.2 package[52] with the AMBER-f99SB-ILDN force field[53] and the water model TIP3P.[54] Based on the general amber force field (GAFF),[55] the parameters for the remdesivir atoms were generated using Antechamber[56] and Acpype.[57] A simple harmonic function form for bonds and angles and the AM1-BCC[58] charge model were used to calculate atomic point charges. The names, types, masses, and charges of the remdesivir atoms are given in Table S2 in the Supporting Information.

The complexes of remdesivir with Mpro and RdRp were solvated in rectangular boxes with dimensions of 8.3 × 9.4 × 12.8 nm and 12.6 × 13.4 × 15.2 nm, respectively. The total charge of the Mpro–remdesivir complex is −4 e. Complexes RdRp (with nonterminal missing residues)–remdesivir and RdRp (with all residues missing)–remdesivir have a total charge of −15 and −10 e, respectively. Na+ counterions were added to neutralize the system.

To calculate vdW forces, a cutoff of 1.4 nm was adopted, while the particle-mesh Ewald summation method was used to calculate the long-range electrostatic interaction with the same cutoff.[59] The leapfrog algorithm was used to solve the Langevin equations with a time step of 2 fs. After the energy minimization using the steepest descent method, MD simulations with position-restrained Cα atoms were performed to equilibrate the system in the NVT and NPT ensembles of 0.5 and 5 ns, respectively. The temperature and pressure of the system were maintained at 300 K and 1 bar, respectively, using the V-rescale and Parrinello–Rahman algorithms.[60,61]

Before carrying out SMD simulations, the systems have been equilibrated in an MD run of 100 ns without position restraints at 300 K and 1 bar. The last obtained structure was used as the initial conformation for SMD simulations.

Steered Molecular Dynamics (SMD)

To prevent the target from drifting under the external force, the Cα atoms were restrained using a harmonic potential with a spring constant of k = 1000 kJ/nm/mol. In the SMD simulation, a spring is attached from one side to a dummy atom and from the other side to the center of mass of remdesivir. The dummy atom is then pulled from its initial position along the direction defined by the MSH (minimal steric hindrance) method[62] (Figure ) with constant speed v. Hence, the elastic force experienced by the ligand is F = kS(Δz – vt), where Δz is the displacement of the ligand’s center of mass connected with the spring in the pulling direction. As in the AFM experiment[63] and our previous works,[62,64] we chose kS = 600 kJ/mol/nm2 and pulling speed v = 5 nm/ns. Since the rupture force[65] and non-equilibrium work[62,66] depend on v but the relative binding affinities are not sensible to it,[62,66] we restricted to this choice of v. We performed 200 SMD trajectories for each protein–ligand complex.

Figure 3

Remdesivir is pulled out from the 6LU7 and 7BTF binding site in the direction determined by the MHS method (arrows). In the green part of the arrow, the space window used in umbrella sampling is 0.08 nm, while the 0.18 nm window was used in the red part.

Umbrella Sampling Method

Combining SMD and the Jarzynski’s equality (JE),[67] in principle, we can calculate the equilibrium binding free energy, but this task is not feasible because an enormous number of SMD runs are required.[68] Other molecular dynamics-based methods such as molecular mechanics Poisson–Boltzmann or generalized Born surface area (MM-PB/GB SA),[69] linear interaction energy (LIE),[70] free energy perturbation,[71] linear thermodynamics integration,[72] and umbrella sampling[73] can be used to estimate the binding free energy. However, we chose umbrella sampling[73] in combination with the WHAM analysis to estimate ΔGbind as umbrella sampling is one of the exact methods.

We calculated the potential of mean force (PMF), which describes the interaction of remdesivir with the receptor, along the line indicated by the arrow in Figure . This line is aligned along the pulling direction Z determined by the MSH method.[62] One end of it is the starting point of remdesivir in the SMD simulation, while the second corresponds to the end of the SMD simulation of 1000 ps, which corresponds to a distance of 5.06 nm for both complexes. To carry out umbrella sampling, the line was divided into two parts (Figure ). For the green segment, corresponding to a distance between 0 and 2 nm, where the interaction of the ligand with the target is still strong (see Figure ), we used a window of 0.08 nm. In the red segment between 2 and 5.06 nm (Figure ), the receptor–ligand interaction is weak; a 0.18 nm window was chosen. In total, we have 2/0.08 + 3.06/0.18 = 42 windows for Mpro and RdRp.

Figure 5

(A) Force–time profiles of a representative SMD trajectory of Mpro (green) and RdRp (red). The rupture force appears at tmax. (B) Same as in (A), but the force is plotted as a function of the ligand displacement. Arrows refer to a distance of 2 and 5.06 nm. In the 0–2 nm range, we used a window of 0.08 nm in umbrella sampling, while a 0.18 nm window was chosen for a distance between 2 and 5.06 nm for both complexes. (C) Free energy profile obtained by using JE and SMD data. TS, bound, and unbound denote a transition, bound, and unbound state, respectively.

The bias harmonic potential was used to keep remdesivir near the center of each window:

with the umbrella force constant k = kS = 600 kJ/mol/nm2, z is the center of umbrella i. For each space window, a 100 ns MD simulation was performed at 300 K and 1 bar with an initial configuration selected from the SMD trajectory at the middle of the window. The weight histogram analysis method[74] was utilized to analyze the results using the WHAM tool in GROMACS package.[75] Errors were calculated using the bootstrap method.[75]

Quantities Used in Data Analysis

The experimental binding free energy ΔGexpbind was obtained from the EC50 value using the formula ΔGexpbind = RTln(EC50), where RT = 0.597 kcal/mol at 300 K, and EC50 is measured in M. The backbone root mean square deviation (RMSD) was used to measure the deviation of the receptor structure with respect to its initial configuration. A hydrogen bond (HB) was formed if the distance between donor D and acceptor A is less than 3.5 Å, the H–A distance is less than 2.7 Å, and the D–H–A angle is larger than 135 degrees. A sidechain contact between remdesivir and the receptor residue is formed if the distance between their centers of mass is less than 0.65 nm. The hydropathy index of residues was obtained from Kyte and Doolittle.[76] The 2D contact network of remdesivir interacting with the target was constructed using Ligplot+ software package.[77]

Using the force-displacement profile obtained in SMD simulation, the pulling work W was calculated using the trapezoidal rule:[66]where N is the number of simulation steps, and F and x are the force experienced by the ligand and position at step i, respectively. To estimate the non-equilibrium binding free energy from SMD simulations, we used JE equality in the presence of external force with constant pulling speed v:[67,78]where ⟨...⟩ is the average over N trajectories, zt is time-dependent displacement, and Wt is the non-equilibrium or pulling work at time t, i.e., Wt = W(t), where W is defined by eq .

Eq means that we can extract an equilibrium quantity by assembling the external work of infinite number of non-equilibrium processes. In this study, when the transformation is not slow enough and the number of SMD runs is finite, we can obtain only the Jarzynski’s non-equilibrium binding free energy ΔGneqJar[66] Therefore, ΔGneqJar is defined by eq but for the non-equilibrium case.

Results and Discussions

Docking Simulation: Binding Sites of Remdesivir in Mpro and RdRp

The configurations obtained in the best docking modes of remdesivir in complex with Mpro 6LU7[27] and RdRp 7BTF[23] are presented in Figure . The docking binding energies of remdesivir with Mpro and RNA polymerase are −7.9 and – 6.5 kcal/mol, respectively, implying that remdesivir binds to Mpro more strongly than to RdRp. The in vitro experiment showed that the EC50 value of remdesivir for SARS-CoV-2 is 0.77 μM.[44] Using the formula ΔGexpbind = RTln(EC50), where gas constant R = 1.987 × 10–3 kcal/mol, T = 300 K, and EC50 is measured in M, we obtained ΔGexpbind = −8.4 kcal/mol, which is roughly consistent with our docking result for Mpro. Thus, docking simulations suggest that remdesivir likely binds to Mpro stronger than to its commonly accepted target RdRp. However, this result is an artefact of the crude docking method since, as shown below, more accurate SMD and umbrella sampling provide the opposite answer.

Figure 4

Binding site of Remdesivir in complex with Mpro and RNA polymerase. (Upper panel) Remdesivir is shown in stick, Mpro is in orange, while nsp12, nsp7, and nsp8 are in blue, red, and magenta, respectively. Residues of the target at the binding site are highlighted in green. The active site of nsp12 is shown in seven motifs A-G. (Bottom panel) Boxed areas are rendered in 2D charts. HBs and non-bonded contacts are highlighted in green and red lines, respectively. Letter A in parentheses refers to chain A. Mpro has only one chain, while RdRp has three chains, and nsp12 is designated as chain A.

The binding site of remdesivir in Mpro is located between domains I and II (Figure ). Remdesivir forms 3 HBs and 12 non-bond contacts with this target. The Mpro residues that form HB with remdesivir are His163, Ser144, and Leu141, and non-bonded contacts are associated with Glu166, Cys145, Met165, Gln189, Arg188, Asp187, His41, Met49, Thr26, Leu27, Thr45, and Thr25. This result indicates that non-bonded contacts govern the interaction between remdesivir and Mpro.

According to Jin et al., the key residues in the binding pockets of inhibitor N3 in Mpro are His41, Tyr54, Met49, Phe140-Cys145, His163–Pro168, His172, and Asp187–Gln192 regions.[27] These areas are similar to the remdesivir binding pocket, indicating that both N3 and remdesivir bind to His41 and Cys145 of the active site of Mpro.

In the case of RNA polymerase, the binding site of remdesivir is close to the active site of nsp12 (Figure ), indicating that remdesivir can affect the function of nsp12. The active site of nsp12 comprises seven motifs A-G. Motifs A, B, and F have residues that are located at the remdesivir binding site. Nsp12 and remdesivir form four HBs and eight non-bond contacts (Figure ). HBs are formed at residues Thr677, Asp757, and Asn688 and non-bond contacts at Asp620, Ser79, Tyr452, Arg550, Lys618, Cys619, Asp615, Thr684, and Ser678 of nsp12. Hence, as in the Mpro case, in the docking simulation, more non-bond contacts are involved in the RdRp–remdesivir stability than HBs.

SMD Results

Remdesivir Binds to RdRp Stronger than to the Main Protease

The profiles of pulling force and position are shown in Figure . The work spent on pulling remdesivir from the Mpro binding site is 106.2 ± 11.6 kcal/mol and Fmax = 716.2 ± 75.7 pN. In the RdRp case, W = 144.6 ± 19.2 kcal/mol and Fmax = 812.5 ± 102.1 pN. Using eq and the ΔG-displacement/time profile (Figure ), we can estimate the non-equilibrium binding free energy ΔGneqJar = ΔG(tend),[66] which is equal to −71.7 ± 1.2 and −89.5 ± 1.2 kcal/mol for Mpro and RdRp, respectively. The large value of ΔGneqJaris due to the fact that the pulling speed is much higher than that used in experiment.[79] Within error bars, Fmax is the same for the two targets, but W of RdRp is greater than that of Mpro, indicating that remdesivir binds to RdRp more strongly than to the main protease. This conclusion is also supported by the result obtained for ΔGneqJar, which is lower for RdRp than Mpro. The SMD result appears to be consistent with experiment, which suggested that remdesivir inhibits corona and Ebola viruses via binding to RdRp.[41,43] However, it contradicts the docking prediction, which shows that the binding affinity for RdRp is lower than for Mpro. This is probably due to the well-known fact that the docking method is not sufficiently accurate in estimating the binding energy.

Although the pulling work of the RdRp complex is higher than that of Mpro, the difference between the two systems is small, suggesting that remdesivir can also bind to Mpro. We will verify this by performing additional umbrella sampling.

Binding and Unbinding Free Energy Barriers

Ligand binding and unbinding are barrier-crossing processes as bound and unbound states are separated by a transition state (TS) (Figure ).[66] The binding barrier ΔG‡bind = ΔGTS – ΔGunbound, while the unbinding barrier ΔG‡unbind = ΔGTS – ΔGbound, where ΔGTS, ΔGbound, and ΔGunbound are the free energy of the transition and bound and unbound states, respectively (Figure ). ΔGTS is the maximum in the free energy profile represented as a function of displacement/time, ΔGbound is determined at zero displacement/time, and ΔGunbound corresponds to the large displacement (the end of simulation) at which the ligand becomes free. Thus, ΔGbound = ΔG(t = 0) and ΔGunbound = ΔG(tend).[66]

Because the number of SMD runs is limited and the pulling speed is high, we can calculate only non-equilibrium binding and unbinding barriers, but they are still useful for predicting relative binding and unbinding times.[66] For both complexes, ΔG‡bind > ΔG‡unbind (Table ), suggesting that remdesivir binds to the target at a longer time scale compared to unbinding. This is consistent with a general experimental trend[80] showing that the ligand exits the binding pocket faster than it joins from the bulk.

Table 1

SMD Results from 200 Independent Trajectoriesa

target	Fmax (pN)	tmax (ps)	W (kcal/mol)	ΔGneqJar (kcal/mol)	ΔG‡unbind (kcal/mol)	ΔG‡bind (kcal/mol)
Mpro	716.2 ± 75.7	207.8 ± 26.5	106.2 ± 11.6	–71.7 ± 1.2	26.3 ± 7.2	132.2 ± 16.9
RdRp	812.5 ± 102.1	235.3 ± 33.1	144.6 ± 19.2	–89.5 ± 1.2	30.8 ± 9.5	175.1 ± 26.1

Errors represent standard deviations.

Since ΔG‡unbind of Mpro (26.3 kcal/mol) is lower than RdRp (30.8 kcal/mol) (Table ), remdesivir should escape from the Mpro binding site faster than from RdRp. To support this conclusion, we calculated the difference of the solvent-accessible surface area (SASA) between the complex and total SASA of the unbound protein and remdesivir, ΔSASA = SASAcomplex – SASAprotein – SASAremdesivir (Figure S3). Obviously, ΔSASA of Mpro–remdesivir reaches 0 faster than RdRp–remdesivir, which indicates that remdesivir leaves the Mpro binding site faster. This conclusion is further confirmed by the time dependence of the number of contacts between remdesivir and the target (Figure S4) since the contacts disappear at about 450 and 800 ps for Mpro and RdRp, respectively.

The binding barrier of Mpro–remdesivir (132 kcal/mol) is less than that of RdRp–remdesivir (175 kcal/mol) (Table ), implying that the remdesivir binding process to RdRp is slower than to Mpro. However, this result should be treated with caution as fast SMD does not produce equilibrium results, and it is unclear whether the relative binding barrier in equilibrium remains unchanged.

Umbrella Sampling Results

Remdesivir Strongly Binds to Mpro and RdRp

Because the results obtained by using SMD are not valid for equilibrium, we utilized umbrella sampling to calculate the equilibrium binding free energy ΔGbind. The potential of mean force (PMF) was determined as a function of a distance to the binding site along the pulling direction as described in the section of Materials and Methods (Figure ). The presence of local minima reflects the complexity of the binding/unbinding process.

Figure 6

Dependence of the potential of mean force on the reaction coordinate Z of two systems. The vertical arrow represents the binding free energy of Remdesivir to Mpro and RdRp. Configurations below the blue arrow at the bottom of the global minimum of PMF were used for data analysis.

The binding free energy ΔGbind was defined as the difference between the maximum and minimum in the PMF profile (Figure ), and we obtained ΔGbind = −8.69 ± 0.36 and −9.34 ± 0.38 kcal/mol for the Mpro–remdesivir and RdRp–remdesivir complexes, respectively. Thus, within the margin of error, remdesivir shows the same binding affinity to both targets.

Relative binding affinity ΔGbind (RdRp)/ΔGbind (Mpro) = 9.34/8.69 = 1.07, which is not so far from the non-equilibrium ratio ΔGneqJar(RdRp)/ΔGneqJar(Mpro) = 89.5/71.7 = 1.25 obtained using JE and SMD. Although the difference in absolute binding affinity is huge (ΔGneqJar(RdRp)/ΔGbind (RdRp) = 89.5/9.34 ≈ 9.6 and ΔGneqJar(Mpro)/ΔGbind (Mpro) ≈ 8.3), the results obtained at non-equilibrium are useful for comparing the binding affinity of remdesivir to various targets. This result is in the line with previous works.[47,62,64,66]

As mentioned above, an in vitro experiment showed[44] that remdesivir binds to novel coronavirus (2019-nCoV) with EC50 ≈ 0.77 μM, which corresponds to ΔGexpbind ≈ −8.4 kcal/mol. This value is very close to our theoretical estimate. The experiment has revealed that remdesivir binds to RdRp,[43] and this is consistent with our in silico results. More importantly, we have shown that this ligand can also associate with Mpro with EC50 in the range of μM. In other words, together with RdRp, Mpro can be a target for remdesivir. We anticipate that the ability to bind to two drug targets is related to the fact that remdesivir is effective for combating COVID-19.

Stability of the Mpro–Remdesivir Complex is Mainly Controlled by the vdW Interaction, while the Electrostatic Interaction is More Important for the RdRp–Remdesivir Complex

We calculated the non-bonded interaction energy between remdesivir and two targets at the global minimum of the PMF curves (Figure ) because in equilibrium, most of the time, the system will be near this minimum. The results were averaged over the configurations located below the blue line shown in Figure . The electrostatic and vdW energies of remdesivir and Mpro are −25.53 ± 0.27 and −32.73 ± 0.09 kcal/mol, respectively, which shows that vdW interaction rules the stability of this complex. The similar molecular mechanism was also observed for lopinavir and ritonavir interacting with Mpro of SARS-CoV-2.[81] For the remdesivir–RdRp complex, in contrast to the Mpro case, the electrostatic interaction (−42.99 ± 0.43 kcal/mol) is stronger than the vdW interaction (and −30.72 ± 0.31 kcal/mol). Thus, the binding mechanism is sensitive to the target.

Most Important Residues

To find residues of the target that make important contribution to the stability of the studied complexes, we calculated the non-bonded energy of the interaction of protein residues with remdesivir at the minimum of the PMF curve (Figure ). The per-residue interaction energy profiles (Figures S5 and S6) show that the vdW and electrostatic interaction plays a key role in the binding of remdesivir to Mpro and RdRp, respectively. Assuming that the most important residues contribute to the complex stability more than 2 kcal/mol, we can demonstrate that Mpro and RdRp have 8 and 11 residues (Figures and 8), respectively. These residues are close to remdesivir.

Figure 7

Non-bonded interaction energy between residues of Mpro and Remdesivir at the global minimum of PMF. Residues with energies below or equal to −2 kcal/mol are labeled.

Figure 8

Non-bonded interaction energy between the residues of RdRp and Remdesivir at the minimum of PMF. Residues that have the absolute energy larger than or equal to 2 kcal/mol are labeled. The residue index range of nsp12 is from 1 to 929, nsp8 is from 930 to 1055, and nsp7 is from 1056 to 1236.

For Mpro, the regions M165–P168 and Q189–Q192 preserve strong interactions. The negatively charged residue E166 of Mpro has the lowest interaction energy of about −18.8 kcal/mol (Figure ). In the case of RdRp, both positively charged (K542, K548, R550, R552, K618, and K795) and negatively charged residues (D449, D615, D620, and D757) make an important contribution to the binding affinity (Figure , Table ). Thus, we have 10 discernible residues (four negatively charged and six positively charged) versus one negatively charged E166 residue in the case of Mpro (Table ). The difference is also clearly visible on the charge surface at the binding pocket of remdesivir (Figure S7), indicating that the charge distribution is denser in RdRp. This further confirms the fact that the electrostatic interaction is dominant in the RdRp–remdesivir interaction.

Table 2

Most Important Residues of the Two Targetsa

targets	most important residues
Mpro	M165, E166, L167, P168, Q189, T190, A191, and Q192
RdRp	D449, K542, K548, R550, A551, R552, D615, K618, D620, D757, and K795

The energy of their interaction with Remdesivir is less than −2 kcal/mol. The red color refers to residues that contribute to the complex stability above 8 kcal/mol.

The total charge and total hydropathy of significant Mpro residues are −1 e (E166) and −5.3, respectively. In RdRp, these values are 2 e and −36.8, respectively, where only residue A551 is neutral (Table ). Since the overall hydropathy of potent residues in Mpro is higher than in RdRp (see also the hydrophobicity surfaces of the two systems in Figure S8), the interaction between remdesivir and Mpro is dominated by the vdW interaction, while in the RdRp case, the electrostatic rules the interaction.

Per-Atom Interaction Energy of Remdesivir

The per-atom distributions of the non-bonded interaction energy of remdesivir with two targets were also calculated at the bottom of the global minimum of PMF (Figure S9). They are similar for both complexes, as the 10–50 atoms (Figure S1) have strong interactions, and the atom P with index 32 has the lowest energy. The electrostatic energy per atom dominates the vdW term for both systems (Figure S9), and this seems to contradict the above analysis, showing that the electrostatic interaction is dominant only in the case of RdRp. However, the contributions of the electrostatic interaction of remdesivir atoms cancel each other out, and, consequently, for Mpro, the average electrostatic and vdW energies of remdesivir are −0.3 and −0.4 kcal/mol, respectively. In the case of RdRp, these values are −0.7 and – 0.5 kcal/mol, respectively. Therefore, consistent with the previous section, the vdW term is more important for Mpro–remdesivir, while the opposite is true for RdRp–remdesivir. From a drug development standpoint, the important role of the 10–50 region suggests that this region can be modified to improve the binding affinity of drug candidates.

To access the contribution of the different parts of remdesivir, we have divided it into six blocks, as shown in Figure and Figure S1. The names of these blocks are shown in Figure S1. The electrostatic interaction of block 1 (pyrrolo[1,2-f][1,2,3,4] triazin-4-amin) destabilizes Mpro–remdesivir but stabilizes RdRp–remdesivir (Figure ). Blocks 2 (2-metyl-3,4-dihydroxy-5-xianhua tetrahydrofuran), 4 (benzene), and 6 (isopentane) promote the stability of remdesivir and two targets. Block 3 (phosphate group) disfavors the remdesivir binding to both targets via electrostatic repulsion. Block 5 (acid propylic) destabilizes the Mpro–remdesivir complex, and it has a weak stabilization effect on the RdRp–remdesivir complex (Figure ).

Figure 9

Non-bonded energies of Remdesivir blocks interacting with Mpro and RdRp. The structure of the six blocks is shown on the top. Results were obtained at the global minimum of PMF.

Thus, our result indicates that the presence of benzene, [2-metyl-3, 4-dihydroxy-5-xianhua tetrahydrofuran], and isopentane will enhance the binding affinity toward the Mpro and RdRp targets. This information is likely useful for the development of drugs for COVID-19.

Hydrogen Bonds

We studied the receptor–ligand HB network formed during umbrella sampling simulations. This network involving 19 and 24 residues of Mpro and RdRp and the population of each HB are shown in Figure . In the Mpro case, the HB is significantly populated with residues E166 (49.2%), Q189 (19.7%), Q192 (18.6%), and N142 (5.7%). In RdRp, similar residues are D615 (12.7%), R553 (10.6%), K618 (10.3%), D449 (9.8%), D620 (9.5%), R552 (8.9%), D757 (8.7%), and R550 (6.2%). These residues also have a low energy of non-bonded interaction, except N142 of Mpro and R553 of RdRp (Figures and 8). Our result is reasonable because residues that have strong interactions with remdesivir are located close to it, increasing[66] the likelihood of HB formation.

Figure 10

HB network formed by Remdesivir with two targets. Carbon atoms of Remdesivir and protein are highlighted in gold and gray, respectively. The acceptor HB atoms in Remdesivir are shown by balls. The number under the residue name refers to the population (%) of HB formed by the residue and Remdesivir. Results were obtained at the global minimum of PMF.

Effect of Terminal Missing Residues on Binding Affinity of Remdesivir to RdRp

The PDB structure 7BTF of RdRp lacks several segments (Table S1 and Figure S2). In the previous sections, we presented the results obtained for a model in which nonterminal missing residues were added. Here, we investigate the effect of terminal missing residues on the binding affinity of remdesivir to RdRp. Most of these residues are located in nsp8 (Table S1).

In the presence of terminal missing residues, the docking binding energy slightly increases from −6.5 to −6.3 kcal/mol. We performed SMD simulations with the same conditions described above and obtained W = 144.2 ± 14.1 kcal/mol and Fmax = 793.80 ± 88.1 pN, which are close to the reported above W = 144.6 ± 19.2 kcal/mol and Fmax = 812.5 ± 102.1 pN for the system without terminal missing residues.

Using the umbrella sampling method, we obtained a free binding energy of −9.49 ± 0.59 kcal/mol, which is equivalent to the model without terminal missing residues (−9.34 ± 0.38 kcal/mol). Consequently, terminal missing residues have little effect on the stability of the RdRp–remdesivir complex. This result is reasonable because these residues are located far enough from the binding site (Figure S2).

Conclusions

Combining various computational tools, we have studied the association of remdesivir with RdRp and Mpro. Molecular docking shows that remdesivir is associated with RdRp weaker than with RdRp, but this finding has not been supported by more advanced SMD and umbrella sampling techniques. Using the latter method, we showed that, in accordance with an in vitro experiment, remdesivir inhibits the 2019-nCoV activity with an EC50 of μM. Importantly, we predict that together with RdRp, Mpro is also a target for remdesivir, which can be used to understand the high efficacy of this repurposed drug. It would be interesting to verify this prediction by in vitro and in vivo experiments.

Our study revealed that the binding of remdesivir to Mpro and RdRp occurs via different molecular mechanisms that the vdW interaction plays a primary role in the stability of the Mpro–remdesivir complex, while the electrostatic interaction is dominant for the RdRp–remdesivir case. Information on the strongest residues of the two targets in association with remdesivir may be useful for the development of potential drugs for COVID-19.

We have studied the binding of remdesivir to Mpro and RdRp, and the extension of this work to other targets will be important from the point of view not only of basic research but also of medical applications.