Targeting multiple conformations of SARS-CoV2 Papain-Like Protease for drug repositioning: An in-silico study.

Papain-Like Protease (PLpro) is a key protein for SARS-CoV-2 viral replication which is the cause of the emerging COVID-19 pandemic. Targeting PLpro can suppress viral replication and provide treatment options for COVID-19. Due to the dynamic nature of its binding site loop, PLpro multiple conformations were generated through a long-range 1 micro-second molecular dynamics (MD) simulation. Clustering the MD trajectory enabled us to extract representative structures for the conformational space generated. Adding to the MD representative structures, X-ray structures were involved in an ensemble docking approach to screen the FDA approved drugs for a drug repositioning endeavor. Guided by our recent benchmarking study of SARS-CoV-2 PLpro, FRED docking software was selected for such a virtual screening task. The results highlighted potential consensus binders to many of the MD clusters as well as the newly introduced X-ray structure of PLpro complexed with a small molecule. For instance, three drugs Benserazide, Dobutamine and Masoprocol showed a superior consensus enrichment against the PLpro conformations. Further MD simulations for these drugs complexed with PLpro suggested the superior stability and binding of dobutamine and masoprocol inside the binding site compared to Benserazide. Generally, this approach can facilitate identifying drugs for repositioning via targeting multiple conformations of a crucial target for the rapidly emerging COVID-19 pandemic.

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the relevant virus for the Coronavirus disease (COVID-19) and was originally discovered in Wuhan, China [[1], [2], [3], [4]]. Based on the situation report of the World Health Organization (WHO) - November 2020, COVID-19 is extremely spreading worldwide over 230 countries and accountable so far for > 111 million cases and >2 million fatalities (https://covid19.who.int/, accessed on Feb 25, 2021). The coronavirus is part of a large family of enveloped single stranded RNA genome (ssRNA) that belong to the Coronaviridae family. They are categorized into four genera: alpha, beta, gamma, and delta coronaviruses [5]. Some of coronaviruses prompted several respiratory diseases, such as SARS-CoV [6], middle east respiratory syndrome coronavirus (MERS-CoV) [7] and the pandemic COVID-19 [3]. SARS-CoV-2 are beta coronaviruses [3,8] with symptoms commonly like other respiratory viruses infection, such as influenza and rhinovirus [9].

Papain-like protease (PLpro) and 3C-like/main protease (3CLpro/Mpro) [10,11] are critical for the release of 16 non-structural proteins (nsps1-16) after processing the two large polyproteins, pp1a and pp1ab [12]. These polyproteins are created following the virion entry to the host cell where their production is initiated via translation of 5′-terminal open reading frames (ORF1a and ORF1ab) [13]. It is known that the establishment of the replicase complex necessary for viral genome replication is conditional on nsps [14]. Nonetheless, PLpro performs a key role for the release of nsp 1-3 from the viral polyprotein which is critical for viral replication. Furthermore, PLpro has been reported to negatively regulate the host innate immune response towards the viral infection by its deubiquinating and deISGylating effect [15,16]. Therefore, PLpro has been identified as an imperative target for viral replication suppression attempts in SARS-CoV and SARS-CoV-2 [15,17].

Structure-based virtual screening (SBVS) remains an essential technique in modern drug discovery [[18], [19], [20], [21]]. Molecular docking is widely employed in SBVS campaigns and computational drug repositioning for COVID-19 [[22], [23], [24], [25], [26], [27], [28], [29], [30]], utilizing the structural information of the protein targets to gauge molecular databases and predict the favored binding of compounds. However, the docking tool should be evaluated utilizing benchmarking molecular sets [31,32]. Also, ensemble docking has emerged as a popular approach which incorporates the protein flexibility and tackles the problem of false positive results of rigid docking [33]. This can be achieved via Molecular dynamics (MD) simulations that has shown to be a powerful approach for sampling different protein conformations. Usually, screening this ensemble of conformations gives better results than single crystallographic structure [34].

Accordingly, it is the target of this study to investigate samples of the conformational space of SARS-CoV2 PLpro. One reason for this is the limited number of X-ray structures deposited in the Protein Data Bank (PDB) for PLpro (around 24 structures) compared to the main protease (Mpro, around 200 structures). Additionally, compared to Mpro [[35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45]], fewer reports about targeting PLpro is presented in literature [[46], [47], [48], [49], [50], [51], [52]]. In the context of drug repurposing to treat COVID-19, Balasubramaniam et al. [46] reported potential targeting of three SARS-CoV-2 proteins, RNA-dependent RNA polymerase, papain-like proteinase and helicase, by the antiviral drug elbasvir through virtual screening of 54 FDA-approved antivirals and 3300 investigational drugs. Kandeel et al. proposed some FDA-approved drugs and supplements to target PLpro. Their methodology was based on molecular dynamics (MD) simulations followed by molecular mechanics/generalized Born surface area (MM/GBSA) binding energy calculations [51]. Also, Cavasotto and Di Filippo presented a docking‐based screening using a quantum mechanical scoring of a library built from approved drugs and compounds undergoing clinical trials, against three SARS‐CoV‐2 target proteins, including PLpro [52]. Kouznetsova et al. [48] repurposed FDA-approved drugs along with inhibitors of Hepatitis C Virus (HCV) and Cytomegalovirus (CMV). De Vita et al. [49] adopted an FDA-approved drugs screening approach combining MM-GBSA calculations for ranking drugs along three different SARS-CoV-2 targets, including PLpro, and MD simulations for the best ligand in each target member. Delre et al. [50] suggested known drugs as both covalent and non-covalent inhibitors for PLpro depending on covalent docking, non-covalent docking and MM-GBSA calculations.

The mentioned approaches lacked to objectively consider different conformations of the highly flexible loop (T265-H272). This encouraged us to explore the conformational space of PLpro for screening purposes. For this we conducted a long-range molecular dynamics (MD) simulation of PLpro for 1 micro-second. Then the simulation trajectory was clustered to extract representative structures for ensemble docking purposes. Furthermore, guided by our recent benchmarking study of SARS-CoV2 PLpro [13], we employed FRED [53,54] for ensemble screening of DrugBank database against the cluster representatives. The results highlight potential consensus binders to many of the MD clusters as well as the newly introduced X-ray structure complexed with small molecule (PDB ID: 7JRN). Generally, this approach can help in identifying drugs for repositioning via targeting multiple conformations of a crucial target for the rapidly evolving COVID-19 pandemic.

Results and discussion

Structural aspects and molecular dynamics

The 24 structures of PLpro comprises the following PDB IDs: 6W9C, 6WRH, 6YVA, 6XA9, 6XAA, 6WZU, 6XG3, 7CJD, 7D6H and 7D47 for the apo structures, while 6WUU and 6WX4 for complexes with a peptide-like inhibitors. Also, PDB IDs of 7JRN, 7JN2, 7JIR, 7JIT, 7JIV, 7JIW, 7CJM, 7KOK, 7KOL, 7KOJ, 7D7T and 7CMD have been recently introduced for protein - small molecule complexes (Resolution and R-Value Free are shown in Table T1 in the Supplementary Material).

Based on previous reports and the recently deposited X-ray PDB structures of PLpro, Tyr268 and Gln269 were proposed to be the key residues for small molecule recognition [13,17]. They encompass the main part of a flexible loop that can adapt various backbone and side chain conformations, as seen in Fig. 1 . However, their unbound X-ray structures appeared to show close conformations (pairwise RMSD < 1A) and their side chains mostly appear to point outwards to the solvent exposed area [13]. Unlike the unbound conformation, the structures complexed with small molecules appeared to show more distinct and ordered conformations for the key Tyr268 and Gln269 where their side chains pointed inwards to offer a hydrophobic wall for interactions with the small molecules [13]. Therefore, based on the data available from the PDB, two main X-ray structures can be extracted representing distinct conformations of the binding site loop (see Fig. S1 in Supplementary Material), e.g., PDB ID: 6W9C and 7JRN for both the unbound and complexed PLpro structures, respectively.

Fig. 1

Superposition of some X-ray structures of SARS-CoV2 PLpro from the PDB (A) illustrating the flexibility and various conformations of the key Tyr268 and Gln269 residues of the binding site loop (B). The color code is assigned randomly.

Hence, we aim in this section to investigate samples of the conformational space of SARS-CoV2 PLpro with a focus on the conformational changes of the key binding site residues. Molecular dynamics simulations were carried out to explore the protein conformational space of the apo/native form (PDB ID: 6W9C) and give insights into the dynamic nature of the binding site. Long-range molecular dynamics simulations for 1 micro-second ware performed to give robust sampling of the protein conformations. The simulation system was subjected to a preliminary minimization step followed by an equilibration step of 6 ns. The system was then subjected to production run of 1 micro-second. Fig. 2 shows Radius of Gyration (Rg) and Root Mean Square Deviation (RMSD) for the protein backbone atoms through the 1 micro-second simulation, where radius of gyration is a measure of the change in the compactness of the protein during the simulation. Its low range of 1 Å (48-49 Å) indicates the stability of the protein during the simulation and the absence of major conformational change in protein structure. While RMSD is a measure of the deviation of the backbone atoms from the initial protein structure. It is observed that RMSD has a range of 3 Å with lowest and highest values of 1 Å and 4 Å, respectively, which gives an indication of the stability of the protein throughout the simulation. It is also observed that the trajectories are divided into three main ensembles. The first ensemble has RMSD of 1–3 Å (0–300 ns), the second and third parts are interchangeable from 300 to 1000 ns showing fluctuation in their RMSD from 2 to 4 Å. Clustering for the simulation trajectory was performed to extract representative structures of the protein different conformations. The clustering algorithm used calculates the average RMSD between trajectory frames and clusters close frames together. In a first run, we assigned the RMSD measurements to be calculated on the backbone atoms. Fig. 2 shows bar plot of timeline frame distribution of the three resulting clusters with a frame distribution between clusters similar to the distribution of backbone RMSD through the simulation run shown in Fig. 2, which supports the aforementioned observation.

Fig. 2

Top: Protein radius of gyration (R) through the 1 micro-second simulation time showing low fluctuation of 48–49 Å with low range of 1 Å indicating the protein compactness during the simulation. Middle: protein RMSD of backbone atoms during the simulation time showing three clusters for the protein with the first cluster from 0 to 300 ns, and the second and the third clusters from 300 ns to 1 micro-second. Bottom: Timeline bar plot of protein cluster distribution during the 1 micro-second simulation assigning the RMSD matrix to be measured on the backbone atoms. It shows three clusters corresponding to the three clusters concluded from protein RMSD measurement, with the first cluster in red and second and third clusters in blue and green, respectively.

Fig. 3 shows per residue Root Mean Square Fluctuation (RMSF) measured for each residue in the protein. RMSF appears to be relatively high in the regions of amino acids: 3–83, 182–200, 223–235 and 306–315, however, such regions are distant from the binding site. Also, RMSF analysis shows the highest fluctuation (>5 Å) in amino acids region T265-H272 comprising the binding site loop. This is not surprising due to the high flexibility nature of this loop, as inspected in Fig. 1 and reported earlier [13]. The rest of the binding site amino acids (L162-E167, P247–P248, Y264–Y273 and T301) showed lower RMSF values < 1 Å indicating their stability during the simulation. The rest of the protein amino acids displayed low RMSF values highlighting their stability as well (Fig. 3).

Fig. 3

Top: PLpro apo structure in khaki color (PDB ID: 6W9C) showing the binding site residues with the highly flexible loop (T265-H272) in red and the rest (L162-E167, P247–P248, Y264–Y273 and T301) in blue. The highly flexible amino acids (3–83, 182–200, 223–235 and 306–315) are shown in yellow. Bottom: per residue RMSF with color notation corresponding to the protein (top) color notation. It shows the binding site residues RMSF (red and blue rectangles) and the highly flexible residues (yellow rectangles).

The MD simulation resulted in a trajectory file of 14,285 protein frames. Considering all the frames in an ensemble docking approach is certainly time-consuming, and redundant outcome can be produced due to high frame similarities in some regions. Therefore, partitioning MD frames into clusters of similar frames addresses the problem of their large number. TTClust program [55] was used for this clustering task, where it calculates the average RMSD between frames and cluster similar frames together. In the first run, the RMSD was assigned to be calculated on the backbone atoms (Fig. S2 in the Supplementary Material). In the second run, heavy atoms of the flexible loop (T265-H272) were assigned for RMSD calculations. It resulted in 5 clusters. Nonetheless, to provide higher conformational sampling, a larger arbitrary number of clusters were pursued. Hence, a third run was carried out with cluster number of 10 as a default output. The results of the third clustering approach are shown in Fig. 4 .

Fig. 4

TTClust results applied to the frames of 1micro-second MD. The heavy atoms of T265-H272 residues are assigned for the calculation of average RMSD between frames with 10 as a default number of output clusters. (A) Hierarchical clustering dendrogram. (B) 2D projection plot of the relative distances between clusters calculated based on the representative structure of each cluster showing the minimum and maximum RMSD and minimum and maximum average RMSD within clusters (spread) (C) Timeline barplot showing the frame distribution along the simulation time. (D) Heatmap showing the RMSD distance matrix between frames. (E) Histogram of number of frames in each cluster. The color code of the cluster is the same in the dendrogram, 2D projection plot, timeline barplot and histogram (a-c, e).

Eight out of the ten produced clusters were present in the first half of the simulation time (0–500 ns). This helped us to sample the first part of the simulation rigorously. It was observed that clusters 1, 2, 3 and 4 are the closest to the X-ray structures (PDB ID: 6W9C and 7JRN) with respect to the loop conformation, especially Tyr268 and Gln269 key binding residues, as inspected in Fig. 5 .

Fig. 5

(A) The ten cluster representative structures superimposed to 6W9C as reference structure (ribbon) with Tyr268 and Gln269 key binding residues (sticks). (B) Depictions of the flexible loop with Tyr268 and Gln269 in sticks of different clusters. (C) Pairwise root mean square deviation (RMSD) between the C1–C10 (numbers 1–10) and X-ray structures PDB ID: 6W9C and 7JRN (numbers 11–12). Blue colors indicate close conformations (low RMSD values), while red colors indicate distinct conformations (high RMSD values).

Virtual screening of DrugBank

Twelve protein structures comprising the ten MD clusters representatives, the starting X-ray structure (PDB ID: 6W9C) and the X-ray structure of PLpro-small molecule complex (PDB ID: 7JRN) were utilized for ensemble docking for screening the DrugBank database [56]. Such set up is to represent various states of the conformation space of PLpro. This approach would enable us to extract a consensus list of drugs that is able to target PLpro in different conformations simultaneously. Based on the results of our previous benchmarking study [13], we applied FRED docking tool for this ensemble screening approach since it showed promising screening performance in recognizing bioactive molecules in a pool of challenging decoys for both the bound and unbound states of PLpro [13,31]. Also, we biased our consensus selection for drugs enriched in the newly introduced X-ray structure complexed with a small molecule (PDB ID: 7JRN), because it showed the best benchmarking results among other PLpro X-ray structures [13]. Top enriched 1% of the screened DrugBank database [56] (25 compounds) were retrieved for each docking run and concatenated together to extract commonly enriched drugs in a consensus list. Drugs with 1% enrichment and in consensus rank for 5, 6, 7 and 8 different conformations are listed in Table 1 . As a consensus, ten drugs were shown to commonly target different conformations of PLpro.

Table 1

Consensus list of the top 10 out of the best enriched 1% FDA approved (DrugBank) against the 12 conformations of PLpro.

Drug	Targeted MD cluster number/X-ray structure of PLpro (Docking score)a	Number of PLpro conformations	Structure
Amifostineb	C2 (−8.14)	Five conformations	Image 2
	C5 (−7.21)
	C7 (−8.20)
	C8 (−7.34)
	6w9c (−8.72)
Benserazide	C3 (−10.28)	Five conformations	Image 3
	C7 (−8.02)
	C10 (−7.97)
	7jrn (−11.15)
	6w9c (−8.67)
Carglumic acid	C2 (−8.87)	Five conformations	Image 4
	C7 (−7.89)
	C8 (−6.74)
	C9 (−7.62)
	6w9c (−9.52)
Dobutamine	C1 (−9.84)	Five conformations	Image 5
	C2 (−8.79)
	C3 (−9.92)
	7jrn (−10.66)
	6w9c (−9.34)
Pemirolast	C2 (−8.60)	Five conformations	Image 6
	C3 (−10.96)
	C4 (−9.38)
	C5 (−9.04)
	C7 (−8.90)
Penciclovir	C1 (−9.84)	Five conformations	Image 7
	C2 (−8.86)
	C5 (−8.41)
	C7 (−7.85)
	6w9c (−8.93)
Masoprocol	C1 (−10.93)	Six conformations	Image 8
	C2 (−8.52)
	C3 (−10.79)
	C4 (−9.31)
	C5 (−7.37)
	7jrn (−10.38)
5-O-phosphono-alpha-d-ribofuranosyl diphosphate	C2 (−9.26)	Seven conformations	Image 9
	C3 (−9.53)
	C4 (−11.52)
	C6 (−8.43)
	C7 (−8.77)
	C10 (−6.86)
	6w9c (−9.02)
Pyridoxal phosphate	C1 (−10.80)	Seven conformations	Image 10
	C2 (−10.83)
	C4 (−9.57)
	C5 (−7.42)
	C6 (−8.77)
	C7 (−8.58)
	C10 (−7.03)
Adenosine phosphate	C1 (−10.14)	Eight conformations	Image 11
	C2 (−8.59)
	C3 (−9.90)
	C4 (−9.27)
	C5 (−8.07)
	C6 (−8.11)
	C7 (−7.89)
	6w9c (−8.91)

The docking score is based on FRED docking outcome.

Prodrug.

Three drugs were enriched for the X-ray structure complexed with small molecule (PDB ID: 7JRN), namely Benserazide, Dobutamine and Masoprocol, with docking scores of −11.15, −10.66 and −10.38, respectively. Interestingly, their docking score is observed to be the best in the X-ray structure (7JRN) compared to other conformations, therefore, we inspected their postulated binding poses in such X-ray structure, as shown in Fig. 6 and Fig. 7 .

Fig. 6

(A) Benserazide docking pose (blue sticks) and the co-crystallized naphthalenyl benzamide derivative (violet sticks) showing the key H-bonding interactions with Asp164 side chain and Gln 269 backbone in orange and black lines, respectively. (B) Metoclopramide, Isocarboxazid, Procainamide and Remoxipride docking poses (white sticks) displaying the two key H-bonding interactions with Asp164 and Gln269 residues (orange lines). Metoclopramide, Isocarboxazid and Procainamide possess extra H-bonds with Lys 157, Tyr273 and Leu162, respectively (black lines). The 2D interaction representation for Benserazide can be found in the Supplementary Material in Fig. S3.

Fig. 7

(A) Dobutamine docking pose illustrating H-bonding interaction with Tyr273, pi-pi stacking with Tyr264 and pi-anion interaction with Asp164. (B) Masoprocol docking pose revealing H-bonding interactions with Leu162, Arg166 and Gln269, and pi-pi stacking with Tyr268, pi-anion with Asp164, pi-sigma with Tyr264, and hydrophobic interactions with both Pro248 and Tyr273. The 2D interaction representation for the ligand-protein interactions can be found in the Supplementary Material in Fig. S4 and Fig. S5.

Benserazide appeared to be enriched for 5 conformations involving both X-ray structures (PDB IDs: 7JRN and 6W9C). It is an inhibitor of levodopa (l-DOPA) decarboxylation that acts peripherally and does not cross the blood-brain barrier. It is used in combination with the dopamine precursor l-DOPA in Parkinson's disease management to prevent its decarboxylation and facilitate its blood-brain barrier crossing. Upon investigation of Benserazide postulated binding mode and comparing it with the co-crystallized ligand in (PDB ID: 7JRN), two hydrogen bonds with Asp164 side chain and Gln269 backbone atoms are conserved in both poses. These two hydrogen bonds are formed by the hydrazide group of Benserazide and the amide group of the co-crystallized ligand. Two additional hydrogen bonds with Tyr268 and Tyr273 side chain atoms are formed with Benserazide. Additionally, the T-shaped pi-pi stack between the co-crystallized ligand naphthyl ring and Tyr 268, and its pi-alkyl interaction with Pro248, are also conserved with the trihydroxy phenyl ring of Benserazide (Fig. 6A).

Interestingly, the top-enriched four drugs for the X-ray PLpro structure (PDB ID: 7JRN) showed the best docking scores globally across all docking runs. These drugs are Metoclopramide, Isocarboxazid, Procainamide and Remoxipride with docking scores of −12.70, −12.18, −12.02 and −11.24, respectively. Amide group is observed to be shared for all of them conserving two main H-bonding interactions with Asp164 and the key Gln269. This observation can be valuable in prospective drug design projects targeting PLpro (Fig. 6B).

Dobutamine is observed to be enriched for 5 conformations including the X-ray structures (PDB IDs: 7JRN and 6W9C). It is a β1 agonist used in the treatment of congestive heart failure. As shown in Fig. 7, its dihydroxy phenyl group maintains a pi-pi stacking with Tyr264 phenyl group, and a hydrogen bond with Tyr273. Additionally, the phenol substituent forms a pi-anion interaction with Asp164 amino acid (Fig. 7A).

Masoprocol is enriched for 6 conformations including the X-ray structure (PDB ID: 7JRN). It is a lipoxygenase inhibitor used as antineoplastic drug for precancerous skin growths. It possesses a structure similarity with Dobutamine, hence, its docking pose reproduced similar binding pi-anion interactions with Asp164 by both phenyl rings. Two additional hydrogen bonds were also predicted between both Leu162 backbone and Gln269 side chain and the dihydroxy phenyl ring, which also forms a T-shaped pi-pi stack with Tyr268. The other dihydroxyphenyl ring also forms a hydrogen bond with Arg166 side chain.

The proposed docking poses of the three drugs (benserazide, dobutamine and masoprocol) were subjected to molecular dynamics simulations for 50 ns to test their stability inside the binding site. The values of Rg, RMSD, and RMSF indicate stable simulation for these ligand-PLpro complexed forms (Fig. S6 in the Supplementary Material).

Fig. 8 A displays the RMSD of the three ligands complexed with PLpro through the 50 ns simulation time. The benserazide-PLpro complex exhibits high RMSD values and obvious fluctuations reflecting its relative weak binding to PLpro compared to the other two drugs. On the other hand, dobutamine and masoprocol showed substantial better binding and stability in the binding site mirrored from their low RMSD values during the 50 ns simulation time. In general, the dynamics of the complexes converged after 10 ns of simulation, implying the idea that the structural changes present in the complexes of dobutamine or masoprocol with PLpro converged to a stable structure. Fig. 8B depicts the number of hydrogen bonds formed between each ligand and the protein. The hydrogen bonds of dobutamine with the binding site residues were the most conserved compared to the other two drugs during the MD simulation, while both benserazide and masoprocol were partially conserved during the simulation. Collectively, from the results of these three drug-PLpro MD runs, both dobutamine and masoprocol are proposed for inhibition of PLpro, with benserazide having the lowest probability of PLpro inhibition due to its relative unstable binding.

Fig. 8

(A) RMSD of the three ligands showing benserazide with the highest RMSD, (benserazide: green, dobutamine: blue and masoprocol: yellow). (B) number of hydrogen bonds formed between each ligand and the protein with dobutamine having the highest number of hydrogen bonds and highest proposed stability.

It is noteworthy mentioning that seven other drugs were predicted to target different conformations of PLpro (Table 1), namely: Amifostine, Carglumic, Pemirolast, Penciclovir, 5-O-phosphono-alpha-d-ribofuranosyl diphosphate, Pyridoxal and Adenosine phosphate. Amifostine is a prodrug which is dephosphorylated by alkaline dephosphatases to yield the free thiol active drug. It is used as cytoprotective in cancer chemotherapy and radiotherapy. Carglumic acid is N-acetylglutamate (NAG) synthetic analog which is used in treatment of hyperammonaemia in NAG-deficient patients through activation of carbamoyl phosphate synthetase 1 (CPS1) liver enzyme to start the urea cycle and convert ammonia into urea. Pemirolast is a mast cell stabilizer used in the treatment of allergic conjunctivitis-induced eye itching. Penciclovir is an antiviral drug used for management of herpes simplex virus types 1 (HSV-1) and 2 (HSV-2). It competes with deoxyguanosine triphosphate, hence, inhibiting DNA polymerase responsible for viral replication. 5-O-phosphono-alpha-d-ribofuranosyl diphosphate is a carbohydrate derivative possessing a pentose ring with phosphate substituents. It is considered the key substance in the biosynthesis of histidine, tryptophan, and purine and pyrimidine nucleotides. Pyridoxal phosphate is the active form of vitamin B6. It serves as a coenzyme for synthesis of amino acids and neurotransmitters such as serotonin and norepinephrine. It is used as a dietary supplement to treat dietary imbalance. Adenosine phosphate is also used as a dietary supplement for treating dietary shortage. The clinical indications of the proposed drugs are tabulated in Table T2 in the Supplementary Material.

Methodology

Molecular dynamics

The PLpro apo crystal structure (PDB ID: 6W9C) was subjected to 1 μs molecular dynamics simulation at NPT ensemble. The simulation system was prepared using OpenMM Setup utility [57]. AMBER14 [58] force field was used for the protein parameters and TIP3P-FB water model [59] was used for water molecules. A preliminary clean up step was done for the protein PDB file by adding missing OXT atom for LYS315 amino acid. Hydrogen atoms were added at pH value of 7.0. The protein was immersed in a cubic water box with a padding distance of 1.3 nm along each axis of the protein molecule. NaCl ions were added with ionic strength of 0.1 Molar to neutralize the model. The final simulation system ended up with 88,222 atoms. The PLpro crystal structure (PDB ID: 7JRN) used for MD simulation with benserazide, dobutamine and masoprocol was treated by the same parameters as the apo protein. Generalized Amber Force Field (GAFF) [60] parameters were used for the three drugs. Both GAFF parameters and ELF10 charges were generated on the fly by the OpenEye toolkits 2019.4.1 [61]. Particle Mesh Ewald (PME) method [62] was used for the calculation of long-range electrostatic interactions with a 0.0005 error tolerance for truncating the Ewald summation. A 1.2 nm cut off was used for both PME direct space interactions and Lennard-Jones interactions. The length of all bonds that involve a hydrogen atom, and the water molecules, involving their bond length and angles, were constrained. Langevin integrator was utilized for integrating force equations with step size of 2 ns, temperature of 310 K and friction coefficient of 1 ps−1. Mote Carlo barostat [63] was used for pressure coupling with assigning 310 K as constant temperature value, 1 atm as constant pressure value and 25 steps as pressure update frequency. The dynamics simulations were carried out using OpenMM python application layer [57]. We started the simulation with a minimization step with 10 kJ/mol as energy convergence criteria, then, 3,000,000 equilibration steps (6 ns). Finally, the simulation was carried out for 500 million steps (1 μs). A trajectory snapshot was reported every 35,000 steps (70 ps), resulting in a final trajectory file of 14,285 total frames. For the three ligand-PLpro MD runs with benserazide, dobutamine and masoprocol, equilibration was run for 500,000 steps (1 ns), the production was run for 25 million steps (50 ns), and trajectory reporting was set to one snapshot per 1,000 steps resulting in a trajectory file of 25,000 snapshots. RMSD, RMSF and radius of gyration was calculated out using ProDy python library [64,65].

Trajectory clustering

As a first step of trajectory clustering, the trajectory frames were aligned on the protein backbone atoms of the starting structure. Then, pairwise RMSD between frames was calculated and stored in a matrix, for backbone in the first clustering run, and heavy atoms of residues T265-H272 in the second and third runs. This matrix was then used in hierarchical clustering of frames implementing the ward method in calculating the distance between the newly formed clusters. The third clustering run was carried out assigning 10 clusters as output cluster numbers. Clustering was carried out using TTClust program [55].

DrugBank virtual screening

DrugBank database [66] was downloaded from www.drugbank.ca (DrugBank - release March 2020) with 2454 ligands. It was prepared by MOE [67] with retrieving one conformer of one protonation state at pH 7.0 implementing Amber: 10EHT forcefield. Multi-conformer file generation for the database was carried out using OMEGA conformer generation module of OpenEye software [68,69] with a default number of 200 conformers generated for each ligand. The 10 representative protein structures, along with the downloaded protein structures (6W9C and 7JRN), were prepared using MakeReceptor GUI module of OpenEye software [70]. With respect to the 10 representative structures, possible cavities on the protein surface was detected by cavity detection utility of MakeReceptor implementing a molecular probe. Then, the search box dimensions were assigned according to the detected cavity volume of the binding site surrounding the binding site amino acids: L162-E167, P247–P248, Y264–Y273 and T301. Virtual screening was carried out using FRED module of OpenEye software [53,54], by docking the multi-conformer DrugBank database file into the prepared protein structures. Protein and protein-ligand 3D interaction figures were generated using UCSF Chimera [71].

Conclusion

Drug repurposing of DrugBank database was carried out by adopting ensemble docking approach for multiple conformations of PLpro SARS-CoV-2 molecular target. Ten PLpro different conformations sampled from 1 μs MD simulation, in addition to two reported PDB PLpro structures, were used for ensemble virtual screening. Based on consensus ranking of the top 1% of the screening hits, three drugs showed to target different conformations with high docking scores and interaction with Tyr268 and Gln269 key binding residues, namely: Benserazide, Dobutamine and Masoprocol. Further MD simulations for these drugs complexed with PLpro indicated higher stability and binding of dobutamine and masoprocol inside the binding site compared to benserazide. Moreover, Metoclopramide, Isocarboxazid, Procainamide and Remoxipride showed the best docking scores over all docking runs combined despite their lower consensus over different conformations. These proposed drugs can be biologically evaluated as potential SARS-CoV-2 treatment candidates. Globally, this protocol can help in identifying drugs for repositioning via targeting multiple conformations of a crucial target for the rapidly evolving COVID-19 pandemic.

Author contributions

TMI conceived the idea and supervised the project. MII caried out all experiments. MII and TMI wrote the manuscript. All authors have given approval to the final version of the manuscript. All authors agree to be accountable for the content of the work.

Declaration of competing interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.