Precision Therapeutics in Triple-Negative Breast Cancer: Boolean Network Control and Pharmacological Mapping of Individualized Models

1. Introduction

Cancer, particularly breast cancer, remains an escalating global health emergency due to rising incidence rates [1] and pronounced inter- and intra-tumoral heterogeneity [2]. This profound molecular complexity frequently drives therapeutic resistance to standard cytotoxic regimens, highlighting an urgent clinical need to transition toward precision oncology [3]. While high-throughput transcriptomics such as RNA sequencing [4] offer comprehensive profiling of aberrant gene expression, these static “snapshots” are insufficient to unravel the non-linear, dynamic behaviors of Gene Regulatory Networks (GRNs). Translating massive omics datasets into actionable clinical insights requires a paradigm shift toward Systems Biology.

Mathematical modeling of GRNs is essential for predicting pharmacological responses in silico. Continuous-time quantitative frameworks, such as Ordinary Differential Equations (ODEs), provide precise kinetic insights but are severely bottlenecked by the lack of experimental parameters for large-scale networks [5,6]. To circumvent this, qualitative Boolean network modeling [7] offers a highly effective, parameter-free logical framework. By abstracting molecular states into binary variables, Boolean modeling thoroughly explores the system’s state space to identify “attractors”—stable phenotypes like proliferation or apoptosis [8]. Consequently, this scalable framework enables investigators to simulate how perturbing specific logical nodes can force a global phenotypic shift, even without detailed kinetic data [9].

In this study, we applied this methodology to Triple-Negative Breast Cancer (TNBC), a highly aggressive subtype lacking actionable targets (ER, PR, HER2) and characterized by rapid chemoresistance and early metastasis [10]. Building upon the generic breast cancer Boolean GRN modeled by Sgariglia et al. [11], we substantially expanded the network to capture TNBC-specific pathology. Concretely, the baseline topology was integrated with a curated subset of critical regulatory nodes that have historically been implicated in TNBC oncogenesis, maintenance of an aggressive phenotype, and therapeutic evasion. This targeted expansion involved the incorporation of IRF2BP2 [12,13], which acts as a key transcriptional cofactor modulating immune evasion and survival pathways; MAPK4 [14], an atypical mitogen-activated protein kinase known to promote oncogenic signaling independently of the canonical RAF/MEK/ERK cascade; BCL11A [15], a potent transcriptional repressor driving cancer stemness and tumor progression; USP7 [16,17,18], a deubiquitinase crucial for stabilizing key oncogenic substrates and hindering apoptotic signaling; and TRIP13 [19,20,21], a mitotic checkpoint chaperone frequently overexpressed and associated with chromosomal instability and drug resistance in refractory breast malignancies.

Our core methodological innovation lies in the patient-specific contextualization of the network topology and drug therapy inference. By integrating individual TNBC RNA-seq data, we dynamically “pruned” biologically silent interactions. This yielded distinct, highly specific sub-network architectures that capture unique inter-tumoral heterogeneity while maintaining a standardized analytical scaffold across the cohort.

To overcome the computationally unfeasible combinatorial explosion of possible system states within these personalized networks, we introduced a Genetic Algorithm (GA). Utilizing a “state-constraining” approach, the GA efficiently converged upon precise, minimal combinations of genes that, when forcibly fixed to specific logical states, irreversibly disrupt the pathological TNBC equilibrium and guide the network toward the apoptotic attractor.

Finally, to bridge computational abstraction and clinical applicability, we mapped these mathematical targets to existing FDA-approved antineoplastic agents, facilitating a cross-indication drug-repositioning strategy. The resulting combinatorial regimens were prioritized based on simulated apoptotic efficacy, maximized druggability, and minimized structural network fragmentation (as a proxy for reducing systemic functional impact). This integrated pipeline provides a robust, scalable blueprint for the automated discovery of tailored, minimally disruptive targeted therapies in TNBC.

Ultimately, the core hypothesis of this study is that integrating patient-specific transcriptomic profiles with Boolean network dynamics and heuristic optimization can uncover minimal therapeutic interventions to overcome TNBC resistance. Our findings robustly validate this approach: the developed computational pipeline successfully personalized regulatory networks across the clinical cohort, identifying highly specific, multilayered combinations of driver nodes that irreversibly induce an apoptotic transition in 92.7% of the analyzed profiles. This report highlights the potential of Boolean modeling to translate malignant regulatory networks into actionable therapeutic targets, thereby facilitating cancer control through clinically rationalized treatment strategies. In other words, this integrated systems biology framework establishes a scalable paradigm for precision oncology, demonstrating that complex tumor networks can be safely and rationally reprogrammed through targeted, computationally guided polypharmacology.

2. Results

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

2.1. Patient-Specific Network Personalization via Topological Pruning

The primary outcome of our structural customization pipeline was the successful generation of highly individualized Boolean regulatory networks for each TNBC patient. A fundamental achievement of this phase was the strategic decision to elaborate the patient-specific topology by acting selectively on the regulatory interactions (edges) rather than deleting nodes (genes).

Unlike standard reductionist approaches that remove nodes based on low expression levels, which often disrupt the structural integrity of the model by fragmenting the scale-free biological topology [22], our approach preserved the overarching Prior Knowledge Network (PKN) scaffold [23]. By systematically identifying and severing only the specific edges that mathematically contradicted the patient’s empirical transcriptomic profile (e.g., paradoxical activations or defunct repressions), we effectively rewired the flow of regulatory information. This resulted in refined sub-networks that perfectly mirror the unique, pathologically altered signaling cascades of each tumor, without compromising the basal biological architecture.

The formalization of this data-driven filtering mechanism culminated in the development of the Topological Network Pruning algorithm (Algorithm 1). Engineered as a deterministic and scalable computational tool, this algorithm represents a foundational result of our framework. It robustly automates the translation of static, continuous RNA-seq data into coherent, individualized Boolean topologies by systematically evaluating both outgoing and incoming regulatory edges against strict logical coherence thresholds. The complete programmatic logic governing this personalization phase is detailed below.

Algorithm 1. Pseudocode for Patient-Specific Topological Network Pruning. The algorithm systematically customizes the generic Boolean regulatory graph ($G$) by integrating individual transcriptomic data (logFC array) against a defined tolerance threshold ($\tau$). Rather than indiscriminately deleting entire nodes, the procedure evaluates the logical coherence of each regulatory interaction. For every significantly dysregulated node (classified as UP or DOWN), the algorithm iteratively verifies both outgoing (Step 3A) and incoming (Step 3B) edges. Regulatory interactions are selectively removed if their theoretical function (+1 for activation, -1 for inhibition) flagrantly contradicts the empirical expression states of the connected node pair. Conversely, edges involving nodes with missing data or non-significant expression changes are conservatively retained to preserve the basal structural integrity of the resulting personalized network ($G_{pruned}$).

The execution of this algorithm successfully yielded a customized therapeutic network for all evaluated patient samples. Comprehensive results for all individualized model personalizations are available in the Zenodo repository cited in the Data Availability Statement below.

2.2. Identification of Therapeutic Driver Nodes via Hybrid Genetic Optimization

Following topological personalization, our objective was to identify specific Driver Nodes capable of actively reprogramming cellular fate. By “Driver Node”, we intend a node that may influence the behavior of the entire network or a substantial part of it [24]. Unperturbed TNBC Boolean models naturally converge into stable attractors that sustain tumor proliferation. Our framework solves this control problem by discovering the precise combinatorial perturbations necessary to force the network’s trajectory out of this pathological trap and irreversibly into a predefined “Apoptosis Attractor” [25].

To navigate the immense combinatorial search space inherent in large-scale network control, we engineered the Unfiltered Strict Hybrid Optimizer (Algorithm 2). This algorithm integrates deterministic Boolean simulations with a heuristic evolutionary search [26]. By permanently locking candidate nodes into fixed logical states—mimicking continuous pharmacological interventions—the optimizer rigorously evaluates the global phenotypic impact of localized structural constraints. Driven by a biological fitness evaluation, the algorithm efficiently converged on minimal, patient-specific sets of Driver Nodes (evaluating up to 8 simultaneous interventions). These combinations maximized alignment with a strict 26-node apoptotic signature while explicitly rewarding the stabilization of terminal executioner pathways, such as the Caspase cascade. Across the entire cohort, this strategy consistently identified the minimal interventions required to secure the apoptotic phenotype. The complete operational architecture of this targeted genetic search is delineated in the pseudocode of Algorithm 2.

The sequential application of the topological pruning and hybrid optimization algorithms to our clinical cohort of 41 TNBC samples yielded highly robust quantitative outcomes. Overall, the computational pipeline successfully identified minimal combinations of driver nodes capable of irreversibly forcing an apoptotic state transition in 38 out of 41 patient-specific networks, achieving a system control rate of 92.7%.

By analyzing the terminal configuration of specific apoptosis nodes within these dynamically controlled networks, we evaluated the underlying biological mechanisms driving these successful transitions. To structurally contextualize these functional outcomes, Figure 1 delineates the foundational architecture of the two primary apoptotic cascades: the receptor-mediated extrinsic pathway and the mitochondria-mediated intrinsic pathway.

Algorithm 2. Pseudocode for the Unfiltered Strict Hybrid Optimization. This algorithm identifies the minimal set of therapeutic driver nodes required to reprogram the patient-specific network toward an apoptotic phenotype. Following a baseline evaluation of the unperturbed pathological state (Phase 1), it employs an iterative genetic search (Phase 2) to systematically navigate the combinatorial perturbation space. Candidate interventions are evaluated via deterministic Boolean simulations and scored using a dual fitness function that rewards exact target alignment and Caspase cascade stabilization. The evolutionary process utilizes elitism, crossover, and mutation, terminating immediately upon achieving a target match to guarantee therapeutic minimality.

As illustrated, while these signaling trajectories process independent upstream stimuli through distinct initiator nodes—such as the Fas-Associated Death Domain (FADD) for the extrinsic route and the Mitochondrion/Caspase-9 axis for the intrinsic route—they ultimately converge upon a common terminal axis: the activation of executioner caspases required to orchestrate programmed cell death. The genetic optimizer dynamically leveraged these distinct structural routes to overcome patient-specific topological resistance.

By mapping the in silico target configurations onto this canonical biological scaffold, we categorized the individualized therapeutic responses into distinct mechanistic profiles. The overall quantitative distribution of these outcomes across the analyzed cohort is comprehensively summarized in Figure 2, revealing a stratification into three distinct response behaviors based on the underlying network dynamics. Specifically, in the majority of successfully controlled networks (n = 24), an extrinsic pathway dominance was observed; the intrinsic mitochondrial pathway remained fundamentally compromised due to the inability to properly configure the initiator Caspase-9, but the optimizer successfully salvaged the apoptotic cascade by routing programmed cell death exclusively through FADD activation. Conversely, a substantial subset of the cohort (n = 14) exhibited a combined intrinsic and extrinsic activation, where the identified driver nodes achieved a highly robust biological response by simultaneously securing the correct logical configuration of both mitochondria-associated initiators and death-receptor adaptors. Finally, a marginal fraction of the cohort (n = 3) yielded a non-significant response, as severe topological distortion and structural fragmentation of these specific patient-specific models prevented the algorithm from discovering a viable combination of driver nodes capable of forcing the apoptotic attractor.

Following the quantitative assessment of the cohort’s therapeutic apoptotic profiles, we evaluated the qualitative nature of the identified interventions to resolve the specific molecular mechanisms driving network control. To this end, we constructed an integrated target matrix mapping each patient-specific model against its required compositional driver nodes and their specific directional perturbation requirements (Figure 3). This qualitative framework distinguishes between the therapeutic necessity of targeted gene activation (+) versus targeted inhibition (-) within the context of network topology. The structural arrangement of the matrix, organized by response groups, immediately contextualizes molecular interventions against phenotypic outcomes. Specifically, the cohort stratification into three distinct subgroups—Combined Activation, Extrinsic Dominance, and Non-Significant Response—reflects the predominant signaling cascade leveraged by the optimization algorithm to trigger apoptosis, distinguishing between cases requiring dual-pathway engagement versus those driven exclusively by the extrinsic apoptotic module. The upper marginal barplot mirrors the number of genes to be activated or repressed for driving the network toward the apoptosis attractor per patient, which illustrates the precise number of regulatory checkpoints required for simultaneous network control. Concurrently, the right marginal barplot defines the Total Frequency of each driver node across the cohort.

The identified patient-specific driver node configurations are dynamically validated to steer each individualized gene regulatory network toward a stable apoptotic attractor. Consequently, the combination of therapeutic targets identified for each patient dictates a precise configuration of the 26 target nodes utilized to define the cellular phenotype [11], successfully establishing the specific steady-state profile that drives the cell toward programmed cell death.

2.3. Aggregate Topological and Structural Analysis of the Interventions

Following the identification of patient-specific apoptotic trajectories, we evaluated the collective topological properties of the computational interventions across the entire cohort of 41 TNBC samples. This aggregate analysis characterized the structural profile of the specific driver nodes selected by the genetic optimizer, revealing how the algorithm navigated the network topology to enforce the target state transitions.

First, we assessed the relationship between a node’s structural prominence and its actual regulatory efficacy across the cohort (Figure 4a). As illustrated in the distribution, the selected driver nodes are predominantly characterized by low-to-moderate connectivity scores, while highly connected structural hubs were consistently bypassed by the selection pipeline. This outcome aligns with established literature, which demonstrates that network control can be effectively achieved through non-hub nodes with low centrality [24,27].

Furthermore, analyzing the structural proximity of the selected driver nodes to the 26-node apoptotic target module revealed a multi-layered cascading control organization (Figure 4b). The spatial distribution of the interventions demonstrates that the algorithm primarily identified and favored optimal control nodes in the immediate structural vicinity of the target group (Distance 1). However, in network topologies where direct local intervention was dynamically insufficient or logically constrained, the algorithm expanded its search space into progressively deeper upstream layers. Consequently, targets were identified at distances up to a maximum of 4 jumps away. This spatial distribution reflects a functional cascading mechanism, wherein initial distal perturbations are transmitted indirectly through intermediate regulatory layers to ultimately force the downstream apoptotic module into its required configuration.

2.4. Pharmacological Therapeutic Optimization and Drug Repositioning

Following the computational identification and structural validation of patient-specific driver nodes, the final outcome of our analytical pipeline was the successful translation of these abstract mathematical perturbations into clinically actionable pharmacological regimens.

To navigate the complex polypharmacological landscape—where single drugs may hit multiple targets and multiple drugs may target the same gene—we engineered the Multi-Objective Pharmacological Optimization framework (Algorithm 3). This declarative, solver-based algorithm represents a critical translational achievement. By imposing an absolute clinical safety constraint (Hard Constraint), the algorithm automatically invalidates any pharmacological combination whose biological mechanism of action contradicts the requisite Boolean regulatory signals needed to force the apoptotic transition.

The execution of this optimization module yielded highly specific, synergistic drug regimens tailored to each patient’s unique network topology. The system systematically extracted these optimal combinations by satisfying a strict, clinically rationalized hierarchy: it primarily maximized phenotypic efficacy (comprehensive driver node coverage), subsequently minimized the overall drug cardinality to reduce polypharmaceutical toxicity, and finally minimized off-target interference to ensure highly localized intervention and systemic safety. The exact logical architecture and programmatic constraints governing this automated pharmacological selection are delineated in Algorithm 3.

Algorithm 3. The pseudocode details the logical constraint satisfaction framework used to pair mathematical Driver Nodes with clinical drugs. The solver guarantees biological coherence by strictly filtering out contradictory compounds (Hard Constraint). An optimal synergistic regimen is then deterministically derived through a multi-objective hierarchy: primarily maximizing the coverage of required therapeutic targets (Priority 3), secondarily minimizing the total number of administered drugs to reduce complexity (Priority 2), and finally resolving redundancies by minimizing off-target side effects (Priority 1).

The optimization procedure (Algorithm 3) identified seven distinct pharmacological agents capable of executing the required directional perturbations on the prioritized regulatory targets within the individualized networks (Figure 5). Specifically, Cladribine, Diethylstilbestrol, Vincristine Liposomal, Dactinomycin, Dexamethasone, Bexarotene, and Tretinoim were selected as the optimal compounds to engage their respective molecular targets across the patient cohort.

To evaluate the broader pharmacological footprint of these interventions, we quantified the sample-specific Side Effect Score displayed in the upper marginal profile of the matrix (Figure 5). For each analyzed patient, this metric explicitly indicates the total number of molecular targets modulated by the selected drug combination that are distinct from the specific control targets identified by the algorithm. By isolating therapeutic driver engagements from external network perturbations, this score provides an objective measurement of the patient-specific off-target burden associated with each combinatorial intervention.

3. Discussion

The present report introduces and validates an innovative computational framework for the dynamic control of gene regulatory networks in TNBC. The cornerstone of this approach lies in the therapeutic personalization, achieved through a systematic topological pruning procedure [28] applied individually to each analyzed clinical case. This pruning strategy, relying on specific transcriptomic profiles of individual samples, enabled the identification of highly specific therapeutic targets tailored to each patient. By integrating the topological mapping of these personalized networks with a multi-objective genetic search framework, we identified minimal sets of Driver Nodes capable of orchestrating the cascading activation of the apoptotic module. Following this algorithmic identification, the implemented procedure successfully translated these computational targets into actionable therapeutic protocols; wherever feasible, the framework mapped specific pharmacological agents to the selected driver genes, strictly adhering to the precise multi-objective optimization criteria embedded within the algorithm. These findings raise critical methodological and clinical implications that we discuss below.

3.1. Interpretative Synthesis of Results

Beginning with the quantitative efficacy of the framework, the programmed phenotypic transition was successfully completed in 38 out of the 41 analyzed personalized regulatory network models (approximately 92.7%). From a computational perspective, this high success rate validates the intrinsic synergy of the implemented hybrid method, which couples the rigorous deterministic nature of Boolean logic rules with the stochastic exploratory efficiency of the genetic algorithm. While Boolean equations guarantee biological fidelity and adherence to the biochemical constraints of molecular interactions, the stochastic engine effectively navigates the immense combinatorial state space, identifying optimal control configurations that would otherwise remain inaccessible through a purely exhaustive or deterministic search. The identification of the optimal Driver Nodes offers an important partial confirmation of fundamental theories on the control of complex systems, with specific reference to the seminal contributions [24] regarding the structural controllability properties of graphs. Following the results of our model, the nodes endowed with the absolute highest topological connectivity do not necessarily coincide with the most effective control node for inducing the transition toward the desired attractor.

This functional selection can be contextualized by the algorithm’s hierarchical preference for topological proximity, paired with the capacity to penetrate deep into the network architecture when neighboring nodes are logically constrained. Rather than distributing perturbations indiscriminately, the optimization framework operates via a systematic search priority relative to the 26-node apoptotic target module. The simulation data indicate that the algorithm prioritizes a direct control strategy, selecting driver nodes in immediate contact with at least one of the target genes to enforce the required binary configuration. Importantly, when such direct local paths are dynamically insufficient or locked by inhibitory Boolean rules, the framework does not fail; instead, it successfully explores deeper regulatory strata. It navigates upstream to identify nodes with indirect connectivity—operating through a single intermediate layer—and, if further constrained, progressively extends its reach through two or more intermediate regulatory steps. Instead of merely exploiting superficial or peripheral network layers, the algorithm demonstrates an advanced structural penetration capability, resolving complex logical bottlenecks by mobilizing distal upstream regulators.

Within this framework, the evaluation of the implemented interventions highlights a highly selective therapeutic approach tailored to each individualized patient network [22,29]. As a validation step, we conducted a literature review focusing on the most frequently activated node identified in this report.

Analyzing the frequency of the algorithm’s choices highlights which specific genes represent the most reliable target nodes across the entire patient cohort.

A quantitative examination of the optimization outputs reveals that the most frequently selected therapeutic targets across the investigated cohort are TRIP13 and CDKN1A, which were identified as optimal driver nodes in 12 and 9 patient cases, respectively. The recurrence of CDKN1A (encoding the p21 protein) underscores its complex context-dependent role in TNBC pathology [30]. Although p21 traditionally acts as a cell cycle regulator within the nuclear compartment, its elevated expression in tumors deprived of hormone receptors serves as a predictor of decreased overall survival and heightened metastatic recurrence [30]. This clinical aggressiveness is primarily driven by the cytoplasmic translocation of p21, which actively suppresses apoptotic cascades through direct interactions with pro-caspase-3 and ASK1 [30], thereby promoting chemoresistance in refractory TNBC cells. Consequently, the computational selection of CDKN1A as a target aligns with therapeutic strategies aimed at disrupting this cytoplasmic survival mechanism.

Given its central position at the intersection of cell cycle progression, DNA repair, and programmed cell death, the systematic targeted inhibition of CDKN1A highlights a promising, network-guided approach to ablate a pivotal oncogenic evasion mechanism in TNBC.

TRIP13 emerged as the most frequently selected driver node across the analyzed cohort, being identified in 12 independent patient cases. This computational recurrence is consistent with established oncological literature characterizing TRIP13 as a potent oncogenic driver highly overexpressed in breast malignancies and consistently associated with poor clinical prognosis, accelerated tumor progression, and metastatic dissemination [31]. Functionally, the systematic silencing or targeted inhibition of TRIP13 has been demonstrated to impair the viability, proliferation, and migratory capacity of TNBC cells [32]. At the molecular level, suppressing TRIP13 signaling triggers programmed cell death, as evidenced by the biochemical cleavage of poly(ADP-ribose) polymerase (PARP) proteins [32], while simultaneously downregulating the PI3K/AKT signaling cascade to suppress downstream cellular survival mechanisms [31]. Furthermore, TRIP13 plays a critical role in driving the epithelial-to-mesenchymal transition (EMT), a key mechanism underlying the invasive nature of refractory subtypes [33]. The targeted ablation of this node effectively reverses the EMT phenotype by upregulating epithelial markers such as E-cadherin and reducing mesodermal factors, including vimentin and Snail, thereby hindering metastatic potential [33]. Given that high baseline expression of this gene is heavily implicated in conventional drug resistance, including insensitivity to standard PARP inhibitors, the network-guided selection of TRIP13 provides a rational therapeutic strategy to resensitize refractory TNBC phenotypes to standard cytotoxic treatments [32].

Beyond identifying the main target genes, the framework actively matches these driver nodes with specific pharmacological treatments.

Building upon the patient-specific target configurations illustrated in Figure 5, we next examine the mechanistic rationale of the selected therapies against TNBC based on existing medical literature below.

Cladribine was identified from a screening of over 1,800 compounds for its unique ability to target the most aggressive TNBC subtypes, specifically those deficient in the CREB3L1 protein, which represents approximately 75% of these tumor cases [34]. It demonstrates high selective cytotoxicity toward metastatic cancer cells, effectively inducing their death while remaining practically harmless to healthy mammary cells, thereby ensuring low systemic toxicity [34]. As a nucleoside analog, it is metabolically activated by the DCK enzyme; since this enzyme is overexpressed precisely in tumors with the poorest prognoses, TNBC cells selectively accumulate the drug and are consequently eliminated, while normal tissues remain spared [34].

Diethylstilbestrol is not considered an effective treatment for TNBC because this tumor subtype is characterized by the absence of estrogen (ER), progesterone (PR), and HER2 receptors [35]. Unfortunately, it has been identified as a significant risk factor for the development of breast cancer, including cases that lack hormone receptor expression, a condition typical of the TNBC phenotype [36].

Vincristine is a chemotherapy drug derived from an alkaloid that binds to tubulin and prevents its polymerization into microtubules, halting cell division in metaphase and driving TNBC cells toward apoptosis [37]. Its encapsulation in liposomes (such as sphingomyelin and cholesterol) is essential for TNBC treatment as it reduces the severe dose-dependent neurotoxicity of the conventional drug. This system allows the administration of higher doses with higher specificity toward tumor tissue, leveraging the enhanced permeability and retention effect [38]. A pegylated liposomal formulation combining vincristine and doxorubicin drastically increases therapeutic efficacy against TNBC cells in both in vitro and in vivo models, representing a more potent combinatorial strategy compared to the single free drugs [39,40]. TNBC cells frequently develop resistance to vincristine via overexpression of P-glycoprotein (P-gp), which actively effluxes the drug from the cell [39]. The use of liposomes for the co-delivery of vincristine alongside sensitizing agents (such as quinine or quercetin) inhibits P-gp and reduces ATP depletion, restoring tumor sensitivity and overcoming multidrug resistance [40].

Dactinomycin emerged among the most effective (“most sensitive”) compounds against TNBC cell lines, in large-scale screenings of FDA-approved drugs, leading to complete cell death [41]. This antineoplastic antibiotic halts cell proliferation by binding to DNA and inhibiting RNA synthesis (transcription). It acts by stabilizing the tumor suppressor protein P53 through direct binding to MDM2 and MDMX, thereby preventing P53 degradation and triggering programmed cell death (apoptosis) [42]. This drug exhibits a potent synergy with doxorubicin, a standard-of-care treatment for TNBC. This combination drastically increases apoptotic rates and significantly reduces tumor mass compared to single-agent treatments by operating through the P53/PUMA/BAX signaling pathway [43].

Dexamethasone is a synthetic glucocorticoid with strong anti-inflammatory and immunosuppressive effects. In cancer therapy, it is used alongside chemotherapy because it may contribute to reducing inflammation around tumors, preventing nausea and vomiting, as well as improving appetite and well-being [29]. It acts through the activation of the Glucocorticoid Receptor (GR), which serves as a survival mediator within TNBC epithelial tumor cells [43]. GR activation induces the expression of pro-survival genes—such as SGK1, MKP1, and KLF5—which can inhibit chemotherapy-induced apoptosis, raising concerns regarding the potential induction of chemoresistance [43,44]. Despite preclinical data suggesting pro-tumorigenic effects, clinical analyses of TNBC patients indicate that dexamethasone use does not have a statistically significant impact on pathologic complete response rates, recurrence-free survival (RFS), or overall survival [44]. However, high GR expression remains a potent predictive biomarker for poor prognosis and an increased risk of metastasis in the triple-negative subtype [45].

Bexarotene selectively activates retinoid X receptors (RXRs) [46]. Bexarotene derivatives have been developed that exhibit significantly greater antitumor potency than the parent compound, even at very low doses. The most critical feature of these new derivatives is their high selectivity; they specifically recognize and attack tumor cells while leaving healthy mammary tissue virtually unaffected [47].

Tretinoin is a retinoid, a natural derivative of vitamin A [48]. Although TNBC is generally resistant to retinoids, Tretinoin acts as a targeted therapy for a specific molecular niche (approximately 3% of cases) characterized by genetic aberrations in the NOTCH1 gene [49]. In these cells, the constitutive production of N1ICD (the active form of NOTCH1) occurs, rendering the tumor dependent on this pathway for growth [49]. Tretinoin intervenes by reducing NOTCH1 mRNA levels, leading to a decrease in N1ICD protein and subsequent proliferation arrest [49].

Interferon alfa-2b (IFN-α2b) is a recombinant form of interferon alpha, a naturally occurring cytokine produced by immune cells in response to viral infections and other immune stimuli [49]. This biopharmaceutical reduces TNBC cell viability by inducing the expression of several members of the Schlafen (SLFN) gene family, including SLFN5, SLFN11, SLFN12, SLFN12-Like, SLFN13, and SLFN14 [50]. SLFN protein induction is directly linked to decreased tumor survival [51]. IFN-α2b’s action is not mediated exclusively by a single gene but operates through a complex and redundant network of cytotoxic signals; notably, the silencing of SLFN14 completely blocks the effects of IFN-$\alpha 2$ on cell viability [50]. In addition, IFN-$\alpha 2$b promotes the differentiation of monocytes into dendritic cells and drives macrophages toward an antineoplastic M1 phenotype. Furthermore, it enhances tumor immunogenicity by upregulating antigen presentation on the surface of neoplastic cells, rendering TNBC more visible to the immune system [52].

Based on this evidence, we can state that the multi-objective optimization framework successfully identified a diverse cohort of pharmacological agents with distinct, validated roles in the management of TNBC. Among these, Cladribine, Vincristine Liposome, Dactinomycin, and IFN-2b act as highly effective, direct antineoplastic agents capable of selectively inducing apoptosis, overcoming multidrug resistance, and remodeling the tumor microenvironment to halt metastatic progression. Conversely, the inclusion of Tretinoin and Bexarotene demonstrates the model’s capacity to exploit patient-specific molecular niches—such as NOTCH1 aberrations—and leverage enhanced derivatives to achieve precision targeting with minimal systemic toxicity. Finally, the identification of Diethylstilbestrol and Dexamethasone highlights critical nuances within the network selection: rather than offering direct therapeutic efficacy, these two agents reveal conflicting evidence and contradictions in the existing literature regarding their roles in TNBC.

Moreover, the patient-dependent variation in the drugs assigned to target TGFB1 and RXRA (as depicted in Figure 5) underscores the practical efficacy of the methodology. This variation confirms that the implemented algorithm does not merely apply uniform protocols but actively operationalizes its optimization specificities.

3.2. Practical and Translational Implications

The findings emerging from this study extend beyond theoretical advancements in systems biology, offering concrete and immediately relevant perspectives for translational medicine and drug development in TNBC treatment.

The primary and most significant practical implication lies in the framework’s capability to guide the rational design of personalized combinatorial therapies. Triple-negative breast cancer is notoriously refractory to therapeutic interventions [53,54]. Our model drastically reduces the biochemical search space, transforming a potentially infinite combinatorial problem into a circumscribed therapeutic target. From an applied perspective, this enables the mapping of the identified driver nodes onto existing drug databases via drug combination approaches. Instead of administering massive and non-specific combinations of chemotherapeutic agents, this approach allows the prescription of cocktails of low-dose compounds, thereby mitigating the risks associated with polypharmacology while maximizing therapeutic synergy [55,56,57].

Secondly, the discovery of the Multi-Layered Cascading Control mechanism provides a solution for targets traditionally considered “undruggable” [58,59]. Many genes directly responsible for apoptosis are structurally difficult to target with currently available molecules. Demonstrating, through our algorithm, that state transitions can be induced by acting on nodes located at topological distances of 2, 3, or 4 from the target module implies that precision medicine can bypass the need for direct interaction with the gene of interest. Clinically, this enables the use of FDA-approved drugs acting as upstream regulators, leveraging the natural signal propagation within the network to achieve the desired cascading effect on downstream apoptotic targets.

Finally, the translational utility of this approach is realized through its capacity to resolve patient-specific interventions into actionable, optimized pharmacological profiles. For each analyzed TNBC sample, the methodology systematically evaluates the individualized driver configuration and, whenever available, identifies the precise therapeutic agents required to modulate those specific control points. This selection pipeline operates under strict multi-objective optimization criteria designed to maximize the directional efficacy of the intervention and minimize the impact on the patient while accounting for the inherent constraints of the network. By aligning personalized drug availability with topological requirements, the computational framework ensures that the selected multi-drug combinations effectively satisfy the logic rules necessary to transition the personalized network out of pathological stability and firmly into a validated apoptotic attractor state.

3.3. Limitations and Future Prospects

While demonstrating substantial computational efficacy, the proposed methodology exhibits intrinsic limitations related to the input data, specifically regarding the structure of the analyzed gene regulatory network and the limited coverage of the pharmacological database utilized for therapeutic mapping. In the former case, the obtained results are strictly constrained by the specific characteristics of the adopted model, namely the number and type of included genes, the network topology, and the logical rules (Boolean functions) governing its dynamics. This underscores a clear necessity to validate the proposed framework across divergent biological contexts and on networks of larger scale and complexity to rigorously test its scalability and generalizability.

Secondly, a critical bottleneck emerged during the matching phase between the identified targets and existing therapeutics, as a significant number of Driver Nodes selected by the algorithm were not associated with any drug in the employed reference database. Consequently, the utilization of continuously updated pharmacological datasets and the inclusion of newly emerging drug-target interactions constitute crucial factors to fully evaluate the actual clinical and translational efficacy of the implemented framework.

Nevertheless, while these points raise legitimate methodological limitations, they concurrently provide direct and constructive inspiration for the future developmental trajectories of the framework. From the perspective of achieving deeper therapeutic personalization, a highly promising evolution would entail integrating new and more granular relational specifications among the constituent elements of the gene regulatory network. This contextual enrichment could be achieved by leveraging advanced technologies, such as Knowledge Graphs coupled with Large Language Models (LLMs), which would enable the extraction and structuring of a significantly broader and more diversified informational payload drawn directly from the entire biomedical literature [60,61]. This would be achieved simply by modulating and recalibrating the parameters within the fitness function of the implemented genetic algorithm, thereby rendering the system an intrinsically adaptive tool.

4. Materials and Methods

The methodology developed for this study was structured as a modular computational pipeline, engineered to systematically process patient transcriptomic profiles and compute optimized pharmacological interventions. The overall system architecture is depicted in Figure 6 and is organized into five sequential analytical modules.

4.1. Module 1: Triple-Negative Breast Cancer Data Acquisition

To translate the theoretical network topology into a biologically accurate and disease-specific computational model, empirical gene expression data were acquired from the publicly accessible Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) repository. Specifically, the transcriptomic dataset GSE65194 was selected for the present study [62].

To ensure the computational analysis was strictly representative of the target pathology, we exclusively extracted and utilized the transcriptomic profiles of the 41 samples clinically classified as TNBC from this cohort. These 41 tumor-specific profiles constituted the sole empirical input utilized in our pipeline. By leveraging this strictly filtered pathological subpopulation, we ensured that the subsequent topological personalization of the Boolean network accurately reflected the unique molecular signatures and the intrinsic heterogeneity characteristic of the TNBC subtype.

4.2. Module 2: Baseline of the Boolean Gene Regulatory Network

To model the complex regulatory architecture of TNBC, we employed a standard Boolean network formalism. In this qualitative framework, critical molecular entities are represented as binary nodes whose temporal evolution is governed by logical regulatory functions mapping established biological crosstalk.

Network dynamics were simulated utilizing a deterministic synchronous updating scheme. While asynchronous updating can capture varied biological timescales, the synchronous approach was selected for its high computational efficiency in large-scale state-space exploration and its ability to systematically minimize biologically irrelevant transient transitions [62]. Under this deterministic paradigm, the finite state space inevitably converges into invariant stable sets (attractors), which in our oncological framework mathematically encode the terminal cellular phenotypes of interest: uncontrolled tumor proliferation and targeted apoptosis.

To enforce biological realism and maintain dynamic robustness, the network’s regulatory rules were strictly constrained to the algebraic class of Nested Canalizing Functions (NCFs) [63]. Unlike arbitrary logical combinations, NCFs accurately model the hierarchical dominance inherent in genuine genomic interactions, wherein specific “master” regulators decisively dictate downstream node states. Master regulators are proteins, typically transcription factors, that sit at the top of gene regulatory hierarchies. Structurally confining the network to NCFs proactively prevents chaotic dynamic regimes, ensuring an intrinsically ordered system characterized by highly convergent trajectories and a stable, well-defined attractor landscape.

Finally, to translate continuous patient transcriptomic data into discrete initial network states, we applied a statistical discretization pipeline [64]. For each TNBC sample, an individual differential expression profile was generated by subtracting the mean transcriptomic baseline of the healthy control cohort. Assuming a normal distribution of these differential values, a Critical Value Cutoff (CVC) was established at the 95th percentile. Nodes exhibiting over-expression above this stringent CVC on the positive side of the normal distribution centered on zero were initialized in an active Boolean state (1), while all others were set to a basal inactive state (0). The same exercise was performed for the down-regulated genes on the negative side of the normal distribution. This thresholding effectively filters out physiological background noise, initializing the computational topology exclusively with the patient’s statistically significant molecular signature.

4.3. Module 3: Customization of Gene Regulation Network

The integration of generic biological knowledge with patient-specific empirical observations represents a cornerstone of our precision medicine pipeline. While the baseline Boolean network formalizes the universal regulatory topology of the system, it lacks the contextual specialization required to model an individual oncological profile. To resolve this limitation, Module 3 executes a data-driven topological pruning procedure, transforming the generalized network scaffold into a personalized, patient-specific dynamic model.

The fundamental objective of the pruning process is the systematic elimination of regulatory interactions (edges) that exhibit flagrant contradictions with the patient’s actual transcriptomic state. An intuitive representation of this filtering mechanism is illustrated in Figure 7, which delineates the computational processing of a simplified three-node sub-network comprising Gene A (upregulated), Gene B (downregulated), and Gene C (upregulated).

The pruning process is performed through the following steps:

• Forward Discrepancy Filtering: When an upstream source node is pathologically active, its downstream regulatory targets must reflect this state according to the literature-derived interaction type. In Figure 6, Gene A is significantly upregulated ($lo{g}_2\ge 1.0$), establishing an active regulatory command. The activating edge (+1) toward Gene B is removed because Gene B is empirically downregulated, representing a “Broken Activation.” Similarly, the inhibiting edge (-1) toward Gene C is dismantled because Gene C evades repression and remains highly expressed, representing a “Broken Inhibition”;
• Backward Causality and Loss-of-Function Validation: The algorithm interprets severe downregulation as a functional silencing of the corresponding regulatory pathway. When a source node is repressed ($lo{g}_{2}FC\ge -1.0$), it cannot logically serve as the primary driver for a highly expressed target. In Figure 2, the activating edge from Gene B (downregulated) to Gene C (upregulated) is stripped. Because the known activator is pathologically repressed, the target’s upregulation must be sustained by alternative, unmapped molecular pathways. Retaining this interaction would introduce false mathematical dependencies into the dynamic model.

Systemically, the pipeline maps the continuous logFC matrix onto the network graph, isolating dysregulated genes based on the predefined CVC. Interactions between non-significant nodes are preserved to maintain basal homeostasis. The algorithm then evaluates every network edge for biological coherence by applying a tolerance zone of $[-CVC,+CVC]$, it accommodates minor physiological noise while strictly eliminating logical paradoxes. Additionally, interactions lacking complete transcriptomic visibility for both adjacent nodes are discarded to prevent the propagation of unverified signals. The comprehensive programmatic workflow governing this customization phase is detailed in Algorithm 1 (see Results section).

4.4. Module 4: Driver Node Search and Network Control

Having successfully tailored the generic gene regulatory network to the specific transcriptomic profile of the TNBC patient, the subsequent and most clinically relevant objective of our computational pipeline is the identification of therapeutic targets. To quantitatively measure the success of the therapeutic optimization, it was imperative to establish a rigorous baseline representing the desired biological endpoint.

The therapeutic objective was defined as the irreversible transition of the network into an apoptotic steady state. This state was mathematically codified by a predefined configuration of 26 critical target nodes (genes) intricately involved in the cell death machinery. As detailed in Figure 8, these 26 critical genes constitute a specialized apoptotic sub-module extracted directly from the broader 136-node global network analyzed in this study. The selection of this specific subset, alongside the precise Boolean activation or repression states required to successfully trigger the apoptotic cascade, was adopted from a previously validated methodological framework [11].

Having defined the exact therapeutic destination through the 26-node apoptotic configuration (Figure 8), the framework must subsequently determine the optimal topological interventions required to force the patient’s personalized network into this desired state. This complex transition is achieved through the computational identification and manipulation of Driver Nodes.

As conceptualized in Figure 9, the dynamic behavior of the Boolean network can be visualized as an attractor landscape [65]. The unperturbed, patient-specific network naturally settles into the “Cancer Attractor,” a stable dynamic basin representing the active disease state and its underlying pathological configuration. To achieve complete state space reprogramming, a targeted computational intervention is necessary. By applying a sustained perturbation force to a highly specific, minimal subset of upstream control genes-defined as the Driver Nodes $(D_1,....,D_K)$[24], the system state is dynamically driven out of the pathological basin. This strategic perturbation overrides the basal regulatory logic, irreversibly redirecting the network’s trajectory toward the “Apoptosis Attractor” and successfully locking the downstream 26 target nodes into their requisite therapeutic configuration.

The methodology is grounded in the evaluation of a specific therapeutic phenotype. Let $T=\{t_1,t_2,....,t_m\}$ be the predetermined set of $m=26$ target genes associated with the apoptotic response, and $V_T=\{v_1,v_2,....v_m\}$be their corresponding desired Boolean steady-state configuration (where $v_i\in \left\{\mathrm{0,1}\right\}$).

The control problem consists of identifying a minimal subset of Driver Nodes $D$ (strictly excluding the target nodes themselves, such that $D\cap T=0$) and determining a corresponding control vector $U_D$. By applying $U_D$ as constant, unyielding external inputs (perturbations), the network’s dynamics are forced to converge toward a steady-state attractor $A$ where the target condition is strictly satisfied according to formula 1.

$$\forall t_i\in T,A\left(t_i\right)=v_i$$

(1)

Initially, the algorithm performs a baseline evaluation of the natural, unperturbed network. Utilizing the patient’s transcriptomic profile as the initial state vector $S\left(0\right)$, the network evolves deterministically toward its natural attractor using a strict synchronous updating scheme following formula 2.

$$S(t+1)=F\left(S\right(t\left)\right)$$

(2)

To systematically discover the optimal minimal set of Driver Nodes $\left(D\right)$ capable of reprogramming the patient-specific TNBC network from malignancy to apoptosis, we developed an Unfiltered Strict Hybrid Optimizer. This computational approach departs from purely algebraic space-reduction techniques by leveraging an iterative genetic search algorithm, heavily constrained by rigorous, deterministic Boolean simulations.

Unlike heuristic models that may tolerate minor update contradictions, our simulator operates under strict logical fidelity; any undefined regulatory rule immediately raises a fatal exception, guaranteeing absolute topological integrity. The resulting baseline attractor is then evaluated against the target configuration $V_T$ to establish an initial comparative reference score. Because exhaustive evaluation is computationally intractable, to navigate the immense combinatorial search space of possible perturbations, the optimizer employs an iterative genetic algorithm [66], scaling the intervention size $K=\left|D\right|$ from 1 up to a maximum predefined threshold (e.g., max perturbations = 8). For each level of $K$ a population of $N$ candidate interventions is randomly initialized from the search space. Each individual (candidate control vector) undergoes rigorous simulation until its corresponding attractor is reached. The performance of each candidate is quantified using a predefined fitness function ($\Phi$), designed to reward both overall stability and the activation of critical biological pathways. The primary fitness metric is the raw count of target genes successfully matching their desired state, as described in formula 3

$${\Phi }_{raw}(D,U_D)={\sum }_{i=1}^{m}I\left(A\right(t_i)=v_i)$$

(3)

where $I$ is the indicator function. Furthermore, to biologically prioritize the execution of the terminal executioner pathways, a specific weighted bonus is appended for the correct stabilization of the Caspase cascade (e.g., CASP3, CASP7, CASP8, CASP9), effectively defining the total fitness score (${\Phi }_{total}$).

Inspired by the principles of Trap Space analysis—where subsets of nodes are irreversibly locked into specific states to restrict the dynamic trajectory—the optimizer applies the candidate interventions $U_D$ as permanent, unalterable Boolean state constraints during the simulation. The evolutionary progression iteratively refines the population through crossover and mutation. At each generation, top-performing candidates are recombined to inherit successful gene-perturbation pairs. Random point mutations (such as inverting a perturbation state or swapping a targeted gene) are introduced at a rate of 20% to prevent premature convergence on local maxima. This iterative deepening continues until an absolute stable state is achieved $({\Phi }_{raw}=m)$, definitively identifying the minimal combination of Driver Nodes required to lock the specific patient’s network into the therapeutic apoptotic phenotype.

4.5. Module 5: Pharmacological Therapeutic Optimization

Translating the mathematically derived Driver Nodes into clinically actionable therapies requires mapping topological perturbations to specific pharmacological compounds. To ensure dynamic compatibility, we developed an automated data enrichment pipeline utilizing the ChEMBL web API (https://www.ebi.ac.uk/chembl/api/data/mechanism) to retrieve the precise Mechanism of Action (MoA) for established oncological drugs [67] (File S1). Compounds were strictly codified to match our network logic: activating agents (e.g., agonists) were mapped to a positive Boolean perturbation (1), and inhibiting agents (e.g., antagonists) to a negative perturbation (0). Ambiguous or uncharacterized compounds were explicitly excluded, yielding a high-confidence, directionally defined drug-target database.

To navigate the resulting polypharmacological landscape—where drugs exhibit multiple targets and targets respond to multiple drugs [68]—we modeled the therapeutic selection as a Constraint Satisfaction Problem (CSP), solved via Answer Set Programming (ASP) [69]. This declarative solver identifies the globally optimal, synergistic drug combination by rigorously evaluating a predefined multi-objective hierarchy consisting of:

• Directional Consistency (absolute constraint): The algorithm enforces a zero-tolerance policy against biological contradictions. A drug is selected only if its pharmacological effect strictly matches the requisite Boolean state of the target Driver Node;
• Phenotypic Maximization (primary objective): If total target coverage is pharmacologically precluded, the solver minimizes the number of missed Driver Nodes, ensuring the closest possible approximation to the ideal apoptotic phenotype;
• Minimal Cardinality (Secondary Objective): To mitigate polypharmaceutical toxicity and unpredictable drug interactions, the system minimizes the total number of administered drugs, prioritizing the simplest effective regimen [70];
• Minimization of Side Effects (Tertiary Objective): Among equally effective minimal regimens, the solver selects the configuration that inadvertently affects the fewest non-Driver nodes, ensuring highly localized intervention and reducing systemic collateral impact.

By computationally satisfying these cascading priorities, the ASP framework automatically generates a targeted, highly specific, and clinically rationalized drug regimen designed to force the patient’s network into the therapeutic apoptotic attractor. The exact programmatic constraints governing this combinatorial selection are detailed in Algorithm 3.

5. Conclusions

This study validates a three-stage computational framework for precision therapeutics in TNBC, designed to drive a phenotypic transition toward apoptosis. The pipeline begins with the ex-ante personalization of gene regulatory networks utilizing individual patient transcriptomic profiles. Subsequently, a hybrid GA, constrained by Boolean determinism, identifies minimal Driver Nodes. This phase revealed a distinct divergence between node centrality and actual dynamic control efficacy and a systemic preference for multi-layered cascading control. Finally, these theoretical targets are mapped to clinically actionable therapies through a rigorous optimization procedure that minimizes the prescription regimen and prioritizes the rational selection of established oncological therapies.

Translationally, this framework suggests that cellular fate can be reprogrammed [71] through highly targeted, non-destructive topological perturbations. This approach facilitates minimally disruptive interventions capable of engaging critical cell death pathways while actively mitigating systemic toxicity.

While highly effective, the current model exhibits limitations inherent to static network topologies and the restricted coverage of existing pharmacological databases. Future enhancements will integrate Knowledge Graphs and LLMs to dynamically extract updated relational specifications from the biomedical literature. By natively incorporating this evidence to recalibrate the genetic algorithm’s fitness function, the system will evolve into an intrinsically adaptive tool.

Ultimately, the synergistic integration of network personalization, genetic optimization, and pharmacological constraint satisfaction establishes a rigorous and scalable platform for the safe, patient-specific therapeutic reprogramming of TNBC.