npj Systems Biology and Applications

Cellular responses to perturbations and drugs are determined by interconnected networks, rather than linear pathways. Individually, the JNK, p38 and p53 stress and DNA-damage response networks are well understood and regulate critical cell-fate decisions, such as apoptosis, in response to many chemotherapeutical agents, such as doxorubicin. To better understand how interactions between these pathways determine the dynamic behaviour of the entire network, we constructed a data-driven mathematical model. This model contains mechanistic details about the kinase cascades that activate JNK, p38, AKT and p53, and free parameters that describe possible interactions between these pathways. Fitting this model to experimental time-course perturbation data (five time-courses with six time-points under five different conditions), identified specific network interactions that can explain the observed network responses. JNK emerged as an important control node. JNK exhibited a positive feedback loop, was tightly controlled by negative feedback and crosstalk from p38 and AKT, respectively, and was the strongest activator of p53. Compared to static network reconstruction methods, such as modular response analysis, the model-based approach identifies biochemical mechanisms and explains the dynamic control of cell signalling.

Treatment with drug combinations carries great promise for personalized therapy. We have previously shown that drug synergies targeting cancer can manually be identified based on a logical framework. We now demonstrate how automated adjustments of model topology and logic equations can greatly reduce the workload traditionally associated with logical model optimization. Our methodology allows the exploration of larger model ensembles that all obey a set of observations. We benchmark synergy predictions against a dataset of 153 targeted drug combinations. We show that well-performing manual models faithfully represent measured biomarker data and that their performance can be outmatched by automated parameterization using a genetic algorithm. The predictive performance of a curated model is strongly affected by simulated curation errors, while data-guided deletion of a small subset of edges can improve prediction quality. With correct topology we find some tolerance to simulated errors in the biomarker calibration data. With our framework we predict the synergy of joint inhibition of PI3K and TAK1, and further substantiate this prediction with observation in cancer cell cultures and in xenograft experiments.

Type 1 diabetes mellitus (T1D) is a systemic disease triggered by a local autoimmune inflammatory reaction in insulin-producing cells that disrupts the glucose-insulin-glucagon system and induces organ-wide, long-term effects on glycolytic and nonglycolytic processes. Mathematical modeling of the whole-body regulatory bihormonal system has helped to identify intervention points to ensure better control of T1D but was limited to a coarse-grained representation of metabolism. To extend the depiction of T1D, we developed a whole-body model using a novel integrative modeling framework that links organ-specific regulation and metabolism. The developed framework allowed the correct prediction of disrupted metabolic processes in T1D, highlighted pathophysiological processes common with neurodegenerative disorders, and suggested calcium channel blockers as potential adjuvants for diabetes control. Additionally, the model predicted the occurrence of insulin-dependent rewiring of interorgan crosstalk. Moreover, a simulation of a population of virtual patients allowed an assessment of the impact of inter and intraindividual variability on insulin treatment and the implications for clinical outcomes. In particular, GLUT4 was suggested as a potential pharmacogenomic regulator of intraindividual insulin efficacy. Taken together, the organ-resolved, dynamic model may pave the way for a better understanding of human pathology and model-based design of precise allopathic strategies.

The pancreatic islet (islet of Langerhans) is a mini-organ comprising several thousand endocrine cells, functioning jointly to maintain normoglycemia. Cellular networks within an islet were shown to influence its function in health and disease, but there are major gaps in our quantitative understanding of such architecture-function relations. Comprehensive modeling of an islet architecture and function requires the integration of vast amounts of information obtained through different experimental and theoretical approaches. To address this challenge, our lab has recently developed Bayesian metamodeling, a general approach for modeling complex systems by integrating heterogeneous input models. Here, we further developed metamodeling and applied it to construct a metamodel of a pancreatic islet. The metamodel relates islet architecture and function by combining a Monte-Carlo model of architecture trained on islet imaging data; and an ordinary differential equations (ODEs) mathematical model of function trained on calcium imaging, hormone imaging, and electrophysiological data. These input models are converted to a standardized statistical representation relying on Probabilistic Graphical Models; coupled by modeling their mutual relations with the physical world; and finally, harmonized through backpropagation. We validate the metamodel using existing data and use it to derive a testable hypothesis regarding the functional effect of varying intercellular connections. Since metamodeling currently requires substantial expert intervention, we also develop an automation tool for converting models to PGMs (step I) using feedforward neural networks. This automation is a first step towards automating the entire metamodeling process, working towards collaborative science through sharing of expertise, resources, data, and models.

The ability to reliably predict and infer cellular responses to environmental exposures would offer a major advance in the investigation of immune regulation in health and disease. One possible approach is the use of in silico modelling. Design of such a mathematical kinetic model would be based on existing knowledge of a biological system and utilise a partial data set to parameterise. However, the process of parameter estimation, key for the accuracy of the model, is difficult to conduct by hand, and thus a computational alternative is necessary. We report the utility of Genetic Algorithm with Rank Selection (GARS) as a parameter estimation tool on multiple biological models, including heat shock, signal transduction via ERK, circadian rhythm and NF{kappa}B systems, where it showed strong accuracy and superiority to the Extended Kalman Filter method, Algebraic Difference Equations, and MATLAB fminsearch approaches. GARS parameter estimation is a valuable tool for biological data because it reliably infers system behaviour from partial data sets, allowing for the prediction of cellular responses to environmental exposures.

The immune system is a complex and dynamic network, crucial for combating infections and maintaining health. Developing a comprehensive digital twin of the immune system requires incorporating essential cellular components and their interactions. This study presents the first blueprint for an immune system digital twin, consisting of a comprehensive and simulatable mechanistic model. It integrates 51 innate and adaptive immune cells, 37 secretory factors, and 11 disease conditions, providing the foundation for developing a multi-scale model. The cellular-level model demonstrates its potential in characterizing immune responses to various single and combinatorial disease conditions. By making the model available in easy-to-use formats directly in the Cell Collective platform, the community can easily and further expand it. This blueprint represents a significant step towards developing general-purpose immune digital twins, with far-reaching implications for the future of digital twin technology in life sciences and healthcare, advancing patient care, and accelerating precision medicine.

The proper balance and transition between cellular quiescence and proliferation are critical to tissue homeostasis, repair, and regeneration. The likelihood of quiescence-to-proliferation transition is inversely correlated with quiescence depth, and deep quiescent cells are less likely to exit quiescence and reenter the cell cycle than shallow quiescent cells. The regulatory mechanisms of quiescence depth are poorly understood but essential for developing strategies against hypo- or hyper-proliferation diseases such as aging and cancer. Our earlier studies have demonstrated that the activation threshold of the bistable Rb-E2F gene network switch (ThE2F) controls quiescence depth. We have also identified coarse- and fine-tuning ThE2F modulators in rat embryonic fibroblasts. To examine whether other quiescent cells (including most adult stem and progenitor cells) under different environmental conditions use the same or different modulators of quiescence depth, here we studied the behaviors of 30,000 theoretical quiescent cell models that each support a functional Rb-E2F bistable switch with a unique parameter set. We found that although the vastly heterogeneous quiescent cell models exhibited no apparent parameter patterns, they converged at two alternative groups of strong quiescence-depth modulators (G1 cyclin/cdk-related and Rb/E2F complex-related). Our further machine learning (decision tree) analysis suggested that the Rb protein level and dephosphorylation rate in quiescent cells determine which modulator group to use to regulate quiescence depth.

Difficult-to-treat (D2T) rheumatic disease affects approximately 12% of rheumatoid arthritis patients and resists sequential biologic therapy, yet no mechanistic model explains this refractoriness as a system-level phenomenon. Here we present the 3-Axis Integrative Framework (3-AIF), a six-variable ordinary differential equation system integrating mucosal tolerance, energy-gated neuroimmune danger sensing, and integrated stress response pathways coupled through Hill-function metabolic gating. Stability analysis reveals bistable dynamics with two co-existing attractors separated by a saddle point. Bifurcation analysis demonstrates fold catastrophe with hysteresis: recovery requires greater therapeutic effort than disease prevention. Sensitivity analysis identifies three dominant parameters mapping to neuroimmune activation, energy drain, and recovery capacity. Cross-disease transcriptomic consistency analysis across six public datasets (n=310, five rheumatic diseases, four tissue types) reveals compartment-specific axis dysregulation -- circulating cells show integrated stress response activation while target tissues show pathway exhaustion -- and disease-specific axis dominance patterns consistent with model predictions.

Macrophages are crucial immune regulators as they can either trigger or resolve inflammation. These properties rely on defined inflammatory states and make macrophages valuable therapeutic targets. Identification of stable steady states and bistability in inflammatory transitions provides deeper insights into immune regulation and facilitates the development of novel therapeutic strategies. We present a multiomic and systems biology approach for the top-down identification of bistable circuits in human macrophages polarized towards pro- and anti-inflammatory phenotypes. Using differential gene expression profiles, we identified three criteria to suspect bistability: two potential attractors, a hysteresis behavior between their transitions, and the presence of potential modules of coregulated genes. This was further confirmed by proteomics data pointing to mutually exclusive and time-dependent profiles of gene expression. By network simplification and creation of a novel pipeline for parameter estimation in bistable models, we obtained a minimal model of inflammatory transitions in which we identified ultrasensitivity and hysteresis. Our minimal model genes establish a regulatory circuit switching miR-155 expression, which in turn regulates the expression of inflammatory marker genes during inflammatory transitions.

Symmetry refers to properties that remain invariant upon mathematical transformations. The principles of symmetry have guided numerous important discoveries in physics and chemistry but not in biology and medicine. Here, we aim to explore the presence of symmetry relationships at the gene expression level as a mean to distinguish between healthy and disease states. We deployed Learning-Based Invariant Feature Engineering - LIFE, a hybrid machine learning approach implemented with two symmetric invariant feature functions (IFFs) to identify Invariant Feature Genes (IFGs), which are gene pairs whose IFF single-value outputs remain invariant across individual samples in a given biological phenotype. Our multiclass classification results across the transcriptomes of 25 normal organs, 25 cancer types, and blood samples obtained from 4 different types of neurodegenerative diseases revealed the presence of unique phenotype-specific IFGs. We constructed networks using these IFGs (IF-Nets) and intriguingly, we demonstrated that the hubs could serve as information encoders, capable of reconstructing sample-wise expression values in relation to their counterpart genes. More importantly, we found that hubs of cancer IF-Nets were enriched with both approved and clinical trial drugs, highlighting "symmetry breaking" as a novel approach for treating diseases.

Sequential drug combinations can significantly enhance therapeutic efficacy, yet the general principles governing when and why sequential administration outperforms concurrent treatment remain poorly understood. While empirical evidence demonstrates that the order and timing of drug exposure can be critical, a mechanistic framework to predict which regulatory architectures are primed for sequential benefit is currently lacking. Here, we systematically enumerated and dynamically analysed 59,040 four-node network topologies to identify the structural design principles that dictate sequential efficacy. Our analysis reveals that only a small fraction of network architectures robustly confer a sequential advantage and identifies a minimal structural requirement for this benefit: a positive feedback loop between the primary drug target and its downstream oncogenic output, coupled with antagonistic crosstalk from a secondary drug target. We demonstrate that this architecture enables bistability, allowing the first drug to reconfigure the network into a suppressed attractor state that is inaccessible through concurrent administration. The treatment schedule determines which of two coexisting stable states the system ultimately occupies, with the gap time between doses defining a critical therapeutic window. Only when the first drug is given sufficient time to displace the system past a threshold does the sequential regimen achieve superior suppression. Our findings establish bistability-enabling network motifs as predictive determinants of sequential drug efficacy and provide a topology-based framework for the rational design of time-dependent combination therapies.

Dose-regimen design requires choosing how much drug to give, when to give it, and how treatment should vary across patients. Mechanistic pharmacokinetic-pharmacodynamic (PK/PD) and quantitative systems pharmacology (QSP) models can predict treatment responses, but optimizing dosing inputs depends on model-specific sensitivity derivations or derivative-free search. Here, we introduce DiffDose, a differentiable programming framework for mechanistic open-loop dose-regimen optimization that uses automatic differentiation (AD) to handle clinically interpretable dose amounts and administration times as differentiable controls. We evaluate our method in three settings: fixed-schedule dose-amplitude optimization in OptiDose PK/PD benchmarks; dose-timing optimization in a chemotherapy-induced neutropenia model with state-dependent delay; and individualized mosunetuzumab dosing in a QSP virtual population. Across these examples, AD produced gradients consistent with references, reduced model-specific derivative work, and shortened benchmark time to solution. DiffDose thereby turns mechanistic PK/PD and QSP models from tools that evaluate prespecified regimens into gradient-based engines for individualized regimen design.

Yes-associated protein (YAP) and its homolog TAZ are transducers of several biochemical and biomechanical signals, serving to integrate multiplexed inputs from the microenvironment into higher-level cellular functions such as proliferation, differentiation, apoptosis, migration, and hemostasis. Emerging evidence suggests that Ca2+ is a key second messenger that closely connects microenvironmental input signals and YAP/TAZ regulation. However, studies that directly modulate Ca2+ have reported contradictory YAP/TAZ responses: In some studies, a reduction in Ca2+ influx increases the activity of YAP/TAZ, while in others, an increase in Ca2+ influx activates YAP/TAZ. Importantly, Ca2+ and YAP/TAZ exhibit distinct spatiotemporal dynamics, making it difficult to unravel their connections from a purely experimental approach. In this study, we developed a network model of Ca2+-mediated YAP/TAZ signaling to investigate how temporal dynamics and crosstalk of signaling pathways interacting with Ca2+ can alter YAP/TAZ response, as observed in experiments. By including six signaling modules (e.g., GPCR, IP3-Ca2+, Kinases, RhoA, F-actin, and Hippo-YAP/TAZ) that interact with Ca2+, we investigated both transient and steady-state cell response to Angiotensin II and thapsigargin stimuli. The model predicts stimuli, Ca2+ transient, and frequency-dependent relationships between Ca2+ and YAP/TAZ primarily mediated by signaling species like cPKC, DAG, CaMKII, and F-actin. Model results illustrate the role of Ca2+ dynamics and CaMKII bistable response in switching the direction of changes in Ca2+-induced YAP/TAZ activity for different stimuli. Frequency-dependent YAP/TAZ response revealed the competition between upstream regulators of LATS1/2, leading to the YAP/TAZ non-monotonic response to periodic GPCR stimulation. This study provides new insights into the underlying mechanisms responsible for the controversial Ca2+-YAP/TAZ relationship observed in experiments.

Assessing metabolic variations and flux quantities enable systematic understandings of metabolic shifts, reprogramming, adaptation and interactions in human diseases. However, omics-based estimation of metabolic flux and its variation remains challenging due to several fundamental limitations: the need for disease and tissue context specific metabolic model; nonlinear enzyme kinetic model that links enzyme and substrate changes to reaction flux; partial, unpaired, and snap-shot measurements across omics modalities; and uncertainty in computational prediction. Here, we present Michaelis-Menten model-based Flux Estimation Analysis (mmFEA), a Monte Carlo framework for estimating condition-specific flux changes by integrating paired or unpaired metabolomics and transcriptomics (or proteomics) data. mmFEA separates each reaction-rate change into enzyme- and substrate-associated components and assess reaction rate using Michaelis-Menten kinetics equation. A baseline metabolite saturation rate is introduced by integrating protein language model predicted kinetic parameters and human baseline level metabolic concentration to enable kinetic-aware integration of unpaired substrate and enzyme level measurements. Distribution of metabolic flux and variations between conditions are further computed using MCMC sampling by treating Michaelis-Menten-derived marginal flux distribution as prior and coherency in flux balance as likelihood. To benchmark mmFEA, we generated an in-house multi-omics data set including transcriptomics, metabolomics, metabolic activity functional assay, and CRISPR screening data using pancreatic cancer cell line system treated by APEX1 inhibitors. We demonstrated that mmFEA could accurately capture experimentally observed metabolic changes and achieved a better performance than all baseline methods. Our analysis revealed the necessity in using both substrate and enzyme modality and kinetic aware model in metabolic flux assessment. Further analysis using independent pancreatic cancer cohorts further validated the robustness of mmFEA, supporting integration of condition-linked unpaired data. Pan-cancer and spatial multi-omic applications demonstrated the use of mmFEA for resolving context-dependent metabolic variation when direct flux measurements are unavailable. Together, mmFEA provides a mechanistically grounded framework for estimating relative metabolic flux changes and their uncertainty from heterogeneous omics data.

Modeling of systems for which data is limited often leads to underdetermined model identification problems, where multiple candidate models are equally adherent to data. In such situations additional optimality criteria are useful in model selection apart from the conventional minimization of error and model complexity. This work presents the attractor landscape as a domain for novel model selection criteria, where the number and location of attractors impact desirability. A set of candidate models describing immune response dynamics to SARS-CoV infection is used as an example for model selection based on features of the attractor landscape. Using this selection criteria, the initial set of 18 models is ranked and reduced to 7 models that have a composite objective value with a p-value < 0.05. Additionally, the impact of pharmacologically induced remolding of the attractor landscape is presented.

Understanding how distinct neuronal subtypes contribute to Alzheimers disease (AD) pathology remains a major challenge. Patient-derived induced pluripotent stem cell (iPSC) studies have shown neuronal subtype-specific molecular and pathological signatures, yet the underlying metabolic shifts driving this selective vulnerability are not completely understood. Here we present iNeuron-GEM, the first manually curated, genome-scale metabolic network of human neurons that integrates transcriptomic and metabolic knowledge to resolve subtype-specific metabolic states. By coupling iNeuron-GEM with single nucleus RNA sequencing data from post-mortem human cohort studies, ROSMAP and SEA-AD, we capture neuronal subtype-specific metabolic features and fluxes and identify perturbations in lipid and energy metabolism across excitatory and inhibitory neurons. Integrative analysis with NPS-AD data shows overlapping metabolic disruptions in AD and schizophrenia (SCZ), suggesting shared molecular vulnerabilities between neurodegenerative and neuropsychiatric disorders. We also developed a computational pipeline to infer transcriptional regulation of metabolic pathways and identify NR6A1 and NR3C1 as important regulators of lipid dysregulation in AD neurons. Our study establishes iNeuron-GEM as a framework to identify neuronal subtype-specific metabolic vulnerabilities in complex brain disorders.

Computational models of signaling networks provide frameworks for predicting how molecular cues guide cell decisions. But they are typically limited by manual curation from incomplete literature. Here, we test whether general-purpose large language models (LLMs) generate accurate models of signaling networks. We find that general purpose LLMs generate 24-58% of the reactions of literature-curated networks for cardiomyocyte hypertrophy, myofibroblast activation, and mechano-signaling, and predicting network responses to perturbations with accuracies of 5-26%. While current general-purpose LLMs generate signaling networks with limited accuracy, this study provides a pipeline and benchmarks to guide future improvements.

Disease reprograms cells through changes in gene regulation, yet identifying these changes remains a major challenge. We introduce NetDes-Duo, a computational method that jointly infers transcription factor regulatory network models for two related conditions using scRNA-seq data. The networks are optimized to have minimal topological differences, while the associated ODE models recapitulate single-cell gene expression trajectories for both conditions. On synthetic benchmarks, NetDes-Duo outperformed methods that infer each network independently. NetDes-Duo was applied to neutrophil reprogramming in naive and tumor-bearing mice, and the network-simulated dynamics reproduced the observed cell state transitions. The naive landscape had two well-separated basins, whereas the tumor-bearing landscape was more continuous, with three shallower basins. Perturbation and driving simulations also identified Cebpb as a key driver of the tumor-bearing transition, consistent with emergency granulopoiesis literature. We expect NetDes-Duo to be a broadly applicable framework for uncovering the regulatory logic of disease-associated cell state transitions.

Although endocrine therapies and Cdk4/6 inhibitors have produced significantly improved outcomes for patients with estrogen receptor positive (ER+) breast cancer, continuous application of these drugs often results in resistance. We hypothesized that cancer cells acquiring drug resistance might increase their dependency on negative regulators of the cell cycle. Therefore, we investigated the effect of inhibiting WEE1 on delaying the development of resistance to palbociclib and fulvestrant. We treated ER+ MCF7 breast cancer cells with palbociclib alternating with a combination of fulvestrant and a WEE1 inhibitor AZD1775 for 12 months. We found that the alternating treatment prevented the development of drug resistance to palbociclib and fulvestrant compared to monotherapies. Furthermore, we developed a mathematical model that can simulate cell proliferation under monotherapy, combination or alternating drug treatments. Finally, we showed that the mathematical model can be used to minimize the number of fulvestrant plus AZD1775 treatment periods while maintaining its efficacy.

Biological networks that are formed through amalgamation of signaling pathways include recurrent configurations called network motifs. These statistically over-represented subgraphs are often formed through interconnected enzyme-substrate relationships that are known to result in highly dynamic downstream behavior, including oscillatory output. Such signals are abundant in biology: heartbeats, circadian rhythms, and cell cycles all exhibit characteristic frequencies. Though there has been great emphasis on how oscillations can be generated through network dynamics, little is known about oscillatory information processing or transduction capacity of network motifs. We employ ordinary differential equations-based dynamical modeling to understand how different network topologies impact oscillatory signal propagation through a multi-enzyme network. We model enzyme-substrate interactions of 20 commonly observed motifs using Michaelis-Menten kinetics. We then perform deterministic Monte Carlo simulations using a range of biologically relevant enzymatic parameters and input frequencies. From these simulations, we quantify signal propagation characteristics using cluster analysis, categorize different motif responses based on output characteristics, and identify potential mechanisms for oscillatory signal processing using parameter sensitivity analysis. We see that the input-output responses depend on network topology and enzyme kinetic parameters. Enzymatic motifs show median oscillatory suppression of 30-135 decibels, with three-node coherent feedforward loops showing the lowest propensity for oscillatory signal suppression. Motifs that contained negative feedback or four-node coherent feedforward loops had the biggest potential to act as AC-to-DC converters, translating oscillatory input signals into transient impulses or sustained continuous outputs, respectively. We conclude that enzyme networks can process and decode information within oscillatory inputs in a frequency- and network-dependent manner.