Gliomas are the most common type of brain tumor in adults and are associated with poor prognosis and high mortality. Despite technological advances, their classification remains challenging in both clinical practice and research, highlighting the need for improved molecular subtyping to enhance diagnosis and treatment. In this study, we systematically evaluate the robustness of current glioma classification using a reproducible pipeline for unsupervised patient stratification. We apply multiple multi-view clustering methods, namely, CIMLR, intNMF, and moCluster, to cluster glioma patients based on transcriptomics and methylomics data. All methods consistently distinguished glioblastoma (GBM) from lower-grade gliomas, with survival analysis confirming poorer outcomes for GBM clusters, in line with clinical expectations. The best partition was obtained with CIMLR, closely aligning with the 2021 WHO classification of central nervous system tumors. In contrast, intNMF identified three clusters with distinct survival distributions, suggesting additional biological heterogeneity that warrants further investigation. To enhance robustness and better capture tumor complexity, we implemented consensus clustering to integrate results across methods. This integrative approach yielded well-defined groups with improved concordance to established glioma subtypes compared to individual methods. Notably, the consensus solution identified more clusters than the traditional classification, revealing deeper molecular characterization and providing a framework to refine glioma classification. Furthermore, we identified key molecular features supported by the literature, including both established and emerging biomarkers. Overall, our findings underscore the value of multi-omics integration and consensus clustering in refining glioma classification, supporting the discovery of novel biomarkers to improve diagnostics, prognostics, and patient outcomes.
Kappa is a graph-rewriting language originally developed for modeling molecular interactions in cellular processes. A key feature of Kappa models is that, inspired by organic chemistry, interaction rules operate on patterns i.e., partially specified molecular species, thus allowing compact descriptions of otherwise large or even infinite models. In this paper, we propose and implement a framework for statistical model checking for stochastic Kappa models against properties written in bounded linear temporal logic (BLTL). A key feature of our approach is that temporal properties operate over Kappa patterns, making the specification language native to Kappa and avoiding the expensive—and sometimes theoretically impossible—translation of the model into an equivalent chemical reaction network with a finite number of species. Concretely, given a Kappa model and a property, we first instrument the model by introducing additional variables and observables that track the truth value of each atomic proposition. For each simulated trace, the BLTL formula satisfaction is evaluated using an offline monitoring procedure. Finally, the satisfaction probability is statistically estimated from repeated model simulations. The proposed model checking framework enables a systematic exploration of behavioral properties of Kappa models, that is computationally efficient while able to capture stochastic and finite size effects with any predefined desired accuracy/precision. As such, it can serve in applications for e.g. model selection, property-driven parameter exploration, or robustness analysis. The framework is illustrated on representative case studies from systems biology and swarm robotics.
In this work, we study the source-target control problem of Boolean networks, which has important applications for cellular reprogramming. More specifically, we want to perturb the input Boolean network such that it leaves a predetermined source attractor to eventually find itself in the target attractor. We use edge perturbations, which modify the update functions of a Boolean network without necessarily setting them to constants. This work improves an existing method in which edge perturbations are used to eliminate all but the target attractor. Both approaches are based on Thomas’ first rule: the existence of at least one positive cycle in the interaction graph of a dynamical system is a necessary condition for the existence of multiple steady states. Thus, if all positive cycles are removed, the resulting dynamical system has a single steady state. This requires removing a subset of edges so that at least one edge of every cycle in the graph is perturbed: the feedback edge set. Information on the source and target attractors allows reducing the number of relevant feedback edge sets and, subsequently, candidate perturbations. We propose a carefully designed strategy that combines exhaustive search and a heuristic inspired by Simulated Annealing to effectively reduce the number of required perturbations with respect to the feedback edge set. We evaluate our method on a wide array of Boolean networks from the literature to demonstrate its efficacy and efficiency.
Aging and Alzheimer’s disease are complex multifactorial processes driven by interacting genetic and metabolic mechanisms, for which single-gene approaches often fail to capture system-level behaviour. Systems biology addresses this challenge by integrating omics data with genome-scale metabolic models (GSMMs) to identify candidate interventions. The Metabolic Transformation Algorithm (MTA) and its robust variant rMTA use constraint-based modelling to predict perturbations that shift a system from a disease-associated toward a healthier metabolic state. To move beyond single-gene prioritization, we applied EA-rMTA, an evolutionary extension of rMTA that searches the combinatorial knockout space by optimizing the robust Transformation Score. We evaluated the framework in two case studies: lifespan-extending unc-62 knockdown in Caenorhabditis elegans and early-onset Alzheimer’s disease (EOAD) in human cortex. In both settings, rMTA found biologically coherent targets and pathway-level signatures consistent with known hallmarks of aging and neurodegeneration. EA-rMTA further identified compact and interpretable multi-gene intervention strategies that outperformed single-gene deletions in shifting the metabolic state toward the desired phenotype. Together, these results show that evolutionary search extends metabolic transformation analysis beyond single-gene perturbations and enables the discovery of experimentally tractable combinatorial intervention hypotheses. More broadly, the framework provides a computational strategy for generating systems-level therapeutic hypotheses in aging, neurodegeneration, and other multifactorial disorders.
Gene Regulatory Networks (GRNs) constitute a useful abstraction of complex molecular mechanisms responsible for cell behaviors, cell differentiation, and various diseases. Advances in high-throughput transcriptomics, such as single cell RNA sequencing data, have enabled unprecedented access to large-scale gene expression data. However, the manual reconstruction of GRNs from such data and literature is a hard task, and the automatic reconstruction by GRN inference algorithms remains a major challenge to analyze and validate. Previous work revealed that the performances of GRN inference algorithms vary significantly according to network topology, degree distribution and motif preponderance. The analysis of sensitivity to those properties requires the generation of synthetic GRNs of various forms, and expression data obtained by simulation. Since the seminal works of Barabási and Albert on the scale-free power-law distribution of connectivity degrees in GRNs, and of Alon and Milo on the particular distributions and significance of small network motifs, various random graph models have been developed to try to capture those particular features. In this paper, we show that a relaxed version of the directed configuration model (RDCM) does allow us to generate random GRNs fitting the graph properties and motif distributions of several large GRNs of the literature: Human-TRRUST, h-ESC, m-ESC, m-DC, and to a lesser extent Yeast and E. coli, thanks to the variance of the results, and despite the discrepancies previously observed in terms of the mean in that random graph model. This is shown with GRNgen, a Python package which takes as input the number of nodes with for each node its in-degree and out-degree, the average path length, diameter, number of arcs and 7 important motif counts, and generates as output a set of random GRNs ranked by their fit to the input properties. We present our results for the generation of 1,000 samples for each of our reference GRNs, and provide elements of comparison with other generators.
Microbial communities play a central role in many bioprocesses with key applications in food fermentation, waste treatment, human and animal well-being, plant protection or metabolite transformation in industrial bioprocesses. However, the metabolic microbial interactions driving the community dynamics remain difficult to characterize because of their complexity and their temporal variability. Recent advances in sequencing and analytical technologies now provide time-resolved multi-omics data at the community scale providing key insights into the mechanisms shaping the community dynamics. However, integrating these heterogeneous data in an interpretable way to decipher species-specific metabolic activity and microbial interactions remains a major challenge in the study of microbial communities. We introduce the community metabolic flux inference (comFI) method, a mathematical framework for inferring the metabolic fluxes of individual microorganisms from community-level longitudinal data. The method formulates flux estimation as a biology-informed constrained inference problem that combines observed microbial abundances and extracellular metabolite exchange data, with metabolic constraints encoded in a metabolic model and transcriptomic-based lasso regularization terms. We evaluated comFI on synthetic datasets generated from dynamic models of microbial communities involving three Escherichia coli mutant strains. The comFI method showed a very good reconstruction accuracy for exchange fluxes, intracellular metabolic fluxes distribution, metabolic pathway activation patterns and strain contribution. We also applied the method to experimental cheese fermentation data involving three bacteria (Lactococcus lactis, Lactobacillus plantarum, and Propionibacterium freudenreichii), combining abundance measurements, targeted metabolomics and metatranscriptomics data. The comFI framework enabled to recover previously identified interaction patterns, and to reconstruct latent intracellular flux states for individual microorganisms alongside with their respective metabolic contributions within the community, consistently with the omics data and known physiology. All together, we demonstrate that comFI provides a practical framework for recovering the metabolic activity of individual microorganisms from community-scale multi-omics time-resolved data.
Kappa offers a modeling environment to describe, simulate, and reason about rule-based models. It has been used to model protein-protein interaction networks, especially models of signaling pathways. Kappa comes with a static analyzer, KaSa, to assist the modeler and assess the consistency of the models. Although efficient, KaSa is sometimes too slow to reason while modifying large models or when editing models within the user interface. Here, we propose an incremental version that updates the result of the current analysis at each model modification. Our approach relies on the use of an abstraction of the relationships between the rules of the model and the properties that they induce. Partial evaluation is used when some rules are removed, to exclude the results that derive from the removed rules. Adding rules is done classically by resuming the iterations of the analysis algorithm. This incremental analysis is available on the command-line or as an electron app, and it is evaluated on examples from the literature.
Reaction Systems (RSs) provide a successful qualitative modelling framework inspired by biochemical reactions. In a RS a computation starts from an initial state given by a set of entities and each following computation state is determined by the application of all the enabled reactions to the previous state. RSs can also model the interaction with the environment. Each entity can either be present or absent in a computation state, as a crisp boolean condition, and also reactions are (or not) enabled under crisp conditions. This framework has proved to have many applications for modelling biomedical and computer science systems, but it can become restrictive when laboratory measurements exhibit graded concentrations, partial inhibition, and noise. We thus introduce Mamdani-driven Fuzzy Reaction Systems (M-FRS), as a graded conservative extension of RSs. Each reaction in the style of RSs is now interpreted as a Mamdani rule, and we formalise a single four-stage fuzzy inference cycle (fuzzification, rule evaluation, aggregation, optional defuzzification) which defines a deterministic discrete-time graded update operator. Fuzzy inference yields a discrete dynamical system. As a first case study, we develop a compact M-FRS model of the hypothalamic–pituitary–thyroid axis.
Predicting gene expression dynamics is challenging due to the complex regulatory interactions within high-dimensional datasets. We evaluate predictive models that integrate temporal patterns with gene–gene networks, comparing a state-of-the-art approach based on Protein-Protein Interaction (PPI) networks from STRING with models utilizing data-driven network inference. Our results show that inferred networks can enhance accuracy over static biological priors. However, simpler models treating genes independently often achieve comparable performance. This suggests that for the considered datasets, the added complexity of explicit gene–gene interactions does not always translate into superior predictive power, opening to further investigations on the most effective ways to represent and leverage biological connectivity in forecasting tasks.
Bursty transcription in single cells typically produces over-dispersed, skewed, and sometimes heavy-tailed expression distributions that are explained by two-state Markov models of the promoters. While the gold standard for simulation is exact stochastic sampling with Gillespie’s algorithm, obtaining thousands of timed traces is computationally costly. Surrogate models based on stochastic differential equations (SDEs) are widely used to speed up this simulation process. An example is the Chemical Langevin Equation based on Gaussian noise, which, however, does not capture heavy-tailed noise. In this work, we present a unified SDE framework that combines deterministic drift, Gaussian fluctuations, and additive sporadic jumps of arbitrary distributions, and provide an open-source Python implementation, bcrnnoise. The framework subsumes standard surrogate models and allows for vectorized generation of batches of transcription traces. We assess computational speed and accuracy of common surrogate models along with new models, showing that high accuracy can be obtained while reducing computational cost up to two orders of magnitude.
Discovering reliable cause-and-effect relationships in real-world medical data is an open challenge. Classical Causal Discovery (CD) algorithms used to solve this task rely on strict assumptions that are rarely met in complex real-world scenarios with limited expert knowledge - the functional form of the causal relationships, the data distribution, the causal sufficiency. Thus, the reliability of CD algorithms can significantly drop, compromising the interpretability of the results and the trustworthiness of downstream decision-making. To overcome these limitations, we introduce the concept of consensus causal model to combine various CD algorithms and enhance their accuracy. Our consensus model can be efficiently constructed from a set of heterogeneous causal graph objects through a homogenisation step, ensuring semantic compatibility with the original edge definitions and enabling meaningful information exchange. To showcase the proposed method, we analyze a lung cancer dataset combining patient-level information such as smoking habits and age, and we study their effect on the onset and development of the disease, the tumor stage, and cellular pathway mutations. By applying multiple classical CD algorithms, we observe significant structural inconsistencies and heterogeneity across individual graphs. We demonstrate that the consensus causal model, unlike the individual models, effectively aggregates the strengths of each algorithm while mitigating their uncertainties. The resulting model reveals biologically validated causal relationships between risk factors, mutations, and pathways that isolated algorithms fail to capture, thereby underscoring the value of consensus causal modelling as a robust alternative to single-model selection for causal discovery.
Inferring culture media that enable specific metabolic functions is a challenging problem due to the vast combinatorial search space induced by all the possible subsets of compounds and reactions within a metabolic network. Existing scalable approaches based on flux optimization or combinatorial enumeration do not integrate prior biological knowledge to tailor the media to the cell physiological context. We introduce a method that solves culture media inference constrained by biological observation and knowledge as a combinatorial optimization problem, implemented using Answer Set Programming (ASP). This method combines four heuristics that progressively restrict the search space toward biologically supported solutions: (i) identify the target synthesis subnetwork achieving target compounds production, (ii) compute the parsimonious network merging the minimal reaction sets contextualized with the provided constraints, (iii) enumerate minimal media from this reduced search space, and (iv) detect and filter out the self activating internal compounds from the solution space, required for algorithmic reasons but not informative as environmental inputs. We evaluate this approach on different genome-scale metabolic networks and test varying biological constraints and heuristic combinations to assess the method’s scalability. This method provides a knowledge-driven and scalable enumeration of minimal culture media, combining logical reasoning, minimality optimizations, as well as knowledge and topology-based constraints to address the combinatorial complexity of the reverse ecology problem.
Biochemical systems involve both the flow of matter, in which entities transform into one another via reactions, and the flow of information, in which entities regulate which reactions may occur. Boolean networks capture the latter; reaction networks capture the former. Yet no unified qualitative formalism treats regulated reactions as its principal objects of study, despite their prominence in standards such as the Systems Biology Graphical Notation Process Description (SBGN-PD) language. We introduce modulation-reaction networks (MR-networks), a mathematical framework in which entities modulate reactions through activations and inhibitions, and study their synchronous Boolean semantics. To reason about MR-networks we develop Modulation-Reaction Logic (MRL), a hybrid modal μ-calculus whose modalities reason about the structure of the network and whose fixed-point operators capture temporal evolution of the computation. We establish a collection of validities, including a complete characterisation of the one-step update rule, and demonstrate the expressive power of MRL by formalising properties of biological interest such as reachability, sustained production, and presence of attractors. We show that MRL admits model-checking via an evaluation game, and introduce a bisimulation relation for MR-networks, which is proved to be invariant for all MRL-formulas. As a step towards a biologically more realistic computational model, we sketch the asynchronous semantics of MR-networks, and outline how the developments for the synchronous case transfer to the study of the asynchronous one.
The Chemical Reaction Networks (CRN) interpreted through the differential semantics, even when restricted to elementary reactions with mass action law kinetics, form a Turing-complete language. This means that any computable real function can thus be programmed, and in fact compiled, in an abstract CRN that will compute it with an arbitrarily high precision. In this computational framework, the information carriers are the molecular concentrations, the required precision is given as input, and the output concentration is guaranteed to satisfy the required precision. On the other hand, one can be interested in estimating the derivative of an unknown input signal or in reading the concentration value of an input molecular species. By nature, such problems can only be approximated with a finite precision. Hence, the computation framework proposed previously cannot be applied and we need to design and analyze custom CRNs to perform these tasks. In this paper, we present an analog-dyadic converter CRN which takes as input one molecular concentration (in [0, 1] but not necessarily computable), and produces as output a sequence of ”on” and ”off” spikes corresponding to some extent to the sequence of bits in the dyadic representation of the input concentration. We provide a detailed analysis of the source of errors and their behavior when varying the reactions rate constants. We conclude by sketching a possible design for a reader module that takes as input an arbitrary concentration and a desired precision and outputs a dyadic encoding approximating the value of the concentration with the desired precision. We leave as an open question to prove the correctness of our construction.
Biological regulatory networks can be represented by computational models, which allow the study and analysis of biological behaviours, therefore providing a better understanding of a given biological process. However, as new information is acquired, biological models may need to be revised in order to also account for this new information. Current model revision tools are scarce and often lack the flexibility to integrate with broader analysis workflows. Here, we present pyModRev, an enhanced iteration of the model revision tool ModRev, capable of verifying the consistency of Boolean regulatory models, and finding minimal repairs in case of inconsistency. pyModRev supports model validation against both steady state observations as well as time-series data, being able to consider different update schemes simultaneously. pyModRev supports different model formats, and is available as a Python package in PyPI, for easy integration with other model analysis tools, significantly improving accessibility and utility for the logical modelling community.
Qualitative models provide crucial instruments for modelling complex biological systems. While advances in automated reasoning and symbolic encodings have enabled rigorous inference of these models from data, the process remains highly fragile. First, biological measurement errors inevitably propagate into formal model specifications. Second, when a specification becomes unsatisfiable, distinguishing between fundamental design flaws and minor technical errors is notoriously difficult. This uncertainty often leads to under-specification, as it is unclear which observations are still “safe” to incorporate. To overcome these challenges, we introduce a robust inference method based on weighted MaxSMT. By encoding uncertain biological observations as weighted soft constraints, our approach enables the solver to identify a model best reflecting the observations, even with some conflicting constraints. Our method allows for Boolean and multi-valued variable domains, alongside observations derived from discretisation (level constraints) and differential expression (ordering constraints). We show our approach can be used to successfully infer neural cell differentiation models from prior-knowledge networks with approximately 200–1,300 genes using ordering constraints on all included genes.
Transcriptional networks represent one of the most extensively studied types of systems in synthetic biology. Although the completeness of transcriptional networks for digital logic is well-established, analog computation plays a crucial role in biological systems and offers significant potential for synthetic biology applications. While transcriptional circuits typically rely on cooperativity and highly nonlinear behavior of transcription factors to regulate protein production, they are often modeled with simple linear degradation terms. In contrast, general analog dynamics require both positive and negative nonlinear terms, seemingly necessitating control over not just transcriptional (i.e., production) regulation but also the degradation rates of transcription factors. Surprisingly, we prove that controlling transcription factor production (i.e., transcription rate) without explicitly controlling degradation is mathematically complete for analog computation, achieving equivalent capabilities to systems where both production and degradation are programmable. We demonstrate our approach on several examples including oscillatory and chaotic dynamics, analog sorting, memory, PID controller, and analog extremum seeking. Our result provides a systematic methodology for engineering novel analog dynamics using synthetic transcriptional networks without the added complexity of degradation control and informs our understanding of the capabilities of natural transcriptional circuits. We provide a compiler, in the form of a Python package that can take any system of polynomial ODEs and convert it to an equivalent transcriptional network implementing the system exactly, under appropriate conditions.
Stochastic dynamical systems like gene regulatory networks (GRNs) often exhibit behavior characterized by metastable sets (representing cellular phenotypes), in which trajectories remain for long times, whereas switches between these sets in the phase space are rare events. One way to capture these rare events is to infer the system's long-term behavior from the spectral characteristics (eigenvalues and eigenvectors) of its Koopman operator. For GRNs, the Koopman operator is based on the chemical master equation (CME), which provides a precise mathematical modeling framework for stochastic GRNs. Since the CME is typically analytically intractable, methods based on discretizing the CME operator have been developed. However, determining the number and location of metastable sets in the phase space as well as the transition rates between them remains computationally challenging, especially for large GRNs with many genes and interactions. A promising alternative method, called ISOKANN (invariant subspaces of Koopman operators with artificial neural networks) has been developed in the context of molecular dynamics. ISOKANN uses a combination of the power iteration and neural networks to learn the basis functions of an invariant subspace of the Koopman operator. In this paper, we extend the application of ISOKANN to the CME operator and apply it to two small GRNs: a genetic toggle switch model and a model for macrophage polarization. Our work opens a new field of application for the ISOKANN algorithm and demonstrates the potential of this algorithm for studying large GRNs.
The behaviour of microorganisms and microbial communities can be abstracted by models combining a description of their metabolic capabilities as metabolic networks, and suitable computational or mathematical paradigms that further integrate simulation conditions. A major component of the latter is the composition of the environment or growth medium that can be referred to as seeds. Predicting the seeds from the metabolic network and an expected behaviour is an inverse problem that can be addressed with linear programming or logic paradigms such as Answer Set Programming (ASP). Here, we formalise seed prediction for microbial communities, taking into account that their members may interact positively through metabolite transfers, which may reduce the need for external seed metabolites. We address the problem with ASP and add a hybrid component ensuring the satisfiability of linear constraints. We explore the subset-minimality solving heuristic of the Clingo solver and develop two heuristics supporting priority of seeds over transfers. We present a proof of concept of seed inference in small-scale communities, and assess the scalability of the three heuristics at genome-scale. Overall, our work introduces a hybrid logic-linear model for seed inference in interacting microbial communities, and new heuristics for the exploration of the solution space with subset minimality optimisations.
