RNA-dependent RNA polymerase (RdRp), is an enzyme essential component in the RNA replication within the life cycle of the severely acute respiratory coronavirus-2 (SARS-CoV-2), causing the deadly respiratory induced sickness COVID-19. Remdesivir is a prodrug that has seen some success in inhibiting this enzyme, however there is still the pressing need for effective alternatives. In this study, we present the discovery of four non-nucleoside small molecules that bind favorably to RdRp over adenosine-triphosphate (ATP) and active-form remdesivir-triphosphate (RTP) using high-throughput virtual screening (HTVS) coupled with extensive (total 4800 ns) molecular dynamics (MD) simulations with using the ZINC compounds database against SARS-CoV-2 RdRp (PDB: 7BV2). We found that the simulations with both ATP and RTP remained stable for the duration of their trajectories, and it was revealed that the phosphate tail of RTP was stabilized by a positive amino acid pocket near the entry channel of RTP and magnesium ions containing residues K551, R553, R555 and K621. It was also found that residues D623, D760, and N691 further stabilized the ribose portion of RTP with U10 on the template RNA strand forming hydrogen pairs with the adenosine motif. Using these models of RdRp, we employed them to screen the ZINC database of ~17 million molecules. Using docking and drug properties scoring, we narrowed down our selection to fourteen candidates. These were subjected to 200 ns simulations each underwent free energy calculations. We identified four hit compounds from the ZINC database that have similar binding poses to RTP while possessing lower overall binding free energies, with ZINC097971592 having a binding free energy two times lower than RTP.
Research Article
Working Mechanism of Remdesivir and Discovering New Potential Inhibitors to SARS-CoV-2 RNA-Dependent RNA Polymerase (RdRp) using High-Throughput Virtual Screening and Molecular Dynamics Simulations .
https://doi.org/10.21203/rs.3.rs-1119894/v1
This work is licensed under a CC BY 4.0 License
You are reading this latest preprint version
SARS-CoV-2
RNA-dependent RNA Polymerase
Remdesivir
Virtual screening
Molecular dynamics simulation
The highly contagious coronavirus induced illness COVID-19 is now a global threat to human well-being, with the current case numbers rising near 220 million and over 4.5 million deaths as of September 2021. With the increase in rollouts of vaccinations around the world, there are now some signs of relief to even the most affected regions. However, there is still the potential for new vaccine resistant strains of the coronavirus are evolving often as seen in the recent deadly outbreak of the new B.1.617.2 “Delta” variant in Nepal [1]. Additionally, there is the possibility for some vaccinated individuals to not develop an adequate immune response to create neutralizing antibodies. Therefore, there is still an urgent need for effective drug therapy in order to slow down the progression of sickness in affected patients. The SARS-CoV-2 virus is a positive-strand RNA virus whose genomic structure expresses high similarity with previous severe acute respiratory syndrome coronavirus (SARS-CoV). Like other RNA viruses, the replication of SARS-CoV-2 requires RNA-dependent RNA polymerase (nsp 12), a key enzyme which regulates the viral genome replication and transcription within the viral life cycle. It is thus a key validated target for the development of COVID-19 disease small molecule therapeutics.
Currently, Gilead’s Remdesivir (RDV) otherwise known as Veklury® has been approved by the FDA as an antiviral prodrug which targets the nucleotide addition process in RdRp. This prodrug is bio transformed in the body into the active-form remdesivir-triphosphate (RTP), a nucleoside mimic to adenosine triphosphate (ATP), which then undergoes nucleotide addition onto the growing RNA chain [2]. It is theorized that RDV inhibits RdRp by chain termination during the RNA translocation step [3]. Its action on treating SARS-CoV-2 infection has seen some success in reducing the time needed for recovery in very sick patients that require hospitalization [4]. However, this drug's efficacy is also controversial, as its administration was recently halted due to the lack of significant improvement expected in COVID-19 infected patients. Therefore, because RDV is the only FDA approved drug for treating COVID-19 patients, there is interest in the discovery of alternative treatments. Additionally, there is interest in developing non-nucleoside mimetics with high RdRp affinity for this purpose, as these are thought to circumvent the polymerization reaction.
There is a remarkable effort underway from researchers around the globe to help shed more light on potential alternatives for the inhibition of RdRp. There have been several studies (Table 1) [3, 5–42] that have reported the discovery of potential small molecule inhibitors of RdRp, mostly deriving from the in-silico approach coupled with virtual drug screening of numerous compound libraries. As of this study, few in-silico studies have considered the approach of using the massive ZINC library of approximately 17 million compounds using high-throughput virtual screening workflow (Figure 1) paired with long, multi-system molecular dynamics (MD) simulations [43]. A previous study reported the discovery of two potential isoformic RdRp inhibitors: ZINC09128258 and ZINC09883305, however these were observed to have lower docked binding affinity than the compared standards [21].
In our present study, we have utilized long MD simulations (1.0 µs) each of SARS-CoV-2 RdRp in complex with free form active RTP and ATP and template-primer RNA starting from the recently solved crystal complex (PDBID: 7BV2) [44]. Using the most abundant conformations of RTP complex from 1.0 µs simulations, we performed pharmacophore-based screening of the ZINC 15 library against the binding site of RdRp, leading to fourteen compounds with top docking scores [43]. From these top fourteen compounds, we performed classical MD simulations (200 ns each) to observe the stability for each complex system. Lastly, we performed Molecular Mechanism-General Born Surface Area (MM-GBSA) energy predictions to observe their relative binding affinities to RdRp after simulation [39, 45–47]. Four ZINC compounds showed stable binding to the RdRp-RNA complex. The structures of these top four compounds are shown in Figure 3. Additionally, these four systems showed improved MM-GBSA energies both the RTP and ATP complex systems as well as in reference to their initial docked complexes. This discovery helps introduce crucial knowledge of small molecule binding to RdRp in addition to the effort in discovering and development of anti-SARS-CoV-2 RdRp agents using in-silico methods.
2.1. Structure Building and Ligand Preparation
The recent electron microscopy solved crystal structure of SARS-COV-2 nsp12-nsp7-nsp8 RdRp in complex with template-primer RNA and covalently bound RDV named PDB ID 7BV2 was retrieved from the RCSB protein database [44, 48]. This structure was then prepared with the default charge state being optimized at pH=7 using Maestro’s built-in Protein Preparation Wizard [49–51]. The complex was then preprocessed, optimized, and minimized using the Protein Preparation Wizard within Maestro software [49]. The ligand RTP and ATP were then prepared for docking by generating the ionization/tautomeric states at pH=7±2 using Maestro’s Epik tool [49, 52, 53]. Using the prepared RdRp structure, a receptor grid file was generated around covalently the bound RDV position. RTP along with ATP were then re-docked into the receptor grid box using GlideXP scoring [50, 51].
2.2. MD System Setup
Multiple MD simulations were performed: two 1000 ns simulations for both RTP and ATP containing RdRp complexes. The selected compounds from HTVS amounted to an additional 2800 ns simulation or 200 ns for each of the 14 complexes [49, 54]. The MD simulations were performed applying the OPLS2005 force field for all protein, RNA and ligands. When creating the systems for RDV-TP, ATP and the 14 ZINC compounds, a salt concentration of 0.15 M NaCl was added to neutralize each system [55, 56]. There were then placed into an orthorhombic water box with a buffer distance of 10 Å using the SPC water model [57].
2.3. Relaxation and Production Runs
Using the Desmond module present in Maestro software, each system was relaxed using the default relaxation protocol for non-membrane proteins [58]. There are eight stages which constitute the protocol: 1). Minimization of the system with a restraint of solute heavy atoms; 2). Minimization without any restraints; 3). Simulation with heating from 0 K to 300 K, H2O barrier and gradual restraining; 4). Simulation under the NPT ensemble which is constant number of particles, pressure of 1 bar and constant temperature of 300K with an H2O barrier and heavy atoms being restrained: 5). Simulation under the NPT ensemble with equilibration solvent; 6). Simulation under the NPT ensemble with protein heavy atoms annealing from 10.0 kcal/mol to 2.0 kcal/mol; 7). Simulation under the NPT ensemble with Cα atoms restrained at 2.0 kcal/mol; and 8). Simulation of 1.5 ns under the NPT ensemble with no restraints. After this relaxation procedure, two 1000 ns production runs for RTP and ATP systems and fourteen 200 ns production runs were conducted under the NPT ensemble using the default protocol.
2.4. Conformational Clustering of RdRp Complexes
The Desmond trajectory clustering tool was used to group the complex structures from the last 100 ns of simulation for the ATP and RTP systems [57]. The merging distance cutoff was set to be 2.5 Å. The centroid representative (i.e. the structure having the largest number of neighbors in the structural family) was used to display the structural family. The most abundant conformations are identified as the cluster with the most occupancy. This was performed for the complexes containing ATP and RTP. Clustering was then performed on the top fourteen ZINC compounds selected from HTVS after 200 ns simulation.
2.5. Pharmacophore Screening and Hit Prioritization
The pharmacophore screening was performed using the active Remdesivir against ZINC product molecules using LigandScout 4.3 [54]. In our model, we determined the seven pharmacophore features to include the two negative ionizable groups, one aromatic ring, and four hydrogen bonds acceptors with their corresponding positions. Before further experiments, we determined the ionizable/tautomeric states of each compound under pH=7±2 using Maestro’s Epik tool [52, 53]. The lowest ionization/tautomeric states were selected, and their bond geometries were minimized to the lowest energetic state. We then employed our prepared RdRp complex to generate a receptor grid file around the RTP binding site, or i site, and the ionization, tautomeric and bond states of ATP were generated at pH=7±2 using Maestro’s Epik tool [52, 53]. This receptor grid file was additionally used for our docking of the ZINC database compounds using the Glide XP docking scoring function within Maestro software [49, 54]. High-Throughput Virtual Screening (HTVS) was also performed using the Glide module on the most abundant cluster representative of the RTP system after simulation utilizing the ZINC15 compound database [43, 50, 51]. This process consists of screening at three levels with gradually increasing the computational cost. The first level is based on ligand ADME/Toxicity prediction using QikProp of Schrodinger and docking using Glide [52, 53]. The second is performed using MD simulations of the docked complexes with explicit solvents (water) [58]. Lastly, the third is determined by the binding energy calculations. After HTVS, we narrowed down the library of compounds to the top fourteen based on their superior Glide scores and #star value [50, 51]. We then subjected these fourteen to 200 ns simulations and performed MM-GBSA calculations to determine their relative binding affinities after the simulation [39, 45]. From the fourteen compounds, four compounds were selected based on their favorable binding energies.
2.6. Confirmation of Stability of Hits and MM-GBSA Calculation
The Desmond SID module was employed to analyze the interaction between proteins and ligands as well as to determine their stability in each MD simulation [58]. This includes the Root-Mean Square Deviation (RMSD) of the protein-ligand complexes and Root-Mean Square Deviation (RMSF) of both the protein and ligands. Additionally, secondary structure elements (SSE) were explored. Lastly, the protein-ligand contacts containing H-bond, ionic, hydrophobic interactions were explored.
In order to check the convergence of the simulations, we investigated the root mean square deviation (RMSD) values of protein Cɑ, RNA O5’ atoms, and ligand main atoms for each trajectory using the SID tool within Maestro software [58]. The root mean square fluctuation (RMSF) of each individual amino acid, RNA and ligand atom were also calculated to characterize the local movement of individual structural components of each complex.
The Molecular Mechanism-General Born Surface Area (MM-GBSA) binding energies were calculated on the frames for the whole duration of both systems as well as for each ZINC system [45]. The OPLS3 force field, VSGB 2.0 solvation model and the default Prime procedure was used for the MM-GBSA calculations [46, 47, 55, 56, 59, 60]. The default procedure consists of three steps: Receptor alone (minimization), Ligand alone (minimization), Receptor-ligand complex (minimization). The total binding free energy equation is: ΔGBind = Ecomplex (minimized) - (Eligand (minimized) + Ereceptor(minimized) ). To gain a more detailed understanding of binding nature, the original interaction terms (Coulombic + H-bond + GB solvation+ van der Waals + π-π packing + self-contact + lipophilic) were merged into three components: Eelectrostatics, EvdW, and Elipophilic, where Eelectrostatics =(Hbond + Ecoulomb +EGBsolvation), EvdW = (EvdW+ Eπ−π +Eself−contact) and Elipophilic.
3.1. The Binding of Remdesivir Triphosphate and ATP to RdRp
Adenosine triphosphate was docked against the crystal structure of SARS-CoV-2 RdRp (PDBID: 7BV2) at the i-binding site which formed hydrogen pairing with the template strand of 14 nucleotide bases and the primer strand of 11 nucleotide bases [44]. The conformational pose of the docked RTP compared to the covalently bound RDV in the crystal structure was nearly identical (Figure S2). We then generated a docking grid from the docking pose of RTP to dock ATP into the same catalytic site of RdRp. The conformational pose of ATP compared to the covalently bound crystal RDV was nearly identical. It was also seen that the triphosphate tail of RTP was positioned near one of the Mg+ ions. These conformations were then used to perform long (1.0 µs) simulations for each complex to further elucidate more about the relaxed conformations of each and to obtain a stable conformation for our high-throughput virtual screening portion of the study.
Table 1. Survey of the current studies on the discovery and/or repurposing of approved drugs for the inhibition of SARs-COV-2 RdRp.
The trajectory convergence of the simulation for both systems was confirmed through analyzing the RMSD during the 1000 ns (1.0 µs) simulations shown in Figure 2A and Figure 2B. The trajectory RMSD is of the protein receptor backbone Cɑ atoms and the ligand main atoms for both the ATP and RTP complexes. It is clearly seen in Figure 2A that the RTP system reached stability very early in the simulation (~50 ns) while the ATP system in Figure 2B reached stability around 550 ns. Importantly, both protein complexes maintained average RMSD values of approximately 2.0 Å. This is indicative that the differences between ATP and RTP binding had no significant change in the overall RdRp conformation through the 1.0 µs simulation.
The ligand RMSD average values for the ATP system are within 1 Å of the protein RMSD values and reached stability at approximately 550 ns. The ligand RMSD for the RTP system also matched the average distance as in the protein with a convergence very early in the simulation; approximately 50 ns. It is clearly seen that the complex containing ATP possessed higher RMSD values compared to the RTP complex indicating that the binding of RTP was more stable during the simulation.
We additionally wanted to confirm the stability of the template-primer RNA in the complex through simulation. The RMSD of the RNA main atoms is shown in Figure S3, both ATP and RTP systems show clear stability and early convergence at around 50 ns with values both approximately 2.0 Å which is also consistent with the protein Cɑ RMSD. Likewise seen in the ligand main atoms RMSD plot, the complex containing RTP showed slightly lower average values compared to the ATP system indicating more stable binding occurred during the simulation. To observe the fluctuation of each individual amino acid in the protein structure, we recorded the RMSF of the RdRp complexes during the simulation.
Ligand main atoms RMSF (Figure S5) was also recorded to observe the fluctuation of the individual atoms of each ligand in order to better understand the regions of stability of ATP and RTP. Clearly, RTP in complex saw less fluctuation per atom over the simulation as compared to ATP. Short peaks at atoms 1-4 indicate the slight instability seen in the gamma phosphate on the phosphate tail which is not abnormal. Interestingly, small peaks in atoms 16 and 30 indicate a slight increase in fluctuation seen at the 3’ and 4’ hydroxyl groups on the ribose ring of RTP. This specificity at only these two polar functional groups is also seen in atoms 29 and 30 of the ATP system and likely indicates positioning for nucleophilic attack during nucleotide addition.
3.2. Clustering Analysis
We then performed cluster analysis of both MD systems after the simulation to generate the most abundant conformational pose during the 1000 ns duration. In the most abundant conformational pose of the simulation (Figure 5), it was observed that the adenine motif of both ATP and RTP maintained H-Bond pairing with U10 on the template strand and partial pi-pi stacking of A11 on the primer strand. The RTP system was more able to sustain salt bridge interactions with positive residues R555, R553, R624, K621. The RTP system also saw interaction with S759 on the nitrile group while Y619 was seen to form H-Bonds with negative oxygen on the gamma phosphate of ATP. In addition, the RTP system also saw stabilization of the 3’ hydroxyl group on the ribose ring by H-Bonds with residues N691 and T680.
During the 1.0 µs simulation of both the ATP and RTP systems, several amino acid contacts were observed to be maintained throughout with high occupancy (Figures S7). For both the ATP and RTP systems, there was significant occupancy in the H-Bond interactions between the triphosphate tails and a positive region of amino acids such as K551, K621, K798, R553 and R555 of ≥ 30% of the simulation time. For the RTP complex, residues K621, R553 and R555 formed H-bonds with negatively charged oxygen on alpha, beta and gamma phosphate with all possessing an occupancy of above 90%, indicating that these contacts could be crucial for the stability and positioning of RTP in the catalytic site of RdRp. These interactions were seen in previous study where both ATP and RTP were docked to RdRp [61]. Both systems also maintained the maximum occupancy interaction with Mg+ ions with the triphosphate tails, suggesting that at least one ion is responsible for the stabilization and positioning of incoming nucleotides. There was modest occupancy of hydrogen bond formation with hydroxyl groups located on the ribose ring portions with residues D623, D760, and N691, indicating secondary stabilizing roles for these amino acids potentially related to the initiation of nucleotide addition. Negative residue D623 saw an average occupancy of 56% bonding with both hydroxyl groups of the ribose ring portion of RTP. In the ATP complex, D760 saw 34% occupancy with the 3’ hydroxyl group located on the ribose ring which is involved in nucleotide addition.
In Figures S6 are the histogram of protein-ligand contacts occupancy during the simulations for ATP and RTP complexes. These show the smaller contacts seen with occurrence of ≤ 30% of the simulation as well as contacts that occurred at over 100% which indicates multiple ligand atoms contacting same amino acid. In the ATP complex, the ligand contacted several more protein residues with brief occupancy compared to the RTP complex which is likely due to lower relative ligand bound stability seen in the former. Also, more water bridge interactions occurred with the ATP complex indicating weaker overall affinity compared to the RTP complex. Furthermore, RTP saw higher occupancy of H-bonds compared to others confirming the relative favorability of bonding to RdRp.
Previous studies have also reported similar findings on the critical amino acids involved in RDV binding. Koulgi et al. described an ensemble approach of the free form binding of prodrug RDV in complex with SARS-CoV-2 RdRp contacted amino acids Y451, T540, M542, R548, K551, R553, R555, A558, D618, S674, D761 and E811 [62]. Amino acid R555 was observed to form H-bonds with the alpha phosphate group in RDV in all five ensemble representatives [62]. Other studies have aimed to observe the active triphosphate form interaction upon RdRp binding. Zhang et al. found that the nitrile group was positioned in a pocket including K545, Y546 and A547 with the phosphate tail interacting with K551, R553 and R555. Additionally, positive amino acids near the palm subdomain were found to be crucial for RTP activity [63].
3.3. High-Throughput Virtual Screening finds 14 compounds.
Once we validated that our simulations with ATP and RTP bound to RdRp met convergence during the simulation, we then used these complexes to perform our initial screening. Firstly, pharmacophoric screening of the ZINC database was performed by selecting seven pharmacophoric features, i.e. two negative ionizable, one aromatic ring, and four H-bond acceptors of RTP. Shown in Table 2 is the results from HTVS docking against RdRp to retrieve the top fourteen hits with the highest GlideXP scores and lowest #stars as the determining factors. The GlideXP scores measured in kcal/mol are tabulated based on decreasing favorability. A more negative GlideXP score typically indicates a good initial binding affinity to RdRp. The initial conformational pose of each ligand in complex with RdRp and RNA as well as their detailed 2D ligand contacts is shown in Table S2. These fourteen compounds were then subjected to classical 200 ns MD simulations and further MM-GBSA energy calculations for each.
|
Molecule |
Docking score (kcal/mol) |
# Stars |
|---|---|---|
|
ZINC000014651456 |
-13.5 |
0 |
|
ZINC000257306096 |
-12.9 |
0 |
|
ZINC000238950253 |
-12.8 |
0 |
|
ZINC000299798705 |
-12.7 |
0 |
|
ZINC000067790716 |
-12.4 |
0 |
|
ZINC000089920955 |
-11.9 |
0 |
|
ZINC000097971592 |
-11.8 |
1 |
|
ZINC000065742965 |
-11.7 |
0 |
|
ZINC000016040970 |
-11.6 |
0 |
|
ZINC000408592119 |
-11.6 |
0 |
|
ZINC000237948681 |
-11.5 |
0 |
|
ZINC000069492350 |
-11.4 |
0 |
|
ZINC000002146610 |
-11.2 |
0 |
|
ZINC000084651559 |
-11.2 |
0 |
| * Docking score is empirically calculated in kcal/mol from rigid receptor GlideHTVS protocols which help define compounds with good binding affinity. The scoring function is comprised of lipophilic, hydrogen bonding, hydrophobic terms as well as a rotable bond penalty. The #stars is a parameter defined by QikProp module and scores compounds based on their similarity to known medicines. | ||
The MM-GBSA energies and average protein-ligand RMSD of all fourteen MD simulations are summarized in Table S3. Within Figures 6A and 6B show the ligand RMSD and RMSF values of all fourteen ZINC candidates with ATP and RTP systems calculated over the last 100 ns of simulation. The free energy of binding gives a more accurate insight into both the simulation efficacy and flexible binding affinity of ligands to their target as opposed to rigid ligand docking. The summary of the 200 ns simulation runs and MM-GBSA energy calculation in Table S3 shows that several compounds possessed more favorable free energy of binding compared to the ATP and RTP systems. From this subset of ligands, we selected the top four with the most favorable MM-GBSA calculated energies for further convergence analysis to ensure that the conformational stability was reached and maintained. Figures 2C-F shows the RMSD of protein Cα atoms and ligand main atoms over the entire duration for these selected compounds. The protein portions reached convergence early in the simulation for each system containing ZINC097971592, ZINC002146610, ZINC069492350 and ZINC408592119. It was observed that the RMSD values of the protein Cα of the four ZINC systems were lower than that of the RTP system. Also, the ligand RMSD values for ZINC097971592 and ZINC408592119 were higher than that of the RTP complex.
We additionally wanted to observe the RMSF of the RNA main atoms during simulation to gain some more insight into the individual fluctuation of each nucleotide in the RNA sequence. In Figure 4A, this displays the template-primer RNA RMSF obtained from the O5’ atoms on the RNA backbone for ATP and RTP systems and the four ZINC systems with Figure 4B indicating the nucleotide sequence indexing used during the recording. The RMSF values were high at nucleotides G1-U3 on the 5’ end of the primer strand as well as high peaks around A22-C25 towards the 3’ end of the template strand for all 6 systems. This indicates high fluctuation induced from the exposure to solvent. In the middle of the plot is a slight peak around U12 on the template strand indicating that this portion did not interact with either ligand. Nucleotides U13-A15 on the template strand show a sharp decrease in fluctuation which is most likely due to the binding and hydrogen pairing of both ATP and RTP.
The average RMSD values over the last 50 ns of simulation for receptor and ligand is also summarized in Table 3. The four ZINC systems maintained similar average receptor RMSD values were approximately 3.3 Å which is near the value of the ATP system of 2.4 Å. This result suggests the RdRp complex system remained in similar conformational state regardless of their structure dissimilarity against ATP. Table 3 is also a breakdown of the individual MM-GBSA terms such as the Van der Waals, electrophilic and lipophilic energies of the top four selected ZINC compounds as well as ATP and RTP for comparison. The summation of these term equates to the relative free energy of binding measured in kcal/mol. Clearly, the four selected compounds from the ZINC database resulted in lower overall MM-GBSA values over the ATP and RTP systems of -35.9 ± 3.1 kcal/mol and -21.3 ± 5.9 kcal/mol, respectively. The RdRp complex containing ZINC097971592 had the most negative MM-GBSA value of -54.5 ± 13.4 kcal/mol, with ZINC002146610, ZINC069492350 and ZINC408592119 having values of -47.5 ± 3.4, -43.9 ± 3.8 and -37.6 ± 4.6 kcal/mol, respectively. The four ZINC systems maintained similar average receptor RMSD values were approximately 3.3 Å which is near the value of the ATP system of 2.4 Å. This result suggests the RdRp complex system remained in similar conformational state regardless of their structure dissimilarity against ATP.
|
Ligand |
ΔGVDW |
ΔGELE |
ΔGLIP |
MM-GBSA (kcal/mol) |
Receptor RMSD (Å) |
Ligand RMSD (Å) |
|---|---|---|---|---|---|---|
|
ATP |
-14.4 ± 5.0 |
-15.3 ± 4.4 |
-6.2 ± 0.7 |
-35.9 ± 3.1 |
2.4 |
4.5 |
|
RTP |
-19.1 ± 4.1 |
6.9 ± 6.2 |
-9.1 ± 1.4 |
-21.3 ± 5.9 |
1.9 |
0.9 |
|
ZINC097971592 |
-42.3 ± 5.4 |
6.2 ± 15.9 |
-18.3 ± 3.1 |
-54.5 ± 13.4 |
3.9 |
4 |
|
ZINC002146610 |
-52.8 ± 3.0 |
23.1 ± 3.0 |
-17.8 ± 1.3 |
-47.5 ± 3.4 |
3 |
1.5 |
|
ZINC069492350 |
-42.2 ± 3.1 |
14.4 ± 3.3 |
-16.1 ± 2.1 |
-43.9 ± 3.8 |
3.1 |
2.8 |
|
ZINC408592119 |
-55.6 ± 3.1 |
35.9 ± 3.7 |
-17.9 ± 1.9 |
-37.6 ± 4.6 |
3.3 |
4.4 |
| *VDW energies was calculated with the summation of ligand ΔGvdW energies, pi-pi packing correction energies and the self-contact correction energies. | ||||||
| *Electrophilic terms are calculated by the summation of coulombic, hydrogen-bonding and generalized born electrostatic solvation energies. Hydrophilic term comprised of lipophilic energies. Total MM-GBSA calculated in kcal/mol is summation of ΔGvdW, ΔGele and ΔGhyd terms. | ||||||
The detailed 2D receptor-RNA-ligand contacts diagrams of the most abundant conformational pose for each of the top four ZINC candidates are shown in Figure 7A-L. The 2d ligand-protein interactions diagram for the most abundant conformational pose of ZINC097971592, or estriol 3 glucoronide, in complex with RdRp is shown in Figure 7C. There were several residue contacts conserved from the ATP and RTP pose. Positive residues R624, R555, and R553 formed H-bonds with mainly the polar negative carboxylate portion of the molecule, indicating that this compound is situated with the polar end in the positive catalytic region of RdRp. Additionally, ZINC097971592 also came into H-Bond contact with S682 on the hydroxyl groups on the aromatic portion. The ligand-protein contacts summary diagrams for these four ZINC candidates are shown in Table S7. In this, several of the residues seen in contact with the most abundant pose were also maintained with large occupancy through the simulation. Residues R555 and R624 established interactions with ZINC097971592 with 89% and 74% occupancy, respectively. Interestingly, K798 formed a salt-bridge with the negative oxygen with a relatively high occupancy value of 85%. Contacts with lesser occupancy were S682 with a 44% occupancy and R551, R553 of 50% or less. This compound’s enantiomers were identified in a previous VirtualFlow in that there was activity towards plpro, nucleoprotein, nsp10, nsp14 and orf7a in addition to nsp12 [64].
In Figure 7F, the psoralen derivative ZINC002146610 also shows some similarity with the RTP pose and with the other ZINC hits in that the polar carboxylate group is facing inwards towards the positive amino acid region where triphosphate tail usually lays. Also, this compound formed pi-pi stacking interactions on the psoralen ring with nucleotide U20 on the primer strand and A11 of the template strand of RNA. Additionally, the polar end formed salt bridge and H-Bond interactions with R555 and the nearby Mg+ ion. The ligand-protein contacts summary for ZINC002146610 shows similar contacts seen in the most abundant pose. Positive amino acid R555 maintained an occupancy of 113% with the polar region of the molecule, meaning there were additional contacts seen with residue in addition to the full simulation time. Like the most abundant pose, the magnesium ion maintained full contact with the polar end via salt-bridge interactions which was seen to be further stabilized by D618. Additional stabilizing water-bridges were established with occupancies of 47% and 48% with D623 and R555, respectively. Interestingly, this compound as well as its derivatives have seen previous attention in being low toxicity potential selective and reversable inhibitors of (B5i) or the chymotrypsin-like subunit of human immunoproteasome which is associated with the treatment of autoimmune diseases and various types of cancer [65].
Table 4. Ligand-protein contacts during simulation from ATP, RTP and the four ZINC compound systems.
*Ligand-protein contacts may consist of hydrogen bonding, lipophilic, hydrophobic and electrostatic interactions. Contacts shown are present in at least 30% of trajectory time. Green highlight are residues that interacted with at least 3 of the six tested ligands.
In Figure 7I, conserved residue R555 is seen in contact with the polar end of ZINC69492350 in similar fashion to the previous ZINC compounds in addition to the RTP system. Amino acid R555 has been seen to be crucial in the activity of small molecule binding to RdRp in this study as well as others [7–11, 31, 39–41, 62, 66, 67]. In this binding pose of ZINC69492350, R555 forms a salt bridge to the negatively charged oxygen on the carboxyl group, with an additional hydrogen bond contact with the partially negative doubly bonded oxygen on the carboxyl end. Likewise, the aromatic portion of the molecule faces in direction with the hydrogen pairing nucleotide U10 on the template strand RNA. Surrounding amino acids are also conserved which includes D760, D623, S682 and K545. Lastly, the magnesium ions are seen to stabilize ZINC69492350 with pi-cation interactions with the thiophene portion. In the ligand-protein contacts summary diagram in the supporting information, R555 maintained full occupancy with the addition of other interactions less seen with a value of 150% with the polar carboxyl end of the molecule. Like the other selected ZINC compounds, D623 formed stabilizing water-bridge interactions with the same polar end with a value of 33%. Polar residue T687 contacted ZINC69492350 55% of the simulation, unlike the other compounds where it served as a spectator.
In Figure 7L, ZINC408592119 lies in a similar binding pose to ATP and RTP. According to the amino acid contacts seen in the most abundant conformational pose of the simulation. Positive residues K551, K545, K621 and R555 were all seen to form hydrogen bonds and pi-cation interactions with partially positive groups on the polar ends of the molecule. Negative residues such as D618 and D760 are seen to form hydrogen bonding with partially positive portions such as on the imidizolidine group. The significance of D618 is consistent with previous studies [39–41, 62, 66, 67]. In the ligand-protein contacts summary diagram in the supporting information, positive residue K551 maintained an occupancy of 34% in the same fashion as the most abundant pose. Magnesium also contributed full occupancy with the polar carboxyl end of the molecule consistent in Table S7 and other ZINC candidates. Lastly, in addition to forming hydrogen bonds with the imidazolidine portion, also formed a water-bridge 76% of the simulation. Table 4 is the receptor-ligand contacts with an occurrence of greater than 30% tabulated for the ATP, RTP and four ZINC hit complexes. Clearly, positive residues K621, R555, R553 and K551 all induced contact with more than two-thirds of the ligands selected, with R555 meeting all 5 ligands. Likewise, negative residues D623, D798 and D760 saw interaction with at least half of the ligands. Lastly, hydrophobic residue N691 interacted with all the selected ligands apart from ZINC408592119.
The absorption, distribution, metabolism and excretion (ADME) properties were also explored for these four compounds in addition to ATP and prodrug RTP using the SwissADME webserver [68] (Table 5). Within the four non-nucleoside analogs, three of them were reported to have high predicted GI absorption and cytochrome P450 inhibition for CYP1A2, CYP2C19, CYP2C9 and CYP3A4. With these compounds being very early drug hits selected using in silico methods, obviously there is the need for experimental evaluation needed to evaluate their biological attributes. Additionally, it is still unclear how these compounds would behave and be metabolized when introduced to a biological system. Based on the predicted bioavailability, three of the four compounds possess high (56%) availability in rats.
Table 5. Summary of drug property predictions of ATP, RTP and top four ZINC candidates obtained from the SwissADME webserver.
In this study, we aimed to determine the binding of the active-form of remdesivir to SARS-CoV-2 RdRp in complex with template-primer RNA and magnesium ions. Using the model constructed from PDB 7BV2 [44], we then constructed a complex with ATP for further validation. We confirmed that these models were stable during the 1.0 µs MD simulation by observing the receptor/ligand RMSD and RMSF values. We then performed HTVS of 1179 compounds from the ZINC database against the most abundant conformation obtained from the RTP system. Based on their best docking score and low #stars, we selected the best fourteen compounds form this list. We then subjected these fourteen to an additional 200 ns relaxation MD simulation as well as MM-GBSA binding energy calculations to further elucidate more about their binding behavior and comparability to the classical nucleotides ATP and RTP. From these results, we ultimately selected the four ZINC candidates: ZINC097971592, ZINC002146610, ZINC69492350 and ZINC408592119 based on their lower MM-GBSA binding energy values than that of ATP and RTP. We also observed the detailed ligand-receptor contacts for each ZINC compounds pre and post simulation to compare their similarity in conformational pose to ATP and RTP. It was observed that R551, R553 and R555 have significant importance relating to negative motif hydrogen bond stabilization of both nucleotide phosphate tails and the ZINC polar ends. Additionally, negative amino acids D618 and D760 were seen to establish metal coordination bridges with Mg+ ions to further stabilization, as well as other water bridges and hydrogen bonding interactions with nearly full occupancies. These discoveries suggest that the four selected ligands from HTVS of the ZINC database could display similar or better binding affinity and inhibition activities of towards RdRp in-vitro compared to traditional nucleoside analogues.
Data Availability statements
All data generated or analysed during this study are included in this published article and its supplementary information document.
Acknowledgements
C.W acknowledges the support by the New Jersey Health Foundation (PC 76-21) and the National Science Foundation under Grants NSF ACI-1429467/RUI-1904797, and XSEDE MCB 170088. The Anton2 machine at the Pittsburgh Supercomputing Center (PSCA17017P) was generously made available by D. E. Shaw Research.
- Lopez Bernal J, Andrews N, Gower C, Gallagher E, Simmons R, Thelwall S, Stowe J, Tessier E, Groves N, Dabrera G, Myers R, Campbell CNJ, Amirthalingam G, Edmunds M, Zambon M, Brown KE, Hopkins S, Chand M, Ramsay M (2021) N Engl J Med 385(7):585
- Kokic G, Hillen HS, Tegunov D, Dienemann C, Seitz F, Schmitzova J, Farnung L, Siewert A, Hobartner C, Cramer P (2021) Nat Commun 12(1):279
- Jockusch S, Tao CJ, Li XX, Chien MC, Kumar S, Morozova I, Kalachikov S, Russo JJ, Ju JY (2020) Scientific Reports 10(1)
- Beigel JH, Tomashek KM, Dodd LE, Mehta AK, Zingman BS, Kalil AC, Hohmann E, Chu HY, Luetkemeyer A, Kline S, Lopez de Castilla D, Finberg RW, Dierberg K, Tapson V, Hsieh L, Patterson TF, Paredes R, Sweeney DA, Short WR, Touloumi G, Lye DC, Ohmagari N, Oh MD, Ruiz-Palacios GM, Benfield T, Fatkenheuer G, Kortepeter MG, Atmar RL, Creech CB, Lundgren J, Babiker AG, Pett S, Neaton JD, Burgess TH, Bonnett T, Green M, Makowski M, Osinusi A, Nayak S, Lane HC, Members A-SG (2020) N Engl J Med 383(19):1813
- Chien MC, Anderson TK, Jockusch S, Tao CJ, Li XX, Kumar S, Russo JJ, Kirchdoerfer RN, Ju JY (2020) Journal of Proteome Research 19(11):4690
- Ju JY, Li XX, Kumar S, Jockusch S, Chien MC, Tao CJ, Morozova I, Kalachikov S, Kirchdoerfer RN, Russo JJ (2020) Pharmacology Research & Perspectives 8(6)
- Aftab SO, Ghouri MZ, Masood MU, Haider Z, Khan Z, Ahmad A, Munawar N (2020) Journal of Translational Medicine 18(1)
- Ahmad J, Ikram S, Ahmad F, Rehman IU, Mushtaq M (2020) Heliyon 6(7)
- Aktas A, Tuzun B, Aslan R, Sayin K, Ataseven H (2020) Journal of Biomolecular Structure & Dynamics
- Alexpandi R, De Mesquita JF, Pandian SK, Ravi AV (2020) Frontiers in Microbiology 11
- Al-Masoudi NA, Elias RS, Saeed B (2020) Biointerface Research in Applied Chemistry 10(5):6444
- Ao S, Han D, Sun L, Wu YH, Liu S, Huang YJ (2020) Frontiers in Genetics 11
- Calligari P, Bobone S, Ricci G, Bocedi A (2020) Viruses-Basel 12(4)
- da Silva FMA, da Silva KPA, de Oliveira LPM, Costa EV, Koolen HHF, Pinheiro MLB, de Souza AQL, de Souza ADL (2020) Memorias Do Instituto Oswaldo Cruz 115
- da Silva JKR, Figueiredo PLB, Byler KG, Setzer WN (2020) International Journal of Molecular Sciences 21(10)
- Elfiky AA (2020) Life Sci 248:117477
- Elfiky AA (2020) Journal of Biomolecular Structure & Dynamics
- Indu P, Rameshkumar MR, Arunagirinathan N, Al-Dhabi NA, Arasu MV, Ignacimuthu S (2020) Journal of Infection and Public Health 13(12):1856
- Kandeel M, Kitade Y, Almubarak A (2020) Peerj 8
- Mohamed TA, Elshamy AI, Ibrahim MAA, Zellagui A, Moustafa MF, Abdelrahman AHM, Ohta S, Pare PW, Hegazy MEF (2020) Rsc Advances 10(57):34541
- Parvez MSA, Karim MA, Hasan M, Jaman J, Karim Z, Tahsin T, Hasan MN, Hosen MJ (2020) International Journal of Biological Macromolecules 163:1787
- Pokhrel R, Chapagain P, Siltberg-Liberles J (2020) Journal of Medical Microbiology 69(6):864
- Ruan ZJ, Liu C, Guo YT, He ZQ, Huang XH, Jia X, Yang T (2021) Journal of Medical Virology 93(1):389
- Tchesnokov EP, Gordon CJ, Woolner E, Kocinkova D, Perry JK, Feng JY, Porter DP, Gotte M (2020) Journal of Biological Chemistry 295(47):16156
- Wu CR, Liu Y, Yang YY, Zhang P, Zhong W, Wang YL, Wang QQ, Xu Y, Li MX, Li XZ, Zheng MZ, Chen LX, Li H (2020) Acta Pharmaceutica Sinica B 10(5):766
- Zhao Z, Bourne PE (2020) Journal of Proteome Research 19(11):4698
- Ahmed S, Mahtarin R, Ahmed SS, Akter S, Islam MS, Al Mamun A, Islam R, Hossain MN, Ali MA, Sultana MUC, Parves MR, Ullah MO, Halim MA (2020) Journal of Biomolecular Structure & Dynamics
- Barage S, Karthic A, Bavi R, Desai N, Kumar R, Kumar V, Lee KW (2020) Journal of Biomolecular Structure & Dynamics
- Borquaye LS, Gasu EN, Ampomah GB, Kyei LK, Amarh MA, Mensah CN, Nartey D, Commodore M, Adomako AK, Acheampong P, Mensah JO, Mormor DB, Aboagye CI (2020) Biomed Research International 2020
- Dwarka D, Agoni C, Mellem JJ, Soliman ME, Baijnath H (2020) South African Journal of Botany 133:273
- El Hassab MA, Shoun AA, Al-Rashood ST, Al-Warhi T, Eldehna WM (2020) Frontiers in Chemistry 8
- Gul S, Ozcan O, Asar S, Okyar A, Baris I, Kavakli IH (2020) Journal of Biomolecular Structure & Dynamics
- Gutierrez-Villagomez JM, Campos-Garcia T, Molina-Torres J, Lopez MG, Vazquez-Martinez J (2020) Journal of Physical Chemistry Letters 11(19):8008
- Kar P, Sharma NR, Singh B, Sen A, Roy A (2020) Journal of Biomolecular Structure & Dynamics
- Khan A, Khan M, Saleem S, Babar Z, Ali A, Khan AA, Sardar Z, Hamayun F, Ali SS, Wei DQ (2020) Interdisciplinary Sciences-Computational Life Sciences 12(3):335
- Mutlu O, Ugurel OM, Sariyer E, Ata O, Inci TG, Ugurel E, Kocer S, Turgut-Balik D (2020) Journal of Biomolecular Structure & Dynamics
- Narayanan N, Nair DT (2020) Iubmb Life 72(10):2112
- Ribaudo G, Ongaro A, Oselladore E, Zagotto G, Memo M, Gianoncelli A (2020) Journal of Biomolecular Structure & Dynamics
- Sharma A, Vora J, Patel D, Sinha S, Jha PC, Shrivastava N (2020) Journal of Biomolecular Structure & Dynamics
- Singh J, Malik D, Raina A (2020) Journal of Biomolecular Structure & Dynamics
- Singh S, Sk MF, Sonawane A, Kar P, Sadhukhan S (2020) Journal of Biomolecular Structure & Dynamics
- Elkarhat Z, Charoute H, Elkhattabi L, Barakat A, Rouba H (2020) Journal of Biomolecular Structure & Dynamics
- Sterling T, Irwin JJ (2015) J Chem Inf Model 55(11):2324
- Yin WC, Mao CY, Luan XD, Shen DD, Shen QY, Su HX, Wang XX, Zhou FL, Zhao WF, Gao MQ, Chang SH, Xie YC, Tian GH, Jiang HW, Tao SC, Shen JS, Jiang Y, Jiang HL, Xu YC, Zhang SY, Zhang Y, Xu HE (2020) Science 368(6498):1499
- Kollman PA, Massova I, Reyes C, Kuhn B, Huo S, Chong L, Lee M, Lee T, Duan Y, Wang W, Donini O, Cieplak P, Srinivasan J, Case DA, Cheatham TE (2000) Accounts of Chemical Research 33(12):889
- Jacobson MP, Friesner RA, Xiang Z, Honig B (2002) Journal of Molecular Biology 320(3):597
- Jacobson MP, Pincus DL, Rapp CS, Day TJF, Honig B, Shaw DE, Friesner RA (2004) Proteins: Structure, Function, and Bioinformatics 55(2):351
- Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, Shindyalov IN, Bourne PE (2000) Nucleic acids research 28(1):235
- Madhavi Sastry G, Adzhigirey M, Day T, Annabhimoju R, Sherman W (2013) Journal of Computer-Aided Molecular Design 27(3):221
- Friesner RA, Banks JL, Murphy RB, Halgren TA, Klicic JJ, Mainz DT, Repasky MP, Knoll EH, Shelley M, Perry JK, Shaw DE, Francis P, Shenkin PS (2004) Journal of Medicinal Chemistry 47(7):1739
- Friesner RA, Murphy RB, Repasky MP, Frye LL, Greenwood JR, Halgren TA, Sanschagrin PC, Mainz DT (2006) Journal of Medicinal Chemistry 49(21):6177
- Greenwood JR, Calkins D, Sullivan AP, Shelley JC (2010) Journal of Computer-Aided Molecular Design 24(6):591
- Shelley JC, Cholleti A, Frye LL, Greenwood JR, Timlin MR, Uchimaya M (2007) Journal of Computer-Aided Molecular Design 21(12):681
- Wolber G, Langer T (2005) Journal of Chemical Information and Modeling 45(1):160
- Harder E, Damm W, Maple J, Wu C, Reboul M, Xiang JY, Wang L, Lupyan D, Dahlgren MK, Knight JL, Kaus JW, Cerutti DS, Krilov G, Jorgensen WL, Abel R, Friesner RA (2016) Journal of Chemical Theory and Computation 12(1):281
- Shivakumar D, Harder E, Damm W, Friesner RA, Sherman W (2012) Journal of Chemical Theory and Computation 8(8):2553
- Mark P, Nilsson L (2001) The Journal of Physical Chemistry A 105(43):9954
- Bowers KJ, Chow E, Xu H, Dror RO, Eastwood MP, Gregersen BA, Klepeis JL, Kolossvary I, Moraes MA, Sacerdoti FD, Salmon JK, Shan Y, Shaw DE. Scalable algorithms for molecular dynamics simulations on commodity clusters. Proceedings of the 2006 ACM/IEEE conference on Supercomputing. Tampa, Florida: Association for Computing Machinery, 2006:84
- Li J, Abel R, Zhu K, Cao Y, Zhao S, Friesner RA (2011) Proteins: Structure, Function, and Bioinformatics 79(10):2794
- Jorgensen WL, Maxwell DS, Tirado-Rives J (1996) Journal of the American Chemical Society 118(45):11225
- Wu Y, Chang KY, Lou L, Edwards LG, Doma BK, Xie ZR (2020) Inform Med Unlocked 21:100461
- Koulgi S, Jani V, Uppuladinne MVN, Sonavane U, Joshi R (2020) Rsc Advances 10(45):26792
- Zhang L, Zhou R (2020) J Phys Chem B 124(32):6955
- Gorgulla C, Padmanabha Das KM, Leigh KE, Cespugli M, Fischer PD, Wang Z-F, Tesseyre G, Pandita S, Shnapir A, Calderaio A, Gechev M, Rose A, Lewis N, Hutcheson C, Yaffe E, Luxenburg R, Herce HD, Durmaz V, Halazonetis TD, Fackeldey K, Patten JJ, Chuprina A, Dziuba I, Plekhova A, Moroz Y, Radchenko D, Tarkhanova O, Yavnyuk I, Gruber C, Yust R, Payne D, Näär AM, Namchuk MN, Davey RA, Wagner G, Kinney J, Arthanari H (2021) iScience 24(2):102021
- Sosič I, Gobec M, Brus B, Knez D, Živec M, Konc J, Lešnik S, Ogrizek M, Obreza A, Žigon D, Janežič D, Mlinarič-Raščan I, Gobec S (2016) Angewandte Chemie International Edition 55(19):5745
- Arba M, Nur-Hidayat A, Usman I, Yanuar A, Wahyudi ST, Fleischer G, Brunt DJ, Wu C (2020) Indonesian Journal of Chemistry 20(6)
- Arba M, Wahyudi ST, Brunt DJ, Paradis N, Wu C (2021) Comput Biol Med 129:104156
- Daina A, Michielin O, Zoete V (2017) Sci Rep 7:42717
No competing interests reported.