In silico study to identify novel potential thiadiazole-based molecules as anti-Covid-19 candidates by hierarchical virtual screening and molecular dynamics simulations

In the present study, a new category of 1,3,4-thiadiazoles was developed by submitting methyl 2-(4-hydroxy-3-methoxybenzylidene) hydrazine-1-carbodithioate to react with the appropriate hydrazonoyl halides in presence of few drops of diisopropyl ethyl amine. The chemical structures of the newly synthesized derivatives were inferred by means of their micro-analytical and spectral data. Utilizing combined molecular docking and molecular dynamics techniques, the binding affinities and features of the synthesized compounds were evaluated against four SARS-CoV-2 target enzymes, namely, main protease (Mpro), papain-like protease (PLpro), RNA-dependent RNA polymerase (RdRp), and receptor-binding domain (RBD) of the spike protein. Compound 7 demonstrated promising binding affinities with the target enzymes Mpro, PLpro, RdRp, and RBD with docking scores of −11.4, −9.4, −8.2, and −6.8 kcal/mol, respectively. In addition, compound 7 exhibited MM-GBSA//100 ns MD docking score of −35.9 kcal/mol against Mpro. Structural and energetic analyses revealed the stability of the 7-Mpro complex over 100 ns MD simulations. In addition, compound 7 obeyed Lipinski’s rule of five, as it has acceptable absorption, distribution, and oral bioavailability inside the body. Therefore, compound 7 is considered as a promising starting point for designing potential therapeutic agents against Covid-19. Supplementary Information The online version contains supplementary material available at 10.1007/s11224-022-01985-1.


Introduction
A new strain for SARS-CoV-1 identified recently as SARS-CoV-2 in late December 2019 resulted in serious physical and psychological damages to the human health; a massive outbreak initially in Wuhan, China, and spread rapidly in different nations around the global in a short time [1,2]. The World Health Organization (WHO) declared this highly infectious respiratory disease Covid-19 as a pandemic [3]. It is acceptable to think that a sufficient understanding of SARS-CoV-2 and the full clinical picture of the resulting Covid-19 disease will take some time. However, the first detected clinical sign of Covid-19 was pneumonia [4]. Recently, asymptomic infections and gastrrointestinal symptoms were also reported especially among young children [5,6]. Pneumonia mostly appeared in the second or third week of the infection. Decreased oxygen saturation, blood gas deviations, and changes visible through chest X-rays are prominent signs of viral pneumonia. In addition, lymphopenia documented to be common, and inflammatory markers (Proinflammatory cytokines and CRP) are elevated. Consequently, investigation of anti-Covid-19 therapeutic agents became an urgent demand and attracted more interest recently owing to the lack of specific drugs for the treatment of Covid-19 [7,8]. Nevertheless, several existing drugs are available only to overcome the clinical symptoms of  On the other hand, 1,3,4-thiadiazole moieties have been reported for their pharmaceutical properties. Many antiviral drugs like acetazolamide, besaglybuzole (glybuzole), and furidiazine (triafur) were reported to append the 1,3,4-thiadiazole in their constructors [9][10][11][12][13].
Urgent needs for develoment novel anti-Covid19 agents have directed us to synthesize some new bioactive heterocyclic molecules. In the present study, we aimed to identify potential Covid-19 inhibitors through a computer-based molecular docking and molecular dynamics techniques [14][15][16]. In addition, ADMET (absorption, distribution, metabolic, excretion, and toxicity) and pharmacokinetics parameters of the prepared ligand molecules were performed to identify their drug-likeness properties [17].
The chemical structures of all newly prepared molecules are affirmed by spectral and elemental data. For instance, IR spectrum of the target molecule 7 revealed a strong broad absorption band at v 3337 assigned for NH group. Additionally, it showed a strong sharp absorption band at v 1681 attributed to the carbonyl group. Meanwhile, 1 H-NMR spectrum exhibited singlet signal at δ 3.83 ppm represented the methoxy group along with multiplet signal at δ 6.86-7.85 ppm for aromatic protons. Also, it revealed doublet signal at δ 7.75 ppm attributed to the aromatic hydrogen and doublet signal at δ 8.15 ppm represented the aromatic hydrogen. Moreover, it showed three singlet signals at δ 8. 36, 9.65, and 10.68 ppm for CH = N, OH and NH, respectively as illustrated in Fig. 1. Figure 2 revealed the significant signals of the target molecule 7 which confirmed the formation of the compound.
The chemical composition of the target molecule 7 was affirmed also by the mass spectrum (m/z 445) [M + ], which agrees with its molecular formula C 23 H 19 N 5 O 3 S, as represented in Fig. 3.

Molecular docking calculations
The molecular docking technique [11,18] was used to predict the binding modes and affinities of the newly synthesized compounds with SARS-CoV-2 targets M pro , PL pro , RdRp, and RBD of S-protein. The predicted docking scores are tabulated in Table 1. 2D (two dimensional) and 3D (three dimensional) representations of binding modes of best docked compound 7 inside the active site of M pro , PL pro , RdRp, and RBD are displayed in Fig. 4. The representations of the rest docked compounds against the targets are shown as Figs. S1-S4, respectively, in the Supplementary file section.

Scheme 1 Synthetic procedures of the desired 1,3,4-thiadiazoles 3-7
It is observed from the data in Table 1 that compound 7 exhibited binding affinities against all the selected targets better than the reference drug (Darunavir) (Fig. 4).
In addition, most of the synthesized compounds demonstrated promising binding affinities against M pro with binding energies ranged from −11.4 to −6.4 kcal/mol. The high docking scores of the studied compounds with M pro would be returned to their ability to form hydrogen bonds, hydrophobic and van der Waals interactions with the amino acid residues of active sites (Fig. 5).
Compared to M pro , the examined compounds showed relatively weak binding affinities with PL pro , RdRp, and RBD, with docking scores ranged from −9.4 to −5.6, −8.2 to −5.2, and −6.8 to −5.3 kcal/mol, respectively. Interestingly, compound 7 displayed the highest binding affinity against all targets with docking scores of −11.4, −9.4, −8.2, and −6.8 kcal/mol, respectively. The promising binding affinity of 7 against the target M pro is attributed to its ability to exhibit two hydrogen bond interactions with Thr24, and Ser144 at 2.17 and 2.91 Å, respectively (Fig. 4). Besides, inspecting the binding modes of compound 7 with PL pro , RdRp, and RBD unveiled its potentiality to form hydrogen bonds and Pi-stacking interactions, as presented in Fig. 4. The ligand molecule 7 docked with PL pro amino acid residues (Lys218, Tyr251, and Phe258) through pi-stacked interactions at distances 4.27, 5.51, and 5.58 Å, respectively. Further, compound 7 interacted with the residues Arg116 and Phe35 of the target RdRp through pi-stacked interactions at distances 5.89, and 3.92 Å, respectively. Finally, it docked with RBD through two hydrogen bonds and one pi-pi interaction with the amino acid residues Tyr385, Arg393, and Phe40.

Molecular dynamics (MD) simulations
Towards more reliable binding affinities, the molecular dynamics simulations were performed for all synthesized compounds in complex with SARS-CoV-2 targets. The binding energies (ΔG binding ) were then calculated using the molecular mechanics-generalized born surface area (MM-GBSA) approach based on the collected snapshots for M pro , PL pro , RdRp, and RBD of spike protein over the production stage of 25 ns. The calculated MM-GBSA binding energies are listed in Table 2.
It is apparent from Table 2 that the examined compounds with M pro showed higher binding affinities over 25 ns MD simulations than with PL pro , RdRp, and RBD. The calculated MM-GBSA binding energies were in line with the predicted docking scores, demonstrating the high potency of the examined ligand molecules with M pro over the other SARS-CoV-2 targets.
Among the examined compounds, 7 exhibited the lowest binding energy with M pro with a ΔG binding value of −39.2 kcal/mol. Moreover, it showed weak binding energies of −15.8, −15.5, and −14.9 kcal/mol with PL pro , RdRp, and RBD, respectively. These results declared the selectivity of the compound 7 towards M pro over PL pro , RdRp, and RBD. However, compound 4 demonstrated the lowest binding with PL pro , RdRp, and RBD with a ΔG binding value of −27.9, −26.1, and −21.7 kcal/mol,   (Table 2). MD simulation for compound 7 in complex with M pro and compound 4 complexed with PL pro , RdRp, and RBD were then elongated to 100 ns. Additionally, the corresponding MM-GBSA binding energy was calculated and was compared to the reference drug Darunavir (Fig. 6). MM-GBSA binding energy of compound 7-M pro , compound 4-PL pro , compound 4-RdRp, and compound 4-RBD complexes was decomposed to explore the predominant interactions between the investigated inhibitors and target. According to the data, it was found that the docking energy was calculated by E vdw interactions with an average value of −47.2, −38.4, and −36.8 kcal/mol for investigated inhibitor with M pro , PL pro , and RdRp, respectively (Fig. 6). For compound 4 complexed with RdRp, the docking energy was dominated by E ele interactions with an average value of −65.5 kcal/mol which was three times higher than that of Lopinavir and curcumin, with an average value of −42.5 kcal/mol (Fig. 6). Together these results demonstrated the promising binding affinity of compounds 7 and 4 with SARS-CoV-2 targets.

Post-dynamics analysis
The interaction nature and stability of compound 7 and Darunavir within the active site of M pro was estimated using structural and energetic analyses. Structural and energetic analyses including energy per-frame, centre-of-mass distance (CoM), root-mean-square deviation (RMSD), and root-mean-square fluctuation (RMSF) were performed over 100 ns MD simulations.

Docking energy per frame
The stability of compound 7 and Darunavir in complex with SARS-CoV-2 M pro was estimated using the correlation between the binding energy per-frame and time. MM-GBSA binding energy was subsequently evaluated per-frame for the most promising compound with each target and displayed in Fig. 7. The most interesting aspect of this graph is the overall stability of two identified compounds over 100 ns MD simulations with average values of −35.9, and −34.8 for compound 7-M pro , and darunavir-M pro complexes, respectively.

Center-of-mass distance
Interestingly, investigating the center-of-mass (CoM) distance between the compound 7, and Darunavir and the residue Glu166 through the 100 ns MD simulations would reflect a strong indication of the high stability of the identified compounds inside the M pro active site. The CoM distances were inspected over the 100 ns MD simulations and  Fig. 8. What stands out in Fig. 8 is the average CoM distance between the identified compounds and the key amino acid residue Glu166 was approximately constant, with average CoM distances of 5.7, and 12.1 Å, respectively. The current data revealed that compound 7 bound more tightly to the M pro complex compared to Darunavir.

Root-mean-square deviation
The structural changes of 7-M pro and darunavir-M pro complexes were evaluated using the root-mean-square deviation

ADMET and drug-likeness properties of the molecules
ADMET studies of the prepared molecules exhibited that they have acceptable absorption and distribution properties in the range of (91-80-98.76%) and (0.53-0.70), respectively. The physiochemical properties of the compounds exhibited acceptable values, as they have molecular weights and partition coefficients in the range of (256.35-445.49 g/ mol), and < 5, respectively. Moreover, the molecules have no toxicity and carcinogenicity. All tested compounds showed Fig. 8 Centre-of-mass (CoM) distances between 7 (in black), and darunavir (in red) with GLU166 of SARS-CoV-2 M. pro throughout 100 ns MD simulations Fig. 9 Root-mean-squaredeviation (RMSD) of the backbone atoms from the initial structure for 7 (in black) and darunavir (in red) with the SARS-CoV-2 main protease (M. pro ) over 100 ns MD simulations good oral bioavailability within the body as they obeyed Lipinski's rule of five (Table 3).

Conclusion
In this study, a new series of 1,3,4-thiadiazole derivatives was synthesized, characterized, and theoretically evaluated as Covid-19 inhibitors against four SARS-CoV-2 targets namely, main protease (M pro ), papain-like protease (PL pro ), RNAdependent RNA polymerase (RdRp), and receptor-binding domain (RBD) of the spike protein. The molecular docking studies and molecular dynamics simulations exhibited the promising binding affinity of compound 7 with all targets. Therefore, it could be select as promising chemical moiety for designing of future inhibitors as anti-Covid-19 agents.

Instrumentation
All melting points were uncorrected and measured using electrothermal device. The IR spectra were recorded (KBr discs) using Shimadzu FT-IR 8201 PC spectrophotometer. 1 H-and 13 C-NMR spectra were recorded in (CD3) 2 SO solutions on a BRUKER 500 FT-NMR system spectrometer, and chemical shifts are expressed in ppm units using TMS as an internal reference. Mass spectra were recorded on a GC-MS QP1000 EX Shimadzu. Elemental analyses were carried out at the Microanalytical Center of Cairo University.

Synthetic procedures of the target molecules (3-7)
A mixture of compound 2 (1.28 gm, 5 mmol) and the selected derivative of the hydrazonoyl halides (5 mmol) and 2-3 drops of DIPEA as a catalyst, were ground well in an open mortar with a pestle for 5-7 min. at RT till the mixture turned into melt. The grinding was continued for approximately 5-10 min, and the reaction was monitored by TLC. The solid was collected and washed with (water/ethanol) the recrystallized from the proper solvent to give the desired derivatives 3-7, respectively.

Inhibitor preparation
The chemical structures of the synthesized compounds were manually constructed, and their 3D structures were generated using Open Babel 2.4.1 tool [25][26][27]. All ligand molecules were then energetically minimized using the CHARMM Force Field [28].

Molecular docking
Molecular docking calculations [10,12,14,18,24,[29][30][31][32][33] were carried out using PyRx -virtual screening software [34]. The pdbqt files of M pro , PL pro , RdRp, and RBD of S-protein targets were prepared according to PyRx protocol. The docking algorithms were conserved to their default values, except the number of genetic algorithms (GA) run and the maximum number of energy evaluation (eval). In the current study, GA and eval were set to 250 and 25,000,000, respectively. The docking grid was set to 25 Å × 25 Å × 25 Å with a spacing value of 0.375 Å [17,35]. The grid center was positioned at the center of the active site of M pro , PL pro , RdRp, and RBD of S-protein. The partial atomic charges of the examined compounds were estimated using the Gasteiger method [36]. The prediction of binding modes for each compound was handled using the built-in clustering analysis with an RMSD tolerance of 1.0 Å. Also, the lowest energy conformation from the largest cluster was picked out as a representative binding pose.

Molecular dynamics simulations
Molecular dynamics simulations were performed for the examined compounds in complex with the four studied SARS-CoV-2 targets using AMBER16 software. In MD simulations, the General AMBER force field (GAFF2) [37] and AMBER force field 14SB [38] were employed to describe the studied compounds and SARS-CoV-2 targets, respectively. The atomic partial charges of the examined compounds were calculated using the restrained electrostatic potential (RESP) approach at the HF/6 −31G* level with the help of Gaussian software [39]. Prior to RESP charge calculations, the studied compounds were first geometrically optimized at the B3LYP/6-31G* level of theory. The docked compoundtarget complexes were solvated in a cubic water box with 15 Å distances between the edges of the box and any atom of compound or compound-target complexes. The solvated compound-target systems were subsequently energy minimized for 5000 steps, gently annealed from 0 to 300 K over 50 ps, and equilibrated for 1 ns. The equilibrated systems were then simulated for 100 ns using periodic boundary conditions and NPT ensemble. The non-bonded cut-off distance was placed at 12 Å, and particle mesh Ewald (PME) method was applied to process long-range electrostatic interactions. The Langevin dynamics with the collision frequency gamma_ln set to 1.0 was used to conserve the temperature of the examined systems at 298 K. The pressure of the system was controlled using Berendsen barostat with a relaxation time of 2 ps. A time step of 2 fs and the SHAKE option to constrain all bonds involving hydrogen atoms were utilized. Coordinates and energy values were collected every 10 ps over the production stage for binding energy calculations and post-dynamics analyses. All MD simulations were conducted with the GPU of pmemd (pmemd.cuda) in AMBER16 on the CompChem GPU/CPU cluster (hpc.compchem.net). 2D and 3D visualization of the compound-targets interactions were performed using the Discovery studio software.

MM-GBSA binding energy
The binding free energies of the examined compounds with SARS-CoV-2 targets were estimated using molecular mechanical-generalized Born surface area (MM-GBSA) approach [40,41]. For MM-GBSA calculations, uncorrelated snapshots were collected every 10 ps over the production stage. The MM-GBSA binding energy ( ΔG binding ) can be conceptually summarized as: where the energy term (G) is estimated as: where E vdw and E ele are the van der Waals and electrostatic energies, respectively. G GB is the electrostatic solvation free energy calculated from the generalized Born equation, and G SA is the nonpolar contribution to the solvation free energy from the solvent-accessible surface area (SASA). Solute entropy contributions to binding energies were neglected, and a single-trajectory approach was employed, in which the coordinates of each molecule, receptor, and complex were isolated from a single trajectory.

ADMET analysis
The freely accessible online softwares such as admetSAR, SwissADME, and Mol inspiration are used to predict ADMET and drug-likeness properties of compounds.
Author contribution H.R.M.R. and A.H.A. made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis, and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Data availability
The data that support the findings of this study are included within the article and supplementary file.

Conflict of interest
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.