IFAC-PapersOnLine

Parameter estimation in nonlinear biological dynamical systems is a difficult inverse problem because the governing equations are often stiff or oscillatory, the data are sparse and noisy, and the objective landscape is non-convex. Physics-informed neural networks (PINNs) offer an alternative to purely simulation-based calibration by representing state trajectories with neural networks while penalizing violations of the governing equations. This paper studies the empirical reliability of PINNs for recovering the parameters of the repressilator, a synthetic genetic oscillator formed by three cyclically repressive genes. We use synthetic time-series generated from the standard ordinary differential equation model and train inverse PINNs to estimate the production parameter {beta} and the Hill coefficient n. The study varies observation noise, partial observation of repressors, sampling density, sensitivity to initial parameter guesses, and the difference between stable and oscillatory regimes. The results show that PINNs can reconstruct trajectories accurately when the model structure is correct and the three repressors are observed, but parameter recovery is more fragile than trajectory fitting. Noise, sparse sampling, unobserved variables, and unfavorable initial guesses increase the risk of biased estimates. The stable regime is easier to reconstruct, whereas the oscillatory regime provides richer information but also exposes optimization sensitivity. These findings support PINNs as a useful reverse-engineering tool for small gene-regulatory ODE models, while highlighting the need for repeated runs, uncertainty reporting, and experimental designs that improve identifiability.

Spatio-temporal interactions shape microbial community dynamics. Metabolism, through competition and cross-feeding, is a foundational mechanism of these interactions. Flux balance analysis enables efficient simulation of steady-state metabolism. Integrating these simulations through time, using dynamic flux balance analysis, provides temporal predictions of growth and metabolism. Incorporating spatial context, through partial differential equations, enables spatio-temporal simulation of microbial communities. In this chapter, we step through this sequential process, moving from steady-state, to temporal, to spatio-temporal simulation of microbial community metabolism. We provide an illustrative example using the modeling software COMETS (Computation of Microbial Ecosystems in Time and Space) to simulate interacting bacterial colonies of Bifidobacterium longum subsp. infantis and Anaerobutyricum hallii (previously Eubacterium hallii). Within this simulation, both competition and cross-feeding influenced the production of butyrate leading to an intermediate optimal interaction distance for metabolite production. We outline each step and provide open-source code such that this simulation can serve as a template for future spatio-temporal simulations of microbial community metabolism.

Standard target-mediated drug disposition (TMDD) models are widely used to describe nonlinear pharmacokinetics driven by high-affinity drug-target interactions. However, their reliance on instantaneous binding limits their ability to capture delayed and history-dependent dynamics observed in vivo. Here, we introduce a fractional TMDD model that incorporates memory effects through a fractional derivative, thereby generalizing the standard TMDD (sTMDD) framework. Although this fractional TMDD (fTMDD) formulation increases modeling flexibility, it also exacerbates parameter identifiability challenges under typical experimental conditions where only drug concentration data are available. To address this limitation, we derive a fractional quasi-steady-state approximation (fQSSA) that reduces model dimensionality while preserving essential nonlinear and memory-dependent pharmacokinetic dynamics. We further establish an explicit validity condition that quantifies the approximation error of both fTMDD and fQSSA without requiring numerical simulation. This condition reveals that the initial drug-to-target ratio is the primary determinant of QSSA validity, whereas the fractional order has a comparatively minor influence. Application of the proposed framework to recombinant human erythropoietin (rhEPO) data demonstrates that fractional dynamics play a population-dependent role, improving model performance in adults but not in infants. Together, this work provides the first systematic derivation of a QSSA framework for fractional TMDD models, along with rigorous and computable applicability conditions. Our results establish a principled foundation for incorporating memory effects into pharmacokinetic modeling and offer a generalizable framework for nonlinear PK-PD systems involving binding-mediated dynamics. Author summaryMany drugs interact strongly with their biological targets, leading to complex and nonlinear pharmacokinetics that are commonly described using target-mediated drug disposition (TMDD) models. However, these models assume that drug-target interactions occur instantaneously, which limits their ability to capture delayed and history-dependent behaviors observed in real biological systems. In this study, we develop a new modeling framework that incorporates such memory effects by extending TMDD models using fractional calculus. To make the model more practical and computationally efficient, we derive a simplified version based on a quasi-steady-state approximation (QSSA) and provide a clear mathematical condition that determines when this simplification is valid. Our analysis shows that the accuracy of the simplified model is primarily controlled by the initial ratio of drug to target, while the influence of memory effects is comparatively smaller. When applied to experimental data for erythropoietin, our model reveals that memory effects are important in adults but negligible in infants, suggesting that these effects may reflect underlying physiological differences. Overall, this work provides a systematic and interpretable framework for incorporating memory effects into pharmacokinetic modeling, with potential applications to a wide range of drug systems involving complex binding dynamics.

The optimization and control of bioprocesses require robust in silico models that can accurately capture the complex and dynamic behavior of living cells. While hybrid models that combine machine learning with mechanistic equations have emerged as a powerful tools, they often require relatively large datasets and might yield inconsistent predictions that violate the stoichiometry of metabolism. In this study, we introduce FBA-Hyb, a multi-scale hybrid modeling framework that tightly integrates genome-scale metabolic networks via flux balance analysis (FBA) into its architecture. In our FBA-Hyb framework, artificial neural networks predict key FBA inputs (substrate uptake rates and cellular objectives) while a surrogate FBA module translates them into the metabolic fluxes that govern the bioprocess. A key novelty is that the FBA optimization step is replaced by a surrogate generated with symbolic regression, which encapsulates the FBA model into a compact analytical expression. This allows easy backpropagation through the integration of the neural controlled differential equationbased FBA-Hyb bioprocess model. We validated FBA-Hyb against a standard hybrid model (Std-Hyb) using two Escherichia coli fedbatch case studies. In the first study, FBA-Hyb achieved a 42 % average improvement in predictive accuracy (R2) during a leave-one-process-out cross validation. Crucially, FBA-Hyb maintains strict stoichiometric feasibility even during extrapolation. Meanwhile, an alternative approach based on standard architectures leads to stoichiometrically inconsistent solutions in 22 % of the cases analyzed. In the second case study, we demonstrate how FBA-Hyb effectively simulates unmeasured chemical species and discovers a metabolic shift in sulfate-limited regimes during bioprocessing. By providing a modular, biologically consistent, and computationally efficient architecture, FBA-Hyb offers a robust foundation for the next generation of bioprocess models and sustainable process optimization. Graphical Abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=81 SRC="FIGDIR/small/720062v1_ufig1.gif" ALT="Figure 1"> View larger version (28K): org.highwire.dtl.DTLVardef@16f011eorg.highwire.dtl.DTLVardef@b25b5borg.highwire.dtl.DTLVardef@18bd178org.highwire.dtl.DTLVardef@65274e_HPS_FORMAT_FIGEXP M_FIG C_FIG HighlightsO_LIFBA-Hyb integrates flux balance analysis (FBA) into hybrid bioprocess models. C_LIO_LISymbolic regression discovers a simple closed-form FBA surrogate model. C_LIO_LIThe FBA surrogate ensures accurate reaction stoichiometry. C_LIO_LIA neural network predicting the FBA objective keeps the model flexible. C_LIO_LIFBA-Hyb has superior capabilities and accuracy compared to the current standard. C_LI

Standard differentiation assays sample cell states only at discrete time points, while the underlying progression unfolds continuously and heterogeneously across cells. As a result, different combinations of proliferation, commitment, and maturation dynamics can converge to similar endpoint measurements. This many-to-one mapping between latent trajectories and observable readouts constitutes a partially observed inverse problem that limits mechanistic interpretation. Although this ambiguity is inherent to many experimental systems, it is rarely examined using models that connect cell-state dynamics to assay-level quantities. We present OsteoMin, a coarse-grained cellular automaton that links stochastic transitions between pre-osteoblast and osteoblast states to experimentally measurable readouts of alkaline phosphatase activity, collagen deposition, and mineralization. Model parameters were constrained using literature-reported kinetics and evaluated against dexamethasone and menaquinone-4 perturbations. The frame-work reproduces qualitative assay trends and enables systematic analysis of how cell-state progression, matrix maturation, and external perturbations shape differentiation outcomes. Using this framework, we quantify the identifiability limits of endpoint assays and test whether standard measurements can distinguish underlying differentiation mechanisms. Distinct perturbation families often produce similar endpoint responses (macro-F1 {approx} 0.42), indicating limited discriminative power. Incorporating temporal trajectories improves separability (macro-F1 {approx} 0.78), demonstrating that most identifiable information resides in marker dynamics rather than terminal measurements. Sobol analysis shows early markers depend on proliferation timing, whereas late mineralization is governed by nonlinear matrix maturation and parameter interactions. Together, these results show that endpoint assays constrain overall progression but do not uniquely identify underlying mechanisms. OsteoMin provides a framework linking differentiation dynamics to assay observables and a basis for assessing identifiability in endpoint-driven systems.

Human walking is intrinsically variable. For example, there is considerable stride to stride variability even when walking speed is constant. This variability is due to uncertainty in the sensorimotor system and the environment, and is shaped by both musculoskeletal dynamics (e.g. joint stiffness and damping originating from muscles) and the control strategy used to mitigate the effects of uncertainty. Yet, insight into how sensorimotor noise shapes walking variability is limited due to a lack of experimental methods to assess sensorimotor noise and control strategies during walking. Simulations that account for uncertainty can elucidate how sensorimotor noise affects movement variability but due to numerical challenges, accounting for sensorimotor noise is not common in simulations of walking. Existing simulations have hugely simplified musculoskeletal dynamics (e.g. no muscles), the control policy (e.g. pre-defined feedback loops), or sensorimotor noise sources (e.g. only motor noise). Here, we performed stochastic optimal control simulations of walking based on a model with 9 degrees of freedom and 18 muscles to study how the level of sensory and motor noise influences walking. We solved for feedforward muscle excitations and full-state time-varying feedback gains that minimised expected effort while generating periodic, and hence stable, gait patterns. To enable these simulations, we approximated the state distribution with a Gaussian and used an unscented transform to propagate the state covariance. Resulting optimisation problems were solved with direct collocation. Sensorimotor noise level had a small effect on the mean kinematics but shaped kinematic and muscle activity variability as well as expected effort. Although simulations underestimated the magnitude of experimental positional variability, they captured its structure. In agreement with experimental results, the control policy prioritised limiting variability of centre of mass kinematics and minimal swing foot clearance over limiting joint angle variability. Hence, our simulations suggest that effort minimisation underlies these observations. Author summaryWhen performing a movement multiple times, each repetition will be slightly different due to random disturbances in the neural signals used to control movement, i.e. sensorimotor noise. Because it is difficult to measure inside the nervous system of a moving person, computer simulations are used to study movement control. They found that both sensorimotor noise and musculoskeletal mechanics determine how people control arm movements and standing. However, there are no simulations of walking that systematically evaluated how sensorimotor noise level influences walking kinematics because they pose computational challenges. Here, we proposed and used an approach for minimal effort simulations of walking in the presence of uncertainty. We imposed forward speed and stability but not kinematics. We found that the level of sensorimotor noise had little effect on the mean movement but a strong effect on the variability and the expected effort. The control strategy prioritised reducing the variability of the centre of mass position and swing foot clearance over reducing the variability of individual joint angles, which is also observed in experiments. Interestingly, strict control of centre of mass position and foot clearance in our simulations emerged from minimising effort.

We present the Toroidal Search Algorithm (TSA), a novel population-based metaheuristic optimization method inspired by the topology of a torus. Conventional metaheuristics frequently suffer from boundary stagnation, a phenomenon that severely degrades performance in bounded and high-dimensional search spaces. TSA addresses this limitation by embedding the search domain into a toroidal geometry, thereby eliminating artificial boundaries and enabling continuous cyclic exploration. Beyond boundary handling, TSA uses winding numbers to capture the history of agent movement across periodic dimensions, which are exploited to adaptively refine local search. A modified sigmoid control function regulates the transition between global and local search. Performance of TSA is evaluated on a collection of unimodal and multimodal benchmark functions at various dimensions. It consistently outperforms established metaheuristics. Notably, TSA demonstrates exceptional robustness to increasing dimensionality, maintaining fast convergence and low variance where competing methods deteriorate. To assess real-world applicability, we apply TSA to an inverse problem from mathematical oncology. With both synthetic and clinical data, TSA reliably recovers physiologically plausible parameters with greater stability and predictive accuracy than competing algorithms. These results demonstrate that TSA is a powerful and robust tool for large-scale global optimization in computational modelling applications. Striking Image O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=200 SRC="FIGDIR/small/709766v1_ufig1.gif" ALT="Figure 1"> View larger version (131K): org.highwire.dtl.DTLVardef@1b2c994org.highwire.dtl.DTLVardef@d045a8org.highwire.dtl.DTLVardef@18d296corg.highwire.dtl.DTLVardef@9a972d_HPS_FORMAT_FIGEXP M_FIG C_FIG Image generated with Google Gemini.

The influence of arbitrary randomness in cell division times on the variability of protein copy numbers within a lineage ensemble has been recently studied, going beyond the contributions of noisy gene expression and partitioning error. However, variability of protein concentrations need separate study, since cell size growth between cell divisions dilute protein concentrations at the same rate as size growth, which also determines mean division times. Here for a model of bursty protein production, we present exact moments (of all orders) of protein concentrations in the cyclo-stationary state, comparing: (i) population and lineage cell ensembles, and (ii) statistics at different cell ages. Two interesting results emerge. While the variance of protein concentration changes with the degree of division time heterogeneity at any cell age, the age-averaged variance is independent of it within lineage ensemble but stays dependent within population ensemble. The skewness within population ensemble is higher in younger cells than within lineage ensemble, and this behavior reverses at older ages. Such a feature vanishes for the age-averaged distribution, with population based skewness always dominating over that of lineage. We also show that mother-daughter correlations in generation times, do not add any significant difference to the results.

The preferred stride ratio (PSR), defined as the ratio of step length to cadence, is approximately invariant across a wide range of walking speeds in healthy adults but breaks down at slow speeds. The lower speed boundary at which this constancy is broken was estimated by Murakami and Otaka (2017) to be approximately 62 m min-1 ({approx} 1.03 m s-1) on the basis of unstandardised K-means cluster analysis applied to data from 21 healthy adults at five speed conditions. The present report re-examines this estimate using the digitised individual-level scatter of Fig. 1-A and the published group-level statistics of Table 1 of that study, applying three breakpoint estimators in parallel: (i) unstandardised K-means (replicating the original method), (ii) a Gaussian mean-and-variance changepoint estimator, and (iii) a piecewise-linear regression on PSR. Applied directly to the digitised scatter (n = 84 resolved markers from a total of 105; 44 of 44 slow-walk markers, 40 of 61 normal-walk markers), the unstandardised K-means estimator returned 62.0 m min-1, matching the originally reported value to the reported precision; the mean-and-variance changepoint estimator returned 55 m min-1; and the piecewise-linear estimator was numerically unstable on the raw heteroscedastic data. To quantify uncertainty, 5 000 Monte Carlo realisations of synthetic individual-level data were generated from a bivariate truncated-normal model conditioned on the published condition means and standard deviations and on the published within-cluster speed-PSR correlations. The Monte Carlo distributions gave median estimates of 61 m min-1 (95 % MC interval 55-67) for unstandardised K-means, 39 m min-1 (29-53) for the mean-and-variance changepoint estimator, and 35 m min-1 (19-49) for piecewise-linear regression. Under a log-normal sensitivity model the corresponding intervals were 60 [55, 66], 34 [20, 58], and 19 [5, 42] m min-1. The likelihood-based estimator placed the central tendency substantially below 62 m min-1, and its Monte Carlo intervals did not include the original boundary under either marginal-distribution model. An additional robust heteroscedastic segmented profile-likelihood analysis on log-PSR yielded lower Monte Carlo median breakpoints across all model specifications, although the full-variance intervals overlapped the original K-means boundary. The qualitative finding of Murakami and Otaka -- that PSR constancy breaks down at slow walking speeds -- is supported by the present reanalysis. The original 62 m min-1 boundary is reproduced under the unstandardised K-means estimator, where it reflects the location of the largest density gap in the published five-condition speed sampling rather than a formally estimated changepoint; estimators formally designed for changepoint detection localise the joint PSR mean-and-variance transition substantially below this value. O_FIG O_LINKSMALLFIG WIDTH=162 HEIGHT=200 SRC="FIGDIR/small/720900v2_fig1.gif" ALT="Figure 1"> View larger version (41K): org.highwire.dtl.DTLVardef@2bb53dorg.highwire.dtl.DTLVardef@187d9bborg.highwire.dtl.DTLVardef@1e7a6a0org.highwire.dtl.DTLVardef@16c587b_HPS_FORMAT_FIGEXP M_FIG O_FLOATNOFigure 1.C_FLOATNO Reproduction and likelihood-based extension of the boundary reported in Murakami and Otaka [5]. (A) Digitised individual-level scatter from Fig. 1-A of [5] (n = 84 resolved markers from a total of 105: 44 of 44 slow-walk markers and 40 of 61 normal-walk markers). The dashed vertical line marks the value 62.4 m min-1 as drawn in the original figure. (B) PSR variance amplification across the five speed conditions, expressed as Var(PSR)/Var(PSR)Preferred, on a logarithmic vertical axis. (C) Distributions of the breakpoint estimates over N = 5 000 Monte Carlo realisations under the bivariate truncated-normal model with cluster-specific within-cluster correlations: unstandardised K-means (median 61 m min-1), the Gaussian mean-and-variance changepoint estimator (median 39 m min-1), and piecewise-linear regression on PSR (median 35 m min-1). The dashed vertical line marks 62.4 m min-1. (D) Sensitivity of each estimator to the choice of marginal-distribution model (truncated normal vs. log-normal); error bars are 95 % Monte Carlo simulation intervals. (E) PSR mean {+/-} SD across the five speed conditions (Table 1 of [5], height-adjusted). C_FIG O_TBL View this table: org.highwire.dtl.DTLVardef@24fe39org.highwire.dtl.DTLVardef@ae8fdborg.highwire.dtl.DTLVardef@66a473org.highwire.dtl.DTLVardef@b6ad84org.highwire.dtl.DTLVardef@139bca7_HPS_FORMAT_FIGEXP M_TBL O_FLOATNOTable 1.C_FLOATNO O_TABLECAPTIONSource data reproduced from Murakami and Otaka [5], height-adjusted, n = 21 per condition. C_TABLECAPTION C_TBL

As the United Kingdom (UK) targets net-zero emissions by 2050, anaerobic digestion (AD) has become a cornerstone of renewable energy infrastructure. However, mathematical models, such as the Anaerobic Digestion Model No. 1 (ADM1), often struggle with high-solids agricultural feedstocks because they rely on Chemical Oxygen Demand (COD), a metric that introduces significant experimental error. To overcome this, this study applies an established mass-based ADM1 framework tailored for the co-digestion of maize silage and cow manure sourced from a UK AD site. This study uses a parallel reactor framework, using two identical laboratory-scale reactors to physically replicate the dynamic conditions of the full-scale site. A Global Sensitivity Analysis was first conducted, identifying biomass decay and carbohydrate breakdown rates as the most influential factors affecting system stability and model accuracy. The model was calibrated using data from the first reactor and then tested against an independent second reactor subjected to significant organic loading stress. Results show high predictive capabilities, with the model achieving a R2 of 0.81 for biogas production during calibration. The model maintained high predictive accuracy during the validation test of the second physical twin, achieving an R2 of 0.85, proving that the framework is robust and not overfitted to a single dataset. While predicting rapid fluctuations in pH and alkalinity remains challenging, the mass-based approach effectively forecasts gas yields and process stability. This methodology provides a reliable foundation for robust process modelling, offering a scalable tool for the UK biogas sector to optimise AD. Graphical Abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=93 SRC="FIGDIR/small/721061v1_ufig1.gif" ALT="Figure 1"> View larger version (32K): org.highwire.dtl.DTLVardef@92c7e2org.highwire.dtl.DTLVardef@80d723org.highwire.dtl.DTLVardef@ac3d24org.highwire.dtl.DTLVardef@1e21a51_HPS_FORMAT_FIGEXP M_FIG C_FIG

Calibration of closed-loop lumped-parameter cardiovascular models remains a major bottleneck for scalable digital-twin generation because inverse estimation is ill-conditioned and typically requires computationally expensive iterative forward simulation. This study investigates whether a supervised neural network (NN) can provide a fast inverse estimator for a paediatric sepsis cardiovascular ODE model by learning a direct mapping from prescribed haemodynamic target vectors to calibrated parameter sets. Training data are generated by sampling model parameters at random, forward-simulating the closed-loop system to steady state, and pairing the resulting target summaries with the corresponding parameters; the same target definitions and evaluation populations are used throughout for consistency. We evaluate NN inference by forward re-simulation to steady state and benchmark performance against a simulator-constrained calibration reference (Embedded Gradient Descent, EGD) using relative-error statistics, distributional similarity of achieved outputs and inferred parameters (median shift, IQR ratio, Wasserstein distance, KS statistic), and target-space localisation of parameter-space disparity (cosine distance). The NN reproduces the prescribed targets with predominantly small errors for most samples, while the largest discrepancies are confined to a well defined set of target configurations that also yield high residuals under the reference method, indicating feasibility limits of the target/model combination. Overall, NN-guided calibration provides a computationally efficient accelerator for virtual-twin generation and target-space screening, with simulator-based refinement and forward re-simulation retained to handle infeasible regimes and enforce mechanistic plausibility.

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 Gillespies 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.

Successful development of multicellular organisms requires cells to transition between states with precise timing. Distinct cell states are often understood as being maintained by stabilizing regulatory networks, such that a complete cell-state transition requires network rewiring through partial dismantling of the current state and concurrent reconfiguration into a new one. Empirically, these transitions are often investigated by quantifying the gain or loss of expression of a small number of state-specific markers, frequently a single proxy. A general quantitative framework for comparing the kinetics of such transitions across experimental conditions is lacking. Here, we show that the delayed Weibull distribution provides a natural description of cell-state transition kinetics when transition is viewed as the cumulative consequence of many molecular events, whose timing may vary between cells and conditions, analogous to system failure in reliability theory. This formulation yields a compact model with interpretable parameters describing the delay before transition onset, the characteristic timescale of transition, the temporal form of the transition hazard, and the fraction of cells competent to respond. Together, this framework provides a practical and interpretable approach for quantifying the kinetics of cell-state transitions and how they are altered by perturbation.

The challenge of stratifying patients for combination therapy is both technically demanding and clinically crucial. In previous work, we introduced a multi-objective optimization framework for identifying optimally synergistic combination protocols that are robust to competing definitions of additivity. This manuscript extends this methodology to quantify how inter-individual variability in drug sensitivity influences the combination doses that optimally balance the competing objectives of synergy of efficacy and synergy of potency (a proxy measure of toxicity). For this methodology, we introduce a voxel-based stratification approach to characterize individuals (model parameterizations) into subgroups based on sensitivity to each drug as a monotherapy and in combination. As a case study, we apply the method to a preclinical dataset of murine response to the combination of an immune checkpoint inhibitor and an antiangiogenic agent. We demonstrate that the algorithm can quantify how the robustly optimal combination therapies vary across different treatment response subgroups and how the algorithm can identify subpopulations for which no meaningfully efficacious combination exists. As applying the methodology requires knowledge of specific parameter values for which measurable biomarkers may be unavailable, we also propose an initiation protocol that permits identification of the parameters necessary to place an individual in a subgroup. This methodology is a step in the direction of determining the right combination therapy for a subgroup and finding the right subgroup for an existing therapy.

SummaryRecent advances in analysis of biochemical networks have contributed the identification of their modular structure based on the concept of multi reaction dependencies and kinetic coupling of reaction rates (Kuken et al., 2022; Langary et al., 2025). Existing implementations of the algorithms to study modular structure do not scale well with the size of the networks, prohibiting their application with genome-scale networks. Here, we introduce COCOA.jl, a multithreaded Julia package for identification of concordant and kinetic modules, with applications in the study of concentration robustness. Availability and implementationCOCOA.jl is implemented in Julia 1.12.2 and is freely available under the MIT license at https://github.com/antoniofranky/COCOA.jl. It runs on Linux, macOS, and Windows; installation is supported via the Julia package manager. COCOA.jl can be called from Python via JuliaCall. Contactantonschaf@posteo.de; ankueken@uni-potsdam.de

Biological data is prone to both intrinsic and extrinsic noise and variability between experimental replicas. That same stochasticity and heterogeneity can carry information about underlying biochemical mechanisms but, if not incorporated in modeling and probabilistic inference, can also bias parameter estimates and misguide predictions and, subsequently, experiment design. Mechanistic inference typically requires lengthy simulations (e.g., the Stochastic Simulation Algorithm (SSA)); approximations to chemical master equation (CME) solutions that lack rigorous error tracking; or deterministic averaging that lacks the complexity necessary to reflect the data. We introduce the Stochastic System Identification Toolkit (SSIT) - a fast, flexible, and open-source software package available on GitHub that makes use of MATLABs efficient and diverse computational architecture. The SSIT is designed for building, simulating, and solving chemical reaction models using ODEs, moments, SSA, Finite State Projection truncations of the CME, or hybrid methods; sensitivity analysis and Fisher information quantification; parameter fitting using likelihood- or Bayesian-based methods; handling of experimental noise and measurement errors using probabilistic distortion operators; and sequential experiment design that empowers users to save time and resources while gaining the most information possible out of their data. The SSIT also offers advanced modeling tools, including model reduction methods for increased efficiency and joint fitting of models and datasets with overlapping reactions or parameters. To facilitate the ease and speed of use, the SSIT provides a graphical user interface and ready-made, adaptable pipelines that can be run in the background from commandline or high-performance computing clusters. We demonstrate features of the SSIT on two experimental datasets: the first consists of published mRNA count data that reflect Saccharomyces cerevisiae yeast cell response to osmotic shock using single-cell single-molecule fluorescence in situ hybridization; the second consists of single-cell RNA sequencing measurements of 151 activating genes in breast cancer cells following treatment with dexamethasone. Author summaryWe present the Stochastic System Identification Toolkit (SSIT) to model, fit, and predict any data that can be interpreted as changing populations or counts through time, including but not limited to single-cell experiments, economics, epidemiology, ecology, sociology, agriculture, and biotechnology. The SSIT was constructed particularly for stochastic modeling, which is important for systems whose states may experience significant fluctuations from mean behavior, thus affecting the inference of the underlying rate parameters and predictions of subsequent behavior. The SSIT provides statistical inference tools for parameter estimation; sensitivity analysis and information calculation; handling of distortions to probability distributions caused by experimental or measurement processes (e.g., dropout in single-cell RNA sequence data and total fluorescence intensities versus spot counting/puncta analysis); and quantitative design of experiments. The SSIT also offers a variety of complex modeling tools, including model reduction methods and fitting of combined models/datasets that share some behavior but remain distinct (e.g., different genes responding a single stimulus). The SSIT generates pipelines for easy, efficient analyses to run in the MATLAB environment, in the background on commandline, or on high-performance computing clusters, thus facilitating users to make informed, time- and cost-effective decisions about their next set of experiments.

Oscillatory behaviour is important in multiple biological contexts. However, the inherent nonlinearity and high dimensionality of mathematical models in biology makes proving the existence and the localization of limit cycle oscillations challenging. Here, we provided an elementary proof for the existence and a method for rigorously localizing the oscillatory solutions in a class of benchmark biomolecular oscillators. To construct the proof, we used a geometric approach based on Brouwers Fixed Point theorem. We constructed a toroidal-like manifold within a positively invariant set by removing the hypervolume containing the fixed point and the trajectories converging to it along its stable manifold. We showed that the vector field describing the system dynamics maps a cross section of the toroidal-like manifold onto itself. The existence of a limit cycle solution in this manifold was guaranteed by Brouwers Fixed Point theorem. For different sets of initial conditions in these cross-sections, we used an interval-based Reachability Analysis to localize the oscillatory behaviour that complements the Brouwers Fixed Point theorem approach. These results add a simple and elegant approach to demonstrating the existence of limit cycles in biomolecular systems as well as a method for rigorous localization of the region of existence.

The endometrium and menstrual disorders, such as endometriosis and adenomyosis, are difficult to study, partly because menstruation depends on interactions between multiple cell types through complex molecular mechanisms. To help understand this system, researchers need humanized experimental and computational models that can interrogate how endometrial cell populations impact each other in both physiological and pathological conditions. Here, we use ordinary differential equations (ODEs) to model changes in the rates of endometrial cell proliferation and death in response to hormones, cytokines, and the specific cell types present. To calibrate this computational model, we used previous-published experimental datasets from a 3D co-culture platform supporting primary human endometrial epithelial organoids and endometrial stromal cells. Our ODE-based model can simulate the size of endometrial epithelial organoids and the density of stromal cells over time under multiple hormone/cytokine conditions in mono- and co-cultures. We further created a second, partial differential equation (PDE)-based model that simulates the diffusion of molecules added to these 3D cultures and their uptake by proliferating endometrial cells using the predicted cell densities from the ODE model as inputs to the PDE simulations. We show that the exposure to hormones and cytokines used in the experiments is reasonably homogenous throughout the 3D culture and identify conditions where this would not be true. Altogether we use these models to quantify the influence of stromal cells on epithelial cell proliferation and vice versa, to identify differences across cells from different donors, and to provide a quantitative assessment of experimental designs.

Thin structures such as heart valves and aortic dissection flaps interact dynamically with blood flow in human vessels. Their flexibility and capacity for large deformations generate complex, highly transient hemodynamic patterns over the cardiac cycle. Accurately resolving these interactions remains challenging for conventional boundary-fitted fluid-structure interaction approaches. We present an immersed boundary method for simulating thin structures in incompressible flow on unstructured grids. The method couples a stabilized finite element fluid solver with a nonlinear, rotation-free shell formulation through a direct forcing immersed boundary approach. The framework supports both weak (explicit) and strong (implicit) time-coupling strategies, enabling stable simulations over a wide range of solid-to-fluid density ratios. Hydrodynamic forces acting on thin structures are computed from fluid solutions sampled on both sides of the structure, allowing accurate force reconstruction for zero-thickness shells. To our knowledge, this is the first immersed boundary formulation that couples an unstructured finite element fluid solver with a two-dimensional, rotation-free shell model to simulate interactions between thin structures and incompressible flow. Fluid-structure coupling is achieved using predefined finite element shape functions, which provide consistent projection between Eulerian and Lagrangian fields without additional interpolation procedures. The framework is validated using three-dimensional benchmark problems involving thin structures. Then, valve-like model is used to compare strong and weak coupling strategies. Finally, the method is applied to an idealized type-B aortic dissection model. The proposed approach is implemented within the open-source software CRIMSON, a finite element platform for cardiovascular simulation.

Variability in gene expression among single cells and growing cell populations can arise from the stochastic nature of protein synthesis, which often occurs in random bursts. This study investigates the variability in the expression of a growth-sustaining protein, whose concentration is regulated by a negative feedback loop due to cell growth-induced dilution. We model the distribution of protein concentration using a Chapman-Kolmogorov equation for single cells and a population balance equation for growing cell populations. For single cells, we derive an explicit solution for the protein concentration distribution in state space and represent it as a Bessel function in Laplace space. For growing populations, we find that the distribution satisfies a Heun differential equation with singular boundary conditions. By addressing the central connection problem for the Heun equation, we quantify the population-level protein distribution and determine the Mathusian parameter, which characterizes population growth. This work provides a comprehensive analytical framework for understanding how stochastic protein synthesis impacts gene expression variability and population dynamics.