Molecular Dynamics Simulation-Assisted siRNA Design for Dual-Ubiquitinated SKP2 Silencing via Ago2 Anchoring in Breast Cancer.

1
Introduction
Amidst the rising global
burden of oncological diseases, breast
cancer has emerged as a major contributor to cancer-related morbidity
and mortality. In the United States, it is projected to account for
nearly 30% of all new cancer diagnoses in 2025, with over 316,000
invasive cases anticipated. Internationally, its incidence continues
to climb at an annual rate of 1–5% across diverse populations,
underscoring an urgent need for deeper molecular insights and the
development of targeted therapeutic strategies. Recently, among an array of oncogenic drivers implicated
in tumorigenesis, S-phase kinase-associated protein 2 (SKP2) has been
identified as frequently overexpressed across diverse human malignancies,
with breast cancer standing out as a particularly prominent case.
For example, elevated levels of both SKP2 mRNA and protein have been
consistently observed in breast cancer cell lines and primary tumor
tissues.
,
It is a type of F-box protein that helps
form the Skp1-Cullin-F-box (SCF) complex, a molecular
machine that tags other proteins for breakdown or signaling changes.
,
In this SCF setup, which includes parts like Cul1, Rbx1, Skp1, and
SKP2, the protein acts as the ″scout″ that picks out
specific targets. SKP2 has been firmly
established as a cancer-promoting gene, and it is often overproduced in 40–60% of breast tumors.
This overproduction is linked to more aggressive features, like higher
tumor grades, spread to lymph nodes, and shorter time before the cancer
worsens, making SKP2 a promising target for new therapies in breast
cancer. Beyond breast cancer, SKP2 overexpression has also been reported
in liver, prostate, colorectal, and pancreatic
cancers, underscoring its broader oncogenic
relevance. These observations suggest that SKP2 functions as a broadly
relevant oncogenic regulator rather than a tissue-restricted factor.
Importantly, the molecular mechanisms governed by SKP2, including
ubiquitin-mediated degradation of cell-cycle inhibitors and modulation
of survival pathways, are conserved across diverse tumor contexts.
Nevertheless, the present work specifically focuses on breast cancer,
in which SKP2 exhibits distinct dual ubiquitination functions that
are closely linked to therapeutic resistance and disease progression,
thereby providing a strong biological rationale for its targeted silencing
in this setting.
SKP2 is an E3 ubiquitin ligase complex within
the ubiquitin–proteasome
system (UPS), which orchestrates the
ATP-dependent ubiquitination of target proteins to regulate cellular
homeostasis, cell-cycle progression, and oncogenesis, exemplified
by its K48-linked polyubiquitination of the cyclin-dependent kinase
inhibitor p27 to promote its proteasomal degradation
,
and its K63-linked ubiquitination of the kinase Akt1 to enhance
its activation and signaling;
,,
both mechanisms dysregulated in breast cancer to drive proliferation
and survival.

A well-known substrate
of SKP2 is the cyclin-dependent kinase inhibitor
p27̂Kip1. Specifically in breast cancer, SKP2 acts as the substrate-recognition
subunit of the SKP2 ubiquitin ligase complex to mediate p27̂Kip1
proteasomal degradation. This mechanism disrupts normal G1/S cell-cycle
progression and fosters malignant transformation.
,
The degradation process begins with CDK2/cyclin E-dependent phosphorylation
of p27 at threonine 187. This modification promotes Cks1-assisted
binding to SKP2, which in turn assembles K48-linked polyubiquitin
chains to direct p27 toward 26S proteasome destruction. Furthermore, SKP2 overexpression defines aggressive
breast cancer subtypes such as triple-negative and HER2-enriched tumors.
,
It inversely correlates with p27 abundance, thereby accelerating
proliferation, invasion, epithelial–mesenchymal transition,
and resistance to drugs including doxorubicin, paclitaxel, and tamoxifen.
As a result, these alterations predict unfavorable clinical outcomes
and reduced survival rates. Consequently,
targeting the SKP2-p27 axis with small-molecule inhibitors or proteolysis-targeting
chimeras (PROTACs) preserves p27 levels, triggers G1 arrest, and boosts
chemosensitivity. These strategies thus highlight SKP2 as a promising
therapeutic target in breast cancer.

In addition to p27 degradation, SKP2 promotes oncogenic signaling
by modulating Akt1, a serine/threonine kinase that is essential for
cell survival, metabolism, and growth pathways. In various cancers, SKP2 goes beyond proteolytic degradation
and instead elevates Akt1 signaling through K63-linked nonproteolytic
ubiquitination, promoting Akt1 localization to the membrane and enhancing
kinase signaling.
,
Amplification of Akt1 kinase
signaling activation can partially explain the survival advantage
and chemoresistance in breast cancer cells during drug exposure, highlighting
the multiple oncogenic functions of SKP2.

Functioning as a pivotal node in oncogenic signaling networks,
SKP2 acts as a dual E3 ubiquitin ligase, mediating K48-linked ubiquitination
of p27̂Kip1 to promote its proteasomal degradation and K63-linked
ubiquitination of Akt1 to facilitate its membrane localization and
activation. This bifunctional activity positions SKP2 as a central
regulator of cell-cycle progression and survival signaling in breast
cancer. Accordingly, therapeutic strategies aimed at attenuating SKP2
expression or function such as by designing RNA interference molecules
or developing structure-based small-molecule inhibitors have garnered
considerable interest as potential interventions to counteract tumor
growth and resistance to apoptosis.

Small interfering RNAs (siRNAs) are 21–23 nucleotides long
with 2-nucleotide 3′ overhangs and drive gene silencing through
the conserved RNA interference pathway in eukaryotic cells. In the
cytoplasm, siRNAs arise either from longer double-stranded RNA processed
by the RNase III enzyme Dicer or as directly introduced synthetic
duplexes.
−

These duplexes integrate into the RNA-induced silencing
complex (RISC) where the guide (antisense) strand persists while the
passenger (sense) strand degrades.
,
The TAR RNA-binding
protein (TRBP) supports Dicer by stabilizing the siRNA and aiding
its transfer to Argonaute 2 (Ago2), the primary catalytic effector
in RISC. Among human Argonaute proteins (Ago1–Ago4), only Ago2
possesses the endonucleolytic activity required for mRNA cleavage.
Its PAZ domain anchors the siRNA’s 3′ overhang, and
the MID domain secures the 5′ phosphate ensuring accurate guide
strand positioning. This precise alignment
allows Ago2 to pair with and cleave target mRNA effectively, halting
translation and reducing protein expression.

Despite their precision in gene silencing, siRNAs encounter
significant
therapeutic delivery challenges. Their susceptibility to nuclease
degradation in serum and their polyanionic, hydrophilic nature hinders
cellular internalization. To address these, siRNAs are chemically
modified or encapsulated in lipid nanoparticles or polymeric carriers,
enhancing serum stability, cellular uptake, endosomal escape, and
cytosolic delivery to maximize RNA interference efficacy while minimizing
off-target effects. Compared to small-molecule
inhibitors or monoclonal antibodies, siRNAs offer superior specificity
with reduced systemic toxicity when carefully designed. Their therapeutic potential, exemplified by
targeting genes like SKP2 in breast cancer, depends on optimizing GC content, thermodynamic stability, and mRNA
secondary structure for efficient silencing. Mitigating off-target
interactions and immune activation is essential to ensuring precision
and safety.

Designing effective
siRNA molecules through in silico approaches
offers a cost-efficient and high-throughput strategy for gene silencing
applications. In this study, we present the first computational framework
for RNA interference-based targeting of the SKP2 oncogene in breast
cancer. Candidate siRNAs were systematically designed and screened
using established selection criteria, including the Ui-Tei, Amarzguioui,
and Reynolds rules, GC content optimization, melting temperature (T
m), and off-target minimization. Thermodynamic
profiling and secondary structure prediction were employed to assess
duplex stability and guide strand accessibility. To ensure functional
compatibility with the RNA-induced silencing complex (RISC), molecular
docking and molecular dynamics (MD) simulations were conducted to
evaluate siRNA interactions with human Argonaute 2 (hAgo2), focusing
on domain-specific anchoring within the MID and PAZ regions. This
integrative computational strategy enabled the identification of siRNA
candidates with high specificity, structural stability, and favorable
binding energetics, laying the groundwork for RNAi-based precision
therapeutics targeting SKP2 in breast cancer.

2
Materials and Methods
The workflow
adopted in this study, along with the algorithms and
tools employed at each computational tier, is systematically summarized
in Table
.
2.1
Foundational Data Acquisition for Structure-Based
siRNA Design Targeting SKP2 in Breast Cancer
SKP2 (S-phase
kinase-associated protein 2) is an F-box protein involved in ubiquitin-intermediated
cell-cycle control and acts as an element of the cyclin A CDK2 kinase
complex. The full coding DNA sequence (CDS) of human SKP2 mRNA was
retrieve from the NCBI Nucleotide database. The retrieve SKP2 mRNA sequence was employed in computational siRNA
design tools to identify seeker small interfering RNAs (siRNAs) with
high particularity for targeting and silencing SKP2 expression. Furthermore,
to perform molecular docking and dynamics, the three-dimensional structure
of the Argonaute 2 (Ago2) protein was stained from AlphaFold Protein
Structure Database. This database is
a slice-edge AI system developed by DeepMind that uses deep learning
neural networks to prognosticate protein 3D structures from amino
acid sequence. The model was trained on thousands of experimentally
determined structures from the Protein Data Bank (PDB), enabling it to learn the complex mapping between
sequence patterns and 3D conformation.

2.2
Designing Candidate siRNAs for SKP2 Silencing
To ensure a rational and biologically meaningful selection of siRNA
candidates, a stepwise exclusion strategy was applied based on widely
accepted computational filtering criteria reported in the literature
for maximizing silencing efficiency and specificity. The selection
process incorporated multiple parameters, including established sequence-based
design rules (Ui-Tei, Amarzguioui, and Reynolds), optimal GC content
thresholds, seed duplex stability, off-target filtering through BLAST
searches, mRNA secondary structure accessibility, and thermodynamic
properties. These parameters are consistently emphasized in recent
computational siRNA design studies targeting disease-associated genes,
thereby providing a reliable framework for identifying potent and
specific SKP2-targeting siRNAs.
,

siRNA candidates
were initially predicted using siDirect version 2.0, with the full-length human SKP2 mRNA sequence as input.
The tool applied three well-established design rules Ui-Tei, Reynolds, and Amarzguioui to ensure high specificity and silencing efficiency.
After filtering based on thermodynamic stability and off-target potential,
a final set of siRNAs was selected. These guide strands met all key
criteria and are considered strong candidates for further experimental
validation due to their high predicted efficacy and minimal risk of
off-target effects. To minimize potential off-target effects, the
maximum melting temperature (T
m) of the
seed duplex was restricted to 21.5 °C. The GC content of the
siRNAs was constrained within the range of 30–55% ensuring
an optimal balance between duplex stability and efficient strand separation
during RISC loading. This range allows the designed siRNAs to maintain
sufficient binding affinity with the target mRNA while preserving
specificity, thereby enhancing cleavage efficiency. While higher GC
content can strengthen mRNA binding, it simultaneously increases the
risk of unintended hybridization with nontarget transcripts, which
was carefully avoided in the selection process.

2.3
Computational Evaluation of Structural Flexibility
and Thermodynamic Stability in siRNA Candidates Targeting SKP2
The guanine–cytosine (GC) composition of each siRNA guide
strand was quantitatively assessed using the Oligonucleotide Properties
Calculator provided by Bio-Synthesis Inc. Nucleotide sequences were input in single-stranded RNA format, and
the tool computed the GC percentage to evaluate duplex stability and
optimize base-pairing fidelity. Additionally, the thermodynamic stability
of each guide strand was assessed through Gibbs free energy calculations
under standard physicochemical conditions. The Gibbs free energy (ΔG) values reflective
of the intrinsic spontaneity of siRNA guide strand folding were computed
under canonical conditions 1 M NaCl, 25 °C, and pH 7 as defined
configuration of the Oligonucleotide Properties Calculator. However,
the ΔG values (reported in kcal/mol) were utilized
exclusively for comparative thermodynamic profiling across siRNA candidates.
Their purpose was to inform relative structural stability during the
selection phase, rather than to replicate physiological environments.
To elucidate the secondary structures of single-stranded siRNA
candidates, the MaxExpect algorithm in the RNAstructure web server
(Version 6.0.1, Mathews Lab, University of Rochester Medical Center)
,
was utilized. This computational tool predicts RNA and DNA secondary
configurations by prioritizing base-pairing interactions with the
highest probabilistic accuracy. In this analysis, each 23-nucleotide
siRNA guide strand was submitted under defined parameter settings.
MaxExpect generates a structure that maximizes the cumulative probability
of correctly predicted base pairs, thereby yielding the most statistically
reliable conformation. Notably, this probabilistic modeling strategy
has demonstrated superior concordance with experimentally validated
RNA structures compared to conventional minimum free energy-based
approaches.
The DI-Nucleic Acid hybridization and melting prediction
DINAmelt
modules, accessible via the UNAFold Web
Server, was employed to characterize
the thermal dissociation behavior of siRNA–mRNA duplexes. This
predictive framework utilizes nearest-neighbor thermodynamic algorithms,
which account for the contextual stability of each base pair relative
to its flanking nucleotides. By simulating the temperature-dependent
denaturation process, the tool generates equilibrium melting profiles
and computes key thermodynamic metrics including melting temperatures T
m (Conc) and T
m (Cp),
which serve as critical indicators of duplex robustness and silencing
potential. For each siRNA candidate, the linear nucleotide sequence
alongside its complementary mRNA target was submitted as input. Melting
temperature values were derived under defined simulation parameters
using an energy minimization approach, thereby facilitating comparative
assessment of duplex stability and guiding the selection of the most
thermodynamically favorable siRNA constructs.
The secondary
structures of individual siRNA guide strands were
predicted using the Mfold algorithm. This
tool employs dynamic programming techniques and empirically validated
thermodynamic parameters to identify the most energetically favorable
conformations for a given nucleotide sequence. Each siRNA strand was
submitted in a linear, single-stranded RNA format, appropriate for
structural modeling. Simulations were conducted under default conditions
37 °C folding temperature and 1 M NaCl ionic strength to calculate
the Gibbs free energy (ΔG) associated with
secondary structure formation. These ΔG values
served as indicators of structural stability, aiding in the comparative
evaluation of siRNA candidates.
To evaluate the hybridization
dynamics and thermodynamic stability
of siRNA–mRNA interactions, the DuplexFold algorithm from the
RNAstructure web server (Version 6.0.1) was employed. This computational tool predicts the most energetically
favorable intermolecular RNA duplexes by exhaustively analyzing all
possible base-pairing configurations between two distinct nucleotide
strands, while explicitly excluding intramolecular folding events.
Each designed siRNA guide strand was aligned to its cognate target
region within the SKP2 mRNA sequence to simulate duplex formation.
Simulations were conducted under default physicochemical conditions
37 °C and 1 M NaCl representing standard ionic strength and temperature
parameters. The output comprised the predicted intermolecular base-pairing
architecture along with the corresponding minimum Gibbs free energy
(ΔG), thereby enabling comparative assessment
of duplex stability and silencing potential across siRNA candidates.
The silencing potential of the designed siRNA guide strands targeting
the full-length SKP2 mRNA transcript was systematically evaluated
using the OligoWalk module. This tool
facilitates quantitative assessment of siRNA–mRNA interactions
by analyzing key determinants of silencing efficacy, including thermodynamic
duplex stability, target site accessibility, and the overall binding
energy landscape. Each siRNA guide strand was individually submitted
under default simulation parameters. OligoWalk computes an efficacy
score expressed as a percentage that reflects the probabilistic likelihood
of successful gene silencing. These scores were subsequently employed
to rank the siRNA candidates, enabling the identification of high-performing
guide strands for downstream functional validation.

2.4
Molecular Docking of siRNA Guide Strands with
Ago2 to Assess Domain-Specific Binding and RISC Compatibility
Molecular docking was carried out to examine how the designed siRNA
candidates interact with the Argonaute 2 (Ago2) protein. The Ago2
structure was obtained from the AlphaFold Protein Structure Database, which uses deep learning to predict protein
conformations. Each siRNA guide strand was modeled using RNA Composer and docked with Ago2 using the HDOCK web server, a hybrid algorithm designed for protein–RNA
interactions.
Each docking simulation was set up by assigning
the Ago2 protein as the receptor and the siRNA guide strand as the
ligand. Docking was performed in RNA–protein mode in HDOCK,
which applies a hybrid strategy that combines template-based modeling
and docking. The algorithm evaluates binding poses using a knowledge-based
scoring function that incorporates shape complementarity, electrostatics,
and statistical potentials derived from known complexes. For each
siRNA-Ago2 pair, the top-ranked docking pose was selected based on
its docking score, which reflects predicted binding affinity and interaction
strength. The resulting complexes were visualized in PyMOL to examine how each siRNA was positioned within
the Ago2 binding pocket, particularly across the PAZ and MID domains
that anchor the guide strand. The most promising complexes were further
analyzed using PLIP (Protein–Ligand Interaction Profiler, 2025
version), which identifies noncovalent
interactions such as hydrogen bonds, hydrophobic contacts, π-stacking,
salt bridges and water-mediated interactions using geometric and physicochemical
rules.

2.5
Target mRNA Accessibility and Hybridization
Analysis
Structural accessibility of the target transcript
was investigated to determine the suitability of candidate siRNAs
for hybridization. Two regions were analyzed: nucleotides 1089–1125,
targeted by siRNA-10, siRNA-11, and siRNA-12, and nucleotides 3001–3023,
targeted by siRNA-92. Because RNA interference efficacy is influenced
not only by siRNA thermodynamic stability and Ago2 affinity but also
by the folding state of the target mRNA, it was essential to characterize
whether these segments were structurally exposed or constrained within
stable secondary structures.
Accessibility was quantified using
RNAup, which implements a thermodynamic
framework for RNA–RNA interaction prediction. The algorithm
decomposes the binding process into two energetic contributions: the
opening energy, representing the cost of destabilizing local secondary
structures in the mRNA, and the duplex formation energy, reflecting
the stabilizing gain achieved through siRNA–mRNA hybridization.
RNAup applies a nearest-neighbor energy model in combination with
dynamic programming recursion to evaluate all possible local structural
contexts. This approach ensures that both the energetic penalty of
structural disruption and the stabilizing effect of duplex formation
are incorporated into the final interaction score.
By integrating
these parameters, RNAup provides a mechanistic measure
of net binding favorability, allowing precise characterization of
whether candidate siRNA sites reside within structurally accessible
regions of the transcript. This algorithmic evaluation establishes
the structural readiness of the selected mRNA segments for effective
siRNA engagement and supports the rational prioritization of siRNA
candidates for downstream validation.

2.6
Dynamic Profiling of siRNA–Ago2 Complexes
Using MD Simulations to Assess Structural Integrity and Silencing
Potential
Molecular dynamics (MD) simulations were conducted
using GROMACS version 2023.1 to assess
and compare the structural stability of selected siRNA–Ago2
complexes. Input files were prepared using the Solution Builder modules
of CHARMM-GUI,
,
a web-based platform that streamlines
the biomolecular system setup across various force fields. The CHARMM36
force field was applied throughout the study.

System preparation began with the PDB Reader module, which
processed the input structures by assigning protonation states, resolving
missing atoms, and converting files into CHARMM-compatible format.
The Solvator module then embedded each complex in a cubic box containing
TIP3P water molecules, with a 10 Å buffer and 0.15 M KCl to simulate
physiological conditions. Ion placement was determined using Monte
Carlo simulations to ensure electroneutrality. CMAP corrections were
included as part of the CHARMM36 parameters. Protonation states of
titratable residues were assigned based on predicted pK
a values at pH 7.0, without manual adjustments. Once solvation
was complete, CHARMM-GUI generated all required input files for GROMACS,
enabling direct execution of energy minimization, equilibration, and
production runs. These files included protocols for steepest descent
minimization, NVT and NPT equilibration, and unrestrained production
dynamics. Each complex underwent a 200 ns production run, a duration
selected to ensure adequate sampling of binding interactions and conformational
behavior, in line with previous Ago2–siRNA studies. Postsimulation
analysis was performed using standard GROMACS version 2023.1 software,
with the modules gmx_rms, gmx_gyrate, gmx_rmsf, and gmx_potentialenergy,
which generate the respective values for RMSD (root-mean-square deviation),
radius of gyration (Rg), and RMSF (root-mean-square fluctuation),
and potential energies were extracted and compared across the siRNA–Ago2
complexes to evaluate their dynamic stability and interaction profiles.
The XMGRACE software and the CHIMERA 1.15 software were used to create the graphs and images, respectively.
The visual inspection of the structures was performed using Chimera
software, which allowed for the detailed visualization and interpretation
of molecular interactions.

3
Results and Discussion
3.1
Retrieval and Preparation of SKP2 mRNA and
AGO2 Structural Data for siRNA Design
NCBI’s Nucleotide
database was used to obtain the FASTA format of the coding DNA sequence
of human SKP2 mRNA, which corresponds to a linear transcript of 3715
nucleotides under accession number NM_005983.4. This sequence served
as the starting point for designing small interfering RNAs aimed at
silencing SKP2 expression, making it essential to model their incorporation
into the RNA-induced silencing complex. For this purpose, the structure
of Argonaute2 (Ago2), which is the central catalytic protein of the
complex, was obtained from the UniProt database with accession Q9UKV8,
encoding 859 amino acids. The AlphaFold-predicted 3D structure complemented
this information and was used to model how siRNA candidates interact
with the PAZ and MID domains that mediate recognition, anchoring,
and cleavage. To ensure the accuracy of the model, residue-wise confidence
scores (pLDDT) were evaluated and visualized, revealing that the PAZ
and MID domains achieved very high scores above 90, while lower scores
below 50 were limited to the N terminal and flexible loop regions.

3.2
Rational Design and Optimization of siRNA
Candidates for Targeted SKP2 Silencing in Breast Cancer
Using
the siDirect 2.0 platform, we generated 127 potential siRNA duplexes
against the human SKP2 mRNA. We then selected 13 duplexes that aligned
with the key principles of three validated siRNA design rules: Ui-Tei,
Reynolds, and Amarzguioui. Refinement proceeded with additional filtersa
GC content of 30–55% and a seed duplex melting temperature
(T
m) ≤21.5 °Cto maximize
silencing efficiency and limit off-target effects. This T
m limit, informed by Ui-Tei et al., minimizes miRNA-like repression by destabilizing the seed
region at physiological temperatures (∼37 °C).
These
13 siRNAs balanced thermodynamic properties for optimal RISC loading
and reduced nonspecific binding, while avoiding the pitfalls of high
GC content that can foster aberrant hybridization. Transcriptome alignments
verified high specificity: Each duplex had at least one partial mismatch
outside its target and no more than a single perfect match elsewhere,
thereby minimizing unintended silencing a vital attribute for experimental
and therapeutic precision.

Supporting Information (Table S1) lists
the target sites, guide and passenger sequences, and complementary
regions for all 127 candidates, while Table
details the 13 finalists, including their
guide and passenger T
m values. By targeting
diverse coding regions across the SKP2 transcript, these siRNAs overcome
potential mRNA secondary structures to ensure reliable knockdown.
This multifaceted approach integrating algorithmic,
thermal, and
specificity criteria provides a robust template for siRNA optimization.
The resulting candidates are well-suited for future SKP2 knockdown
studies in functional genomics and therapeutics, with planned validation
to confirm their efficacy and selectivity.

3.3
In Silico Structural and Thermodynamic Profiling
Identifies Stable siRNA Candidates for Efficient SKP2 Silencing
We performed an integrated computational analysis to evaluate the
functional potential of 13 siRNA duplexes designed to target human
SKP2 mRNA, focusing on key biophysical and thermodynamic parameters
that influence RNA interference efficacy. As summarized in Table
, each duplex was
assessed for guide strand secondary structure, thermodynamic stability,
target-binding affinity, melting temperature profiles, and predicted
silencing efficiency. These metrics collectively inform the duplexes’
structural accessibility, hybridization strength, and resilience under
physiological conditions. By synthesizing these diverse features,
we gained a comprehensive understanding of each siRNA’s suitability
for effective incorporation into the RNA-induced silencing complex
(RISC), highlighting candidates with optimal stability, specificity,
and gene knockdown potential for SKP2 silencing.
A fundamental
starting point in this analysis was the assessment of GC content for
each duplex. Every candidate maintained a GC composition squarely
within the recommended range of 30–55%. This particular window
is crucial because it allows for the formation of stable double-stranded
hybrids with the target mRNA, ensuring precise and efficient recognition
of the SKP2 sequence during the silencing process. At the same time,
it prevents the duplex from becoming overly rigid, which could otherwise
complicate the unwinding step mediated by the RNA-induced silencing
complex (RISC) or increase the risk of unintended interactions with
nontarget transcripts. Certain duplexes particularly exemplified this
optimal tuning, including siRNA 4 with 43% GC, siRNA 10 at 38%, siRNA
11 at 33%, siRNA 12 at 33%, siRNA 15 at 33%, siRNA 20 at 38%, siRNA
92 at 33%, and siRNA 124 at 48%. In these cases, the balanced nucleotide
makeup not only strengthens the reliability of base-pairing interactions
but also supports the smooth extrusion and release of the guide strand
once loaded into RISC, thereby streamlining the overall pathway toward
target cleavage.
Building upon this compositional foundation,
we next quantified
the folding free energies (ΔG) for the single-stranded
guide components of each siRNA. These calculations, performed using
the Oligonucleotide Properties Calculator, produced values ranging
from −22.5 to −26.1 kcal/mol across the set. Such moderately
negative ΔG figures are indicative of a strategically
balanced secondary structure for the guides. On one hand, they provide
the intrinsic stability needed to protect the strands from rapid degradation
by ubiquitous cellular nucleases, which could otherwise shorten their
functional lifespan in vivo. On the other hand, the structures avoid
excessive compactness, preserving the flexibility essential for quick
and seamless incorporation into the RISC complex during the initial
loading phase of RNAi. This dual property resilience paired with adaptability
reinforces the suitability of these duplexes as robust candidates
for practical SKP2 knockdown applications in both research and therapeutic
contexts.
To gain deeper insights into the structural nuances
of the guide
strands, we proceeded with predictions of their secondary conformations
using the MaxExpect algorithm implemented in the RNAstructure web
server. The resulting models, which are depicted in Figure
to visually represent the
predicted folding patterns for each siRNA, revealed a consistent pattern
of low overall stability. In particular, the strands showed only sparse
instances of internal base pairing, accompanied by free energy minima
between 1.6 and 1.9 kcal/mol. This subdued propensity for self-folding
represents a significant advantage in the design of effective siRNAs.
It effectively eliminates the potential for steric clashes or entropic
barriers arising from pronounced intramolecular hydrogen bonds, which
might otherwise impede the efficient partitioning of the guide strand
into RISC. Among the candidates, siRNA 4, siRNA 10, and siRNA 20 were
especially noteworthy in this regard. Their highly flexible and predominantly
linear conformations, as illustrated in the figure, suggest they are
exceptionally well-configured for direct and unhindered engagement
with the RISC machinery, potentially accelerating the transition to
active silencing complexes. During RNA interference, siRNAs are first
incorporated into the RNA-induced silencing complex (RISC) as double-stranded
molecules; however, productive activation of the complex requires
selective unwinding of this duplex. As this unwinding proceeds, the
passenger strand is actively removed, whereas the guide strand remains
bound to Argonaute in a single-stranded state. Once retained, the
guide strand’s intrinsic folding behavior becomes particularly
important, since excessive self-folding can hinder its stable accommodation
within the Argonaute binding channel and disrupt proper seed-region
alignment. Because the passenger strand is rapidly cleaved or discarded
during RISC maturation, its folding characteristics are comparatively
inconsequential. Consequently, evaluating the single-stranded folding
propensity of the guide strand provides mechanistic insight into RISC
activation and, in turn, serves as a meaningful predictor of downstream
gene-silencing efficiency.
−

For greater confidence in these structural predictions,
we conducted
an independent validation using the mfold tool, a widely used algorithm
for RNA folding analysis. This cross-check fully corroborated the
MaxExpect findings, confirming the absence of any prominent or highly
stable secondary elements that could compromise the functionality
of the guide strands. Instead, the siRNAs consistently favored primarily
extended and linear architectures, occasionally interrupted by minor,
low-energy loops or short stem-loop motifs. These subtle structural
features are energetically modest and transient, ensuring that the
thermodynamic landscape remains versatile enough to accommodate the
dynamic requirements of gene silencing, such as rapid hybridization
and dissociation events within the cellular milieu.
Shifting
focus to the critical step of target recognition, we investigated
the formation of duplexes between each guide strand and its corresponding
site on the SKP2 mRNA using the DuplexFold module from the RNAstructure
suite. Figure
provides
a graphical overview of these hybridization simulations, plotting
the free energy landscapes for duplex assembly across all candidates.
The analysis demonstrated that every siRNA formed thermodynamically
favorable hybrids, as evidenced by consistently negative ΔG values that denote spontaneous association and long-lasting
stability. The most thermodynamically robust interactions were those
involving siRNA 4 (−35.7 kcal/mol), siRNA 124 (−35.1
kcal/mol), siRNA 10 (−34.8 kcal/mol), siRNA 20 (−34.2
kcal/mol), siRNA 11 (−33.5 kcal/mol), siRNA 12 (−32.7
kcal/mol), siRNA 15 (−31.2 kcal/mol), and siRNA 92 (−31.1
kcal/mol). These highly exergonic binding profiles highlight a superior
capacity for tight and persistent engagement with the mRNA target,
laying a mechanistic groundwork for potent endonucleolytic cleavage
and subsequent degradation of the SKP2 transcript.
In parallel, we evaluated the thermal resilience
of these siRNA–mRNA
hybrids by deriving detailed melting curves through computations on
the DINAMelt server. The predicted melting temperatures under concentrated
conditions (T
m (conc)) ranged broadly
from 82.9 to 90.9 °C, values that substantially surpass typical
physiological temperatures. This elevated thermal threshold underscores
the inherent durability of the duplexes, affirming their ability to
remain intact and functional within the warmer confines of cellular
environments. It is important to note that while RNAi mechanisms do
not involve heating to these extremes, the high T
m serves as a reliable proxy for strong binding affinity.
At 37 °C, such interactions would thus maintain structural coherence,
enabling dependable target site recognition and precise cleavage by
the RISC-associated Argonaute protein.
Within this cohort, siRNA
4 emerged as the leader in guide strand
thermal stability, with a T
m (conc) of
88.9 °C. It was closely followed by siRNA 20 at 87 °C, siRNA
12 at 85.5 °C, siRNA 124 at 85.3 °C, siRNA 10 at 84.9 °C,
siRNA 11 at 84.4 °C, and siRNA 92 at 83.3 °C. This clear
hierarchy in melting temperatures reflects varying degrees of enhanced
thermal fortitude among the candidates, which correlates with their
potential for prolonged occupancy on the mRNA target and, consequently,
more sustained knockdown effects over time in biological systems.
To cap our evaluation, we forecasted the overall silencing efficacy
of the duplexes using the OligoWalk algorithm, an advanced tool that
holistically integrates binding free energies, the accessibility of
target sites within the mRNA’s native fold, and the thermodynamics
of duplex formation. This multifaceted scoring predicted siRNA 12
as the standout performer, achieving an efficacy of 96.35%. Trailing
just behind were siRNA 11 and siRNA 10, both surpassing the 94% threshold,
a benchmark commonly linked to strong and reproducible knockdown in
cellular assays. Although siRNA 20 and siRNA 80 recorded slightly
more modest predictions, they still resided comfortably within ranges
deemed functionally relevant, making them attractive for scenarios
involving combinatorial therapies or tissue-specific adaptations.
Collectively, this multilayered computational framework successfully
identified eight premier siRNA candidates siRNA 4, siRNA 10, siRNA
11, siRNA 12, siRNA 15, siRNA 92, and siRNA 124 for advancing to experimental
stages. Table
consolidates
key performance metrics for the finalist siRNAs, including GC content,
free energies (Delta G), and various T
m values, facilitating direct comparisons. These metrics provide a
comprehensive assessment of the duplexes’ structural integrity
and predicted silencing efficacy. Specifically, the table details T
m calculations using the Basic, Salt Adjusted,
and Nearest-Neighbor models. Furthermore, it presents T
m variants (Conc and Cp), which characterize the melting
profile under specified strand concentrations. These metrics were
calculated assuming physiological conditions (100 mM Na+). The resulting high T
m values (typically
75–85 °C) confirm robust structural stability at physiological
temperatures (37 °C). These duplexes collectively demonstrate
pliant yet resilient architectures, potent target hybridization capabilities,
steadfast thermodynamic profiles, and elevated projections for silencing
success. The supplementary candidates, namely, siRNA 19, siRNA 40,
siRNA 78, siRNA 80, and siRNA 107, also aligned with several core
benchmarks and could serve as valuable alternatives or adjuncts in
downstream validation studies. Ultimately, this rigorous, data-informed
selection process provides a robust platform for subsequent in vitro
and in vivo testing, with direct implications for developing SKP2-targeted
interventions in oncology and beyond.

3.4
Docking and Interaction Mapping of siRNA Candidates
with Ago2 Reveals Domain-Specific Binding and Structural Anchoring
for RISC Activation
The eight prioritized siRNA candidates
targeting SKP2 mRNA (siRNA 4, siRNA 10, siRNA 11, siRNA 12, siRNA
15, siRNA 20, siRNA 92, and siRNA 124) were subjected to rigid-body
molecular docking simulations against human Argonaute 2 (Ago2; UniProt
Q9UKV8) to quantify their prospective incorporation into the RISC.
Simulations employed the HDOCK web server, a template-based algorithm
that generates binding poses and assigns docking scores (more negative
values indicating higher predicted affinity) alongside confidence
scores (0–1 scale; >0.95 denoting high-probability binding).
As shown in Table
, all candidates achieved docking scores from −305.83 to −387.99,
markedly exceeding the empirical threshold of ∼−200
kcal/mol for stable RNA–protein interfaces derived from crystallographic
benchmarks. siRNA 10, siRN12, siRNA 11, and siRNA 92 had the top four
most negative docking scores of −387.99, −377.05, −343.68,
and −338.75, while each having confidence greater than 95%,
implying preferential guide strand sequestration via 3′-end
clamping in the PAZ domain and 5′-end nucleation in the MID
domain, thereby optimizing slicer activation.
Resultant Ago2–siRNA complexes were rendered
in PyMOL for
qualitative inspection of guide strand topology within the central
nucleic acid-binding channel (Figure
). This revealed conformational distortions in siRNA
4, siRNA 15, siRNA 20, and siRNA 124, manifesting as aberrant intramolecular
base pairing or helical perturbations that misaligned the seed region
(nucleotides 2–8) proximal to the PIWI active site. Consequently,
these were deprioritized due to projected inefficiencies in Argonaute-mediated
unwinding and target scanning. The retained candidates (siRNA 10,
siRNA 11, siRNA 12, and siRNA 92) exhibited fidelity to the canonical
A-form helix, warranting quantitative dissection of noncovalent interactions
via the Protein–Ligand Interaction Profiler (PLIP).
PLIP interrogation delineated residue-specific
contacts across
Ago2 domains (N-terminal, PAZ, MID, and PIWI), revealing domain-distributed
stabilization motifs. siRNA 10, with the most favorable docking score
(−387.99 kcal/mol), engaged 12 hydrogen bonds in the N-terminal
domain predominantly with Asn42, Asn197, Asn200, His47, Ser29, and
Ser208 supplemented by two salt bridges (Asp–Glu pairs) and
three cation−π interactions (Arg/Lys-π stacks).
MID-domain hydrogen bonds (e.g., Gln824 and Tyr529) and PIWI contacts
(Gly953 and Asp1001) further buttressed allosteric clamping, as quantified
by a total interaction energy of −45.2 kcal/mol.
Analogously,
siRNA 11 (−343.68 kcal/mol) accrued 10 N-terminal/MID
hydrogen bonds (Asn42, Asn197, Asn200, Ser29, Ser208, Asn355, and
Asn351), three salt bridges (Lys350–Arg437 network), and eight
PIWI hydrogen bonds, yielding a cumulative −41.8 kcal/mol and
evincing broad-domain cooperativity.
siRNA 12 (−377.05
kcal/mol) emphasized N-terminal dominance
(nine hydrogen bonds: Asn42, Asn197, Asn200, Ser29, Ser208, and His47),
with attenuated PAZ/MID engagements but compensatory PIWI salt bridges
and hydrogen bonds (−39.1 kcal/mol total), suggestive of a
PIWI-biased loading trajectory.
siRNA 92 (−338.75 kcal/mol)
distributed contacts equitably:
seven N-terminal hydrogen bonds (Asn42, Ser29, Asn197), four PAZ (Gly232–Lys236),
five MID (Asn351–Lys350), and six PIWI (Arg635, Glu637, Lys660)
hydrogen bonds, plus dual-domain salt bridges (−37.4 kcal/mol). Figure
focalizes PIWI interactions,
underscoring their catalytic relevance via proximity to the Mg2+-coordinated slicer triad (Asp579, Asp593, Asp767).
These profiles collectively substantiate high-fidelity
Ago2 binding,
with siRNA 10 and siRNA 11 privileging delocalized networks for enhanced
kinetic stability, versus the more modular affinities of siRNA 12
and siRNA 92. Table
corroborates terminus fidelity: all four aligned the guide 3′/5′
ends within 2.5 Å RMSD of PAZ/MID crystal templates (PDB 4W5O), ensuring positional
accuracy for RISC maturation.
Owing to their synergistic attributes superior docking
energetics,
domain-spanning noncovalent arrays (hydrogen bonds, salt bridges),
and seed-region optimality siRNA 10, siRNA 11, siRNA 12, and siRNA
92 advanced to all-atom molecular dynamics simulations (200 ns trajectories)
for temporal assessment of complex dynamics, thereby refining predictions
of RNAi fidelity.

3.5
siRNA 10 and siRNA 11 Identified as the Most
Effective Candidates for the 1089–1125 mRNA Region
RNAup analysis of the 1089–1125 region showed that siRNA 10
had a total free energy of −25.69 kcal/mol, with a duplex formation
energy of −34.82 kcal/mol and opening costs of 9.13 kcal/mol
for the mRNA and 0.87 kcal/mol for the siRNA. siRNA 11 displayed a
similar profile, with total free energy −25.40 kcal/mol, duplex
energy −33.71 kcal/mol, and opening energies of 7.61 kcal/mol
(mRNA) and 0.69 kcal/mol (siRNA). These results, summarized in Table
, indicate that both
siRNAs combine strong duplex stability with moderate energetic requirements
for mRNA accessibility, which is further illustrated in the Supporting Information (Figure S1)
.

This observation aligns with previous reports emphasizing
the critical
role of target accessibility in RNAi efficacy. In particular, Sfold-based
analyses have demonstrated that the “disruption energy”,
the energetic cost required to unwind local target structures, serves
as a key determinant of knockdown efficiency across diverse siRNAs. In another large-scale study, siRNAs targeting
regions of native mRNAs that were predicted to be unpaired (i.e.,
structurally accessible) achieved stronger silencing than those targeting
highly structured regions. Moreover,
early experimental work using designed mRNA constructs found that
local base pairing at the target site can significantly impair gene
silencing activity.

Taken together,
the thermodynamic profiles of siRNA 10 and siRNA
11 support their designation as the strongest candidates not only
because they form stable hybrids but also because their target region
is sufficiently accessible, aligning with the well-established importance
of mRNA structure in effective RNAi.

3.6
Molecular Dynamics Reveals siRNA 10 and 11
as Structurally Stable and Functionally Compatible RNAi Candidates
for SKP2 Silencing in Breast Cancer
To evaluate the dynamic
stability and interaction behavior of the top-performing siRNA candidates
identified from molecular docking, molecular dynamics (MD) simulations
were carried out using the CHARMM36 force field for a total duration
of 200 ns. Out of eight initially docked siRNAs (siRNA 4, 10, 11,
12, 15, 20, 92, and 124), only siRNA 10, 11, 12, and 92 demonstrated
favorable docking scores and stable binding conformations. All four
selected siRNAs (siRNA 10, 11, 12, and 92) were first subjected to
individual molecular dynamics (MD) simulations to examine their intrinsic
structural stability and dynamic behavior in the solvated environment
prior to complex formation. These simulations were performed for 200
ns using the CHARMM36 force field, and trajectory analyses were conducted
to assess key stability parameters, including root-mean-square deviation
(RMSD), root-mean-square fluctuation (RMSF), radius of gyration (R
g), and potential energy (Table
and Figure
A,B). The objective was to identify which siRNAs retained
the most stable duplex architecture under physiological-like conditions,
as stable siRNA structures are known to exhibit improved persistence,
base-pairing integrity, and silencing efficiency in the RNA-induced
silencing complex (RISC).
,

Throughout the simulation period, siRNA 10 and
siRNA 11 displayed
pronounced structural stability, as reflected by their consistently
low RMSD profiles, averaging 0.645 and 0.658 nm, respectively. In
contrast, siRNA 12 (1.401 nm) and siRNA 92 (0.697 nm) exhibited larger
fluctuations, indicating conformational instability and deviations
from their native duplex conformations. The lower RMSD values observed
for siRNAs 10 and 11 suggest that both retained their A-form helical
geometry with minimal backbone perturbations throughout the trajectory.
This behavior indicates robust internal hydrogen bonding and stable
Watson–Crick base pairing, which are critical determinants
of siRNA functionality and efficient target recognition.
Furthermore,
the RMSF analysis that provides residue-level insights
into atomic flexibility demonstrated that nucleotides in siRNA 10
(average RMSF 0.309 nm) and siRNA 11 (0.264 nm) were significantly
less mobile than those in siRNA 12 (0.758 nm) and siRNA 92 (0.426
nm). This reduced local flexibility indicates enhanced duplex rigidity,
consistent with efficient stacking interactions and minimized strand
fraying. Such structural rigidity enhances siRNA resistance to degradation
and ensures accurate positioning within Argonaute’s guide-binding
groove during RISC loading.
,

The compactness
of each siRNA structure was examined via the radius
of gyration, which measures mass distribution about the molecular
centroid and reflects the degree of structural folding. siRNA 10 (1.67
nm) and siRNA 11 (1.54 nm) showed lower R
g values relative to siRNA 12 (1.93 nm) and siRNA 92 (2.01 nm), confirming
tighter molecular packing and stable tertiary organization. Compact
structures are energetically favorable and resist conformational unwinding,
facilitating improved retention in the silencing complex. Consistently,
potential energy analysis revealed highly negative mean energy values
for siRNA 10 (−2,310,694 kJ/mol) and siRNA 11 (−2,310,611
kJ/mol), indicating thermodynamically stable configurations, whereas
siRNA 12 and 92 occupied higher energy states corresponding to less
stable geometries.
Based on these findings, siRNA 10 and siRNA
11 were selected for
subsequent MD simulations in complex with human Argonaute 2 (Ago2)
to explore their compatibility and stability within the silencing
machinery. The simulations were again performed for 200 ns under identical
conditions to compare the dynamic behavior of Ago2 wild type (WT)
alone and in complex with each siRNA (Table
and Figure
A,B). Argonaute proteins are central to RNA interference,
orchestrating siRNA binding and mRNA cleavage.
−

In the Ago2–WT simulation, the RMSD stabilized
around 0.756
nm, consistent with the intrinsic flexibility of unbound Argonaute
proteins. However, upon siRNA loading, the RMSD values decreased to
0.563 nm for Ago2–siRNA 10 and 0.576 nm for Ago2–siRNA
11, indicating reduced backbone motion and enhanced global stability
of the complexes. This stabilization suggests that siRNA binding constrains
the conformational flexibility of key structural domains. The RMSF
analysis confirmed this observation, showing reduced fluctuations
in the PAZ and MID domains regions critical for 3′-end recognition
and seed-sequence anchoring. The reduced atomic displacement in these
regions implies that both siRNAs established strong intermolecular
interactions that reinforced Argonaute’s structural integrity
and functional configuration.
The R
g profiles of the complexes (3.10
nm for hAgo2-siRNA 10 and 3.19 nm for hAgo2–siRNA 11) were
slightly smaller than that of free hAgo2 (3.33 nm), suggesting that
siRNA engagement induced a more compact, ordered conformation. Such
compactness has been previously associated with catalytically active
Argonaute conformations capable of precise target cleavage. Consistent with this, the potential energy trajectories
revealed that the complexes attained more stable energy states compared
to the unbound protein, reflecting energetically favorable protein–RNA
interactions.
Altogether, the MD simulations demonstrate that
siRNA 10 and siRNA
11 are structurally robust, are thermodynamically stable, and form
highly compatible complexes with hAgo2. These siRNAs exhibited reduced
atomic fluctuations, compact duplex geometry, and favorable interaction
energetics both as free duplexes and in protein-bound form. Their
ability to stabilize Ago2’s conformation supports their potential
as efficient guide strands for RISC-mediated gene silencing. Biologically,
this implies that siRNA 10 and siRNA 11 could effectively target and
silence SKP2. Thus, these candidates emerge as the most promising
siRNAs from our computational analysis, potentially capable of disrupting
SKP2’s oncogenic signaling and attenuating breast cancer proliferation
and drug resistance.

4
Conclusions
SKP2 is a key protein that
helps tag other proteins for modification
inside cells. In breast cancer, it plays two harmful roles: one by
marking tumor-suppressing proteins for destruction (through K48-linked
ubiquitination) and another by switching on cancer-promoting signals
(through K63-linked ubiquitination). These dual actions make SKP2
a major driver of cancer progression, promising biomarker, and therapeutic
target in solid tumors, especially breast cancer. To counter its oncogenic
activity, we developed siRNA-based therapeutics using a high-throughput
computational pipeline guided by established design rules. Starting
with 127 candidate sequences, we applied validated screening parameters
including thermodynamic stability, secondary structure prediction,
and off-target analysis based on our previously published work, to identify eight high-confidence siRNAs.
Subsequent molecular docking with human Argonaute 2 (hAgo2), which
is known to a catalytic core of the RNA-induced silencing complex
(RISC), revealed four siRNAs with binding
confidence scores exceeding 96%, indicating robust compatibility with
the silencing machinery. Importantly, beyond duplex incorporation,
the efficiency of RISC loading is strongly influenced by strand folding
and accessibility. Consistent with prior reports, regions of the guide
strand that remain relatively unstructured facilitate effective engagement
with Ago2, while excessive folding in either strand can hinder duplex
unwinding and bias strand selection. This mechanistic consideration
further supports the high-confidence loading potential of our shortlisted
siRNAs.
,

Domain-resolved interaction analysis
confirmed that siRNA 10 and
siRNA 11 anchored within the MID and PAZ domains of Ago2 regions critical
for guide strand stabilization and RISC loading.
,
To further validate structural integrity and dynamic behavior, we
conducted two tiers of molecular dynamics simulations: first, unbound
siRNAs in a solvated CHARMM36m environment, and second, siRNA–Ago2
complexes compared against wild-type Ago2. Across 200 ns trajectories,
siRNA 10 and siRNA 11 consistently demonstrated superior RMSD, RMSF,
radius of gyration, and potential energy profiles, confirming their
conformational stability and binding resilience. However, this study
focuses on bioinformatics and molecular dynamics–guided siRNA
design. Several targeted delivery strategies for breast cancer have
been reported, such as folate receptor ligands,
,
RGD peptides,
,
aptamer–siRNA chimeras
against HER receptors,
,
antibody–siRNA conjugates,
,
and ligand-decorated liposomal or LNP platforms.
,
These approaches provide translational pathways to advance the designed
siRNAs into in vivo applications. As a next step, we plan to validate
the silencing efficiency of lead candidates and investigate their
delivery in breast cancer models using receptor-targeted ligands and
nanoparticle-based encapsulation systems to strengthen therapeutic
potential.
While the present study focused on sequence selection
and functional
validation, future investigations should extend to the evaluation
of our lead candidates (siRNAs 10 and 11) under clinically established
chemical modification frameworks. Chemical modifications such as 2′-O-methyl and 2′-fluoro substitutions, phosphorothioate
linkages, and GalNAc conjugation have been shown to enhance nuclease
resistance, reduce immunostimulatory responses, and enable tissue-specific
delivery. Integrating these strategies
into subsequent analyses will provide a more accurate assessment of
translational potential and align our findings with the current standards
of oligonucleotide therapeutics. Comprehensive reviews by Egli have
highlighted how structural re-engineering of RNA molecules has transformed
them into clinically approved drugs, underscoring the importance of
coupling sequence design with chemical tailoring.

Collectively, this study introduces siRNA 10 and
siRNA 11 as domain-anchored
RNAi candidates targeting the dual ubiquitination axis of SKP2, advancing
mechanistic insight into SKP2-mediated oncogenic regulation in breast
cancer. Although derived from in silico analyses, these findings provide
a strong foundation for experimental validation. The convergence of
optimized sequence design with emerging siRNA chemical modification
and delivery strategies supports the translational potential of these
candidates as precision RNAi-based therapeutic leads.

Supplementary Material