Computational and Experimental Verification of Cabozantinib Targeting DDX11 to Inhibit DNA Damage Repair in Liver Cancer.

1
Introduction
Hepatocellular carcinoma
(HCC), the most prevalent form of primary
liver cancer, represents approximately 80% of all liver cancer cases
and remains a significant global health burden. Originating from hepatocytes, the functional cells of the
liver, HCC often arises in the context of chronic liver diseases such
as hepatitis B virus (HBV) and hepatitis C virus (HCV) infections,
cirrhosis, and nonalcoholic fatty liver disease (NAFLD). The global
incidence of HCC has increased markedly in recent decades, particularly
in regions with high prevalence rates of HBV and HCV, such as East
Asia and Sub-Saharan Africa. In 2012 alone, there were approximately
782,000 new cases and 745,000 deaths attributed to liver cancer worldwide,
positioning HCC as the second leading cause of cancer-related mortality
globally. Even in regions with relatively
lower HCC incidence, such as the United States, its occurrence has
been steadily increasing, largely due to rising obesity, metabolic
syndrome, and chronic viral hepatitis rates. These trends underscore the critical need for deeper insights into
the molecular mechanisms underpinning HCC development and progression
to develop more effective diagnostic and therapeutic strategies.
Genetic Mutations often lead to the progression of diverse cancers
including breast cancer, lung cancer and liver cancer. Likewise, HCC pathogenesis involves a complex interplay
of genetic, epigenetic, and environmental factors. Key risk factors
include chronic HBV and HCV infections, excessive alcohol consumption,
obesity, metabolic syndrome, and exposure to carcinogens like aflatoxins.
These factors create a pro-inflammatory and fibrotic hepatic microenvironment
that fosters the transformation of hepatocytes into malignant cells.
Chronic liver injury leads to repeated cycles of cell death, inflammation,
and regeneration, during which genetic and epigenetic alterations
accumulate, driving the initiation and progression of HCC. Among the molecular alterations frequently implicated
in HCC are dysregulated signaling pathways that control cellular proliferation,
apoptosis, angiogenesis, and DNA repair. Specifically, pathways involving
Wnt/β-catenin, PI3K/Akt/mTOR, Ras/MAPK, and DDR mechanisms play
central roles in HCC pathobiology. Moreover,
the DDR system is crucial for maintaining genomic integrity by detecting
and repairing DNA damage, thereby preventing the accumulation of mutations
that can lead to oncogenesis. Among the key mediators of DDR is the
ATM kinase, which orchestrates the repair of double-strand breaks
by phosphorylating serine residues in SQ/TQ motifs of DDR proteins.
This post-translational modification activates downstream repair pathways,
including HR repair, a high-fidelity repair mechanism critical for
resolving DNA double-strand breaks.

One of the proteins regulated by ATM is the DDX11 helicase, a multifunctional
DNA-dependent ATPase involved in DNA replication, sister chromatid
cohesion, and HR-mediated DNA repair. DDX11
plays an essential role in alleviating replication stress and ensuring
proper chromosome segregation during cell division. Dysregulation
or mutations in DDX11 have been implicated in various diseases, including
cancer, due to their impact on DDR pathways. Recent studies have identified a specific mutation in DDX11, Q238H,
which profoundly affects its function in DDR. The mutation alters
the SQ motif comprising serine 237 (S237) and glutamine 238 (Q238),
impairing ATM kinase recognition and preventing phosphorylation at
S237. This phosphorylation event is essential for DDX11’s interaction
with BRCA2, a protein that facilitates the recruitment of RAD51 to
DNA damage sites during HR repair. The inability to phosphorylate
S237 disrupts this interaction, leading to the failure of RAD51 recruitment
and the accumulation of unrepaired DNA lesions. This impairment in DDR compromises genomic stability and
can exacerbate cancer progression. Beyond its impact on DDR, the Q238H
mutation sensitizes HCC cells to poly­(ADP-ribose) polymerase (PARP)
inhibitors.
,
PARP inhibitors exploit synthetic
lethality by targeting tumor cells with defective DDR pathways, further
impairing their ability to repair DNA damage and inducing apoptosis.
The Q238H mutation, by compromising HR repair, enhances the vulnerability
of liver cancer cells to PARP inhibition. This presents a promising
therapeutic opportunity for HCC patients with tumors harboring the
Q238H mutation, allowing for personalized treatment approaches that
exploit these molecular vulnerabilities. Despite its potential, the
rarity of the Q238H mutation highlights the need for more detailed
investigations to fully understand its functional and clinical significance.
This study aims to comprehensively investigate the role of the
Q238H mutation in DDX11 within the context of HCC pathogenesis and
therapeutic sensitivity. The putative inhibitors, specifically targeting
potential DDX11 residues including S237 and Q238, were designed and
validated through comprehensive virtual screening and Molecular dynamics
(MD) simulation-based strategies. Through an integrative approach
combining advanced computational modeling, and drug repurposing assays,
this research aims to deepen our understanding of DDR dysfunction
in HCC and its exploitation for precision medicine. This approach
could revolutionize the drug repurposing process by significantly
reducing the time and cost associated with traditional drug development
methods. Furthermore, this approach enables for further preclinical
and clinical testing of our proposed drugs and also highlights the
comprehensive use of computer-aided strategies in the field of drug
repurposing. By addressing the unique aspects of the S237 and Q238
mutations, this work contributes to the broader goal of improving
treatment outcomes for patients with liver cancer and advancing the
field of personalized oncology.

2
Experimental Section
2.1
Structure Prediction and Refinement
The three-dimensional structure of human DDX11 was obtained from
the AlphaFold Protein Structure Database (https://alphafold.ebi.ac.uk/entry/AF-Q96FC9-F1), which provides open-access structural predictions generated by
the AlphaFold 2 algorithm. The retrieved
model was first inspected for stereochemical quality and refined using
the GalaxyRefine server to optimize side-chain conformations and relieve
steric clashes. The refined structure was subsequently protonated
and energy-minimized in the Molecular Operating Environment (MOE) software to ensure optimal hydrogen-bonding
geometry and stable folding. The final minimized structure was exported
in PDB format for downstream docking and molecular dynamics analyses. The three-dimensional structures of FDA-approved
small-molecule drugs were retrieved from publicly accessible repositories,
including DrugBank (https://go.drugbank.com), ChEMBL (https://www.ebi.ac.uk/chembl/), ZINC20 (https://zinc20.docking.org/), and PubChem (https://pubchem.ncbi.nlm.nih.gov/). A curated library of 2367 bioactive
compounds was assembled for screening. All ligands were geometry-optimized
using the MOE energy minimization module, ensuring proper valence
states, protonation, and polar hydrogen assignment prior to docking
simulations.

2.2
Site-Specific Docking of FDA-Approved Drugs
To identify small-molecule inhibitors of DDX11, AutoDock 4.2 was
employed for molecular docking simulations. The Lamarckian Genetic Algorithm (LGA) was used for conformational
searching, with a population size of 150, a maximum of 2.5 ×
106 energy evaluations, and 100 independent runs per ligand.
The AutoDock scoring function, which combines van der Waals, electrostatics,
hydrogen bonding, desolvation, and torsional free energy terms, was
used to rank binding poses.
,
Site-specific docking
focused on the Ser237 (S237) and Gln238 (Q238) residues, which are
critical for ATM kinase binding. The docking grid was defined using
AutoGrid 4.2, centered on the S237–Q238 region. Grid dimensions
were set to 60 × 60 × 60 points with a grid spacing of 0.375
Å, corresponding to a box size of ∼22.5 Å per axis,
large enough to encompass the active site pocket. The grid was centered
at coordinates (x = 45.231, y =
62.876, z = 38.554 Å). To further validate and
visualize the docking complexes, the results were analyzed using PyMOL
(http://www.pymol.org/pymol), and Discovery Studio. These tools facilitated detailed examination
of the inhibitor–receptor complexes, providing insights into
the molecular interactions underlying the high binding affinities
of the selected drug candidates.

2.3
ADMET Analysis and Drug-Likeness Evaluation
From the initial screening of potential DDX11 inhibitors, the most
promising compounds were selected based on the Lipinski Rule of Five
(Ro5). This analysis ensured that the compounds met the following
criteria: molecular weight below 500 Da, no more than 5 hydrogen bond
donors, no more than 10 hydrogen bond acceptors, and a logP (water
partition coefficient) value below 5. To carry out this analysis, the physicochemical properties of the
shortlisted inhibitors were evaluated using AdmetSAR2 (http://lmmd.ecust.edu.cn/admetsar2/) and SwissADME (http://www.swissadme.ch/). This comprehensive ADMET (absorption, distribution,
metabolism, excretion, and toxicity) analysis was conducted to predict
the safety profile, side effects, and drug-like properties of the
selected compounds in humans. Additionally, the drug-likeness scores
of the compounds were calculated using MolSoft software. This multistep evaluation ensured that the inhibitors
demonstrated favorable characteristics for potential therapeutic development.

2.4
Bioactivity Analysis
The bioactivity
of the final shortlisted compounds was evaluated using the SMILES
(Simplified Molecular Input Line Entry System) canonical data through
the PASS (Prediction of Activity Spectra for Substances) online tool
(https://www.way2drug.com/PASSOnline/predict.php). This platform calculates the Probability
of Activity (Pa) for a given compound, representing the likelihood
of interaction with the target protein. To be considered bioactive,
the Pa value of a compound must exceed its Probability of Inhibition
(Pi) value. A Pa value greater than 0.30 is indicative of medium confidence
in computational predictions and is deemed reliable for bioactivity
assessment. This analysis provided critical insights into the interaction
potential of the selected compounds, supporting their viability as
effective DDX11 inhibitors.

2.5
Molecular Dynamic Simulation
MD simulations
were performed to evaluate the stability and conformational behavior
of the DDX11–ligand complexes. All simulations were carried
out using AMBER22 with the ff19SB force field for protein atoms. Ligand topologies were generated using Antechamber
with the GAFF2 force field and AM1-BCC charges, and systems were prepared
using tleap under physiological protonation conditions (pH 7.4). Each
system was solvated in an orthorhombic TIP3P water box, extending
10 Å from the protein surface, and neutralized with Na+/Cl– counterions. Energy
minimization was conducted in two stages: (i) 2500 steps with restraints
on the protein (1000 steepest descent +1500 conjugate gradient) to
relax solvent molecules, followed by (ii) 2500 unrestrained minimization
steps. Systems were heated from 0 to
300 K over 50 ps under NVT conditions using a Langevin thermostat
(collision frequency 2 ps–1). Equilibration was
performed for 1 ns under NPT, using the Berendsen barostat (τp = 2 ps). Production runs were conducted for 100 ns using
Monte Carlo barostat (NPT) with the Langevin thermostat. This dual-step minimization ensured a well-relaxed
system free from structural anomalies. Following minimization, the
systems were gradually heated to 300 K over 50 ps. A Berendsen barostat
was used only during equilibration to stabilize pressure, while a
Langevin thermostat-maintained temperature control.
,
For production simulations, the Parrinello–Rahman barostat
was employed to ensure accurate pressure coupling and reliable free
energy calculations. Covalent bonds involving hydrogen atoms were
constrained using the SHAKE algorithm. Equilibration was performed for 1000 ps under an NPT ensemble to
achieve system stability before proceeding to the production phase.
,
For each protein–ligand complex, 100 ns MD simulations were
conducted using the GPU-accelerated PMEMD.cuda module in AMBER22. The simulation trajectories were analyzed using
PTRAJ and CPPTRAJ to extract structural stability, atomic fluctuations,
and interaction dynamics.
,
These analyses provided
critical insights into the long-term behavior of the complexes, confirming
their stability and reliability under physiological conditions.

2.6
Free Energy Calculations
The binding
free energies of all simulated DDX11–drug complexes were calculated
using two well-established methods: MM-GBSA (Molecular Mechanics/Generalized
Born Surface Area) and MM-PBSA (Molecular Mechanics/Poisson–Boltzmann
Surface Area). These methods are widely
recognized for determining protein–ligand affinities. Calculations
were performed using the MMPBSA.py module on 100 ns simulation trajectories,
consisting of 10,000 frames sampled at 50-frame intervals. To enhance
structure–function correlation, frames representing different
conformational states were extracted at 500 ps intervals. The ΔG_bind for the complexes
was computed using the equationWhere ΔE_MM represents the total gas phase energy, calculated as the sum of
electrostatic (ΔE_electrostatic),
internal (ΔE_internal), and van
der Waals (ΔE_vdw) energies. The
term ΔG_SOL (PB/GB) is the
sum of polar and nonpolar contributions to solvation-free energy,
as estimated by the GB or PB methods, respectively. The conformational
entropy (TΔS) accounts for
entropy changes during the binding process. The gas phase energy (ΔE_MM) was derived from the molecular mechanics
(MM) force field, with ΔE_internal reflecting bond, angle, and dihedral interactions. Electrostatic
energy contributions (ΔE_electrostatic) were calculated using the MM force field, while van der Waals interactions
(ΔE_vdw) were also included. The
solvation-free energy (ΔG_solvation (PB/GB)) was further divided into polar (ΔG_(PB/GB)) and nonpolar (ΔG_SA‑(PB/GB)) components. The nonpolar solvation term was empirically estimated
using a suitable model, while the polar solvation term was calculated
using either the Poisson–Boltzmann (PB) or Generalized Born
(GB) approach. This comprehensive analysis provided a detailed understanding
of the binding affinities and energy contributions of individual residues,
shedding light on the energetic landscape governing the interactions
between DDX11 and the selected drug candidates.

2.7
Principal Component Analysis
Principal
component analysis (PCA), an unsupervised machine learning approach,
was utilized to analyze the MD trajectories and elucidate the functional
mobility of macromolecules. This analysis
focused on internal protein motions, offering insights into the dynamic
behavior of the complexes. The Cartesian coordinates of carbon (C)
atoms from 5000 snapshots of each simulation trajectory were processed
using the cpptraj package. Positional covariance matrices were generated
to examine atomic motion, followed by the computation of eigenvector
covariance matrices. Diagonal transformation of the positional covariance
matrix produced eigenvalues and eigenvectors, forming the foundation
for PCA representation. The first two principal components, PC1 and
PC2, which represent the most significant collective motions within
the protein structure, were calculated and plotted to track the dynamic
behavior over the course of the simulation. This approach enabled
the identification of major conformational transitions, capturing
essential features of the protein’s functional dynamics. The
relationship between the principal components and the system’s
free energy was expressed mathematically using the equationWhere G is the total free
energy, K
B is the Boltzmann constant, T is the system temperature, and P­(PC1/PC2) represents the
probability function of the principal components. The reaction coordinates,
denoted as (PC1/PC2), provided a quantitative link between the conformational
states and the system’s energy landscape. This methodology
allowed for the quantification and visualization of significant collective
motions within the protein, offering a deeper understanding of conformational
changes and functional dynamics. By integrating PCA with MD trajectories,
the study provided valuable insights into the macromolecular behavior
of the simulated systems.

2.8
Differential Proteomics Based on Data-Independent
Acquisition (DIA) Quantitative Technology
To elucidate the
molecular mechanisms underlying cabozantinib therapeutic effects in
PLC/PRF/5 cells, we conducted differential proteomics analysis using
three biologically independent replicates. Protein extraction was
performed by PBS washing followed by lysis buffer treatment (8 M urea,
100 mM Tris HCl [pH 8.0], protease inhibitors). Samples underwent
three 10-s ice-cold ultrasonication cycles, followed by centrifugation
at 12,000 rpm for 15 min at 4 °C. The supernatant was collected
for protein quantification using the BCA method. For protein digestion,
samples were reduced with 5 mM DTT (56 °C, 1 h), alkylated with
20 mM IAA (room temperature, 45 min in dark), and subjected to overnight
tryptic digestion via FASP method. The digested peptides were collected
by centrifugation (13,000g, 10 min), vacuum-dried,
and desalted using C18 columns. DIA mass spectrometry was performed
in diaPASEF mode with a mass range of 349–1229 m/z and 40 Da isolation windows, employing collision
energy ramping based on ion mobility during PASEF MS/MS scans. Data were processed using Spectronaut 18.0 with
default parameters, including automated retention time alignment and
mass calibration using built-in iRT standards, with intelligent extraction
window optimization. Protein identification criteria: Precursor Threshold
1.0% FDR, Protein Threshold 1.0% FDR. The SwissProt human database served as the reference proteome for
this analysis.

2.9
Protein Purification, Site-Directed Mutagenesis,
and Protein–Ligand Interactions
For structural and
functional studies, a truncated DDX11 construct (residues 227–774)
was cloned into pET-28a­(+) to exclude the N-terminal intrinsically
disordered region (IDR) present in the wild-type protein. The pET-28a­(+)-DDX11 construct was commercially
synthesized (GentleGen, China). The DDX11 S237A-Mut mutant plasmid
was constructed using the multisite mutagenesis kit Mut Express-MultiS
Fast Mutagenesis Kit from Vazyme Biotech. Site-directed mutagenesis
PCR primers for DDX11 S237A were designed, and PCR amplification was
performed using the recombinant plasmid pET-28a­(+)-DDX11 as a template,
followed by plasmid circularization. Successful introduction of the
corresponding amino acid mutation was verified by DNA sequencing.
Both wild-type (WT) and DDX11 S237A mutant proteins were expressed
in Escherichia coli BL21­(DE3). Subsequently,
bacterial cells were harvested by centrifugation at 5000g for 10 min at 4 °C and washed twice with 1× PBS. The bacterial
suspension was then subjected to three freeze–thaw cycles in
liquid nitrogen and ultrasonically lysed on ice (200–300 W
power, 3–5 bursts of 10–20 s each with 30 s intervals).
The lysate was centrifuged at 12,000g for 30 min
at 4 °C, and inclusion body proteins were solubilized using 8
M urea or 6 M guanidine hydrochloride, followed by renaturation with
1× PBS. The target protein was concentrated using a 3 kDa ultrafiltration
tube, and its concentration was finally determined by the BCA method.

The Open-SPR instrument (nicoya, Canada)
enables the detection of protein–ligand interactions. Using 1× PBS buffer as the mobile phase,
a 1:1 mixture of 50 mM N-ethyl-N′-(3-diethylaminopropyl) carbodiimide and 0.2 M N-hydroxysuccinimide was prepared to activate the COOH-modified gold
nanoparticle sensor chip. The mixture was diluted to a final concentration
of 50 μg/mL with 10 mM sodium acetate buffer (pH 5.0) and injected
at a flow rate of 35 μL/min for 500 s. Subsequently, diluted
DDX11 protein was coupled to the activated chip using the mobile phase,
followed by a 5 min injection of 1 M ethanolamine (pH 8.5) as blocking
buffer to deactivate remaining reactive groups. Drug molecules were
diluted to different concentration gradients, and the binding interactions
between each sample and chip-immobilized DDX11 generated distinct
response curves. The dissociation constant (K
D) was determined using an A + B ↔ AB binding model
for data fitting.

3
Results and Discussion
3.1
Structural Modeling and Confidence Assessment
The AlphaFold-derived model represents the full-length human ATP-dependent
DNA helicase DDX11 (970 amino acids) and comprises two principal regions:
an N-terminal helicase core and a C-terminal zinc-binding domain.
Model reliability was assessed using predicted Local Distance Difference
Test (pLDDT) scores and Predicted Aligned Error (PAE) metrics, which
indicated high confidence across the helicase core and moderate flexibility
in peripheral regions (Supporting Figure S1a–c). The DDX11 model exhibited very high confidence (pLDDT > 90)
across
the helicase core, ATP-binding cleft, and DNA-binding motifs, visualized
in dark blue on the confidence map. Moderate confidence regions (70
< pLDDT < 90) were mainly localized to flexible loops and linker
segments connecting subdomains, consistent with their inherent conformational
variability. In contrast, the N- and C-terminal tails displayed lower
confidence (pLDDT < 70), suggesting flexibility and possible disorder,
which is typical of peripheral helicase extensions. The PAE heatmap
indicated an expected positional error below 5 Å within each
structured domain, confirming accurate intradomain residue alignment.
Slightly elevated PAE values (10–15 Å) between domains
reflected moderate interdomain uncertainty, a common feature in large,
multidomain helicases that undergo conformational transitions during
ATP hydrolysis. Collectively, these findings confirm that the helicase
catalytic domain and nucleotide-binding pocket are modeled with high
accuracy, validating their suitability for subsequent docking and
molecular dynamics (MD) simulations.

3.2
Structural Refinement and Preparation for
Docking
To further enhance stereochemical quality, the AlphaFold
2 model was refined using the GalaxyRefine Web server, which iteratively
optimizes side-chain conformations through relaxation and energy minimization.
Among the five refined models generated, models d and f exhibited
the most favorable stereochemical statistics and lowest overall energy
values (Supporting Figure S1d,e). The best
refined model which demonstrated the best MolProbity and Ramachandran
metrics, was subsequently selected for all downstream analyses. This
optimized structure was protonated and energy-minimized in the Molecular
Operating Environment (MOE) to achieve optimal hydrogen-bond geometry
and backbone conformation. The final refined model exhibited reduced
steric clashes, improved backbone dihedral distributions, and favorable
stereochemical parameters, confirming its reliability for receptor-based
simulations. The model also retained high pLDDT confidence across
the Walker A, Walker B, and helicase C motifs, reinforcing its functional
integrity. Notably, residues S237 and Q238, situated in the high-confidence
catalytic region of the helicase core, were clearly resolved and structurally
stable, providing a robust foundation for assessing cabozantinib binding
interactions and the structural consequences of the Q238H mutation
during MD simulations.

3.3
Virtual Screening
The predicted three-dimensional
structure of DDX11, obtained from the AlphaFold 2 model and refined
through energy minimization and stereochemical optimization, was employed
for virtual screening to identify potential inhibitors. A curated
library of 2367 FDA-approved pharmacological compounds was docked
against DDX11 using AutoDock 4.2. To ensure comprehensive identification
of potential binding sites, an initial global docking approach was
used to explore the entire protein surface without predefined binding
constraints. Subsequently, focused redocking was performed around
the S237 and Q238 residues, which are catalytically relevant and implicated
in DDX11 functional regulation, to verify specific hydrogen-bonding
interactions. Binding affinities (kcal/mol) were analyzed to evaluate
interaction strengths, while binding conformations revealed the three-dimensional
orientations of ligands within the target pocket. From this two-stage
screening, seven compounds exhibiting binding affinities between −6.5
and −7.7 kcal/mol were shortlisted as promising candidates
(Table
). For this
study, a docking score threshold of approximately −6.5 kcal/mol
was adopted, as compounds below this cutoff typically exhibit stable
and favorable protein–ligand interactions. The selected compounds
were subjected to subsequent interaction and stability analyses to
identify key inhibitors of DDX11.
Among the shortlisted compounds, Afatinib displayed
the strongest
predicted binding affinity (−7.674 kcal/mol). However, detailed
binding mode analysis revealed that Afatinib did not establish stable
hydrogen-bond interactions with the key catalytic residues Ser237
and Gln238, which are essential for DDX11 inhibition. Consequently,
despite its higher docking score, Afatinib was deprioritized due to
the lack of critical residue engagement. In contrast, dacomitinib
(−7.108 kcal/mol), ergotamine (−6.814 kcal/mol), and
cabozantinib (−6.567 kcal/mol) exhibited both favorable binding
affinities and strong hydrogen-bond interactions with Ser237 and Gln238.
Among the screened compounds, cabozantinib was prioritized based on
its unique interaction with Ser237, a known phosphorylation site of
ATM kinase, and its established clinical use in hepatocellular carcinoma.
This combination of favorable binding features and clinical relevance
highlighted cabozantinib as a strong candidate for downstream validation.
Collectively, these results supported the selection of dacomitinib,
ergotamine, and cabozantinib as the top candidates for subsequent
molecular dynamics simulations and experimental validation. The binding
interactions of these compounds with DDX11 at the Ser237–Gln238
site are shown in Figure
A–C. We presented the potential cocrystallization regions
for the binding of these three ligands to the protein (Figure
D).

3.4
ADMET Analysis and Drug-Likeness Evaluation
The drug-likeness (DL) properties of the shortlisted compounds
were evaluated using ADMET-based screening tools, with Lipinski’s
Rule of Five (Ro5) serving as a benchmark. According to these criteria,
drug candidates are generally expected to have fewer than five hydrogen
bond donors, fewer than ten hydrogen bond acceptors, a molecular weight
below 500 Da, and a Log P value not exceeding
five. In addition, a drug-likeness score threshold of greater than
0.18, determined using the Molsoft toolset, was applied to further
refine the selection. Although all compounds under investigation are
FDA-approved, their repurposing for liver cancer therapy necessitates
assessment of physicochemical and pharmacokinetic properties, particularly
given the potential for administration at higher doses. As shown in Table
, dacomitinib fully
complies with Lipinski’s Ro5, whereas ergotamine (581.66 Da)
and cabozantinib (501.51 Da) slightly exceed the molecular weight
threshold. These deviations highlight that not all selected compounds
strictly adhere to the classical Ro5. However, it is important to
note that such exceptions are common among anticancer drugs, where
higher molecular weight and complexity can still be compatible with
oral bioavailability and therapeutic activity. Beyond Lipinski’s
criteria, ADMET profiling indicated acceptable absorption, distribution,
metabolism, and toxicity profiles for all three compounds. This balance
between efficacy and safety is a critical consideration in cancer
drug repurposing, where therapeutic benefit often outweighs minor
rule violations. Furthermore, Prediction of Activity Spectra for Substances
(PASS) analysis confirmed their classification as biologically active
compounds. Collectively, these findings support the continued consideration
of dacomitinib, ergotamine, and cabozantinib as viable drug candidates
despite minor deviations from strict Ro5 compliance.

3.5
Molecular Dynamic Simulation
MD simulations
were performed using the Amber software suite to investigate the time-dependent
structural changes and stability of DDX11–drug complexes over
a 100 ns simulation period. The dynamic behavior of the complexes
and any degenerative changes in the DDX11 enzyme upon drug binding
were assessed by analyzing key parameters, including root-mean-square
deviation (RMSD), root-mean-square fluctuation (RMSF), protein–ligand
interactions, β-factor, and the radius of gyration. These metrics
provided a comprehensive evaluation of the structural stability and
flexibility of the complexes. Simulations were conducted on three
selected drug candidates: dacomitinib, ergotamine, and cabozantinib.
Although cabozantinib uniquely interacts with Ser237, dacomitinib
and ergotamine were also included in molecular dynamics simulations
to assess the robustness of the computational workflow and to enable
comparative evaluation of alternative binding modes involving Gln238
(Q238). Such Q238-centered interactions may be relevant in mutation
contexts (e.g., Q238H) or inform future structure–activity
relationship (SAR) studies for DDX11 inhibitors. Notably, all three
complexes exhibited high structural stability throughout the 100 ns
simulation, underscoring their potential to act as effective inhibitors
of DDX11 activity. By inhibiting DDX11, these compounds may interfere
with pathways critical for tumor progression, offering promising therapeutic
avenues for preventing liver cancer.

3.6
Root Mean Square Deviation Analysis
RMSD is a key metric in postsimulation trajectory analysis used to
assess the stability and binding conformation of protein–ligand
complexes during MD simulations. The
dynamic stability and conformational behavior of DDX11 upon ligand
binding were examined through 100 ns molecular dynamics (MD) simulations
using RMSD, RMSF, and radius of gyration (RoG) as key structural metrics
(Figure
). The backbone
RMSD trajectories for the dacomitinib, ergotamine, and cabozantinib
DDX11 complexes (Figure
A–C) showed an initial relaxation phase within the first 20–30
ns, during which the complexes adapted from their docked conformations
to more energetically favorable states. After this equilibration period,
the RMSD values plateaued between approximately 3.5 and 4.0 Å,
indicating that no major conformational drift or unfolding occurred.
Such moderate RMSD values are characteristic of local rearrangements
within the binding pocket rather than global structural instability,
reflecting the complexes’ attainment of equilibrium within
the protein’s active site.
The residue-wise RMSF profiles (Figure
D–F) further revealed
that flexibility
was primarily confined to loop regions and the N- and C-terminal ends
of the protein, whereas the helicase core domain remained comparatively
rigid throughout the trajectory. Importantly, the residues surrounding
Ser237 and Gln238 which constitute the catalytic and ligand-binding
region of DDX11 displayed minimal fluctuations (<2 Å), confirming
the structural integrity of the binding pocket and the stability of
protein–ligand contacts during the simulation. These restrained
fluctuations suggest that ligand binding effectively stabilized the
local environment near the active site.
The RoG plots (Figure
G–I) demonstrated
that the global compactness of the
DDX11–ligand complexes was maintained over the 100 ns simulation,
with values fluctuating narrowly between 54.5 and 55.6 Å. The
absence of significant deviations in RoG indicates that none of the
ligands induced expansion or collapse of the protein structure, corroborating
the RMSD results. Collectively, the RMSD, RMSF, and RoG analyses confirm
that all three complexes remained structurally stable and well-equilibrated
during the simulation period. The data suggest that dacomitinib, ergotamine,
and cabozantinib preserve the conformational integrity of DDX11 while
allowing subtle adaptive flexibility conducive to stable and persistent
binding interactions within the helicase core.

3.7
Hydrogen Bond Analysis
Hydrogen bonds
and hydrophobic interactions at the protein–ligand interface
play a pivotal role in strengthening complex stability and enhancing
binding interactions. To evaluate these
interactions at the atomic level, we quantified all hydrogen bonds
formed in each system based on established criteria. A hydrogen bond
was considered formed when the donor–acceptor distance was
less than 0.35 nm, and the donor-hydrogen-acceptor angle was within
30°. The temporal analysis of hydrogen
bond formation, presented in Figure
, illustrates robust hydrogen bonding networks across
all three complexes throughout the simulation. This finding highlights
the substantial binding affinities of the proposed drug candidates,
enabling precise and stable binding to DDX11. Such interactions are
critical for disrupting the DDX11-ATM kinase interaction, offering
potential therapeutic benefits against liver cancer progression. Moreover,
the consistent presence of hydrogen bonds during the simulation underscores
the strong and stable nature of these complexes. The analysis confirmed
that each system maintained a considerable number of hydrogen bonds,
which is integral to preserving the structural integrity of the protein–ligand
complexes. These robust hydrogen bonding networks are not only essential
for stabilizing the complexes but also enhance their therapeutic efficacy,
further supporting their potential as effective inhibitors of DDX11-ATM
kinase interactions.

3.8
Free Energy Calculations
Binding
free energy is composed of several key components, including van der
Waals energy (vdW), electrostatic energy (EEL), surface energy (ESURF),
Poisson–Boltzmann energy (EPB), polar energy (EnPolar), Generalized
Born energy (EGB), and surface energy (ESURF). The MM-GBSA methodology
was employed to re-evaluate ligand-docked conformations, enabling
identification of the most accurate ligand orientations and corresponding
binding free energies. For the dacomitinib, ergotamine, and cabozantinib
complexes with DDX11, the calculated MM-GBSA binding free energy values
were −36.18 ± 7.32, −26.35 ± 3.36, and −34.26
± 6.79 kcal/mol, respectively (Table
). These findings highlight that the proposed
drug–receptor complexes exhibited favorable binding free energy
values, indicative of high binding affinity and structural stability.
In addition to cabozantinib, analysis of dacomitinib and ergotamine
enabled comparison of energetics associated with Q238-mediated binding
modes. Importantly, these three ligands were not selected solely on
the basis of docking energy but rather for their biologically relevant
interactions with key residues and therapeutic potential. The results
of MM-GBSA and MM-PBSA analyses consistently identified all three
ligands as strong binders to DDX11 (Table
and Figure
), reinforcing the stable interactions observed in
MD simulations and supporting their potential to disrupt the DDX11–ATM
kinase interaction, a pathway critical to liver cancer progression.

3.9
Per-Residue Energy Decomposition Analysis
To further dissect the molecular determinants of binding, per-residue
energy decomposition was performed using MM-GBSA on the 100 ns MD
trajectories Table
. This approach allowed quantification of the energetic contributions
of individual amino acids within the binding pocket of DDX11. The
analysis revealed that Ser237 and Gln238 contribute substantially
to ligand stabilization, particularly in the cabozantinib complex,
consistent with docking predictions and our Surface Plasmon Resonance
(SPR) validation.

3.10
Principal Component Analysis
PCA
was performed on 100 ns MD trajectories to investigate and compare
the conformational changes induced by protein–ligand interactions.
The analysis revealed distinct differences among the complexes, alongside
shared cluster-like motion patterns (Figure
). All three complexes exhibited relatively
compact conformational behavior, with average PC1 values ranging from
−400 to +300 and PC2 spanning −200 to +300. These findings
indicate that ligand binding induces restricted protein dynamics in
the complexes, suggesting stabilization of specific conformational
states. Distinct motion patterns associated with Q238-centered interactions
were observed for dacomitinib and ergotamine relative to the Ser237-targeting
cabozantinib complex. These Q238-mediated dynamics may become especially
relevant in the context of mutations such as Q238H, which could alter
binding preferences. The observed variations in principal component
distributions therefore not only highlight the stabilizing effect
of cabozantinib but also broaden our understanding of alternative
conformational landscapes that may guide future structure–activity
relationship (SAR) studies.

3.11
Differential Proteomic Validation of Cabozantinib
Treatment in PLC/PRF/5 Liver Cancer Cells
In our previous
study, we confirmed that the phosphorylation of DDX11 Ser237 by ATM
kinase is a key factor in enhancing HR repair capability in HCC cells. Through molecular docking, we predicted that
three drugs dacomitinib, ergotamine, and cabozantinib could interact
with DDX11. However, dacomitinib and ergotamine are primarily used
to treat metastatic nonsmall cell lung cancer with EGFR-sensitive
mutations and acute migraines, respectively, while cabozantinib can
be applied in HCC therapy. Based on the virtual screening results,
cabozantinib was selected for experimental validation in PLC/PRF/5
liver cancer cells. Its interaction with the Ser237–Gln238,
the key phosphorylation site, and has established therapeutic relevance
in hepatocellular carcinoma. More importantly, only cabozantinib can
target the DDX11 Ser237 site. Therefore, in this study, we focused
more on the role of cabozantinib in treating HCC and its relationship
with DDX11.
Subsequent differential proteomics analysis was
performed on three independent biological replicate samples using
DIA quantitative technology (Figure
A). Compared with the NC group, 818 differentially
expressed proteins were identified in the cabozantinib treatment group,
including 139 upregulated proteins and 679 downregulated proteins
(Figure
B). GO enrichment
analysis of these 818 differentially expressed proteins showed that
they were significantly enriched in biological processes, cellular
components and molecular functions related to ribosome biogenesis,
rRNA metabolic process, transcription coactivator activity, and histone
binding (Figure
C).
Further KEGG pathway analysis revealed that the differentially expressed
proteins were significantly enriched in Carbon metabolism and Autophagy
processes, which are associated with various diseases including Alzheimer’s
disease, Parkinson’s disease, and cancers (Figure
D). These results collectively
indicate that cabozantinib can disrupt translation and metabolic pathways
in PLC/PRF/5 cells, which may inhibit DNA damage repair.

3.12
Cabozantinib Inhibits the DNA Damage Repair
Ability of PLC/PRF/5 Liver Cancer Cells
Our cluster analysis
of the 818 differentially expressed proteins revealed that 26 DNA
damage repair-related proteins were significantly downregulated (Figure
A) (Table S1). This may be attributed to cabozantinib inhibiting
the active site of DDX11 Ser237, thereby preventing its phosphorylation
and subsequently impairing the DNA damage repair capacity of PLC/PRF/5
cells. We then constructed an interaction network of these 26 DNA
damage repair-related proteins, identifying six proteins that interact
with DDX11: FEN1, CDCA5, FANCA, MCM8, RAD51, and FANCD2 (Figure
B). Within this network,
RAD51a core DNA repair proteinexhibited extensive
interactions with multiple proteins. Western blot analysis confirmed
that, compared to the NC group, the expression levels of FEN1, CDCA5,
FANCA, MCM8, RAD51, and DDX11 were downregulated in the cabozantinib-treated
group (Figure
C).
Since ATM kinase and BRCA2 were not detected in our proteomics data,
we assessed the mRNA levels of FEN1, CDCA5, FANCA, MCM8, RAD51, FANCD2,
DDX11, ATM, and BRCA2 via qRT-PCR (Figure
D) (Table S2).
The results demonstrated that cabozantinib significantly suppressed
the transcriptional levels of these genes (P <
0.05). These findings validate the reliability of our proteomics data
and confirm that cabozantinib targets the DNA damage repair pathway,
particularly DDX11.
Notably, we evaluated the HR repair efficiency
in PLC/PRF/5 cells
and found that both shDDX11#1 and shDDX11#2 exhibited significantly
reduced HR repair capacity (P < 0.01), while the
cabozantinib-treated group also showed a decline in HR repair activity
(P < 0.05) (Figure
E). In contrast, no significant differences were observed
in nonhomologous end joining (NHEJ) repair efficiency among the shDDX11#1,
shDDX11#2, and cabozantinib-treated groups compared to controls (Figure
F). These results
demonstrate that cabozantinib suppresses HR repair in PLC/PRF/5 cells
by inhibiting the expression of DNA damage repair-related proteins,
highlighting its potential therapeutic application in liver cancer.

3.13
Cabozantinib Targets and Binds to the Ser-237
Residue Site of DDX11
To validate the binding between cabozantinib
and DDX11, we expressed and purified the core region of DDX11 (residues
227–774) in vitro. SDS-PAGE results showed that DDX11 had an
approximate molecular weight of 57 kDa with high purity (Figure
A). We then tested
the interaction using cabozantinib at concentrations of 0.0625 μM,
0.125 μM, 0.25 μM, 0.5 μM, 1.0 μM, 2.0 μM,
and 5.0 μM. Fitting the data from three independent experiments
revealed that cabozantinib exhibited strong binding affinity to WT
DDX11, with a K
D value of 0.726 ±
0.03 μM (Figure
B). Next, we generated the DDX11 S237A mutant protein by site-directed
mutagenesis, substituting serine 237 with alanine. Binding assays
with cabozantinib at concentrations of 0.25 μM, 0.5 μM,
1.0 μM, 2.0 μM, 4.0 μM, 8.0 μM, 16.0 μM
demonstrated a significant decrease in affinity for the mutant protein,
with a K
D value of 536.315 ± 0.70
μM (Figure
C).
These results demonstrate that cabozantinib primarily targets the
Ser237 site, which would severely inhibit its activity. In conclusion,
our biochemical characterization is consistent with the reliability
of the molecular docking predictions.

4
Conclusions
This study employed a structure-based
drug repurposing approach
to identify inhibitors targeting the critical S237 and Q238 residues
in DDX11, key phosphorylation sites for ATM kinase involved in liver
cancer progression. Screening 2367 FDA-approved drugs provided seven
candidates, with three showing strong hydrogen bonding to the target
residues. Molecular dynamics simulations confirmed the stability of
these compounds, supported by consistent RMSD, RMSF, RoG, and robust
hydrogen bonding. MM-GBSA/MM-PBSA analyses revealed high binding affinity,
while PCA indicated ligand-induced stabilization of DDX11 conformations.
These findings highlight the potential of the proposed compounds to
disrupt DDX11-ATM kinase interactions and inhibit liver cancer progression.
Given that cabozantinib possesses anticancer activity and this candidate
drug has the potential to directly target DDX11 Ser237, subsequent
differential proteomics data supported that cabozantinib causes downregulation
of DNA damage repair-related proteins in hepatocellular carcinoma
cells, thereby inhibiting HR repair capacity. Further in vitro expression
and biochemical characterization experiments demonstrated that cabozantinib
may specifically target DDX11 Ser237. Future studies are encouraged
to validate these results through clinical investigations, advancing
these compounds toward therapeutic applications in targeted cancer
treatment.

Supplementary Material