*Article* **In Silico Design, Synthesis and Biological Evaluation of Anticancer Arylsulfonamide Endowed with Anti-Telomerase Activity**

**Giulia Culletta 1,2 , Mario Allegra 2 , Anna Maria Almerico 2 , Ignazio Restivo <sup>2</sup> and Marco Tutone 2, \***


**Abstract:** Telomerase, a reverse transcriptase enzyme involved in DNA synthesis, has a tangible role in tumor progression. Several studies have evidenced telomerase as a promising target for developing cancer therapeutics. The main reason is due to the overexpression of telomerase in cancer cells (85–90%) compared with normal cells where it is almost unexpressed. In this paper, we used a structure-based approach to design potential inhibitors of the telomerase active site. The MYSHAPE (Molecular dYnamics SHared PharmacophorE) approach and docking were used to screen an *in-house* library of 126 arylsulfonamide derivatives. Promising compounds were synthesized using classical and green methods. Compound 2C revealed an interesting IC<sup>50</sup> (33 ± 4 µM) against the K-562 cell line compared with the known telomerase inhibitor BIBR1532 IC<sup>50</sup> (208 ± 11 µM) with an SI ~10 compared to the BALB/3-T3 cell line. A 100 ns MD simulation of 2C in the telomerase active site evidenced Phe494 as the key residue as well as in BIBR1532. Each moiety of compound 2C was involved in key interactions with some residues of the active site: Arg557, Ile550, and Gly553. Compound 2C, as an arylsulfonamide derivative, is an interesting hit compound that deserves further investigation in terms of optimization of its structure to obtain more active telomerase inhibitors

**Keywords:** sulfonamides; arylsulfonamide; anticancer compounds; telomerase inhibitors; structurebased drug design; pharmacophore modeling; docking; molecular dynamics

#### **1. Introduction**

The nuclear protein complex, Telomere, defends the terminal ends of chromosomes from degradation, end-to-end fusion, and recombination [1–3]. The telomere structure is subject to several changes. After each cell division cycle, the telomere gradually shortens until the chromosomal DNA is exposed, inducing a DNA damage response [4,5]. This event helps to maintain the stability of genetic information and protects the genome in a "timebomb" manner [6]. When the length of telomeres reaches a critical point, cells reach the cycle of termination, aging, and death [5,6]. Normal cells cannot survive this progressive shortening. Sometimes, cells can extend telomeres by reactivating telomerase activity or through a telomere replacement elongation mechanism (ALT) to help cells survive the crisis [7]. The reactivation of telomerase is observed in 85–90% of human tumor cells [8]. Telomerase is a reverse transcriptase that contains an RNA template (TER) with its binding domain (TRBD) and reverse transcriptase unit (TERT). Telomerase is a challenging but appealing target, because the inhibition of its activity can be achieved by acting at different stages and using various mechanisms [9–11]. A first but now less popular used approach is based on inhibiting telomerase access to DNA by stabilizing G-quadruplexes formed by a 3' DNA overhang, which has the task to block telomere elongation [11–15]. Another

**Citation:** Culletta, G.; Allegra, M.; Almerico, A.M.; Restivo, I.; Tutone, M. In Silico Design, Synthesis and Biological Evaluation of Anticancer Arylsulfonamide Endowed with Anti-Telomerase Activity. *Pharmaceuticals* **2022**, *15*, 82. https:// doi.org/10.3390/ph15010082

Academic Editor: Osvaldo Andrade Santos-Filho

Received: 6 December 2021 Accepted: 5 January 2022 Published: 10 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

approach includes the use of inhibitors of the enzymatic active site of telomerase including antisense compounds to TER [16,17]. But the most promising approach is possibly the use of compounds that are able to block the active center of the enzyme in the catalytic subunit, directly. BIBR1532, which selectively inhibits telomerase activity, was first proposed in 2001 and has been patented as a potential anticancer drug [18,19]. As reported in a recent review, structure-based drug design strategies can be used to design potential inhibitors of the telomerase active site [20]. The first structure showing a co-crystallized inhibitor was published in 2015, revealing the *Tribolium castaneum* full-length catalytic subunit of telomerase in complex with the compound BIBR1532 (PDB 5CQG), which showed that the studied ligand interacts outside of the active center [21]. This experimental crystal structure changed the understanding of the expected mechanism of activity of this compound and possibly other active site-directed molecules. BIBR1532 binds to a highly conserved hydrophobic pocket (FVYL) motif on the outer surface of the thumb domain of telomerase.

Several other compounds designed using in silico approaches have been tested as potential inhibitors of the catalytic subunit of telomerase such as benzylidene–hydrazone analogs [22], dihydropyrazole [23–29], dibenzopyrroles [30], flavone pyridines [31], 1,3,4-oxadiazole derivatives [32–36], pyrazole-5-carboxamides and pyrazole-pyrimidines [37], spiroketals [37], celastrol derivatives [38], myricetin derivatives [39], indolyl-2'-deoxynucleotide analogs [40], flavonoid derivatives [41], and chrolactomycin derivatives [42]. Other classes of compounds developed as telomerase inhibitors are reported in two recent reviews [43,44]. Thus, considering the direct inhibition of TERT as the most promising approach for blocking the action of telomerase and the lack of an approved drug (Imetelstat is the only one to reach clinical trials [20])makes the search for new and effective inhibitors a hot topic of the pharmaceutical sciences. To the best of our knowledge, there is just one piece of evidence of arylsulfonamide derivatives (SEW05920) as an inhibitor of the TERT [45]. Exploiting our previous experience and outcomes in the use of computational approaches [46–50], we developed a structure-based computational approach to performing a virtual screening of an *in-house* arylsulfonamides library. Some arylsulfonamides have been identified as hits, also synthesized using green chemistry approaches and tested against three cancer cell lines K-562, HCT-116, and MCF-7. The results determined compound 2C as more active than the reference compound BIBR1532. The identification of this novel scaffold, which will undergo optimization, could help to identify new potent telomerase inhibitors.

#### **2. Results**

#### *2.1. In Silico Modeling and Virtual Screening*

In this study, we performed a structure-based virtual screening by using MD pharmacophore modeling and docking studies to identify potential telomerase inhibitors in our *in-house* library of arylsulfonamide compounds. Firstly, the *Tribolium castaneum* full-length catalytic subunit of telomerase in complex with the compound BIBR1532 (PDB 5CQG) was selected. The crystal structure was optimized by completing and refining the entire structure and optimizing amide groups of asparagine (Asn) and glutamine (Gln) as well as the imidazole ring in histidine (His); then, we predicted the protonation states of histidine, (His), aspartic acid (Asp), glutamic acid (Glu), and the tautomeric states of histidine. Starting from the PDB coordinates set of BIBR1532, a static pharmacophore model was created by using LigandScout containing six features (Figure 1A):


This model was used to screen the *in-house* library of arylsulfonamide derivatives.

The MYSHAPE pharmacophore model had six features (Figure 1C):

*Pharmaceuticals* **2022**, *15*, 82 3 of 22

ferent conformations.

 Three hydrophobic features; Two H-bond acceptors; One anionic feature.

The pharmacophore features were decreased to five, the two hydrophobic features on the naphthyl ring of BIBR1532 were fused into one, and the tolerance radius for the new pharmacophore feature was increased by 0.15 Å to compensate for small deviations (Figure 1B). To improve the performance of the virtual screening (VS) process, the recent MYSHAPE (Molecular dYnamics SHared PharmacophorE) approach was used [51,52]. According to this approach, the exploration of the protein conformations by molecular dynamics coupled with the pharmacophore modeling improved the result of the VS concerning the corresponding model generated from the PDB coordinates set. To build the MYSHAPE model, 20 ns of molecular dynamics simulation of the BIBR1532–protein complex was run. A new interaction was retrieved such as a hydrogen bond interaction of the carbonyl group of BIBR1532 with the Met482 by a water bridge. The new MY-SHAPE pharmacophore feature was added to the original pharmacophore model. In addition, in this case, the tolerance radius for the added pharmacophore feature was increased by 0.15 Å to compensate for small deviations in the 3D coordinates of the dif-

**Figure 1.** Pharmacophore models derived from the X-ray pose of BIBR1532 (PDB ID: 5CQG): (**A**) static pharmacophore model from the PDB with six features—4 hydrophobic features (yellow spheres), 1 H-bond acceptor (red sphere), and one anionic feature (red star); (**B**) static pharmacophore model after fusing with 5 features—3 hydrophobic features (yellow spheres), 1 H-bond acceptor (red sphere), and 1 anionic feature (red star); (**C**) MYSHAPE pharmacophore model—3 hydrophobic features (yellow spheres), 2 H-bond acceptors (red spheres), and 1 anionic feature (red star). **Figure 1.** Pharmacophore models derived from the X-ray pose of BIBR1532 (PDB ID: 5CQG): (**A**) static pharmacophore model from the PDB with six features—4 hydrophobic features (yellow spheres), 1 H-bond acceptor (red sphere), and one anionic feature (red star); (**B**) static pharmacophore model after fusing with 5 features—3 hydrophobic features (yellow spheres), 1 H-bond acceptor (red sphere), and 1 anionic feature (red star); (**C**) MYSHAPE pharmacophore model—3 hydrophobic features (yellow spheres), 2 H-bond acceptors (red spheres), and 1 anionic feature (red star).

At the same time as the pharmacophore modeling, standard precision (SP) molecular docking studies were performed using Glide [53], considering the conserved hydrophobic pocket (FVYL motif) where BIBR1532 binds. The docking studies were performed centering the docking boxes on the 3D coordinates of BIBR1532. The RMSD of BIBR1532 was calculated showing a value of 0.2 Å. The residues that characterized the The pharmacophore features were decreased to five, the two hydrophobic features on the naphthyl ring of BIBR1532 were fused into one, and the tolerance radius for the new pharmacophore feature was increased by 0.15 Å to compensate for small deviations (Figure 1B). To improve the performance of the virtual screening (VS) process, the recent MYSHAPE (Molecular dYnamics SHared PharmacophorE) approach was used [51,52]. According to this approach, the exploration of the protein conformations by molecular dynamics coupled with the pharmacophore modeling improved the result of the VS concerning the corresponding model generated from the PDB coordinates set. To build the MYSHAPE model, 20 ns of molecular dynamics simulation of the BIBR1532–protein complex was run. A new interaction was retrieved such as a hydrogen bond interaction of the carbonyl group of BIBR1532 with the Met482 by a water bridge. The new MYSHAPE pharmacophore feature was added to the original pharmacophore model. In addition, in this case, the tolerance radius for the added pharmacophore feature was increased by 0.15 Å to compensate for small deviations in the 3D coordinates of the different conformations.

The MYSHAPE pharmacophore model had six features (Figure 1C):


This model was used to screen the *in-house* library of arylsulfonamide derivatives.

At the same time as the pharmacophore modeling, standard precision (SP) molecular docking studies were performed using Glide [53], considering the conserved hydrophobic pocket (FVYL motif) where BIBR1532 binds. The docking studies were performed centering the docking boxes on the 3D coordinates of BIBR1532. The RMSD of BIBR1532 was calculated showing a value of 0.2 Å. The residues that characterized the binding site were Phe478, Met482, Met483, Arg486, Phe494, Gly495, Ile497, Trp498, Ile550, Tyr551, Gly553, Lys552, Leu554, and Arg557. BIBR1532 established an H-bond interaction and a salt bridge

between the oxygen of the amide group and Arg486 and aromatic H-bond interactions between the naphthyl ring and Phe494 and Ile550 (Figure 2). interaction and a salt bridge between the oxygen of the amide group and Arg486 and aromatic H-bond interactions between the naphthyl ring and Phe494 and Ile550 (Figure 2).

binding site were Phe478, Met482, Met483, Arg486, Phe494, Gly495, Ile497, Trp498, Ile550, Tyr551, Gly553, Lys552, Leu554, and Arg557. BIBR1532 established an H-bond

*Pharmaceuticals* **2022**, *15*, 82 4 of 22

**Figure 2.** Binding mode after re-docking of BIBR1532: the yellow, dashed lines represent the H-bond interaction and, the purple, dashed lines represent the salt bridge between the O of the carboxyl group and Arg486, and the light blue dashed lines represent the aromatic H-bond interactions between the naphthyl ring and Phe494 and Ile550. **Figure 2.** Binding mode after re-docking of BIBR1532: the yellow, dashed lines represent the H-bond interaction and, the purple, dashed lines represent the salt bridge between the O of the carboxyl group and Arg486, and the light blue dashed lines represent the aromatic H-bond interactions between the naphthyl ring and Phe494 and Ile550.

The *in-house* library was used to perform the VS by docking. MM–GBSA-binding

free energy calculations were performed for the molecules obtained from the VS and compared with BIBR1532. Hit compounds that exhibited ∆G binding values (1A, 1B, 1C, 1E, 1G, 2B, and 2C, −70.75/−62.97 kcal/mol) lower than that obtained for BIBR1532 (−62.76 kcal/mol) were selected for synthesis and in vitro tests together with hit compounds retrieved using the MYSHAPE approach (1D, 1F, and 2A) (Figure 3). The *in-house* library was used to perform the VS by docking. MM–GBSA-binding free energy calculations were performed for the molecules obtained from the VS and compared with BIBR1532. Hit compounds that exhibited ∆G binding values (1A, 1B, 1C, 1E, 1G, 2B, and 2C, −70.75/−62.97 kcal/mol) lower than that obtained for BIBR1532 (−62.76 kcal/mol) were selected for synthesis and in vitro tests together with hit compounds retrieved using the MYSHAPE approach (1D, 1F, and 2A) (Figure 3).

#### *2.2. Synthesis*

For the preparation of sulfonamide compounds of type 1 and 2 (Figure 3), several procedures are present in the literature [54–61], and a few of them were also proposed recently in light of using more environmentally friendly green chemistry approaches [62–66].

Our approach to targeting derivatives of type 1 and 2 involved reactions of the sulfonyl chloride 3 and suitable benzylamines 4 or amine 5, according to Schemes 1 and 2. The benzylamines were commercially available, whereas the 1-(4-aminophenyl)-3,5-dimethylpyrazole (5) was prepared in two steps from acetylacetone and 4-nitrophenylhydrazine and the subsequent reduction of the nitro group with H<sup>2</sup> and Pd/C, using standard literature methodology (see Section 4).

**Figure 3.** Structures of sulfonamides 1 and 2 prepared in this work (R = Me, NHAc; NO2; NH2), (R' = 4 Cl; 4-OMe; 2,5-diOMe). **Figure 3.** Structures of sulfonamides 1 and 2 prepared in this work (R = Me, NHAc; NO<sup>2</sup> ; NH<sup>2</sup> ), (R' = 4 Cl; 4-OMe; 2,5-diOMe).


**Scheme 1.** The general method for the synthesis of benzyl-sulfonamides of type 1. **Scheme 1.** The general method for the synthesis of benzyl-sulfonamides of type 1.

> S Cl

5 3 2

cause of the solubility problem of the starting material (Scheme 1).

O O

N N N In all cases, initially, the reactions were carried out under classical literature conditions, according to Methods X or Y (see Section 4) affording derivatives 1 in good to high isolated yields (64–97%) and slightly lower (49%) only in the case of 1E, probably because of the solubility problem of the starting material (Scheme 1).

**Scheme 2.** General method for the synthesis of sulfonamide 2, R = Me (2A); NO2 (2B); NH2 (2C).

In all cases, initially, the reactions were carried out under classical literature conditions, according to Methods X or Y (see Section 4) affording derivatives 1 in good to high isolated yields (64–97%) and slightly lower (49%) only in the case of 1E, probably be-

The aminophenyl-3,5-dimethylpyrazole 5 also reacted with sulfonylchloride 3 in anhydrous THF in the presence of an equimolar amount of triethylamine, affording sul-

The possibility of avoiding the use of dangerous/dry solvents was also considered. The pilot reaction of benzylamine and para-toluensulfonyl chloride, in Schot-

HN

O S O

R

N

+

fonamides 2A,B (Scheme 2).

R

NH<sup>2</sup>

R

**Scheme 1.** The general method for the synthesis of benzyl-sulfonamides of type 1.

*Pharmaceuticals* **2022**, *15*, 82 6 of 22

3 4

O

1

O

H2N

+ R'

R'

S

O O

Cl

R S NH

**Scheme 2.** General method for the synthesis of sulfonamide 2, R = Me (2A); NO2 (2B); NH2 (2C). **Scheme 2.** General method for the synthesis of sulfonamide 2, R = Me (2A); NO2 (2B); NH2 (2C).

In all cases, initially, the reactions were carried out under classical literature conditions, according to Methods X or Y (see Section 4) affording derivatives 1 in good to high isolated yields (64–97%) and slightly lower (49%) only in the case of 1E, probably be-The aminophenyl-3,5-dimethylpyrazole 5 also reacted with sulfonylchloride 3 in anhydrous THF in the presence of an equimolar amount of triethylamine, affording sulfonamides 2A,B (Scheme 2).

cause of the solubility problem of the starting material (Scheme 1). The aminophenyl-3,5-dimethylpyrazole 5 also reacted with sulfonylchloride 3 in anhydrous THF in the presence of an equimolar amount of triethylamine, affording sulfonamides 2A,B (Scheme 2). The possibility of avoiding the use of dangerous/dry solvents was also considered. The pilot reaction of benzylamine and para-toluensulfonyl chloride, in Schot-The possibility of avoiding the use of dangerous/dry solvents was also considered. The pilot reaction of benzylamine and para-toluensulfonyl chloride, in Schotten–Baumann conditions under pH control with Na2CO3, as discussed in [67], was explored, but the result was the incomplete conversion to the desired sulfonamide. Recently, a new environmentally safe methodology was reported for the preparation of a wide range of aliphatic- and aryl-sulfonamides [68], and although the methodology was not extended to investigate the reactivity of benzylamines, we decided to explore whether this kind of amine resulted in a suitable substrate to prepare sulfonamides via this route. The green method afforded yields from good to quantitative, giving in some cases a similar efficient conversion of the reactants into products, if compared to Method X (Table 1). The major concern of this methodology is the competitive reactions leading to bis-sulfonylation which, however, in the case benzylamines were not observed, suggesting that this procedure has preferable results, especially when electron-withdrawing groups are present in both the reactants (as in the case of entry 1C).


**Table 1.** Comparison of the yields obtained using classical and green chemistry methods.

\* According to [68], the amine (2 mmol) and anhydrous NaHCO<sup>3</sup> (approximately 1 g) were ground together into a fine powder, and arylsulfonyl chloride (2 mmol) was added under vigorous stirring at room temperature. The progress of the reaction was monitored by TLC until the conversion of the amine was completed.

Compound 2C was prepared by reduction of the corresponding nitro compound with H2/Pd in ethanol a quantitative yield.

#### *2.3. Cytotoxic Activity and Selectivity Index*

To select the most promising cytotoxic agents, the newly synthesized compounds were initially evaluated on human colorectal carcinoma (HCsT-116) [69,70], human breast adenocarcinoma (MCF-7) [70,71], and human chronic myeloid leukemia (K-562) [71,72] cell lines that express the active telomerase.

As shown in Table 2, except for compound 2C, the antiproliferative efficacy of the compounds was modest and only evident at the highest micro-submillimolar concentrations. Moreover, compounds 1B, 1D, 1E, 1F, 1G, and 2B exhibited solubility issues that prevented the evaluation of their cytotoxicity.

*Pharmaceuticals* **2022**,*15*, 82

**Table 2.** Cytotoxicity of the synthesized compounds on human tumor cell lines HCT-116, MCF-7, and K-562. Cells were treated for 24 h in the absence or the presence of the compound tested at the indicated concentration, and cell viability was assessed by an MTT test as reported in Section 4. Values are expressed as the mean ±SD of three separate experiments conducted in triplicate. \* Sol: solubility issue.


On the other hand, compound 2C revealed an interesting activity, and it was possible to calculate its IC<sup>50</sup> (Table 3). Interestingly, this value on K-562 cells appeared 3.6-fold lower than that on HCT-116 and 4.2-fold lower than that on MCF-7 cells. Moreover, in the K-562 cell line tested, compound 2C showed an IC<sup>50</sup> 6.8-fold lower than the reference compound BIBR1532, and an IC<sup>50</sup> comparable to BIBR1532 in the MCF-7 cell line.

**Table 3.** IC<sup>50</sup> of 2C against HCT-116, MCF-7, K-562, and BALB/3-T3 cell lines. Values were calculated by plotting the percentage viability versus concentration on a logarithmic graph and are expressed as the mean ± SD of three separate experiments conducted in triplicate.


Selective cytotoxicity is a pivotal requirement for anticancer drugs. To determine the selectivity of compound 2C, its cytotoxicity against the cancer cell lines employed (i.e., HCT-116, MCF-7, and K-562) was compared with that against the non-cancerous, murine, embryonic, fibroblast cell line BALB/3-T3. As shown in Table 3, the IC<sup>50</sup> of compound 2C on BALB/3-T3 cells was remarkably higher than those on the tumoral cell lines tested. Moreover, the calculated selectivity indexes (SIs) of compound 2C for HCT-116, MCF-7, and K-562 cells were 2.9, 2.5, and 9.8, respectively. Relevantly, these values are above the accepted threshold for antitumor drugs (SI = 2.0) [72,73].

Although selectivity for cancer cells cannot be easily derived from the comparison of toxicity parameters in different cell cultures, these data indicate that compound 2C shows preferential toxicity towards cancer cells. Along these lines, the potential use of compound 2C as a novel lead molecule for the development of more potent and selective antiproliferative agents can, therefore, be envisaged.

#### *2.4. Docking of Compound 2C*

Analysis of the best docking pose of compound 2C showed a pi–pi stacking interaction between Phe494 and the central aromatic ring as observed in BIBR1532, where the same residue is involved in two H-aromatic bonds with the naphthyl ring (Figure 2). The same residues Phe494, Asp493, and Gly495 establish positive van der Waals (vdW) contacts with the aniline group of 2C. Other positive vdW contacts are formed by the sulfonamide moiety of 2C and Ile550, Gly553, and Leu554. The major part of positive vdW contacts are established by the dimethylpyrazolyl moiety with Met482, Met483, Arg486, Phe494, and Ile550 (Figure 4).

The same residues Phe494, Asp493, and Gly495 establish positive van der Waals (vdW) contacts with the aniline group of 2C. Other positive vdW contacts are formed by the sulfonamide moiety of 2C and Ile550, Gly553, and Leu554. The major part of positive vdW contacts are established by the dimethylpyrazolyl moiety with Met482, Met483,

*Pharmaceuticals* **2022**, *15*, 82 9 of 22

Arg486, Phe494, and Ile550 (Figure 4).

**Figure 4.** Binding mode of the best docking pose of 2C after the docking study. The light blue, dashed lines represent the pi–pi stacking. The displayed residues establish positive vdW contacts with the compound. **Figure 4.** Binding mode of the best docking pose of 2C after the docking study. The light blue, dashed lines represent the pi–pi stacking. The displayed residues establish positive vdW contacts with the compound.

#### *2.5. Molecular Dynamics Simulation*

*2.5. Molecular Dynamics Simulation*  The running of dynamics simulations of protein–ligand complexes over time could be considered a major step toward accuracy in computer-assisted drug design. In this study, an unbiased molecular dynamics simulation was performed to investigate the conformational stability and the time-dependent binding capability of 2C in the active site of the telomerase. Additionally, we tried to understand if the protein target undergoes conformational alteration after interacting with 2C. Therefore, starting from the previous docking, 100 ns of MD simulation was carried out. Various analyses, such root mean square deviation (RMSD), root mean square fluctuation (RMSF), and determination of the number and types of protein–ligand contacts, were performed to obtain a The running of dynamics simulations of protein–ligand complexes over time could be considered a major step toward accuracy in computer-assisted drug design. In this study, an unbiased molecular dynamics simulation was performed to investigate the conformational stability and the time-dependent binding capability of 2C in the active site of the telomerase. Additionally, we tried to understand if the protein target undergoes conformational alteration after interacting with 2C. Therefore, starting from the previous docking, 100 ns of MD simulation was carried out. Various analyses, such root mean square deviation (RMSD), root mean square fluctuation (RMSF), and determination of the number and types of protein–ligand contacts, were performed to obtain a more detailed analysis of the 2C–target complex.

#### more detailed analysis of the 2C–target complex. 2.5.1. Stability Analysis

2.5.1. Stability Analysis The RMSD was selected as a criterion to evaluate the dynamic stability of the ligand-bound system [73]. The RMSD values of the protein's atoms and ligand are reported The RMSD was selected as a criterion to evaluate the dynamic stability of the ligand-bound system [73]. The RMSD values of the protein's atoms and ligand are reported in Figure 5. The system reached equilibrium quickly and fluctuated around the average RMSD value of <3 Å. The average value of ligand vs. protein RMSD of ~6.4 Å indicated a strong stability of 2C in the binding pocket compared to the ligand vs. ligand RMSD of ~1.81 Å.

~1.81 Å

~1.81 Å

Å.

*Pharmaceuticals* **2022**, *15*, 82 10 of 22

*Pharmaceuticals* **2022**, *15*, 82 10 of 22

**Figure 5.** Compound 2C and protein RMSD during MD simulation (100 ns). The *x*-axis is expressed by the number of snapshots extracted, and 1 snapshot = 100 ps; the *y*-axis is expressed in **Figure 5.** Compound 2C and protein RMSD during MD simulation (100 ns). The *x*-axis is expressed by the number of snapshots extracted, and 1 snapshot = 100 ps; the *y*-axis is expressed in Å. Å.

in Figure 5. The system reached equilibrium quickly and fluctuated around the average RMSD value of <3 Å. The average value of ligand vs. protein RMSD of ~6.4 Å indicated a strong stability of 2C in the binding pocket compared to the ligand vs. ligand RMSD of

in Figure 5. The system reached equilibrium quickly and fluctuated around the average RMSD value of <3 Å. The average value of ligand vs. protein RMSD of ~6.4 Å indicated a strong stability of 2C in the binding pocket compared to the ligand vs. ligand RMSD of

2.5.2. Residue Mobility and Protein–Ligand Contact Analyses 2.5.2. Residue Mobility and Protein–Ligand Contact Analyses To examine the structural flexibility effect of 2C upon TERT per residue, the main

2.5.2. Residue Mobility and Protein–Ligand Contact Analyses To examine the structural flexibility effect of 2C upon TERT per residue, the main chain average of the root mean square fluctuation (RMSF) of the complex was calculated for the entire 100 ns of simulation. The residue-wise fluctuation of the complexes was plotted and is presented in Figure 6. As reported, the RMSF plot was low for the identi-To examine the structural flexibility effect of 2C upon TERT per residue, the main chain average of the root mean square fluctuation (RMSF) of the complex was calculated for the entire 100 ns of simulation. The residue-wise fluctuation of the complexes was plotted and is presented in Figure 6. As reported, the RMSF plot was low for the identified residues involved in the interactions with 2C in the docking analysis. chain average of the root mean square fluctuation (RMSF) of the complex was calculated for the entire 100 ns of simulation. The residue-wise fluctuation of the complexes was plotted and is presented in Figure 6. As reported, the RMSF plot was low for the identified residues involved in the interactions with 2C in the docking analysis.

**Figure 6.** Protein RMSF during MD simulation (100 ns). The *x*-axis reports the number of the resi-**Figure 6.** Protein RMSF during MD simulation (100 ns). The *x*-axis reports the number of the residues of TERT, and the *y*-axis is expressed in Å. **Figure 6.** Protein RMSF during MD simulation (100 ns). The *x*-axis reports the number of the residues of TERT, and the *y*-axis is expressed in Å.

dues of TERT, and the *y*-axis is expressed in Å. Evaluation of the protein interactions provides a measure of the interaction power between the ligands and the protein and can be categorized by type and summarized as represented in Figure 7. In Figure 7A, the interactions that occurred for more than 10% of the simulation time are reported. In Figure 7B, a timeline representation of the interactions and contacts (H-bonds, hydrophobic, ionic, and water bridges) summarize the total number of specific contacts the protein makes with the ligand throughout the simulation. The bottom panel of Figure 7B shows which residues interact with the ligand in each trajectory frame. Some residues make more than one specific contact with the lig-Evaluation of the protein interactions provides a measure of the interaction power between the ligands and the protein and can be categorized by type and summarized as represented in Figure 7. In Figure 7A, the interactions that occurred for more than 10% of the simulation time are reported. In Figure 7B, a timeline representation of the interactions and contacts (H-bonds, hydrophobic, ionic, and water bridges) summarize the total number of specific contacts the protein makes with the ligand throughout the simulation. The bottom panel of Figure 7B shows which residues interact with the ligand in each trajectory frame. Some residues make more than one specific contact with the lig-Evaluation of the protein interactions provides a measure of the interaction power between the ligands and the protein and can be categorized by type and summarized as represented in Figure 7. In Figure 7A, the interactions that occurred for more than 10% of the simulation time are reported. In Figure 7B, a timeline representation of the interactions and contacts (H-bonds, hydrophobic, ionic, and water bridges) summarize the total number of specific contacts the protein makes with the ligand throughout the simulation. The bottom panel of Figure 7B shows which residues interact with the ligand in each trajectory frame. Some residues make more than one specific contact with the ligand, and they are represented by a darker shade of orange, according to the scale to the right of the plot. Protein–ligand contacts are categorized into four types: hydrogen bonds, hydrophobic, ionic, and water bridges. The stacked bar charts are normalized throughout the trajectory: a value of 1.0 suggests that 100% of the simulation time, the specific interaction is maintained.

specific interaction is maintained.

*Pharmaceuticals* **2022**, *15*, 82 11 of 22

and, and they are represented by a darker shade of orange, according to the scale to the right of the plot. Protein–ligand contacts are categorized into four types: hydrogen bonds, hydrophobic, ionic, and water bridges. The stacked bar charts are normalized throughout the trajectory: a value of 1.0 suggests that 100% of the simulation time, the

**Figure 7.** Protein–compound 2C contacts: (**A**) the interactions that occur for more than 10% of the simulation time are reported; (**B**) a timeline representation of the interactions and contacts (H-bonds, hydrophobic, ionic, and water bridges) during the simulation; (**C**) histogram of protein–compound 2C interactions fractions: hydrogen bonds (green), hydrophobic (purple), ionic (red), and water bridges (blue). **Figure 7.** Protein–compound 2C contacts: (**A**) the interactions that occur for more than 10% of the simulation time are reported; (**B**) a timeline representation of the interactions and contacts (H-bonds, hydrophobic, ionic, and water bridges) during the simulation; (**C**) histogram of protein–compound 2C interactions fractions: hydrogen bonds (green), hydrophobic (purple), ionic (red), and water bridges (blue).

The analysis of the trajectory evidenced that the pi-stacking interaction between the aniline ring of 2C and Phe494 occurred at 57%. This output evidences a shifting in the pi-stacking interaction with respect to the docking where the central phenyl ring was involved, but it clearly defines Phe494 as a key residue in the inhibition pattern. As also reported in Figure 7B,C, Phe494 established at least one contact for the entire duration of the simulation and, for a major part of the time, 2–3 contacts with the phenyl ring of 2C. Gly553 showed for 57% and 17% of the time two H-bonds interactions with the oxygen of the sulfonamide group emerging after ~18 ns of simulation. Interestingly, after 22 ns of the simulation a new in pi–cation interaction appeared for 14% of the time between Arg557 and the pyrazole ring, even though it was in a spotted fashion. During the first 20 ns (16% of the total time), the protein–ligand interaction stabilized with an H-bond between Ile550 and the NH of the sulfonamide group. It is interesting evidence of the role of a water molecule as a water bridge between Trp498 and the pyrazole ring. Other residues involved in the protein–ligand interaction but with a minor role were Arg486, Ile497, Trp498, Tyr551, Lys552, and Leu554. The analysis of the trajectory evidenced that the pi-stacking interaction between the aniline ring of 2C and Phe494 occurred at 57%. This output evidences a shifting in the pi-stacking interaction with respect to the docking where the central phenyl ring was involved, but it clearly defines Phe494 as a key residue in the inhibition pattern. As also reported in Figure 7B,C, Phe494 established at least one contact for the entire duration of the simulation and, for a major part of the time, 2–3 contacts with the phenyl ring of 2C. Gly553 showed for 57% and 17% of the time two H-bonds interactions with the oxygen of the sulfonamide group emerging after ~18 ns of simulation. Interestingly, after 22 ns of the simulation a new in pi–cation interaction appeared for 14% of the time between Arg557 and the pyrazole ring, even though it was in a spotted fashion. During the first 20 ns (16% of the total time), the protein–ligand interaction stabilized with an H-bond between Ile550 and the NH of the sulfonamide group. It is interesting evidence of the role of a water molecule as a water bridge between Trp498 and the pyrazole ring. Other residues involved in the protein–ligand interaction but with a minor role were Arg486, Ile497, Trp498, Tyr551, Lys552, and Leu554.

#### 2.5.3. ADME Calculation for Compound 2C 2.5.3. ADME Calculation for Compound 2C

ADME calculation for compound 2C was performed using the SwissADME web tool [74]. In the hexagon drug-likeness graph (Figure 8), each vertex represents a param-ADME calculation for compound 2C was performed using the SwissADME web tool [74]. In the hexagon drug-likeness graph (Figure 8), each vertex represents a parameter that defines a bioavailable drug. The pink regions represent the optimum range of the following six properties: lipophilicity = 2.53 (XLOGP3 between −0.7 and +5.0), size = 342.42 g/mol (MW between 150 and 500 g/mol), polarity = 98.39 Å<sup>2</sup> (TPSA between 20 and 130 Å<sup>2</sup> ), solubility (log S not higher than 6), saturation = 0.12 (fraction of carbons in the sp3 hybridization not less than 0.25), and flexibility (no more than nine rotatable bonds). It was found that compound 2C was slightly outside the pink area on one side

due to the inconformity of its insaturation. The pharmacokinetics analysis showed that compound 2C is probably not a P-glycoprotein (P-gp) substrate and not BBB permeant. It has potentially good gastrointestinal (GI) absorption and ABS, up to nearly 70%. The %ABS, a very functional physiochemical variable, defines a drug's transport properties. It was calculated according to the equation %ABS = 109 − (0.345 × TPSA) [75,76]. TPSA values below 98.39 Å<sup>2</sup> characterize a significant permeability in the cellular plasma membrane. It has a bioavailability score of 0.55, which means good pharmacokinetic properties according to the Rule of 5 by Lipinski [77]. compound 2C is probably not a P-glycoprotein (P-gp) substrate and not BBB permeant. It has potentially good gastrointestinal (GI) absorption and ABS, up to nearly 70%. The %ABS, a very functional physiochemical variable, defines a drug's transport properties. It was calculated according to the equation %ABS = 109 − (0.345 × TPSA) [75,76]. TPSA values below 98.39 Ǻ<sup>2</sup> characterize a significant permeability in the cellular plasma membrane. It has a bioavailability score of 0.55, which means good pharmacokinetic properties according to the Rule of 5 by Lipinski [77].

eter that defines a bioavailable drug. The pink regions represent the optimum range of the following six properties: lipophilicity = 2.53 (XLOGP3 between −0.7 and +5.0), size =

and 130 Å2), solubility (log S not higher than 6), saturation = 0.12 (fraction of carbons in the sp3 hybridization not less than 0.25), and flexibility (no more than nine rotatable bonds). It was found that compound 2C was slightly outside the pink area on one side due to the inconformity of its insaturation. The pharmacokinetics analysis showed that

(TPSA between 20

*Pharmaceuticals* **2022**, *15*, 82 12 of 22

342.42 g/mol (MW between 150 and 500 g/mol), polarity = 98.39 Å<sup>2</sup>

**Figure 8.** The bioavailability radar (the pink area exhibits the optimal range of a particular property) for compound 2C (LIPO = lipophilicity as in XLOGP3; SIZE = size as in molecular weight; POLAR = polarity as TPSA (topological polar surface area); INSOLU = insolubility in water by log S scale; INSATU = insaturation per fraction of carbons in the sp3 hybridization; FLEX = flexibility per rotatable bonds). **Figure 8.** The bioavailability radar (the pink area exhibits the optimal range of a particular property) for compound 2C (LIPO = lipophilicity as in XLOGP3; SIZE = size as in molecular weight; POLAR = polarity as TPSA (topological polar surface area); INSOLU = insolubility in water by log S scale; INSATU = insaturation per fraction of carbons in the sp3 hybridization; FLEX = flexibility per rotatable bonds).

#### **3. Discussion 3. Discussion**

In the past decade, there has been tangible progress in the definition of the role of telomerase in tumor progression. Several studies evidenced telomerase as a promising target for developing cancer therapeutics. The main reason is due to the overexpression of telomerase in cancer cells (85–90%) compared with normal cells where it is almost unexpressed. Despite the fact that efforts have allowed for the identification of several small molecules, oligonucleotides, natural products, and immunotherapeutics, no telomerase-based cancer drugs have yet been approved by the FDA due to the long lag time between administration of the drug and clinical response. For these reasons, the identification and then the optimization of new small molecules, such as telomerase in-In the past decade, there has been tangible progress in the definition of the role of telomerase in tumor progression. Several studies evidenced telomerase as a promising target for developing cancer therapeutics. The main reason is due to the overexpression of telomerase in cancer cells (85–90%) compared with normal cells where it is almost unexpressed. Despite the fact that efforts have allowed for the identification of several small molecules, oligonucleotides, natural products, and immunotherapeutics, no telomerasebased cancer drugs have yet been approved by the FDA due to the long lag time between administration of the drug and clinical response. For these reasons, the identification and then the optimization of new small molecules, such as telomerase inhibitors, remain a crucial point that deserves to be further explored. A structure-based drug design approach could be used to design potential inhibitors of the telomerase active site. Exploiting the experimental structure of TERT showing a known inhibitor co-crystallized, we developed a combined structure-based strategy to screen an *in-house* library of 126 arylsulfonamide derivatives. The combined computational strategy was performed using the MYSHAPE approach, first, to identify hit compounds. Successively, docking and calculation of the binding free energy using the MM-GBSA method identified other hits. These approaches guided the selection of the most promising compounds for synthesis and in vitro tests.

As described in Section 2.2., the compounds were initially prepared using classical methodologies. To explore the possibility of avoiding using dangerous and/or dry solvents, a reported green procedure was successfully utilized in the case of derivatives listed in Table 1, affording yields from good to quantitative for a major part of the compounds.

With the aim of identifying the most promising cytotoxic agents, the selected hits, newly synthesized, were evaluated on three cancer lines: human colorectal carcinoma (HCT-116), human breast adenocarcinoma (MCF-7), and human chronic myeloid leukemia (K-562) cell lines, considering the known telomerase inhibitor BIBR1532 as the reference compound. One of the major problems faced during the in vitro test was the poor solubility for some of the tested compounds, which prevented the evaluation of the cytotoxicity and gave results in micro-submillimolar concentrations. In general, the pyrazole derivatives (2A, 2B, and 2C) seems to show a higher activity on the tested cell lines, probably even due to the major solubility

One of the selected hits (2C) revealed an interesting IC<sup>50</sup> in all three cell lines. In particular, the IC<sup>50</sup> results were comparable with BIBR1532 for the HCT-116 cell line and the MCF7 cell line, but several times lower in the K-562 cell lines. These findings were confirmed by the evidence of the selective cytotoxicity against the non-cancerous, murine, embryonic, fibroblast cell line BALB/3-T3. The IC<sup>50</sup> of compound 2C on BALB/3-T3 cells was remarkably higher than those on the tumoral cell lines tested. Another interesting concern regarding the data was the absence of the solubility issue encountered during the in vitro test for the compound 2C in addition to the calculated ADME values that seemed to exclude a binding with the P-gp, an important parameter for the insurgence of resistance issues. In the end, compound 2C showed a good calculated pharmacokinetic profile. Further, the MD simulation of compound 2C at the telomerase active site showed good stability and evidenced Phe494 as the key residue also in BIBR1532. But the more interesting evidence of the simulation regarded that each moiety of compound 2C was involved in key interactions with some residues of the active site: the pyrazole moiety in a cation–pi-stacking interaction with Arg557, the two phenyl rings in the pi-stacking interaction with Phe494, and the sulfonamide moiety in the H-bond interaction with Ile550 and Gly553. Compound 2C as arylsulfonamide derivative is an interesting hit compound that deserves, for the reasons stated above, further investigation in terms of the optimization of its structure to obtain more active telomerase inhibitors.

#### **4. Materials and Methods**

#### *4.1. Library and Protein Preparation*

The library used for this experimental work comprised a series of 126 sulfonamide derivates. The sulfonamide and BIBR1532 structures were built using the builder panel in Maestro, and ligand preparation was carried out by LigPrep (Schrödinger, LLC, New York, NY, USA, 2021). The force field adopted was OPLS4 [78] and Epik 5.5 (Schrödinger, 2021-1) was selected as an ionization tool at pH 7.2 ± 0.2. Tautomer generation was unflagged, and the maximum number of conformers generated was set at 32. For the purpose of this study, the 2.30 Å resolution crystal structure of *Tribolium castaneum* telomerase in complex with the highly specific inhibitor BIBR1532 (PDB ID: 5CQG) [21] was used. The structure was optimized using the Protein Preparation Wizard in Maestro (Schrödinger, 2021-1) adding bond orders and hydrogen atoms to the crystal structure using the OPLS4 force field. The remaining settings are reported in [79].

#### *4.2. Docking Studies*

The library was submitted to a docking study using Glide v9.0 [80] in standard precision (SP) with the OPLS4 force field. The grid box was built considering the BIBR1532 as the centroid of the grid. The study was performed using no constraints. Van der Waals radii were set at 0.8, and the partial cutoff was 0.15 with flexible ligand sampling. Bias sampling torsion penalization for amides with non-planar conformation and Epik state penalties were added to the docking score.

#### *4.3. MYSHAPE Approach*

MYSHAPE [51] is an approach in which the pharmacophore model is built using only the common pharmacophore feature patterns that the ligands exhibit during MD simulations. The evaluation of the ligand–protein interaction patterns during the MD simulation was investigated using the data recorded employing the MD simulation (exactly 1001 frames for each simulation). The MYSHAPE model was used to distinguish this type of model from the "default" pharmacophore models that were generated using the crystal structure of the ligand–protein complex. For the newly added pharmacophore features, the tolerance radius was increased by 0.15 Å to compensate for small deviations in the 3D coordinates. The computational details are reported in previously published papers [51,52]

#### *4.4. MM-GBSA Free Energy Calculations*

The output of docking was used to calculate the ∆G values of the BIBR1532 and molecules using MM-GBSA (molecular mechanics generalized-born surface area) (Prime, Schrödinger, LLC, New York, NY, USA, 2021). Protein–ligand binding free energy using MM-GBSA was calculated as the difference between the energy of the bound complex and the energy of the unbound protein and ligand. The entropy term −T∆S was not calculated to reduce computational time as previously reported [79,81–83]. The VSGB solvation model was chosen using OPLS4 force field with the minimized sampling method.

#### *4.5. Molecular Dynamic Simulations*

A 100 ns MD simulation was carried out using a Desmond 6.5 (Desmond Molecular Dynamics System, D. E. Shaw Research, New York, NY, USA, 2021) using the OPLS4 force field for the complex TERT/2C. The complex was solvated in a cubic box using the TIP3P. Ions were added to neutralize charges. The systems were minimized and equilibrated at a temperature of 303.15 K and a pressure of 1.013 bar. The system was simulated as an NPT ensemble; a Nose–Hover thermostat and Martyna–Tobia–Klein barostat were used. The integration time step was chosen to be 2 fs. To keep the hydrogen-heavy atom bonds rigid, the SHAKE algorithm was used. A 9 Å cutoff radius was set for the short-range Coulomb interactions, and smooth particle mesh Ewald was used for the long-range interactions. The detection ranges for energy were 1.2 ps, and 5.0 ps for the trajectory.

#### *4.6. ADME Prediction*

Drug-likeness, physicochemical properties, lipophilicity, solubility, and pharmacokinetics properties were analyzed by the SwissADME web tool (http://www.swissadme.ch/ index.php (accessed on 30 November 2021)). SwissADME is a new comprehensive tool run by the Swiss Institute of Bioinformatics (SIB) enabling estimation of ADME (absorption, distribution, metabolism, and excretion) parameters of drug candidates. 2D structural models of the compound were drawn in the molecular sketcher into the ChemAxon's Marvin JS window and transferred to the SMILES (simplified molecular-input line-entry system) format to predict suitable properties [74].

#### *4.7. Chemistry*

All melting points were taken on a Büchi melting point M-560 apparatus and were uncorrected. IR spectra were recorded in Bromoform with a JASCO FT-IR spectrophotometer. 1H and 13C NMR spectra were measured at 200 and 50.0 MHz, respectively, in DMSO-d<sup>6</sup> solution and TMS as an internal standard, using a Bruker Avance II Series 200 MHz spectrometer or at 300 and 75 MHz (APT) with a Bruker AC-E spectrometer. Column chromatography was performed with Merck silica gel (230–400 mesh ASTM). Elemental analyses (C, H, and N) were within ±0.4% of theoretical values. The substituted benzylamine derivatives were commercially available and were used without further purification. For all the compounds already cited in the literature, the IR and NMR spectra results were identical to those reported (1A–E). The 1H and 13C NMR spectra of new synthesized compounds (1F,G and 2A–C) are reported in the Supplementary Materials.

4.7.1. General Methods for the Preparation of *N*-(R'-Benzyl)-4-R-benzenesulfonamides of Type 1

Method X: To a stirred solution of amino derivative (1.5 mmol) in anhydrous THF (20 mL), an equimolar amount of triethylamine (0.21 mL) and the suitable sulfonylchloride were added. The reaction mixture was stirred for the appropriate time-lapse until the disappearance of the starting amine (TLC monitoring). The ammonium salt was collected and taken up with 20 mL H2O. The aqueous solution was extracted using DCM. The organic layers were combined with the organic mother liquor, dried over Na2SO4, and concentrated under reduced pressure to yield the desired sulfonamide of type 1.

Method Y: The amino derivative (5 mmol) was dissolved in dry pyridine, and the suitable sulfonylchloride (5 mmol) was added. The reaction mixture was stirred for the appropriate time-lapse until the disappearance of starting amine (TLC monitoring). The mixture was poured onto ice/water and the solid precipitate was collected by filtration, air-dried, and recrystallized from ethanol.

#### *N*-(4-Chlorobenzyl)-4-methylbenzenesulfonamide (1A)

This compound was prepared according to Method X. The reaction mixture was stirred for 24 h at room temperature. Yield 70%. Mp 98 ◦C ([84] 107–108 ◦C; [85] 106.1 ◦C). 1H NMR ppm: 2.43 (s, 3H, Me), 4.06 (d, J = 6.4 Hz, 2H, CH2), 5.18 (t, J = 6.0 Hz, 1H, NH), 7.11 (d, J = 8.4 Hz, 2H, Ar-H), 7.20 (d, J = 8.4 Hz, 2H, Ar-H), 7.27 (d, J = 8.0 Hz, 2H, Ar-H), 7.70 (d, J = 8.4 Hz, 2H, Ar-H). 13C NMR ppm: 21.6, 46.6, 127.2, 128.8, 129.3, 129.9, 133.7, 135.1, 136.9, and 143.8.

#### *N*-(4-Chlorobenzyl)-4-acetylaminobenzenesulfonamide (1B)

This compound was prepared according to Method X. The reaction mixture was stirred for 24 h at room temperature and under reflux for an additional 4 h. Yield 64%. Mp 152◦C [56] 172–174 ◦C). 1H NMR ppm: 2.10 (s, 3H, Me), 3.08 (bs, 1H, NH), 3.96 (s, 2H, CH2), 7.30–8.11 (m, 8H, Ar-H), 10.42 (bs, 1H, NH). 13C NMR ppm: 24.1, 45.5, 118.6, 127.6, 128.5, 129.4, 131.6, 134.1, 136.9, 142.5, and 169.0.

#### *N*-(4-Chlorobenzyl)-4-nitrobenzenesulfonamide (1C)

This compound was prepared according to Method X. The reaction mixture was stirred under reflux for 4.5 h. Yield 78%. Mp 174–176 ◦C ([86] 168–170 ◦C). 1H NMR ppm: 3.80 (d, J = 6.0 Hz, 2H, CH2), 6.61 (d, J = 8.0 Hz, 2H, Ar-H), 6.85 (d, J = 8.0 Hz, 2H, Ar-H), 7.15 (d, J = 8.0 Hz, 2H, Ar-H), 7.45 (d, J = 8.0 Hz, 2H, Ar-H), 7.56 (t, J = 6.0 Hz, 1H, NH). 13C NMR ppm: 45.5, 112.6, 113.5, 125.5, 128.4, 128.9, 129.8, 152.4, and 158.3.

#### *N*-(4-Methoxybenzyl)-4-methylbenzenesulfonamide (1D)

This compound was prepared according to Method X. The reaction mixture was stirred for 3 h at room temperature. Yield 71%. Mp 122–123 ◦C ([87] 114–117 ◦C; [88] 122–123 ◦C from EtOAc). 1H NMR (300 MHz, DMSO) ppm: 2.38 (s, 3H, Me), 3.72 (s, 3H, Me), 3.89 (s, 2H, CH2), 6.85 (d, J = 8.4 Hz, 2H, Ar-H), 7.17 (d, J = 8.4 Hz, 2H, Ar-H), 7.38 (d, J = 8.0 Hz, 2H, Ar-H), 7.71 (d, J = 8.4 Hz, 2H, Ar-H), 8.03 (bs, 1H, NH). 13C NMR ppm: 20.9, 45.8, 55.0, 113.6, 126.5, 128.9, 129.4, 129.5, 137.8, 142.5, and 158.4.

#### *N*-(4-Methoxybenzyl)-4-acetylaminobenzenesulfonamide (1E)

This compound was prepared according to Method Y. The reaction mixture was stirred for 1 h at room temperature and under reflux for an additional 2 h. Yield 49%. Mp 192– 193 ◦C ([87] 164–167 ◦C). 1H NMR ppm: 1.84 (s, 2H, CH2), 2.12 (s, 3H, Me), 4.46 (sa, 1H, NH), 7.74 (d, J = 8.5 Hz, 2H), 7.79 (d, J = 8.5 Hz, 2H), 9.41 (sa, 1H, NH). 13C NMR ppm: 24.2, 45.6, 55.2, 113.6, 118.6, 127.6, 129.1, 129.5, 134.0, 142.9, 158.5, and 169.4

#### *N*-(2,5-Dimethoxybenzyl)-4-methylbenzenesulfonamide (1F)

This compound was prepared according to Method Y. The reaction mixture was stirred for 10 h at room temperature. Yield 80%. Mp 110–111 ◦C. IR ν: 3268 (NH) cm−<sup>1</sup> . 1H-NMR ppm: 2.37 (3H, s, Me), 3.65 (6H, s, Me), 3.91 (2H, t, J = 3.9 Hz, CH2), 6.77–6.83 (3H, m, ArH), 7.10 (2H, d, J = 8.5 Hz), 7.36 (2H, d, J = 7.6 Hz), 7.68 (2H, d, J = 7.6 Hz), 7.89 (1H, t, J = 3.9 Hz, NH). 13C NMR ppm: 21.4, 41.3, 55.8, 56.1, 111.8, 113.0, 115.0, 126.7, 127.0, 130.0, 138.2, 143.0, 150.9, and 153.4.

#### *N*-(2,5-Dimethoxybenzyl)-4-nitrobenzenesulfonamide (1G)

This compound was prepared according to Method X. The reaction mixture was stirred for 24 h at room temperature. Yield 97%. Mp 129–132◦C. ir ν: 3268 (NH), 1523 and 1348 (NO2) cm−<sup>1</sup> . 1H NMR ppm: 3.35 (s, 3H, Me), 3.63 (s, 3H, Me), 4.03 (s, 2H, CH2), 6.65–6.86 (m, 3H, Ar), 7.96 (2H, d, J = 8.8 Hz), 8.34 (2H, d, J = 8.8 Hz). 13C NMR ppm: 46.0, 55.7, 56.1, 111.8, 113.2, 115.5, 124.7, 125.9, 128.4, 147.0, 149.8, 151.0, and 153.2.

#### 4.7.2. Preparation of 1-(4-Aminophenyl)-3,5-dimethylpyrazole

A solution of acetylacetone (10 mmol) and 4-nitrophenylhydrazine (10 mmol) in acetic acid was heated under reflux for 9 h. The reaction mixture was cooled to room temperature and poured onto ice/water to give a solid precipitate. Purification by column chromatography (eluant DCM/EtOAc 4:1) gave the nitro compound. Yield 80%, mp 100 ◦C ([89] 102 ◦C from EtOH). 1-(4-Nitrophenyl)-3,5-dimethylpyrazole (0.56 mg) dissolved in EtOH (20 mL) Pd/C was added, and the mixture was reduced overnight at room temperature in an H2 atmosphere (50 psi) in a Parr apparatus. The catalyst was filtered off, and the solution was evaporated under reduced pressure. The amino derivatives were isolated as white crystals. Yield 100%. Mp 84 ◦C ([90] 82–84 ◦C from benzene).

4.7.3. General Method for the Preparation of *N*-[4-(3,5-Dimethyl-1H-pyrazol-1-yl)phenyl]- 4- R-benzenesulfonamide of Type 2

To a stirred solution of 1-(4-aminophenyl)-3,5-dimethylpyrazole (1.5 mmol) in anhydrous THF (20 mL), an equimolar amount of triethylamine (0.21 mL) and the suitable sulfonylchloride were added. The reaction mixture was stirred for the appropriate timelapse until the disappearance of the starting amine (TLC monitoring). The ammonium salt was collected and taken up with 20 mL H2O. The aqueous solution was extracted using DCM. The organic layers were combined with the organic mother liquor, dried over Na2SO4, and concentrated under reduced pressure to yield the desired sulfonamide.

#### *N*-[4-(3,5-Dimethylpyrazol-1-yl)phenyl]-4-methylbenzenesulfonamide (2A)

The reaction mixture was stirred for 30 h at room temperature. Yield 96%. Mp 167 ◦C. ir ν: 3392 (NH) cm−<sup>1</sup> . 1H NMR ppm: 2.09 (s, 3H, Me), 2.13 (s, 3H, Me), 2.28 (s, 3H, Me), 6.00 (s, 1H, pyrazole-CH), 7.20 (dd, 4H, ArH), 7.27 (d, 2H, ArH), 7.64 (d, 2H, ArH), 10.41 (s, 1H, NH). 13C NMR ppm: 12.3, 13.4, 21.3, 107.5, 120.7, 125.5, 125.8, 127.2, 128.8, 130.3, 136.7, 134.0, 144.2, and 148.6.

#### *N*-[4-(3,5-Dimethylpyrazol-1-yl)phenyl]-4-nitrobenzenesulfonamide (2B)

The reaction mixture was stirred for 3 h at room temperature. The solid residue was recrystallized from ethanol. Yield 57%. Mp 219 ◦C. ir ν: 3392 (NH), 1532, and 1348 (NO2) cm−<sup>1</sup> . 1H NMR ppm: 2.12 (s, 3H, Me), 2.21 (s, 3H, Me), 6.02 (s, 1H, pyrazole-CH), 7.27 (dd, 4H, ArH), 8.19 (dd, 4H, ArH), 10.81 (s, 1H, NH). 13C NMR ppm: 12.0, 13.2, 107.0, 120.9, 124.7, 125.0, 128.3, 132.2, 135.5, 136.2, 139.0, 147.7, and 149.9.

#### Preparation of *N*-[4-(3,5-Dimethylpyrazol-1-yl)phenyl]-4-aminobenzenesulfonamide (2C)

To the nitro derivative, 2B dissolved was in EtOH (20 mL), Pd/C was added, and the mixture was reduced overnight at room temperature in an H<sup>2</sup> atmosphere (50 psi) in a Parr apparatus. The catalyst was filtered off, and the solution was evaporated under reduced pressure. The amino derivative was isolated as white crystals. Yield 100%. Mp 106–110 ◦C. ir ν: 3482 and 3380 (NH2), 3234 (NH) cm−<sup>1</sup> . 1H NMR ppm: 1.92 (s, 3H, Me), 2.33 (s, 3H, Me), 6.00 (s, 1H, pyrazole-CH), 7.30–7.35 (m, 4H, Ar-H), 7.42–7.46 (m, 4H, Ar-H), 13C NMR ppm: 12.0, 13.2, 112.6, 119.5, 122.6, 124.8, 125.0, 128.3, 130.3, and 152.9.

#### *4.8. Cell Proliferation Assay*

#### Cell Culturing and MTT Assay

Unless stated otherwise, all reagents were from Merck (Milan, Italy) and of the highest purity grade commercially available. All synthesized compounds were dissolved in dimethyl sulfoxide (DMSO) and then diluted in a culture medium so that the effective DMSO concentration did not exceed 0.1% (*v*/*v*). HCT-116, MCF-7, K-562, and BALB/3-T3 cell lines were purchased from the American Type Culture Collection, Rockville, MD, USA. Except for BALB/3-T3 cells, which were grown in DMEM, all other cells were cultured in the RPMI-1640 medium. Both DMEM and RPMI-1640 were supplemented with L-glutamine (2 mM), 10% fetal bovine serum (FBS), penicillin (100 U/mL), streptomycin (100 µg/mL), and gentamicin (5 µg/mL). Cells were maintained in the log phase by seeding them twice a week at a density of 3 × 10<sup>5</sup> cells/mL, in humidified 5% CO<sup>2</sup> atmosphere at 37 ◦C.

The cytotoxic activity of the synthesized compounds against all the cell lines employed was determined by the MTT colorimetric assay as previously reported [91]. The assay is based on the reduction of 3-(4,5-dimethyl-2-thiazolyl)bromide-2,5-diphenyl-2-*H*tetrazolium to purple formazan by the mitochondrial dehydrogenases of living cells. Briefly, cells at the passage that did not exceed the number 20, were seeded into 96-well plates (Corning, New York, NY, USA) at a density of 2.0 × 10<sup>4</sup> cells/cm<sup>2</sup> , incubated overnight, and then treated with either the compounds or the vehicle (control) for 24 h. Afterward, the medium was carefully removed, and 200 µL of 5 mg/mL MTT was added. The supernatant was discarded after a 2 h of incubation at 37 ◦C, and the formazan blue formed dissolved in DMSO. The absorbance at 565 nm of the formazan product was measured using a microplate reader (LTek, INNO, Seongnam, South Korea), and the value of control cells was taken as 100% of viability. Each experiment was repeated three times in triplicate to obtain the mean values. No differences were found between cells treated with DMSO 0.1% and untreated cells in terms of cell number and viability.

The growth inhibition activity of the tested compounds was defined as the IC<sup>50</sup> value that represents the concentration of the compound that inhibits 50% of cell viability. IC<sup>50</sup> values were calculated using the dose–response inhibition model in Prism 8 (GraphPad Software, San Diego, CA, USA). The SI, a measure of the therapeutic potential of the tested compound, was calculated by dividing the IC<sup>50</sup> for normal BALB/3-T3 cells by the IC<sup>50</sup> for HCT-116, MCF-7, and K-562 cancer cells.

Two hundred milliliters (6.0 × 104 cells/cm<sup>2</sup> ) of a K-562 cell suspension were plated in each well of a 96-well plate and treated with different concentrations of compounds. An equal volume of DMSO was added to the control well, and the cells were cultured for an additional 24 h; then, 20 µL MTT (5 mg/mL) in the growth medium was added per well, and the plates were incubated at 37 ◦C for 4 h as reported in [91] with some modifications. Plates were then centrifuged at 400× *g* for 10 min. Supernatants were removed from the wells, and the reduced MTT dye in each well was solubilized in 200 µL DMSO. Absorbance was measured in a microplate reader (LTek, INNO, Seongnam, Korea), and the value of the control cells was taken as 100% of viability.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ph15010082/s1, 1H and 13C NMR spectra of new derivatives synthetized and tested in this work.

**Author Contributions:** Conceptualization, M.T.; methodology, M.T., G.C. and M.A.; software, G.C. and M.T.; synthesis G.C. and A.M.A.; biological evaluation M.A. and I.R.; data curation, M.T. and

M.A.; writing—original draft preparation, M.T.; writing—review and editing, M.T., M.A. and A.M.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available in article and Supplementary Material.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **A Review on Parallel Virtual Screening Softwares for High-Performance Computers**

**Natarajan Arul Murugan 1, \* ,† , Artur Podobas 1 , Davide Gadioli 2 , Emanuele Vitali 2 , Gianluca Palermo 2 and Stefano Markidis 1, \***


**Abstract:** Drug discovery is the most expensive, time-demanding, and challenging project in biopharmaceutical companies which aims at the identification and optimization of lead compounds from large-sized chemical libraries. The lead compounds should have high-affinity binding and specificity for a target associated with a disease, and, in addition, they should have favorable pharmacodynamic and pharmacokinetic properties (grouped as ADMET properties). Overall, drug discovery is a multivariable optimization and can be carried out in supercomputers using a reliable scoring function which is a measure of binding affinity or inhibition potential of the drug-like compound. The major problem is that the number of compounds in the chemical spaces is huge, making the computational drug discovery very demanding. However, it is cheaper and less time-consuming when compared to experimental high-throughput screening. As the problem is to find the most stable (global) minima for numerous protein–ligand complexes (on the order of 10 6 to 10 12 ), the parallel implementation of in silico virtual screening can be exploited to ensure drug discovery in affordable time. In this review, we discuss such implementations of parallelization algorithms in virtual screening programs. The nature of different scoring functions and search algorithms are discussed, together with a performance analysis of several docking softwares ported on high-performance computing architectures.

**Keywords:** computational drug discovery; virtual screening; molecular docking; chemical space; parallelization; high-performance computers and accelerators

#### **1. Introduction**

Drug discovery is one of the most highly challenging, time-consuming and expensive projects in the healthcare sector. The usual time involved in bringing a drug from basic research to market is 12–16 years, and the cost associated is about 2.5 billion dollars [1–4]. To meet one of the EU Sustainable Development Goals [5] aimed at the good health and wellbeing for everyone, drugs should be made available to common people at an affordable price, and the current protocols in drug development need to be redesigned to make the discovery process economically sustainable. One of the most promising techniques to accelerate the drug discovery process, and to make it more cost-effective, is to perform in silico virtual screening, and to exploit the computational power of large high-performance computing (HPC) systems.

One of the major contributing factors to the cost and time associated with the discovery is that it has been reported [6,7] that only one in 10,000 compounds subjected to research and development (R&D) turns out to be successful. The drug discovery involves various steps such as target discovery, lead identification, lead optimization, ADMET (absorption, distribution, metabolism, excretion, toxicity) properties optimization, and clinical trials [8].

**Citation:** Murugan, N.A.; Podobas, A.; Vitali, E.; Gadioli, D.; Palermo, G.; Markidis, S. A Review on Parallel Virtual Screening Softwares for High-Performance Computers. *Pharmaceuticals* **2022**, *15*, 63. https://doi.org/10.3390/ph15010063

Academic Editor: Osvaldo Andrade Santos-Filho

Received: 21 November 2021 Accepted: 28 December 2021 Published: 4 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Once a valid target is known for a disease, compounds from different chemical libraries are subjected to high-throughput screening against this target. If the number of compounds used for screening can be narrowed down to a few hundred, the cost and time associated with a drug discovery process can be drastically reduced. Using computational approaches, many of the steps involved in the drug discovery projects can be made to be cost-effective and less time-consuming. For example, in the case of protein tyrosine phosphatase-1B [9], the experimental high-throughput screening of a chemical library with 400,000 compounds yielded a success rate of 0.021% in identifying the ligands that can inhibit the enzyme with IC<sup>50</sup> values less than 100 µM. However, with the use of a preliminary screening phase using a computational approach, the success rate turned out to be 34.8% starting from a chemical library of 235,000 compounds.

To summarize, the experimental high-throughput screening is not suitable to deal with modern chemical spaces as they are composed of up to billions of molecules. To solve this problem, it is common to use computational approaches on HPC systems. In this review, we highlight various currently available implementations of virtual screening softwares suitable for high-performance computers. Below, we provide general introduction to virtual screening (VS) problems and discuss the possibilities for the parallelization so that it can be effectively implemented for computing facilities offered by HPCs.

The paper is organized as follows. Section 2 introduces the computational VS with details on scoring functions and search algorithms. Section 2.2 presents details on the major breakthroughs obtained in VS and Section 3 presents the main parallelization techniques used in VS and why they target HP systems. In Section 4, we provide an overview about the implementations of different VS softwares. Finally, we discuss the opportunities offered by reconfigurable architectures such as FPGAs.

#### **2. The In Silico Virtual Screening Problem**

In general, the computational approaches for molecular docking have two main components: sampling and scoring. Sampling refers to generation of various conformations and orientations for the ligand within a target binding site (defined usually by a grid box). Scoring refers to evaluating the binding/docking energies for various configurations of the ligand within the binding site. The most stable configuration of the ligand is referred to as binding pose. The VS protocol where molecular dockings are carried out for all the ligands from a chemical library includes a third component referred to as ranking, where different ligands are ranked with respect to their binding potential. Overall, the VS identifies the ligand with the topmost binding affinity (which is based on the docking energies of different ligands) for a given biomolecular target. In addition, the most stable binding mode/pose for each of the ligands within the binding site is found (which is based on the relative docking energies of different configurations of the same ligand). Figure 1 shows the general workflow of computer-aided drug discovery where the VS approach is used to identify lead compounds. It shows the steps involved in the binding pose identification of ligands within the binding site, and the ranking of different ligands is subsequently carried out to identify the lead compounds. Most of the VS schemes do not include the flexibility for the target protein, and only the sampling over translational, rotational, and torsional degrees of freedom of ligand is accounted for.

**Figure 1.** Workflow for computer-aided drug discovery where the lead compounds are identified using VS. It shows that the various configurations of the ligands within binding sites are generated using sampling, which is scored using a scoring function to identify the most stable binding mode/pose. The docking energies of the most stable configurations of all the ligands are used in ranking them to identify a list of lead compounds which are taken for further experimental validation. As the sampling and scoring need to be performed for all the ligands in the chemical library, these steps are shown within a loop.

#### *2.1. Scoring Functions*

The reliability and accuracy of the scoring functions used for screening compounds are the most important parameters that dictate the success rate of the computational screening approaches. The scoring functions are mostly defined to be proportional to the binding affinity of the ligand towards a target. The scoring functions are often classified as *physics-based*, *knowledge-based,* and *empirical*.


As discussed above, there are different scoring functions developed, and this section mainly focuses on implementations available in open-source softwares such as Dock [11], Autodock4.0 [12,13], Autodock Vina [14,15], and Gnina [16]. The docking energy defined to rank protein–ligand complexes (*s f*) in Autodock4.0 is classified as physics-based and is defined as the sum of van der Waals, electrostatic, hydrogen bonding, and desolvation energy, as shown in Equation (1).

$$sf = \mathcal{W}\_{\text{adv}} \sum\_{i,j} (\frac{A\_{ij}}{r\_{ij}^{12}} - \frac{B\_{ij}}{r\_{ij}^6}) + \mathcal{W}\_{\text{HB}} \sum\_{i,j} (\frac{\mathcal{C}\_{ij}}{r\_{ij}^{12}} - \frac{D\_{ij}}{r\_{ij}^{10}}) + \mathcal{W}\_{\text{dcc}} \sum\_{i,j} (\frac{q\_i q\_j}{\epsilon(rij) r\_{ij}}) + \mathcal{W}\_{\text{ds}} \sum\_{i,j} (\mathcal{S}\_i V\_j + \mathcal{S}\_j V\_i)(\exp\frac{-r\_{ij}^2}{2\sigma^2}) \tag{1}$$

In addition, the entropic contribution which is proportional to number of rotatable bonds is also added to the docking energy. In the equation, *rij* refers to the distance between the two atoms, *i* and *j*, centered on protein and ligand subsystems. Similarly, *q<sup>i</sup>* and *q<sup>j</sup>* refer to charges on these atoms. *Aij* and *Bij* are the coefficients of the potential energy functions describing van der Waals interaction. *Cij* and *Dij* are the coefficients of the potential energy functions describing hydrogen bonding interaction. The terms *S<sup>i</sup>* and *V<sup>i</sup>* refer to the solvation parameter and fragmental volume of atom i, respectively.

Because the protein–ligand complexes are considered to be in an aqueous environment, the binding free energies need to account for this, and solvation energy adds the binding free energy differences due to vacuum to aqueous-like environments. In particular, the last term in the equation accounts for this solvation effect (refer to Equation (1)).

In general, the entropic contributions can be due to translational, rotational, and torsional degrees of freedom. However, the docking energies implemented in the aforementioned molecular docking softwares only account for the contributions due to torsional degrees of freedom. The contributions due to other degrees of freedom are assumed to be negligible in ranking different ligands. It is worth noting that certain free energy calculation tools such as MMPBSA.py estimate the entropic contributions due to all degrees of freedom based on normal mode analysis [17]. The entropic contributions due to torsional degrees of freedoms in molecular docking softwares are oversimplified and they are estimated from the number of rotatable bonds (each bond contributes with 0.2–0.5 kcal/mol) [18]. In the case of Autodock Vina, the scoring function can be majorly classified as empirical in nature and it is a sum of various distance-dependent pairwise interactions [14]. It includes terms for describing steric, hydrophobic, and hydrogen bond interactions. The values for different parameters and weights for different terms of potential functions were obtained from nonlinear regression of the PDBbind 2007 dataset. In other words, the empirical scoring functions can have the same mathematical expression as in Equation (1), but the weights/coefficients for different types of interactions are obtained from fitting to experimental binding potentials. The knowledge-based scoring functions [19] (*s f kb*) have the following form:

$$sf\_{kb} = \sum\_{i}^{L} \sum\_{j}^{R} -k\_B T \ln[\mathcal{g}(r\_{ij})] \tag{2}$$

Here, the summation runs over the ligand and receptor atoms, and *g*(*rij*) refers to the relative probability distribution of distances of a specific types of protein–ligand atom pairs in the docked complex structure when compared to reference experimental complex structure.

Recently, deep learning networks have been proposed to provide scoring functions. For instance, Gnina uses convolutional neural network (CNN)-based scoring function to rank the protein–ligand complexes [16]. The neural networks were trained using threedimensional protein–ligand complex structures from the PDBbind database. In particular, the dataset contained two sets of poses for the ligands within the binding site. The group of positive poses had root mean square deviation (RMSD) value below 2 Å, when compared to the crystallographic poses, while the rest were considered as a group of negative poses. Here, RMSD is the root mean square deviation in atomic positions in the predicted pose when compared to reference pose, as in the crystal structure. The positive and negative poses were generated using the experimental protein–ligand three-dimensional structures by adopting a random conformation generation algorithm. The CNN model was trained using the 4D grid (which was constructed using the protein–ligand coordinates within the grid box and atom types) to classify the poses.

#### *2.2. Search Algorithms*

The search algorithms aim at finding the protein–ligand structure corresponding to the global minimum in a potential energy surface. However, this is a very challenging problem and many search algorithms end up in a local minima. Therefore, molecular docking software uses several techniques, such as deterministic search [20], genetic algorithm, Monte Carlo with simulated annealing [21], particle swarm optimization [22], or Broyden– Fletcher–Goldfarb–Shanno (BFGS) [14]. Deterministic approaches apply techniques such as gradient descent, and they focus on a reproducible sequence of pose evaluations. Monte Carlo algorithms generate numerous poses by using random values for translational, rotational, and torsional degrees of freedom of the ligand. In Monte Carlo-simulated annealing, the heating step allows the system to escape from the local minimum (during which it can sample high-energy regions of the potential energy surface). In genetic algorithms, each pose is defined by a vector of genes that correspond to translational, rotational, and torsional degrees of freedom. By varying these values in the genes, new poses can be generated. The fitness function aims at finding the minimum energy of the pose.

#### *2.3. Validation of Molecular Docking Approaches*

As discussed above, molecular docking approaches employ different types of scoring functions, and before implementation they were validated rigorously against available experimental data. In particular, two properties obtained from molecular docking can be considered in general for benchmarking:


The RMSD in the above list is computed from the experimental and predicted protein– ligand complex structures and provides an estimate of how well the molecular docking software is capable of producing the most stable binding mode and binding pose of the ligand within the target biomolecule. An RMSD value of <2 Å is considered as a threshold value for the correct prediction of complex structures [23].

A benchmark study using Autodock4.0 and Autodock Vina on the complex structures (190 in number) from PDBbind showed that the latter could predict structures within the threshold value (i.e., <2 Å) for about 78% complexes while the former one achieved 42% [14]. A recent study compared the performance of Autodock4.0 and Autodock Vina in discriminating the active compounds and decoys using a dataset of 102 protein targets, 22,432 active compounds, and 1,380,513 decoy molecules [24], The study showed that Autodock4.0 was better in discriminating actives and decoys in the targets having hydrophobic binding sites, while the Autodock Vina's performance was better for those targets having binding sites with polar and charged residues [24]. There are other studies reported in the literature which compared the performance of various molecular docking softwares such as AutoDock, DOCK, FlexX, GOLD, and ICM, and the ICM turned out to be the superior performer, with structures predicted for 93% complexes within the acceptable accuracy [23].

The other set of quantities used for benchmarking the molecular docking approaches are the inhibition constants, dissociation constants, IC<sup>50</sup> and pIC50, which are available from experimental binding assay studies. All these quantities refer to the binding potential or inhibiting potential of ligands to a specific biomolecular target. The dissociation constants and binding free energies are related to each other through the following equation:

$$
\Delta G = RT \ln \mathbf{K}\_d \tag{3}
$$

where *R* is the gas constant (equal to 1.987 cal K−<sup>1</sup> mol−<sup>1</sup> ) and *T* refers to temperature (set to 298.15 K). Thanks to this equation, the computed docking energies can be directly validated using the experimental binding assay results.

#### *2.4. Computational Cost Associated with Virtual Screening*

In a computational drug discovery project, high-affinity lead compounds against a target biomolecule (can be an enzyme, membrane protein, DNA, RNA, Quadruplex, membrane, or fibril aggregates) are identified using a reliable scoring function. One can identify lead compounds for a target from various chemical spaces. The popular chemical spaces are ZINC, Cambridge, Chemspider, ChEMBL, Pubchem, Pubmed, DrugBank, TCM, IMPPAT, and GDB13-17. Refer to Table 1 for a list of different chemical libraries with their properties [25]. The estimate for the size of chemical space for carbon-based compounds with molecular mass <500 daltons is 1060, which clearly indicates that there are limitless possibilities for designing a therapeutic compound. As can be seen, even the use of the top supercomputers with exascale computing speed (which can perform 10<sup>18</sup> floating point operations per second running the high-performance Linpack benchmark) for screening these compounds will take a universal lifetime. Therefore, it is reasonable to use certain filters to reduce the number of compounds before subjecting to screening. Recently, a filtering procedure based on Bayesian optimization algorithm, referred to as MolPAL, was used to identify the top 50% compounds by developing a model with data of explicit screening of less than 5% compounds of the chemical space [26].

Otherwise, it is reasonable to use other chemical libraries having compounds that are easy to synthesize and having favourable pharmacokinetic (ADMET) properties. The chemical library with the largest number of compounds is GDB17, which contains 166 billion organic molecules made of just 17 atoms of C, O, N, S and halogens [27]. Most of the drug discovery applications use the DrugBank database, Enamine database, ZINC15, [28], and Cambridge, and the number of compounds in these chemical libraries are listed in Table 1. As can be seen, the compounds range from tens of thousands to billions, and the computational cost associated with screening is enormous, which requires use of the HPCs and accelerators.

To demonstrate the computational demand associated with virtual screening, we describe an application below. In the case of the AmpC target, Lyu et al. docked 99 million molecules. For each compound, 4000 orientations on average were generated. Further, for each orientation, about 280 conformations were generated [29]. Therefore, for each ligand, 1.1 M docking energy calculations were carried out, and given that the number of compounds considered were 99 M, a total of 10<sup>13</sup> of such calculations were carried out. This will be further increased in the flexible receptor docking where the sampling over side-chain conformations of residues needs to be accounted for [29]. As can be seen, the computational demand is really huge with such virtual screening applications, and so it is inevitable to develop parallel algorithms and to use HPCs to accomplish such screenings within affordable time.


**Table 1.** Top chemical spaces available for VS.

#### **3. Milestones in Virtual Screening**

The first virtual screening using 3D structures of chemical compounds was carried out in 1990 against the target, dopamine D2 agonists based on which agonist with pKi 6.8 was successfully found [34]. A similar search against other targets, such as alpha-amylase, thermolysin, and HIV-I protease, yielded inhibitors with significant inhibition potential [34,35]. The top compounds obtained from screening were found to be potential inhibitors, and with further optimization of the lead compounds, the inhibition potential increased considerably. With the use of more accurate scoring functions and large-sized chemical libraries and powerful computers, the computer-aided drug discovery can become a potential route to narrow down the search space before subjecting to experimental high-throughput screening. Thanks to the currently available Petaflop/s computing facilities, one can screen a billion compounds in a day. The latest record on the high-speed virtual screening is reported by Jens Claser et al. [36], where the authors screened a billion compounds within 21 h.

The number of compounds screened against various targets keeps increasing with time. As the chance for a compound with better binding affinity increases with size of the chemical library, such screening procedures result in identifying highly potent compounds. We here list a few example cases (refer to Table 2): The screening study by J. Lyu et al. [29] yielded an active compound (339204163), which had 20 times more potency than the known inhibitors for AmpC target. The most potent inhibitor for the same target was designed by the same research group by optimizing the lead compound (the end compound was referred to as 549719643). Similarly, 10-fold more potent agonists (465129598, 270269326, and 464771011) were identified for the D4 receptor [29]. The most potent compound with 180 pM binding affinity, 621433144, for the same target has been identified by the same group. Using the multistage docking workflow, referred to as Virtual Flow for Virtual Screening (VFVS) (where different docking softwares such as Quickvina2, Smina vinardo, and Autodock Vina were used in sequence), the inhibitors for KEAP1 target were recently identified from the chemical space of 1.4 billion compounds (made of Enamine Real database and ZINC15) [37]. The first round of scoring was carried out using Quickvina2, while the rescoring was carried out for the top 3 million compounds using Smina vinardo and Autodock Vina. An inhibitor, iKEAP1 with dissociation constant 114 nM, was identified which is shown to interrupt the interaction between the KEAP1 and transcription factor NRF2 [37].

Thanks to parallel implementations of molecular docking software, not only has the number of compounds used for screening increased drastically, but the time required to accomplish the high-throughput screening has also reduced considerably. The inhibitors for the targets Purine Nucleoside Phosphorylase and Heat Shock Protein 90 were identified from the REAL enamine database (1.4 B compounds) using CPU-enabled Orion software [38]. The recent screening of compounds from Enamine database against the enzyme targets from SARS-CoV-2 was carried out in 21 h instead of 43 days with the use of parallel implementation of Autodock4.0 (referred to as Autodock-GPU) on the Summit supercomputer [39,40]. Currently, the largest experiment ever run was achieved using the EXSCALATE platform based on LIGEN software. It virtual-screened a library of 71.6 B compounds against 15 docking sites of 12 viral proteins of SARS-CoV-2. The experiment was carried out on the CINECA-Marconi100 and ENI-HPC5 supercomputers, and it lasted 60 h [41], performing more than 1 trillion docking evaluations overall. A list of megadocking and gigadocking screening calculations on large-sized chemical libraries is presented in Table 2.


**Table 2.** Topmost mega or gigadocking applications reported in the literature.

#### **4. High-Performance Computing**

Molecular docking for a single chemical compound involves sampling and scoring. This refers to sampling over the configurational phase space for the compound in the target binding site and computing the scoring functions for each of these configurations generated. In the case of virtual screening, compounds from a chemical library are ranked against a single target with respect to their binding affinities, and this additional component is referred to as ranking (refer to Figure 1). As discussed in the previous section, this is computationally very demanding and can highly benefit from the use of the parallel algorithms which can run on high-performance computers with different (shared memory and/or distributed) architectures [42,43]. Nowadays, the massively parallel computing units, referred as graphical processing units (GPUs) and field programmable gate arrays (FPGAs), are accessible to research groups or even individuals, and with the development of parallel virtual screening software, drug discovery projects can be offloaded to such groups, making the drug discovery economically sustainable. Below, we highlight the opportunities for the parallelization of VS protocol.

#### *4.1. Parallelization Strategies of Virtual Screening for High-Performance Computers*

As we discussed above, the virtual screening involves three key steps:


Each of these key steps can be parallelized, as the computations can be carried out independently.

The estimate of scoring function for each binding mode/pose within the protein binding site involves the computation of docking energies or binding free energies or any other empirical potentials. We need to find the most stable binding mode for each of the ligands in the target binding site (which corresponds to a global minimum in the protein– ligand potential energy surface). Therefore, for each ligand, millions of configurations (each configuration is a point in the ligand configurational phase space) are generated and the scoring functions are estimated for these structures. As discussed above, the configurations can be generated by changing the translational, rotational, or torsional degrees of freedom. These changes can be performed with a deterministic methodology or by using randomdriven approaches such as Monte Carlo simulations or genetic algorithms. The estimate of energies for each of these configurations can be carried out independently; these can be distributed to different computing units.

In the virtual screening, the aforementioned procedure is repeated for all the ligands in the chemical library, and even this can be carried out independently, and so can be distributed to different computing units. Finally, even the calculation of the energy of a single conformation can be parallelized, as it needs to estimate the interactions of all the ligand–protein atoms couples. Overall, the parallelization of the virtual screening approach can be implemented in the following three steps, as we also show in Figure 2.


(iii) **High-level parallelization (HLP)**: Parallelizing the ligands evaluations on different computing units.

These three different levels have different efficiencies, as the amount of data transfer between the computing units is very different and may become a strong overhead. Moreover, it is also possible to combine more of these strategies in a single application to address different level of parallelism available (such as, for example, HLP for multi-node and MLP for multi-core architectures). The low-level parallelism approach is, however, usually not appealing, since the computed docking energies for different docking poses need to be compared, and this involves frequent data transfer between different computing units. The docking energy calculation involves summation over the pairwise interactions of the atoms centered on protein and ligand subsystems. In case of flexible targets, the intramolecular energies also need to be computed, which is obtained from a double summation over the number of atoms in the target. The non-bonded interaction calculation (sum of electrostatic and van der Waals) is the most time-consuming part of the energy calculation (by 80%) [42]. Stone et al. showed a 100-times speed-up for the non-bonded energy calculation using GPUs, while Harvey et al. showed 200-times increased performance [44,45]. A more sophisticated algorithm which effectively distributed the tasks to CPUs and GPUs for computing non-bonded interactions was developed by Gine's D. Guerrero et al. [46]. In this case, each atom of the receptor was assigned to a single thread in GPUs, which handle computing its interaction with all the ligand atoms [46]. Each thread was provided with necessary ligand and receptor atom coordinates and charges. The speed-up due to this algorithm was up to 280 times.

The parallelization of scoring function calculation is the most problematic approach, as it requires to sum all the atom contributions that are evaluated by different compute units (see Figure 2). This introduces frequent small data transfers that may be limiting the scaling behavior For this reason, the other two techniques are the most used in parallelizing VS software. Indeed, computing the energy for different configurations of the ligand can be carried out independently. Similarly, finding the global minimum structure for each of the ligands as well can be carried out independently. Each molecular docking step involves generation of millions of configurations, and the scoring functions for all the conformers can be calculated independently by assigning the tasks to different computing units. The energies for conformers can be gathered and checked for the least energy configuration corresponding to global minimum in the protein–ligand free energy surface.

**Figure 2.** Various parallelization opportunities in virtual screening. The terms nlig, nconf, nlig\_atoms, and nprot\_atoms refer to number of ligands in a chemical library, number of configurations (specific to each ligand and dictated by number of rotatable bonds), number of ligand atoms, and number of protein atoms, respectively.

#### *4.2. HPCs and Accelerator Technology as Problem Solvers*

Parallel computing architectures can be classified based on the memory availability to computing units:

(i) In the shared memory architectures, a set of processors use the same memory segment.

(ii) In the case of distributed memory architectures, each computing unit has its own memory.

In general, a high-performance computing environment is usually made of both architectures where the shared-memory multi-core CPUs are connected through highperformance network connections. The memory in different nodes are local and are available to those cores within the specific node. The parallel execution of tasks in a shared memory architecture can be achieved by compiler directives such as OpenMP pragmas or libraries such as pthreads.

On the other hand, the parallel implementation of tasks in distributed memory architectures requires to handle the communication manually, using libraries such as Message Passing Interface (MPI) between different nodes of the machine.

As most of the supercomputers have both shared memory behavior (on a single node) and distributed memory behavior, as they are multi-node machines, it is easy to see the usage of both programming libraries for parallelizing the VS. In addition to multiple CPUs, recent HPCs also are integrated with accelerators such as graphical processing units, which have their own memory unit. The programming languages for CPU+GPU supercomputers are CUDA, OpenACC, and OpenCL. To use the accelerators, it is mandatory to move the data between the host and the device memory, which means that the developer needs to copy the data on the device before performing the computation and the result from the device after the computation is done. In addition, the tasks distribution to GPU cores and synchronization of the task execution in CPU and GPUs are achieved.

In virtual screening, the docking of different ligands can be carried out independently from each other. The procedure can be effectively parallelized on supercomputers with thousands of multiprocessor nodes, as well as in multi-CPU+GPU computers. As discussed above, there are different possibilities for running the virtual screening in HPCs. The molecular docking procedure for each ligand can be parallelized by distributing the docking energy calculations for different conformers of the ligand in a target binding site. In this case, the receptor coordinates and ligand coordinates and their charges need to be provided to all computing units (i.e., server nodes or GPU cards). In the second scenario, the molecular docking workload for different ligands is distributed to different computing units. In this case, the receptor coordinates and charges should be made available to all the computing units, while the ligand coordinates can be made available to the specific computing unit handling this specific ligand. The main disadvantage in this way of distributing tasks is that the ligands can have different sizes and so a rank assigned with the smaller sized ligand completes the task and waits until the molecular docking procedure is completed for all the ligands in different ranks to accept the tasks of second-round screening. This way, the computing time in this specific node is wasted. Thus, it is usually recommended to sort the ligands in terms of their size, and then the distribution of ligands can be carried out from the list. In this way, all the ranks will have more-or-less equivalent-sized tasks, and the wall time in different ranks is efficiently used. The next section provides an overview of the different software targeting HPC systems, also highlighting the type of parallelization employed.

#### **5. Current Implementation of VS Available for Workstations, Accelerators, and HPCs**

Currently, there are many parallel implementations for performing virtual screening in multicore machines, clusters, and accelerators. As discussed above, different strategies were adopted for parallel virtual screening. Only a few virtual screening softwares, such as Flexscreen, use the parallel implementation of docking energy calculation. The remaining perform the parallelization by distributing the conformational sampling/scoring or the ligands ranking segments over different computing units.

Table 3 shows a comparison of the analyzed software under some key features. As we can notice, most of the available softwares have the possibility to scale on multiple nodes, thus are able to exploit the whole HPC machine. All of them are able to fully utilize a node, while only a few have the support for the GPU acceleration. Moreover, as we already anticipated, all of them use an MLP (parallel conformational search) or HLP (parallel ligand

evaluations) or both strategies to parallelize the computation, while none of them uses a low-level parallelism (parallel energy evaluations).

We now analyze, one by one, more in detail, these softwares: we focus on their performance in massively parallel architectures and accelerators when compared to CPU implementation. Wherever possible, the accuracy of the docking results (in terms of reproducing the experimental protein–ligand complex structure and experimental inhibition constants) will be discussed.


**Table 3.** Different parallel virtual screening softwares and important features.

#### *5.1. Dock5,6*

Dock5 and 6 [47,48] are the parallel implementations of virtual screening program written in C++ with MPI libraries. They are based on the original version referred to as Dock1, which was developed by Kuntz and coworkers [11]. The Dock1 used a geometric shape-matching approach for identifying the lead compounds for a given protein target. The subsequent versions adopted physics-based scoring functions for ranking and had improvement over the thoroughness of sampling and accounted for the ligand flexibility. The recent version allowed using multiple scoring functions where the solvation energies are computed using different implicit models, such as Zou GB/SA [49], Hawkins GB/SA [50], PB/SA [51], and generalized solvent models as implemented in Amber16. The benchmarking study to evaluate the performance of recent Dock6 (V6.7) used SB2012 dataset [52] which is a collection of crystallographic structures for about 1043 receptor– ligand complexes. It was able to reproduce the crystallographic poses in this dataset to an extent of 73.3% (i.e., RMSD between the experimental and predicted binding poses was <2 Å), which is due to reduction in the sampling failure of previous Dock versions. There is also a report in the literature on the offloading of Dock6 to CPU+GPU architectures using CUDA [53]. Only the ranking using amber scoring was offloaded to GPU architectures. In this offloading, the coordinates, gradients, and velocities are copied to host (CPU) memory to device (GPU) memory and the results are copied back from GPU to CPU memory. As the GPU could only handle single-precision numbers, the original data in double-precision were converted to single-precision numbers before transfer. Overall, the study reported about a 6.5-speedup for the amber scoring in Dock6 in GPU (NVIDIA GeForce 9800 GT) when compared to AMD dual-core CPU [53].

#### *5.2. DOVIS2.0 VSDocker2.0*

Both these parallel VS softwares use Autodock4.0 as docking engine. The former one has been developed for the multi-CPU systems with Linux operating system (OS), while the latter one is for parallel computer system running with MS Windows OS. DOVIS2.0 [54] uses multithreading, SSH, and batch system for parallelization, while C++ library, Perl, and Python are used for ligand format conversion, virtual screening workflow, and receptor grid files preparation, respectively. Instead of preparing the receptor grid files for each ligand, the program reuses them, and so the I/O file operation is reduced drastically. As the program employs dynamic job controlling, the load balancing due to different size of the ligands is effectively handled. Once the docking of block of ligands assigned to a CPU is completed, the results are written to project directory. Subsequently, the CPU

requests for newer assignment for screening. This way, the CPUs are continuously in action until all the ligands are ranked and the results are written to shared directory. During a virtual screening for 2.3 M ligands from ZINC database on 256 CPUs, the DOVIS2.0 achieved a ranking speed of 670 ligands/CPU/day. The VSDocker program is developed to carry parallel virtual screening in both multiprocessor computing clusters and work stations running MS windows OS. The program is written in C++ and parallel library mpi.h [55]. In this VS implementation, the receptor map file generation and analysis of docking results for different ligands are carried out sequentially while the docking of different ligands are carried out in parallel. The performance of VSDocker has been tested using a chemical library of 86,775 ligands from ZINC database against a target made of 145 amino acids. The speed-up was comparable to that of DOVIS2.0, which is similar to VSDocker in implementation, but is suited for Linux running computer clusters.

#### *5.3. Autodock Vina*

Autodock4.0 is the most widely used molecular docking software based on physicsbased scored function, but the original version is sequential in nature. However, the Autodock Vina, which is also from the same developers at Scripps Institute, can use multiple cores simultaneously for carrying out the docking. The calculations can be executed with single-threaded or multithreading options. The implementation uses C++ with Boost thread libraries. In the multithreading version, mutiple Monte Carlo simulations are initiated with different random number seed to explore different areas of conformational space of the ligand within the binding site. In a benchmarking study, for the same protein– ligand complex, the Autodock Vina [14] with single thread ran 62 times faster than the Autodock4.0. Further, running Autodock Vina with multithreading option in eight CPU machines yielded a 7.3-times-faster (when compared to single-threaded option) completion of molecular docking. The performance of single-threaded Autodock Vina compared to Autodock should be attributed to the difference in the computational cost of scoring functions. The more than 7-times speed-up with multithreading option in an eight-core machine shows that the molecular docking calculation scales well with the number of cores. Autodock Vina relies on OpenMP for distributing tasks to different threads and so is suitable for workstations with multiple cores. However, for the supercomputers with distributed memory, this version of Autodock Vina is not suitable, and rather more robust programs that allow the data transfer and communication between different nodes need to be used. It is also worth mentioning that there are updated versions of Autodock Vina, namely, Qvina 1 and Qvina 2, which showed some speedup due to the improvement in local search algorithm. SMINA [56] is a fork of Autodock Vina with a number of additional features such as user-specified scoring function, creating grid box for docking based on the coordinates of ligand bound to target, improved minimization, feasibility to include residues for flexible docking, and possibility to print more than 20 poses. In terms of speed-up in HPCs, this did not contribute to any improvement.

#### *5.4. MPAD4*

MPAD4 [57] is a parallel implementation of Autodock4.0, and the important features when compared to parent code are listed as follows: (i) It uses MPI to distribute docking jobs across the cluster; (ii) The grid maps generated for receptor are reused for all the docking calculations while in the default version, and these files are generated for each ligand docking with the target receptor and loaded into memory and released at the completion of docking. In MPAD4, The maps are loaded into memory of the node at the beginning of tasks and are used for all remaining docking calculations. This greatly reduces I/O usage, contributing to speed-up in the performance; (iii) The OpenMP is used for the node-level parallelization in executing the LGA for finding the global minimum. Overall, it allows system-level and node-level parallelization. The performance analysis of MPI version and MPI+OpenMP version can be studied by setting OMP to 1 and 4, respectively. The computers used for the initial performance analysis were IBM Blue gene/P and sharedmemory 32-core POWER7 p755 server. The dataset used for performance analysis was

HIV protease target and 9000 compounds from a druglike subset of ZINC8 database (which has 34,481 ligands in total). The performance analysis in Blue gene/P with 2048 (8192) node (core) showed that grid map reuse has reduced single-threaded execution time by 17.5%. Multithreaded execution of the code yielded 25% improvement in the overall performance. The execution of the code on nodes ranging from 512–4096 showed near-linear scaling behavior for symmetric multiprocessing (SMP) mode (OMP = 4). In particular, for 16,384 core system, the speed gained was 92% compared to that of the ideal case. The virtual node (VN) mode (OMP = 1) with the grid map "reuse" option showed, however, 72% speed-up when compared to the ideal case. The gain in the performance of SMP mode should be attributed to multithreading. The node utilization can be further improved with the use of preordering ligands with decreasing number of torsional angles.

#### *5.5. VinaLC*

VinaLC [58] is an Autodock Vina extension to use MPI library for parallelization in large supercomputers. This implementation is very similar to VinaMPI (which is described below) and uses MPI and multithreading hybrid scheme for parallelization across nodes and within node, respectively. The computational cost of each docking calculation depends on the size of the ligand, receptors, and the grid box. If the calculations are distributed on all MPI processes due to this uneven size of the input systems, there can be MPI load imbalance where the processes are waiting for other processes to complete the task. This load imbalance is tackled effectively by the master–slave MPI scheme where the master process handles the inputs, and outputs job allocation to the slave processes. The tasks in the master slaves are handled by three loops: (i) the first loop is over each combination of receptor–target ligand which assigns docking task to a free slave; (ii) The second loop collects the docking results from slave processes and the new tasks are assigned in case of unfinished calculations; (iii) The third loop frees the slave processes. In the slave processes, an infinite whole loop is initiated which ends when the "job finished" signal is received from the master. The ideal slave processes are identified from the MPI\_ANY\_source tag and docking tasks are assigned upon the completion of previous task. In this way, the computing resources available in different processes are utilized efficiently. To make the communication effective, all the inputs needed for a single docking calculation are sent by single MPI\_Send call. Therefore, coordinates of the receptor and ligand and grids (which are computed on the fly in the master process) are sent to the slave process. The outputs from the slave processes are collected by the master process into a few files instead of generating file for each ligands (which will generate a million or billion files depending upon the size of the chemical library). The benchmarking study was carried out using the two datasets, namely ZINC and DUD (directory of useful decoys). The target was chosen as Thermus thermophilus gyrase B ATP-binding domain. The benchmarking calculations were carried out on HPC machines at Lawrence Livermore National Laboratory, and the number of cores used were in the range 600–15,408. The study showed that the average CPU time per docking was closer to ideal average CPU time. The VinaLC was shown to scale well up to 15 K cores. The percentage of I/O activity was reported to be negligible when compared to the total computing time. For aforementioned target using 15 K cores, VinaLC could screen about one million compounds from Zinc15 database in 1.4 h. This can be extrapolated to 17 million compounds per day which suggests the suitability of VinaLC for the most time-taking mega- or gigadocking screening applications.

#### *5.6. VINAMPI*

VinaMPI is another implementation of Autodock Vina for distributed computing architectures [59]. It is written in C and for communication between the nodes, it uses MPI libraries. In order to avoid poor scaling behavior of the parent-child (or masterslave) distribution scheme in massively parallel supercomputers, this implementation uses all-worker scheme. It is worth recalling that rather VinaLC used Master-slave scheme for distributing tasks. In this code, each worker (or each MPI rank) deals with its own protein-ligand complex and within each rank the computation (related to search of global

minimum) is carried out using multithreading. Due to this reason it is also suitable for the virtual screening for more than one targets. Further the computations are sorted out in terms of complexity so that the work loads at a given round of screening can be distributed equivalently. The computational complexity is measured based on the number of rotatable bonds and size of the ligands which is used to sort the tasks in the beginning of the screening. As a benchmark, the dockings were carried out for targets namely ACE (angiotensin-converting enzyme), ER AGONIST (estrogen receptor agonist), VEGFR2 (vascular endothelial growth factor receptor kinase), and PARP (poly(ADP-ribose) polymerase) with a chemical library of 98,164 compounds (comprised of ligands and decoys). Running this screening in a 516 cores supercomputer costed about 103 s [59]. It is expected that the same implementation in a supercomputer with 0.3 M cores can be used to screen 250 M compounds in 24 h.

#### *5.7. LiGen Docker-HT*

LiGen [60] is a VS software that leverages CPU and GPU to perform the required computation. Several versions of the tool have been developed, starting from a CPU-only application [20,60,61], then the main kernels have been ported to the GPU by [62] using OpenACC. Finally, it has been optimized using CUDA for the GPU kernels and this last version has been used to perform a large VS experiment in the search of a therapeutic cure for COVID-19 screening 71.6 B compounds against 15 binding sites from 12 Sars-Cov2 viral proteins on two supercomputers, accounting for 81PFlops [41]. LiGen uses deterministic algorithms to generate the different conformers of a ligand, and an empirical scoring function to select the best molecule. The Docker-HT application is the version of LiGen that is designed targeting a large VS campaign, and it is able to leverage multi-node, multi-core, and heterogeneous systems. In particular, it uses MPI to perform the multi-node communication, which is limited as much as possible by the algorithm to avoid large communication overheads [63]. Indeed, the amount of data that needs to be processed by every node is divided beforehand, and it may create load-balancing issues as there is no mechanism to rebalance it during the execution of the application. On the single node, it leverages the C++ thread library. Finally, it uses CUDA to support GPU acceleration. Within each node, the program uses pipeline parallelism and work-stealing to process the ligands.

#### *5.8. Geauxdock*

This is a parallel implementation of virtual screening available for multicore CPU, GPU, and Xeon Phi-computers. The software uses a common code for front-end computations in all these computers [64]. However, the back-end codes have one version for CPU and Xeon Phi architectures, and another for GPU. The code for CPU and Xeon Phi architectures written in C++, OpenMP and IntelSIMD pragmas. The GPU version is written in C++ and CUDA. The program uses the Monte Carlo approach for conformational search and for identifying the global minimum of the protein–ligand complex. The scoring function is based on physics-based energy terms combined with statistical and knowledge-based potentials. The performance of the code has been tested in multi-core CPUs and massively parallel architectures, namely Xeon Phi and NVIDIA GPUs. The testing using CCDC/Astex dataset showed a 1.9-times increase in performance for Xeon Phi when compared to 10-core Xeon CPU. Further, on the GeForce GTX 980 GPU accelerator, the performance was 3.5 times higher when compared to the CPU version.

#### *5.9. POAP*

Poap is a GNU parallel-based multithreaded pipeline for preprocessing ligands, for virtual screening, and for postprocessing the docking results [65]. It also allows the minimal use of memory through optimized dynamic file-handling protocol. It has also been optimized in a way that erroneous ligand input does not affect the workflow. It can be integrated with any of other sequential or multithreaded molecular docking softwares, such as Autodock4.0, Autodock Vina, or AutodockZn. In the case of Autodock-based virtual

screening, the map files are generated for each ligand in the datasets. In the case of POAP, the map files are directed to a common hub directory and so occupancy of space due to redundant atom types in ligands is overcome. The number of threads to be used should be mentioned when the Autodock Vina is used as a docking software. In the case of Autodock and AutodockZn, the number of parallel jobs to be executed (which can be equal to the number of CPU threads) should be defined by the user. The performance of POAP has been tested using the virtual screening for the targets namely Human ROCK I, HTH-type transcriptional regulator, Polyketide synthase, and PqsA (Anthranilate-coenzyme A ligase) using the ligands from the chemical library, DrugBank. Since POAP does not have any serial code, the speed-up (theoretical estimate) is directly proportional to the number of processors used. The performance analysis of Autodock in a T5510 DELL workstation with Intel Xeon(R) CPU E5-2620V2, 2.10 GHz clock speed (12 Cores, 24 threads) with 62.9 GB RAM showed a 12.4-times speed-up when compared to serial mode (the number of parallel jobs specified is 24). Similarly, the Autodock Vina showed 2.4-times speed when compared to default mode (which is already multithreaded) and here, the number of jobs was set to three.

#### *5.10. GNINA*

GNINA is a fork of SMINA [56] and Autodock Vina [16]. When compared to the hybrid scoring function employed in Autdock Vina, it provides options to use various built-in scoring functions (such as Vina, Vinardo) along with user-customized scoring functions. More importantly, it provides a machine learning-based scoring function to rank the complexes. The default scoring function (called "none" option) is the same as used in Vina or Smina, while the rescoring (called "rescoring" option) allows the topmost ligand poses to be ranked using machine learning-based scoring functions.

In particular, this scoring function is based on convolutional neural networks trained using 3D protein–ligand complex structures (as reported in PDBBind or BindingDB) and corresponding experimental inhibition constants. They were trained to reproduce binding pose and the binding affinity. There are multiple machine learning functions (namely, crossdock\_ default2018, dense, general\_default2018, redock\_default2018, and default2017) provided by GNINA, and these have been developed using different datasets. The CNN-based scoring function outperforms the scoring function implemented in Vina in reproducing the binding poses. The RMSDs for the predicted poses in the unseen examples are below 2 Å in as many cases as 56%. Further, the binding prediction of poses within this cutoff improves to 79% if the redocking is performed.

When compared to Autodock Vina, the grid box center can be provided with the help of a ligand file. For the conformational search, GNINA uses Monte Carlo sampling scheme. The sampling is carried out over the ligand translational, rotational, and torsional degrees of freedom. In the case of flexible docking, the sampling is also carried out over the residue side-chain conformations. Unlike Smina, GNINA performs computing in singleprecision (32 bit) which allows the possibility of offloading the CNN scoring tasks to GPUs. Even though GNINA can be used in massively parallel supercomputers and HPCs with accelerators, there is no performance analysis or profiling when compared to other docking softwares reported in the literature.

#### *5.11. AUTODOCK-GPU*

Autodock4.0 is one of the most widely used molecular docking softwares, but it is a serial code which runs on a single thread so cannot be effectively used in high-performance computing environments with multiple CPUs and GPUs. Autodock-GPU [66] is the version of Autodock developed for multiple node parallel computers with GPU accelerators. It is worth recalling that the above discussed MPAD4 was developed for multi-CPU architectures. This program has been developed using the application programming interface, OpenCL, as it allows portability to hybrid platforms with CPUs and GPUs. When compared to Autodock4.0, the local search algorithm uses derivatives of energies with respect to translations, rotations, and torsions (this implementation of gradient-based local

search is referred to as ADADELTA). In the case of CPU+GPU architectures, the workflow consists of a sequence of host and device functions. In analogy to biological gene, the state of the protein–ligand complex is represented by a sequence of variables. In the case of a rigid docking (where the protein framework is treated as a rigid body), the variables represent positions, orientations, and conformation of the ligand. The number of variables are 6+Nrot, where Nrot is the number of rotatable bonds. The docking is aimed at finding the genotype which corresponds to the global minimum in the protein–ligand potential energy surface and the ranking is dictated by the scoring function.

The performance of Autodock-GPU has been tested using the diverse data set of 140 protein–ligand complexes from Astex Diversity Set [67], CASF-2013 [68], and protein databank. The reference docking calculations were performed using the single-threaded Autodock4.2.6. The speedup performance was dependent on the minimization algorithm used for the local search (whether Solis–Wets or ADADELTA), GPU type, and the type of protein used in the docking. With the use of Solis–Wets local search algorithm, the speed-up was 30 to 350 times in GPUs, with the M2000 showing the least performance and with TITAN V showing the best. However, with the use of ADADELTA local search, the speedup was only 2 to 80 times improved, which has to be attributed to the computationally expensive calculation of gradients and difficulties associated with the parallelization of this local minimization step. In general, TITAN V cards showed 10 times higher speed-up when compared to M2000 versions. The performance analysis in multiple-core CPUs showed a similar trend where for the Solis–Wets search, the speed-up was in the range 5 to 33 times (the number of cores employed 8–36), while for the ADADELTA local search the speed-up was 2–20 times better.

#### *5.12. Other VS Tools*

The focus of this review was mostly about the open-source parallel VS softwares which are summarized in Table 3 along with some important features. The details about the year they were introduced and source URL pages are listed in Table 4. Many of these softwares such as GeauxDock, Autodock-GPU, and GNINA, were introduced in recent years and so their capacity in the lead compounds identification from huge chemical libraries needs to be validated extensively. Meanwhile, many already existing virtual screening softwares have contributed to successful lead compound identification and lead optimization over the years. In particular, the softwares such as FlexX, DOCK (the sequential version of above discussed Dock6), SLIDE, Fred (OpenEye), GOLD, LigandFit, PRO\_LEADS, ICM, GLIDE, LUDI, and QXP are worth mentioning [69]. Among these, LigandFit and QXP employ Monte Carlo for sampling, while SLIDE and Fred employ conformational ensembles approach. GOLD, ICM, and GLIDE, respectively, adopt genetic algorithm, pseudo-Brownian sampling/local minimization, and exhaustive search for sampling. Finally, Dock and FlexX use incremental build approach for identifying the most stable binding mode/pose for the ligand. There are other VS softwares, such as RosettaDock [70], Surflex [71], and LIDAEUS [72], which are not discussed here as the review focuses on those with parallellism capability. It is also worth mentioning Spark-VS [73] software, which uses Google's MapReduce to run parallel virtual screening in distributed cloud resources. The parallel efficiency of Spark-VS against a chemical library of 2.2 M compounds is reported to be 87% when compared to a public cloud environment [73]. This opens up another possibility of using cloud computing resources for parallel virtual screening without a need to buy our own HPCs and accelerators.


**Table 4.** Timeline for different parallel virtual screening software and source URLs (with accessed dates in bracket).

#### **6. Emerging Reconfigurable Architectures for Molecular Docking**

In prior sections, we have focused exclusively on reviewing methods of virtual screening and molecular docking that target modern central processing unit (CPU) and graphics processing unit (GPU) solutions (refer to Figure 3). At the same time, we know that Moore's law (transistor scaling) is terminating, which could motivate (or even necessitate) the search for alternative computing platforms that can continue the performance trend that molecular docking has come to rely upon. Among the many (so-called) post-Moore technologies [74], reconfigurable architectures are perhaps the most noticeable, partially because they are readily available today. A reconfigurable architecture, such as a fieldprogrammable gate array (FPGA) or coarse-grained reconfigurable array (CGRAs) [75] is a system which aspires to retain some of the silicon plasticity that is lost when manufacturing an application-specific integrated circuit (ASIC). In turn, users can leverage reconfigurable systems to perfectly match the hardware to the application, which in turn can lead to improvement in performance and reduction in energy costs. For example, the expensive von Neumann bottleneck associated with the decoding of instructions in CPUs can be virtually eliminated. Traditionally, reconfigurable architectures such as FPGAs have been programmed using complex, low-level hardware description languages (HDLs) such as VHDL or Verilog. This, in turn, has limited exposure of using these devices to specialized hardware and is thus out of reach for typical HPC users. However, with the increase in maturity of high-level synthesis (HLS) [76] tools in the past decade, today, programmers can describe their hardware in abstract languages such as C/C++ and directive-driven models (e.g., OpenCL [77] or OpenMP [78]) and automatically translate the code down to specialized hardware. Modern HLS has, in turn, facilitated the accelerated use and research of FPGAs in other HPC applications such as computational fluid dynamics, neuroscience, and molecular docking.

**Figure 3.** A conceptual picture of different processors and accelerators, showing (**a**) CPUs, which are latency-focused architectures with few processing units and large (and deep) memory hierarchies, (**b**) a GPU, which is a throughput-focused architecture with more processing units (contra CPUs) and a shallower memory hierarchy, and (**c**) FPGAs, which offer much more parallelism compared to both CPUs and GPUs, with finer control over individual unit types (here shown in different controls), but is harder to use.

Pechan et al. [79] evaluated and compared the use of FPGAs against both GPU and CPU solutions of the popular Autodock software. They created a custom RTL-based threestage FPGA accelerator that computes the performance-critical sections of the Autodock

algorithm. More specifically, the custom accelerator has four modules capable of exploiting MLP (see Section 4.1), while LLP is exploited inside each module (through pipelining); the accelerator relies on other methods to exploit HLP. They compared their solution against a custom (CUDA-based) GPU solution (GT220 and GTX260) and a CPU (Intel Xeon 3.2 GHz) version on the 1hvr and 2cpp protein pairs (from the Protein Data Bank). The overall results showed that FPGAs outperformed the CPUs for both use-cases independent of the number of dockings that were used. The GPU, however, had a clear advantage when a large number of dockings were executed, and the FPGA was only preferable when a few number of docking runs were executed. Recent work by Solis-Vasquez et al. [80,81] focused particularly on using OpenCL HLS tools to create custom accelerators that target FPGAs. Aside from disseminating their design process, they also vary several different architectural properties in their accelerator. For example, they consider both floating-point and fixed-point representation for various phases of the computation, which demonstrates an advantage that FPGAs can provide over more general-purpose approaches. The accelerator runs at a fairly high frequency (between 172 MHz and 215 MHz) on a Intel Arria 10, and consumes a varying amount of resources (subject to their design-space exploration). They compare their accelerator against the single-threaded Autodock software on five protein targets, and show that they reach between 1.73× and 2.77× speed up.

Today, there is a remarkably small number of published work that leverages FP-GAs in the Autodock software (for surveys using FPGAs on other molecular algorithms, see [82,83]). Even more surprising is that (to the authors' knowledge) CGRAs have been largely unexplored in this domain. With both FPGAs and CGRAs emerging as performance (and, more importantly, *greener*) alternatives to traditional CPUs and GPUs, we believe that these systems will come to play a much larger role in molecular docking and virtual screening in the future than they have so far.

#### **7. The Advent of Quantum Computing for Molecular Docking**

With the advent of publicly available quantum computers via cloud computing, such as the IBM, Righetti, Google, and D-Wave quantum systems [84], quantum computing is becoming a promising approach to support and accelerate molecular docking computations.

A study of a molecular docking implementation on a photonic quantum computer was presented in [85]. The authors used a Gaussian boson sampler (GBS) which is a special model of photonic quantum computer where the computation is realized via the interference of identical photons that are passing through a circuit or a network of beam splitters and phase shifters. In this work, the binding interaction graph between ligands and receptor is used to generate the ligand orientations within the protein pocket. A simplified pharmacophore representation is used, limiting the graph size from all-molecular model of the ligand and receptor to a set of points having large influence on the interactions, i.e., negative/positive-charged atoms, hydrogen bond donor/acceptor atoms, hydrophobic characteristics, and aromatic ring positions. The docking problem has been formalized by mapping it to the identification of large clusters in a weighted graph. The GBS device was used to search for the largest cliques while considering the graph weights. The method shows very good results compared to solving the same problem in a classical way; however, it cannot be used alone in a virtual screening process unless to pair it with classical data postprocessing techniques (scoring) thus generating a hybrid-quantum approach.

The usage of quantum annealers to understand the capabilities of these devices to improve the quality and the throughput of molecular docking methods is presented in [86]. In particular, the paper focused on a specific phase of the molecular docking, consisting of ligand manipulation in terms of its rotatable bonds. The authors propose a quantum annealing approach to molecular docking by formulating it as a high-order unconstrained binary optimization (HUBO), which was possible to solve on the latest D-Wave annealing hardware (2000Q and Advantage). The work demonstrated how a lot of simplifications have to be taken into consideration during the problem formulation and embedding phase, even with small molecules. The results show that despite that the current hardware is

not yet mature to solve the molecular docking problem on real-life scales, there is a clear positive trend in that direction.

#### **8. Conclusions**

The parallel implementation virtual screening algorithms in massively parallel computers with multiple CPUs and/or GPUs have the high potential to speed up the exploration of gigantic chemical spaces (having compounds in the range 10<sup>9</sup> to 1012) in real time. In a serial version of virtual screening software, it may take many years of CPU hours for such tasks. The current regard for gigantic docking is the screening of billions of compounds from ZINC15 and Enamine databases with the use of Autodock GPU in Summit HPC computers in less than a day. The parallel implementations and reliable scoring functions will increase the success rates in the lead compounds identification for drug discovery. This makes the drug discovery less time-consuming and economically sustainable. Further, as the chemical spaces are really huge, the drugs with entirely different scaffold geometry can be identified. The speed-up of the virtual screening software is found to be dependent on the number of factors: energy minimization algorithm, scoring function, biomolecular target, and computer architecture. More elaborate studies will allow us to develop highly optimized virtual screening software in the future. The implementation of VS for FPGAs and quantum computing is still in its infancy, and a dedicated research is needed for adopting such architectures for drug discovery projects.

**Author Contributions:** Conceptualization, N.A.M. and S.M.; writing—original draft preparation, N.A.M., A.P. and S.M.; writing—review and editing, N.A.M., A.P., D.G., E.V., G.P. and S.M.; funding acquisition, S.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** Funding for the work is received from the European Commission H2020 program, Grant Agreement No. 801039 (EPiGRAM-HS). We also acknowledge funding from EuroHPC-JU under grant agreement No 956137 (LIGATE) and European Commission H2020 program.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data has been presented in main text.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


*Article*
