The RNN was implemented using Tensorflow (v2.0²) and Keras (v2.3³) in Python (v3.7⁴) and RDkit⁵ was used for the processing of the molecules.

To accelerate the synthesis of potential inhibitors against 3CL^pro, a transfer learning approach was applied. Transfer learning is a machine learning technique where a pre-trained model is used as a starting point for training a new model with a similar task or domain. This allows the model to leverage the knowledge and experience gained from the pre-training to adapt to the new data and tasks more quickly and efficiently. In other words, transfer learning allows for faster and more accurate model development by building on top of previously learned representations.

Here, a publicly available model named LSTM_Chem (License: CC BY-NC-ND 4.0) was used (Gupta et al., 2017). The model consists of two LSTM layers with a 256-sized hidden state vector. It is regularized, having dropout layers. The two layers are followed by a final dense output layer with the softmax activation function. The model input is a bit array sequence of the molecule in the simplified molecular-input line-entry system (SMILES) format. This model was initially trained to produce novel TRPM8 inhibitors.

For the LSTM_Chem model to generate potential inhibitors for 3CL^pro, it was necessary to retrain and optimize the model on a ligand dataset that exhibits a certain degree of activity against 3CL^pro. The training process involved two stages with distinct datasets. The first stage-trained the model on a training dataset to learn the latent space features of chemical molecules. The second stage involved fine-tuning the model using a separate dataset to enable it to generate ligands that possess the chemical features of protease inhibitors for COVID-19.

The training dataset consists of a large volume of diverse ligands from which the RNN learns to produce valid ligands with high accuracy (Figure 1). The dataset was obtained from ChEMBL22⁶ and contained 556,134 SMILES strings, which were processed to remove duplicates, salts, and stereochemical information, resulting in a collection of unique ligands. Furthermore, only SMILES strings that had lengths between 34 and 74 tokens were retained, leading to a final size of 439,217 SMILES strings. This methodology was chosen following the work done by Gupta et al. (2017). The SMILES string length was constrained as having very long strings would result in the vanishing gradient problem, and the network would not learn anything. Although LSTMs are good at tackling the vanishing gradient, they are not completely immune to it (Moret et al., 2019). Additionally, the LSTM_Chem model accepts a bit array sequence as input, which was obtained by converting the SMILES strings using the Morgan algorithm in RDKit Open-source cheminformatics; (Open-Source Cheminformatics, see Footnote 5). The Morgan algorithm is a graph relaxation algorithm used for molecule canonicalization, which assigns a unique identifier to a molecule regardless of its representation. However, the Morgan algorithm has known issues that can result in noncanonical atom orderings and can be problematic when used with large molecules such as proteins. Therefore, restricting the length of the SMILES strings to fall within a range of 34–74 tokens limits the size of the molecules to small ligands and reduces the likelihood of encountering issues with the Morgan algorithm (Schneider et al., 2015; Moret et al., 2019; Schneider, 2019).

To create a dataset for fine-tuning the model, 845 drugs undergoing clinical trials and drugs that demonstrated activity in different biological assays for COVID-19 were collected from PubChem.⁷ This dataset included compounds that showed activity against not only 3CL^pro but other drug targets of SAR-CoV-2 as well. In PubChem BioAssay, “PUBCHEM_ACTIVITY_OUTCOME” is a column that reports the outcome of a specific assay run for a given compound. It describes whether the tested compound showed activity (i.e., produced a measurable effect) against the target of interest or not. Only compounds that were termed “active” were selected. The ligand structures were retrieved in SDF format using compound IDs, via PUG_REST, an application program interface (API) for accessing PubChem (Kim et al., 2018), and SMILES strings were generated using the MolToSmiles function from RDKit. This dataset was preprocessed similarly to the training dataset, resulting in 639 ligands that were used for fine-tuning the model. However, unlike the training dataset, no restriction was set on the string length of the SMILES in the fine-tuning dataset.

The model underwent training for 50 epochs on the training dataset, followed by 25 epochs of fine-tuning on the fine-tuning dataset. During fine-tuning, 1,000 SMILES strings were generated and evaluated for their validity. The validated ligands were then subjected to docking onto 3CL^pro. The ligands that exhibited a binding affinity greater than that of the native ligands were incorporated back into the fine-tuning dataset. The fine-tuning dataset underwent preprocessing and shuffling, with this entire fine-tuning process being repeated for three cycles. This process of adding validated ligands with high-binding affinity back into the dataset and repeating the fine-tuning process is a way to iteratively improve the performance of the model and helps it identify the characteristics of ligands that contribute to their strong binding affinity with 3CL^pro. A schematic overview of the model training and optimization process is provided in Figure 1.

Virtual screening involves docking ligand libraries to a target macromolecule to discover a lead that would confer a biological function. The virtual screening was done using AutoDock Vina in PyRx (Dallakyan and Olson, 2014).

The generated ligands were converted from SMILES to SDF using openBabel-GUI (O’Boyle et al., 2011). To obtain the lowest free energy of the ligand, the Merck molecular force field (mmff94) parameter was used in PyRx. Finally, the ligands were converted to PDBQT format, preparing the ligand for molecular docking (Morris et al., 1998; Huey et al., 2007).

The 3D crystalline structure of SARS-CoV-2 main protease or 3CL^pro (PDB: 6 LU7) was obtained from PDB⁸; this served as the target for docking. The target was prepared by removing the native ligand present (N3 Inhibitor) and water molecules using Biovia Discovery Studio (Biovia, 2017).

The native ligand was docked onto the target molecule, and the binding affinity was found to be −7.9 kcal/mol. The amino acid residues involved in binding with the native ligand were obtained using the 2D structure in Biovia Discovery Studio. The amino acid residues are Thr24, Thr26, Phe140, Asn142, Gly143, Cys145, His163, His164, Glu166, His172. The grid box was then positioned over the binding site (center: x = −10.606, y = 17.214, z = 64.716, total size: x = 24.074 Å, y = 24.134Å, z = 19.174Å).

Further post-docking analysis and visualization of the ligand-target complex were carried out in Biovia Discovery Studio (Biovia, 2017).

To evaluate the molecular properties of the ligands, an online tool, Molinspiration⁹ was used by uploading the ligands in SMILES format. Molinspiration also provides bioactivity scores for drug targets such as Ion channel modulators, GPCR (G protein-coupled receptor) ligands, kinase inhibitors, nuclear receptor inhibitors, protease inhibitors, and enzyme inhibitors. The bioactivity score of Molinspiration is calculated by a machine learning-based model that predicts the probability of a molecule being active against a particular target. The model is trained on a large database of known active and inactive compounds and uses various molecular descriptors, such as physicochemical properties and substructure information, to make predictions. These bioactivity scores provide an additional metric for evaluating the drug-like properties of the ligands (Vardhan and Sahoo, 2020). Another online tool Molsoft¹⁰ was used to evaluate the drug-likeness score of the ligands (Prabhavathi et al., 2020).

The absorption, distribution, metabolism, elimination, and toxicity (ADMET) are some of the important pharmacokinetic properties that must be evaluated. An online tool called admetSAR¹¹ was used to obtain the ADMET profile of the ligands (Yang et al., 2018). Some of the properties calculated include ames mutagenesis, blood-brain barrier penetration, BSEP inhibition, Caco-2, Carcinogenicity, cytochrome p450 substrate and inhibitors, glucocorticoid receptor binding, hepatotoxicity, human ether-a-go-go (hERG) inhibition, p-glycoprotein inhibitors and substrate, human intestinal absorption, and human oral bioavailability.

To evaluate the capability of the generated ligands as potential anti-COVID drugs, a reference dataset comprising 20 drugs currently in clinical trials for COVID-19 treatment was obtained from PubChem. This included drugs such as Remdesivir, Ritonavir, Galidesivir, etc. (For full list of drugs—Appendix A). The properties of the generated ligands were compared to those of the reference drugs to assess their potential as anti-COVID agents. This comprehensive comparison of the ligands and clinical trial drugs facilitates the assessment of the ligand’s properties.

Molecular Dynamics simulation is a sophisticated computational tool for predicting and analyzing the dynamic behavior of molecules (Verdonk et al., 2003; Radinnurafiqah et al., 2016; Girdhar et al., 2019; Choubey et al., 2022; Mishra et al., 2022). The stabilities of six selected protein-ligand complexes were assessed using GROMACS 2021 package through Molecular Dynamics (MD) simulations (Van Der Spoel et al., 2005). The complexes included SARS_COV2_MOL_1, SARS_COV2_MOL_3, SARS_COV2_MOL_9, SARS_COV2_MOL_10, SARS_COV2_MOL_17, and SARS_COV2_MOL_20. The ligand topology parameter for CHARMM forcefield (Vanommeslaeghe et al., 2010) was created using the CGenFF server.¹² A cubic box of TIP3P water models was used to solvate all the complexes. To maintain the periodic boundary conditions, the distance between the protein and the box edge was kept at 1 nm. The systems were neutralized by adding 0.15 M NaCl. Energy minimization was performed using the steepest descent method followed by the conjugate gradient method with maximum number of minimization 50,000 per alogorithm. The Particle Mesh Ewald (PME) method was employed to calculate long-range interactions (Abraham and Gready, 2011). The first phase of equilibration was carried out with an NVT ensemble, where the temperature was equilibrated using 50,000 iterations of 2 fs each. In the second phase, the pressure was equilibrated at 300 K with an NPT ensemble using Parrinello-Rahman, a pressure coupling method. The temperature inside the system was regulated using V-rescale, a modified Berendsen thermostat. Finally, a production run of 100 ns was established to gain insights into the dynamic behavior of the complex.

The obtained trajectories after the MD simulations were analyzed for calculations such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and inter-molecular hydrogen bond using the in-built tools of the GROMACS package. To compute the RMSD in the protein backbone, the rms module of GROMACS was employed. RMSD of the ligands were also calculated using the same module, whereas the rmsf module was used to determine the RMSF in the atomic positions of the protein Cα backbone. In addition, modules like h-bond, gyrate, and SASA were used to calculate the number of hydrogen bonds, Rg, and SASA, respectively.

After the first stage of training on the training dataset, a final loss of 0.427 on the training set and 0.567 on the validation set (20% of data from the training dataset) was obtained. Additionally, the model had an accuracy of 82% in generating valid ligands, i.e., out of every 100 ligands the model produces, 82 are valid molecules.

The model then underwent three cycles of fine-tuning on the fine-tuning dataset, and after each cycle, the binding affinity of 30 randomly selected ligands was evaluated. Figure 2 depicts the distribution of the binding affinities across the three cycles, indicating that the 3rd generation of molecules had a significantly better binding affinity to 3CL^pro than the previous two generations (Mann Whitney U Test, p < 0.001). The average binding affinity for the 3rd generation was −8.406 ± 0.087 kcal/mol, and 26 (86.67%) of the ligands had a value higher than the binding affinity of the native Ligand (−7.9 kcal/mol).

After evaluating the binding affinity of 30 ligands generated by the model, the top 20 ligands were chosen for additional investigations. To confirm the novelty of these molecules, a search was conducted in the PubChem database, which did not yield any results for these ligands. Therefore, it can be inferred that the generated molecules are novel. The molecules were named SARS_COV2_MOL_1 - SARS_COV2_MOL_20 as an identifier.

Various molecular and ADMET properties of the generated ligands and drugs in the reference dataset were calculated and contrasted to assess the efficacy of the generated ligand to serve as a potential inhibitor to the 3CL^pro protease (Tables 1, 2).

The drug-likeness score of the generated ligands predicted by MolSoft are shown in Figure 8A. The average drug-likeness score was found to be 0.663 ± 0.118, and the scores of the generated ligands were found to be in the range with FDA reference drugs used by MolSoft. Drugs that have drug–likeness scores greater than 0 are considered to have drug-like properties (Prabhavathi et al., 2020). In Figure 8B, ligands that fall in the green region have a binding affinity greater than −7.9 kcal/mol and have a drug-likeness score greater than 0. 16 (80%) ligands fall under the green region, indicating that these ligands have a good potential to be developed as lead molecules.

From the 20 generated ligands, 6 were selected to study their interactions with 3CL^pro using Discovery Studio. Figures 9A–F, displays the interaction of the six ligands with the protease, and the amino acid residues interacting with the ligand have been annotated. Hydrogen bonds are an essential factor that determines the stability of the docked complex. Most of the ligands can be seen interacting with receptor residues Leu 141, Asn 142, Gly 143, and Ser 144 via a hydrogen bond. SARS_COV2_MOL_17 (Figure 9E) shows the highest binding affinity of −9.4 kcal/mol and a drug-likeness score of 0.29. The ligand showed hydrogen bond interaction with receptor residues at Leu 141, Gly 143, and Ser 144. However, there is an unfavorable donor-donor interaction at Cys 145. Unfavorable bonds greatly hinder the stability of the protein-ligand complex. SARS_COV_MOL_20 (Figure 9F) has the highest drug-likeness score of 1.57 and a binding affinity of −8.3 kcal/mol. This ligand can be seen interacting with hydrogen bonds at receptor residues Thr 26, Glu 166, and Gln 189. However, among the six ligands, SARS_COV2_MOL_1, 3, 9, and 10 (Figures 9A–D respectively) have shown high binding affinity with a good drug-likeness score and less than 1 violation of Lipinski’s Rule. SARS_COV2_MOL_9 has shown a high drug-likeness score of 1.27 and is also a good protease inhibitor (Protease inhibitor bioactivity score - 0.26 (refer Appendix-A)). It also exhibits a high number of hydrogen bond interactions at the receptor (Thr 24, Thr 26, Gly 143, Glu 166, and Gln 189) and interacts with a pi-pi stacked bond at His 41 and pi-alkyl bond at Pro 168, and Met 49. In conclusion, the generated ligands are seen to be forming stable complexes with 3CL^pro protease having high binding affinities.

To further validate the structural stability of the generated ligands, MD simulations was performed on the docked complexes for 100 ns. The trajectories obtained after the simulations were analyzed to calculate RMSD, RMSF, hydrogen bonds, Rg, and SASA to assess the stability of the simulated systems.

RMSD is commonly used to evaluate docked complex stability (Martínez, 2015; Sargsyan et al., 2017). It measures the difference between the initial position and the final conformation of the protein backbone. From Figure 10, the SARS_COV2_MOL_3 (red) complex showed the lowest average RMSD value, around 0.18 nm, among all six complexes. The average RMSD values of the other five complexes with SARS_COV2_MOL_1 (black), SARS_COV2_MOL_9 (green), SARS_COV2_MOL_10 (blue), SARS_COV2_MOL_17 (yellow), and SARS_COV2_MOL_20 (brown) was estimated to be ~0.2 nm same as that of the apoprotein. The RMSD values remained nearly constant over the 100 ns period, indicating that the protein-ligand complexes were structurally stable (Figure 9A). The low RMSD values suggest that the protein-ligand interactions were energetically favorable and contributed to the stability of the complexes. The SARS_COV2_MOL_20 and SARS_COV2_MOL_1 complexes showed slight deviations of ~0.4 and ~0.38 nm around 50 and 66 ns, respectively. These fluctuations were further confirmed by observing local changes at the residue level using the RMSF plot.

The ligand RMSD ranged between 0.16–0.39 nm as shown in Figure 9B. SARS_COV2_MOL_17 had the least average RMSD which suggest the stability of the protein-ligand system when bound with SARS_COV2_MOL_17. Rests of the ligands also showed no sharp deviation and were stable throughout the period of simulation. It indicates the stability of the complexes. Alignment of the post-MD complexes with their respective initial-docking poses corroborate the low deviations observed, as illustrated in Figure 11.

The fluctuations in the protein can be determined by calculating RMSF, which measures the flexibility of each residue over time. The stability of the protein-ligand complexes can be inferred from the RMSF scores, with higher values indicating less stability and more flexibility. The RMSF of Cα atoms was calculated for all complexes, and the resulting average values for the six ligands were between 0.103 and 0.150 nm, as shown in Figure 12. The RMSF of the apoprotein was around 0.134 nm with no major fluctuation with respect to the protein-ligand complexes. However, the residues in the range of 46–50 showed slight fluctuation of average 0.37 nm. These low RMSF values suggest that the protein-ligand complexes are relatively stable and exhibit a moderate degree of flexibility. The close proximity of the complexes with the apoprotein also indicates about the stability of the complexes. This indicates that the ligands can bind to the protein without causing significant changes in its conformation. So, the predicted system appears to be stable (Farmer et al., 2017).

The compactness and stability of protein structures can be measured using the Rg, which represents the mass-weighted root mean square distance of the atomic distribution from their mutual center of mass. The Rg values depict the inclusive dimensions of the protein and protein-ligand complexes and reflect their appropriate interactions. The protein-ligand complexes that displayed the least radius of gyration are considered to be more compact and stable. The Rg values of all the complexes were analyzed, and the solvent-accessible surface area (SASA) was also computed for all the proteins for 100 ns. The SASA is an important measure to determine the area of the receptor exposed to the solvents during the simulation. As shown in Figure 13A, all the systems exhibited similar Rg values, ranging from 2.2 to 2.6 nm throughout the simulation, indicating their stability. The estimated SASA values also displayed a similar pattern, varying between 150.59 and 153.93 nm², as depicted in Figure 13B, with the highest value observed for the SARS_COV2_MOL_3 complex. These observations confirm the stability and compactness of all the protein-ligand complexes, as smaller deviations in average Rg and SASA values (Figure 13) Suggesting a stronger binding between the protein and ligand (Shaji, 2016).

The stability and molecular recognition process of a protein-ligand complex is affected by the intermolecular hydrogen bonds (H-bonds) between interacting atom pairs. The number of H-bonds formed between the receptor protein and selected ligands was determined during the 100 ns MD simulations to ascertain the dynamic stability of each complex. The binding strength and specificity of the protein-ligand complex are determined by hydrogen bonds. Figure 14, represents the number of hydrogen bonds formed between the receptor protein and selected ligands throughout the MD simulation. The complex formed with SARS_COV2_MOL_3 (red) showed a higher number of hydrogen bonds, while the rest of the complexes showed a stable number of hydrogen bonds throughout the 100 ns simulation (Figure 12). The results further suggest the stability of the studied 3CL^pro inhibitors (Pereira et al., 2019; Zhu et al., 2022).