IFAC-PapersOnLine

Process models are increasingly used to support upstream process development in the biopharmaceutical industry for process optimization, scale-up and to reduce experimental effort. Parametric unstructured models based biological mechanisms are highly promising, since they do not require large amounts of data. The critical part in the application is the certainty of the parameter estimates, since uncertainty of the parameter estimates propagates to model predictions and can increase the risk associated with those predictions. Currently Fisher-Information-Matrix based approximations or Monte-Carlo approaches are used to estimate parameter confidence intervals and regularization approaches to decrease parameter uncertainty. Here we apply profile likelihood to determine parameter identifiability of a recent upstream process model. We have investigated the effect of data amount on identifiability and found out that addition of data reduces non-identifiability. The likelihood profiles of non-identifiable parameters were then used to uncover structural model changes. These changes effectively alleviate the remaining non-identifiabilities except for a single parameter out of 21 total parameters. We present the first application of profile likelihood to a complete upstream process model. Profile likelihood is a highly suitable method to determine parameter confidence intervals in upstream process models and provides reliable estimates even with non-linear models and limited data.

We present an automated model reduction algorithm that uses quasi-steady state approximation to minimize the error between the desired outputs. Additionally, the algorithm minimizes the sensitivity of the error with respect to parameters to ensure robust performance of the reduced model in the presence of parametric uncertainties. We develop the theory for this model reduction algorithm and present the implementation of the algorithm that can be used to perform model reduction of given SBML models. To demonstrate the utility of this algorithm, we consider the design of a synthetic biological circuit to control the population density and composition of a consortium consisting of two different cell strains. We show how the model reduction algorithm can be used to guide the design and analysis of this circuit.

This paper concerns the inverse problem of characterising the state of a bioreactor from observations. In laboratory settings, the bioreactor is represented by a device called a chemostat. We consider a differential description of the evolution of the state of the chemostat under environmental fluctuations. First, we model the state evolution as a stochastic process driven by Brownian motion. Under this model, our best knowledge about the state of the chemostat is described by its probability distribution in time, given the distribution of the initial state. The corresponding probability density function solves a deterministic partial differential equation (PDE), the Kolmogorov forward equation. While this provides a probabilistic description, incorporating an observation process allows for a more refined characterisation of the state. More formally, we are interested in obtaining the distribution of the state conditional on an observation process as the solution to a filtering problem, with the corresponding conditional probability density function solving a non-linear stochastic PDE, the Kushner-Stratonovich equation. This paper focuses on the pathwise formulation of this filtering problem in which inferences about the state are obtained conditional on a fixed stream of observations. We establish the existence and uniqueness of solutions to the governing differential equations, ensuring well-posedness before presenting numerical approximations. We approximate the pathwise solution to the filtering problem by combining the finite difference and splitting methods for solving PDEs, and then compare the approximated solution with results from a linearisation method and a classical sequential Monte Carlo method.

The paper addresses the critical challenge of accurately characterising steady states in biomolecular systems, which are often complex, nonlinear, multistable and subject to significant uncertainties. Traditional numerical methods often fail to provide complete or guaranteed solutions under these conditions. To overcome these limitations, the research proposes and evaluates the application of interval analysis methodologies. We provided algorithms for interval Newton and interval Krawczyk methods for rigorously bounding all possible steady states (both stable and unstable) in multistable, multidimensional nonlinear systems. This study involves a comparative analysis of these two methods in conjunction with interval bisection and interval constraint propagation. We addressed numerical examples for an array of biologically plausible models, involving both feedback and feedforward gene networks. The work recommends the choice of the most suitable method for various types of biomolecular systems, ultimately offering a robust computational framework to understand cellular functions and design synthetic biological circuits.

Modeling dynamical biological systems is key for understanding, predicting, and controlling complex biological behaviors. Traditional methods for identifying governing equations, such as ordinary differential equations (ODEs), typically require extensive quantitative data, which is often scarce in biological systems due to experimental limitations. To address this challenge, we introduce an approach that determines biomolecular models from qualitative system behaviors expressed as Signal Temporal Logic (STL) statements, which are naturally suited to translate expert knowledge into computationally tractable specifications. Our method represents the biological network as a graph, where edges represent interactions between species, and uses a genetic algorithm to identify the graph. To infer the parameters of the ODEs modeling the interactions, we propose a gradient-based algorithm. On a numerical example, we evaluate two loss functions using STL robustness and analyze different initialization techniques to improve the convergence of the approach.

Ordinary differential equations are frequently employed for mathematical modeling of biological systems. The identification of mechanisms that are specific to certain cell types is crucial for building useful models and to gain insights into the underlying biological processes. Regularization techniques have been proposed and applied to identify mechanisms specific to two cell types, e.g., healthy and cancer cells, including the LASSO (least absolute shrinkage and selection operator). However, when analyzing more than two cell types, these approaches are not consistent, and require the selection of a reference cell type, which can affect the results. To make the regularization approach applicable to identifying cell-type specific mechanisms in any number of cell types, we propose to incorporate the clustered LASSO into the framework of ordinary differential equation modeling by penalizing the pairwise differences of the logarithmized fold-change parameters encoding a specific mechanism in different cell types. The symmetry introduced by this approach renders the results independent of the reference cell type. We discuss the necessary adaptations of state-of-the-art numerical optimization techniques and the process of model selection for this method. We assess the performance with realistic biological models and synthetic data, and demonstrate that it outperforms existing approaches. Finally, we also exemplify its application to published biological models including experimental data, and link the results to independent biological measurements. Contactadrian.hauber@fdm.uni-freiburg.de Author SummaryMathematical models enable insights into biological systems beyond what is possible in the wet lab alone. However, constructing useful models can be challenging, since they both need a certain amount of complexity to adequately describe real-world observations, and simultaneously enough simplicity to enable understanding of these observations and precise predictions. Regularization techniques were suggested to tackle this challenge, especially when building models that describe two different types of cells, such as healthy and cancer cells. Typically, both cell types have a large portion of biological mechanisms in common, and the task is to identify the relevant differences that need to be included into the model. For more than two types of cells, the existing approaches are not readily applicable, because they require defining one of the cell types as reference, which potentially influences the results. In this work, we present a regularization method that is independent from the choice of a reference. We demonstrate its working principle and compare its performance to existing approaches. Since we implemented this method in a freely available software package, it is accessible to a broad range of researchers and will facilitate the construction of useful mathematical models for multiple types of cells.

Here we present details of the calculation necessary to estimate the impact of avidity in a mathematical model of a bispecific antibody with two membrane-bound targets. The calculation is used to reproduce the results reported in Rhoden et al. (2016) and implemented in Kareva et al. (2018). We reproduce the impact of difference in relative concentration of the two targets on projections of free and bound concentrations of both targets and the antibody and highlight the applicability of this approach for supporting model informed decision making particularly in the early stages of drug discovery and development.

The dynamics of heartbeat intervals provide important insights into cardiovascular and autonomic nervous system function. Conventional analytical approaches often use fixed-window averaging, which can obscure rapid changes and reduce temporal resolution. Point process models address this limitation by operating in continuous time, enabling more precise characterization of heartbeat variability. A landmark example is the history-dependent inverse Gaussian (IG) point process model of Barbieri et al. (2005), which captures temporal dependencies in heartbeat timing. However, the nonconvex likelihood of the IG model complicates parameter estimation, requiring careful initialization and adding computational burden. In this work, we introduce a convex alternative: a history-dependent gamma generalized linear model (GLM) for heartbeat dynamics. Applied to a tilt-table dataset, our approach yields accurate and robust heart rate estimation. We further extend the model to two more applications: (1) sequential prediction of interbeat intervals, outperforming common machine learning algorithms, and (2) computation of information-theoretic measures demonstrating its utility in quantifying the influence of cardiac medications on heartbeat dynamics.

The long-term behaviors of biochemical systems are described by their steady states. Deriving these states directly for complex networks arising from real-world applications, however, is often challenging. Recent work has consequently focused on network-based approaches. Specifically, biochemical reaction networks are transformed into weakly reversible and deficiency zero networks, which allows the derivation of their analytic steady states. Identifying this transformation, however, can be challenging for large and complex networks. In this paper, we address this difficulty by breaking the complex network into smaller independent subnetworks and then transforming the subnetworks to derive the analytic steady states of each subnetwork. We show that stitching these solutions together leads to the the analytic steady states of the original network. To facilitate this process, we develop a user-friendly and publicly available package, COMPILES (COMPutIng anaLytic stEady States). With COMPILES, we can easily test the presence of bistability of a CRISPRi toggle switch model, which was previously investigated via tremendous number of numerical simulations and within a limited range of parameters. Furthermore, COMPILES can be used to identify absolute concentration robustness (ACR), the property of a system that maintains the concentration of particular species at a steady state regardless of any initial concentrations. Specifically, our approach completely identifies all the species with and without ACR in a complex insulin model. Our method provides an effective approach to analyzing and understanding complex biochemical systems. Author summarySteady states describe the long-term behaviors of biochemical systems, which are typically based on ordinary differential equations. To derive a steady state analytically, significant attention has been given in recent years to network-based approaches. While this approach allows a steady state to be derived as long as a network has a special structure, complex and large networks rarely have this structural property. We address this difficulty by breaking the network into smaller and more manageable independent subnetworks, and then use the network-based approach to derive the analytic steady state of each subnetwork. Stitching these solutions together allows us to derive the analytic steady state of the original network. To facilitate this process, we develop a user-friendly and publicly available package, COMPILES. COMPILES identifies critical biochemical properties such as the presence of bistability in a genetic toggle switch model and absolute concentration robustness in a complex insulin signaling pathway model.

SummaryThe dynamics of ordinary differential equation (ODE) models strongly depend on the model structure, in particular the existence of positive and negative feedback loops. LoopDetect offers user-friendly detection of all feedback loops in ODE models in three programming languages frequently used to solve and analyze them: MATLAB, Python, and R. The developed toolset accounts for user-defined model parametrizations and states of the modelled variables and supports feedback loop detection over ranges of values. It generates output in an easily adaptable format for further investigation. Availability and ImplementationLoopDetect is implemented in R, Python 3 and MATLAB. It is freely available at https://cran.r-project.org/web/packages/LoopDetectR/, https://pypi.org/project/loopdetect/, https://de.mathworks.com/matlabcentral/fileexchange/81928-loopdetect/ (GPLv3 or BSD license). Contactkatharina.baum@hpi.de

MotivationMetabolic kinetic models are widely used to model biological systems. Despite their widespread use, it remains challenging to parameterize these Ordinary Differential Equations (ODE) for large scale kinetic models. Recent work on neural ODEs has shown the potential for modeling time-series data using neural networks, and many methodological developments in this field can similarly be applied to kinetic models. ResultsWe have implemented a simulation and training framework for Systems Biology Markup Language (SBML) models using JAX/Diffrax, which we named jaxkineticmodel. JAX allows for automatic differentiation and just-in-time compilation capabilities to speed up the parameterization of kinetic models. We show the robust capabilities of training kinetic models using this framework on a large collection of SBML models with different degrees of prior information on parameter initialization. Finally, we showcase the training framework implementation on a complex model of glycolysis. These results show that our framework can be used to fit large metabolic kinetic models efficiently and provides a strong platform for modeling biological systems. ImplementationImplementation of jaxkineticmodel is available as a Python package at https://github.com/AbeelLab/jaxkineticmodel. Author summaryUnderstanding how metabolism works from a systems perspective is important for many biotechnological applications. Metabolic kinetic models help in achieving understanding, but there construction and parametrization has proven to be complex, especially for larger metabolic networks. Recent success in the field of neural ordinary differential equations in combination with other mathematical/computational techniques may help in tackling this issue for training kinetic models. We have implemented a Python package named jaxkineticmodel that can be used to build, simulate and train kinetic models, as well as compatibility with the Systems Biology Markup Language. This framework allows for efficient training of kinetic models on time-series concentration data using a neural ordinary differential equation inspired approach. We show the convergence properties on a large collection of SBML models, as well as experimental data. This shows a robust training process for models with hundreds of parameters, indicating that it can be used for large-scale kinetic model training.

Time-lapse microscopy has become increasingly prevalent in biological experimentation, as it provides single-cell trajectories that unveil valuable insights into underlying networks and their stochastic dynamics. However, the limited availability of fluorescent reporters typically constrains tracking to only a few network species. Addressing this challenge, the dynamic estimation of hidden state-components becomes crucial, for which stochastic filtering presents a robust mathematical framework. Yet, the complexity of biological networks often renders direct solutions to the filtering equation intractable due to high dimensionality and nonlinear interactions. In this study, we establish and rigorously prove the well-posedness of the filtering equation for the time-evolution of the conditional distribution of hidden species. Focusing on continuous-time, noise-free observations within a continuous-time discrete state-space Markov chain model, we develop the Filtered Finite State Projection (FFSP) method. This computational approach offers an approximated solution by truncating the hidden species state space, accompanied by computable error bounds. We illustrate the effectiveness of FFSP through diverse numerical examples, comparing it with established filtering techniques such as the Kalman filter, Extended Kalman filter, and particle filter. Finally, we show an application of our methodology with real time-lapse microscopy data. This work not only advances the application of stochastic filtering to biological systems but also contributes towards more accurate implementation of biomolecular feedback controllers. Author SummaryThe aim of this paper is to introduce a novel computational approach for numerically solving high-dimensional filtering problems associated with stochastic reaction network models in intracellular processes. This method, termed the Filtered Finite State Projection (FFSP) method, can reliably predict the dynamics of hidden species in reaction systems based on time-course measurements of the stochastic trajectories of certain species. While stochastic filtering is extensively utilised in engineering, its application in biology has been limited, primarily due to the nonlinear nature of biological interactions and the discrete, non-Gaussian nature of state variables. Traditional filtering techniques, such as the Kalman filter, often encounter difficulties under these conditions. We demonstrate that the FFSP method provides an accurate solution to the stochastic filtering problem, complete with a computable error bound. We present several numerical examples to showcase the effectiveness of FFSP and its superior performance compared to other filtering methodologies. Additionally, we apply FFSP to biological data, successfully reconstructing the hidden dynamics of a yeast transcription system from partial measurements obtained through time-lapse microscopy. We believe that FFSP could be a valuable tool for elucidating hidden intracellular dynamics and understanding stochastic cellular behaviours.

Cells, in order to thrive, make efficient use of metabolites, proteins, energy, membrane space, and time. How, for example, should they allocate the available amount of protein to different metabolic pathways or cell functions? To model metabolic behaviour as an economic problem, some flux analysis model, kinetic models, and cell models apply optimality principles. However, due to their different assumptions these models are hard to compare and combine. Benefits and costs of metabolic pathways - e.g. favouring high production fluxes and low metabolite and enzyme cost - can be derived from general fitness objectives such as fast cell growth. To define pathway objectives, we may assume "optimistically" that, given a pathway state, any cell variables outside the pathway will be chosen for maximal fitness. The resulting fitness defines an effective pathway objective as a function of the pathway variables. Here I propose a unified theory that considers kinetic models, describes the set of feasible states as a state manifold and score each state by cost and benefit functions for fluxes, metabolite concentrations, and enzyme levels. To screen the state manifold and to find optimal states, the problem can be projected into flux, metabolite, or enzyme space, where effective cost and benefit functions are used. We reobtain existing modelling approaches such as enzyme cost minimisation or nonlinear versions of Flux Balance Analysis. Due to their common origin, the different approaches share mathematical optimality conditions of the same form. A general theory of optimal metabolic states, as proposed here, provides a logical link between existing modelling approaches and can help justify, interconvert, and combine metabolic optimality problems.

Genetic feedback loops can be used by cells to regulate internal processes or to keep track of time. It is often thought that, for a genetic circuit to display self-sustained oscillations, a degree of cooperativity is needed in the binding and unbinding of actor species. This cooperativity is usually modeled using a Hill function, regardless of the actual promoter architecture. Furthermore, genetic circuits do not operate in isolation and often transcription factors are shared between different promoters. In this work we show how mathematical modelling of genetic feedback loops can be facilitated with a mechanistic fold-change function that takes into account the titration effect caused by competing binding sites for transcription factors. The model shows how the titration effect facilitates self-sustained oscillations in a minimal genetic feedback loop: a gene that produces its own repressor directly without cooperative transcription factor binding. The use of delay-differential equations leads to a stability contour that predicts whether a genetic feedback loop will show self-sustained oscillations, even when taking the bursty nature of transcription into account.

The antibiotic trimethoprim targets the bacterial dihydrofolate reductase enzyme and subsequently affects the entire folate network. We present an expanded mathematical model of trimethoprims action on the Escherichia coli folate network that greatly improves upon Kwon et al. (2008). The improvement upon the Kwon Model lends greater insight into the effects of trimethoprim at higher resolution and accuracy. More importantly, the presented mathematical model enables drug target discovery in a way the earlier model could not. Using the improved mathematical model as a scaffold, we use parameter optimization to search for new drug targets that replicate the effect of trimethoprim. We present the model and model-scaffold strategy as an efficient route for drug target discovery.

Systems biology tackles the challenge of understanding the high complexity in the internal regulation of homeostasis in the human body through mathematical modelling. These models can aid in the discovery of disease mechanisms and potential drug targets. However, on one hand the development and validation of knowledge-based mechanistic models is time-consuming and does not scale well with increasing features in medical data. On the other hand, more data-driven approaches such as machine learning models require large volumes of data to produce generalizable models. The integration of neural networks and mechanistic models, forming universal differential equation (UDE) models, enables the automated learning of unknown model terms with less data than the neural network alone. Nevertheless, estimating parameters for these hybrid models remains difficult with sparse data and limited sampling durations that are common in biological applications. In this work, we propose the use of physiology-informed regularization, penalizing biologically implausible model behavior to guide the UDE towards more physiologically plausible regions of the solution space. In a simulation study we show that physiology-informed regularization not only results in a more accurate forecasting of model behaviour, but also supports training with less data. We also applied this technique to learn a representation of the rate of glucose appearance in the glucose minimal model using meal response data measured in healthy people. In that case, the inclusion of regularization reduces variability between UDE-embedded neural networks that were trained from different initial parameter guesses. Author summarySystems biology concerns the modelling and analysis of biological processes, by viewing these as interconnected systems. Modelling is typically done either using mechanistic differential equations that are derived from experiments and known biology, or using machine learning on large biological datasets. While mathematical modelling from biological experiments can provide useful insights with limited data, building and validating these models takes a long time and often requires highly invasive measurements in humans. Efforts to combine this classical technique with machine learning have resulted in a framework termed universal differential equations, where the model equations contain a neural network to describe unknown biological interactions. While these methods have shown success in numerous fields, applications in biology are more challenging due to limited data-availability, high data sparsity. In this work, we have introduced physiology-informed regularization to overcome these instabilities and to constrain the model to biologically plausible behavior. Our results show that by using physiology-informed regularization, we can accurately predict future unseen observations in a simulated example, with much more limited data than a similar model without regularization. Additionally, we show an application of this technique on human data, applying a neural network to learn the appearance of glucose in the blood plasma after a meal.

We present a method for modeling the dynamics of enzymatic chemical reaction networks using species-reaction graphs. These bipartite graphs contain two sets of vertices, one representing species and the other representing reactions, connected by directed edges that indicate relationships between them. Our approach starts by assigning appropriate edge weights to the bipartite graph, which are then used to compute the weighted graph Laplacian. This Laplacian allows us to reformulate the standard system of ordinary differential equations governing the network dynamics, emphasizing the flow of information - such as influences or dependencies between species and reactions - throughout the chemical reaction network considered as causal network. As an application of this framework, we introduce a novel model reduction technique based on the Kron reduction of the weighted Laplacian matrix associated with the species-reaction graphs. Our systematic approach involves identifying nodes for deletion while preserving the bipartite structure, followed by constructing the Kron-reduced model. To demonstrate the effectiveness of our method, we apply it to complex biochemical networks, showing how model simplification facilitates analysis and interpretation of these systems.

The Izhikevich artificial spiking neuron model is among the most employed models in neuromorphic engineering and computational neuroscience, due to the affordable computational effort to discretize it and its biological plausibility. It has been adopted also for applications with limited computational resources in embedded systems. It is important therefore to realize a compromise between error and computational expense to solve numerically the models equations. Here we investigate the effects of discretization and we study the solver that realizes the best compromise between accuracy and computational cost, given an available amount of Floating Point Operations per Second (FLOPS). We considered three fixed-step solvers for Ordinary Differential Equations (ODE), commonly used in computational neuroscience: Euler method, the Runge-Kutta 2 method and the Runge-Kutta 4 method. To quantify the error produced by the solvers, we used the Victor Purpura spike train Distance from an ideal solution of the ODE. Counterintuitively, we found that simple methods such as Euler and Runge Kutta 2 can outperform more complex ones (i.e. Runge Kutta 4) in the numerical solution of the Izhikevich model if the same FLOPS are allocated in the comparison. Moreover, we quantified the neuron rest time (with input under threshold resulting in no output spikes) necessary for the numerical solution to converge to the ideal solution and therefore to cancel the error accumulated during the spike train; in this analysis we found that the required rest time is independent from the firing rate and the spike train duration. Our results can generalize in a straightforward manner to other spiking neuron models and provide a systematic analysis of fixed step neural ODE solvers towards an accuracy-computational cost tradeoff.

Bioprocess development is commonly characterized by long development times, especially in the early screening phase. After promising candidates have been pre-selected in screening campaigns, an optimal operating strategy has to be found and verified under conditions similar to production. Cultivating cells with pulse-based feeding and thus exposing them to oscillating feast and famine phases has shown to be a powerful approach to study microorganisms closer to industrial bioreactor conditions. In view of the large number of strains and the process conditions to be tested, high-throughput cultivation systems provide an essential tool to sample the large design space in short time. We have recently presented a comprehensive platform, consisting of two liquid handling stations coupled with a model-based experimental design and operation framework to increase the efficiency in High Throughput bioprocess development. Using calibrated macro-kinetic growth models, the platform has been successfully used for the development of scale-down fed-batch cultivations in parallel mini-bioreactor systems. However, it has also been shown that parametric uncertainties in the models can significantly affect the prediction accuracy and thus the reliability of optimized cultivation strategies. To tackle this issue, we implemented a multi-stage Model Predictive Control (MPC) strategy to fulfill the experimental objectives under tight constraints despite the uncertainty in the parameters and the measurements. Dealing with uncertainties in the parameters is of major importance, since constraint violation would easily occur otherwise, which in turn could have adverse effects on the quality of the heterologous protein produced. Multi-stage approaches build up scenario tree, based on the uncertainty that can be encountered and computing optimal inputs that satisfy the constrains despite of such uncertainties. Using the feedback information gained through the evolution along the tree, the control approach is significantly more robust than standard MPC approaches without being overly conservative. We show in this study that the application of multi-stage MPC can increase the number of successful experiments, by applying this methodology to a mini-bioreactor cultivation operated in parallel.

Since its introduction by Briat, Gupta and Khammash, the antithetic feedback controller design has attracted considerable attention in both theoretical and experimental systems biology. The case in which the plant is a two-dimensional linear system (making the closed-loop system a nonlinear four-dimensional system) has been analyzed in much detail. This system has a unique equilibrium but, depending on parameters, it may exhibit periodic orbits. An interesting open question is whether other dynamical behaviors, such as chaotic attractors, might be possible for some parameter choices. This note shows that, for any parameter choices, every bounded trajectory satisfies a Poincare-Bendixson property. The analysis is based on the recently introduced notion of k-cooperative dynamical systems. It is shown that the model is a strongly 2-cooperative system, implying that the dynamics in the omega-limit set of any precompact solution is conjugate to the dynamics in a compact invariant subset of a two-dimensional Lipschitz dynamical system, thus precluding chaotic and other strange attractors.