pyPESTO: a modular and scalable tool for parameter estimation for dynamic models

A stochastic modelling framework for cancer patient trajectories: combining tumour growth, metastasis, and survival

Cancer is a major burden of disease around the globe and one of the leading causes of premature death. The key to improve patient outcomes in modern clinical cancer research is to gain insights into dynamics underlying cancer evolution in order to facilitate the search for effective therapies. However, most cancer data analysis tools are designed for controlled trials and cannot leverage routine clinical data, which are available in far greater quantities. In addition, many cancer models focus on single disease processes in isolation, disregarding interaction. This work proposes a unified stochastic modelling framework for cancer progression that combines (stochastic) processes for tumour growth, metastatic seeding, and patient survival to provide a comprehensive understanding of cancer progression. In addition, our models aim to use non-equidistantly sampled data collected in clinical routine to analyse the whole patient trajectory over the course of the disease. The model formulation features closed-form expressions of the likelihood functions for parameter inference from clinical data. The efficacy of our model approach is demonstrated through a simulation study involving four exemplary models, which utilise both analytic and numerical likelihoods. The results of the simulation studies demonstrate the accuracy and computational efficiency of the analytic likelihood formulations. We found that estimation can retrieve the correct model parameters and reveal the underlying data dynamics, and that this modelling framework is flexible in choosing the precise parameterisation. This work can serve as a foundation for the development of combined stochastic models for guiding personalized therapies in oncology.

Integrated mathematical and experimental modeling uncovers enhanced EMT plasticity upon loss of the DLC1 tumor suppressor

Epithelial-mesenchymal transition (EMT) plays an essential role in embryonic development, wound healing, and tumor progression. Partial EMT states have been linked to metastatic dissemination and drug resistance. Several interconnected feedback loops at the RNA and protein levels control the transition between different cellular states. Using a combination of mathematical modeling and experimental analyses in the TGFβ-responsive breast epithelial MCF10A cell model, we identify a central role for the tumor suppressor protein Deleted in Liver Cancer 1 (DLC1) during EMT. By extending a previous model of EMT comprising key transcription factors and microRNAs, our work shows that DLC1 acts as a positive regulator of TGFβ-driven EMT, mainly by promoting SNAIL1 expression. Our model predictions indicate that DLC1 loss impairs EMT progression. Experimental analyses confirm this prediction and reveal the acquisition of a partial EMT phenotype in DLC1-depleted cells. Furthermore, our model results indicate a possible EMT reversion to partial or epithelial states upon DLC1 loss in MCF10A cells induced toward mesenchymal phenotypes. The increased EMT plasticity of cells lacking DLC1 may explain its importance as a tumor suppressor.

Modeling diurnal rumen metabolism dynamics in dairy cattle: An update to a mechanistic model representing eating behavior, rumen content, rumination, and acid-base balance

Eating and rumination activities influence rumen environmental characteristics, such as volume, digesta particle size, microbial activity, and the acidity of rumen fluid. The aim of this study was to update an existing mechanistic, steady-state rumen model to simulate diurnal patterns of rumen functioning, contributing to a better understanding of rumen metabolism. Daily patterns of eating behavior were defined as model inputs based on meal characteristics and feeding regimens. Mechanisms controlling the dynamics of rumen content and fluid volume, particle diminution of roughage, rumination behavior, lactic acid metabolism, and rumen fluid pH were mathematically represented. The integration of these mechanisms into the existing steady-state rumen model enabled simulation of the relationship between eating behavior and rumen function. Scenarios of frequent eating (i.e., 12 equal meals) resulted in small fluctuations of rumen conditions during the day, whereas scenarios of restricted feeding with large meals (i.e., 4 daily meals, with 1 meal larger than the others) resulted in pronounced diurnal fluctuations in simulated rumen volume, lactic acid concentration, and rumen pH but not in simulated overall rumination activity. Also, the simulated responses of rumen function to changes in eating behavior showed interactions with dietary characteristics. Short, intense meals moderately increased daily mean rumen volume in high-roughage diets but not in high-concentrate diets. In contrast, the simulation of short and intense meals in concentrate-rich diets induced lactic acid formation and reduced the minimum pH, although having little effect in high-roughage diets. The model was calibrated and evaluated using a database including 78 experimental treatments from 25 studies where TMR diets were fed to lactating dairy cows. The prediction performance was assessed for rumen volume, VFA concentration, mean and minimum pH, NDF digestibility, and daily rumination time both for the original model and the presented updated model. Incorporating diurnal dynamics in the updated model allowed for the explanation of variability of the experimental data, and further independent evaluation is required to confirm its adequacy in predictive performance. In summary, the updated model integrates multiple factors involved in the diurnal patterns of rumen functioning, advancing quantitative understanding of this complex system.

The Dawn of High-Throughput and Genome-Scale Kinetic Modeling: Recent Advances and Future Directions

Researchers have invested much effort into developing kinetic models due to their ability to capture dynamic behaviors, transient states, and regulatory mechanisms of metabolism, providing a detailed and realistic representation of cellular processes. Historically, the requirements for detailed parametrization and significant computational resources created barriers to their development and adoption for high-throughput studies. However, recent advancements, including the integration of machine learning with mechanistic metabolic models, the development of novel kinetic parameter databases, and the use of tailor-made parametrization strategies, are reshaping the field of kinetic modeling. In this Review, we discuss these developments and offer future directions, highlighting the potential of these advances to drive progress in systems and synthetic biology, metabolic engineering, and medical research at an unprecedented scale and pace.

Assessment of uncertainty quantification in universal differential equations

Scientific machine learning is a new class of approaches that integrate physical knowledge and mechanistic models with data-driven techniques to uncover the governing equations of complex processes. Among the available approaches, universal differential equations (UDEs) combine prior knowledge in the form of mechanistic formulations with universal function approximators, such as neural networks. Integral to the efficacy of UDEs is the joint estimation of parameters for both the mechanistic formulations and the universal function approximators using empirical data. However, the robustness and applicability of these resultant models hinge upon the rigorous quantification of uncertainties associated with their parameters and predictive capabilities. In this work, we provide a formalization of uncertainty quantification (UQ) for UDEs and investigate key frequentist and Bayesian methods. By analyzing three synthetic examples of varying complexity, we evaluate the validity and efficiency of ensembles, variational inference and Markov-chain Monte Carlo sampling as epistemic UQ methods for UDEs.This article is part of the theme issue 'Uncertainty quantification for healthcare and biological systems (Part 2)'.

Quantifying cell cycle regulation by tissue crowding

The spatiotemporal coordination and regulation of cell proliferation is fundamental in many aspects of development and tissue maintenance. Cells have the ability to adapt their division rates in response to mechanical constraints, yet we do not fully understand how cell proliferation regulation impacts cell migration phenomena. Here, we present a minimal continuum model of cell migration with cell cycle dynamics, which includes density-dependent effects and hence can account for cell proliferation regulation. By combining minimal mathematical modeling, Bayesian inference, and recent experimental data, we quantify the impact of tissue crowding across different cell cycle stages in epithelial tissue expansion experiments. Our model suggests that cells sense local density and adapt cell cycle progression in response, during G1 and the combined S/G2/M phases, providing an explicit relationship between each cell-cycle-stage duration and local tissue density, which is consistent with several experimental observations. Finally, we compare our mathematical model's predictions to different experiments studying cell cycle regulation and present a quantitative analysis on the impact of density-dependent regulation on cell migration patterns. Our work presents a systematic approach for investigating and analyzing cell cycle data, providing mechanistic insights into how individual cells regulate proliferation, based on population-based experimental measurements.

CrossLabFit: A Novel Framework for Integrating Qualitative and Quantitative Data Across Multiple Labs for Model Calibration

The integration of computational models with experimental data is a cornerstone for gaining insight into biomedical applications. However, parameter fitting procedures often require a vast availability and frequency of data that are challenging to obtain from a single source. Here, we present a novel methodology "CrossLabFit" designed to integrate qualitative data from multiple laboratories, overcoming the constraints of single-lab data collection. Our approach harmonizes disparate qualitative assessments-ranging from different experimental labs to categorical observations-into a unified framework for parameter estimation. By using machine learning algorithms, these qualitative constraints are represented as dynamic "qualitative windows" that capture significant trends to which models must adhere. For numerical implementation, we developed a GPU-accelerated version of differential evolution to navigate in the cost function that integrated quantitative and qualitative data. We validate our approach across a series of case studies, demonstrating significant improvements in model accuracy and parameter identifiability. This work opens a new paradigm for collaborative science, enabling a methodological road to combine and compare findings between studies to improve our understanding of biological systems and beyond.

Mathematical basis and toolchain for hierarchical optimization of biochemical networks

Biological signalling systems are complex, and efforts to build mechanistic models must confront a huge parameter space, indirect and sparse data, and frequently encounter multiscale and multiphysics phenomena. We present HOSS, a framework for Hierarchical Optimization of Systems Simulations, to address such problems. HOSS operates by breaking down extensive systems models into individual pathway blocks organized in a nested hierarchy. At the first level, dependencies are solely on signalling inputs, and subsequent levels rely only on the preceding ones. We demonstrate that each independent pathway in every level can be efficiently optimized. Once optimized, its parameters are held constant while the pathway serves as input for succeeding levels. We develop an algorithmic approach to identify the necessary nested hierarchies for the application of HOSS in any given biochemical network. Furthermore, we devise two parallelizable variants that generate numerous model instances using stochastic scrambling of parameters during initial and intermediate stages of optimization. Our results indicate that these variants produce superior models and offer an estimate of solution degeneracy. Additionally, we showcase the effectiveness of the optimization methods for both abstracted, event-based simulations and ODE-based models.

Benchmarking methods for computing local sensitivities in ordinary differential equation models at dynamic and steady states

Estimating parameters of dynamic models from experimental data is a challenging, and often computationally-demanding task. It requires a large number of model simulations and objective function gradient computations, if gradient-based optimization is used. In many cases, steady-state computation is a part of model simulation, either due to steady-state data or an assumption that the system is at steady state at the initial time point. Various methods are available for steady-state and gradient computation. Yet, the most efficient pair of methods (one for steady states, one for gradients) for a particular model is often not clear. In order to facilitate the selection of methods, we explore six method pairs for computing the steady state and sensitivities at steady state using six real-world problems. The method pairs involve numerical integration or Newton's method to compute the steady-state, and-for both forward and adjoint sensitivity analysis-numerical integration or a tailored method to compute the sensitivities at steady-state. Our evaluation shows that all method pairs provide accurate steady-state and gradient values, and that the two method pairs that combine numerical integration for the steady-state with a tailored method for the sensitivities at steady-state were the most robust, and amongst the most computationally-efficient. We also observed that while Newton's method for steady-state computation yields a substantial speedup compared to numerical integration, it may lead to a large number of simulation failures. Overall, our study provides a concise overview across current methods for computing sensitivities at steady state. While our study shows that there is no universally-best method pair, it also provides guidance to modelers in choosing the right methods for a problem at hand.

A sphingolipid rheostat controls apoptosis versus apical cell extrusion as alternative tumour-suppressive mechanisms

Evasion of cell death is a hallmark of cancer, and consequently the induction of cell death is a common strategy in cancer treatment. However, the molecular mechanisms regulating different types of cell death are poorly understood. We have formerly shown that in the epidermis of hypomorphic zebrafish hai1a mutant embryos, pre-neoplastic transformations of keratinocytes caused by unrestrained activity of the type II transmembrane serine protease Matriptase-1 heal spontaneously. This healing is driven by Matriptase-dependent increased sphingosine kinase (SphK) activity and sphingosine-1-phosphate (S1P)-mediated keratinocyte loss via apical cell extrusion. In contrast, amorphic hai1afr26 mutants with even higher Matriptase-1 and SphK activity die within a few days. Here we show that this lethality is not due to epidermal carcinogenesis, but to aberrant tp53-independent apoptosis of keratinocytes caused by increased levels of pro-apoptotic C16 ceramides, sphingolipid counterparts to S1P within the sphingolipid rheostat, which severely compromises the epidermal barrier. Mathematical modelling of sphingolipid rheostat homeostasis, combined with in vivo manipulations of components of the rheostat or the ceramide de novo synthesis pathway, indicate that this unexpected overproduction of ceramides is caused by a negative feedback loop sensing ceramide levels and controlling ceramide replenishment via de novo synthesis. Therefore, despite their initial decrease due to increased conversion to S1P, ceramides eventually reach cell death-inducing levels, making transformed pre-neoplastic keratinocytes die even before they are extruded, thereby abrogating the normally barrier-preserving mode of apical live cell extrusion. Our results offer an in vivo perspective of the dynamics of sphingolipid homeostasis and its relevance for epithelial cell survival versus cell death, linking apical cell extrusion and apoptosis. Implications for human carcinomas and their treatments are discussed.

Computational Model Predicts Patient Outcomes in Luminal B Breast Cancer Treated with Endocrine Therapy and CDK4/6 Inhibition

PURPOSE

Development of a computational biomarker to predict, prior to treatment, the response to CDK4/6 inhibition (CDK4/6i) in combination with endocrine therapy in patients with breast cancer.

EXPERIMENTAL DESIGN

A mechanistic mathematical model that accounts for protein signaling and drug mechanisms of action was developed and trained on extensive, publicly available data from breast cancer cell lines. The model was built to provide a patient-specific response score based on the expression of six genes (CCND1, CCNE1, ESR1, RB1, MYC, and CDKN1A). The model was validated in five independent cohorts of 148 patients in total with early-stage or advanced breast cancer treated with endocrine therapy and CDK4/6i. Response was measured either by evaluating Ki67 levels and PAM50 risk of relapse (ROR) after neoadjuvant treatment or by evaluating progression-free survival (PFS).

RESULTS

The model showed significant association with patient's outcomes in all five cohorts. The model predicted high Ki67 [area under the curve; AUC (95% confidence interval, CI) of 0.80 (0.64-0.92), 0.81 (0.60-1.00) and 0.80 (0.65-0.93)] and high PAM50 ROR [AUC of 0.78 (0.64-0.89)]. This observation was not obtained in patients treated with chemotherapy. In the other cohorts, patient stratification based on the model prediction was significantly associated with PFS [hazard ratio (HR) = 2.92 (95% CI, 1.08-7.86), P = 0.034 and HR = 2.16 (1.02 4.55), P = 0.043].

CONCLUSIONS

A mathematical modeling approach accurately predicts patient outcome following CDK4/6i plus endocrine therapy that marks a step toward more personalized treatments in patients with Luminal B breast cancer.

Quantitative Online Monitoring of an Immobilized Enzymatic Network by Ion Mobility-Mass Spectrometry

The forward design of in vitro enzymatic reaction networks (ERNs) requires a detailed analysis of network kinetics and potentially hidden interactions between the substrates and enzymes. Although flow chemistry allows for a systematic exploration of how the networks adapt to continuously changing conditions, the analysis of the reaction products is often a bottleneck. Here, we report on the interface between a continuous stirred-tank reactor, in which an immobilized enzymatic network made of 12 enzymes is compartmentalized, and an ion mobility-mass spectrometer. Feeding uniformly 13C-labeled inputs to the enzymatic network generates all isotopically labeled reaction intermediates and products, which are individually detected by ion mobility-mass spectrometry (IMS-MS) based on their mass-to-charge ratios and inverse ion mobilities. The metabolic flux can be continuously and quantitatively monitored by diluting the ERN output with nonlabeled standards of known concentrations. The real-time quantitative data obtained by IMS-MS are then harnessed to train a model of network kinetics, which proves sufficiently predictive to control the ERN output after a single optimally designed experiment. The high resolution of the time-course data provided by this approach is an important stepping stone to design and control sizable and intricate ERNs.

Processivity and specificity of histone acetylation by the male-specific lethal complex

Acetylation of lysine 16 of histone H4 (H4K16ac) stands out among the histone modifications, because it decompacts the chromatin fiber. The metazoan acetyltransferase MOF (KAT8) regulates transcription through H4K16 acetylation. Antibody-based studies had yielded inconclusive results about the selectivity of MOF to acetylate the H4 N-terminus. We used targeted mass spectrometry to examine the activity of MOF in the male-specific lethal core (4-MSL) complex on nucleosome array substrates. This complex is part of the Dosage Compensation Complex (DCC) that activates X-chromosomal genes in male Drosophila. During short reaction times, MOF acetylated H4K16 efficiently and with excellent selectivity. Upon longer incubation, the enzyme progressively acetylated lysines 12, 8 and 5, leading to a mixture of oligo-acetylated H4. Mathematical modeling suggests that MOF recognizes and acetylates H4K16 with high selectivity, but remains substrate-bound and continues to acetylate more N-terminal H4 lysines in a processive manner. The 4-MSL complex lacks non-coding roX RNA, a critical component of the DCC. Remarkably, addition of RNA to the reaction non-specifically suppressed H4 oligo-acetylation in favor of specific H4K16 acetylation. Because RNA destabilizes the MSL-nucleosome interaction in vitro we speculate that RNA accelerates enzyme-substrate turn-over in vivo, thus limiting the processivity of MOF, thereby increasing specific H4K16 acetylation.

Bayesian modelling of time series data (BayModTS)-a FAIR workflow to process sparse and highly variable data

MOTIVATION

Systems biology aims to better understand living systems through mathematical modelling of experimental and clinical data. A pervasive challenge in quantitative dynamical modelling is the integration of time series measurements, which often have high variability and low sampling resolution. Approaches are required to utilize such information while consistently handling uncertainties.

RESULTS

We present BayModTS (Bayesian modelling of time series data), a new FAIR (findable, accessible, interoperable, and reusable) workflow for processing and analysing sparse and highly variable time series data. BayModTS consistently transfers uncertainties from data to model predictions, including process knowledge via parameterized models. Further, credible differences in the dynamics of different conditions can be identified by filtering noise. To demonstrate the power and versatility of BayModTS, we applied it to three hepatic datasets gathered from three different species and with different measurement techniques: (i) blood perfusion measurements by magnetic resonance imaging in rat livers after portal vein ligation, (ii) pharmacokinetic time series of different drugs in normal and steatotic mice, and (iii) CT-based volumetric assessment of human liver remnants after clinical liver resection.

AVAILABILITY AND IMPLEMENTATION

The BayModTS codebase is available on GitHub at https://github.com/Systems-Theory-in-Systems-Biology/BayModTS. The repository contains a Python script for the executable BayModTS workflow and a widely applicable SBML (systems biology markup language) model for retarded transient functions. In addition, all examples from the paper are included in the repository. Data and code of the application examples are stored on DaRUS: https://doi.org/10.18419/darus-3876. The raw MRI ROI voxel data were uploaded to DaRUS: https://doi.org/10.18419/darus-3878. The steatosis metabolite data are published on FairdomHub: 10.15490/fairdomhub.1.study.1070.1.

Long-term monitoring of SARS-CoV-2 seroprevalence and variants in Ethiopia provides prediction for immunity and cross-immunity

Under-reporting of COVID-19 and the limited information about circulating SARS-CoV-2 variants remain major challenges for many African countries. We analyzed SARS-CoV-2 infection dynamics in Addis Ababa and Jimma, Ethiopia, focusing on reinfection, immunity, and vaccination effects. We conducted an antibody serology study spanning August 2020 to July 2022 with five rounds of data collection across a population of 4723, sequenced PCR-test positive samples, used available test positivity rates, and constructed two mathematical models integrating this data. A multivariant model explores variant dynamics identifying wildtype, alpha, delta, and omicron BA.4/5 as key variants in the study population, and cross-immunity between variants, revealing risk reductions between 24% and 69%. An antibody-level model predicts slow decay leading to sustained high antibody levels. Retrospectively, increased early vaccination might have substantially reduced infections during the delta and omicron waves in the considered group of individuals, though further vaccination now seems less impactful.

Iterative design of training data to control intricate enzymatic reaction networks

Kinetic modeling of in vitro enzymatic reaction networks is vital to understand and control the complex behaviors emerging from the nonlinear interactions inside. However, modeling is severely hampered by the lack of training data. Here, we introduce a methodology that combines an active learning-like approach and flow chemistry to efficiently create optimized datasets for a highly interconnected enzymatic reactions network with multiple sub-pathways. The optimal experimental design (OED) algorithm designs a sequence of out-of-equilibrium perturbations to maximize the information about the reaction kinetics, yielding a descriptive model that allows control of the output of the network towards any cost function. We experimentally validate the model by forcing the network to produce different product ratios while maintaining a minimum level of overall conversion efficiency. Our workflow scales with the complexity of the system and enables the optimization of previously unobtainable network outputs.

Data integration strategies for whole-cell modeling

Data makes the world go round-and high quality data is a prerequisite for precise models, especially for whole-cell models (WCM). Data for WCM must be reusable, contain information about the exact experimental background, and should-in its entirety-cover all relevant processes in the cell. Here, we review basic requirements to data for WCM and strategies how to combine them. As a species-specific resource, we introduce the Yeast Cell Model Data Base (YCMDB) to illustrate requirements and solutions. We discuss recent standards for data as well as for computational models including the modeling process as data to be reported. We outline strategies for constructions of WCM despite their inherent complexity.