The complexes were first sketched using the 2D sketcher of Maestro (Schrödinger), by drawing the triphenylphosphine, bisdiphenylphosphines, and seven ligands as building blocks for the final assembly of the mono and binuclear Phosphinogold(I) thiocarbohydrate complexes as illustrated in Figure 1. The sketched molecules in Figure 1 were then assembled as the 20 mono and binuclear complexes described by (Adokoh et al., 2017) utilizing the single complex builder within Maestro (Schrödinger, 2023a). This was done by selecting gold(I) as the central atom(s), and linear geometry for the final complex around the central atom (Au). The complex builder was run to join the gold atom(s) to the phosphine’s phosphorus atom(s) at one end, and to the sulfur atom(s) of the thiocarbohydrate ligand(s) at the opposite end. The final assembled complexes are illustrated in Figure 1. The LigPrep tool in (Schrödinger, 2023b) was used to convert the 2D structures into 3D structures. Subsequently, each complex’s geometry was optimized using OPLS4 force field minimization by running the complex cleanup tool in Maestro. The final library of prepared complexes was then used to perform the subsequent molecular docking and MM-GBSA free energy calculation studies.

The 20 complexes have previously shown remarkable anticancer activities in 3 cell lines of the breast (MCF7), prostate (PC3), and colon (HCT116) cancers (Adokoh et al., 2017). However, the mechanism of action of these complexes is not yet known. To investigate the possible mechanisms for the observed anticancer activities, a thorough literature review was initiated to look for the available known targets involved in the three types of cancer cells that have significant effects on their life cycle. The review discovered 17 targets that could be investigated through an in silico approach to evaluate their binding interactions and to estimate binding affinities toward the Phosphinogold(I) Thiocarbohydrate Complexes, which will make it possible to hypothesize possible mechanisms for the observed anticancer activities. Table 1 summarizes the selected targets and their involvement in various types of cancer. The crystallized structures of the targets were downloaded as pdb files from the protein data bank, accessible at https://www.rcsb.org/, then they were loaded into the protein preparation workflow in Maestro (Schrödinger, 2023c) at the default settings. Protein preparation for molecular docking simulations was conducted through a series of sequential steps, involving filling in missing side chains, assigning bond orders, placing hydrogens, generating het states at pH 7.4 ± 2, deleting bulk waters, optimizing hydrogen-bond assignments at pH 7.4, and lastly, energy minimization of the protein using the OPLS4 force field. The co-crystallized ligands or inhibitors (Table 1) from the protein targets were extracted, prepared in LigPrep (Schrödinger, 2023d), and subsequently used to validate docking protocols and they were used as reference to compare binding affinities across the 20 phosphinogold(I) thiocarbohydrate complexes using Glide (Schrödinger). It is important to mention that, the target protein β-catenin with PDB ID: 1JDH had no cocrystallized ligand or inhibitor, therefore, the binding site was determined utilizing the SiteMap function within Maestro software, and for subsequent docking studies, PRI-724 (also known as Foscenvivint); a known potent inhibitor for β-catenin signaling pathway (Schmidtova et al., 2021), was selected as the standard ligand for 1JDH to compare binding affinity between ligands to this target.

Molecular docking was performed using Glide (Schrödinger, 2023a) software. The procedure was initiated by selecting each prepared target protein as the macromolecule and the 20 prepared gold(I) complexes as the ligands. A receptor grid box using glide was previously generated, centered on the position of the cocrystallized ligand, and a midpoint box of 10 Å diameter in all three coordinates. Flexible ligand sampling with extra precision (XP) was selected to run the molecular docking simulations. The output settings of the docking results were kept at their default values to report the best pose with the highest docking score for each ligand.

The reliance on redocking in validating docking protocols (Table 2) and assessing accuracy will help assess benchmarking performance. When there is a high level of repetition between the experimentally determined and docked poses, it gives confidence that the docking results reflect, to a high degree, what is happening. In the present study, the native ligands were docked using the same grid box generated for ligand docking, and the poses of the co-crystallized and docked ligands were compared by calculating root mean square deviation (RMSD) (Figure 2). Following this outcome, if the RMSD is small (preferably less than 2.0), then the docking protocol is considered valid (Hevener et al., 2009; Jain, 2008).

The MM-GBSA method combines molecular mechanics force fields (MM) with implicit solvation models (Generalized Born and Surface Area, (GBSA)) to estimate the binding free energy (binding affinity). MM-GBSA is computationally efficient and often used as a first-pass estimate of binding affinity. Prime MM-GBSA can be performed without running a full-fledged MD simulation, utilizing the docked poses as a starting point.

Molecular docking relies on scoring functions to assess the quality of the predicted binding poses generated during conformational search. These functions utilize approximations to streamline the calculations, enabling high-throughput screening of potential drug candidates. Due to the use of approximate scoring functions, molecular docking often falls short of accurately predicting binding energies when compared to experimental measurements. Although, numerous docking programs effectively identify potential ligand binding conformations, a universal scoring function that accurately predicts binding energies for all molecules and protein families remains elusive. Consequently, rescoring steps after molecular docking are often essential to refine the initial predictions (Zhang et al., 2017). On the other hand, MM-GBSA utilizes molecular mechanics and free energy calculations to consider entropic contributions and thus provide more accurate binding affinity predictions (Sgobba et al., 2012). MM-GBSA is computationally more expensive but offers greater accuracy and insights into binding mechanisms, making it ideal for ranking and refining promising ligands (Genheden and Ryde, 2015; Sgobba et al., 2012). Instead of simulating the system’s molecular dynamics over time, the Prime MM-GBSA module in the Schrodinger suite, performs a series of calculations on a single optimized pose from Glide’s docking output, using a Generalized Born solvation model (VSGB 2.0) (Li et al., 2011) to estimate the binding free energy. These calculations involve minimization to relax any strained interactions, energy calculation of the complex and its components using the MM-GBSA method, and finally, determining the binding affinity by calculating the difference in energies (E) of the protein-ligand’s complex and its components (Muddagoni et al., 2021): Δ G bind = E _ complex   minimized   –   E _ ligand   minimized   –   E _ receptor   minimized

This approach offers significant time efficiency compared to standard MM-GBSA, which relies on lengthy MD simulations. However, it is important to note that this method is limited by the single, optimized docked pose, potentially missing the full range of conformations the ligand can adopt in the binding site and potentially neglecting entropic effects. Overall, Prime MM-GBSA without MD can be a useful tool for quick binding affinity estimation. However, for more accurate and complete understanding of the binding process, a full MD simulation is recommended. Therefore, in our study, the poses obtained with the highest docking score for each ligand across all 17 targets were re-scored using the Prime (Schrödinger, 2023a) MM-GBSA to calculate binding free energy (ΔG) in kcal/mol. Then, MD simulations were performed for the best scoring complex.

Building upon the promising MM-GBSA and docking results, Complex 9, which emerged as the top contender from a group of twenty Phosphinogold(I) Thiocarbohydrate Complexes, was subjected to detailed molecular dynamics (MD) simulations. The 100 ns simulations, performed using Desmond software (Schrödinger, 2018), aimed to explore the dynamic interactions between Complex 9 and its target proteins: Poly (ADP-ribose) polymerase (PARP-1, PDB ID 3L3M) and Protein Kinase B (AKT2, PDB ID 8Q61). To prepare the simulations, the protein-ligand complexes were placed in minimized, solvated orthorhombic boxes, ensuring a 10 Å buffer of TIP3P water molecules. The simulation environments were brought to physiological ionic strength by adding counter ions for neutralization and 0.15 M NaCl. The MD simulations were performed under NPT conditions (constant particle number, pressure, and temperature) at 300 K and 1.01325 bar, with the OPLS3e force field. Prior to the 100 ns production run, all systems were relaxed using Desmond’s default 160 picosecond protocol. Following the simulation, Maestro (Schrödinger, 2023b) was employed for comprehensive analysis of simulation data. Key parameters assessed included root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), the dynamic protein-ligand contact profiles, and finally, the calculation of post-simulation MM-GBSA binding free energy (ΔGbind in kcal/mol).

The reliability of the molecular docking protocol was established through redocking of native ligands into their respective protein structures, confirming accurate reproduction of experimentally observed binding poses (Figure 2; Table 2). Root mean square deviation (RMSD) was calculated to quantify the difference in position between co-crystallized and re-docked native ligands. An RMSD value equal to or less than 2.0 Å was used as a criterion to validate the docking protocol (Hevener et al., 2009; Jain, 2008). Following validation, molecular docking was performed using Glide XP, and subsequently, the more rigorous Prime MM-GBSA method was employed to refine binding affinity estimations. A detailed breakdown of observations for each complex can be found in the supplementary materials (Supplementary Tables S1–20; Supplementary Figurea S1–14).

The analysis identified several key protein targets, the most significant are Human Double Minute 2 protein (HDM2), DNA methyltransferase-1 (DNMT1), Human AKT2, and Poly (ADP-ribose) polymerase (PARP-1), which frequently exhibited strong interactions with multiple complexes. Number of complexes exhibited strong binding affinities, with several of them surpassing the binding strength of native ligands (Tables 3, 4). This suggests potential multi-targeting mechanisms and broad-spectrum activity. A general trend observed was that binuclear complex (those numbered 8 and above) consistently exhibited higher binding affinities than mononuclear complexes, as illustrated in Figures 3, 4. This suggests that the presence of two gold centers within the complex may be associated with enhanced binding. This relationship was confirmed by a statistically significant (α = 0.05) negative correlation between complex type (mononuclear or binuclear) and docking scores with two targets; DNMT-1 (Spearman’s correlation coefficient r = −0.837, p = 0.0061) and HDM2 (r = −0.717, p = 0.0242). The negative correlation indicates that binuclear complexes, which correspond to lower docking scores, bind more tightly to the proteins (Supplementary Tables S21, 22). The noted increase in binding affinity with binuclear complexes is entirely consistent with the results of Adokoh et al., who found that dinuclear gold(I) complexes (8-20) displayed significantly enhanced growth inhibition of cancer cells when compared to their mononuclear counterparts (1-7). Acetylation has been found to have a significant effect on binding affinity, with acetylated complexes (complexes 5-7 and 12-20) consistently exhibiting higher docking scores (lower binding affinity) (as detailed in Supplementary Tables S1, 20) than their non-acetylated analogs (complexes 1-4 and 8-11). This observation is supported by statistically significant positive correlations between acetylation and docking scores for DNMT-1 (Spearman’s r = 0.8452, p = 0.0016), HDM2 (r = 0.8563, p = 0.0001), PARP-1 (r = 0.8452, p = 0.0016), and AKT2 (r = 0.7681, p = 0.0081) (Supplementary Tables S22, 23). These positive correlations indicate that acetylation is associated with higher (less negative) docking scores, suggesting a weaker binding affinity. While this reduction in binding affinity can be explained by reduced hydrogen bond forming capability, Adokoh et al. observed that acetylation paradoxically enhanced the activity of gold(I) complexes (specifically complexes 5-7 in their study). This difference between the in silico finding and the experimental in vitro activity underscores that acetylation has a complex and multifaceted effect, where its impact on binding affinity does not necessarily correlate directly with its impact on overall cellular activity due to factors such as cell penetration, which Adokoh et al. attributed to increased lipophilicity as a result of reduced hydrogen bonding potential. The following illustrative examples underscore the significant interactions observed in this study. Table 4 provides a complete listing of the strongest interactions.

Complex 9: This complex stood out as a top performer against Human AKT2 isoform 8Q61, displaying an MM-GBSA binding score of −82.4 kcal/mol (Supplementary Table S9; Figure 5C), which exceeded even the native ligand’s score of −59.3 kcal/mol. This strong interaction aligns well with the strong anti-proliferative effects that have been seen with complexes 8-11 in vitro, and suggests that disruption of AKT2 may be a key factor in its success. Complex 9 also demonstrated strong affinity for PARP1 with MM-GBSA score of −75.7 kcal/mol and was therefore, considered for further MD study. Unfortunately, complex 9 was selected for in vitro studies at the time and however, it is imperative to further investigate complex 9. Complex 5 showed the highest binding affinity towards the HDM2 with a binding free energy of −76.0 kcal/mol. Several other complexes also showed promising interactions with HDM2, this suggests that inhibiting 5HMH/p53 interaction plays a crucial role in the mechanism of action of these complexes. Complex 10, demonstrated multi-target binding interactions, in particular HDM2, DNMT-1, and a particular affinity towards β-catenin (1JDH), were exceeding the MM-GBSA score of the reference inhibitor by a notable margin (Table 4). The tight binding of complex 10 with beta-catenin is crucial in the Wnt signaling pathway (Zhao et al., 2022) which supports its potent activity in colon cancer (HCT116 cell line) in vitro (IC₅₀ = 0.90 µM). Similarly, complex 11, also displayed a diverse array of binding targets, notably with DNMT-1, PARP-1, and the strongest binding affinity observed for the Alpha-Beta-tubulin dimer (1JFF) which could also account for strong in vitro anti-proliferative effects against colon and prostate cancer with IC_50S of 0.63 µM and 0.22 µM respectively (Adokoh et al., 2017). Complex 8 being the shorter chain of 9 and 11 demonstrated a significant interaction with both DNMT-1 and HDM2, surpassing the respective native ligand of PARP1 in line with in vitro data, which was the most potent growth inhibitor of prostate (PC3) cell line (IC₅₀ = 0.003 μM) (Adokoh et al., 2017).

Complex 9’s binding stability with PARP-1 and AKT2 was investigated via 100 ns MD simulations. As shown in Figure 6A and Table 5, the PARP-1-Complex 9 system achieved equilibrium after approximately 10 ns, exhibiting a stable protein Cα RMSD averaging 2.0 Å (range: 1.8–2.4 Å). This suggests minimal overall protein structural fluctuation. While the ligand RMSD relative to the protein (ligand fitted on protein) showed some fluctuations, reaching up to ∼6 Å in few instances, it was generally stabilized between 2.0 and 3.0 Å, averaging 2.72 Å. This suggests sustained binding despite some ligand mobility within the pocket, likely due to conformational adjustments. The ligand’s internal structure remained rigid, as evidenced by the low average RMSD of 2.21 Å (range: 1.6–2.5 Å) when fitted on itself. In contrast, the AKT2-Complex 9 simulation (Figure 6B) revealed a less stable protein Cα RMSD profile, although equilibrium was reached around 20 ns with an average of 3.03 Å (Table 5). Ligand mobility was more pronounced in this system, with an average ligand RMSD (relative to the protein) of 4.13 Å (Table 5). A period of increased fluctuation between 46 and 76 ns suggests substantial ligand conformational changes during this timeframe, potentially impacting the binding site and correlating with the observed jump in protein RMSD (Figure 6B). Despite this increased mobility relative to the protein, the ligand’s internal structure remained relatively rigid, averaging 2.74 Å RMSD when fitted on itself. This indicates the ligand’s conformational changes involve shifts in position or orientation within the binding pocket rather than significant internal rearrangements.

Figure 7 summarizes the properties of complex 9 bound to (A) PARP-1 and (B) AKT2 during MD simulations, including radius of gyration (rGyr), molecular surface area (MolSA), solvent accessible surface area (SASA), and polar surface area (PSA). Complex 9 displayed a more compact average structure and reduced solvent exposure with AKT2 (average rGyr ≈7.0 Å, range 6.2–8.6 Å; average SASA ≈458 Å², range 330–783 Å²) compared to PARP-1 (average rGyr ≈7.8 Å, range 6.8–8.9 Å; average SASA ≈428 Å², range 300–590 Å²), although the PARP-1 system showed greater rGyr fluctuation. Intramolecular hydrogen bonding was similarly low and variable in both simulations (average <2, range 0–4 H-bonds), indicating it is unlikely to be a key stabilizing factor. Polar surface area values were comparable (≈411 Å² for AKT2 and ≈423 Å² for PARP-1). These results indicate Complex 9 binds more dynamically to AKT2 than to PARP-1. Table 5 also summarizes the average values for these properties.

Analysis of Cα RMSF for PARP-1 bound to complex 9 (Figure 8A) suggests that ligand binding modulates protein flexibility. While the N- and C-termini appear generally rigid, ligand contacts at various locations, including residues 27-28, 55-60, 92-109, 200-275, and 322, may contribute to localized changes in flexibility. For example, increased RMSF values around residues 57-64 could indicate ligand-induced conformational changes. It is important to acknowledge that these observations are based on a single MD simulation and require further validation, including comparison to apo PARP-1 dynamics and experimental studies, to definitively establish the effects of complex 9 binding. On the other hand, complex 9 binding differentially affects the flexibility of AKT2. PARP-1’s N-terminus becomes rigid upon binding, whereas AKT2’s N-terminus shows a substantial increase in flexibility (i.e., higher RMSF) (Figure 8B). Overall, complex 9 has a more dramatic effect on AKT2’s flexibility, particularly at the N- and C-termini. In contrast, complex 9’s impact on PARP-1 is more localized. AKT2 also displays a greater number of contact points with complex 9. This broader interaction interface may suggest a more extensive allosteric effect of complex 9 on AKT2 compared to PARP-1. Further studies are needed to validate how these distinct dynamic responses relate to the functional regulation of each protein. Figures 9, 10 illustrate the interaction profiles of complex 9 with PARP-1 and AKT2 binding site residues, respectively. The prominent mode of interaction in both cases is through H-bonding and to a lesser extent, water bridges, which is a special form of H-bonds mediated by surrounding water molecules. This explain the higher affinities of non-acetylated series of complexes such as complex 9, compared to those of the per-acetylated complexes such as 12-20 as we have discussed earlier.

Prime MM-GBSA calculations, using 31 frames extracted at ∼3.3 ns intervals from 100 ns MD trajectories (Figure 11), provide a more comprehensive understanding of the binding profile of complex 9 with PARP-1 and AKT2 compared to previous calculations based on the single docked poses in the molecular docking study. These post-MD calculations, incorporating protein and ligand flexibility by sampling representative snapshots across the simulation, offer a more realistic representation of the binding dynamics. Complex 9 displays a more favorable and stable binding profile with PARP-1, exhibiting an average binding free energy of −71.82 kcal/mol and fluctuations mostly within a −50 to −91 kcal/mol range. In contrast, complex 9 binding to AKT2 is characterized by a lower average binding free energy (−63.31 kcal/mol) and wider fluctuations (−38 to −91 kcal/mol), suggesting a weaker or more dynamic interaction.