Computational characterization of the interactions between novel Boronic Acid derivatives with urokinase-type plasminogen activator (uPA) and their binding energy

Urokinase type plasminogen activator is expected to play a significant role in metastasis therefore various inhibitors are being prepared for this target protein. However, the binding site with residues that are involved in binding and inhibition is unidentified. Hence, comprehensive computational techniques are applied for finding the binding pocket, important amino acid residues and for the characterization of the binding energy of the best ligand among seven novel boronic acid derivative inhibitors within the binding pocket. Among seven test compounds, C14H21BN2O2S showed best results in structure based molecular docking through Molecular Operating Environment (MOE) and GOLD suit with −3.2481 kcal/mol binding affinity and 46.4523 GOLD Score. C14H21BN2O2S also showed high binding affinity within the binding pocket in DFT (Density Functional Theory) studies. DFT was carried out using hybrid functional B3LYP in combination with basis set LANL2DZ level of density functional theory on the extracted geometry of bound ligand C14H21BN2O2S to the binding pocket of uPA with a −2 charge on amino-acid residue ASP189. Computational analysis values on Geometric Optimization (opt), Single Point Energy (SPE) and Self-Consistent Reaction Field (SCRF) were 53.9, −66.3 and −49.0 respectively. Hence it is concluded that C14H21BN2O2S shows better binding with uPA binding pocket when there is a negative two charge on it ASP189 amino acid residue in the binding pocket. These seven ligands were also used for generating pharmacophore model through random selection with genetic algorithm by MOE having sensitivity of 79% towards the test set, specificity of 78% towards test set and 51% calculated Matthews coefficient correlation. 2. Author Summary Boronic-acid based proteasome inhibitor like Bortezomib and Ixazomib are Food and Drug Administration (FDA) approved drugs, which are being used for fighting cancer. They can be considered as a template for understanding the pharmacokinetics and role of Boronic-acid ligands in the process. Boron-based warheads with stabilised functionality along with reduced toxicity are beneficent therapeutically. We have utilized computational quantum mechanical techniques in predicting binding free energies for ligands and proteins in a solvent environment. Instead of providing precise estimations, these techniques are more suitable for prediction purposes. The main challenge is developing inhibitors for uPA sub-sites that have high selectivity, potency, and improved pharmacokinetic properties. We have used Molecular docking and ligand-based techniques to analyze the binding interactions between seven ligands and uPA. Among these ligands, C14H21BN2O2S is identified as the most appropriate inhibitor based on scores and its interactions with specific receptor amino acid residues. Computational quantum mechanical studies are conducted using electron density and hybrid functional B3LYP to determine the binding energy. A pharmacophore model is designed to identify crucial descriptors and search for compounds that can effectively inhibit uPA. The model’s accuracy is assessed through QSAR analysis, which reveals favorable hydrogen bond donor and acceptor groups as well as aromatic hydrophobic rings in proximity to the ligands. The designed model demonstrates good sensitivity, specificity, and calculated Matthews coefficient correlation.


Introduction
Cancer metastasis is the start of last stages of tumorigenesis; it is the dissemination of diseased cells from the site of origin through detachment, followed by the movement towards other sites for invasion by blood vessels or lymphatic vessels. It is responsible for 90% mobidity of cancer patients(1), and new prognostic markers are needed to correctly predict whether a patient will develop metastases later (2). The trypsin-like serine protease (TLSP) urokinase-type plasminogen activator (uPA) may play a role in the propagation of transformed cancerous cells that cause invasive spread, metastasis, and angiogenesis in aggressive cancers like prostate, triple-negative breast, pancreatic, gastric, and colorectal (3).
ELISA and other methods show that uPA plays a major role in metastases (4). Results show high levels of uPA in tumour cells with poor patient prognosis and are seen in different aggressive cancers, making it the highest-graded biomarker for cancer metastases with the highest level of evidence (according to US National Cancer Institute scale of evidence, i.e. 1) in clinical application by performing pooled analysis with a sample of 8377 primary breast cancer patients (5). This biomarker is the best predictor of progression-free and overall survival in many tumours (6). Novel uPA inhibitors were developed due to a deep commitment to tumour suppression. After much research, antibodies and peptidomimetics inhibited metastasis and tumour growth in mice (7)(8) (9). These successes inspired researchers to find less hazardous and more effective drug-like molecules. After its role in angiogenesis and metastasis was discovered, researchers began searching for drugs that could inhibit uPA and its endogenous receptor (uPAR). Peptide aldehydes 5 were initially utilised to block cellular proteases; however, bioavailability, stability, and medication toxicity failed due to off-target toxicity (10). In 2010, Julian Adams produced peptide boronates, which were more potent against uPA (11). Boron's vacant p-orbital accepts a lone pair of oxygen on the binding pocket's serine residue to increase potency (12). These compounds are known for their strong selectivity toward serine proteases because they engage with residues in the S3 and S4 binding subsites of enzymes like elastase and chymotrypsin to produce maximum inhibitory activity. Thrombin prefers the catalytic region's P1 basic residues. Thus, leucine boronic-acid inhibitors do not interrupt action (12). First described by Xue in 2017, flavonoid inhibitors inhibit uPA molecularly. The phenolic hydroxyl ring of flavonoids binds to the S1 pocket with an IC50 of 7μM (13). Pharmacists want to highlight effective boronic acid inhibitors that generate reversible covalent connections with target proteins (11) are analysed based on their multiplicity, which is one of the confidence scores that indicate the goodness of the pocket based on the frequency with which the pocket appeared in the template structure. The first predicted pocket with the highest multiplicity 127 having the amino acid residues (H46, D192, S193, C194, Q195, G196, S198, V216, S217, W218, G219, G221, and C222) was most favourable for inhibition of uPA as suggested in the literature. was docked in uPA with ten generated conformations. The third generated conformation showed the best score of -3.3591 kcal/mol, along with the best interactions within the binding pocket. One of the OH groups showed interaction with the residue ASP189, the Sulphur present in the ligand C 8 H 11 BN 2 O 2 S showed interaction with SER214 and NH 2 group showed interaction with HIS57. C 8 H 11 BN 2 O 2 S-a (right in Figure 3) was docked in uPA, and the first conformation showed the best score of -4.4899 kcal/mol, along with the best interactions within the binding pocket. The two OH groups showed interaction with the residues HIS57 and SER214, the NH 2 group showed interaction with SER146 and GLN192, and the NH group of the ligand showed interaction with ARG217. For C₁₈H₃₁N₃O₂S (left in Figure 4), the first conformation showed the best score of -5.3983 kcal/mol. Only the NH group showed interaction with ARG217 within the binding pocket. C 9 H 9 F 3 N 2 S (right in Figure 4) was docked, and the first conformation was selected as it gave the best score of -3.7709 kcal/mol. The Sulphur group showed interaction with the residue ASP189, and the NH 2 group showed interaction with GLY219 and SER190.
10 C 14 H 21 BN 2 O 2 S (left in Figure 5) was docked in uPA with ten generated conformations. The first conformation showed the best score of -3.2481 kcal/mol, along with the best interactions within the binding pocket. The Sulphur present in the ligand C 14 H 21 BN 2 O 2 S shows interaction with GLY219, and the NH 2 group shows interaction with ASP189 and SER190. C₁₁H₁₈BNO₂S (righ in Figure 5) was docked in uPA with ten generated conformations. The first conformation showed the best score of -4.8438 kcal/mol showing no favourable interactions within the binding pocket. The Sulphur present in the ligand C₁₁H₁₈BNO₂S shows interaction with GLY219 and ARG217, and the NH group also shows interaction with ARG217.
12 Figure 6 Visual representation of ligand C 9 H 10 N 2 O 2 S docked in uPA binding pocket.  Figure 12 provides the structures of the seven test compounds along with their scores, that is, the binding free energy 13 (kcal/mol), electrostatic interaction energy (kcal/mol) and their van der wall interaction energy (kcal/mol) within the binding pocket of the protein of interest uPA.

Molecular docking studies using GOLD Suit
Molecular docking was performed using a GOLD suit to verify the results obtained from MOE, and each pose was generated using the scoring function GOLD fitness score. Ten poses were generated for each test compound. The first conformation for C 8 H 11 BN 2 O 2 S (left in Figure 7) was suggested as the best fit solution with a score of 58.3489, along with favourable interactions within the binding pocket. The NH2 group showed interaction with ASP189 and Lus224, and the Sulphur atom present in the test compound C 8 H 11 BN 2 O 2 S showed interaction with TRP215. The second conformation for C 8 H 11 BN 2 O 2 S-a (right in Figure 7) was suggested as the best fit solution with a score of 65.1328, along with favourable interactions within the binding pocket. The NH 2 group showed interaction with ASP189, ARG217, and LYS224, and the Sulphur atom present in the test compound C 8 H 11 BN 2 O 2 S-a showed interaction with TRP215.

Figure 7 Visual representation of ligand C 8 H 11 BN 2 O 2 S (left) and C 8 H 11 BN 2 O 2 S-a (right) docked in uPA binding
pocket.
The first conformation of C₁₈H₃₁N₃O₂S (left in Figure 8) was suggested as the best fit solution with a score of 47.9117, along with no favourable interactions within the binding pocket. The NH 2 group showed interaction with TRP215. C 9 H 9 F 3 N 2 S (right in Figure 8) was docked in uPA, and the first conformation showed the best fit solution with a score of 57.0528, along with the best interactions within the binding pocket. The Sulphur present in the ligand C 9 H 9 F 3 N 2 S shows interaction with ASP189 and TRP215.
15 shows interactions with ASP189 of the binding pocket. For C₁₁H₁₈BNO₂S (right in Figure 9), the first conformation showed the best fit solution with a score of 53.7381. However, the test compound does not show any favourable interactions within the binding pocket. The OH present in the ligand C₁₁H₁₈BNO₂S shows interaction with SER190.  the protein ligand complex needs to be truncated and reduced to the level of only ligand and the amino acid residues, ASP189, GLY216, GLY219, SER190, SER214, HIS57, TRP215 and LYS224 in the target protein located at the binding site and takes part in the binding interactions.
The extracted model is given below in Figure 11  For geometry optimization hybrid density functional method B3LYP was used in combination with a basis set LANL2DZ. The optimized geometry was then further utilised in the calculation of the single point energies in both gas and solvent phase using B3LYP/LANL2DZ level of DFT.

Binding pocket geometry optimization
To understand the protein binding interaction energy, the optimized geometry of protein-ligand complex were used. To obtain minimized energy for protein binding pocket without ligand, ligand structure was removed and then protein binding pocket geometry model was optimized using same level of DFT, B3LYP/LANL2DZ. Figure 13 given below depicts the selected binding pocket without the ligand.
21 Figure 13 The optimized binding pocket of uPA

Ligand Geometry optimization
For ligand model geometry optimization, ligand C 14 H 21 BN 2 O 2 S was extracted from X-ray crystal structure, PDB ID: 1W10 (16). Figure 14 shows the image of optimized ligand C 14 H 21 BN 2 O 2 S which was optimized using B3LYP/LANL2DZ level of DFT studies.

Frequency calculation for all optimized geometries
To ensure that all the model geometries are fully optimized, frequency calculations were performed on the optimized ligand, protein, and protein-ligand complex using DFT method (B3LYP/LANL2DZ). Result shows that there is no imaginary frequency present in all three model geometries and hence they are the energy minimised structures.

Single point energy calculation for all optimized geometries
Single point energies were calculated on all the optimized geometries in both gas phase as well as

Ligand Binding Affinity
The ligand binding affinity for ligand C 14 H 21 BN 2 O 2 S with protein uPA was calculated using Equation 1 given below Where ∆ is the ligand binding affinity is the energy of the ligand bound to the binding pocket is the energy calculated for the binding pocket is the energy calculated for optimized ligand All these G values are calculated in a.m.u (atomic mass unit), therefore, to convert it to kcal/mol, it needs to be multiplied with 627.51. Therefore, the equation 1 can also be written as follows: The computed relative energy value (Kcal/mol) for binding affinity of C 14 H 21 BN 2 O 2 S is presented in the Table 1 given below The results indicate that the binding interaction between C 14 H 21 BN 2 O 2 S and uPA is stable/feasible as the relative energy is exothermic with the release of -49.0 kcal/mol energy.

D. Pharmacophore modelling
The intention for building a pharmacophore model was to show a concept which explains the importance of different selected pharmacological features in the test compounds and not to optimize the inhibitory effect of these test compounds that act as inhibitors. descriptors that had no effect on the differentiation of active from inactive. The selected pharmacophoric features were revised through modifications or by altering the Gaussian radius so that the maximum number of actives could be selected as hits by the generated pharmacophoric model. The pharmacophore model for the seven ligands was 78% accurate, suggesting that the generated model in this analysis is well-predicted, making it efficient in distinguishing between active and inactive compounds, which indicates it is able to classify actives as hits selectively.

26
In conclusion, the model shown in figure 15 is finalised for the selected test set of 263 compounds that were able to select all active compounds as hits except for 13 active compounds. The pharmacophore model selected for screening the boronic acid derivative inhibitors of uPA comprises five distinct features that are: The model has 79% sensitivity to distinguishing between active and inactive test compounds.

Specificity = 78%
The model has 78% specificity towards active compounds as hits.

Discussion
Computational quantum mechanical techniques are unable to produce accurate results for binding free energies specifically for ligands and proteins in a solvent environment that mimic the real situation. Consequently, the quantum mechanical techniques would be appropriate for prediction purposes instead of the precise estimation of binding free energies for the ligands and proteins.
The development of inhibitors containing moieties which can interact with different uPA sub-sites with high selectivity, potency, and improved pharmacokinetic properties are the main challenges

29
A pharmacophore model was designed because of significant descriptors. These descriptors can be used to search for compounds that may act as efficient inhibitors that fit the model, as these features are characterised as essential regarding the biological activity of inhibitors against uPA.
The biological activity was correlated with the effect of 3-dimensional (3D) properties of ligands.

A. Prediction of the protein binding pocket
The protein data bank (PDB) provided the structure chosen for binding site prediction. The PDB contains multiple crystal structures for the same protein, and choosing the best one relies on the type of study. The criteria used to choose the crystal structure uPA(1W10) (16) was RaptorX needs the protein's FASTA sequence as an input file for the prediction of the binding pocket. With various multiplicity scores, four predicted binding pockets for the target protein uPA were predicted, and the pocket with the highest multiplicity score was selected. The amino acid residues within the predicted binding pockets were visualised using PyMOL (18). The residues were labelled as H46, D192, S193, C194, Q195, G196, S198, V216, S217, W218, G219, G221, C222 (19).

B. Molecular docking
Two software tools, molecular operating environment (MOE) and genetic optimization for ligand (GOLD) suit, were used for docking studies. The DOCK module of MOE was used to generate ten poses for each of the seven test ligands, and the steps for the ligand conformational search were placement, scoring, and refinement by energy minimization under a defined force field. The scoring function used was LondonDG (20) in combination with the placement method Alpha triangle for docking optimization protocol, and the forcefield applied for energy minimization was AMBER99 (21). Genetic optimization and ligand binding (GOLD) suit version 5.6.1 was used to generate ten poses for each of the seven test ligands. The genetic algorithm that implements quantitative prediction of binding energies was used to dock the ligand within the target protein uPA with partial flexibility (22). Ligands were docked using the scoring function GoldScore into the binding pocket of the target protein with PDB ID: 1W10 with resolution 2˚A. X, Y and Z coordinates were used for the allocation of the binding site in the docking studies within a 15˚A radius (23). The generated poses were analysed through MOE to select the most favourable conformation based on the binding interaction between the ligands and the residues of the binding pocket.

C. Quantum mechanical studies/calculations
A molecular visualisation tool, Swiss-Pdb viewer, was used to extract the binding pocket of uPA from within the whole biomolecule, which was utilized for the generation of model complex geometry (24). The Quantum mechanical (QM) calculations were performed by Density functional theory (DFT) (25) method to find the binding affinity of the test ligands within the binding pocket of the target receptor protein uPA using the commercial software GAUSSIAN 09 (26). Binding energies, single point energies, self-consistent reaction field (SCRF) and frequencies for the molecular complex were calculated, and the results generated were observed with GaussView showing the atomic density, electron density, molecular density, and molecular orbitals (29).

D. Pharmacophore model
For the target protein, uPA pharmacophore queries were generated by using MOE software. The model was built on the basis of those descriptors that could efficiently differentiate the active from the inactive. The selected pharmacophoric features were revised through modifications or by altering the Gaussian radius so that the maximum number of actives could be selected as hits by the generated pharmacophoric model (30)(31).

Acknowledgments
All the praises be to Almighty Allah the most compassionate and the most merciful, who has bestowed upon me the power and ability to think and grow, empowering me to play my role in conveying a little share of my knowledge. I will forever be grateful to my dear