Discovery of inhibitors of the SARS-CoV-2 S-protein through drug computational drug repurposing
We present here a computational multistep procedure that tries to find compounds from the Selleck FDA approved drugs and the Selleck database of Natural Products able to interact with the S-protein of SARS-CoV-2 and inhibit its interactions with human ACE2. The multistep procedure uses a sieving approach, starting from low cost computational methods that filter out compounds to pass to next, more computationally expense, stages. A standard docking is applied together with a pharmacophore biased one, and protein flexibility is taken into account from a clustering of four 3 𝜇s trajectories. After visual inspection of the final candidates we propose a list of 10 compounds to be essayed experimentally.
Introduction<!>METHODS 2.1. System preparation<!>Equilibration<!>Molecular dynamics simulations<!>Clustering<!>Pharmacophore determination<!>Virtual screening<!>Pharmacophores<!>Virtual screening
<p>Coronaviruses are positive single-stranded RNA viruses, and are classified into four genera: 𝛼, 𝛽, 𝛾 and 𝛿. 1 They infect a wide range of hosts, including birds, pigs, and humans. 2 Although most family members produce mild respiratory effects in humans, new members have recently appeared that produce serious illness in affected individuals. We can mention as previous variants to the one that has been studied here, 229E, NL63, OC43 and HKU1, which cause symptoms of the common cold, and two more problematic, such as SARS-CoV-1 and MERS-CoV. SARS-CoV-1 was identified as the pathogen responsible for an outbreak of severe acute respiratory syndrome (SARS) in Guangdong Province, China, in 2002, and ten years later MERS-CoV was identified in the sputum of a patient who was retrospectively diagnosed with Middle East Respiratory Syndrome (MERS) in Jordan. 3,4 Both pathogens produced an epidemic that spread to several countries due to international travel by infected people, epidemics that ended about a year after each outbreak after taking strict infection control measures. 5 SARS-CoV-2 was first reported in Wuhan, Hubei Province, China, in December 2019. Severe acute respiratory syndrome caused by this virus, whose symptoms include respiratory illness and pneumonia 6 was named COVID-19, and the rapid spread of SARS-CoV-2 around the world led the World Health Organization to declare the pandemic on March 11th, 2020. As of June the 21th, 2022, there have been 536590224 confirmed cases and 6316655 deaths from COVID-19. 7 Humans are infected commonly and throughout the globe with seasonal coronaviruses, which usually cause a respiratory disease with mild symptoms. Hence, they have been never treated as an important public health threat, and development of a specific anti-viral treatment has not been considered a priority. Consequently, when SARS-CoV-2 emerged, there were no specific antiviral treatments for COVID-19.</p><p>The development of a drug requires on average more than two billion dollars and can take up to 10-15 years. 8 This is why drug repurposing was considered a valuable strategy to deal with the pandemic caused by SARS-CoV-2 in a reasonable time. 9 In addition, computational discovery of repurposed drugs helps avoid costly trial-and-error experiments, which require cell cultures, biochemical screening, and living systems. 10 Thus, baricitinib (rheumatoid arthritis drug) was proposed, using artificial intelligence techniques, as a repurposed drug 11 and was subsequently approved as an emergency treatment by the FDA (US Food and Drug Administration) for the treatment of COVID-19 in combination with redemsivir. 12 Other antivirals already marketed such as favipiravir 13 and EIDD2801 14 show mixed evidence, while drugs such as lopinavir, ritonavir, hydroxychloroquine and chloroquine eventually proved ineffective. 15 Currently, the clinical treatment of COVID-19 is mainly symptomatic, using anti-inflammatories such as dexamethasone 16 or cytosine inhibitors, combined with antibiotics to treat secondary infections. There is therefore an urgent need to develop specific antiviral therapies against SARS-CoV-2.</p><p>SARS-CoV-2 has a sequence identity of 79.6% with SARS-CoV-1 and 96.2% with bat coronavirus. 3 From this high similarity in nucleotide sequences it was inferred that, like SARS-CoV-1, SARS-CoV-2 uses a highly glycosylated protein found in its envelope, called S-protein, to bind to the cell through interaction with the angiotensin-converting enzyme, ACE2 3,17 . The S-protein, made up of the S1 and S2 subunits, plays a crucial role during the first contact between the virus and the ACE2 receptor. The S1 subunit binds to ACE2 and the S2 subunit is activated by the protein serine protease TMPRSS2, which promotes membrane fusion and cell infection. 18 Both SARS-CoV-1 and SARS-CoV-2 use ACE2 to enter host cells, so blocking this critical interaction would prevent the entry of the virus. 19 The S1 subunit contains the so-called receptor binding domain, RBD, a fragment of it that is sufficient to bind to the peptidase domain of ACE2. 17,20 Most studies have focused on this domain, because it binds to cell receptors and, in addition, is the main target of neutralizing antibodies during infection. [17][18][19][20]21 In the present study we have carried out a virtual screening of S-protein inhibitors, focusing on the area of interaction between the RBD and the ACE2 protein. The screening has been done on two databases of molecular structures. On the one hand, Selleck FDA approved drugs 22 (FDA from now on), a database of FDA-approved drugs, and on the other Selleck database of Natural Products, (NAT from now on) database of natural compounds. 23 We have included the flexibility of the protein through the study and subsequent clustering of classical molecular dynamics trajectories, and in addition, different docking methodologies have been used, such as those included in the Qvina 24 and MOE programs. 25 The use of the MOE program allows the use of a pharmacophore guided docking, for which two possible pharmacophores were defined: one focused on imitating the best interactions between the RBD of the SARS-CoV-2 protein-S and ACE2 (COV2-ACE2), and another focused on imitating the interactions that improve the most when moving from the RBD system of the SARS-CoV-1 and ACE2 protein-S to the SARS-CoV-2 plus ACE2 system. The second pharmacophore focuses, therefore, on trying to describe what changes occur in the RBD of SARS-CoV-2 versus SARS-CoV-1 so that the first binds with much greater affinity to ACE2. The results of these different dockings were be subjected to successively more computationally complex methodologies, so that a virtual screening process of a total set of 2791 compounds ended in a set 48, of which 10 were subjected to experimental tests.</p><!><p>The crystallographic structure of the RBD SARS-CoV-2 plus ACE2 system, with access code 6M0J, 26 was the starting structure for this part of the study. As the S-protein is glycosylated with N-acetylglucosamine (NAG), and due to the high computational cost of mixing the parameters of the force fields of proteins and saccharides, it was decided to eliminate these NAG residues, since they are far from the binding zone with ACE2. Besides, hydrogens were added to the residues of both proteins according to the physiological pH, and the orientation of the side chains was adjusted using the Protein Preparation Wizard workflow included in the Maestro software. 27 For the treatment of the zinc atom present in the ACE2 protein, the dummy cationic approach was used. 28 Using the LeaP module of AMBER20 29 the system was placed in a cubic box filled with TIP3P water molecules, 30 establishing a minimum distance of 15 Å between the solute and the walls of the box, and those water molecules closer to 2 Å were removed from any atom in the complex. LeaP was also used to add the necessary counterions, in the positions of lower electrostatic potential, to neutralize the system. All calculations were performed using the ff14sb 31 force field with a cutoff of 10 Å for non-covalent interactions. Electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method. 32</p><!><p>Before starting the molecular dynamics (MD) calculations, the structure was relaxed to eliminate possible steric clashes in a multistep minimization procedure already optimized by our group, 33 using the steepest descent method.</p><!><p>After the minimization process, a heating and equilibration protocol was followed that includes first a heating up to 300 K at a rate of 1.5 K ps -1 under the NVT ensemble, followed by a simulation in the NPT ensemble in order to equilibrate density, both steps using the parameters already optimized by our group. 33 After the equilibration protocol we started the so-called production stage, where in addition to the SHAKE 34 algorithm used in the equilibration stage, which allows the integration time to be 2 fs, the hydrogen mass partitioning 35 scheme was applied to use an integration time of 4 fs. In order to improve conformational sampling by using several molecular dynamics instead of a longer 36 one, we opted to obtain and analyze 4 MDs, each of 3 𝜇s.</p><!><p>The four trajectories were joined into one, the water molecules and the ACE2 protein were eliminated and a multi-stage process was followed in order to have a trajectory ready to analyze and cluster the major conformational changes for the zone of interaction of the SARS-CoV-2 RBD with ACE2, using the AMBER20 cpptraj module. It was sought to superimpose all the snapshots of the global trajectory to the first one, using as reference atoms for this superposition those C𝛼 of the RBD of the zone of interaction with ACE2 that had the least fluctuations, thus focusing the clustering precisely in those with the largest fluctuations. In a first stage, the trajectory is superimposed to the first photo and the root mean square factor (RMSF) is calculated, then we select for the next superposition the C𝛼 of residues with RMSF < 2 Å in the interaction zone. We then lowered the cutoff to RMSF < 1 Å, and finally the minimum cutoff for RMSF that allowed a reasonable number of residues obeying the criterium was 0.6 Å.</p><p>Clustering was carried out using the average linkage 37 , obtaining two populated clusters, one with 461732 members (from now on cluster 1), and another with 101971 members (cluster 2), whose centroids will be the two structures of the SARS-CoV-2 RBD on which the docking with Qvina will be carried out. Thus, the workflow for virtual screening in several stages from docking with QVina will be carried out 4 times: FDA database on cluster 1, FDA database on cluster 2, NAT database on cluster 1 and NAT database on cluster 2.</p><!><p>Following the same protocol of preparation, minimization and production already described, 4 MDs of 3 𝜇s were obtained for the system composed of RBD of SARS-CoV-1 and ACE2, starting from the experimental structure with PDB code 2AJF. 17 For the analysis of interactions that lead to the obtainment of a pharmacophore, it would be convenient to use only those sections of the MD run that are considered energetically converged. Therefore, using the MMPSA.py 38 script of AMBER20, a calculation of the evolution of the Gibbs free binding energy between each RBD (i.e. SARS-CoV-1 and SARS-CoV-2) with ACE2 against time was carried out. Specifically, the MMGBSA 39 approach was used using for the generalized Born (GB) method the parameters optimized by Onufriev et al. 40 Both in the case of SARS-CoV-1 and SARS-CoV-2 it was concluded that the interaction analysis should be done on the last 200 ns of each of the 4 MDs obtained for each system.</p><p>Thus, an analysis of the established hydrogen bonds together with their strength was carried out (the strength of each hydrogen bond can be correlated with the percentage of snapshots of the trajectory in which that bond is being established), as well as an MMGBSA energy decomposition calculation, so that we can obtain a list of those ACE2 residues that establish better interactions with SARS-CoV-2 from which we can propose pharmacophore 1, and a list with those residues showing the biggest differences when moving from SARS-CoV-1 to SARS-CoV-2, from which we can propose pharmacophore 2. Therefore, the workflow for the multi-stage virtual screening from pharmacophore biased docking with MOE will be performed 4 times: FDA database with pharmacophore 1, FDA database with pharmacophore 2, NAT database with pharmacophore 1 and NAT database with pharmacophore 2. This sums a total of 8 virtual screening workflows described below.</p><!><p>We propose a sieving methodology which starts from a docking performed with QVina or MOE on both full databases, from which in different stages of increasing computational requirements a decreasing number of compounds is being analyzed. All the scripts used within this methodology are available in TO BE DECIDED. Thus the search box for the docking with QVina was defined manually with the help of the UCSF Chimera program. 41 After the execution of the four molecular docking procedures, the structures found were ordered in descending order of interaction, selecting for the next step those poses with an interaction energy better than -6,8 kcal mol -1 . Docking with MOE was performed in two steps. In the first one, after defining each pharmacophore, a docking directed to the pharmacophore of each database with no subsequent refinement (that is, each docked structure is not minimized once it has been placed in the corresponding binding site) is performed. The structures selected from the docking are subjected to a second step of docking, now asking for minimization and obtaining list of compounds sorted by interaction energy. The structures obtained in the previous step were solvated with TIP3P waters after having obtained the GAFF 42 parameters for each ligand, and underwent a minimization process of 5000 steepest descent steps. For each minimized structure, the free binding energy was calculated with both MMGBSA and MMPBSA and a consensus criterion is used to choose the 100 best structures that are among the 100 best in both MMGBSA and MMPBSA calculations.</p><p>The structures selected in the previous stage followed the heating and equilibration protocol already described to finally obtain for each one a trajectory of 100 ns, which was used to perform an MMGBSA calculation. Graphs of Gibbs free energy (∆G) binding versus time were used in a program created in our group to filter the best structures to move on to the next stage of virtual screening. This program calculates the average ∆G every 10 ns, and compares it with the previous averages, determining if this latter average deviates by more than twice the standard deviation calculated until that moment. If the given trajectory has a ∆G versus time that is stable at least for the final 25% of the trajectory we consider that the evolution of ∆G versus time is good and that structure proceeds to the next stage. The next step consists of an increase of the MD run with additional 100 ns and a re-analysis of the MMGBSA results against time to further filter the structures that will pass to the final stage. Finally, after obtaining a total of 500 ns of MD, an automatic selection of compounds whose profile of ∆G versus time is stable is obtained in the same way, although this time the criterion is hardened to require that it be so at least during the final 50% of the trajectory. The final analysis of the structures was also carried out manually, and the final selection of compounds proposed for experimental evaluation was made by visual inspection and manual comparison of the energy profiles against time, final values of ∆G, LE (ligand efficiency) and strength of established hydrogen bonds as determined with cpptraj.</p><!><p>Figure S1 shows the MMGBSA calculated binding energy versus time for each of the 4 MD runs (3 𝜇s) obtained for RBD of SARS-CoV-2 + ACE2 and RBD of SARS-CoV-1 + ACE2, from which it can be inferred that the analysis of intermolecular interactions between the RBD and ACE2 can be performed for the last 200 ns of simulation. This analysis can be summarized, in the case of the RBD of SARS-CoV-2 + ACE2 system, in Tables 1 and 2, where the best residues of ACE2 as concerns binding energies and hydrogen bonds, respectively, are shown. As the pharmacophore has to include points which are relatively close to each other, we chose residues ASP338, LYS336 and TYR24 of ACE2 to define the pharmacophore. Due to their good van der Waals interactions, the center of mass of each of these three residues is defined as a pharmacophoric point, while the center of mass between OD1 and OD2 of ASP338 defines a hydrogen acceptor point, and O of LYS336 another hydrogen acceptor point. GLU_18@OE2,OE1 GLN_762@HE22 GLN_762@NE2 72.8 2.831 163.625 ASP_21@OD1 GLN_767@HE22 GLN_767@NE2 36.1 2.831 164.065</p><p>On the other hand, Table 3 shows those residues of ACE2 with better than -1 kcal mol -1 MMGBSA binding energy with RBD SARS-CoV-2 as compared with RBD SARS-CoV-1. It is clear from this table that the better interactions of these residues come from the electrostatic component of the Gibbs free energy, which can be attributed to higher number of hydrogen bonds established with RBD SARS-CoV-2 as compared with RBD SARS-CoV-1. Indeed, although the hydrogen bonds established with SARS-CoV-2 are not always established with the same residue, the percentage of the simulation time that a hydrogen bond is formed between ACE2 and RBD SARS-CoV-2 is 89.2%, 72.8% and 67.5% for ASP21, GLU18 and GLU20, respectively. All these three residues act as hydrogen acceptors, so that we add 3 hydrogen acceptor pharmacophoric points to our pharmacophore, adding an option in MOE that it is not needed that all these eight pharmacophoric points have to be fulfilled at the same time, but only 4. Table 3. Residues of ACE2 with differences in MMGBSA binding energies higher than -1.0 kcal mol -1 , with RBD SARS-CoV-2 and RBD SARS-CoV-1, showing also their van der Waals (VDW) and electrostatic components. All units in kcal mol -1 .</p><!><p>As explained in the Methods section, a multistage virtual screening protocol was performed in order to select a relatively small number of poses to be manually analyzed and to obtain a list of compounds to be experimentally tested. Thus, we used QVina to dock the FDA and NAT databases on the 2 centroids of the clusters found from the MD trajectory, and we used MOE to run a pharmacophore guided docking using 2 different pharmacophores and the centroid of cluster 1 as reference structure. Figure 1 allows to visualize how the number of poses diminishes as the computational requirements of each stage increase. One has to take into account that the number of poses does not necessarily equal the number of compounds, as in some cases different poses of the same ligand are considered. As an example, the docking with QVina of the FDA database (919 poses) allowed to select 553 to proceed to the MMGBSA and MMPBSA one-snapshot calculation, from which 133 poses were selected to pass to the next stage. Then, an MMGBSA Gibbs free energy versus time analysis of a 100 ns trajectory of the 133 poses allow to select 83 poses which will be simulated 100 ns more. These 200 ns filtered from 83 to 30 the number of poses subjected to a total of 500 ns of MD simulation, from our script automatically selected 11 poses. We also individually visualized each of the 30 500 ns MD runs to finally select 6 poses. This procedure has allowed to select 48 poses out of 2791 to be further analyzed.  Thus, Tables 4-6 allowed us to obtain a list of ten compounds to be experimentally essayed against SARS-CoV-2. We have manually selected these compounds taking together into account the values of LE, the establishment of strong hydrogen bonds, and also trying to propose candidates from both databases and both docking procedures. Taking all these into account the compounds essayed were glyburide, clofoctol, garenoxacin, telmisartan and pazopanib from the FDA database, and fenretidine, fisetin, glycochenodeoxycholic acid, engeletin and estradiol cypionate from the NAT database.</p>