Remdesivir analog as SARS-CoV-2 polymerase inhibitor: virtual screening of a database generated by scaffold replacement

By the end of 2019, a novel strain of the corona viral family named SARS-CoV-2 emerged in Wuhan, China and started to spread worldwide causing one of the most dangerous lethal pandemics. Researchers utilized various reported inhibitors and drug databases for virtual screening analysis against this novel strain. Later on, they succeeded to fish and repurpose remdesivir, an antiviral nucleotide analogue that inhibits RNA polymerase of the Ebola virus, as a promising candidate against SARS-CoV-2. In this study, we used the interactions of the co-crystallized metabolite of remdesivir with SARS-CoV-2 RdRp isozyme (PDB 7BV2) to design an analog with potential extra activity. This design was based on a scaffold replacement of a pyrrolotriazine moiety. This design was guided by a generated structure-based pharmacophore. The database generated from scaffold replacement was subjected to molecular docking and molecular dynamics simulations within the active site of SARS-CoV-2 RdRp (PDB 7BV2) to suggest HA-130383 and HA-130384 as potential lead compounds.


Introduction
Coronaviruses (CoV) represent a diverse family of enveloped, positive-sense, single stranded RNA viruses that are responsible for many respiratory and enteric diseases in humans. 1 The most notable highly pathogenic human corona viruses (CoVs) were Severe Acute Respiratory Syndrome corona virus (SARS-CoV) and Middle East Respiratory Syndrome corona virus (MERS-CoV). They were capable of causing a severe fatal respiratory disease. 1 Both SARS-CoV and MERS-CoV broke out and spread with mortality rates of 10% and 35%, respectively. 1 By the end of 2019, a novel strain of this family emerged in Wuhan, China and started to spread worldwide causing one of the most dangerous lethal pandemics in years. 2 Rapidly on the way to discovering a cure, trials were performed to determine one of the viral components that could be targeted and effectively inhibited. [2][3][4][5][6] Many research groups focused on targeting a virus-specic protein known as RNA-dependent RNA polymerase (RdRp) which is considered as one of 16 non-structural proteins (nsp) that were generated inside SARS-CoV-2. 1,7 This polymerase (nsp12) is responsible for viral replication. 8 These research groups utilized various reported inhibitors for virtual screening analyses against this viral target. [2][3][4][5]9 On June 2020 Yin, et al. reported and published the crystal structure of the nsp12-nsp7-nsp8 complex bound to the template-primer RNA and triphosphate form of remdesivir active metabolite (RTP). 10 They noticed that the nsp12, RdRp, possesses minimal activity on its own, and it required the addition of co-factors like nsp7 and nsp8 for full polymerase activity. 10 Later on, remdesivir, which was an antiviral nucleotide analogue that inhibit the RNA polymerase of the Ebola virus, was investigated as a promising candidate against SARS-CoV-2. 5, 11 Lo, et al. suggested the potential mechanism of inhibiting the SARS-CoV-2 RdRp for remdesivir. 12 According to that, remdesivir acts as a prodrug that is metabolized in cells to yield an active remdesivir triphosphate that could be used by the RdRp as a substrate, leading to the incorporation of remdesivir monophosphate (RMP) into the growing RNA product. 13 Finally, on October 2020, remdesivir showed promising results against SARS-CoV-2 in cell culture, animals and clinical trials was approved by FDA for the treatment of adult and pediatric patients over 12 years old suffering COVID-19 requiring hospitalization. 14 In spite of the aforementioned data, the use of remdesivir will remain limited by many factors, including the high price of dose regimen related to the low percent of decreasing the mortality rate. 15 In addition, its I.V. dosage form and its contraindication in both, patients suffering elevation of liver enzymes and those who suffer from renal impairment are important factors that limit its use. 16 Recently on 19 December 2020, WHO withdrew the use of remdesivir for the approved purpose. 15 Based on the previously mentioned drawbacks, remdesivir is considered as a proper lead compound and the improvement of remdesivir was an essential need to obtain new analogs with proper activity and minimum drawbacks. This target will be the research interest in this research work following steps demonstrated Fig. 1 as will be discussed later.

Protein preparation
SARS-CoV-2 polymerase (PDB 7BV2) was retrieved from protein databank and pyrophosphate, RNA, and water molecules were removed. Aer removal, the remaining residues were subjected to preparation by quick prep protocol in Molecular operating environment (MOE, 2019.0102) with the use of AMBER force-eld (Amber10: EHT) and R-Field solvation parameter 17 with standard settings. These settings include protonation of protein and ipping of ASN/GLN/HIS in Protonate3D. Also, renement for the protein was followed with RMS gradient of 0.1 kcal mol À1 A À1 .

Structure based pharmacophore
The interaction of the co-crystallized metabolite of remdesivir was the guide to select ve different features which were 1 Don (radius 1.5), 2 Aro (radius 1) and 2 Acc/Don (radius 1.3). This was achieved from compute/pharmacophore, query editor sequence from the main menu of MOE 2019.0102. The protocol was run according to EHT scheme (extended Huckle theory). 17 Application of Hückel theory to pharmacophore discovery. CICSJ Bulletin, 33(2), 33 and the strength of the different features was determined as Don > 0.5 and Acc > 0.8. Each feature was assigned from the interactions of cocrystallized ligand with the active site. The radius was adjusted to give a selective pharmacophore through the validation process. The external validation and selectivity were checked by retting of FDA-US approved database (Antiviral Library from Selleck, Catalog No. L7000) consisting of 302 drug molecules to give only 7 compounds including remdesivir.

Scaffold replacement
Pyrrolotriazine moiety was chosen as the main scaffold to be replaced. This could be done using scaffold replacement panel from compute/fragment/scaffold replacement. Selection of pyrrolotriazine, leads to 9 leaving atoms and 2 exit vectors. The linker database for this process consists of 40 626 fragments varied between aliphatic, aromatic and fused aromatic fragments. The process was constrained by two methods, lead-like descriptor lter and the structure-based pharmacophore. GBVI/WSA dG was the scoring function used in the protocol. Finally, generated database was minimized and clustered according to the score of the scoring function.

Molecular docking
The generated database was prepared to generate bioactive conformers according to the given protocols in MOE 2019.0102. The docking protocol for this process was restricted using the structure-based pharmacophore and renement of the docking process was performed using rigid receptor. The scoring function of the placement was London dG and for the renement was GBVI/WSA dG. The whole protocol was validated by re-docking of the co-crystalized remdesivir fragment to give RMSD value 0.2178Å. The given compounds were clustered by docking score and the best compound was subjected to another docking process with induced t receptor renement.  parameterized using SwissParam 19 while protein was prepared using all-atom CHARMM36 force eld. 20 Simulations were performed in TIP3P explicit water 21 and neutralized by Na + ions. Minimization was done using steepest descent algorithm and maximum force was set to less than 1000 kJ mol À1 nm À1 . Systems were initially equilibrated using rst constant volume and temperature (NVT ensemble) and then constant pressure and temperature (NPT ensemble) for 1 ns each. Parrinello-Rahman barostat 21 and V-rescale algorithm 22 were used to control pressure at 1 atmosphere and temperature at 300 K respectively. The LINear Constraint Solver (LINCS) algorithm 23 was used for bond's length constraints and particle mesh Ewald (PME) method 24 was used for long-range electrostatics calculations. For all simulations 2 fs timestep was used. van der Waals cut-off distance (rvdw) was set to 1.2 nm. Initial coordinates were taken from docking pose.

Rational of the design
Remdesivir is a prodrug which is activated intracellularly by breakage of side chains attached to phosphate group as shown in Fig. 2. It follows the approach to design some antiviral drugs as protides (prodrugs of nucleotides). 25,26 This pathway of design was based on the use of nucleoside monophosphate skeleton capped with an aryl group and an amino acid ester (a phosphoramidate). 27 This process was clear as the metabolic process of remdesivir were extensively studied as reported. 28 In addition, the metabolite appeared co-crystallized with RdRp of SARS-CoV-2 as shown in the PDB deposit 7BV2. 29 Due to this outcome, our work will focus on nding remdesivir analogs from a replacement of a scaffold present in the metabolite. Therefore, we designed a workow for this study Fig. 2 Metabolite of remdesivir and the scaffold of interest that will be replaced to get more active remdesivir analog. Fig. 3 The generated pharmacophore based on the interaction of remdesivir metabolite within the active site of RdRp SARS-CoV-2 (PDB 7BV2). Pharmacophore query consists of 1 Don, 2 Aros and 2 Accs %26 Dons. depending on scaffold replacement of pyrrolotriazine ring of remdesivir metabolite with respect to retaining the ve membered sugar to keep the nucleoside skeleton to follow the protide approach. The side chain of remdesivir is used to ensure the hydrophobicity of the drug to penetrate the cell then it will be removed by metabolic enzymes. Our protocol will proceed through the following steps; we rst create a structure-based pharmacophore which will used to retain the important contacts of remdesivir fragment inside the active site of RdRp. This process will be followed by molecular docking of the generated database form the scaffold replacement process.

Generation of structure-based pharmacophore
The interaction of co-crystallized remdesivir fragment was the guide to choose the annotation points for the pharmacophore query. Remdesivir fragment displayed a hydrogen bond between 5-CN group of the tetrahydrofuran ring with Asn691 and 4-OH from the same ring with Arg555. Different features were selected to retain the interaction (Fig. 3); two hydrogen bond acceptor and donor (Accs & Dons) features to retain the interaction with Arg555. In addition, two aromatic (Aro) features were added to represent the presence of fused system to retain the spatial arrangement of the contacting groups. Finally, a donor (Don) feature was selected due to the 4-amino group on pyrrolotriazine ring as it displayed a hydrogen bond with the dissecting RNA in the 3D structure with U20 as reported. 29 The generated pharmacophore was validated by pharmacophore search for a decoy test set consists of 302 US-FDA approved drugs (Antiviral Library from Selleck, Catalog No. L7000) for the treatment of SARS-CoV2. Fitting of these compounds to the generated pharmacophore would give 7 compounds including remdesivir as shown in Table 1.

Scaffold replacement
Pyrrolotriazine fused system will be the target of the interest in this process. This process will be constrained by the generated structure-based pharmacophore. In this process we will focus on occupying the side area of this fused system as shown in Fig. 4.
The linker database for this step is available in MOE 2019.0102 and consists of 40 626 linkers. This process resulted in a database consisting of 1231 remdesivir metabolite analogs. In the given database, binding free energy was calculated for any given structure according to the scoring function GBVI/WSA dG. 30 The given metabolites were clustered and arranged according to their scores.

Molecular docking
The given database from scaffold replacement was subjected to molecular docking simulations within the active site of RdRp SARS-CoV-2 (PDB 7BV2). The docking process was restricted by the structure-based pharmacophore to ensure the proper orientation with respect to the active remdesivir fragment. The compounds were clustered according to score and it was observed that all of them with an extension to the replaced fused system and the top ve ranked compounds are shown in Table 2 Table 2. Fig. 4 The targeted side area of the pyrrolotriazine fused system of remdesivir within the active site of RdRp SARS-CoV-2 (PDB 7BV2).
Metabolite HA1 was considered as a good candidate as it was characterized by extra interaction due to the side chain of the fused ring system. It was found that the benzyl ring occupies the vacant pocket to form a p interaction with Arg555 with retaining all other interactions of the co-crystallized metabolite of remdesivir as shown in Fig. 5A.
While metabolite HA2 that was ranked as the second good candidate showed a successful occupying of the vacant pocket nearby Arg555 through its aliphatic side chain in addition to retaining all other interactions of the co-crystallized metabolite of remdesivir as shown in Fig. 5B. All poses of other metabolites were illustrated in Fig. 5.
For more investigation about the given poses, molecular dynamic simulations for the ve given docking poses would be taken into considerations in order to assess the stability of these docking poses in the active site.

Molecular dynamic simulations
Molecular dynamics simulations were done to study the stability of the top ve compounds in docking study (Table 2) and compare their stability to the co-crystalized ligand (control) as well as the apoprotein. The complexes of these compounds in the active site of SARS-CoV-2 RdRp (PDB 7BV2) were used and their coordinates were obtained from docking poses. The protein used in the molecular dynamics included all the 3 protein chains (nsp12, nsp7 and nsp8). The full nsp12-nsp7-nsp8 complex is known to have the RNA polymerase activity 15 and hence was included during the MD simulation. Aer initial minimization and equilibration, the complexes were run for 100 ns production run. Data was collected every 5000 steps (10 ps) of the production run giving 10 000 points to analyze.
One of the most important parameters that could be extracted from the production run trajectories is the root mean square deviation (RMSD) of the ligand heavy atoms. A high uctuation in the RMSD of ligand heavy atoms indicate that the ligand is not stable in the active site. Fig. 6 shows a plot of ligands heavy atoms RMSD. For examples, the position of the  ligand heavy atoms of HA1 was found to be highly uctuating during the production run during most of the production run especially aer 80 ns. An analysis of the studied ligands heavy atoms RMSD is shown in Table 3. As seen in the table, control has a uctuation of 10.8% of the RMSD average. Among tested compounds (HA1-HA5), HA2 and HA3 showed least uctuation and hence best stability. Another important factor that contributes to the stability of ligands in the active site is the number of hydrogen bonds between ligands and proteins. The number of hydrogen bonds between the six studied ligands and target proteins are plotted in Fig. 7. Analysis of these results are shown in Table 4. Compound HA3 has shown the maximum number of hydrogen bonds (4.24 hydrogen bonds) followed by HA4 and HA2 which has shown a number of hydrogen bonds comparable to control (2.07 hydrogen bonds).
In addition, root mean square uctuation (RMSF) of SARS-CoV-2 RdRp protein residues showed similar results for apo- protein, co-crystalized ligand as well as tested ligands (Fig. 8). It is important to notice that HA4 has shown slightly higher uctuation in loops around residues 40 and 232 which might indicate movements of these loops to accommodate HA4.
These results suggest the higher stability of HA2 and HA3 in the active site of SARS-CoV-2 RdRp during production run. In addition, total interaction energy between ligand and protein was calculated as the total of short-range Lennard-Jones (LJ) and Coulomb interaction. The total interaction energy of HA2 (À278.00 AE 34.43 kJ mol À1 ) and HA3 (À393.13 AE 20.03 kJ mol À1 ) were found signicantly better compared to other ligands such as HA1 (À46.14 AE 12.35 kJ mol À1 ) and HA5 (À171.53 AE 22.38 kJ mol À1 ). These results further support the superiority of HA2 and HA3. As a result, the suggested prodrugs in Fig. 9 are potential candidates for further investigations against RdRp of SARS-CoV-2.

Conclusion
In this study, we tried to obtain a remdesivir analog through studying of the interaction of co-crystallized remdesivir metabolite inside active site of RdRp. The obtained results were the guide to choose the pyrrolotriazine as a scaffold of interest and then annotation points for the pharmacophore query. Creation of a structure-based pharmacophore was performed, and different features were selected to retain the interaction. The generated pharmacophore was validated and a test dataset was tted to the generated pharmacophore that gave 7 compounds including remdesivir. During this process the linker database supported by MOE 2019.0102 was used and focused on occupying the side area of pyrrolotriazine fused system. The given metabolites were clustered and it was observed that most of compounds with extension of the replaced fused scaffold were having the best score. The given database from scaffold replacement was subjected to molecular docking simulations within the active site of RdRp SARS-CoV-2 (PDB 7BV2). The results showed that metabolite HA2 and HA3 and their corresponding prodrugs HA-130383 and HA-130384 are considered as a good candidate as it has higher stability and lower energy compared to the other studied ligands.