Synthesis, Molecular Docking, and Dynamic Simulation Targeting Main Protease (Mpro) of New, Thiazole Clubbed Pyridine Scaffolds as Potential COVID-19 Inhibitors

Many biological activities of pyridine and thiazole derivatives have been reported, including antiviral activity and, more recently, as COVID-19 inhibitors. Thus, in this paper, we designed, synthesized, and characterized a novel series of N-aminothiazole-hydrazineethyl-pyridines, beginning with a N′-(1-(pyridine-3-yl)ethylidene)hydrazinecarbothiohydrazide derivative and various hydrazonoyl chlorides and phenacyl bromides. Their Schiff bases were prepared from the condensation of N-aminothiazole derivatives with 4-methoxybenzaldehyde. FTIR, MS, NMR, and elemental studies were used to identify new products. The binding energy for non-bonding interactions between the ligand (studied compounds) and receptor was determined using molecular docking against the SARS-CoV-2 main protease (PDB code: 6LU7). Finally, the best docked pose with highest binding energy (8a = −8.6 kcal/mol) was selected for further molecular dynamics (MD) simulation studies to verify the outcomes and comprehend the thermodynamic properties of the binding. Through additional in vitro and in vivo research on the newly synthesized chemicals, it is envisaged that the achieved results will represent a significant advancement in the fight against COVID-19.


Introduction
Recently, antiviral, chemotherapeutic drugs are ineffective in clinic settings. This is a result of the development of a number of significant viral infections, which has resulted in widespread human disease and mortality. Coronaviruses (CoV) are a large group of viruses that affect a wide variety of animals. They have caused serious and deadly respiratory infections in both humans and animals, such as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and the Middle East respiratory syndrome coronavirus (MERS-CoV) [1][2][3][4][5]. The World Health Organization (WHO) reported the development of coronavirus disease 2019 . In terms of the persons afflicted and the geographic scope of the outbreak, COVID-19 has significantly exceeded SARS and MERS [6]. Therefore, new antiviral

Chemistry
Synthesis of thiazole derivatives 6a-e and 13a-c.

Chemistry
Synthesis of thiazole derivatives 6a-e and 13a-c.
Catalytic amounts of TEA were added into a solution of N'-(1-(pyridin-3-yl) ethylidene)hydrazine-carbothiohydrazide (3) (0.209 g, 1 mmol) and the appropriate hydrazonoyl chlorides 4a-e or α-bromoketones 11a-c (1 mmol for each) in DMF (20 mL). The reaction mixture was refluxed for 3-6 h. Finally, the formed precipitate was isolated and recrystallized from the suitable solvent to yield the compounds 6a-e or 13a-c, respectively. A catalytic amount of HCl concentration was added to a solution of 4-methoxybenzaldehyde (7) (1.36 g, 10 mmol), and the appropriate 8a,d or 13a-c (1 mmol for each) was added to DMF (20 mL). The reaction mixture was refluxed for 2-4 h. Finally, the formed precipitate was recrystallized from the suitable solvent to yield compounds 6a-e or 13a-c, respectively.
Physical and spectral data of all synthesized compounds 6a-e, 13a-c, 8a,d and 14a-c are frond in the supporting information file.

Docking Method
The newly synthesized compounds were subjected to docking tests using the molecular operating environment 2019.012 suite (Montreal, QC, Canada) [37] to ascertain how well they bound and to propose their mechanism of action as SARS-CoV-2 Mpro inhibitors in comparison to the co-crystallized inhibitor (N3), which was used as a reference standard. Energy was minimized and a partial charge was added to the freshly synthesized molecules inside the MOE window [38,39]. The synthesized compounds were then combined with (N3) in one database and stored as an MDB file that could be transferred into the ligand icon during the binding stage. The Protein Data Bank provided the X-ray crystallography target, M pro , of SARS-CoV-2 (PDB code: 6LU7) [40]. Additionally, it was ready for docking by carefully following the previously detailed methods [41,42]. Furthermore, the downloaded protein was energy-reduced, 3D-hydrogen-loaded, and error-corrected [43,44]. The newly created molecules were substituted for the ligand location in a general docking approach. The co-crystallized ligand site was chosen as the docking site after adjusting the default program settings that were provided [45]. In a nutshell, the dummy atoms method was used to select the docking point. The placement and scoring procedures that were chosen were the triangle matcher and London dG. Out of a total of 100 poses for each docked molecule, the stiffer receptor was employed as the new refining strategy and the GBVI/WSA dG was employed as the new scoring methodology [46,47]. The optimal site for each ligand with the highest favorable scores, binding modes, and RMSD values was selected for further investigation. In the first step of the program validation method for the MOE program used, the co-crystallized instinctual inhibitor (N3) was redocked at its binding pocket of the generated main protease [48,49]. By obtaining a low root mean square deviation value (1.29) when comparing the freshly synthesized compounds and the redocked N3 ligand, a valid performance was demonstrated.

Molecular Dynamics Simulation (MDs)
MD simulations were performed using the Desmond 2020.1 (Schrödinger, New York, NY. 2017) from Schrödinger, LLC on the docked complex for 6LU7 with the 8a ligand. In this system, the explicit solvent model with the TIP3P water molecules and the OPLS-2005 force field [50][51][52] were applied in a period boundary salvation box with dimensions of 10 Å x 10 Å x 10 Å [53]. Na + ions were supplied to the system to balance the 0.15 M charge, and NaCl solutions were added to mimic physiological conditions. To retrain the system over the protein ligand complexes, the system was initially equilibrated using an NVT ensemble for 10 ns. Following the preceding phase, an NPT ensemble was used to complete a brief run of equilibration and minimization for 12 ns. The Nose-Hoover chain coupling approach [45,54] was used to set up the NPT ensemble and the variable temperature. Throughout all simulations, an active suspension of 1.0 ps and a pressure of 1 bar were maintained. The time step used was 2fs. The Martyna-Tuckerman-Klein chain coupling technique was used to manage pressure [55], employing a barostat method with a 2 ps relaxed period. The long-range electrostatic interactions were estimated using the particle mesh Ewald method [56], with the radius for the coulomb interactions set at 9. A RESPA integrator was used to determine the bonded forces for each trajectory during a time step of 2 fs. Making use of Geo Measures v0.8 ( https://github.com/lkagami/geo_ measures_pymol, accessed on 1 January 2023) [57], the complexes underwent primary component analysis (PCA). Geo Measures is provided with a substantial library of g sham and eigen values, which are represented in a 3D visual using the Python program matplotlib (https://github.com/matplotlib/matplotlib, accessed on 1 January 2023). The final production run lasted 100 ns. The root mean square deviation (RMSD), radius of gyration (Rg), root mean square fluctuation (RMSF), quantity of hydrogen (H-bonds), salt bridges, and SASA were calculated to monitor the stability of the MD simulations.

Binding Free Energy Analysis
The ligand-protein complex binding free energies were calculated using the molecular mechanics combined with the generalized Born surface area (MM-GBSA) method. Over the previous 50 frames, the Prime MM-GBSA binding free energy in the simulation trajectory with a one-step sampling size was calculated using the thermal mmgbsa.py Python script. The binding free energy of the Prime MM-GBSA (kcal/mol) was determined using the additivity concept, which required adding up each individual energy module, such as the columbic, covalent, hydrogen bond, van der Waals, self-contact, lipophilic, solvation of protein, and ligand modules.
The following equation is applied to determine Gbind: In which: -∆G bind specifies the binding free energy; -∆G MM specifies the difference between the free energies of the ligand-macromolecule complex and the total energies of receptor and ligand in isolated forms; -∆G Solv specifies the differences in the GSA solvation energies of the ligand-macromolecule complex and the sum of the solvation energies of the receptor and the ligand in the unbound state; -∆G SA specifies the difference in the surface area energies for the receptor and the ligand.
The thiazole derivatives 6a-e were produced through the reaction of compound 3 with the hydrazonoyl chlorides 4a-e [58] in the presence of Et 3 N. This was achieved by first performing a substitution reaction with the removal of the HCl molecule to produce the substituted intermediate 5, which was then followed by in situ cyclization with the removal of the water molecule (Scheme 1). Elemental analysis and spectral data ( 1 H-NMR, mass, IR) were used to clarify the structure of the products 6a-e. In each case, two stretching bands at 1692 and 3421-3160 cm −1 , attributed to the carbonyl and NH groups, could be seen in the IR spectra of product 6. The singlet signal at δ =10.69 ppm associated with the -NH proton was observed in the 1 H-NMR spectra of compound 6, in addition to the aromatic and alkyl protons. Each mass spectrum of products 6a-e showed a molecular ion peak with the appropriate molecular weight for that molecule. a molecular ion peak with the appropriate molecular weight for that molecule.
It was proposed that the hydrazine carbon of compound 4 is initially attacked by the thiol group of compound 3 to yield intermediate 5, which is then cyclized to products 6. By forming the Schiff bases 8a and 8d as a result of their interactions with 4-methoxybenzaldehyde 7 while being refluxed in acetic acid, the structural integrity of product 6 was further demonstrated. The structures of the isolated products 8a and 8d were elucidated based on their 1 H-NMR, IR, and mass spectra (see experimental section).

Scheme 1. Synthesis of thiazoles 6a-e and Schiff bases 8a,d.
A different synthetic approach might be used to create the real samples of 8a,d. Thus, Schiff base 9 was produced as a result of compound 3 reacting with 4-methoxybenzaldehyde 7 while being heated in EtOH\AcOH. After reacting with hydrazonoyl halides 4a,d in refluxing DMF, compound 9 produced the corresponding authentic samples 8a,d by forming intermediate 10 (Scheme 2).

Scheme 1. Synthesis of thiazoles 6a-e and Schiff bases 8a,d.
It was proposed that the hydrazine carbon of compound 4 is initially attacked by the thiol group of compound 3 to yield intermediate 5, which is then cyclized to products 6. By forming the Schiff bases 8a and 8d as a result of their interactions with 4-methoxybenzaldehyde 7 while being refluxed in acetic acid, the structural integrity of product 6 was further demonstrated. The structures of the isolated products 8a and 8d were elucidated based on their 1 H-NMR, IR, and mass spectra (see experimental section).
A different synthetic approach might be used to create the real samples of 8a,d. Thus, Schiff base 9 was produced as a result of compound 3 reacting with 4-methoxybenzaldehyde 7 while being heated in EtOH\AcOH. After reacting with hydrazonoyl halides 4a,d in refluxing DMF, compound 9 produced the corresponding authentic samples 8a,d by forming intermediate 10 (Scheme 2). Through its reaction with α-bromoketones, the viability of compound 3 as a building block for the synthesis of additional groups of the anticipated biologically active thiazoles was also investigated. Therefore, just one isolable product was produced in each case when compound 3 was allowed to react with the substituted phenacyl bromides 11a-c in DMF\TEA (Scheme 3). Based on the 1 H-NMR, mass and IR spectral data, the latter products were determined to be 13a-c. The 1 H-NMR spectra of product 13a, considered to be a sample example of the product 13, revealed two singlet signals assigned to the NH2 and thiazole-H5 protons, in addition to the predicted signals assigned to the aromatic protons. The mass spectra of products 13a-c demonstrated Scheme 2. Alternative synthesis of Schiff bases derivatives 8a,d.
Through its reaction with α-bromoketones, the viability of compound 3 as a building block for the synthesis of additional groups of the anticipated biologically active thiazoles was also investigated. Therefore, just one isolable product was produced in each case when compound 3 was allowed to react with the substituted phenacyl bromides 11a-c in DMF\TEA (Scheme 3). Based on the 1 H-NMR, mass and IR spectral data, the latter products were determined to be 13a-c. The 1 H-NMR spectra of product 13a, considered to be a sample example of the product 13, revealed two singlet signals assigned to the NH 2 and thiazole-H5 protons, in addition to the predicted signals assigned to the aromatic protons. The mass spectra of products 13a-c demonstrated peaks that matched their molecular ions.

Physiochemical and Pharmacokinetics Profiling
An ADMET investigation on all synthesized compounds was evaluated using the SwissADME platform (http://www.swissadme.ch/, last accessed on 21 December 2022). To ensure that the specified compounds were new drugs, the SwissADME tool was applied [59,60]. Table 1a-c contain data of the physicochemical parameters. The listed data in Table 1a-c, show that all synthesized compounds except 8d and 14b have a molecular weight < 500 [61], which increased the transport and absorption as well as improved the transmissibility of these thiazole compounds to membranes. With the exception of 6e and 13c, all of these thiazole derivatives' topological polar surface areas (TPSA) were less than 140A 0 , according to the reported range (Tables 1b,c) [62]. The lipophilicity (Log P) values of all the thiazole derivative examined lie within the 1.5-4.81 range, being acceptable as per Lipinski's five-factor rule. It is important to note that lipophilicity is associated to toxicity [63]. The findings indicate that the physicochemical properties of the newly synthesized thiazole derivatives are within the acceptable range, as shown by the bioavailability radar ( Figure 2).  Moreover, N-aminothiazoles 13a-c reacted with aldehyde 7 to yield the respective Schiff bases 14a-c (Scheme 3). The structure of compound 14 was deduced from its spectral data (IR, 1 H-NMR, and MS) and elemental analyses (cf. experimental part).
A different strategy was used to provide further support for the structure allocated to products 14a-c. As a result, the reaction between 9 and 11a-c in DMF under reflux produced a product that was identical to the product produced by the reaction between 13a-c and 7 in all aspects (Scheme 3).

Physiochemical and Pharmacokinetics Profiling
An ADMET investigation on all synthesized compounds was evaluated using the SwissADME platform (http://www.swissadme.ch/, last accessed on 21 December 2022). To ensure that the specified compounds were new drugs, the SwissADME tool was applied [59,60]. Table 1a-c contain data of the physicochemical parameters. The listed data in Table 1a-c, show that all synthesized compounds except 8d and 14b have a molecular weight < 500 [61], which increased the transport and absorption as well as improved the transmissibility of these thiazole compounds to membranes. With the exception of 6e and 13c, all of these thiazole derivatives' topological polar surface areas (TPSA) were less than 140A 0 , according to the reported range (Table 1b,c) [62]. The lipophilicity (Log P) values of all the thiazole derivative examined lie within the 1.5-4.81 range, being acceptable as per Lipinski's five-factor rule. It is important to note that lipophilicity is associated to toxicity [63]. The findings indicate that the physicochemical properties of the newly synthesized thiazole derivatives are within the acceptable range, as shown by the bioavailability radar ( Figure 2). A molecule's druglikeness is represented by the bioavailability radar. The pink area corresponds to the optimal range for each property (lipophilicity: Log P between 1.52 and 4.81; size: MW between 209 and 504 g/mol; polarity: TPSA between 92 and 142 0A2; solubility: log S no higher than 6; etc.), with the pink area representing the best range for each property. Carbons must make up at least 0.25;2.4 of the sp3 hybridization's carbon content to be considered saturated.

Molecular Docking Studies
Using a molecular docking simulation, the fifteen thiazole derivatives were tested for their capacity to engage with the main protease of COVID-19 (Pdb ID: 6LU7). Table 2 and Figures 3-7 list the outcomes of this docking investigation. First, a redock was performed on the co-crystallized ligand (N3) for verification. ASN 142, GLY 143, GLU 166, GLN 189, SER 144, and CYS 145 residues formed an H-bond acceptor with the C=O of the co-crystallized ligand, demonstrating a high docking score of -8 kcal/mol (Figure 3).  A molecule's druglikeness is represented by the bioavailability radar. The pink area corresponds to the optimal range for each property (lipophilicity: Log P between 1.52 and 4.81; size: MW between 209 and 504 g/mol; polarity: TPSA between 92 and 142 0A2; solubility: log S no higher than 6; etc.), with the pink area representing the best range for each property. Carbons must make up at least 0.25;2.4 of the sp3 hybridization's carbon content to be considered saturated.

Molecular Docking Studies
Using a molecular docking simulation, the fifteen thiazole derivatives were tested for their capacity to engage with the main protease of COVID-19 (Pdb ID: 6LU7). Table 2 and   thiazole derivatives 6b, 6e, 14a, 14b, and 14c were higher than those for the co-crystallized ligand. As is shown in Table 2 and Figure 4, these thiazole compounds were incorporated into the SARS CoV-2 main   docked compounds 6b, 6e, 14a, 14b, and 14c into the active site of M pro . Figure 4. Two-dimensional interactions of docked compounds 6b, 6e, 14a, 14b, and 14c

Simulation of Molecular Dynamics
The stability and convergence of the 8a-ligand for the COVID-19 major protease (Mpro) (PDB ID: 6LU7) were investigated using molecular dynamics modelling (MD) techniques. The simulation of 100 ns demonstrated a steady conformation when comparing the findings of the root mean square deviation (RMSD). When it was coupled with the 8a-ligand, the 6LU7 C-RMSD backbone produced a deviation of 1.9 A 0 . ( Figure  8A), while the ligand RMSD of 8a was depicted as 2.0 A 0 ( Figure 8A). Good convergence and stable conformations throughout the simulation are indicated by the stable RMSD graphs. The strong binding affinity of the ligand suggests that when 8a is linked to 6LU7, it is quite stable in its complex. Although the root mean square fluctuations (RMSF) plot for the 6LU7 protein showed small fluctuation spikes, no substantial spikes were found, suggesting that the residues may be more flexible. Most of the residues showed minimal fluctuation over the duration of the 100 ns simulation ( Figure 8B), which represents stable amino acid conformations. Since the protein structure is stiff during simulation in ligand-bound conformations, it can be deduced from RMSF plots.
The radius of gyration (Rg) measures a protein's degree of compactness. This experiment saw a decrease in the radius of gyration (Rg) of the 6LU7 C-backbone linked to the 8a-ligand from 22.3 to 22.01 ( Figure 8C). When the gyration (Rg) is noticeably reduced, this indicates that the protein is strongly oriented in a ligand-bound state. The presence of hydrogen bonds between the protein and the ligand points to the stability and strong interaction of the complex. Throughout the 100 ns of the simulation, there were considerable amounts of hydrogen bonds between compound 8a and 6LU7 ( Figure  8D). The average constant number of hydrogen bonds between 6LU7 and 8a-ligand was three on average ( Figure 8D). Salt bridges were formed between the oppositely charged residues close to each other and played a significant role in protein stability [64]. In this study, average single numbers of salt bridges were formed between 6LU7 and the 8a-ligand ( Figure 8E). An Rg analysis was followed by similar patterns being observed in the solvent accessible surface area (SASA), both in the ligand-bound and unbound states. It is evident from ( Figure 8F) that the protein 6LU7 had a high surface area, which was accessible to the solvent when the 8a-ligand was not attached to the receptor ( Figure 8F, red). When using the 8a-ligand to bind, the SASA value decreased in comparison to the unbound state ( Figure 8F, black). According to the overall analysis of the Rg, the All the thiazole compounds under investigation produced docking scores between −5.8 and −8.6 kcal/mol. The docking scores for the thiazole derivatives 6b, 6e, 14a, 14b, and 14c were higher than those for the co-crystallized ligand. As is shown in Table 2 and Figure 4, these thiazole compounds were incorporated into the SARS CoV-2 main protease (Pdb: 6LU7) active site through interactions with the amino acid residues, which were H-acceptor, H-donor, and hydrophobic. Compound 8a showed the highest binding affinity, demonstrating a binding energy value of −8.6 Kcal/mol against the active site of Mpro (Table 2). It showed three H-bonds with MET 17,GLN 19, and GLY 71, six hydrophobic interactions with GLU14, GLY 120, ALA 70, LYS 97, VAL 18, and TRP 31, as well as one electrostatic interaction with LYS 97 (Figures 4-6). In light of these encouraging docking simulation results using the SARS CoV-2 main protease (Mpro) (Pdb ID: 6LU7), we recommend that thiazole derivatives be tested in vitro as antivirals to suppress the SARS CoV-2 main protease.

Simulation of Molecular Dynamics
The stability and convergence of the 8a-ligand for the COVID-19 major protease (Mpro) (PDB ID: 6LU7) were investigated using molecular dynamics modelling (MD) techniques. The simulation of 100 ns demonstrated a steady conformation when comparing the findings of the root mean square deviation (RMSD). When it was coupled with the 8a-ligand, the 6LU7 C-RMSD backbone produced a deviation of 1.9 A 0 . (Figure 8A), while the ligand RMSD of 8a was depicted as 2.0 A 0 ( Figure 8A). Good convergence and stable conformations throughout the simulation are indicated by the stable RMSD graphs. The strong binding affinity of the ligand suggests that when 8a is linked to 6LU7, it is quite stable in its complex. Although the root mean square fluctuations (RMSF) plot for the 6LU7 protein showed small fluctuation spikes, no substantial spikes were found, suggesting that the residues may be more flexible. Most of the residues showed minimal fluctuation over the duration of the 100 ns simulation ( Figure 8B), which represents stable amino acid conformations. Since the protein structure is stiff during simulation in ligand-bound conformations, it can be deduced from RMSF plots. es Mol. Biol. 2023, 1, FOR PEER REVIEW 17 matching proteins were compelled to become more compact and less flexible when the ligands were bound. The radius of gyration (Rg) measures a protein's degree of compactness. This experiment saw a decrease in the radius of gyration (Rg) of the 6LU7 C-backbone linked to the 8a-ligand from 22.3 to 22.01 ( Figure 8C). When the gyration (Rg) is noticeably reduced, this indicates that the protein is strongly oriented in a ligand-bound state. The presence of hydrogen bonds between the protein and the ligand points to the stability and strong interaction of the complex. Throughout the 100 ns of the simulation, there were considerable amounts of hydrogen bonds between compound 8a and 6LU7 ( Figure 8D). The average constant number of hydrogen bonds between 6LU7 and 8a-ligand was three on average ( Figure 8D). Salt bridges were formed between the oppositely charged residues close to each other and played a significant role in protein stability [64]. In this study, average single numbers of salt bridges were formed between 6LU7 and the 8a-ligand ( Figure 8E). An Rg analysis was followed by similar patterns being observed in the solvent accessible surface area (SASA), both in the ligand-bound and unbound states. It is evident from ( Figure 8F) that the protein 6LU7 had a high surface area, which was accessible to the solvent when the 8a-ligand was not attached to the receptor ( Figure 8F, red). When using the 8a-ligand to bind, the SASA value decreased in comparison to the unbound state ( Figure 8F, black). According to the overall analysis of the Rg, the matching proteins were compelled to become more compact and less flexible when the ligands were bound.

Calculations of Molecular Mechanics Generalized Born Surface Area (MM-GBSA)
The binding free energy and additional contributing energy in the form of MM-GBSA were calculated for the 6LU7 + 8a complex using the MD simulation trajectory. This was followed by Rg analysis, which likewise showed a similar trend. The results (Table 3) suggested that the maximum contribution to ∆G bind in the stability of the simulated complexes were due to ∆G bind Coulomb, ∆G bind vdW, ∆G bind H bond , and ∆G bind Lipo, while ∆G bind Covalent and ∆G bind SolvGB were linked to the corresponding complexes' instability. The 6LU7 + 8a complex had comparatively higher binding free energies, higher than other complexes ( Table 3). The potential for 8a to bind to protein with a high affinity, efficiency, and the capacity to assemble a stable protein-ligand complex was substantiated by these findings.

Principal Component Analysis
The outcomes of a study to explain the random, global mobility of the atoms in amino acid residues are displayed in Figure 9's principal component analysis (PCA) of the MD simulation trajectories for 6LU7 + 8a. The more flexible scattered trajectories (0-600 frames) are interpreted by this technique as a result of non-correlated global motion due to the protein structure's randomness. A covariance matrix contained the internal coordinate mobility into three dimensions throughout the spatial time of 100 ns. Orthogonal sets, or eigenvectors, were used to represent the rational motion of each trajectory. The MD simulation trajectory of the Cα atoms of the 6LU7 + 8a protein displayed more unordered orientation in PC1 and PC2 modes and was oriented more toward a negative correlation from the initial 600 frames (Figure 9). Interestingly, for the last 400 frames (from 600-1000), it exhibited a positive correlation motion and clustered into a more oriented manner. As a result, it was obvious that the centering of the frames in a single cluster by 6LU7 + 8a (dark green) indicated that the periodic motion of MD trajectories was caused by steady, structural global motion. Consequently, the frames become more stable at the completion of the simulation (Figure 9).

Conclusions
According to the current work, a novel series of thiazole clubbed pyridines was created by reacting a pyridine thiocarbohydrazone derivative with a variety of hydrazonoyl halides and α-bromoacetophenones. The physicochemical parameters, toxicity assessment, and molecular docking approaches showed that compound (8a) was not toxic, did not violate Lipinski's rule of five, and could fit well with the binding site of the SARS-CoV-2 target protease through interactions with essential amino acids, hence inhibiting the replication function of the virus. Furthermore, the compound (8a) complex with the SARS-CoV-2 target protease showed stability during the 100 ns trajectory of the molecular dynamic simulation analysis. The predicted findings could be a helpful suggestion for us and other scientists to conduct further confirmation via in vitro experimental studies. Nevertheless, the preliminary results of the present study may pave the way for developing more potent agents against SARS-CoV-2 in the near future.
Supplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/xxx/s1, Physical and spectral data of the synthesized compounds and their 1 Hand 13 C-NMR spectra.

Conclusions
According to the current work, a novel series of thiazole clubbed pyridines was created by reacting a pyridine thiocarbohydrazone derivative with a variety of hydrazonoyl halides and α-bromoacetophenones. The physicochemical parameters, toxicity assessment, and molecular docking approaches showed that compound (8a) was not toxic, did not violate Lipinski's rule of five, and could fit well with the binding site of the SARS-CoV-2 target protease through interactions with essential amino acids, hence inhibiting the replication function of the virus. Furthermore, the compound (8a) complex with the SARS-CoV-2 target protease showed stability during the 100 ns trajectory of the molecular dynamic simulation analysis. The predicted findings could be a helpful suggestion for us and other scientists to conduct further confirmation via in vitro experimental studies. Nevertheless, the preliminary results of the present study may pave the way for developing more potent agents against SARS-CoV-2 in the near future.