This is an uncorrected proof.
Figures
Abstract
Microtubules are dynamic cytoskeletal polymers assembled from /-tubulin heterodimers. Microtubules are validated targets for novel therapeutic drugs, yet therapeutic efficacy targeting them is often compromised by multidrug resistance (MDR). 5-(3-chlorophenyl)-N-(3-pyridinyl)-2-furamide (CPPF) is a novel microtubule-targeting anticancer agent which was found to disrupt microtubule growth in cells and inhibit tubulin polymerization in vitro. CPPF could suppress the growth of multidrug-resistant cell lines and demonstrate anti-tumor efficacy in animal models. However, the fundamental mechanism of how CPPF disrupts microtubule assembly remains unclear. To investigate this, we performed structure prediction using Protenix, RoseTTAFold All-Atom (RFAA) and Umol, as well as molecular dynamics (MD) simulations, using the human /-tubulin heterodimer (PDB ID: 5IJ0) for analyzing the interaction between CPPF and tubulin. Our results showed that CPPF binds at the / interface with dominant contributions from -tubulin residues, particularly VAL236 and LEU253. In isolated-monomer comparisons, -tubulin exhibited more favorable binding energetics and deeper, broader free-energy minima than -tubulin. Supplementary simulations on the alternative -tubulin conformational state (PDB: 6E7B; straight microtubule-lattice) indicated that CPPF binding is preserved across the two major -tubulin conformations, with comparable MM-PBSA binding free energies. Taken together, these findings establish a plausible binding mode for CPPF at the /-tubulin interface, providing a structural hypothesis for understanding its anticancer potential against multidrug-resistant cancers. Experimental validation through binding assays or crystallography is warranted to further substantiate these computational insights.
Author summary
Microtubules are essential for cell division and intracellular transport, which makes tubulin a major target for anticancer drugs. However, many tumors become resistant to commonly used microtubule-targeting agents, creating a need for new compounds with distinct binding mechanisms. CPPF is a recently reported microtubule inhibitor that remains active in multidrug-resistant cancer cell lines, but how it interacts with tubulin at the molecular level is still unclear. Here, we used an integrated computational framework that combines AI-based protein–ligand structure prediction (Protenix, RoseTTAFold All-Atom (RFAA), and Umol), binding-pocket analysis, residue-level interaction profiling, and all-atom molecular dynamics simulations to characterize CPPF–tubulin recognition. Across multiple models, CPPF showed more consistent binding poses and recurrent contact residues on -tubulin than on -tubulin, with key hotspots including VAL236 and LEU253. Simulations further suggested that the CPPF–-tubulin complex is dynamically more stable and occupies deeper free-energy minima than the -tubulin complex, with an ensemble-averaged MM-PBSA binding free energy indicating a substantially favorable interaction. Complementary simulations on an alternative -tubulin conformational state (the microtubule-lattice-related straight form) indicated that CPPF binding is preserved across both major -tubulin conformations, suggesting the pocket is accessible at multiple stages of microtubule dynamics. Together, our results propose a plausible binding mode for CPPF at a composite / interface pocket with dominant contributions from -tubulin residues, providing structural hypotheses to guide future biochemical validation and the rational optimization of tubulin inhibitors that may help overcome drug resistance.
Citation: Yang J, Liang L, Zhu D, Yin X, Feng Y, Tang S, et al. (2026) Integrative AI-assisted modeling suggests CPPF binding at a composite /-tubulin interface pocket dominated by -tubulin contacts. PLoS Comput Biol 22(9): e1014804. https://doi.org/10.1371/journal.pcbi.1014804
Editor: Mohammad Sadegh Taghizadeh, Shiraz University, IRAN, ISLAMIC REPUBLIC OF
Received: February 3, 2026; Accepted: September 7, 2026; Published: September 18, 2026
Copyright: © 2026 Yang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All processed data supporting this study, including structural metrics, interaction summaries, MD summary statistics, and free-energy landscape data, are provided in the Supporting Information. The complete computational workflow, analysis scripts, documentation, and reproducible Nextflow pipeline are available under the MIT License at github.com/jasperyeoh/integrative-ai-assisted-modeling-of-cppf-tubulin-interactions. All-atom MD production trajectories, paired run-input files, and SHA256 checksums are available at huggingface.co/datasets/jasperyeoh2/MD-trajectories-CPPF-tubulin-heterodimer-and-monomers.
Funding: The author(s) received no specific funding for this work.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Microtubules are dynamic cytoskeletal polymers playing fundamental roles in critical cellular processes such as intracellular transport, mitotic spindle formation, and directional cell migration [1–3]. Microtubules are assembled through GTP (guanosine triphosphate)-dependent polymerization of -/-tubulin dimers, which are the fundamental building blocks of microtubules [4]. Different - and -tubulin isoforms can confer unique biochemical properties and drug sensitivities, contributing to differential cellular responses to microtubule-targeting agents [5,6]. Some of these differences predominate in certain cancer cells and some are associated with the development of drug resistance [7]. For example, 3-tubulin (TUBB3) expression is up-regulated in several epithelial tumors [8] and, in lung cancer, has been associated with resistance to anti-programmed cell death protein 1 (PD-1) immunotherapy [7]. Disrupting these resistance-related processes therefore represents a promising therapeutic strategy. Microtubules are key targets in this context, with drugs such as taxanes, vinca alkaloids, and colchicine derivatives exerting anticancer effects by modulating microtubule dynamics [9]. Among the arsenal of anticancer agents, paclitaxel, marketed as Taxol, is a cornerstone drug that exerts cytotoxic effects by targeting microtubules, thereby blocking cancer cell division [10]. Paclitaxel binds to -tubulin near the lateral interface between protofilaments [11–13], stabilizing this flexible region. Importantly, this lateral interface is inherently dynamic; experimentally, microtubules can exhibit a range of 12–15 protofilaments [14,15] and tubulin dimers can switch between curved (unpolymerized) and straight (lattice-incorporated) conformations depending on their nucleotide state and binding environment [16–19].
5-(3-chlorophenyl)-N-(3-pyridinyl)-2-furamide (CPPF) was first identified in 2020 as a potent microtubule-targeting anticancer compound [20]. It was demonstrated that CPPF could inhibit microtubule growth in vitro and in vivo [20]. CPPF induced apoptotic cell death in various cancer cell lines from different origins (lung, breast, etc.), and remained effective against multidrug resistance (MDR) cancer cell lines, which are resistant to known microtubule-targeting drugs such as paclitaxel and colchicine [20]. Structural docking suggests the potential binding site of CPPF overlaps with the colchicine binding site at the /-tubulin interface, a region critical for microtubule depolymerization due to its hydrophobic nature [20]. This site’s potential to facilitate CPPF binding, combined with its efficacy in MDR cells, highlights its therapeutic promise. Furthermore, CPPF demonstrated antitumor efficacy by suppressing tumor growth in different animal models, and displayed low toxicity in initial tests [20]. However, the mechanism by which CPPF disrupts microtubule dynamics was unclear.
To understand the principle of how the novel microtubule inhibitor, CPPF, interacts with /-tubulin, we selected the high-resolution cryo-EM structure of the human /-tubulin heterodimer (PDB ID: 5IJ0) as a biologically relevant template [21]. We specifically selected 5IJ0 because it captures a soluble /-tubulin heterodimer conformation relevant to microtubule-destabilizing ligands, which can prevent the conformational transitions required for productive lattice incorporation. This dimer was strategically selected due to the distinct expression profiles and functional roles of TUBA1B (1B-tubulin) and TUBB3 (3-tubulin). TUBA1B is frequently used in structural and biochemical studies of microtubules [6,11], while TUBB3 is predominantly expressed in neurons and has been linked to cancer progression, chemoresistance, and altered drug response [6,8,22].
Due to the lack of high-resolution experimental structures for CPPF–tubulin complexes, we employed a structure-based computational strategy. Deep learning methods for biomolecular structure prediction, such as RoseTTAFold All-Atom (RFAA) [23], enable modeling of assemblies containing proteins and small molecules from amino acid sequence and ligand chemical structure. RFAA extends the RoseTTAFold framework to predict three-dimensional protein–ligand complex conformations [23]. Protenix is a deep learning framework for protein–ligand docking and binding site prediction developed by ByteDance. Protenix is specifically designed to predict how a small molecule binds to a protein. It uses a geometric deep learning approach to model the complex interactions and steric constraints at a protein’s binding pocket, directly outputting the likely 3D binding pose of the ligand [24]. In our study, Protenix and RFAA provide strong capabilities for specifying ligand properties, selecting binding pockets, and customizing residue-level constraints, making them suitable for structure-guided docking tasks.
Building on this, we elucidated the atomic-level binding behavior of CPPF through a structure-based computational approach integrating ligand docking, residue-level interaction analysis, and molecular dynamics (MD) simulations. We first identified potential binding pockets on the 1B/TUBB3 heterodimer using ProteinsPlus (https://proteins.plus [25], followed by protein–ligand complex modeling using Protenix, Umol [26], and RFAA, covering both dimeric and monomeric forms. The predicted binding poses were benchmarked using the known microtubule inhibitor nocodazole for validation.
To further dissect molecular interactions, key interacting residues were identified through spatial contact analysis and the Protein–Ligand Interaction Profiler (PLIP) [27], while Proteins, Interfaces, Structures and Assemblies (PDBePISA) [28] was employed to characterize interface energetics. Finally, we performed all-atom MD simulations in Groningen Machine for Chemical Simulations (GROMACS) [29] to assess the structural stability of the predicted complexes and constructed free energy landscapes (FELs) to probe conformational dynamics. This multi-platform framework enables the structural characterization of CPPF–tubulin interactions, bridging the gap between functional assays and atomic-resolution insight.
The contributions of this work are threefold: (1) methodological, we demonstrate a fit-for-purpose AI-assisted protein–ligand modeling strategy for large, flexible systems such as the tubulin heterodimer, benchmarked against experimental crystal structures with sub-1.0 Å agreement for known tubulin–ligand poses; (2) analytical, we introduce a reproducible, end-to-end Nextflow [30] pipeline integrating pocket analysis, multi-platform pose generation, replicate MD simulations, free-energy landscape construction, and MM-PBSA energetics; with all analysis scripts and trajectories publicly deposited, this workflow provides an open-access template that can be readily applied to other drug–target systems; and (3) scientific, we propose a structural model in which CPPF occupies a composite / interface pocket with dominant contributions from -tubulin residues, spatially distinct from the canonical colchicine-bound pose. Additionally, supplementary MD simulations on the -GTP straight lattice conformation (PDB: 6E7B) suggest this pocket remains accessible across the two major -tubulin states, offering a potential rationale for CPPF’s activity in colchicine-resistant cancer cell lines. Modeling the binding behavior of CPPF offers a structural hypothesis for how it may disrupt microtubule dynamics, helping guide future experimental validation.
Materials and methods
Code and data availability
The computational workflow, scripts, and project documentation are available at github.com/jasperyeoh/integrative-ai-assisted-modeling-of-cppf-tubulin-interactions.
The computational workflow is additionally implemented as a reproducible Nextflow (v26.04.1) [30] pipeline available under workflows/ in the same repository, with conda environment specifications and full documentation for end-to-end reproduction.
All-atom MD production trajectories are available on Hugging Face: huggingface.co/datasets/jasperyeoh2/MD-trajectories-CPPF-tubulin-heterodimer-and-monomers.
Ligand physicochemical analysis and tubulin structure preparation
The cryo-EM structure of the 1B/3-tubulin heterodimer (PDB ID: 5IJ0) was retrieved from the Protein Data Bank as the structural template for all modeling tasks [21]. Structural visualization was performed in PyMOL (https://pymol.org [31], where the -subunit was rendered in pink and the -subunit in blue. Magnesium ions (Mg2+) were shown as purple spheres, and the bound nucleotide ligands, GTP and GDP were depicted in lemon and orange stick representations, respectively. The chemical structure of CPPF [20] was constructed in SMILES format using SwissADME [32], and its physicochemical properties, including lipophilicity, polarity, solubility, and flexibility, were evaluated using radar-based descriptors. Predicted toxicity and off-target activity profiles of CPPF were obtained using ProTox 3.0 [33] and are summarized in S1 Fig. The 2D structure and physicochemical profile of CPPF, along with the tertiary structure of the /-tubulin heterodimer, are shown in Fig 1A–1C.
(A) 2D structure of the small-molecule compound CPPF. (B) Physicochemical profile of CPPF represented as a radar plot, including properties such as polarity, lipophilicity, solubility, and molecular size, computed via SwissADME. (C) Overall tertiary structure of the 1B/3 tubulin heterodimer (PDB ID: 5IJ0), visualized in PyMOL. The -subunit is shown in pink and the -subunit in blue. Magnesium ions (purple spheres), GTP (lemon sticks), and GDP (orange sticks) are included as cofactors in the native structure. (D) Predicted ligand-binding pockets identified by ProteinsPlus. Nine surface-accessible pockets (P0 to P8) are rendered in distinct colours and mapped onto the tubulin heterodimer. Pockets with Drug Scores > 0.7 are considered druggable. P2 exhibited the highest Drug Score (0.9) and is located near the / interface, whereas P0, despite its large volume, is likely occluded in the dimerized state. (See S1 Table for detailed quantitative pocket parameters.).
Binding pocket analysis
To identify potential ligand-binding sites for CPPF on the tubulin heterodimer, we employed ProteinsPlus (https://proteins.plus a structure-based platform that evaluates surface geometry and physicochemical parameters of binding pockets [25]. The software computes properties such as volume, depth, surface curvature, and hydrogen bond donor-to-acceptor ratio, which jointly contribute to the estimation of druggability scores. These structural features, particularly large, enclosed cavities with favorable hydrogen bonding potential, are typically associated with strong ligand-binding affinity. The resulting pockets provided spatial constraints for subsequent docking and validation steps.
Protein-ligand complex structure prediction
We employed multiple structure prediction platforms, including Protenix [24], Umol [26], and RFAA [23], to model the binding conformations of CPPF with the 1B/TUBB3 heterodimer, as well as the individual - and -subunits. For each platform, the SMILES representation of CPPF and the corresponding amino acid sequences of the protein targets were provided as inputs. Protenix and RFAA automatically generated multiple sequence alignments and three-dimensional structural models for protein–ligand complexes.
Nocodazole, a well-characterized microtubule-destabilizing agent with an experimentally determined crystal structure binding mode [34], was used as a positive control to benchmark the modeling pipeline. The predictive accuracy of Protenix was further benchmarked by comparing Protenix-predicted nocodazole–-tubulin complexes against the experimental crystal structure of nocodazole-bound -tubulin (PDB: 5CA1), shown in S2 Fig. The predicted poses showed strong structural agreement with the experimental binding mode, supporting the platform’s reliability for modeling CPPF–tubulin interactions. In standard protein–ligand docking and AI-based modeling workflows, crystallographic waters, ions, and non-target cofactors are routinely simplified or omitted during pose prediction, as their explicit inclusion can vary across platforms; the structural validity of this approach is supported by the nocodazole benchmark described above. While Protenix and RFAA produced reproducible pose geometries across -, -, and /-subunit models, predictions from SwissDock [35] and Chai-1 [36] were notably less stable; SwissDock in particular yielded high root-mean-square deviation (RMSD) values relative to crystallographic references. All output structures were aligned to their respective PDB references (5IJ0 for the /-tubulin dimer, TUBA1B for -tubulin, and TUBB3 for -tubulin) using PyMOL for visual inspection and subsequent structural analysis (S3 Fig).
Key residue identification
Residue-level contact analysis was performed across all predicted CPPF and tubulin complexes using a 4 Å heavy-atom distance threshold to define ligand–protein interactions. Biopython [37] was used to identify contact residues, and PyMOL was employed for visual inspection and polar interaction annotation. Complexes derived from Protenix and RFAA, covering both monomeric and dimeric forms of 1B- and 3-tubulin were included in the ensemble. In addition to distance-based contact identification, we further analyzed representative protein–ligand complexes using the Protein–Ligand Interaction Profiler (PLIP) to detect biochemically meaningful interactions such as hydrogen bonds, – stacking, and salt bridges. This analysis allowed for the annotation of key interacting residues based not only on spatial proximity but also on chemical complementarity. Example PyMOL commands for optional hydrophobicity-colored cavity surfaces are available under scripts/ in the accompanying GitHub repository.
Interface analysis
To evaluate the physical and energetic properties of the CPPF–tubulin binding interface, we used the PDBePISA [28] server to analyze protein–ligand complexes generated by RFAA. Unlike docking methods that infer binding poses, PDBePISA quantifies interactions based on existing structural models. The analysis reports key energetic metrics such as interface area, number of contacting atoms and residues, solvation energy, hydrogen bonding, and salt bridge formation, which together contribute to the Gibbs free energy of dissociation ( G). Models with negative G values are considered stable and likely to form persistent complexes in solution.
Molecular dynamics simulation and free energy analysis
All-atom molecular dynamics simulations were performed using GROMACS 2024.5 [29] to assess the conformational stability of CPPF-bound tubulin systems. For the / heterodimer, production MD was initiated from Protenix pose 1 (Table 1) after aligning the predicted protein–ligand assembly to PDB 5IJ0. The crystallographic nucleotide cofactors (GTP at the -site, GDP at the -site) and the Mg2+ ion present in 5IJ0 were not retained in the production MD topology. This protocol is justified by three considerations: (i) the CPPF pose predicted in this study localizes near the -tubulin N-site at the / interface, while the GTP/GDP cofactors occupy the E-site on the opposite face of -tubulin (~15–20 Å away), making direct steric or electrostatic ligand–nucleotide coupling negligible; (ii) experimental affinity measurements at the overlapping colchicine site indicate that ligand binding to the tubulin dimer is not measurably modulated by the -nucleotide state, suggesting the same insensitivity for CPPF; and (iii) the indirect cofactor effect on backbone conformation is encoded upstream by the Protenix prediction, which was performed with all cofactors explicitly present. Monomer simulations used Protenix-predicted CPPF-bound - or -tubulin complexes aligned to their respective PDB references (TUBA1B or TUBB3). Proteins were parameterized with the AMBER99SB-ILDN force field, and the CPPF ligand was parameterized using GAFF2 with RESP2(0.5) atomic charges derived from a B3LYP/6-31G(d) optimization (Multiwfn + ACPYPE). Each protein–ligand complex was solvated in a cubic box filled with SPC216 water molecules under periodic boundary conditions, and neutralized with sodium counterions. Energy minimization was conducted using the steepest descent algorithm, followed by equilibration under NVT and NPT ensembles at 300 K. Production simulations were carried out for 200 ns for monomer systems (three replicates each for - and -tubulin) and 400 ns for the / heterodimer (three independent replicates), using a 2 fs timestep, particle mesh Ewald (PME) for long-range electrostatics, and LINCS to constrain bond lengths.
Post-simulation analyses included the calculation of RMSD for protein backbone atoms and ligand heavy atoms (least-squares fit on backbone), minimum protein–ligand distance (mindist), and radius of gyration (Rg). Trajectories were first processed under periodic boundary conditions using a two-step procedure (trjconv -pbc nojump followed by trjconv -pbc cluster -center for dimers or -pbc mol -ur compact -center for monomers) to avoid discontinuities that can introduce RMSD/Rg artefacts. Hydrogen bond analysis between CPPF and the protein was performed using gmx hbond-legacy due to compatibility issues of the newer gmx hbond implementation with the ligand typing in these systems.
To quantify whether CPPF remained near the predicted -tubulin pocket in monomer simulations, we computed the minimum distance between CPPF and the canonical pocket residues VAL236, LEU253, and ALA314 during the final 50 ns window (150–200 ns). For each residue, an atom group corresponding to the residue was generated and the per-frame minimum distance to CPPF was computed using gmx mindist; distances were summarized as the mean over frames in the window.
To further explore the conformational landscape sampled during simulation, two-dimensional free energy landscapes were constructed using backbone RMSD and radius of gyration (Rg) as reaction coordinates. The joint probability distribution P(x, y) was converted into a free energy surface using Boltzmann inversion:
where is the Boltzmann constant, T is the temperature (300 K), and P(x, y) is the normalized occupancy. FELs were constructed without biasing potentials. RMSD and Rg were calculated from backbone atoms. All calculations were implemented using in-house Python scripts and GROMACS trajectory tools. For the main FEL comparison, the full 200 ns trajectories from all three replicates were concatenated within each class prior to SHAM to reduce sampling noise, and and FELs were plotted using shared axis limits and a shared free-energy colour scale (kcal/mol; capped at 5 kcal/mol) to enable direct visual comparison.
Supplementary 6E7B Simulations (GTP-bound -tubulin conformational state)
To test whether CPPF retains stable binding in the alternative -tubulin conformational state captured by PDB 6E7B (-GMPCPP, straight microtubule-lattice conformation), we performed an independent set of MD simulations using a protocol designed to mirror the primary 5IJ0 simulations. The CPPF–tubulin complex starting structure was generated by running Protenix with the full 6E7B cofactor complement explicitly present in the input (one GTP at the -site, one GMPCPP/G2P at the -site, two Mg2+ ions, and CPPF), so that the predicted pose reflects the cofactor-conditioned pocket geometry of the lattice conformation. All five Protenix samples passed quality thresholds (pLDDT > 94; ipTM > 0.93); the top-ranked sample (ranking score 0.9478) was selected and aligned to the 6E7B template via Kabsch superposition on common C atoms, yielding a backbone RMSD of 1.82 Å relative to the experimental structure. The predicted 6E7B complex and its structural overlay with the 5IJ0-predicted complex are shown in S4 Fig.
The aligned protein–CPPF complex was then prepared for production MD using the identical toolchain as 5IJ0: AMBER99SB-ILDN for the protein, GAFF2/RESP2 for CPPF (re-using the topology generated for 5IJ0, as the ligand is identical), TIP3P water in a rhombic-dodecahedral box (1.0 nm minimum distance to box edge), neutralizing Na+ ions, and a target physiological ionic strength of 0.15 M. As in the 5IJ0 protocol, the bound cofactors were not retained in the production topology (see rationale above). The system contained approximately 168,000 atoms in total. Energy minimization, NVT equilibration (100 ps, 300 K), and NPT equilibration (100 ps, 300 K, 1 bar) used the same MDP parameter files as 5IJ0. Three independent replicates of 200 ns each were performed. Replicates 2 and 3 were initialized from independent NVT/NPT equilibration runs (different velocity seeds) starting from the common minimized structure, ensuring trajectory independence while preserving the equilibrated geometry.
The same analysis pipeline applied to the 5IJ0 trajectories (backbone RMSD, ligand RMSD, radius of gyration, minimum CPPF–protein distance, hydrogen-bond counts, 2D free energy landscapes via Boltzmann inversion of the (RMSD, Rg) joint distribution, and MM-PBSA-GB binding free energy using gmx_MMPBSA with the same GB-OBC2 settings as the primary 5IJ0 simulations (igb = 5, intdiel = 1.0, extdiel = 78.5, ff99SB+GAFF, 1 ns sampling over the last 50 ns)) was applied to the 6E7B trajectories, with shared axis limits and shared free-energy colour scales used to enable direct cross-system visual comparison.
Results
Identification and druggability analysis of binding pockets
To identify potential ligand-binding sites on the tubulin heterodimer, we employed ProteinsPlus (https://proteins.plus which evaluates surface topography and chemical properties to compute pocket-specific druggability scores. A total of 29 surface pockets were detected on the /-tubulin dimer (PDB ID: 5IJ0), with quantitative metrics including volume, surface area, depth, polarity, and hydrophobicity for each pocket (S1 Table). To prioritize biologically relevant pockets, we applied a filtering threshold of Drug Score > 0.7, a commonly used cutoff in the DoGSiteScorer druggability-scoring scheme [38] implemented within ProteinsPlus to highlight potentially druggable pockets for downstream inspection. This yielded nine candidate pockets (P0–P8). Among these, P2 exhibited the highest Drug Score (0.9), with a moderate volume (559.2 Å3) and a favorable depth (36.79 Å), indicating strong potential to accommodate small-molecule ligands. P1 and P7 also displayed high scores (0.83 and 0.84, respectively) and were located on accessible protein surfaces. In contrast, P0, despite its large volume (2445.4 Å3), was positioned at the / interface and is likely inaccessible under dimerized conditions. This integrative analysis, which combines geometric and physicochemical features, provided a rational basis for screening structurally viable pockets for subsequent CPPF docking simulations. The spatial distribution of the identified binding pockets is visualized in Fig 1D. For the /-tubulin dimer, Protenix identified five candidate poses, with pose 1 selected as having the most favorable interaction geometry (Table 1). In this pose, CPPF inserts into a cleft between the - and -subunits, with its pyridinyl nitrogen forming a polar contact with VAL236 on -tubulin at a distance of 2.9 Å (Fig 2A, bottom). This interaction was corroborated by PLIP analysis, which identified a hydrogen bond between the same CPPF atom and VAL236 across multiple poses (distance range: 2.05–2.95 Å), reinforcing the plausibility of this anchoring site (S2 Table).
The predicted binding orientation between CPPF and the /-tubulin dimer was simulated by Protenix (A). The predicted binding orientation between CPPF and -tubulin was predicted by Protenix (B) and RFAA (C). The predicted binding orientation between CPPF and -tubulin was predicted by Protenix (D) and RFAA (E). Zoom-in views below each overview show semi-transparent surfaces with representative interaction distances. Pocket localization and druggability screening are summarized in Figs 1D and 3, respectively.
Structure prediction and validation
In -tubulin-only models, CPPF binding was more variable across platforms. Protenix placed the ligand near ASN257, forming a hydrogen bond at 3.0 Å (Fig 2B), while RFAA predicted an adjacent but more surface-exposed binding mode (Fig 2C) with a hydrogen bond to ALA270 and halogen bonds to SER241 (S2 Table). These differences may reflect the lower structural constraints in the monomeric -tubulin models. Conversely, -tubulin predictions were more consistent across tools. In both Protenix and RFAA models, CPPF occupied a similar hydrophobic pocket involving residues TYR200, GLU198, LEU253, and VAL236, forming recurrent hydrophobic and hydrogen bonding interactions (Fig 2D–2E). In addition, Umol independently predicted a comparable CPPF binding mode on the -tubulin monomer (S5A Fig), centering around residue VAL236 also recurrent in Protenix and RFAA results, thereby reinforcing the consistency and robustness of this predicted anchoring site across platforms. Taken together, these consistent interaction profiles across Protenix, RFAA, and Umol strengthen the confidence in the -tubulin pocket as a stable and conserved binding site for CPPF.
To assess the reproducibility and convergence of predicted CPPF binding poses across oligomeric states, we compared five Protenix-generated models for the /-tubulin dimer, -tubulin monomer, and -tubulin monomer. Monomeric predictions consistently yielded lower RMSD values (0.578–1.028Å) than dimeric counterparts (1.318–1.950Å), as shown in Table 1 and S5B–S5D Figs. For example, pose 4 achieved the best alignment in the -tubulin monomer (0.578Å), while the corresponding dimer model showed 1.869Å. Averaged across all poses, -tubulin models exhibited the highest structural convergence (mean RMSD = 0.682Å), followed by -tubulin (0.949Å), and dimers (1.766Å).
Confidence metrics (pTM and ipTM) were high across all monomeric models, with -tubulin showing slightly higher average scores (pTM = 0.9726, ipTM = 0.9835) than -tubulin (pTM = 0.9696, ipTM = 0.9814). This may reflect the more regular backbone architecture of -tubulin, as AI predictors tend to assign higher confidence to structurally conserved domains. Nevertheless, the tighter fit and lower RMSD values observed in -tubulin models suggest a more constrained and reproducible binding pocket geometry for CPPF.
When focusing specifically on monomeric subunits, -tubulin consistently outperformed -tubulin in pose reproducibility. RMSD values for -tubulin ranged from 0.578 to 0.841 Å, compared to 0.872 to 1.028 Å for -tubulin. Pose 4 was optimal for -tubulin (0.578 Å), while pose 2 yielded the lowest RMSD for -tubulin (0.872 Å). These findings indicate that although -tubulin predictions were marginally more confident overall, 3-tubulin provided a more consistent spatial environment for CPPF docking across replicate models.
To further evaluate pose-level agreement, cross-platform structural overlays were generated. As illustrated in S5 Fig, CPPF binding orientations were topologically conserved across different models, particularly in -tubulin monomer and /-tubulin dimer contexts. Coloring by pose ID (e.g., pale green for pose 0, purple for pose 1, etc.) revealed consistent pocket occupation, similar ligand orientations, and recurrent contact residues. These convergent features support the view that the 3-tubulin cleft represents a structurally consistent anchoring site for CPPF.
Key residue contributions to ligand binding
In -tubulin (TUBB3), PLIP consistently identified recurrent hydrophobic interactions involving VAL236, LEU253, ALA314, and VAL316, along with hydrogen bonds formed with polar residues such as TYR200, ASN256, and VAL236. Additionally, halogen bonding between CPPF and LEU135 was observed in select models, suggesting a role for CPPF’s halogen substituents in enhancing complex stability. In contrast, for -tubulin (TUBA1B), PLIP detected hydrophobic interactions with residues including LEU241, LEU247, PHE254, and LEU317, and hydrogen bonds with ASN248 and ASN257 (S2 Table). However, these -tubulin interactions appeared less frequently across predicted models (typically in 1–4 out of 10 predictions), and were more dispersed across surface-exposed helices and loops, indicating a more variable and less confined binding environment compared to -tubulin. Mapping the contact probabilities onto the tubulin sequence revealed two primary clusters of CPPF-binding residues (Fig 3A). The first cluster, on 1B-tubulin (TUBA1B), included ALA180, SER237, and THR257, each with moderate contact probabilities of 0.36–0.45, appearing in approximately 4–5 out of 10 evaluated models. These residues are located on a solvent-accessible surface and were associated with fewer hydrophobic interactions in PLIP analysis, indicating a weaker and more variable binding mode.
All panels visualize predicted protein–ligand contacts derived from structural models using a 4 Å distance cutoff between heavy atoms of CPPF and tubulin residues. (A) Per-residue contact probability distribution along tubulin sequences. Bar plot summarizing the predicted contact probabilities of individual residues within TUBA1B (blue) and TUBB3 (red). Contact probabilities were computed across a panel of modeled /-tubulin–CPPF complexes. Peaks highlight residues with increased likelihood of interaction, with annotations indicating residue index and identity. The vertical dashed line separates - and -chain data. (B) 3D spatial distribution of predicted contact residues on the tubulin dimer. Three surface-rendered views (front, left, and plan) of the predicted TUBA1B–TUBB3 dimer, showing residues colored by contact probability bin: very low (0.08–0.09, green), low (0.09–0.36, blue), moderate (0.36–0.45, purple), and high (0.45–0.92, red). Only residues predicted to contact CPPF are colored; remaining protein regions are shown in transparent gray. This spatial mapping reveals clustering of interaction hotspots at the inter-chain interface.
In contrast, the second cluster on 3-tubulin (TUBB3) comprised VAL236, LEU253, and ALA314, which showed higher contact probabilities ranging from 0.75 to 0.92 and were present in 7 of the 10 models analyzed. These residues formed a contiguous patch at the / dimer interface. PLIP analysis supported their functional relevance by identifying consistent hydrophobic and hydrogen bonding interactions, reinforcing their role as anchoring hotspots for CPPF. This region is partially buried and topologically confined, supporting a stable binding configuration.
Three-dimensional projections (Fig 3B) revealed that high-contact -tubulin residues clustered within an interfacial groove between - and -subunits, whereas -tubulin contacts were more dispersed. Collectively, this multi-resolution analysis delineates a structurally and chemically coherent CPPF-binding interface enriched on -tubulin, consistent across modeling platforms and interaction profiling.
Structural comparison with the canonical colchicine binding site
To contextualize the predicted CPPF binding mode relative to the canonical colchicine site, we performed structural overlays between the experimentally determined colchicine-site crystal structure of /-tubulin (PDB: 4O2B) [39] and Protenix-predicted CPPF–tubulin models (S6 Fig). The overall structural alignment revealed that CPPF occupies a region proximal to the colchicine binding site at the / interface (S6A Fig). However, a close-up comparison of the binding pocket (S6B Fig) suggested that the predicted CPPF pose is positioned in a deeper -tubulin cleft that is spatially distinct from the crystallographic colchicine-site-bound conformation. This geometric divergence was observed across both the heterodimer and monomer CPPF models.
Interface stability and energetic evaluation
PDBePISA analysis of RFAA-predicted CPPF–tubulin complexes revealed compact but energetically favorable binding interfaces. In the CPPF–TUBB3 complex, 16 protein atoms (0.5%) and 11 residues (2.4%) contributed to an interface area of 183.9 Å2, while all atoms of CPPF were engaged. The CPPF–TUBA1B complex showed a slightly smaller interface (171.2 Å2) with similar protein residue involvement. In both models, the ligand remained largely solvent-exposed (97.8% and 99.9% in TUBB3 and TUBA1B, respectively), suggesting peripheral binding rather than deep cavity insertion (Table 2).
Thermodynamic analysis indicated weak but favorable interactions, with Gibbs free energy of dissociation ( G) estimated at –0.4 kcal/mol (TUBB3) and –0.2 kcal/mol (TUBA1B). Despite minimal protein surface burial, the consistent full engagement of the ligand in both complexes supports the presence of a structurally coherent but potentially transient interface.
At the residue level, PDBePISA identified LEU253 in TUBB3 as a key contact residue, forming a hydrophobic interaction with CPPF at 1.39 Å. In the TUBA1B complex, multiple residues including THR239 and ILE238 were involved, with contact distances ranging from 1.10 to 1.81 Å. These contacts predominantly involved nonpolar atoms and lacked evidence of hydrogen bonds or salt bridges, indicating a hydrophobic interface without significant polar stabilization.
Molecular dynamics and binding stability
We next assessed binding stability using multi-replicate MD simulations, focusing on the assembled / heterodimer as the physiological binding context. Across three independent heterodimer replicates (400 ns each), CPPF maintained stable interfacial contact distances throughout the trajectories (Fig 4C), with minimum protein–ligand distances consistently in the 0.14–0.24 nm range. Backbone RMSD and Rg remained in physically reasonable ranges after periodic-boundary processing, indicating no spurious jump artefacts (Fig 4A–4B). Notably, the three replicates exhibited different convergence behaviors in backbone RMSD: replicate 1 reached a stable plateau (0.3 nm) after approximately 50 ns, while replicates 2 and 3 showed continued conformational drift reaching 0.55 and 0.65 nm by 400 ns, suggesting ongoing structural exploration at this timescale. Importantly, the stable ligand contact distances across all three replicates indicate that the conformational variability primarily reflects backbone rearrangements rather than loss of the binding interface. Summary statistics for the final analysis windows are provided in S3 Table. To provide a complementary view of binding dynamics in monomeric contexts, we compared - and -tubulin monomer simulations (200 ns; three replicates each). Ligand RMSD traces revealed substantially greater variability and replicate-to-replicate spread in the -tubulin monomer than in the -tubulin monomer, while minimum distance traces indicated more consistent proximity in the systems (Fig 5). Summary boxplots for the final 50 ns quantify these monomeric differences (S7 Fig). Notably, residue-level analysis during the final 50 ns (150–200 ns) indicated that CPPF did not consistently maintain proximity to the canonical -tubulin pocket residues VAL236, LEU253, and ALA314 in isolated monomer simulations: mean minimum distances to LEU253 and ALA314 ranged from 0.64–1.36 nm and 0.90–1.59 nm across replicates, respectively, consistent with partial disengagement from the predicted binding site in the monomeric state. In the -tubulin monomer, replicate 1 showed a marked ligand RMSD transition near 80 ns, consistent with a ligand relocation event and underscoring the lower binding-site confinement in the monomeric context.
Panels show: (A) backbone RMSD, (B) protein radius of gyration (Rg), (C) minimum protein–ligand distance, and (D) protein–ligand hydrogen bond count, each plotted as three replicate traces with an across-replicate min–max band.
Panels show ligand RMSD for (A) -tubulin and (B) -tubulin monomers, and minimum protein–ligand distance for (C) -tubulin and (D) -tubulin monomers.
To further examine the conformational space sampled during simulation, we constructed two-dimensional free energy landscapes (FELs) using backbone RMSD and Rg as collective variables. The full 200 ns trajectories from all three replicates were concatenated within each class ( and ) prior to FEL construction to improve sampling coverage. Consistent with the time series analyses, the combined landscape exhibited a deeper and broader low-energy basin than the combined landscape (Fig 6), indicating a larger ensemble of energetically favorable conformations in the binding context. The global minimum of the combined landscape is located near Rg 2.19 nm and RMSD 0.43 nm, whereas the combined landscape minimum occurs near Rg 2.20 nm and RMSD 0.48 nm, reflecting that the most probable conformations remain more structurally drifted. Both landscapes are plotted with shared axis limits and a shared free-energy colour scale (kcal/mol; capped at 5 kcal/mol) to enable direct visual comparison.
The full 200 ns trajectories from all three replicates were concatenated within each class ( and ) to improve sampling. Panels show the -tubulin landscape as (A) a 3D surface and (B) a 2D contour map, and the -tubulin landscape as (C) a 3D surface and (D) a 2D contour map. and landscapes are plotted with shared axis limits and a shared free-energy colour scale (kcal/mol; capped at 5 kcal/mol) to enable direct visual comparison of basin depth and localization.
Binding free energy was estimated using the MM-PBSA (GB) method applied to the heterodimer trajectories over the final 50 ns (350–400 ns, 51 frames at 1 ns intervals). Across the three replicates, the binding free energy components were: VDW = −41.96 1.33 kcal/mol, EEL = −5.91 2.22 kcal/mol, EGB = +22.29 0.63 kcal/mol, ESURF = −5.61 0.23 kcal/mol, and TOTAL = −31.19 4.04 kcal/mol (mean SD across replicate averages). Bond energy terms that exhibited numerical overflow in individual sander outputs were flagged and excluded from per-frame decomposition; they do not affect TOTAL because bond contributions are identical in the complex and receptor calculations and cancel upon subtraction. These values indicate a substantially more favorable binding interaction than the static PDBePISA estimate (G = −0.4 kcal/mol), which reflects only the interface geometry of the predicted pose without accounting for conformational sampling or solvation dynamics.
CPPF binding stability in the 6E7B (GTP-bound -tubulin) conformational state
To examine whether CPPF binding is preserved across the major -tubulin conformational states, we performed three independent 200 ns MD replicates on the 6E7B template, in which -tubulin adopts the GTP-bound straight microtubule-lattice conformation (Methods and S8–S10 Figs). The protein–CPPF MD protocol was identical to that used for the 5IJ0 simulations to enable direct cross-system comparison.
Across all three 6E7B replicates, CPPF remained stably bound throughout the trajectories (S8 Fig). Over the final 50 ns analysis window, the backbone RMSD across replicates was nm (per-replicate means: 0.270, 0.295, 0.321 nm), which is comparable to the lowest-drift 5IJ0 replicate (rep1, 0.30 nm plateau) and substantially lower than 5IJ0 rep2/rep3 (0.55 and 0.65 nm). The tighter convergence is consistent with the more constrained backbone geometry of the straight lattice conformation relative to the curved soluble dimer. The minimum CPPF–protein distance was nm across replicates (per-replicate means: 0.209, 0.214, 0.188 nm), within the 5IJ0 range of 0.14–0.24 nm, indicating that CPPF maintains continuous van-der-Waals contact with the binding pocket in the 6E7B conformation. The 2D free energy landscape constructed by Boltzmann inversion of the concatenated (RMSD, Rg) joint distribution exhibits a single well-defined global minimum at RMSD 0.27 nm, Rg 2.99 nm (S9 Fig). The larger Rg compared with the 5IJ0 minimum (Rg 2.19 nm) reflects the extended geometry of the straight lattice conformation rather than loss of binding; basin localization is comparable to that of 5IJ0 , indicating equivalent conformational confinement of the CPPF-bound state. MM-PBSA-GB analysis (gmx_MMPBSA, GB-OBC2, last 50 ns at 1 ns sampling, identical protocol to the 5IJ0 simulations) yielded an across-replicate mean binding free energy of kcal/mol (per-replicate: , , kcal/mol), which is comparable to the 5IJ0 value of kcal/mol, with overlapping 1 intervals (S10 Fig).
Together, these supplementary findings indicate that CPPF accommodation is robust to the -tubulin conformational state. The ligand maintains its position at the /-interface with high stability. Its backbone RMSD profile shows less drift than two of the three 5IJ0 replicates, while all other metrics remain comparable. This is consistent with the biochemical expectation that CPPF and the nucleotide cofactors occupy spatially distinct sites at the -tubulin N-site and E-site, respectively, and with published experimental evidence that ligand binding at the overlapping colchicine site is not appreciably modulated by the -tubulin nucleotide state.
Discussion
This study models the predicted recognition site of /-tubulin dimer by CPPF through a structure-guided framework integrating AI-powered structure prediction (Protenix, RFAA, Umol), spatial contact profiling, atomic-resolution interaction annotation (PLIP), and molecular dynamics (MD) simulations. The analysis reveals that CPPF targets a composite / interfacial pocket dominated by 3-tubulin contacts, supported by comparative modeling across /-tubulin dimers and monomers (Figs 2 and S5 and S2 Table). The primary simulations employed the soluble, curved conformation of tubulin (PDB: 5IJ0), which represents the state relevant to microtubule-destabilizing ligands. The supplementary 6E7B simulations (3 200 ns) further indicate that CPPF binding is preserved in the GTP-bound straight microtubule-lattice conformation, with binding-pocket stability metrics comparable to the 5IJ0 replicates and a comparable MM-PBSA binding free energy ( vs kcal/mol, overlapping 1 intervals), supporting that CPPF accommodation is robust to the -tubulin conformational state. Comparative modeling reveals that monomeric 3-tubulin exhibits superior structural convergence compared to 1B-tubulin and dimers. RMSD values for 3-tubulin range from 0.578 to 0.841 Å, lower than 1B-tubulin (0.872–1.028 Å) and dimers (1.318–1.950 Å) (Table 1), despite -tubulin’s slightly higher pTM (0.9726) and ipTM (0.9835) scores versus -tubulin’s (0.9696, 0.9814). PDBePISA analysis further supports this preference, indicating a larger interface area (183.9 Å2 vs. 171.2 Å2), greater atom engagement (0.5% vs. 0.4%), and a more favorable binding energy (G = -0.4 vs. -0.2 kcal/mol) for 3-tubulin (Table 2), with tighter hydrophobic anchoring at LEU253 (1.39 Å) enhancing van der Waals complementarity. These findings indicate a geometrically optimized 3-tubulin site with higher structural consistency. Interaction profiling identifies a predicted -tubulin binding pocket featuring VAL236, LEU253, and ALA314, with high contact probabilities (0.75–0.92) across models (Fig 3A), as supported by Protenix, RFAA, Umol, and PLIP analysis (S2 Table). These residues contribute to favorable hydrophobic and hydrogen-bond interactions (e.g., VAL316, TYR200, ASN256), including a halogen bond with LEU135 in select 3-tubulin complexes. In contrast, the 1B-tubulin region (ALA180, SER237, THR257) shows moderate contact probabilities (0.36–0.45) and occupies more solvent-exposed topographies with reduced polar interactions, suggesting a less stable binding environment. MD simulations further support that the assembled / heterodimer provides the relevant structural framework for stable CPPF accommodation. In isolated -tubulin monomer simulations (200 ns), CPPF did not consistently maintain proximity to the canonical pocket residues over the trajectory. During the final 50 ns, the mean minimum distances to LEU253 and ALA314 ranged from 0.64–1.36 nm and 0.90–1.59 nm across replicates, respectively, indicating partial disengagement from the predicted binding site in the monomeric state. In contrast, within the /-tubulin heterodimer (400 ns; three independent replicates), CPPF maintained stable interfacial contacts throughout the simulations with consistent minimum protein–ligand distances of 0.14–0.24 nm (Fig 4), and a mean MM-PBSA binding free energy of kcal/mol. The deeper and broader free energy basin observed for -tubulin in the FEL analysis further supports that heterodimer formation and the resulting interfacial geometry are necessary for stable CPPF accommodation (Fig 6). These findings suggest that the biological context of the assembled dimer, rather than the isolated subunit, is the relevant structural framework for CPPF binding.
Prior MM-GBSA studies of colchicine-site ligands bound to human -tubulin heterodimers report endpoint binding free energies on the order of tens of kilocalories per mole under implicit-solvent end-state protocols (for example, to kcal/mol for DAMA-colchicine across tubulin isotypes [40]). Absolute magnitudes are not directly comparable across studies, as they depend on the choice of ligand, force field, implicit solvent model (GB versus PB), trajectory length, and entropy approximation. Our ensemble-averaged MM-PBSA TOTAL for CPPF ( kcal/mol) is of the same order of magnitude as, though quantitatively lower than, these DAMA-colchicine benchmarks [40]; the offset is consistent with expected differences in ligand chemistry (CPPF is smaller and less buried than DAMA-colchicine) and protocol details, and this comparison is intended to establish scale rather than direct affinity ranking. Our result nevertheless remains substantially more favorable than the PDBePISA single-structure interface estimate that omits conformational sampling [28]. Together with docking evidence placing CPPF at the colchicine-site region [20], these comparisons contextualize the present energetics relative to prior in silico studies of colchicine-class ligands.
The stable interaction between the 3-tubulin and CPPF is consistent with Han et al.’s (2020) prediction that CPPF binds at the colchicine site [20], where it can depolymerize microtubules and suppress growth in multidrug-resistant cancer cell lines, suggesting a mechanism distinct from paclitaxel’s microtubule stabilization [22]. This -tubulin-dominant contact pattern may help explain CPPF activity in multidrug-resistant cellular contexts, including TUBB3-enriched tumors such as glioblastoma models, where TUBB3 overexpression is associated with chemoresistance [8]. The PDBePISA static estimate (G = -0.4 kcal/mol) reflects only the interface geometry of a single predicted pose; MM-PBSA calculations incorporating conformational sampling yield a substantially more favorable TOTAL of -31.19 4.04 kcal/mol, supporting a stable binding interaction that warrants further experimental validation. Notably, the predicted CPPF binding pocket diverges from vinblastine’s lateral interface engagement, suggesting a unique tubulin interaction mode that may mitigate resistance mechanisms associated with traditional microtubule inhibitors [20]. Given that CPPF occupies a -tubulin cleft partially overlapping with the colchicine site and exhibits potential stabilizing contacts at the / dimer interface (inferred from PDBePISA interface area, Table 2), it is plausible that such binding may interfere with protofilament assembly or destabilize lateral contacts critical for microtubule polymerization. This contrasts with paclitaxel, which binds along the microtubule lumen. These distinctions raise the possibility that CPPF may reverse paclitaxel resistance or enhance synergistic effects in combination therapies. However, functional assays, such as microtubule polymerization kinetics, cell-based synergy tests, or MDR reversal studies, are required to validate these hypotheses. The structural overlay with the crystallized colchicine-site crystal structure of /-tubulin (PDB: 4O2B) [39] further refines this interpretation (S6 Fig). Although CPPF is positioned near the broader colchicine-site region at the / interface, the Protenix-predicted CPPF pose occupies a deeper -tubulin cleft that is spatially distinct from the crystallographic colchicine-site-bound conformation. This divergence may reconcile an apparent discrepancy in the literature: although prior AutoDock modeling suggested partial overlap between CPPF and colchicine binding sites [20], CPPF remains effective in colchicine-resistant cancer cell contexts [20]. The distinct binding geometry identified here suggests that CPPF may engage -tubulin residues outside the canonical colchicine-bound pose, potentially helping to circumvent colchicine-specific resistance mechanisms. This interpretation is consistent with the observed MDR-bypassing activity of CPPF and supports a binding mode related to, but functionally distinct from, classical colchicine-site inhibitors.
The multi-conformational compatibility of the CPPF binding pocket has a further mechanistic implication. Whereas paclitaxel binds a lateral inter-protofilament interface that is only fully formed once tubulin has polymerized into the microtubule lattice [11,12], and colchicine-site occupancy is coupled to the curved-to-straight transition of the /-heterodimer [16,17], the interfacial site occupied by CPPF is present in both the soluble curved dimer (5IJ0) and the straight lattice-related conformation (6E7B). CPPF retains stable binding in both conformational states across the four independent stability metrics reported here: backbone RMSD, minimum protein–ligand distance, FEL basin localization, and MM-PBSA binding free energy ( vs kcal/mol). This cross-state accessibility supports the hypothesis that CPPF may engage tubulin at multiple stages of microtubule dynamics, potentially both antagonizing curved-to-straight polymerization from soluble dimers and destabilizing already-assembled lattice; definitive mechanistic distinction from stage-restricted binders would require dynamic assembly assays or single-molecule imaging.
While these analyses provide useful structural insights, limitations include the finite MD timescale (200 ns monomers; 400 ns heterodimers; 200 ns 6E7B supplementary replicates) and the intrinsic uncertainty of AI-based structure prediction and docking in capturing long-timescale conformational transitions [41,42]. The observed replicate-to-replicate variability in backbone drift in the heterodimer trajectories further indicates ongoing structural exploration at the 400 ns timescale, even when the binding interface remains stable. These discrepancies necessitate experimental validation. Although the supplementary 6E7B simulations indicate that CPPF binding is preserved in the GTP-bound straight -tubulin conformation, a fully explicit GMPCPP-parameterized comparison, which would resolve any direct chemical contribution of the bound nucleotide at the E-site, remains a refinement for follow-up work. In summary, this study presents computational evidence that CPPF engages a composite / interface pocket characterized by primary contacts with 3-tubulin. The enhanced structural convergence, favorable energetics, and compact binding interface of 3-tubulin suggest that it contributes substantially to CPPF recognition and its potential microtubule-depolymerizing and MDR-bypassing effects. These computational predictions provide useful reference points for future analog design and warrant further validation through experimental assays.
Supporting information
S1 Fig. Toxicity and off-target activity profile of CPPF predicted by Pro-Tox 3.0.
Radar plot showing the predicted toxicity probabilities of CPPF across a panel of toxicity endpoints and off-target targets. The blue line with filled area represents the toxicity probabilities of CPPF, while the orange line with filled area shows the average activity probabilities of known molecules from the ProTox 3.0 training dataset.
https://doi.org/10.1371/journal.pcbi.1014804.s001
(TIF)
S2 Fig. Structural alignment and benchmark validation of Protenix-predicted -tubulin–nocodazole complexes relative to the experimental crystal structure (PDB: 5CA1).
Structural alignment was performed between the experimentally determined crystal structure of nocodazole-bound -tubulin (PDB: 5CA1) and corresponding computational models predicted by Protenix. (A) Overview of the experimental -tubulin–nocodazole structure overlaid with Protenix-predicted complexes (poses 0–2). (B) Close-up view of the ligand–receptor binding interface comparing the experimental structure and predicted poses. The crystal structure of -tubulin in complex with nocodazole was derived from the T2R-TTL-nocodazole assembly (PDB: 5CA1). Cofactors and non-target protein components were omitted in PyMOL to enhance visual clarity. Structural alignments were conducted independently between PDB 5CA1 and each of the Protenix-predicted poses 0, 1, and 2. Color coding: experimental -tubulin–nocodazole, magenta; predicted pose 0, cyan; predicted pose 1, blue; predicted pose 2, purple.
https://doi.org/10.1371/journal.pcbi.1014804.s002
(TIF)
S3 Fig. Structural alignment of predicted CPPF–tubulin complexes with experimental PDB structures.
Predicted CPPF–tubulin complexes were structurally aligned to experimentally determined tubulin structures from the Protein Data Bank (PDB). (A) Alignment of the /-tubulin dimer (PDB ID: 5IJ0) with the CPPF–/-tubulin dimer predicted by Protenix. (B) Alignment of -tubulin (TUBA1B) from the PDB with CPPF-bound -tubulin predicted by Protenix. (C) Alignment of -tubulin predicted by RFAA with the experimental structure. (D) Overlay of the experimental -tubulin structure with Protenix- and RFAA-predicted CPPF-bound models. (E) Alignment of -tubulin (TUBB3) from the PDB with the CPPF–-tubulin complex predicted by Protenix. (F) Alignment of -tubulin predicted by RFAA with the experimental structure. (G) Overlay of the experimental -tubulin structure with Protenix- and RFAA-predicted CPPF-bound models. For all alignments, the best representative Protenix pose (pose 1) was selected. Color scheme: PDB structures (pink for -tubulin, blue for -tubulin), Protenix predictions (purple), and RFAA predictions (yellow). RMSD values are reported in Table 1. All structural alignments and visualizations were performed using PyMOL.
https://doi.org/10.1371/journal.pcbi.1014804.s003
(TIF)
S4 Fig. Structural comparison of the Protenix-predicted CPPF–6E7B complex with the 5IJ0-predicted complex.
(A) Overall view of the Protenix-predicted CPPF–/-tubulin complex generated from the 6E7B (-GMPCPP, straight microtubule-lattice) template, used as the MD starting structure for the supplementary 6E7B simulations (Methods and S8–S10 Figs). (B) Zoomed-in view of the predicted CPPF-binding pocket in the 6E7B-derived complex, showing CPPF (yellow) adjacent to -tubulin VAL236 (orange; closest heavy-atom distance 2.48 Å) and LEU253 (teal; closest heavy-atom distance 2.80 Å), the same two residues identified as principal CPPF contacts in the 5IJ0 analysis (Fig 2A and S2 Table). (C) Structural overlay of the 6E7B-predicted complex (blue -tubulin backbone; yellow CPPF; orange VAL236/LEU253 side chains) with the 5IJ0-predicted complex (salmon -tubulin backbone; magenta CPPF; purple VAL236/LEU253 side chains), aligned on the -tubulin backbone C atoms (PyMOL align; RMSD = 0.212 Å over 382 atom pairs). The close superposition of both the protein backbone and the two independently predicted CPPF poses indicates that the predicted binding site is conserved between the two -tubulin conformational states, independent of and complementary to the MD/MM-PBSA-based evidence in S8–S10 Figs. Cofactors and non-target chains hidden for visual clarity.
https://doi.org/10.1371/journal.pcbi.1014804.s004
(TIF)
S5 Fig. Predicted CPPF binding poses on tubulin generated by Umol and Protenix.
Predicted binding orientations of CPPF on tubulin were generated using Umol and Protenix. (A) Binding pose of CPPF on -tubulin predicted by Umol. (B–D) Five predicted poses (pose 0–4) of CPPF in complex with the /-tubulin dimer (B), -tubulin monomer (C), and -tubulin monomer (D) generated by Protenix, aligned to their respective experimental PDB structures. Structures were color-coded as follows: PDB -tubulin (pink), PDB -tubulin (blue), Protenix pose 0 (pale green), pose 1 (purple), pose 2 (lavender blue), pose 3 (wheat), and pose 4 (gray). Root-mean-square deviation (RMSD) values for each pose are reported in Table 1. All structural visualizations and alignments were performed using PyMOL.
https://doi.org/10.1371/journal.pcbi.1014804.s005
(TIF)
S6 Fig. Structural comparison of the crystallized colchicine-site crystal structure of /-tubulin (PDB: 4O2B) with Protenix-predicted CPPF-bound -tubulin models.
(A) Overall view of the experimental colchicine-site crystal structure of /-tubulin overlaid with Protenix-predicted CPPF-bound complexes, including the full /-tubulin heterodimer and the isolated -tubulin monomer. (B) Zoomed-in view of the ligand-binding pocket, comparing the experimental colchicine-site-bound complex with the predicted CPPF–-tubulin interface from both the heterodimer and monomer models. The colchicine-site crystal structure of /-tubulin was obtained from the Protein Data Bank (PDB: 4O2B) [39]. Cofactors and non-target protein components were hidden in PyMOL to enhance visual clarity. Structural alignments were performed between PDB 4O2B and the Protenix-predicted CPPF-bound /-tubulin heterodimer and -tubulin monomer models. Color scheme: experimental /-tubulin–colchicine-site-bound complex, pink; Protenix-predicted /-tubulin heterodimer–CPPF complex, yellow; Protenix-predicted -tubulin monomer–CPPF complex, cyan.
https://doi.org/10.1371/journal.pcbi.1014804.s006
(TIF)
S7 Fig. Monomer MD summary metrics over the final 50 ns (boxplots).
Boxplots comparing - and -tubulin monomer simulations (200 ns; three independent replicates each) over the final 50 ns window (150–200 ns). (A) Minimum protein–ligand distance. (B) Mean binding-site RMSF (nm) averaged over residues VAL236, LEU253, and ALA314. (C) Ligand RMSD. (D) Protein radius of gyration (Rg). For each metric, distributions pool all frames within the window per replicate; boxplots summarize versus classes (three replicate values per class for panel B; pooled frame distributions for panels A, C, and D). These summaries complement the full-trajectory time series in the main text (Fig 5) and support greater binding-site variability and weaker late-stage confinement in -tubulin monomer simulations relative to -tubulin.
https://doi.org/10.1371/journal.pcbi.1014804.s007
(TIF)
S8 Fig. CPPF binding stability in the 6E7B (-GTP/microtubule-lattice straight) conformation — MD time series.
Three independent 200 ns MD replicates of the CPPF–-tubulin complex starting from a Protenix-predicted pose aligned to PDB 6E7B (Methods). (A) Backbone RMSD; (B) protein radius of gyration; (C) minimum CPPF–protein distance; (D) CPPF–protein hydrogen-bond count. Per-replicate traces: rep1 blue, rep2 orange, rep3 green; gray shaded band: across-replicate min–max envelope. The MD protocol mirrors the 5IJ0 simulations (AMBER99SB-ILDN + GAFF2/RESP2; TIP3P water; 0.15 M NaCl; 2 fs timestep; cofactors not retained in the production topology). Across all three replicates, last-50-ns mean values are: backbone RMSD = nm, minimum CPPF–protein distance = nm. Corresponding window statistics are provided in S4 Table.
https://doi.org/10.1371/journal.pcbi.1014804.s008
(TIF)
S9 Fig. 6E7B two-dimensional free energy landscape (FEL).
2D FEL constructed by Boltzmann inversion of the concatenated (Rg, backbone RMSD) joint distribution from all three 6E7B replicates. Energies are plotted in kcal/mol and capped at 5 kcal/mol to match the visualization convention of Fig 6 and S11 Fig. Global minimum (black dot): Rg 2.99 nm, backbone RMSD 0.27 nm. The higher Rg relative to the 5IJ0 minimum (Rg 2.19 nm) reflects the extended geometry of the straight lattice conformation rather than loss of binding; basin localization is comparable to that of 5IJ0 in Fig 6, indicating equivalent conformational confinement of the CPPF-bound state.
https://doi.org/10.1371/journal.pcbi.1014804.s009
(TIF)
S10 Fig. MM-PBSA-GB binding free energy comparison between 5IJ0 and 6E7B.
Bar comparison of CPPF binding free energy () between the 5IJ0 heterodimer simulations (3 400 ns; blue) and the 6E7B simulations (3 200 ns; orange). Bar height: across-replicate mean. Error bars: across-replicate s.d. Black dots: per-replicate values. Both systems were analyzed with the identical MM-PBSA-GB protocol (gmx_MMPBSA v1.5 + ; GB-OBC2 with igb = 5, intdiel = 1.0, extdiel = 78.5; ff99SB + GAFF; last 50 ns at 1 ns sampling, 50 snapshots per replicate). Across-replicate means: 5IJ0 = kcal/mol; 6E7B = kcal/mol. The two 1 intervals overlap completely, indicating comparable CPPF binding energetics across the two -tubulin conformational states.
https://doi.org/10.1371/journal.pcbi.1014804.s010
(TIF)
S11 Fig. Supplementary free energy landscapes (Rg–backbone RMSD) for monomer simulations.
Two-dimensional free energy landscapes (FELs) for CPPF–tubulin monomer simulations constructed using backbone RMSD and Rg as collective variables. The full 200 ns trajectories from all three replicates were concatenated within each monomer class prior to SHAM analysis. Landscapes are plotted in kcal/mol with a colour scale capped at 5 kcal/mol (same convention as the 5IJ0 FEL comparison in Fig 6). (A) Combined -tubulin monomer landscape. (B) Combined -tubulin monomer landscape. The -tubulin landscape exhibits a deeper and broader low-energy basin than the -tubulin landscape, consistent with Fig 6.
https://doi.org/10.1371/journal.pcbi.1014804.s011
(TIF)
S1 Table. ProteinsPlus-predicted pocket properties for the /-tubulin heterodimer.
This table reports the complete raw quantitative output from ProteinsPlus for all 29 surface pockets (P0–P28) detected on the /-tubulin dimer (PDB ID: 5IJ0). For each pocket, geometric and physicochemical descriptors are provided, including pocket volume, surface area, depth, surface-to-volume ratio, SimpleScore, DrugScore, hydrogen-bond acceptor/donor counts, hydrophobic interaction counts, metal-ion involvement, and the relative composition of nonpolar/polar/apolar surface areas. Nine pockets (P0–P8) with DrugScore > 0.7 were considered potentially druggable and were prioritized for downstream inspection. The spatial localization of these pockets on the tubulin structure is visualized in Fig 1D.
https://doi.org/10.1371/journal.pcbi.1014804.s012
(PDF)
S2 Table. Residue-level interaction summary between CPPF and /-tubulin across predicted models identified by PLIP.
This table summarizes the predicted non-covalent interactions between the small molecule CPPF and - and -tubulin residues as identified using the Protein–Ligand Interaction Profiler (PLIP). Multiple binding poses were analyzed across different structural models, including Protenix-based docking predictions and RFAA-predicted tubulin structures. For each interaction, the interaction type (hydrophobic interaction, hydrogen bond, or halogen bond), interacting residue (with chain identifier and residue number), and associated geometric descriptors (interatomic distance and interaction vector angles) are reported.
https://doi.org/10.1371/journal.pcbi.1014804.s013
(PDF)
S3 Table. Summary statistics for MD simulation observables over the final analysis window.
Per-system summary of key MD-derived metrics computed over the final 50 ns window (350–400 ns for heterodimer replicates; 150–200 ns for monomer replicates). Metrics include backbone RMSD, ligand RMSD, minimum protein–ligand distance (mindist_pl), hydrogen bond count (hbond_num), radius of gyration (Rg), solvent-accessible surface area (SASA), and per-residue RMSF. Values represent mean SD across frames within the specified window. All distances are in nm and SASA in nm2.
https://doi.org/10.1371/journal.pcbi.1014804.s014
(CSV)
S4 Table. Per-replicate window statistics for the 6E7B MD simulations.
Per-replicate summary of key MD-derived metrics computed over the final 50 ns analysis window (150–200 ns) of the three 6E7B simulations. Metrics reported for each replicate include backbone RMSD, minimum CPPF–protein distance (mindist_pl), hydrogen bond count (hbond_num), radius of gyration (Rg), and solvent-accessible surface area (SASA). Values represent mean SD across frames within the analysis window. All distances are in nm and SASA in nm2.
https://doi.org/10.1371/journal.pcbi.1014804.s015
(CSV)
Acknowledgments
Part of the computational work reported in this study was performed using computing resources provided by the High Performance Computing Platform of The Hong Kong University of Science and Technology (Guangzhou), with support from the College of Future Technology.
References
- 1. Desai A, Mitchison TJ. Microtubule polymerization dynamics. Annu Rev Cell Dev Biol. 1997;13:83–117. pmid:9442869
- 2. Kapitein LC, Hoogenraad CC. Building the neuronal microtubule cytoskeleton. Neuron. 2015;87(3):492–506. pmid:26247859
- 3. Heald R, Khodjakov A. Thirty years of search and capture: The complex simplicity of mitotic spindle assembly. J Cell Biol. 2015;211(6):1103–11. pmid:26668328
- 4. Olmsted JB, Borisy GG. Microtubules. Annu Rev Biochem. 1973;42:507–40.
- 5. Verdier-Pinard P, Wang F, Burd B, Angeletti RH, Horwitz SB, Orr GA. Direct analysis of tubulin expression in cancer cell lines by electrospray ionization mass spectrometry. Biochemistry. 2003;42(41):12019–27. pmid:14556633
- 6. Ludueña RF. Multiple forms of tubulin: Different gene products and covalent modifications. Int Rev Cytol. 1997;178:207–75.
- 7. Zhao D, Deshpande R, Wu K, Tyagi A, Sharma S, Wu S-Y, et al. Identification of TUBB3 as an immunotherapy target in lung cancer by genome wide in vivo CRISPR screening. Neoplasia. 2025;60:101100. pmid:39671912
- 8. Duly AMP, Kao FCL, Teo WS, Kavallaris M. βIII-tubulin gene regulation in health and disease. Front Cell Dev Biol. 2022;10:851542. pmid:35573698
- 9. Jordan MA, Wilson L. Microtubules as a target for anticancer drugs. Nat Rev Cancer. 2004;4(4):253–65. pmid:15057285
- 10. Rowinsky EK, Donehower RC. Paclitaxel (Taxol). N Engl J Med. 1995;332(15):1004–14.
- 11. Nogales E, Wolf SG, Downing KH. Structure of the alpha beta tubulin dimer by electron crystallography. Nature. 1998;391(6663):199–203. pmid:9428769
- 12. Alushin GM, Lander GC, Kellogg EH, Zhang R, Baker D, Nogales E. High-resolution microtubule structures reveal the structural transitions in αβ-tubulin upon GTP hydrolysis. Cell. 2014;157(5):1117–29. pmid:24855948
- 13. Nogales E, Whittaker M, Milligan RA, Downing KH. High-resolution model of the microtubule. Cell. 1999;96(1):79–88. pmid:9989499
- 14. Chrétien D, Metoz F, Verde F, Karsenti E, Wade RH. Lattice defects in microtubules: Protofilament numbers vary within individual microtubules. J Cell Biol. 1992;117(5):1031–40. pmid:1577866
- 15. Wade RH, Chrétien D, Job D. Characterization of microtubule protofilament numbers. How does the surface lattice accommodate?. J Mol Biol. 1990;212(4):775–86. pmid:2329582
- 16. Fedorov VA, Orekhov PS, Kholina EG, Zhmurov AA, Ataullakhanov FI, Kovalenko IB, et al. Mechanical properties of tubulin intra- and inter-dimer interfaces and their implications for microtubule dynamic instability. PLoS Comput Biol. 2019;15(12):e1007327.
- 17. Nogales E, Wang H-W. Structural mechanisms underlying nucleotide-dependent self-assembly of tubulin and its relatives. Curr Opin Struct Biol. 2006;16(2):221–9. pmid:16549346
- 18. Nogales E, Wang HW, Niederstrasser H. Tubulin Rings: Which Way Do They Curve? Curr Opin Struct Biol. 2003;13(3):256–61.
- 19. Brouhard GJ, Rice LM. The contribution of αβ-tubulin curvature to microtubule dynamics. J Cell Biol. 2014;207(3):323–34. pmid:25385183
- 20. Han HJ, Park C, Hwang J, Thimmegowda NR, Kim SO, Han J, et al. CPPF, a novel microtubule targeting anticancer agent, inhibits the growth of a wide variety of cancers. Int J Mol Sci. 2020;21(13):4800.
- 21. Ti S-C, Pamula MC, Howes SC, Duellberg C, Cade NI, Kleiner RE, et al. Mutations in human tubulin proximal to the kinesin-binding site alter dynamic instability at microtubule plus- and minus-ends. Dev Cell. 2016;37(1):72–84. pmid:27046833
- 22. Banerjee A. Increased levels of tyrosinated-,III-, andIV-tubulin isotypes in paclitaxel-resistant MCF-7 breast cancer cells. Biochem Biophys Res Commun. 2002;293(2):598–601.
- 23. Krishna R, Wang J, Ahern W, Sturmfels P, Venkatesh P, Kalvet I, et al. Generalized biomolecular modeling and design with RoseTTAFold All-Atom. Science. 2024;384(6693):eadl2528. pmid:38452047
- 24. ByteDance AML AI4Science Team, Chen X, Zhang Y, Lu C, Ma W, Guan J, et al. Protenix: Advancing structure prediction through a comprehensive AlphaFold3 reproduction. bioRxiv. 2025;2025.01.08.631967.
- 25. Schöning-Stierand K, Diedrich K, Ehrt C, Flachsenberg F, Graef J, Sieg J, et al. ProteinsPlus: A comprehensive collection of web-based molecular modeling tools. Nucleic Acids Res. 2022;50(W1):W611–5. pmid:35489057
- 26. Bryant P, Kelkar A, Guljas A, Clementi C, Noé F. Structure prediction of protein-ligand complexes from sequence information with Umol. Nat Commun. 2024;15(1):4536. pmid:38806453
- 27. Adasme MF, Linnemann KL, Bolz SN, Kaiser F, Salentin S, Haupt VJ, et al. PLIP 2021: Expanding the scope of the protein-ligand interaction profiler to DNA and RNA. Nucleic Acids Res. 2021;49(W1):W530–4. pmid:33950214
- 28. Krissinel E, Henrick K. Inference of macromolecular assemblies from crystalline state. J Mol Biol. 2007;372(3):774–97. pmid:17681537
- 29. Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, et al. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX. 2015;1–2:19–25.
- 30. Di Tommaso P, Chatzou M, Floden EW, Barja PP, Palumbo E, Notredame C. Nextflow enables reproducible computational workflows. Nat Biotechnol. 2017;35(4):316–9. pmid:28398311
- 31.
Schrödinger LLC. The PyMOL Molecular graphics system, Version 1.8. Schrödinger: LLC. 2015. 2025. https://pymol.org/
- 32. Daina A, Michielin O, Zoete V. SwissADME: A free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Sci Rep. 2017;7:42717. pmid:28256516
- 33. Banerjee P, Kemmler E, Dunkel M, Preissner R. ProTox 3.0: a webserver for the prediction of toxicity of chemicals. Nucleic Acids Res. 2024;52(W1):W513–20. pmid:38647086
- 34. Wang Y, Zhang H, Gigant B, Yu Y, Wu Y, Chen X, et al. Structures of a diverse set of colchicine binding site inhibitors in complex with tubulin provide a rationale for drug discovery. FEBS J. 2016;283(1):102–11. pmid:26462166
- 35. Bugnon M, Röhrig UF, Goullieux M, Perez MAS, Daina A, Michielin O, et al. SwissDock 2024: major enhancements for small-molecule docking with Attracting Cavities and AutoDock Vina. Nucleic Acids Res. 2024;52(W1):W324–32. pmid:38686803
- 36. Chai Discovery Team, Boitreaud J, Dent J, McPartlon M, Meier J, Reis V, et al. Chai-1: Decoding the molecular interactions of life. bioRxiv. 2024;2024.10.10.615955.
- 37. Cock PJA, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, et al. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics. 2009;25(11):1422–3. pmid:19304878
- 38. Volkamer A, Kuhn D, Rippmann F, Rarey M. DoGSiteScorer: A web server for automatic binding site prediction, analysis and druggability assessment. Bioinformatics. 2012;28(15):2074–5.
- 39. Prota AE, Danel F, Bachmann F, Bargsten K, Buey RM, Pohlmann J, et al. The novel microtubule-destabilizing drug BAL27862 binds to the colchicine site of tubulin with distinct effects on microtubule organization. J Mol Biol. 2014;426(8):1848–60. pmid:24530796
- 40. Kumbhar BV, Borogaon A, Panda D, Kunwar A. Exploring the origin of differential binding affinities of human tubulin isotypes αβII, αβIII and αβIV for DAMA-colchicine using homology modelling, molecular docking and molecular dynamics simulations. PLoS One. 2016;11(5):e0156048. pmid:27227832
- 41. Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ronneberger O, et al. Highly accurate protein structure prediction with AlphaFold. Nature. 2021;596(7873):583–9. pmid:34265844
- 42. Agarwal V, McShan AC. The power and pitfalls of AlphaFold2 for structure prediction beyond rigid globular proteins. Nat Chem Biol. 2024;20(8):950–9. pmid:38907110