key: cord-0894923-hpm7q0a6
authors: Tanimoto, Shoichi; Itoh, Satoru G.; Okumura, Hisashi
title: “Bucket brigade” using lysine residues in RNA-dependent RNA polymerase of SARS-CoV-2
date: 2021-07-31
journal: Biophys J
DOI: 10.1016/j.bpj.2021.07.026
sha: 711c103394ebec0ab3dd736f4c65cadd063a62a0
doc_id: 894923
cord_uid: hpm7q0a6

The RNA-dependent RNA polymerase (RdRp) of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a promising drug target for coronavirus disease 2019 (COVID-19) because it plays the most important role in the replication of the RNA genome. Nucleotide analogs such as remdesivir and favipiravir are thought to interfere with the RNA replication by RdRp. More specifically, they are expected to compete with nucleoside triphosphates, such as adenosine triphosphate (ATP). However, the process in which these drug candidates and nucleoside triphosphates are taken up by RdRp remains unknown. In this study, we performed all-atom molecular dynamics simulations to clarify the recognition mechanism of RdRp for these drug candidates and ATP that were at a distance. The ligand recognition ability of RdRp decreased in the order of remdesivir, favipiravir, and ATP. We also identified six recognition paths. Three of them were commonly found in all ligands, and the remaining three paths were ligand-dependent ones. In the common two paths, it was observed that the multiple lysine residues of RdRp carried the ligands to the binding site like a “bucket brigade”. In the remaining common path, the ligands directly reached the binding site. Our findings contribute to the understanding of the efficient ligand recognition by RdRp at the atomic level.

S759 to D761), which are also conserved in many viral RdRps. The motifs E and F play an important role in primer binding, and motif G is mainly involved in RNA template binding. In addition, RdRp sequence shares a high degree of identity with those of other coronaviruses. In particular, SARS-CoV-2 and SARS-CoV nsp12s show more than 96 % sequence identity, and the tertiary structure of SARS-CoV-2 nsp12 is similar to that of SARS-CoV (19, 25) . This data on conservation and homology of SARS-CoV-2 is useful in applying therapeutic agents tested on other RNA viruses.

Nucleotide analogs such as remdesivir and favipiravir (trade name: Avigan) are promising drug candidates targeting the viral RdRps, and have shown potential therapeutic effects. Remdesivir was created by Gilead Sciences as a treatment for Ebola virus disease (26) . Favipiravir was created by Toyama Chemical as an anti-influenza virus agent (27) . These drugs are expected to inhibit replication of RNA. RdRp usually recognizes nucleoside triphosphates (NTPs) such as adenosine triphosphate (ATP) to replicate RNA. In cells, remdesivir is triphosphorylated, and favipiravir is ribosylated and triphosphorylated to become active metabolite form (from now on, they are referred to as RemTP and FavTP, respectively; Fig. S1 ), which are thought to interfere with the RNA replication by RdRp (20, (28) (29) (30) (31) .

Several experimental and computational studies on the activity of these nucleotide analogs have been reported. In vitro and in vivo experiments showed that remdesivir and favipiravir suppress the replication of SARS-CoV, MERS-CoV, and SARS-CoV-2 with the high potential to inhibit the RdRp function (32) (33) (34) (35) (36) (37) (38) . The tertiary structure of an RNA-RdRp complex with remdesivir monophosphate or FavTP was determined by using cryo-electron microscopy (cryo-EM) (20, (28) (29) (30) (31) . Several molecular docking and free energy perturbation (FEP) studies reported the tight binding J o u r n a l P r e -p r o o f of RemTP and FavTP to the binding site of RdRp and the estimated binding poses (24, (39) (40) (41) . The binding mechanisms of remdesivir were investigated using steered molecular dynamics (SMD) and umbrella sampling (42) , and molecular dynamics (MD) simulation in implicit solvent (43, 44) . Quantum mechanics/molecular mechanics (QM/MM) calculations and MD simulations were performed to investigate the inhibition of RdRp activity by RemTP (45, 46) . We recently compared the dynamic properties of RdRps without ligands in SARS-CoV and SARS-CoV-2 using MD simulations (47) .

In this study, we perform MD simulations to clarify how SARS-CoV-2 RdRp takes up and recognizes RemTP, FavTP, and ATP that were at a distance. The MD simulation is a powerful theoretical approach to study the dynamics of drugs (48, 49) and disease-related biomolecules, such as proteins (50) and glycans (51) . Previous studies on the RdRp-systems, including ligands, have focused on the RdRp structure after the drug candidate was recognized (39, 46) . No studies have been reported for the process in which the drug candidates and NTPs have been recognized. This study is the first work to reveal the recognition process of the drug candidates and ATP by RdRp.

We employed the cryo-EM structure of SARS-CoV-2 RdRp (protein data bank (PDB) entry: 7bv2) (20) as an initial one. Since one of the nsp8 in this cryo-EM structure was missing, the nsp8 corresponding to the chain D in the apo RdRp (PDB entry: 7bv1) (20) was added. The missing non-terminal residues in the cryo-EM structure were complemented using the modeller plugin (52) in UCSF Chimera (53) .

Hydrogen atoms were added using the LEaP and reduce plugin (54) in AMBER.

Proteins were described with Amber ff14SB force field (55) . TIP3P parameters (56) were employed for water molecules. The geometries of the ligand molecules (RemTP, FavTP, and ATP) were optimized at Hartree-Fock/6-31+G(d) in vacuo. The restraint electrostatic potential (RESP) method (57) was used to determine the atomic partial charges of the ligand molecules. The geometrical optimization and RESP calculations of the ligand molecules were performed using the Gaussian 16 Revision C.01 package (58) . The Lennard-Jones (LJ) parameters for the ligand molecules were obtained from the general Amber force field (GAFF) parameter version 2 set (59) assigned by Antechamber (60) . The LJ parameters for the divalent and monovalent metal ions were taken from references (61) and (62), respectively. Ligand molecules were distributed using the Generalized-Ensemble Molecular Biophysics (GEMB) program, which was developed by one of the authors (H. O.) (63) . Initial conformations of several protein systems have been prepared, and their simulations have been performed using this program so far (64) (65) (66) . 100 ligand molecules were randomly placed around RdRp at a distance of more than 50 Å from one of the two Mg 2+ ions (residue ID: 1193) at the binding site (Fig. 1) . The distance between the centers of mass of any two ligands was set to 17 Å or more. Using different seeds of random numbers for the ligand arrangement, five initial conformations were prepared for each system. 415 Na + ions were added to neutralize the system, and the total system consisted of about 400,000 atoms (402,188 atoms for the RemTP system, 402,489 for the FavTP system, and 402,465 for the ATP system). The entire system was immersed in a cubic box with a volume of 160 3 Å 3 . The ligand concentration in each system was set to 40.5 mM to observe the spontaneous ligand recognition by RdRp during the simulation time. Note J o u r n a l P r e -p r o o f that this concentration is higher than the physiological concentration of ATP in vivo (about 2-8 mM) (67) and the 50% cytotoxic concentration of remdesivir and favipiravir in vitro (~ 100-400 M) (68).

All MD simulations were performed using the program AMBER18 (69). The system was first minimized with 5,000 steps of the steepest descent method and then with 5,000 steps of the conjugate gradient method. The system with minimized energy was gradually heated from 0 K to 310 K in 310 ps. After the heating phase, short equilibration at 310 K lasted for 690 ps. A production run simulation was then performed for 110 ns. 50 independent MD simulations were conducted for each system with the combinations of five initial conformations and ten initial velocities. All MD simulations were performed in the isothermal-isochoric (NVT) ensemble with T = 310 K except for the heating phase. Langevin heat bath with a damping coefficient of 2 ps -1 was used to keep the temperature constant. The LJ potential cutoff was at 12 Å. Longrange electrostatic interactions were calculated using the Particle Mesh Ewald method (70) with an Ewald coefficient of 0.305 Å -1 , and the grid size set to 160 3 . Periodic boundary conditions were used. Bonds involving hydrogen atoms were constrained using the SHAKE method (71) . The equations of motion were integrated using the leapfrog algorithm with a time step of 2 fs except for the heating phase, where a time step of 1 fs was used. The cpptraj program (72) 

We performed MD simulations for three systems. Each system consists of RdRp with either RemTPs, FavTPs, or ATPs. We observed that the ligands were recognized by RdRp in all three systems. The number of trajectories, 50, and simulation time, 110 ns, in this study are more than those in similar simulations for other proteins with many ligands (49) , and ligand recognition behavior was sufficiently observed within 110 ns.

More specifically, the ligands were taken up at the binding site of RdRp during the simulations, as shown in supplemental movie 1. To investigate the ligand dependence of the ligand recognition of RdRp, we calculated ligand recognition probabilities. Here, when the closest distance between Mg 2+ at the binding site and any atoms of the ligand was less than 3.5 Å, it was considered as a "contact" event. Furthermore, when the contact events were observed in more than 10 snapshots out of the last 1,000 snapshots (corresponding to the last 10 ns) in the simulations, it was considered as a "ligand recognition" event. The ligand recognition probability was obtained by dividing the number of the MD simulations with the ligand recognition events by the total number of the MD simulations (50 MD simulations), as listed in Table 1 . Although within the statistical errors, RemTP shows the highest probability, and FavTP has the secondhighest probability, followed by ATP. These results are qualitatively consistent with previous experimental studies (20, 36) . In addition, our results agree with the MD simulations using the FEP method for the RdRp-RemTP complex, in which it was shown that RemTP binds much stronger to RdRp than ATP (39) . We remark that the ligands sometimes aggregated due to a large number of the counter ions (415 Na + ) in the systems. These aggregates were stabilized by the interaction between the phosphate J o u r n a l P r e -p r o o f groups of the ligands and the counter ions. To eliminate this effect, even if the ligand aggregates were recognized by RdRp, it was not considered as a ligand recognition event. Here, the number of trajectories with aggregation events were only 1 for RemTP, 3 for FavTP, and 2 for ATP.

In order to understand the ligand recognition mechanism by RdRp, trajectories of the recognized ligands were investigated. We found that the ligand recognition paths Fig. 2 ), the ligands are in frequent contact with any residues of D155 to E167 (especially K160) of nsp12, which is shown on the right side of the binding site. In path 3 (cyan line in Fig. 2 ), the ligands enter the binding site from the front and bind to Mg 2+ with almost no interaction with any residues of RdRp. In addition to these common paths, the following three paths (paths 4-6) were also found in one or two of the systems, namely, they are ligand-specific paths. In path 4, which was observed in the RemTP and FavTP systems, the ligands interacted with the histidine residues (H810 and H816 of nsp12) in the vicinity of the binding site. In path 5, which was found only in the FavTP system, FavTP reached the binding site after interacting with Y826 of nsp12. In path 6, which was found only in J o u r n a l P r e -p r o o f the ATP system, ATP bound Mg 2+ after ATP was trapped in the gap between nsp8 and nsp7. We individually discuss the details of these paths below.

We classified the trajectories in which the phosphate group of the ligand had contact with K43 of nsp7 and K438, K551, and K798 of nsp12 before the ligand recognition as path 1. To understand the ligand recognition process in path 1, we analyzed how the ligands interacted with residues of RdRp. The representative trajectory of RemTP in path 1 is shown in Fig. 3 . Snapshots in which RemTP has contact with specific residues are also presented in Fig. 3 . Here, when the closest distance between any atoms of the ligand and those of a residue of RdRp was le ss than 3.5 Å, it was regarded as a "contact" between the ligand and residue. In Figs. 3(b)-(e), the capital letters at the beginnings of the residues (A or C) correspond to the chain labels in the cryo-EM structure of the original PDB (chains A and C are nsp12 and nsp7, respectively). From now on, the residue of RdRp is described as "chain label + residue number + residue name" for convenience.

In this path, interesting behavior of RdRp was observed, in which lysine residues transported the ligand to the binding site as if it were a "bucket brigade" (supplemental movie 2). These lysine residues (C2LYS, C43LYS, A438LYS, and A551LYS) have a positive charge and are arranged in a straight line toward the binding site. In this ligand transport process, the phosphate group of RemTP first interacted with the sidechain of C2LYS (S1; Fig. 3(b) ). Subsequently, C2LYS passed it to C43LYS, which was spatially close (S2; Fig. 3(c) ). C43LYS then passed RemTP to A551LYS J o u r n a l P r e -p r o o f (S3; Fig. 3(d) ), and finally, the phosphate group of RemTP reached the binding site (S4; Fig. 3(e) ). At the binding site, the ligand electrostatically interacted with A621LYS and A798LYS, which was located in the vicinity of the binding site. The same ligandtransport process was observed in the FavTP and ATP systems, except for minor differences (supplemental movies 3 and 4). In the FavTP and ATP systems, C43LYS

instead of C2LYS first contacted with the phosphate group of the ligands. In the FavTP system, A438LYS between C43LYS and A551LYS also participated in the ligand transport. Figure 4 (a) shows that the contact probabilities between RemTP and C2LYS, C43LYS, and A551LYS at the initial stage in path 1 were particularly high. This result also suggests that the electrostatic interaction between the phosphate group of RemTP and these lysine residues plays an important role in the ligand recognition process in path 1. In addition, the probabilities of contact with C43LYS exceeded 0.9 just before that with C2LYS became less than 0.9.

Similarly, the probabilities of contact with A551LYS increased when that with C43LYS decreased. After the probabilities of contact with A551LYS exceeded 0.9, that with J o u r n a l P r e -p r o o f A621LYS and A798LYS exceeded 0.9. These results also indicate that RemTP reached the binding site through these lysine residues in sequence. Several experimental studies have reported that the NTP entry channel of RdRp has positively charged residues placed on a line (18, 20, 73) . However, no experimental evidence has been reported that the lysine residues of RdRp transport the ligands to the binding s ite. This is the first study to show the ligand transport mechanism of RdRp via multiple lysine residues.

To identify the residues involved in the ligand recognition, we set the threshold of the distance between the ligand and RdRp as 3.5 Å. To examine the threshold dependence, in addition to using this threshold, we calculated the average contact probability between the RemTP and each residue before the first contact event between the RemTP and Mg 2+ using the threshold set to 5.0 Å, as shown in Figs. 4 

We can see that the average contact probabilities of the lysine residues, arranged in a straight line toward the binding site, are higher than those of the surrounding residues regardless of the thresholds. Therefore, it is considered that the threshold dependence is small in identifying the residues involved in the ligand recognition.

It has been reported that these positively-charged basic residues, especially those at motif F (residues A544LEU to A555ARG) in nsp12, favor NTP uptake (19, 20) . Moreover, the above lysine residues (A438LYS, A551LYS, A621LYS, A798LYS, C2LYS, and C43LYS) are highly conserved in RdRp of SARS-CoV (21), as shown in Fig. S3 . Therefore, for SARS-CoV, it is assumed that the NTP recognition ability of RdRp is enhanced by carrying NTP to the binding site by these linearly arranged lysine residues.

J o u r n a l P r e -p r o o f Path 2 includes trajectories in which the phosphate group of the ligand was in frequent contact with any residues between D155 and E167 (especially K160) of nsp12 before the ligand was taken up. The analyzed results of the RemTP trajectory are shown in Figs. 5 and 6. The behavior of the lysine residues that carry RemTP to the binding site like a "bucket brigade" was again observed in path 2 (supplemental movie 5). First, the phosphate group of RemTP interacted with A160LYS (S1; Fig. 5(b) ). The ligand was passed to A798LYS (S2; Fig. 5(c) ). After that, the phosphate group electrostatically interacted with the sidechain of A551LYS (S3; Fig. 5 (d)), then with that of A621LYS (S4; Fig. 5 (e)) before reaching the binding site. Such ligand uptake was also observed in the FavTP and ATP systems (supplemental movies 6 and 7).

The time series of contact probabilities between the ligand and each residue in path 2 are shown in Fig. 6(a) . We again focused on the results of RemTP because similar results were obtained in the FavTP and ATP systems. As can be seen in Fig. 6(a) , A160LYS and the ligand were in frequent contact at the beginning of the simulation. It suggests that for this path, the electrostatic interaction between the phosphate group of the ligand and this lysine plays a particularly important role in ligand recognition. Furthermore, as in path 1, when the contact probabilities with A160LYS became less than 0.9, those with A551LYS and A798LYS exceeded 0.9. More specifically, RemTP was passed between these lysine residues and approached the binding site. Because A160LYS is highly conserved in RdRp of SARS-CoV as shown in Fig. S3 , it is speculated that this lysine residue also contributes to the NTP recognition of RdRp in SARS-CoV.

We calculated the average contact probability between the RemTP and each residue of RdRp before the first contact event between the RemTP and Mg 2+ again J o u r n a l P r e -p r o o f using two thresholds of 3.5 Å and 5.0 Å, as shown in Figs. 6(b) and 6(c). We can see that the average contact probabilities of the lysine residues mainly involved in the ligand recognition in path 2 are higher than those of the surrounding residues regardless of the thresholds. Therefore, it is considered that the threshold dependence is also small in identifying the residues involved in the ligand recognition in path 2.

It is worth mentioning that the effect of counter ions was seen in path 2. For example, the phosphate group of RemTP and A167GLU near the binding site interacted with each other via Na + (Fig. S4) . Such a counterion effect was observed because three acidic residues (A161ASP, A164ASP, and A167GLU) are located between A160LYS and the binding site. Na + binds to these negatively charged acidic residues, neutralizing the electrostatic repulsion between the phosphate groups of the ligands and these residues. Such interaction was also observed in the other ligand systems.

This path is the simplest ligand recognition process. The ligands entered the binding site with almost no interaction with any residues of RdRp due to the electrostatic interaction with the two Mg 2+ ions at the binding site. As seen in the time series of the contact probabilities between RemTP and each residue in Fig. S5 , little contact was seen until the ligand reached the RdRp binding site (supplemental movie 8).

Since this path was observed in all ligand systems as well as paths 1 and 2, it is also one of the main ligand recognition processes.

Path 4 was found only in the RemTP and FavTP systems. The analysis results of J o u r n a l P r e -p r o o f the RemTP trajectory in path 4 are shown in Figs. 7 and S6. As a feature of this path, histidine residues in the vicinity of the binding site (A810HIS and A816HIS) mainly contributed to the ligand recognition. As shown in Fig. 7 , the base moiety of RemTP first formed a - stacking with the sidechains of A810HIS and A816HIS (S1; Fig.   7 (b)). Subsequently, the phosphate group of RemTP electrostatically interacted with A438LYS (S2; Fig. 7(c) ). RemTP was passed from A438LYS to A551LYS and A836ARG (S3; Fig. 7(d) ), and then interacted with A551LYS and A798LYS (S4; Fig.   7 (e)). RemTP finally reached the binding site while maintaining interaction with A551LYS. During these processes (S1-S4), the - stacking between the base moiety of RemTP and A810HIS was almost maintained, and the - stacking was broken when RemTP reached the binding site (supplemental movie 9). Figure S6 also shows that before reaching the binding site, RemTP was in contact with A810HIS and A816HIS for a longer period of time compared to other residues. In addition, similar to paths 1 and 2, the basic residues sequentially carried the phosphate group of RemTP to the binding site. As for FavTP, although sequential delivery of the ligand between these basic residues was not observed, the ligand approached the binding site while maintaining - stacking with A810HIS in the same way as RemTP.

This path suggests that residues other than the positively charged basic ones near the binding site also promote the ligand recognition of RdRp. The location of the aromatic residues such as histidine residues in the vicinity of the binding site would make a beneficial contribution to the transport of the NTP and NTP-like inhibitors to the binding site. Although this recognition process mediated by the histidine residues was not observed for ATP, it is expected that a similar path will be obser ved in the ATP system by improving the statistics. Furthermore, these histidine residues are also J o u r n a l P r e -p r o o f conserved in RdRp of SARS-CoV (Fig. S3) , and can also play an important role in the NTP recognition of RdRp in SARS-CoV.

Path 5 and path 6 were observed only in the FavTP system and the ATP system, respectively. Figure 8 shows the trajectory of FavTP in path 5. The time series of contact probabilities between FavTP and each residue are presented in Fig. S7 . As a noticeable difference from the other paths, path 5 is characterized by long-term contact with A819LEU and A826TYR (Fig. S7 ). Since these residues do not have a net electrostatic charge, it is expected that a moiety of the ligand other than the phosphate group interacts with these residues. The analysis of the trajectory in path 5 revealed that the base moiety of FavTP and the sidechain of A826TYR formed - stacking (supplemental movie 10), resulting in high contact probabilities with these residues.

Similarly, we show the trajectory of ATP in path 6 and the time series of contact probabilities with each residue in Fig. 9 and Fig. S8 , respectively. A distinct feature of this path is that the long-term contact of ATP with multiple residues in chains B (nsp8) and C (nsp7) was observed (Fig. S8) . It was attributed to the base moiety of ATP being trapped in the gap between chain B and chain C (see Fig. 9 and supplemental movie 11).

To investigate the key residues of the ligand recognition, we calculated probabilities with which residues were in contact with the ligands before and after the first contact event between the ligands and Mg 2+ . The contact probabilities in the J o u r n a l P r e -p r o o f RemTP system were averaged for all 12 trajectories in which the ligand recognition events were observed. Figure 10 shows the average contact probabilities of the top 20 residues in contact with RemTP before the first contact event between Mg 2+ and RemTP.

The average contact probabilities for FavTP and ATP were also calculated in the same way, as shown in Fig. S9 . In Figs. 10 and S9, the red letters and pink vertical bars showed the residues that were not included in the top 20 ones after the first contact between the ligands and Mg 2+ (Fig. S10) . In other words, the residues shown in red letters indicate that they have been in contact with the ligand only before the ligand recognition. Figure 10 shows that A160LYS and C43LYS, which mainly contribute to the uptake of RemTP in paths 1 and 2, appear only before the ligand recognition. These two residues also appear in the FavTP and ATP systems only before the recognition (Fig.   S9 ). This suggests that these residues are essential in the NTP recognition by RdRp, as discussed in sections 3.2.1 and 3.2.2. Other lysine residues (A438LYS, A551LYS, A621LYS, and A798LYS) identified in paths 1 and 2 also appear in all ligand systems.

Furthermore, in other paths (except for path 3), there are also high contact probabilities between the ligand and these lysine residues before the first contact between the ligand and Mg 2+ (Figs. S6-S8) . Therefore, the sequential electrostatic interaction between the phosphate group of the ligand and these lysine residues is the main driving force for the ligand recognition in paths other than path 3.

Since these lysine residues (A160LYS, A438LYS, A551LYS, A621LYS, A798LYS, and C43LYS) have higher contact probabilities for all ligands, paths 1 and 2 are considered as the main common paths for the NTP recognition by RdRp.

Conversely, the residues identified in path 4 (A810HIS and A816HIS) appear only for J o u r n a l P r e -p r o o f RemTP and FavTP (Figs. 10 and S9(a) ). In addition, the residues identified in path 5 (A819LEU and A826TYR) and those in path 6 (C27LYS and C34GLN) appear only for FavTP and ATP, respectively. Therefore, these results show that paths 4-6 are relatively minor paths.

According to the literature (74) , even divalent ions such as Mg 2+ have an interaction range of only about 10 Å in water at most. Therefore, it is considered that when the ligand happens to approach the binding site while floating in solution, it is attracted to the Mg 2+ ions, rather than a distant ligand is attracted directly to the Mg 2+ ions. In fact, the ligand was suddenly attracted to the Mg 2+ ions when it came near the binding site in path 3 (supplemental movie 8). In addition to the direct interaction with the Mg 2+ ions, the presence of the lysine residues, such as C2LYS (34-36 Å away from the Mg 2+ ions) and A160LYS (21-25 Å away from the Mg 2+ ions), makes it possible to capture ligands that are further distant. We consider that transporting the ligands to the binding site like a "bucket brigade" by these lysine residues also contributes to the recognition of distant ligands by RdRp.

The MD simulations in this study were performed without RNA strands. To date, some structures of the RdRp bound to an RNA duplex were experimentally solved (20, (28) (29) (30) (75) (76) (77) (78) . From these structures, it is clarified that the nascent RNA extends on the opposite side of the ligand entry side. It is also shown that the RNA template strand enters the binding site from the opposite side of the ligand entry side and translates in the same direction as the nascent RNA elongation. Therefore, we consider that the effect of the presence of these RNA strands on the ligand recognition mechanism of RdRp revealed in this study is small.

In this study, we employed a simple model for the divalent metal ions, which is J o u r n a l P r e -p r o o f represented only by the Coulombic and LJ interactions. Other models such as cationic dummy atom models (79) and polarizable models (80) are useful for considering charge transfer, polarization, and covalent interactions. However, our study focuses on the ligand recognition paths before the ligand binding with the metal ions. Therefore, the force field employed in this study is sufficient for our discussion on the ligand recognition of RdRp.

The structure of RdRp was used without RNA as the initial structure in our MD simulations. The binding site is a cavity containing only two catalytic Mg 2+ ions (Fig.   S11 ). Therefore, the ligands can enter the binding site from either front ("Front" in Fig.   S11 ) or back ("Back" in Fig. S11 ) in these MD simulations, although the back side is inherently blocked by the extended RNA. One might expect that the frequency of the ligand entry from the front is the same as that from the back. However, as shown in Table S1 , the number of trajectories approaching from the correct direction (front entry) was more than that from the opposite direction (back entry). Note that the paths we discussed in the previous subsections are only those from the front.

To understand the difference in the spatial accessibility to the binding site, we calculated the solid angle with which outside is seen from the two Mg 2+ ions without being disturbed by any residues of RdRp. Figures S12 and S13 show the parts of the spherical surface created by these solid angles whose vertices are the two Mg 2+ ions.

Here, the initial structure of RdRp was used. The solid angle at the front side is larger than that at the back side for both Mg 2+ ions ( Table 2 ). For Mg 2+ (1) in Fig. S11 , particularly, the solid angle at the front side is approximately twice larger than that at J o u r n a l P r e -p r o o f the back side. This indicates that the ligand is more accessible from the correct direction when it reaches the binding site with little contact with the residues as in path 3. In addition, it is considered that, as seen in other paths, the intermolecular interaction between the ligand and residues, especially the positively charged basic residues, contributes to the ligand recognition of RdRp from the correct direction.

In this study, we performed the MD simulations of SARS-CoV-2 RdRp with the three kinds of ligands, RemTP, FavTP, and ATP, to clarify the recognition mechanism for the drug candidate compounds. It was found that the recognition probability for RemTP was the highest, that for FavTP was intermediate, and that for ATP was the lowest.

We identified six recognition paths (paths 1-6) in this study. Paths 1-3 were commonly found in all ligands, and paths 4-6 were ligand-dependent paths. In path 1, "bucket brigade" behavior was observed, in which the ligands were carried to the binding site by multiple lysine residues. The lysine residues interacting with the phosphate group of the ligand in this path were K2 and K43 of nsp7 and K438, K551, K621, and K798 of nsp12. In path 2, the transport of the ligands from K160 of nsp12 to the lysine residues around the binding site was also observed. These paths are considered as the main ones in the NTP recognition by RdRp together with path 3, in which the ligands directly reach the binding site. Conversely, in paths 4-6, the residues other than the positively charged basic ones near the binding site were mainly involved in the ligand recognition of RdRp. In our simulations, these paths were not common to all ligands and thus may be relatively minor paths. Fig. 2(a) ). The black circles show positions at which RemTP being in stable contact with residues. (b)-(e) Characteristic structure at each state (S1 to S4). In panels (b)-(e), RemTP and the lysine residues involved in the ligand recognition are represented as red and blue stick models, respectively. The first capital letters of residues (A and C) represent the chain labels in the original PDB (residues with "A" and "C" are residues of nsp12 and nsp7, respectively). Fig. 2(d) ). The black circles showed the positions where RemTP was in stable contact with residues. (b)-(e) Characteristic structure at each state (S1 to S4). In panels (b)-(e), RemTP, the lysine residues involved in the ligand recognition, and A167GLU are shown in red, blue, and brown, respectively. The black circles showed positions at which RemTP was in stable contact with residues.

(b)-(e) Characteristic structure at each state (S1 to S4). In panels (b)-(e), RemTP, the lysine residues involved in the ligand recognition, histidine residues, and A836ARG are shown in red, blue, orange, and light blue, respectively. 

The species Severe acute respiratory syndrome-related coronavirus: classifying 2019-nCoV and naming it SARS-CoV-2

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): The epidemic and the challenges

COVID-19: consider cytokine storm syndromes and immunosuppression

Preparing for the Most Critically Ill Patients with COVID-19: The Potential Role of Extracorporeal Membrane Oxygenation

Critical Illness in Patients With COVID-19: Mounting an Effective Clinical and Research Response

The reproductive number of COVID-19 is higher compared to SARS coronavirus

A pneumonia outbreak associated with a new coronavirus of probable bat origin

World Health Organization (WHO). 2021. Weekly Operational Update on COVID-19 -12

A Novel Coronavirus from Patients with Pneumonia in China

The Architecture of SARS-CoV-2 Transcriptome

A new coronavirus associated with human respiratory disease in China

Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding

Pharmacologic Treatments for Coronavirus Disease 2019 (COVID-19): A Review

Analysis of therapeutic targets for SARS-CoV-2 and discovery of potential drugs by computational methods

A SARS-CoV-2 protein interaction map reveals targets for drug repurposing

The Human Coronavirus Disease COVID-19: Its Origin, Characteristics, and Insights into Potential Drugs and Its Mechanisms. Pathogens

RNA-dependent RNA polymerase of SARS-CoV-2 as a therapeutic target

RNA dependent RNA polymerases: Insights from structure, function and evolution

Structure of the RNA-dependent RNA polymerase from COVID-19 virus

Structural basis for inhibition of the RNA-dependent RNA polymerase from SARS-CoV-2 by remdesivir

Structure of the SARS-CoV nsp12 polymerase bound to nsp7 and nsp8 co-factors

One severe acute respiratory syndrome coronavirus protein complex integrates processive RNA polymerase and exonuclease activities

A Structural View of SARS-CoV-2 RNA Replication Machinery: RNA Synthesis, Proofreading and Final Capping

SARS-CoV-2 RNA dependent RNA polymerase (RdRp) J o u r n a l P r e -p r o o f targeting: an in silico perspective

Novel Coronavirus Polymerase and Nucleotidyl-Transferase Structures: Potential to Target New Outbreaks

Therapeutic efficacy of the small molecule GS-5734 against Ebola virus in rhesus monkeys

Mechanism of Action of T-705 against Influenza Virus

Structural Basis for RNA Replication by the SARS-CoV-2 Polymerase

Mechanism of SARS-CoV-2 polymerase stalling by remdesivir

Structural Basis for the Inhibition of the SARS-CoV-2 RNA-dependent RNA polymerase by favipiravir-RTP

Structural basis of SARS-CoV-2 polymerase inhibition by Favipiravir

Coronavirus susceptibility to the antiviral remdesivir (GS-5734) is mediated by the viral polymerase and the proofreading exoribonuclease

Broadspectrum antiviral GS-5734 inhibits both epidemic and zoonotic coronaviruses

Comparative therapeutic efficacy of remdesivir and combination lopinavir, ritonavir, and interferon beta against MERS-CoV

The antiviral compound remdesivir potently inhibits RNA-dependent RNA polymerase from Middle East respiratory syndrome coronavirus

Remdesivir is a direct-acting antiviral that inhibits RNAdependent RNA polymerase from severe acute respiratory syndrome coronavirus 2 with high potency

Rapid incorporation of Favipiravir by the fast and permissive viral RNA polymerase complex results in SARS-CoV-2 lethal mutagenesis

Favipiravir at high doses has potent antiviral activity in SARS-CoV

−infected hamsters, whereas hydroxychloroquine lacks activity

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

Ribavirin, Remdesivir, Sofosbuvir, Galidesivir, and Tenofovir against SARS-CoV-2 RNA dependent RNA polymerase (RdRp): A molecular docking study

Anti-HCV, nucleotide inhibitors, repurposing against COVID-19

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

Remdesivir-bound and ligand-free simulations reveal the probable mechanism of inhibiting the RNA dependent RNA polymerase of severe acute respiratory syndrome coronavirus 2

Revealing the Inhibition Mechanism of RNA-Dependent RNA Polymerase (RdRp) of SARS-CoV-2 by Remdesivir and Nucleotide Analogues: A Molecular Dynamics Simulation Study

RNA-DEPENDENT RNA POLYMERASE FROM SARS-COV-2. MECHANISM OF REACTION AND INHIBITION BY REMDESIVIR

1′-Ribose cyano substitution allows Remdesivir to effectively inhibit nucleotide addition and proofreading during SARS-CoV-2 viral RNA replication

Dynamic properties of SARS-CoV and SARS-CoV-2 RNA-dependent RNA polymerases studied by molecular dynamics simulations

Replicapermutation molecular dynamics simulations of an amyloid-β(16-22) peptide and polyphenols

ColDock: Concentrated Ligand Docking with All-Atom Molecular Dynamics Simulation

Amyloid Fibril Disruption by Ultrasonic Cavitation: Nonequilibrium Molecular Dynamics Simulations

Conformational Change of Amyloid-β 40 in Association with Binding to GM1-Glycan Cluster OPEN

Comparative protein modelling by satisfaction of spatial restraints

UCSF Chimera -A visualization system for and analysis

Asparagine and glutamine: Using hydrogen atom contacts in the choice of sidechain amide orientation

ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB

Comparison of simple potential functions for simulating liquid water

A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model

Gaussian 16, revision C.01

Automatic atom type and bond type perception in molecular mechanical calculations

Rational design of particle mesh ewald compatible lennard-jones parameters for +2 metal cations in explicit solvent

Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations

Temperature and pressure denaturation of chignol Folding and unfolding simulation by multibaric-multithermal molecular dynamics method

Transformation between α-helix and β-sheet structures of one and two polyglutamine peptides in J o u r n a l P r e -p r o o f explicit water molecules by replica-exchange molecular dynamics simulations

Molecular dynamics simulations of amyloidβ(16-22) peptide aggregation at air-water interfaces

Structural and fluctuational difference between two ends of Aβ amyloid fibril: MD simulations predict only one end has open conformations

Physiological concentrations of purines and pyrimidines

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-nCoV) in vitro

Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems

Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes

PTRAJ and CPPTRAJ: software for processing and analysis of molecular synamics trajectory data

RNA-dependent RNA polymerase: Structure, mechanism, and drug discovery for COVID-19

Aqueous solutions of divalent chlorides: Ions hydration shell and water structure

Structural Basis for Helicase-Polymerase Coupling in the SARS-CoV-2 Replication-Transcription Complex

Structure of replicating SARS-CoV-2 polymerase

Architecture of a SARS-CoV-2 mini replication and transcription complex

Remdesivir is a delayed translocation inhibitor of SARS-CoV-2 replication

Magnesium-cationic Dummy Atom Molecules Enhance Representation of DNA Polymerase β in Molecular Dynamics Simulations : Improved Accuracy in Studies of Structural Features and Mutational Effects

Many-body effect determines the selectivity for Ca2+ and Mg2+ in proteins

S.T. performed simulations, analyzed the results, and wrote the paper. All authors discussed the results and revised the paper.