The two DEL libraries are combined as a dataset containing 1.04 billion molecules with corresponding count and enrichment fold (EF) values. We provide the information of DEL in the supplement material. We cut out the data with very low count, the remaining data is 3,702,672. Then it was divided into a training set and a test set (0.8/0.2), of which the training set and test set have 2,962,138 and 740,534 molecules, respectively. First, the training set is used to train the molecule generation model. Then, molecules 18,129) with higher count values in DEL were selected as the positive samples of transfer learning to fine-tune the pre-trained model.

SMILES or molecular graphs are commonly used for molecular representation. SMILES is a text of molecular structures (Weininger, 1988). Molecular graphs represent the structure of molecules as graphs, where the edges of the graph represent bonds and the dots represent atomic structures (Sun et al., 2020). Molecules were represented as SMILES in three models, RNN, VAE, and AAE, to perform molecule generation. We used pytorch, sklearn, numpy, rdkit,umap-learn, and MOSES benchmark platform (Polykovskiy et al., 2020) to complete the experiments.

Recurrent neural network (RNN) (Segler et al., 2018) are designed to learn sequential data such as text or speech. The SMILES in DEL are made into a corpus. RNN can learn the grammatical information about the SMILES to know which parts of molecules tend to be connected. RNN can generate sequences through forward propagation (Bjerrum and Threlfall, 2017). By treating the molecule generation process as a series of steps and sampling the network at each step, generating effective molecules is highly probable and structurally similar to the training molecules. The architecture used in RNN consists of an embedding layer, three LSTM layers, and a linear layer.

Variational autoencoder (VAE) (Bowman et al., 2015) consists of an encoder and a decoder. The former encodes the input data into a latent vector, which obeys the Gaussian Distribution. The decoder restores the latent vector result to the target sample. SMILES are used as the model’s input and output to establish a VAE model. The VAE architecture consists of an embedding layer, an encoder layer, and a decoder layer. The encoder and decoder layers consist of a GRU layer and two linear layers.

Adversarial Autoencoder (AAE) (Makhzani et al., 2015) is similar to VAE principally. The difference is that based on the encoder and the decoder, a discriminator is introduced, which is responsible for distinguishing the calculated latent vector in the encoder from the real sample. The encoder and decoder are still accountable for encoding and reconstructing the data. AAE also uses SMILES as input and output. The AAE’s encoder part includes an embedding layer, an LSTM layer, and a linear layer. The decoder consists of two linear layers, one embedding layer, and one LSTM layer. The discriminator consists of two linear layers where the activation function is ELU (Clevert et al., 2015).

Transfer learning (Amabilino et al., 2020) is a fine-tuning model technique that fixes the original model’s specific parameters while others are still training and updating. This technique aims to streamline expansive chemical space in the generative model, searching for target molecules in the relatively small chemical space. The model is retrained by inputting molecules with high-count values to generate more distribution-similar molecules than those with high-count values. To AAE, we fine-tune the decoder’s last linear layer and the discriminator’s last linear layer. To VAE, the last two linear layers of the decoder’s model are fine-tuned. To RNN, we also fine-tune the last linear layer.

Each model generated 10,000 molecules, which were evaluated using the evaluation metrics provided in Moses (Polykovskiy et al., 2020), including valid, unique, novelty, internal diversity (IntDiv), and scaffold similarity (Scaff).

Validity is the proportion of valid molecules in the generated molecules.

The number of molecules generated by each model is 1,000,000 molecular datasets. First, the dataset was filtered for validity and reproducibility, followed by drug-likeness: 250 ≤ MW ≤ 750, logP ≤ 5, HBD ≤ 5, HBA ≤ 10, RB < 10, and 0.5 < QED. The next step was improving drugability by applying Medicinal Chemistry filters (MCFs) (Kalgutkar et al., 2005) and Pan Assay Interference Compounds (PAINS) filtering (Baell and Holloway, 2010). Without transfer learning, RNN, AAE, and VAE have no significant performance differences and predicted affinity distribution. Therefore, these datasets are merged and divided into groups depending on whether transfer learning is applied or not. 5000 molecules were obtained after pharmacophore filtration. Finally, molecular docking was employed, and molecule with higher docking score was selected for molecular dynamics simulation.

The remained SMILES after validity and repeatability filtration and the original 3,702,672 SMILES from the DEL dataset were transformed into Morgan fingerprints with 1024 dimensions and 2 radius (Rogers and Hahn, 2010). These fingerprints were then used to build a UMAP(Uniform manifold approximation and projection) (metric = “jaccard,” n_components = 2) model for dimensionality reduction visualization (McInnes et al., 2018).

The model for affinity prediction was established according to our previous study (Xiong et al., 2022). First, we sorted the molecules in the DEL dataset by the count value, then oversampled the top 10,000 ranked molecules by ten times. The step was set as 0, and every other step of the remaining molecules was sampled to form a training set.

The 3D molecular similarity was calculated through the shape and color similarity score (SC score), which represents the pharmacophoric feature similarity (Landrum et al., 2006) and the shape similarity (Putta et al., 2005). This score was used for the previously generated dataset. The 3D similarity score is a floating point value in the range of [0, 1], with a higher value indicating higher similarity between candidate and reference molecules. The native ligand in PDB:7L13 from the RSC-PDB database was used as a reference structure (Zhang et al., 2021). 100 conformations were generated for each molecule from the dataset using the RDKit UFF (Universal Force Field) force field. The lowest energy conformation was applied for the next step.

The A-chain of the complex PDB:7L13 (resolution 2.17 Å) of 3CLpro protein was split as a docking template to obtain accurate docking results. Subsequently, the complex was preprocessed using the Protein Preparation Wizard module of the Maestro suite (version: 13.1.141, Schrödinger Inc.) with the default setting, including the addition of hydrogen and side chains, removal of water molecules, and calculation of partial charges and protonation states using the OPLS4 force field (Poltev et al., 1996). Followed by a grid generation module, a similar-sized grid box centered on the native ligand was made to determine the binding pocket. All molecules were preprocessed by the LigPrep module. The ionization states were calculated using Epik (Shelley et al., 2007) at pH = 7.0 ± 2.0. Finally, all molecules were docked into the binding pocket within the grid and evaluated using the standard precision (SP) of Glide-v9.4. The scale factor and partial charge intercept are set to 0.8 and 0.15, respectively. 1000 poses per ligand were generated for docking evaluation. Post-docking binding site analysis and generation of interaction graphs were finished using Maestro.

Amolecular dynamics simulation was carried out to analyze further the dynamic interaction process between protein and ligand and the stability of binding status. Molecular dynamics simulation is a popular technique to study protein motion by tracking its conformational changes over time (Collier et al., 2020). Molecular interaction and visualization analysis based on SP docking results, the top-ranked molecules were used for the molecular dynamics simulation (MD-simulation) study. MD-simulation was performed using the GROMACS software package (version 2021.5) (Rakhshani et al., 2019). The AMBER14SB force field parameter was used for the protein. The ligand atomic charge was calculated using the B3LYP/6-31G* basis set. The ligand topology was computed using the GAFF2 force field parameter. The TIP3P water model was used to add Na⁺ and Cl⁻ ions to neutralize the charge. Electrostatic interactions are handled separately using the Particle Mesh Ewald (PME) and Verlet algorithms. The heavy atoms of the protein are constrained, and the energy minimization is carried out through 50,000 steps using the steepest descent method. The simulated system was equilibrated for 100 ps using a canonical ensemble (NVT) and an isothermal-isobaric ensemble (NPT). Both van der Waals and Coulomb interactions were calculated using a cutoff of 1.4 nm. Afterward, the system was run at constant temperature (300 K) and constant pressure (1 bar) for 100 ns molecular dynamics simulations with a time step of 2 fs and trajectory data saved every 5 ps. Finally, the ligand and protein complex’s root mean square deviation (RMSD, Å) at 100 ns was measured. By examining the interaction of the ligand with active site residues and the structural changes of the complex, the complexes were considered stable.

We evaluate each model’s validity, uniqueness, novelty, intDiv, and Scaff metrics (Tables 1, 2). The validity and uniqueness indicators of all models perform satisfactorily, indicating that the models can learn the grammatical information of the SMILES structure. The performance of novelty and IntDiv indicators is relatively poor, meaning that the model’s generalization ability may not be strong enough. VAE and RNN have higher Scaff values, meaning that the model can generate the same skeleton as the training set, but the ability of generating new skeleton is weak, while AAE is the opposite. In other words, the generated dataset and the training dataset had an apparent overlap. Especially after using transfer learning, Novelty’s metrics dropped further. A more complex model may be beneficial to address such a problem, so the MCMG (Multi-constraint molecular generation) model was also established (Wang et al., 2021). Unfortunately, MCMG performed relatively poorly in affinity prediction, so we decided not to analyze it further (Supplementary Figure S13).

The results visualization using UMAP dimensionality reduction are shown in Figure 2. The molecules generated by each model closely resemble the chemical space distribution of DEL’s. This indicates that the model could learn molecular distribution sufficiently from the source dataset.

Figure 3 shows the distribution of the molecules from each model’s affinity predictions. The molecules distribution without transfer learning is mainly located in the area of 0.1–0.2, while the corresponding values after using transfer learning are mostly above 0.2. Such improved affinity indicates the beneficial effect of transfer learning, which expectedly to improve the success rate and efficiency for further structural optimization. Scheme 2 shows the representative molecules with high affinity scores, which were expected to be potentially bio-active.

According to the calculated SC scores of all molecules in the dataset and native ligand, 5,000 ligands with the highest SC score were selected for the follow-up study.

Molecular docking was utilized to analyze the 3D conformational and pharmacophore-based screens and to study the structural basis of the interaction between 3CLPro and ligands. First, the reliability of the glide docking algorithm (standard precision mode, enhanced conformational sampling by four times) was confirmed by re-docking the native ligand to the receptor.

The re-docked conformation was presented in supporting information. Subsequently, the selected 5000 molecules were preliminarily docked to the revealed binding site using standard precision mode. According to the docking evaluation score and molecular conformation, 500 different conformations were selected, and four times enhanced conformational sampling was used to generate the ligand-binding pose more accurately.

In addition, the docking scores with and without transfer learning for pharmacophoric models were analyzed (Table 3). Applying transfer learning, 4.9% of ligands possess a score greater than 8, while 0.3% have a score higher than 9. In contrast, the corresponding values without transfer learning are 3.3% and 0.1%. This result is consistent with the affinity model, indicating that transfer learning can effectively increase the percentage of positive samples.

Moreover, the molecules filtered by the affinity model also performed molecular docking as a pharmacophoric model. This aims to explore the possibility of replacing an external pharmacophore model with an own-built affinity model. The docking scores of the filtered molecules were combined in Table 3 for comparison. We found an exciting revelation that the proportion of molecules filtered by the affinity model with docking scores of 9–11 was twice as high as that from the pharmacophore model. Therefore, replacing pharmacophore filtering with affinity models alone may be a promising option to utilize DEL’s datasets more fully.

The binding free energy can be used as a reference standard for assessing the activity of molecules. Generally, the lower the binding value means, the more stable the complex formed is. The ligands’ binding free energies were calculated by psp-v6.7 MMGBSA. The self-established ADMET prediction model and Maestro’s QikProp module evaluated the corresponding properties of molecules with better conformation. In fact, our ADMET prediction model is similar to ADMETlab (Dong et al., 2018). Molecules N1-N8 were finally selected for follow-up research considering the above ADMET properties with results of binding energy calculation (Scheme 3).

From Table 4, the binding free energies (dG_Bind) of N1 to N8 indicate their potential biological activities. Van der Waals energy (dG_Bind_vdW) shows that hydrophobic interaction is the main contributor to the ligand binding process. According to the molecular docking conformation, we selected N1 and N2 for subsequent research, in which the position and interaction of N1 and N2 on receptors are consistent with the previous report (Figure 4) (Zhang et al., 2021; Stille et al., 2022).

To further analyze the interaction-related atomic details between molecules and 3CLPro, GROMACS was applied for molecular dynamics simulations using the above docking results. 100 ns run time for MD simulation is considered sufficient for side chain rearrangement. The result will confirm whether or not the complex remains in the most stable association. Since N2 is believed to be more compatible with the receptor pocket (ligand conformation), and the calculation of MMGBSA binding free energy indicates that N2 has a stronger interaction with the receptor, molecular dynamics simulations of the complex formed by N2 were conducted.

RMSD values plotted over the simulation time revealed a stable kinetic equilibrium of the complex. In detail, the 3CLPro protein with N2 and the ligand N2 showed steady kinetics after 30 and 50 ns, respectively (Figure 5). By monitoring the fluctuation of RMSD, each system is in the range of 2Å after 50 ns. This suggests that the complex undergoes a conformational change during the simulation that promotes tight binding between the N2 and receptor, and finally, the system reaches a steady state. The RMSF value showed minimal fluctuation, and it remained in the range of 0.05–0.2 nm throughout the simulation period for most residues, except that a peak in RMSF value was observed only at residue 1. The less fluctuating performance confirms the strong attachment of the ligand to the protein (Figure 6).