Shao-Long Lina,
Yan-Song Chenb,
Ruo-Yu Liua,
Mei-Ying Zhuc,
Tian Zhuc,
Ming-Qi Wang*b and
Bao-Quan Liu*a
aDepartment of Bioengineering, College of Life Science, Dalian Minzu University, Dalian 116600, China. E-mail: lbq@dlnu.edu.cn
bSchool of Pharmacy, Jiangsu University, 212013, Zhenjiang, PR China. E-mail: wmq3415@163.com
cBeijing Bionaxin Biotech Co, 1000000, Beijing, China
First published on 13th March 2024
Prostate-specific antigen (PSA) serves as a critical biomarker for the early detection and continuous monitoring of prostate cancer. However, commercial PSA detection methods primarily rely on antigen–antibody interactions, leading to issues such as high costs, stringent storage requirements, and potential cross-reactivity due to PSA variant sequence homology. This study is dedicated to the precise design and synthesis of molecular entities tailored for binding with PSA. By employing a million-level virtual screening to obtain potential PSA compounds and effectively guiding the synthesis using machine learning methods, the resulting lead compounds exhibit significantly improved binding affinity compared to those developed before by researchers using high-throughput screening for PSA, substantially reducing screening and development costs. Unlike antibody detection, the design of these small molecules offers promising avenues for advancing prostate cancer diagnostics. Furthermore, this study establishes a systematic framework for the rapid development of customized ligands that precisely target specific protein entities.
PSA, an enzyme of the serine protease class, is secreted by prostate epithelial cells and presents as a glycoprotein composed of 237 amino acid residues.4 Prostate-specific antigen is also a natural constituent of semen, where its primary function lies in facilitating the hydrolysis and liquefaction of semen clots, thereby aiding in the release of sperm.8 The Prostate Health Index (PHI) and other composite indices evaluating prostate health, incorporating comprehensive assessments of cPSA, fPSA, and p2PSA, have significantly enhanced the detection rate of prostate cancer.7
Currently, a diverse array of PSA detection chips and assay kits is commercially available, primarily relying on the principle of antigen–antibody specific recognition. However, it is important to note that antibodies used in these assays are associated with considerable costs, necessitate stringent low-temperature storage conditions, and contribute to elevated detection expenses.9 Concurrently, the considerable amino acid sequence homology among cPSA, fPSA, and p2PSA poses a significant challenge, as it can readily lead to antibody cross-reactivity, a phenomenon readily discerned when employing these aforementioned chips and kits.10,11
Interface diversity is crucial for understanding protein functions, designing drug molecules, and developing new antibodies.12,13 It enables proteins to interact with a variety of partners, thus participating in a wide range of biological processes. For instance, ab initio docking methods for antibodies play a significant role in leveraging interface diversity, especially in the recognition of prostate-specific antigens.14 Although antibody recognition is a commonly used approach, exploiting the potential of interface diversity for the combined use of ligands and antibodies in recognizing prostate-specific antigens could offer innovative potential.
Virtual screening has emerged as a pivotal tool within the domain of small molecule development, encompassing applications in drug development and the preparation of molecular probes. This computational approach comprises two fundamental strategies: structure-based design and ligand-based design. Structure-based design capitalizes on the spatial arrangement of proteins and ligands, as well as the interactions with local residues, to construct a lock-and-key model for guiding drug development and molecular probe formulation.15 In contrast, ligand-based design encompasses methodologies such as Quantitative Structure–Activity Relationship (QSAR) and pharmacophore modeling, aiming to identify ligand structures exhibiting structural or functional similarities to known ligands16 The integration of both structure-based and ligand-based virtual screening methodologies has demonstrated notable success in drug development endeavors, thereby enhancing the likelihood of successfully identifying target compounds.17,18 While direct utilization of open source databases like ZINC enables the screening of a vast array of small molecule compounds with binding potential to target proteins, challenges arise when dealing with unknown proteins. In such cases, the establishment of effective pharmacophore models becomes impractical, thereby necessitating substantial efforts for compound property verification. Furthermore, it is worth noting that many compounds within these open databases are characterized by their intricate synthesis requirements and high production costs.19
This research endeavors to elucidate the design and synthesis of a small molecule compound ligand tailored to the distinct binding pocket of prostate-specific antigen (PSA). We employed a comprehensive analytical approach, which incorporated a synergistic integration of Virtual Screening and Machine Learning methodologies, supported by a diverse spectrum of computational and empirical methods, including Surface Plasmon Resonance (SPR), Microscale Thermophoresis (MST), Circular Dichroism (CD), and Molecular Dynamics (MD) simulations, to meticulously evaluate the binding affinity of the synthesized compound with the target protein (refer to Fig. 1).20–23 This endeavor has yielded a promising lead compound characterized by a robust binding affinity, holding significant potential for advancing the diagnosis and treatment of prostate cancer. Moreover, it lays the groundwork for a viable development pathway for designing ligands targeted at specific protein. This multifaceted approach not only facilitates the development of a lead compound for prostate cancer management but also exemplifies a systematic framework for the rational design and synthesis of binding ligands tailored to target proteins.
For the exploration of potential binding sites within the prostate-specific antigen, we harnessed the capabilities of the PlayMolecule BindScope module. This innovative module leverages voxelization techniques to delineate binding pockets and ligand poses based on distinctive pharmacophore class attributes. Notably, it employs a 3D convolutional neural network (CNN) to predict binding affinities, thus enhancing the precision of our virtual screening efforts.26 Furthermore, we utilized drugs targeting prostate-specific antigen developed by scientists like LeBeau AM for the full protein docking of prostate-specific antigen. This aligns with the cavity prediction by the aforementioned software, therefore, we are confident in designating this cavity as the one for virtual screening.27
In the realm of virtual screening, Autodock Vina emerged as our tool of choice. It capitalizes on a global optimization algorithm, a marked improvement over the genetic algorithm employed by AutoDock 4, thereby bolstering search efficiency and result accuracy. Moreover, Autodock Vina offers support for multithreaded parallel processing, augmenting computational throughput.
Our computational endeavors were executed on a robust 192-core Linux server, enabling the completion of all docking simulations within a time frame of 30 days. The simulation parameters were set with a box size of 60 × 60 × 60 Å3, centered at coordinates X = −30.53, Y = −30.375, and Z = −24.965. A search accuracy level of exhaustiveness = 32 was employed, and the analysis focused on the ten best docking poses selected from the screening results.
Initially, the SMILES representations of the selected compounds postdocking are converted into molecular fingerprints using the RDK software. Additionally, a binary coding or counting vector representation is created for the SMILES notation of chemical and physical properties, en-compassing molecular weight, polarity, charge distribution, stereochemistry, and other relevant features.32 Principal Component Analysis (PCA) is subsequently employed to reduce the dimensionality of the data, ultimately resulting in a 166 bit binary evaluation index serving as the input file for K-means cluster analysis.33
Utilizing the generated algorithm input file, a random function is formulated to select the initial cluster centers for the clustering process. The Euclidean distance metric is then utilized to determine the nearest cluster center for each data point, thereby minimizing the sum of distances between each data point and its respective cluster center. The first iteration produces the initial cluster distribution by assigning data points to the nearest cluster centers. Subsequently, the mean value of data points within each cluster is calculated to establish new cluster centers and minimize intracluster variance. This step is iteratively applied until the cluster distribution stabilizes, indicating the completion of the clustering process.30
In this study, the K-means clustering model adopts the Euclidean distance metric for assigning data points to their respective clusters. The calculation formula is as follows:
![]() | (1) |
To evaluate the quality of the clustering process and assess the cohesion and separation of the clustering outcomes, the Silhouette Coefficient and Calinski–Harabasz Index are employed as performance indicators.
The Silhouette Coefficient is calculated as follows: For each data point i, the average Euclidean distance between i and other members within the same cluster is computed as a(i). Additionally, b(i) is determined, representing the average Euclidean distance from data point i to all members within other clusters. Subsequently, the silhouette coefficient for each data point i is calculated as follows:
![]() | (2) |
For the entire cluster, the average silhouette coefficient is computed, with n signifying the total number of cluster members:
![]() | (3) |
The Silhouette Coefficient yields values within the range [−1, 1]. A higher Silhouette Coefficient indicates superior clustering results, implying that data points are tightly clustered within their respective clusters while maintaining clear separation from other clusters. A coefficient closer to 1 signifies optimal clustering outcomes, while a value near 0 suggests potential overlap among clusters. Conversely, a coefficient near −1 indicates potential misclustering of data points.31
The Calinski–Harabasz Index is computed as follows:
The inter-class dispersion B of the clusters is evaluated. Inter-class dispersion measures the difference between each cluster and the sum of squared distances between all cluster members and the overall average
![]() | (4) |
To iteratively expand the subgraph, all CSGs are amalgamated, adhering to the constraint of preserving the molecular connectivity. The parameters -atomcompare- and -bondCompare- serve as pivotal tools for evaluating atom and bond equivalences, respectively. Through the meticulous application of these parameters, the framework succeeds in extracting the Maximum Common Subgraph (MCSG). Ultimately, this process culminates in the extraction of the compound's skeletal characteristics, thereby facilitating in-depth analysis and further scientific inquiry.34
To prepare the CM5 chip for target protein immobilization, the carboxylic acid groups on the chip's surface (Cytiva) were initially evaluated using a solution comprising a mixture of EDC and NHS (Cytiva) at a temperature of 25 °C and a flow rate of 10 μL min−1. This activation process spanned a duration of 7 minutes. Following activation, the chip was subjected to an injection of a sodium acetate soluble buffer (10 mM; pH 4.5) containing the target protein until the protein content immobilized on the chip reached a level of 2500 resonance units (RU). Ultimately, the multicycle kinetics of the chip were quenched using ethanolamine.35
Subsequently, the target protein was securely anchored to the CM5 chip, and the gradient-diluted analytes were consecutively introduced, with each concentration examined during a distinct cycle. During these experiments, solvent corrections were diligently applied using DMSO to ensure precision in the evaluation of affinity and kinetics. The concentration series of small molecules tested included the following values: 0, 1.25, 2.5, 5, 10, and 20 nM. The experimental conditions were maintained at a temperature of 25 °C, a flow rate of 30 μL min−1, a binding duration of 90 seconds, a dissociation duration of 60 seconds, and a running buffer comprising 5% DMSO in PBS-PSA.
The CD experiment was conducted using the Chirascan circular dichroism spectrometer produced by Applied Photophysics. The cuvette had a volume of 400 μL, with a PSA concentration of 2 μM. The buffer solution used was PBS with a concentration of 10 mM and a pH of 7.4. During the experiment, a scanning wavelength range of 180–260 nm was employed. In the cuvette, 1.5 μL of a 100 μM small molecule was added sequentially with thorough mixing between additions. Circular dichroism spectra were recorded after each compound addition to investigate the impact of compound titration on the protein.
The obtained results were subjected to analysis using the Circular Dichroism by Neural Networks (CDNN) software. CDNN utilizes neural networks to predict the secondary structure of proteins based on CD spectroscopy data.37
The Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) method is widely employed for the analysis of binding affinities between proteins and ligands.40 Its binding free energy calculation formula includes contributions from van der Waals interactions, electrostatic interactions, polar solvation, and solvent-accessible surface area (SASA). Free energy landscapes provide insights into the interactions and energy distribution among molecules within the system, aiding in the understanding of molecular interactions, conformations, and stability characteristics. Typically, free energy values are represented using color intensity and contour lines, with the X-axis representing RMSD and the Y-axis representing R(g) (Radius of Gyration). R(g) is a measure of molecular volume and compactness, where smaller R(g) values indicate a compact molecular structure and larger R(g) values indicate a more extended or loose molecular structure.
To streamline the characterization and synthesis efforts following virtual screening, K-means clustering analysis was applied to the 1160 compounds with different cluster numbers (K = 2, 3, 4, 5). The results, as depicted in the Fig. 2, show Silhouette Coefficient values ranging from −1 to 1, with larger values indicating better clustering performance. The Silhouette Coefficient was highest at 0.60 for K = 2, making it the optimal choice among the various clustering options. Subsequently, the compounds of Cluster I, identified through this analysis, were used for further investigation.
To characterize the central compounds within the clusters, the Euclidean distance formula was utilized to calculate the cluster center for Cluster I (ZINC000013683205) and its surrounding compounds, namely ZINC000028333181, ZINC000257278584, and ZINC000257297002, as shown in Fig. 3B. Subsequent sections will provide a detailed analysis of these four structures, including molecular similarity and their interactions with protein cavities.
In Fig. 3, we employed PYMOL for visualization.42 It can be observed that these four compounds bind to the same cavity. The compounds ZINC000013683205, ZINC000028333181, ZINC000257278584, and ZINC000257297002, enclosed in the red box, share similar structural motifs, benzimidazolone. This region is primarily involved in hydrogen bond interactions, driven by the structural similarities of these molecules. The structures enclosed in the blue box exhibit greater structural diversity and pose greater synthesis challenges. The blue box primarily consists of hydrophobic groups based on benzene rings, with hydrophobic interactions being the dominant mode of interaction, as indicated in the accompanying table. The central part of the structure includes amide and other structural elements contributing to fewer interactions, suggesting its role as a linker between hydrogen bonding moieties and variable hydrophobic groups.
Due to the sensitivity of the K-means algorithm to initial centroid selection, there is a risk of converging to local optima. This study combined the features of central compounds and their surrounding compounds, along with insights into the cavity's structural characteristics, to further deduce the skeletons of the 1160 compounds, aiding in guiding compound synthesis. The maximum common substructure (MCS) for this cluster was established by employing rdkit software and the VF2 graph matching algorithm on the results of the K-means clustering (K = 2). Ultimately, two maximum common substructures were identified (Fig. S1A and B†).
![]() | ||
Fig. 4 Maximum common substructures selected by K-means clustering and newly designed molecules, with (A) and (B) showcasing the newly designed compounds LIG1 and LIG2, respectively. |
The synthesized compounds conform to the key functional group requirements, skeleton criteria, binding cavity specifications, and exhibit lower synthetic complexities, as depicted in Fig. S2.† The nuclear magnetic resonance spectra of these molecules are provided in ESI Fig. S3.†
In the SPR experiment, all obtained data underwent kinetics/affinity fitting analysis, employing a 1:
1 binding model. Data analysis was conducted utilizing Biacore evaluation software (version T200 2.0). The results revealed a binding affinity of 59.33 nM between LIG 1 and PSA, whereas the binding affinity between LIG2 and PSA was determined to be 144 nM (Fig. 5A–D). This observation signifies that LIG 1 demonstrates a notably higher binding affinity towards PSA.
In contrast to SPR principles, MST relies on the influence of temperature gradients on the diffusion behavior of molecules in solution. When two molecules bind or dissociate, their diffusion behavior in a temperature gradient undergoes changes. In this study, we validated the binding of PSA with its ligands in a multisystem approach, providing more comprehensive and reliable results. MST results showed that between concentrations of 0.15 nM to 5 μM, as the concentrations of LIG 1 and LIG 2 increased, the fluorescence signal gradually intensified, indicating that both LIG 1 and LIG 2 could bind to PSA (Fig. 5G and H). The binding affinity between LIG 1 and PSA was determined to be 78.91 nM, while the bind-ing affinity between LIG 2 and PSA was 106.31 nM, indicating that both LIG 1 and LIG 2 possess strong binding affinities, with LIG1 exhibiting stronger binding to PSA.
To further explore the impact of synthesized compounds on the spatial structure of PSA, circular dichroism (CD) titration experiments were conducted. As shown in Fig. 5I and J, it can be observed that upon addition of the compounds, the negative peak height of PSA in the range of 200 nm to 210 nm decreased with increasing compound concentration. This change was further analyzed using CDNN software to assess the impact of compound addition on protein structure. As depicted in Fig. S4,† in the absence of compounds, beta sheets constituted the major part of the protein, accounting for approximately 44% ± 1%, while alpha helices were the least prevalent, comprising approximately 23%. Upon the addition of LIG 1, there was a slight decrease in beta sheet structure by approximately 1%, and when the ratio of alpha helix to PSA and LIG 1 was 1:
3, there was a reduction of 3%, resulting in an alpha helix composition of approximately 15%, along with an increase in random coil structure by 5%, accounting for approximately 13.5% of the final structure (as shown in S4A†). In the case of LIG 2 addition, there was a slight increase of approximately 1% in beta sheet structure, and when the ratio of alpha helix to PSA and LIG 2 was 1
:
3, there was a reduction of 5%, resulting in an alpha helix composition of approximately 29%, along with an increase in random coil structure by 5%, accounting for approximately 12.8% of the final structure (as shown in Fig. S4B†). In summary, given the structural similarity between LIG 1 and LIG 2, their impact on protein structure exhibited a consistent trend, as both compounds significantly altered the PSA structure by converting alpha helices into random coil configurations upon binding with PSA.
ID | Docking score (kcal mol−1) | Hydrophobic residues | Hydrogen bond residues |
---|---|---|---|
LIG 1 | −10.8 | His 57, Leu 95, Cys 191, Ser 192, Gly 216, Cys 220 | Ser 189, Thr 190, Ser 217 |
LIG 2 | −10.3 | His 57, Leu 95, Asp 102, Cys 191, Cys 220, Gly 216 | Thr 190, Ser 193, Ser 217 |
Drug 1 | −8.9 | Leu 95, Cys 191, Trp 215 | His 57, Thr 190, Ser 192, Gly 193, Ser 195, Cys 220, Ser 226 |
Drug 2 | −8.1 | Leu 95, Thr 190, Try 215 | ALA39, His 57, Ser 99, Thr 190, Ser 192, Gly 216 |
To further evaluate the stability of compound-PSA-ligand complexes and ligand conformations, various indicators including RMSD (Root Mean Square Deviation), RMSF (Root Mean Square Fluctuation), and MMPBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) were assessed in this study (Fig. 6E–H). Simulations for all structures were carried out over a 300 ns timescale. The RMSD data reveals that the average RMSD values for LIG1 and LIG2 are 0.1293 nm and 0.1321 nm, with standard deviations of 0.0132 nm and 0.0144 nm, respectively. For Drug 1 and Drug 2, the average RMSD values are 0.1628 nm and 0.1621 nm, with standard deviations of 0.0165 nm and 0.0170 nm, respectively. Lower RMSD values indicate greater system stability, and both the average and variance values for LIG 1 and 2 are lower than those for Drug 1 and Drug 2, suggesting better stability. Overall, all RMSD values are below 0.3 nm, indicating the stability of the four compound structures and validating the rationality of the system setup for further in-depth analysis.
MMPBSA results are presented in Table 2. Among the four systems, the PSA-LIG1 complex exhibits the lowest binding free energy of −36.60 kcal mol−1. Although the free energy scores for LIG2 and Drug 1 are relatively close, with LIG2 scoring −35.36 kcal mol−1 compared to Drug 1 scoring 33.89 kcal mol−1, LIG2 still maintains a lower free energy score. To enhance visualization, MMPBSA results were visualized using the builtin plugins of Gromacs 2022.3.38,39
Complex | Binding free (kcal mol−1) | Van der Waals (kcal mol−1) | Electrostatic (kcal mol−1) | Polar solvation (kcal mol−1) | SASA (kcal mol−1) |
---|---|---|---|---|---|
PSA-LIG1 | −36.60 | −46.53 | −23.87 | 38.75 | −4.95 |
PSA-LIG2 | −35.36 | −34.97 | −36.75 | 48.87 | −12.51 |
PSA-Drug1 | −33.89 | −54.72 | −51.95 | 80.02 | −7.24 |
PSA-Drug2 | −30.59 | −54.33 | −31.80 | 63.41 | −7.87 |
In this study, 2D and 3D free energy landscape plots were constructed. The rmsd_gyrate.log, bindex.ndx, and rmsd_gyrate.log files recorded the relationship between indices and energies. By analyzing these files, the conformations corresponding to the lowest energy for each system were identified at 266 ns, 85 ns, 114 ns, and 71 ns for LIG1, LIG2, Drug1, and Drug2, respectively (as shown in Fig. S6†). The RMSD values for these conformations are close to the average, and R(g) values are relatively low, indicating structural similarity to the initial structures and compactness in the low-energy region. This suggests that during the simulation, molecules may remain in this state for an extended period, demonstrating stability and further confirming the reliability of the binding free energy calculations.
High-throughput screening techniques enable the detection of alterations in enzyme or receptor function, the interaction between probes and proteins, as well as the kinetic properties of protein-ligand binding.46 Nevertheless, the scarcity of compound samples for high-throughput screening poses a challenge: a substantial number of compound samples are required, but obtaining these samples can be difficult.27,43–45 For instance, compounds targeting prostate-specific antigen, designed by previous researchers using high-throughput screening, entail considerable experimental costs. In contrast, this study employs an integrated approach, combining computer aided drug design and machine learning techniques, to rapidly and cost-effectively identify and synthesize ligands with high binding affinity for prostate-specific antigen within the diverse ZINC database.
In this study, a comprehensive evaluation of compound-protein binding affinity is conducted using four techniques: molecular docking, molecular dynamics simulations, surface plasmon resonance (SPR), and microscale thermophoresis (MST). SPR is an optical phenomenon that arises when incident light at the interface between a metal and a dielectric medium meets specific energy and momentum matching conditions, leading to the excitation of coherent oscillations of free electrons on the metal surface.20 MST employs infrared lasers for localized heating of molecules, inducing directional movement, and subsequently analyzing the molecular distribution ratio within the temperature gradient field via fluorescence analysis.21 Despite the distinct principles underlying the two systems, the experimental values obtained are comparable. Although the data for LIG 1 and LIG 2 are relatively similar, both systems demonstrate that LIG 1 exhibits a higher binding affinity for prostate-specific antigen than LIG 2, corroborating the molecular docking and molecular dynamics simulation data. This, to some extent, validates the reliability of the computational data and ultimately establishes a dependable, rapid, and cost-effective ligand screening approach for specific proteins.
However, our work also revealed certain challenges. Conventional antibody-based prostate cancer detection methods exhibit limitations related to temperature sensitivity, cost, and cross-reactivity. While our synthesized compounds exhibit lower binding affinities, they address cost and storage issues. Nevertheless, there remains room for optimization. Future research directions could include enhancing compound binding affinities through fragment-based drug design and designing residues for multiple PSA cavities.47 Additionally, cost-effective fluorescent labeling methods, such as incorporating Fmoc groups, could be explored to improve PSA detection.48
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3ra08550c |
This journal is © The Royal Society of Chemistry 2024 |