Interactions of Urea-Based Inhibitors with Prostate-Specific Membrane Antigen for Boron Neutron Capture Therapy.

In this study, molecular interactions of prostate-specific membrane antigen (PSMA) with five chemically distinct urea-based boron-containing inhibitors have been investigated at the atomic level using molecular docking and molecular dynamics simulations. The PSMA-inhibitor complexations have been analyzed by comparing their binding modes, secondary structures, root-mean-square deviations, noncovalent interactions, principal components, and binding free energies. PSMA is a cell surface glycoprotein upregulated in cancerous cells and can be targeted by boron-labeled inhibitors for boron neutron capture therapy (BNCT). The effective BNCT requires the selective boron delivery to the tumor area and highly specific PSMA-mediated cellular uptake by tumor. Thus, a potent inhibitor must exhibit both high binding affinity and high boron density. The computational results suggest that the chemical nature of inhibitors affects the binding mode and their association with PSMA is primarily dominated by hydrogen bonding, salt bridge, electrostatic, and π-π interactions. The binding free energies (-28.0, -15.2, -43.9, -23.2, and -38.2 kcal/mol) calculated using λ-dynamics for all inhibitors (In1-5) predict preferential binding that is in accordance with experimental data. Among all inhibitors, In5 was found to be the best candidate for BNCT. The binding of this inhibitor to PSMA preserved its overall secondary structure. These results provide computational insights into the coordination flexibility of PSMA and its interaction with various inhibitors. They can be used for the design and synthesis of efficient BNCT agents with improved drug selectivity and high boron percentage.

Introduction

Boron neutron capture therapy (BNCT) is a noninvasive two-step cancer treatment modality utilizing the nuclear fission reaction when neutrons are captured by the boron-10 isotope (B10).[1] The first step is the selective delivery of B10-containing chemicals to the tumor cells; the second step is the irradiation of the boron-loaded tumor area with either a thermal or an epithermal neutron beam.[2,3] Irradiation with a high flux neutron beam results in the emission of lithium ions and high linear energy transfer (LET) alpha particles, which travel less than 10 μm in human tissues.[1] The ability of B10 to release high LET alpha particles at such a short distance, approximately the diameter of a cell, enables the focal deposition of ionizing radiation energy selectively to the tumor area while sparing surrounding normal tissues.[4] The boron-containing compounds for BNCT should have these two properties: (1) selective uptake by tumor cells to achieve a high concentration of absorbed neutrons in tumor. (2) Rapid clearance from circulation and normal tissues with low systemic toxicity.[5]

The main challenges for developing an ideal BNCT reagent are to achieve relatively high tumor-specific uptake (more than 20 μg/g tumor) and low normal tissue toxicity.[6] Within the short period of time between infusion and neutron irradiation, boron-containing reagents need to maintain a high concentration in the tumor area, and after that, they should be quickly cleared from normal tissue and circulation.[3,5,6] Currently, 4-borono-l-phenylalanine (BPA) and sodium borocaptate (BSH) are two clinically used BNCT reagents.[7−9] The major limitation with both BPA and BSH is the relatively low tumor-to-normal tissue ratio because of the variability in tumor uptake.[10] Due to these disadvantages, clinical usage of BNCT has been mainly limited to patients with high-grade gliomas or recurrent head and neck cancers which do not respond to conventional therapy.[11−14] For other types of cancers, BNCT is still not applied as the first-line standard of care therapy.[15] Recently, several tumor targeting strategies have been incorporated to improve the efficiency of boron delivery and retention using unnatural amino acids, peptides, proteins, sugars, and nanoparticles.[6,16,17] A couple of cell surface receptors that overexpressed on tumors have been utilized as targets for developing boron-containing ligands.[18] These receptors include the vascular endothelial growth factor receptor (VEGFR), somatostatin receptors, thymidine kinase 1, and the epidermal growth factor receptor (EGFR).[19−22] Targeted delivery of BNCT reagents such as boronated VEGFR, EGFR, and anti-EGFR monoclonal antibodies have been used in tumor vasculature and EGFR-positive glioma.[19−21]

Prostate-specific membrane antigen (PSMA), also called glutamate carboxypeptidase II (GCPII), is a 750-amino-acid-residue-containing class II transmembrane metalloenzyme, which possesses a dinuclear zinc core at the active site.[23−25] Encoded by the humans’ FOLH1 (folate hydrolase (1) gene,[26] PSMA is expressed in many tissues including the prostate, kidney, liver, intestinal epithelium, and central nervous system,[27,28] and it plays distinct biological roles in different tissues. In the intestine, the FOLH1 gene is associated with impaired dietary folate absorption, which can lead to conditions such as a low blood folate level and hyperhomocysteinemia.[29,30] In the central nervous system, N-acetyl-l-aspartyl-l-glutamate (NAAG) is hydrolyzed by PSMA into N-acetylaspartate (NAA) and glutamate (Glu), which is a neurotransmitter that relates to neurotoxicity and neuronal death in the brain at enhanced levels.[31−33] In addition, PSMA expression is upregulated in cancerous cells and used as an effective diagnostic and prognostic indicator for prostate cancer.[34] The expression of PSMA is observed both in primary and metastatic prostate cancer. Its expression levels correlate positively as prostate cancer develops into high-grade and metastatic cancer.[35] Interestingly, PSMA expression is also regulated by the androgen receptor (AR).[36] Androgen deprivation therapy such as enzalutamide can upregulate the PSMA expression level in castration-resistant prostate cancer.[37]

As an independent diagnostic/prognostic marker, PSMA ligands (especially those labeled with Gallium-68 and Fluor-18) have been extensively applied in positron emission tomography/computed tomography (PET/CT) for oligometastatic prostate cancer imaging.[38−40]

In 2021, the U.S. Food and Drug Administration (FDA) has approved 177Lu-PSMA-617, a PSMA-targeted radioligand therapy for the treatment of metastatic castration-resistant prostate cancer (mCRPC). Although therapy with 177Lu-PSMA-617 demonstrated successful tumor control, it also exhibits side effects such as xerostomia and bone marrow suppression.[41,42] Unlike 177Lu-PSMA-617 that allows PSMA-labeled radioisotopes to be concentrated within the cell, boron-10-containing compounds are nonradioactive inhibitors with low toxicity. Important organs such as the heart, kidney, liver, and lungs can be protected by avoiding dose to the tissues outside the irradiation field. Therefore, a treatment of prostate cancer with PSMA-targeting BNCT reagents may be a promising strategy for highly specific tumor delivery as it produces minimum side effects and toxicity for normal tissues. PSMA can be targeted by various molecular agents, which can be categorized into several classes: small molecules,[43] small peptides,[44] polymers,[45,46] nanoparticles,[47,48] and monoclonal antibodies.[49,50] In comparison to other reagents, small molecules and peptides can easily circulate in the blood stream through the major organs including the tumor site within a few hours because of their smaller size.[51] After 24–36 h, they can also be rapidly eliminated from the plasma by excreting into urine through the kidney. In contrast, when monoclonal antibodies against PSMA are injected, the tumor-to-normal tissue signal ratio takes around 24 h to reach its peak and can last as long as 72–120 h[52] because these antibodies possess higher molecular weight, longer stability in bloodstream, and stronger binding affinity. Therefore, smaller molecules are proposed to be the promising reagents for BNCT on account of their quicker delivery to the tumor region and shorter half-lives in the circulation system.

The X-ray structure of PSMA has been resolved at 1.84 Å resolution (PDB ID: 4NGM).[24,53] In this structure, the zinc-containing active site is located at a depth of approximately 20 angstrom (Å) from the surface.[24] The whole cavity around the active site can be separated into two substrate binding sites: S1′ site and S1 site, which accommodate the P1′ and P1 portions of inhibitors, respectively. It has been established that the S1′ site is more specific to glutamate residue or glutamate analogues, whereas the S1 site is more flexible and can accommodate different molecules.[24,25,54−57] PSMA ligands can be classified into two categories (Figure ): the first class is the phosphorus-based ligands mimicking the transition state of hydrolytic reaction and the second class is urea-based inhibitors with the hydrolysis-resistant peptide bond surrogate.[58−60] The urea-based inhibitors can bind to both S1′ and S1 sites like the NAAG, which is the natural substrate of PSMA.[61] They provide several advantages such as an ease of large-scale synthesis, penetration of the blood–brain barrier, and radiolabeling.[62,63] The unique chemical properties of these urea derivatives result in their better tumor uptake and higher binding affinity to a lipophilic pocket located near the active site of PSMA.[64] Moreover, the urea-based scaffold is also very tolerant with regard to structural modification. Therefore, these inhibitors can be exploited not only as diagnostic tools for PSMA-targeted molecular imaging but also as radiolabeled small molecules to treat advanced prostate cancer.[65]

Figure 1

Chemical structures of 2-PMPA, NAAG, urea-based inhibitor scaffold, and boron-containing inhibitors (In1–In5).

Previously reported PSMA-targeted BNCT compounds were mostly based upon the urea scaffold. Due to high boron content, a series of carborane derivatives have been synthesized and evaluated for inhibitory activity, tumor uptake, and biodistribution.[60,66,67] In this study, we aim to computationally analyze the interactions between PSMA and potential boron-containing PSMA-targeted inhibitors, which are proposed according to previous urea-based molecules and a structure–activity relationship (SAR) study.[56] In order to be clinically effective agents for BNCT, the PSMA inhibitors must possess high tumor-to-normal tissue ratios, high binding affinities, and high boron densities. A deeper understanding of interactions between the PSMA and inhibitors at the atomistic level will facilitate this goal.

Results and Discussion

The crystal structure of PSMA and its active site are displayed in Figure a. The bulk part (707 residues) of this enzyme resides in extracellular space, and the remaining part includes a short cytoplasmic N-terminal region (19 residues) and a single membrane-spanning segment (24 residues).[24,68] The extracellular part of this enzyme can be divided into three distinct domains (protease domain, apical domain, and C-terminal domain), all of which are involved in substrate recognition and binding. The protease domain (domain I, 57–116 and 352–590 residues, green color in Figure a) is mainly composed of seven central β sheets with 10 flanking α helices, whereas the apical domain (domain II, 117–351 residues, cyan color in Figure a) forms a deep funnel-shaped tunnel, from which the substrate can enter the active site (20 Å away from the protein surface). The C-terminal domain (domain III, 591–750 residues, orange color in Figure a) is constituted by four α helical bundles with an up-down-up-down order.[24,68] Under the physiological conditions, PSMA tends to exist as an active homodimer, whose dimerization interface is stabilized by the interactions between domain I/II of one monomer and domain III of another.[55]

Figure 2

(a) Crystal structure of PSMA and its binuclear active site and (b) ESP of PSMA.

Harbored in the protease domain, the active site of PSMA (Figure a) contains two zinc metal ions, designated as Zn1 and Zn2. Among them, Zn1 interacts with Asp387, Glu425, and His553, while Zn2 is ligated to His377, Asp387, and Asp453. Glu425 and Asp453 coordinate to Zn1 and Zn2, respectively, in a bidentate binding mode. Additionally, Asp387 and a hydroxyl group (μOH) bridge both metal ions. The metal–metal and metal–ligand distances are summarized in Table . The electrostatic surface potential (ESP) indicates that the funnel-shaped tunnel (black circle in Figure b) is largely composed of positive charges. The chemical structures of urea-based inhibitors are shown in Figure . The urea-based scaffold is linked to distinct boron ligands through an aromatic or an aliphatic chain. The number of boron atoms in these inhibitors gradually increases from 1 to 10 to fulfill a key requirement for effective BNCT (at least 20–50 μg of boron uptake per g of tumor).[3]

Table 1

Metal Coordination Number and Metal–Metal, Metal–Substrate, Metal–OH, and Metal–Ligand Distances (Å) for the First Coordination Shell Residues in the PSMA Crystal and Equilibrated Structures (Inhibitors 1 and 2) and the Calculated Binding Free Energies of Distinct Inhibitors (kcal/mol)a

	PSMA crystal	PSMA–In1	PSMA–In2
coordination number	5(Zn1), 4(Zn2)	6(Zn1), 6(Zn2)	6(Zn1), 6(Zn2)
MZn1–MZn2	3.26	3.45	3.55
inhibitor		2.15(Zn1), 2.19(Zn2)	2.19(Zn1)
μOH	2.08(Zn1), 1.68(Zn2)	1.91(Zn1), 1.90(Zn2)	1.88(Zn1), 1.91(Zn2)
His377	1.98(Zn2)	2.15(Zn2)	2.14(Zn2)
Asp387	1.98(Zn1), 2.03(Zn2)	2.23(Zn1), 2.08(Zn2)	2.24(Zn1), 2.12(Zn2)
Glu424	3.73(Zn2)	5.84(Zn2)	2.20(Zn2)
Glu425	2.17, 2.42(Zn1)	2.12, 2.09(Zn1)	2.21, 2.14(Zn1)
Asp453	2.02, 2.91(Zn2)	2.19, 2.11(Zn2)	2.26, 2.21(Zn2)
His553	2.04(Zn1)	2.11(Zn1)	2.21(Zn1)
binding free energy		–28.0	–15.2

Metal coordination sites (Zn1 and Zn2) are shown in the parentheses.

In this study, we investigated the interactions of PSMA with five chemically distinct urea-based boron-containing inhibitors (In1–In5) at the atomic level using molecular dynamics (MD) simulations. The structures of all complexes (PSMA–In1, PSMA–In2, PSMA–In3, PSMA–In4, and PSMA–In5) and their interactions are discussed using their binding modes, secondary structures, root-mean-square deviations (RMSD), noncovalent interactions (NCIs), principal components, and binding free energies as parameters.

PSMA–In1 Interactions

Since SAR study[56] suggested that P1′ and P1 parts of inhibitors preferred bulky groups, a phenyl moiety is introduced at both positions. As a consequence, the glutamate analogue with a phenyl ring of tyrosine occupies the P1′ segment and the lysine analogue with a phenyl ring of benzene resides in the P1 position. For Inhibitor1 (In1), a boronic acid (shown in the circle in Figure ) is attached to the phenyl ring of the P1 portion to form a phenylboronic acid. As shown in Figure a, the binding of In1 to the active site is predominately propelled by the interactions between two positively charged zinc metal ions and negatively charged oxygen atoms of the inhibitor. Specifically, the carbonyl oxygen of ureido linkage coordinates to Zn1 at a distance of 2.15 Å and one oxygen of the carboxylate group from P1 connects to Zn2 with a distance of 2.19 Å, resulting in a μ-1,6 bidentate mode. The bridging hydroxide μ-OH shifts to a more symmetrical position (1.91 Å (Zn1) and 1.90 Å (Zn2), Table ) and the metal–metal distance is elongated by 0.19 Å upon inhibitor binding, which is consistent with experimental observations.[24,68] All ligands surrounding Zn2+ ions (six-coordinated) adopt a distorted octahedral geometry. The boronic acid that links to the aromatic moiety of P1 points to the opposite direction of the funnel-shaped tunnel, protruding into a pocket enclosed by Arg210 and Tyr234. The RMSD value derived from the superposition (Figure b) of equilibrated and crystallographic structures of the active site is only 0.61 Å, substantiating the validity of our MD simulations. Furthermore, the RMSD plot (Figure a) for the α carbon of the enzyme indicates that the simulation of the PSMA–In1 complex achieves equilibrium at around 40 ns. This equilibrium is further confirmed by the pairwise RMSD plot (Figure S1a) that calculates the RMSD value of each frame to all other frames in the trajectory. Dominating fluctuations focus on the S1′ site and C-terminal loop, whereas the S1 site and the first coordination shell experience less fluctuations (Figure b). It also exhibits the rigidity of the S1′ site and flexibility of the S1 site. The PCA of the PSMA–In1 complex (Figure c) displays a V-shaped graph accompanied with a variety of narrow spread energy basins, manifesting the existence of multiple ensemble conformations. A similar PCA graph is also observed in the MD simulation of glycerolphosphodiesterase (GpdQ) combined with the paraoxon substrate.[69,70]

Figure 3

Most representative structure of the PSMA–In1 complex derived from MD simulations: (a) active site with the inhibitor, (b) superposition of the active site of the equilibrated structure (magenta carbon) and crystal structure (cyan carbon), (c) principal component analysis, (d) NCI plot of the S1′ site with amino acid residues labeled, (e) NCI plot of the S1 site with amino acid residues labeled, and (f) contact map of the S1 site.

Figure 4

(a) RMSD of all complexes and (b) root-mean-square fluctuations of all complexes.

The glutamate analogue of In1 occupies the S1′ pocket, which is composed of 12 amino acid residues. These residues are highlighted as the primary source of interaction by the NCI contour plot (Figures d and S2a). The α-carboxylic group of In1 at the S1′ site forms a salt bridge and a hydrogen bond with the guanidinium moiety (Nη) of Arg210 and the backbone amino group of Asn257 at distances of 2.89 and 1.99 Å, respectively. The aromatic ring interacts with the γ-amino group (Nδ) of Asn257 and the side chain of Phe209 through NH−π and T-shaped π–π interactions, respectively. Overall, the substrate binding mode of In1 in the S1′ site is analogous to those of the PSMA complex with GPI-18431[24] and 2-PMPA,[68] both of which possess the exact glutamate in their structures. Typically, Lys699 and Tyr700 are denominated as the “glutarate sensor” to recognize two carboxylate groups of glutarate in the substrate via the induced-fit mechanism.[24,68] When there are no substrates, these two residues are located outside the S1′ site, while in the presence of the substrate, they enter the S1′ site and principally select glutarate-like structures. However, in the equilibrated structure of the PSMA–In1 complex, although Lys699 is confirmed to clamp the terminal carboxylate group, the hydroxyl group of Tyr700 is approximately 5 Å away from the α-carboxylic group. The reason for this discrepancy could be that the strong π–π interaction between Tyr700 and the aromatic ring of the P1 segment restricts the movement of this residue to a great extent.

The S1 pocket of PSMA exhibits the capability to accommodate various binding modes on the grounds of its high flexibility and accessory binding tunnel. The NCI contour plot of PSMA–In1 in the S1 site (Figure e) illustrates that apart from coordinating to Zn2, the α-carboxylate group of In1 forms a strong salt bridge (2.31 Å) with Arg536 and receives a strong hydrogen bond (1.70 Å) contributed by Tyr552. The phenyl ring of the P1 portion is sandwiched by two phenol rings of Tyr552 and Tyr700 through parallel-displaced π–π stacking. The side chains of Arg210 and Tyr234 function like two gates to encompass the hydroxyl groups of boronic acid, exploiting the electrostatic and OH−π interactions, respectively. Previous studies proposed that the aromatic ring of some inhibitors such as the PEGylated inhibitor[55] and DCIBzL[61] was trapped within an arginine patch (arginines 463, 534, and 536) and stabilized by cation−π interactions. Normally, when PSMA was bound with natural substrate NAAG or in its resting state, the side chain of Arg536 was either in “stacking” or “binding” conformation, accompanied by the “up” or “down” position of the side chain of Arg463. The “binding” conformation of Arg536 was stabilized by the hydrogen bonds with the α-carboxylate group of P1. The opening of this patch or the so-called S1-accessory hydrophobic pocket was only associated with an abnormal combination of Arg463’s “up” position and Arg536’s “binding” conformation.[54,61] Although the full insertion of the terminal phenyl ring into the pocket led to the enhanced inhibition constant, this phenomenon was only observed for limited inhibitors with an appropriate length of the P1 chain and suitable substituent ligands (one halide atom on the para position of the phenyl ring). Nevertheless, it served as an important way to design and synthesize PSMA inhibitors with high binding affinity. However, in our simulation, even though Arg463 and Arg536 stay in an open state, In1 rejects to enter this arginine patch. This is most likely due to the steric clash produced by the comparatively large boronic acid ligand with bifurcated hydroxyl groups. Instead, it is captured by another hydrophobic pocket (Figure S3a) surrounded by Arg210, Tyr234, Phe546, Tyr552, and Tyr700. Their intimate distances are also revealed by the contact map in Figure f. Therefore, the loss of binding affinity for not locating in the arginine patch is greatly compensated by the insertion of the P1 portion into another hydrophobic pitfall. The binding free energy of In1 to PSMA is calculated to be −28.0 kcal/mol (Table ), utilizing the λ particle approach in a thermodynamic cycle. This energy is mainly contributed by the following three interactions: (1) interactions between Zn ions and oxygen atoms of the inhibitor; (2) interactions between the S1′ site and the P1′ segment; 3. interactions between S1, especially the hydrophobic pocket, and the P1 region. Experimental results[67] on a structurally similar compound to In1 demonstrated a half maximal inhibitory concentration (IC50) value of 130.3 nM, which was higher than that of 2-phosphonomethyl pentanedioic acid (strongest PSMA inhibitor ever reported, 2-PMPA, ∼0.3 nM)[71] but lower than that of glutamate (∼30 μM).[72] Although the IC50 value is not a direct indicator of binding affinity, the two can be related to each other using the Cheng–Prusoff equation.[73] Essentially, the lower the IC50 value, the higher the binding affinity and thus the more negative the binding free energy.

PSMA–In2 Interactions

A boronic ester—pinacolborane (4,4,5,5-tetramethyl-1,3,2-dioxaborolane shown in the circle in Figure ) is attached to In2 to investigate the influence of the ester group formed by bifurcated methyl groups. In contrast to In1, In2 binds in the monodentate manner to the active site of PSMA (Figure a) that is characterized by the following factors: the ureido oxygen is engaged in the coordination of Zn1 (2.19 Å) and the P1 carboxylate group drifts away from the second metal Zn2 (4.67 Å). As a result, Zn2 coordinates to a previous second coordination shell residue—Glu424 (2.20 Å) to maintain its six coordination number and distorted octahedral geometry. Except for His377, all first coordination shell residues interact with both metal cations at a relatively longer distance in comparison to PSMA–In1 (Table ). On account of the large size of pinacolborane, the phenylboronic ester escapes from the hydrophobic site that wraps around In1 and intercalates into the funnel-shaped tunnel that has enough space to accommodate large ligands. The equilibrated structure of the active site deviates from the crystal structure by an average of 1.41 Å (Figure b), winding up with significant displacement of nearly every residue and metal ion. Akin to PSMA–In1, the simulation of the PSMA–In2 complex attains equilibrium at approximately 40 ns (Figures a and S1b). Upon In2 binding, the population of conformational ensembles in PSMA–In2 is quite comparable to that in PSMA–In1 because the corresponding PCA (Figure c) reveals an R-shaped graph with a similar amount of energy basins and whole surface area.

Figure 5

Most representative structure of the PSMA–In2 complex derived from MD simulations: (a) active site with the inhibitor, (b) superposition of the active site of the equilibrated structure (magenta carbon) and crystal structure (cyan carbon), (c) principal component analysis, (d) NCI plot of the S1′ site with amino acid residues labeled, (e) NCI plot of the S1 site with amino acid residues labeled, and (f) contact map of the S1 site.

The outline of the tunnel that incorporates the P1 part of In2 is clearly observed in the NCI contour plot of the PSMA–In2 complex (Figure S2b). In comparison to the PSMA–In1 complex, the composition of the S1′ site is analogical but still composed of several different residues. The phenyl ring forms tilted T-shaped π–π interaction with the side chain of Phe209. In contrast to In1, the terminal carboxylate group of In2 forms only two hydrogen bonds (2.34 Å and 1.81 Å, respectively) with Leu261 and Asn262. Besides, this group and adjacent ester oxygen are recognized by the ε-amine group of Lys699 (Nζ) through a weaker salt bridge (3.30 Å) and a weaker hydrogen bond (2.76 Å), respectively. All these differences may correlate with a slight rotation of P1′ in PSMA–In2 induced by the position variance of P1. Analogous to In1, the strong π–π interaction with the aromatic ring of P1 precludes Tyr700 to function as a “glutarate probe sensor”. Alternatively, Asn257 takes the place of Tyr700 to interact with the α-carboxylate group of P1′, and Tyr700’s contribution is largely compensated by π interactions around the aromatic ring of P1′.

Depicted in Figure e, the NCIs in the S1 site of the PSMA–In2 complex are substantially weaker in comparison to PSMA–In1. The phenylboronic ligand breaks away from the clamp of Tyr552 and Tyr700, shifting to the accessory binding tunnel (Figure S3b) that contains sufficient space for its accommodation. In this tunnel, the side chain of Arg463 is included in a cation−π interaction with the phenyl ring of P1 and a hydrogen bonding with one ester oxygen atom. Parallel-displaced π–π stacking manifests itself between the phenyl ring and Tyr700. These interacting residues are clearly displayed as orange square points at the third row of the S1 site’s contact map (Figure f). Since the size of boronic pinacol ester is larger than that of boronic acid, the possibility of its location at the preceding arginine patch can be completely obviated. The association of In2 with PSMA demonstrates a subtly increased (1.7%) α-helix content (39.1%) but the overall similar secondary structure (Figure S4) as the PSMA–In1 complex. The computed binding free energy for In2 is −15.2 kcal/mol (Table ), which is substantially higher than the one (−28.0 kcal/mol) calculated for In1. The decreased binding affinity of In2 is substantiated by the increased IC50 value (318.4 nM) of an In2 analogue[67] and is related to the conversion from the bidentate binding mode to the monodentate mode. The weaker NCIs in S1′ and S1 sites also contribute to this reduction.

PSMA–In3 Interactions

As discussed above, in order to invoke the successful BNCT, at least ∼20 μg of B10 per weight of tumor has to be selectively delivered into the tumor cells. Once these boron atoms absorb enough thermal neutrons, they will produce high LET α particles with pathlengths of 5–9 μm in tissues to destroy the boron-containing tumor cells while sparing the adjacent normal cells.[1,3] To achieve that, boron density should be enhanced by increasing the percentage of boron atoms within an inhibitor. In In3, the phenyl ring of P1 is substituted by a borazine ring, which is a cyclic compound constructed by three BH units and three NH units alternatively. Although isoelectronic and isostructural with benzene, borazine is suggested to be aromatic but possesses less electron delocalization than all-carbon analogues due to different electronegativities of the boron and nitrogen atoms.[74,75]In3 coordinates to the bimetallic center of PSMA in a μ-1,6 bidentate fashion (Figure a) analogous to In1. The relative positions of the first coordination shell residues are preserved (Table ) so that both Zn2+ ions are still six coordinated and adopt the distorted octahedral geometry. Based on the structural similarity, the boron ligand (borazine ring and boronic acid) of In3 inserts into the same hydrophobic pocket that accommodates the phenylboronic acid of In1. The first and second principal components of α carbon atoms construct an inverse V-shaped graph in the PCA of the PSMA–In3 complex (Figure c). In contrast to PSMA–In1 and PSMA–In2 complexes, fewer energy basins in this graph suggest that the In3 binding triggers less structural variations in PSMA.

Figure 6

Most representative structure of the PSMA–In3 complex derived from MD simulations: (a) active site with the inhibitor, (b) superposition of the active site of the equilibrated structure (magenta carbon) and crystal structure (cyan carbon), (c) principal component analysis, (d) NCI plot of the S1′ site with amino acid residues labeled, (e) NCI plot of the S1 site with amino acid residues labeled, and (f) contact map of the S1 site.

Table 2

Metal Coordination Number and Metal–Metal, Metal–Substrate, Metal–OH, and Metal–Ligand Distances (Å) for the First Coordination Shell Residues in PSMA Equilibrated Structures (Inhibitors 3, 4, and 5) and the Calculated Binding Free Energies of Distinct Inhibitors (kcal/mol)a

	PSMA–In3	PSMA–In4	PSMA–In5
coordination number	6(Zn1), 6(Zn2)	7(Zn1), 6(Zn2)	7(Zn1), 6(Zn2)
MZn1–MZn2	3.48	3.40	3.35
inhibitor	2.24(Zn1), 2.24(Zn2)	2.18,2.43(Zn1), 2.32(Zn2)	2.20, 2.50(Zn1), 2.36(Zn2)
μOH	1.89(Zn1), 1.85(Zn2)	1.85(Zn1), 1.87(Zn2)	1.87(Zn1), 1.99(Zn2)
His377	2.15(Zn2)	2.07(Zn2)	1.96(Zn2)
Asp387	2.26(Zn1), 2.11(Zn2)	2.09(Zn1), 2.12(Zn2)	2.39(Zn1), 2.05(Zn2)
Glu424	8.52(Zn2)	6.00(Zn2)	4.48(Zn2)
Glu425	2.25, 2.11(Zn1)	2.27, 2.11(Zn1)	2.20, 2.27(Zn1)
Asp453	2.11, 2.17(Zn2)	2.27, 2.17(Zn2)	2.12, 2.24(Zn2)
His553	2.15(Zn1)	2.05(Zn1)	2.02(Zn1)
binding free energy	–43.9	–23.2	–38.2

Metal coordination sites (Zn1 and Zn2) are shown in the parentheses.

The overall microenvironment of In3 inside PMSA resembles that of In1 (Figure S2c). However, in comparison to PSMA–In1, Ser517 and Gly518 interchange their positions with Gly256 and Glu522 at the S1′ site of PSMA–In3 (Figure d). The α-carboxylate group of P1′ in In3 is recognized by the guanidinium moiety (Nη) of Arg210 and the backbone amino group of Asn257 through a salt bridge (2.43 Å) and a hydrogen bond (1.99 Å), respectively. The borazine ring of P1 approaches the side chains of Phe209 and Arg210, suggesting that In3 intercalates more deeply into the S1′ site than In1. This deeper intercalation can also be confirmed by the shorter hydrogen bonds (1.74 and 1.92 Å) between the terminal carboxylate of P1′ and Leu261/Asn262. As a result, the aromatic ring in P1′ forms cation−π interactions with the ε-amine group of Lys699, while in the previous two cases, Lys699 interacts with the terminal carboxylate group through a salt bridge. Additionally, the phenyl ring of Phe209 further stabilizes the same aromatic ring by tilted T-shaped π–π interactions. Similar to the complexes of In1 and In2, Tyr552 and Tyr700 are restrained by the borazine ring of In3, preventing them from donating hydrogen bonds to the α-carboxylate group. On the other hand, the guanidinium group (Nη) of Arg536 in the S1 site (Figure e) interacts with the α-carboxylic group of P1 by a strong salt bridge (2.16 Å), thus staying in the “binding” conformation. Arg463 can play the role of position mark, and the absence of its interaction with In3 indicates that the inhibitor is engaged in a translational motion toward the S1′ site and the boron ligand is more buried in the hydrophobic pocket (Figure S3c). In the PSMA–In1 complex, the aromatic ring of P1 interacts with Tyr552 and Tyr700, while Arg210 and Tyr234 only form interactions with the hydroxyl groups of boronic acid. In contrast, albeit with a slightly weaker aromaticity, the increased localization of electrons on the BH and NH units of the borazine ring can interact with Arg210 and Tyr234. Moreover, it turns out that all these four residues and Phe209 besiege the borazine ring (close distances shown in Figure f). The phenyl ring of Phe209 and the phenol rings of Tyr552/Tyr700 form π–π interactions with the borazine ring, which in turn associates with Tyr234 through OH−π interaction. A binding free energy of −43.9 kcal/mol for In3 is significantly lower than those for In1 and In2 (Table ). This difference suggests that PSMA prefers to associate with In3. The stronger binding of In3 than of In2 is understandable (bidentate vs monodentate binding). However, its stronger binding in comparison to In1 is trivial. From the viewpoint of NCIs, the slightly less aromatic but more electrostatic borazine ring of In3 is more buried in the hydrophobic pocket, producing more interactions around the boron ligand and pushing P1′ more deeply into the S1′ site. This movement again generates more π–π interactions around the phenyl ring of P1′ and stronger hydrogen bonds on the terminal carboxylate group.

PSMA–In4 Interactions

In4 consists of the same number of boron atoms as In3, but it contains a different aromatic ring—boroxine ring, which is composed of alternating oxygen and singly hydrogenated boron atoms (shown in the circle in Figure ). Thus, the boronic acid can only occupy the meta rather than the para position with respect to the P1 chain. The boroxine is also isoelectronic to benzene, yet it possesses less aromatic character than borazine due to the larger difference in electronegativity between the oxygen and boron atoms.[76] With vacant p-orbitals on boron atoms, boroxine can accept electrons from phenyl or negatively charged groups. In the active site of PSMA–In4 (Figure a), all three oxygens from the ureido linkage and carboxylate side chain of P1 are involved in a tripartite ligation of the bimetallic cations. This appealing observation seems to be a bit gratuitous at first glance but maybe correlated with the movement of the P1 part in the S1 site. Although the coordination between the carboxylate oxygen and Zn1 is loose (2.43 Å, Table ), it alters the molecular geometry at the Zn1 site from the distorted octahedron to the distorted pentagonal bipyramid. First coordination shell residues experienced some fluctuations (Figure b) and display a slightly higher RMSD value (0.76 Å, Figure b) than those for In1 and In3 complexes. The PCA of the PSMA–In4 complex (Figure c) offers an inverse V-shaped graph that has four distinct energy basins. Since PCA of the In3 complex has a similar number of energy basins, these two inhibitors are considered to generate comparable conformations during the simulations.

Figure 7

Most representative structure of the PSMA–In4 complex derived from MD simulations: (a) active site with the inhibitor, (b) superposition of the active site of the equilibrated structure (magenta carbon) and crystal structure (cyan carbon), (c) principal component analysis, (d) NCI plot of the S1′ site with amino acid residues labeled, (e) NCI plot of the S1 site with amino acid residues labeled, and (f) contact map of the S1 site.

As expected, the meta substitution of boronic acid on the boroxine ring causes it to jump out of the hydrophobic pocket that can occupy both In1 and In3. The drift of the boron ligand of In4 is deemed to result in the tripartite coordination to the bimetallic site and the shift-up of the P1′ part, which promotes the loss of contact between P1′ and Phe209/Lys699 (Figures d and S2d). Additionally, the salt bridge between Arg210 and the α-carboxylate side chain of P1′ becomes very weak (3.04 Å), while the hydrogen bond (2.28 Å) donated from Asn257 is conserved. The residues of Asn257, Gly427, Leu428, and Asn519 contribute to the stabilization of the phenyl ring of P1′from both the sides. Furthermore, the terminal carboxylate group is trapped by a hydrogen bonding network created by residues Leu261, Asn262, Ala264, and Gly265. In the S1 site, the boron ligand resides itself in the accessory binding tunnel (Figure S3d) and less NCIs are observed between In4 and environmental residues (Figure e,f). The α-carboxylate group of P1 makes a very weak salt bridge (4.00 Å) with the guanidinium moiety (Nη) of Arg536. Differing from previously inhibitors, Arg534 also exhibits some weak electrostatic interactions with the amide bond of P1. A nearly perfect cation−π interaction manifests itself between the boroxine ring and the guanidinium group of Arg463. It is noteworthy that the boroxine ring is completely devoid of π–π interaction with Tyr700 due to its’ weaker aromaticity. In general, tripartite ligation of the inhibitor to metal ions is suggestive of a high binding potency. Nevertheless, the binding free energy of In4 is computed to be −23.2 kcal/mol, which is 4.8 and 20.7 kcal/mol higher than those of In1 and In3, respectively. This lower binding affinity could be explained by less NCIs around the boron ligand, weaker salt bridges on α-carboxylate groups, and loss of contact with Phe209/Lys699 described previously.

PSMA–In5 Interactions

In order to further increase the boron density within PSMA, a closo carborane ball (1,2-C2B10H11–1, shown in the circle in Figure ) is attached to the phenyl ring of P1 to create In5. 1,2-C2B10H11–1 is an icosahedral electron-delocalized nonclassical bonding cluster composed of hydrogen, boron, and carbon atoms.[77] It is usually characterized by the electron-deficient 3-center 2-electron (3c-2e) chemical bond, and its stability is derived from completely filled bonding molecular orbitals. The carborane clusters and related metallacarboranes are employed in a wide range of applications including catalysis, medicines, electroactive materials, and recovery of radioactive metals from nuclear waste and heat-resistant polymers. The electronic structures of these electron-delocalized polyhedral compounds can be predicted by the Wade-Mingos rules[78,79] that are invoked depending on the number of electrons per vertex. The equilibrated structure of the PSMA–In5 active site (Figure a) displays a tripartite coordination similar to that of the PSMA–In4 complex. The comparatively large-sized carborane ball is completely accommodated in the accessory binding tunnel like the cases of In2 and In4. This is in accordance with the crystal structure of PSMA in complex with a carborane-containing inhibitor.[60] Large deviations in the equilibrated active site of PSMA–In5 are reflected by the high RMSD value (1.38 Å, Figure b). When compared with In1–In4, the conformational variations upon binding of In5 are the largest because the PCA of the PSMA–In5 complex (Figure c) shows a unique U-shaped plot with maximum conformational ensembles. It illustrates that the accommodation of In5 inside PSMA induces maximum structural variations in contrast to all previous inhibitors.

Figure 8

Most representative structure of the PSMA–In5 complex derived from MD simulations: (a) active site with the inhibitor, (b) superposition of the active site of the equilibrated structure (magenta carbon) and crystal structure (cyan carbon), (c) principal component analysis, (d) NCI plot of the S1′ site with amino acid residues labeled, (e) NCI plot of the S1 site with amino acid residues labeled, and (f) contact map of the S1 site.

Just like the In2 complex, the shape of the accessory binding tunnel is clearly shown in the NCI plot of the whole cavity in PSMA–In5 (Figure S2e). Due to similar positions of the boron ligand and identical tripartite ligation, the P1′ part of In5 resembles that of In4 in terms of the NCIs in the S1′ site. The absence of salt bridges contributed by Arg210 and Lys699 could be ascribed to the upward movement of the P1′ portion. The β-amide group of Asn257 and the peptide bond of Gly427/Leu428 are involved in parallel and tilted T-shaped π stacking with the aromatic ring of P1′, respectively. The methylene/methyl groups of Leu428 engaged in CH−π interactions with the same ring, which adopts a tilted T-shaped arrangement relative to the side chain of Phe209. It is worth mentioning that aromatic-amide interaction is one of the most common NCIs in proteins. However, the NH−π contacts are generally outnumbered by the aromatic-amide-stacked structures with the sp2-hybridized nitrogen atoms.[80] This phenomenon can be attributed to the higher number of conventional hydrogen bonds the NH group can form in the stacked orientation.[81] Analogously, the terminal carboxylate side chain interacts with the backbone amino groups of Leu261, Asn262, Gly263, and Ala264 via hydrogen bonds. In the accessory binding tunnel (Figure S3e), the P1 of In5 displays more interactions with residues that wrap around it (Figure e) in comparison to In2 and In4. Two strong hydrogen bonds (2.07 and 2.32 Å) are formed between the carbonyl oxygen of amide in P1 and guanidinium protons (Hθ) of Arg534. The aromatic ring of P1 is arranged in a perfect parallel-displaced π stacking with the phenol ring of Tyr700 on one side. On the other hand, the backbone carbonyl bond, α-methylene group, and phenyl ring of Phe546 add to the stability of this aromatic ring. Being electron-deficient, the carborane ball is predisposed to receive electrons from electron-donating groups to enhance the delocalization effect. The methylene groups of Lys207, guanidinium group of Arg463, peptide bond of Lys545/Phe546, γ-amino group of Asn698, backbone carbonyl oxygen of Tyr700, and α-methyl group of Ala701 are all found to donate electrons to carborane. All these interactions and tripartite coordination are confirmed in the contact map (Figure f), thus predicting a high binding potency for In5. Its binding free energy is calculated to be −38.2 kcal/mol, which is the second lowest among all inhibitors studied. This may result from the lack of three salt bridges that originated from Arg210, Lys699, and Arg536 in the PSMA–In5 complex. The IC50 value of an In5 mimic was measured to be 20.3 nM.[67] Moreover, a series of closo-, nido-, and iodo-C-hydroxy carborane clusters have been reported to bind with PSMA for BNCT.[66] Among them, a radiolabeled 123I analogue exhibited the highest affinity with an IC50 value of 73.2 nM. This indicated that except for the urea-based scaffold, the carborane cluster is also tolerant with respect to structural modification. In this situation, it is beneficial for us to search for the best BNCT compounds by modifying the existing molecules.

Summary and Conclusions

In this study, the molecular docking and MD simulation techniques have been employed to investigate the binding potencies and interactions of five potential urea-based boron-containing inhibitors with PSMA. These inhibitors are derived from the common urea-based PSMA binding skeleton that has been reported previously. They contain chemically distinct boron ligands (Figure ): (1) phenylboronic acid in In1; (2) phenylboronic ester in In2; (3) borazine ring and boronic acid in In3; (4) boroxine ring and boronic acid in In4; and (5) carborane ball in In5. The number of boron atoms within one molecule is gradually increased from 1 to 10 because the effective BNCT not only requires the strong binding affinity of the inhibitor to the target protein but also necessitates the high boron density within tumor cells.[3,5] It is found that different inhibitors bind in a rather distinct mode to generate PSMA–inhibitor complexes, which are mainly dominated by hydrogen bonding, salt bridge, electrostatic, and π–π interactions. Specifically, In1 and In3 bind to the bimetallic active site in a μ-1,6 bidentate mode. Their boron ligands are inserted into a hydrophobic pocket encompassed by Arg210, Tyr234, Tyr552, and Tyr700 (Figure S3). In contrast, In2 is singly coordinated to Zn1, while In4 and In5 are involved in a tripartite ligation to both metal ions. The boron ligands of these three inhibitors are accommodated in the funnel-shaped tunnel due to different reasons: the comparatively large size of boron ligands in In2 and In5 and the meta substitution of boronic acid on the boroxine ring in In4. Nonetheless, they all share a commonality: none of the ligands are trapped in the arginine patch mentioned in the previous studies.[55,61] This may be attributed to the steric hindrance created by those large boron ligands. By and large, these results demonstrate the flexibility of the S1 site of PSMA because the P1 portion of an inhibitor can be located in various positions based on its chemical nature. On the other hand, the P1′ part of all inhibitors is located in the S1′ site of PSMA and interactions around this part are influenced by the position of the P1 portion. For In3, the P1′ intercalates deeper into the S1′ site, exhibiting more interactions with PSMA. However, for In4 and In5, the P1′ shifts up and loses contact with several important residues (Phe209 and Lys699) at the bottom of the S1′ site.

The association of the inhibitor to PSMA does not alter the overall secondary structure of PSMA (Figures S4 and S5), which consists of approximately 38% helical and 15% beta-sheet conformations. The binding free energies are calculated using the λ-particle approach in a thermodynamic cycle and are analyzed based on the NCIs between the inhibitor and PSMA. The computed free energies suggest a binding preferential order (Tables and 2): In3 (−43.9 kcal/mol) > In5 (−38.2 kcal/mol) > In1 (−28.0 kcal/mol) > In4 (−23.2 kcal/mol) > In2 (−15.2 kcal/mol). The measured IC50 values (130.3, 318.4, and 20.3 nM)[67] for analogues of three inhibitors (In1, In2, and In5) are in excellent agreement with the order predicted by binding free energies. Even though the binding property of In3 to PSMA is the best, the boron density of In5 is the highest and its binding potency is only slightly lower than that of In3. By taking into account the balance between binding potency and boron density, In5 should be considered as the best option for BNCT. In summary, our computational study can complement experimental binding studies between PSMA and potential BNCT reagents that are typically performed by isothermal calorimetry (ITC) or Trp fluorescence spectroscopy. It provides the detailed atomistic and thermodynamic properties of each enzyme–inhibitor complex that may be useful for the design of In5 like PSMA inhibitors. This information is expected to stimulate the development of effective PSMA-targeted BNCT reagents with both high boron density and enhanced binding affinity. However, they need to be synthesized and characterized using experimental techniques.

Computational Details

The 1.84 Å resolution crystal structure of PSMA (PDB ID: 4NGM) was obtained from Protein Data Bank (PDB).[55] The structures of various urea-based boron-containing inhibitors were fully optimized without any geometrical constraint at the B3LYP[82]/6-31G(d)[83] theory level by utilizing the Gaussian 09 program.[84] These optimized structures were used to develop their force field parameters utilizing Automated Topology Builder (ATB).[85] These parameters are compatible with the GROMOS force field family in a wide range of formats. The molecular docking procedure was performed using Autodock Vina 1.5.6 software[86] to derive the binding poses of inhibitors to the bimetallic active site of PSMA. Two docking approaches (rigid docking and flexible docking) were used to yield 20 poses with an exhaustiveness value of 20 for each trial, and the size of the grid box was chosen to cover the entire active site of PSMA. The lowest energy poses provided by these docking protocols were used as the starting points for MD simulations, which were performed using the GROMACS program[87] and the Gromos 54A7 force field.[88] The enzyme–substrate complexes were equilibrated in a cubic box with dimensions of 10 × 10 × 10 nm for all simulations to obviate the adverse impacts that maybe induced by the application of periodic boundary conditions.[89] This box was saturated with SPC water molecules,[90] some of which were replaced by sodium and chloride ions to simulate a physiological ion concentration of 154 mM and neutralize the system. Since the active site of PMSA is deeply buried away from the water solvent, SPC is an adequate water model for this system. The energy minimization with 3000 steps was performed on the starting structures through the steepest descent method with the coordinates of the active site and inhibitor frozen. The resulting structures from that energy minimization were first equilibrated for 20 ns by placing distance restraints around the active site and substrate to minimize the energy of the environment. Subsequently, the MD simulations on these equilibrated structures were performed for 100 ns without any restraints. All these simulations were carried out using a constant number of particles (N), pressure (P), and temperature (T), that is, NPT ensemble or isobaric ensemble. The LINCS algorithm[91] was employed to constrain the bond lengths of the peptide, and the SETTLE algorithm[92] was used to constrain the bond lengths and angles of water molecules. The long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) method[93] with a cutoff value of 1.2 nm. Peptides, metal ions, inhibitors, water molecules, and ions were coupled separately in a bath at 1 atm and 300 K with coupling constants of 1.0 and 0.1 ps, respectively. A time step of 2 fs was utilized to compute the trajectory of each model, and a pH of 7.0 was used to set amino acid residues to their normal ionization states.

Several built-in tools within GROMACS were employed to analyze the trajectories obtained from MD simulations. The most representative structures of enzyme–substrate complexes were derived from the cluster analysis, in which the frame with the greatest number of neighbors was selected as the middle structure to represent the cluster that has been constructed by grouping together the structurally similar frames (RMSD cutoff of 0.3 nm). The computed RMSD and pairwise RMSD[94] values confirmed the convergence of all equilibrated structures within the simulation timeframe. The binding free energies between PSMA and inhibitors were calculated by the lambda (λ) particle approach[95−97] in a thermodynamic cycle that describes the bound and unbound states. In this cycle, the relative binding energies between an enzyme and a ligand can be defined as the difference in free energy associated with the chemical changes of the ligand into the enzyme in their bound and solvated states. In this approach, the Coulombic and van der Waals interactions between PSMA and inhibitors were turned off in a sequential and systematic way: first Coulombic interactions were turned off and then the van der Waals interactions. This sequence of operations avoided the interactions of opposite charges at intimate distances, which could lead to unrealistic configurations and imprecise energies. The soft-core interactions[98] were applied to Lennard-Jones and Coulomb potentials to remove singularities in these potentials and circumvent a poor convergence at λ close to 0 or 1. The values of soft-core parameters (α, sc_alpha in mdp), soft-core power (λ, sc_power in mdp), and radius of interaction (σ, sc_sigma in mdp) when either C6 or C12 is zero were set to 0.5, 1.0, and 0.3, respectively. NCIs between inhibitors and PSMA were calculated and visualized using the NCIPLOT program,[99] which identifies the NCIs based on the electron densities and their derivatives. The contact maps and defined secondary structure protocol (DSSP)[100] of all complexes were generated from their compacted trajectory files (xtc files) by inbuilt tools of mdmat and DSSP within the GROMACS software package. The zinc metal ions, inhibitors, and S1 site residues were chosen as references to highlight their relative positions in the contact maps. They were labeled by the residue index in a numerical order as listed below: 1-Zn1, 2-Zn2, 3-inhibitor, 4-μOH, 5-Tyr205, 6-Lys207, 7-Phe209, 8-Arg210, 9-Tyr234, 10-Glu457, 11-Arg463, 12-Arg534, 13-Arg536, 14-Lys545, 15-Phe546, 16-Ser547, 17-Gly548, 18-Tyr552, 19-Asn698, 20-Tyr700, and 21-Ala701. ESP of PSMA was created using Adaptive Poisson-Boltzmann Solver software[101] and visualized with PyMol.[102] In addition, VMD,[103] ChemDraw,[104] YARASA,[105] and Chimera[106] software programs were also utilized to visualize and prepare the diagrams used in this study. The conformational dynamics of all complexes were investigated by performing principal component analysis (PCA)[107,108] of the alpha carbon atoms. The PCA reduces the dimensionality of large data variables, while preserving as much information as possible.