A Kinetic Finite Volume Discretization of the Multidimensional PIDE Model for Gene Regulatory Networks.

In this paper, a finite volume discretization scheme for partial integro-differential equations (PIDEs) describing the temporal evolution of protein distribution in gene regulatory networks is proposed. It is shown that the obtained set of ODEs can be formally represented as a compartmental kinetic system with a strongly connected reaction graph. This allows the application of the theory of nonnegative and compartmental systems for the qualitative analysis of the approximating dynamics. In this framework, it is straightforward to show the existence, uniqueness and stability of equilibria. Moreover, the computation of the stationary probability distribution can be traced back to the solution of linear equations. The discretization scheme is presented for one and multiple dimensional models separately. Illustrative computational examples show the precision of the approach, and good agreement with previous results in the literature.

IDESS: a toolbox for identification and automated design of stochastic gene circuits.

MOTIVATION: One of the main causes hampering predictability during the model identification and automated design of gene circuits in synthetic biology is the effect of molecular noise. Stochasticity may significantly impact the dynamics and function of gene circuits, specially in bacteria and yeast due to low mRNA copy numbers. Standard stochastic simulation methods are too computationally costly in realistic scenarios to be applied to optimization-based design or parameter estimation. RESULTS: In this work, we present IDESS (Identification and automated Design of Stochastic gene circuitS), a software toolbox for optimization-based design and model identification of gene regulatory circuits in the stochastic regime. This software incorporates an efficient approximation of the Chemical Master Equation as well as a stochastic simulation algorithm-both with GPU and CPU implementations-combined with global optimization algorithms capable of solving Mixed Integer Nonlinear Programming problems. The toolbox efficiently addresses two types of problems relevant in systems and synthetic biology: the automated design of stochastic synthetic gene circuits, and the parameter estimation for model identification of stochastic gene regulatory networks. AVAILABILITY AND IMPLEMENTATION: IDESS runs under the MATLAB environment and it is available under GPLv3 license at https://doi.org/10.5281/zenodo.7788692.

Automated Design of Synthetic Gene Circuits in the Presence of Molecular Noise.

Microorganisms (mainly bacteria and yeast) are frequently used as hosts for genetic constructs in synthetic biology applications. Molecular noise might have a significant effect on the dynamics of gene regulation in microbial cells, mainly attributed to the low copy numbers of mRNA species involved. However, the inclusion of molecular noise in the automated design of biocircuits is not a common practice due to the computational burden linked to the chemical master equation describing the dynamics of stochastic gene regulatory circuits. Here, we address the automated design of synthetic gene circuits under the effect of molecular noise combining a mixed integer nonlinear global optimization method with a partial integro-differential equation model describing the evolution of stochastic gene regulatory systems that approximates very efficiently the chemical master equation. We demonstrate the performance of the proposed methodology through a number of examples of relevance in synthetic biology, including different bimodal stochastic gene switches, robust stochastic oscillators, and circuits capable of achieving biochemical adaptation under noise.

Automated design of gene circuits with optimal mushroom-bifurcation behavior.

Recent advances in synthetic biology are enabling exciting technologies, including the next generation of biosensors, the rational design of cell memory, modulated synthetic cell differentiation, and generic multifunctional biocircuits. These novel applications require the design of gene circuits leading to sophisticated behaviors and functionalities. At the same time, designs need to be kept minimal to avoid compromising cell viability. Bifurcation theory addresses such challenges by associating circuit dynamical properties with molecular details of its design. Nevertheless, incorporating bifurcation analysis into automated design processes has not been accomplished yet. This work presents an optimization-based method for the automated design of synthetic gene circuits with specified bifurcation diagrams that employ minimal network topologies. Using this approach, we designed circuits exhibiting the mushroom bifurcation, distilled the most robust topologies, and explored its multifunctional behavior. We then outline potential applications in biosensors, memory devices, and synthetic cell differentiation.

Global Optimization Approach for Parameter Estimation in Stochastic Dynamic Models of Biosystems.

Mechanistic dynamic models have become an essential tool for understanding biomolecular networks and other biological systems. Biochemical stochasticity can be extremely important in some situations, e.g., at the single-cell level where there is a low copy number of the species involved. In these scenarios, deterministic models are not suitable to characterize the dynamics, so stochastic dynamic models should be considered. Here, we address the challenging problem of parameter estimation in stochastic dynamic models. Despite recent advances, this area is considerably less mature than its deterministic counterpart. We present a novel strategy based on two components: (i) global optimization via a hybrid stochastic-deterministic approach, and (ii) stochastic simulation techniques tailored to the sparsity of the available experimental data. Regarding the latter, for cases of dense population data we make use of a novel approach using a Partial Integro-Differential Equation (PIDE) model solved using a semilagrangian method. In order to further speed up the simulations, we also present efficient parallel implementations for multi-core CPUs and also for graphical processing units (GPUs). Importantly, whereas SDE and Fokker Planck approximations of the Chemical Master Equation (CME) apply when the reactant populations are sufficiently large, the PIDE approximation to the CME is valid for very low copy numbers, and therefore they enable us to tackle parameter estimation for systems with large intrinsic molecular noise, (highly stochastic regimes far from the thermodynamic limit). We test our strategy with four challenging problems: a Lotka-Volterra system, a polarization system in S. cerevisiae, a genetic toggle switch, and a genetic circadian oscillator. Our method could successfully solve these problems in very reasonable computation times (often a few minutes for the first two problems) using standard low-cost hardware, showing very significant speedups with respect to recent alternative methods. The code used to obtain the results reported here is available at https://doi.org/10.5281/zenodo.5195408.

Stochastic SIR model predicts the evolution of COVID-19 epidemics from public health and wastewater data in small and medium-sized municipalities: A one year study.

The level of unpredictability of the COVID-19 pandemics poses a challenge to effectively model its dynamic evolution. In this study we incorporate the inherent stochasticity of the SARS-CoV-2 virus spread by reinterpreting the classical compartmental models of infectious diseases (SIR type) as chemical reaction systems modeled via the Chemical Master Equation and solved by Monte Carlo Methods. Our model predicts the evolution of the pandemics at the level of municipalities, incorporating for the first time (i) a variable infection rate to capture the effect of mitigation policies on the dynamic evolution of the pandemics (ii) SIR-with-jumps taking into account the possibility of multiple infections from a single infected person and (iii) data of viral load quantified by RT-qPCR from samples taken from Wastewater Treatment Plants. The model has been successfully employed for the prediction of the COVID-19 pandemics evolution in small and medium size municipalities of Galicia (Northwest of Spain).

Wastewater and marine bioindicators surveillance to anticipate COVID-19 prevalence and to explore SARS-CoV-2 diversity by next generation sequencing: One-year study.

This study presents the results of SARS-CoV-2 surveillance in sewage water of 11 municipalities and marine bioindicators in Galicia (NW of Spain) from May 2020 to May 2021. An integrated pipeline was developed including sampling, pre-treatment and biomarker quantification, RNA detection, SARS-CoV-2 sequencing, mechanistic mathematical modeling and forecasting. The viral load in the inlet stream to the wastewater treatment plants (WWTP) was used to detect new outbreaks of COVID-19, and the data of viral load in the wastewater in combination with data provided by the health system was used to predict the evolution of the pandemic in the municipalities under study within a time horizon of 7 days. Moreover, the study shows that the viral load was eliminated from the treated sewage water in the WWTP, mainly in the biological reactors and the disinfection system. As a result, we detected a minor impact of the virus in the marine environment through the analysis of seawater, marine sediments and, wild and aquacultured mussels in the final discharge point of the WWTP.

OptDesign: Identifying Optimum Design Strategies in Strain Engineering for Biochemical Production.

Computational tools have been widely adopted for strain optimization in metabolic engineering, contributing to numerous success stories of producing industrially relevant biochemicals. However, most of these tools focus on single metabolic intervention strategies (either gene/reaction knockout or amplification alone) and rely on hypothetical optimality principles (e.g., maximization of growth) and precise gene expression (e.g., fold changes) for phenotype prediction. This paper introduces OptDesign, a new two-step strain design strategy. In the first step, OptDesign selects regulation candidates that have a noticeable flux difference between the wild type and production strains. In the second step, it computes optimal design strategies with limited manipulations (combining regulation and knockout), leading to high biochemical production. The usefulness and capabilities of OptDesign are demonstrated for the production of three biochemicals in Escherichia coli using the latest genome-scale metabolic model iML1515, showing highly consistent results with previous studies while suggesting new manipulations to boost strain performance. The source code is available at https://github.com/chang88ye/OptDesign.

Fast biofoundries: coping with the challenges of biomanufacturing.

Biofoundries are highly automated facilities that enable the rapid and efficient design, build, test, and learn cycle of biomanufacturing and engineering biology, which is applicable to both research and industrial production. However, developing a biofoundry platform can be expensive and time consuming. A biofoundry should grow organically, starting from a basic platform but with a vision for automation, equipment interoperability, and efficiency. By thinking about strategies early in the process through process planning, simulation, and optimization, bottlenecks can be identified and resolved. Here, we provide a survey of technological solutions in biofoundries and their advantages and limitations. We explore possible pathways towards the creation of a functional, early-phase biofoundry, and strategies towards long-term sustainability.

A numerical approach for detecting switch-like bistability in mass action chemical reaction networks with conservation laws.

BACKGROUND: Theoretical analysis of signaling pathways can provide a substantial amount of insight into their function. One particular area of research considers signaling pathways capable of assuming two or more stable states given the same amount of signaling ligand. This phenomenon of bistability can give rise to switch-like behavior, a mechanism that governs cellular decision making. Investigation of whether or not a signaling pathway can confer bistability and switch-like behavior, without knowledge of specific kinetic rate constant values, is a mathematically challenging problem. Recently a technique based on optimization has been introduced, which is capable of finding example parameter values that confer switch-like behavior for a given pathway. Although this approach has made it possible to analyze moderately sized pathways, it is limited to reaction networks that presume a uniterminal structure. It is this limited structure we address by developing a general technique that applies to any mass action reaction network with conservation laws. RESULTS: In this paper we developed a generalized method for detecting switch-like bistable behavior in any mass action reaction network with conservation laws. The method involves (1) construction of a constrained optimization problem using the determinant of the Jacobian of the underlying rate equations, (2) minimization of the objective function to search for conditions resulting in a zero eigenvalue, (3) computation of a confidence level that describes if the global minimum has been found and (4) evaluation of optimization values, using either numerical continuation or directly simulating the ODE system, to verify that a bistability region exists. The generalized method has been tested on three motifs known to be capable of bistability. CONCLUSIONS: We have developed a variation of an optimization-based method for the discovery of bistability, which is not limited to uniterminal chemical reaction networks. Successful completion of the method provides an S-shaped bifurcation diagram, which indicates that the network acts as a bistable switch for the given optimization parameters.

Automated Biocircuit Design with SYNBADm.

SYNBADm is a Matlab toolbox for the automated design of biocircuits using a model-based optimization approach. It enables the design of biocircuits with pre-defined functions starting from libraries of biological parts. SYNBADm makes use of mixed integer global optimization and allows both single and multi-objective design problems. Here we describe a basic protocol for the design of synthetic gene regulatory circuits. We illustrate step-by-step how to solve two different problems: (1) the (single objective) design of a synthetic oscillator and (2) the (multi-objective) design of a circuit with switch-like behavior upon induction, with a good compromise between performance and protein production cost.

Automated engineering of synthetic metabolic pathways for efficient biomanufacturing.

Metabolic engineering involves the engineering and optimization of processes from single-cell to fermentation in order to increase production of valuable chemicals for health, food, energy, materials and others. A systems approach to metabolic engineering has gained traction in recent years thanks to advances in strain engineering, leading to an accelerated scaling from rapid prototyping to industrial production. Metabolic engineering is nowadays on track towards a truly manufacturing technology, with reduced times from conception to production enabled by automated protocols for DNA assembly of metabolic pathways in engineered producer strains. In this review, we discuss how the success of the metabolic engineering pipeline often relies on retrobiosynthetic protocols able to identify promising production routes and dynamic regulation strategies through automated biodesign algorithms, which are subsequently assembled as embedded integrated genetic circuits in the host strain. Those approaches are orchestrated by an experimental design strategy that provides optimal scheduling planning of the DNA assembly, rapid prototyping and, ultimately, brings forward an accelerated Design-Build-Test-Learn cycle and the overall optimization of the biomanufacturing process. Achieving such a vision will address the increasingly compelling demand in our society for delivering valuable biomolecules in an affordable, inclusive and sustainable bioeconomy.

Synthetic Gene Circuit Analysis and Optimization.

Synthetic biology aims at engineering synthetic circuits with pre-defined target functions. From a systems (model-based) perspective, the following problems are of central importance: (1) given the model of a biomolecular circuit, elucidate whether it is capable of a certain behavior/functionality; and (2) starting from a pre-defined required functionality and a library of biological parts, find the biomolecular circuit that, built as a combination of such parts, achieves the desired behavior. These two problems, framed, respectively, as nonlinear analysis and automated design problems, are tackled here by efficient optimization methods. We illustrate these methods with case studies considering the analysis and design of biocircuits capable of bistability (bistable switches). Bistability is of particular interest in the context of systems and synthetic biology because it endows cells with the capacity to make decisions.

Solving the Competitive Binding Equilibria between Many Ligands: Application to High-Throughput Screening and Affinity Optimization.

Modern small-molecule drug discovery relies on the selective targeting of biological macromolecules by low-molecular weight compounds. Therefore, the binding affinities of candidate drugs to their targets are key for pharmacological activity and clinical use. For drug discovery methods where multiple drug candidates can simultaneously bind to the same target, a competition is established, and the resulting equilibrium depends on the dissociation constants and concentration of all the species present. Such coupling between all equilibrium-governing parameters complicates analysis and development of improved mixture-based, high-throughput drug discovery techniques. In this work, we present an iterative computational algorithm to solve coupled equilibria between an arbitrary number of ligands and a biomolecular target that is efficient and robust. The algorithm does not require the estimation of initial values to rapidly converge to the solution of interest. We explored binding equilibria under ligand/receptor conditions used in mixture-based library screening by affinity selection-mass spectrometry (AS-MS). Our studies support a facile method for affinity-ranking hits. The ranking method involves varying the receptor-to-ligand concentration ratio in a pool of candidate ligands in two sequential AS-MS analyses. The ranking is based on the relative change in bound ligand concentration. The method proposed does not require a known reference ligand and produces a ranking that is insensitive to variations in the concentration of individual compounds, thereby enabling the use of unpurified compounds generated by mixture-based combinatorial synthesis techniques.

CRNT4SBML: a Python package for the detection of bistability in biochemical reaction networks.

MOTIVATION: Signaling pathways capable of switching between two states are ubiquitous within living organisms. They provide the cells with the means to produce reversible or irreversible decisions. Switch-like behavior of biological systems is realized through biochemical reaction networks capable of having two or more distinct steady states, which are dependent on initial conditions. Investigation of whether a certain signaling pathway can confer bistability involves a substantial amount of hypothesis testing. The cost of direct experimental testing can be prohibitive. Therefore, constraining the hypothesis space is highly beneficial. One such methodology is based on chemical reaction network theory (CRNT), which uses computational techniques to rule out pathways that are not capable of bistability regardless of kinetic constant values and molecule concentrations. Although useful, these methods are complicated from both pure and computational mathematics perspectives. Thus, their adoption is very limited amongst biologists. RESULTS: We brought CRNT approaches closer to experimental biologists by automating all the necessary steps in CRNT4SMBL. The input is based on systems biology markup language (SBML) format, which is the community standard for biological pathway communication. The tool parses SBML and derives C-graph representations of the biological pathway with mass action kinetics. Next steps involve an efficient search for potential saddle-node bifurcation points using an optimization technique. This type of bifurcation is important as it has the potential of acting as a switching point between two steady states. Finally, if any bifurcation points are present, continuation analysis with respect to a user-defined parameter extends the steady state branches and generates a bifurcation diagram. Presence of an S-shaped bifurcation diagram indicates that the pathway acts as a bistable switch for the given optimization parameters. AVAILABILITY AND IMPLEMENTATION: CRNT4SBML is available via the Python Package Index. The documentation can be found at https://crnt4sbml.readthedocs.io. CRNT4SBML is licensed under the Apache Software License 2.0.

BioSwitch: a tool for the detection of bistability and multi-steady state behaviour in signalling and gene regulatory networks.

MOTIVATION: Multi-steady state behaviour, and in particular multi-stability, provides biological systems with the capacity to take reliable decisions (such as cell fate determination). A problem arising frequently in systems biology is to elucidate whether a signal transduction mechanism or a gene regulatory network has the capacity for multi-steady state behaviour, and consequently for a switch-like response to stimuli. Bifurcation diagrams are a powerful instrument in non-linear analysis to study the qualitative and quantitative behaviour of equilibria including bifurcation into different equilibrium branches and bistability. However, in the context of signalling pathways, the inherent large parametric uncertainty hampers the (direct) use of standard bifurcation tools. RESULTS: We present BioSwitch, a toolbox to detect multi-steady state behaviour in signalling pathways and gene regulatory networks. The tool combines results from chemical reaction network theory with global optimization to efficiently detect whether a signalling pathway has the capacity to undergo a saddle node bifurcation, and in case of multi-stationarity, provides the exact coordinates of the bifurcation where to start a numerical continuation analysis with standard bifurcation tools, leading to two different branches of equilibria. Bistability detection in the G1/S transition pathway of Saccharomyces cerevisiae is included as an illustrative example. AVAILABILITY AND IMPLEMENTATION: BioSwitch runs under the popular MATLAB computational environment, and is available at https://sites.google.com/view/bioswitch.

SELANSI: a toolbox for simulation of stochastic gene regulatory networks.

Motivation: Gene regulation is inherently stochastic. In many applications concerning Systems and Synthetic Biology such as the reverse engineering and the de novo design of genetic circuits, stochastic effects (yet potentially crucial) are often neglected due to the high computational cost of stochastic simulations. With advances in these fields there is an increasing need of tools providing accurate approximations of the stochastic dynamics of gene regulatory networks (GRNs) with reduced computational effort. Results: This work presents SELANSI (SEmi-LAgrangian SImulation of GRNs), a software toolbox for the simulation of stochastic multidimensional gene regulatory networks. SELANSI exploits intrinsic structural properties of gene regulatory networks to accurately approximate the corresponding Chemical Master Equation with a partial integral differential equation that is solved by a semi-lagrangian method with high efficiency. Networks under consideration might involve multiple genes with self and cross regulations, in which genes can be regulated by different transcription factors. Moreover, the validity of the method is not restricted to a particular type of kinetics. The tool offers total flexibility regarding network topology, kinetics and parameterization, as well as simulation options. Availability and implementation: SELANSI runs under the MATLAB environment, and is available under GPLv3 license at https://sites.google.com/view/selansi. Contact: antonio@iim.csic.es.

Chemical Reaction Network Theory elucidates sources of multistability in interferon signaling.

Bistability has important implications in signaling pathways, since it indicates a potential cell decision between alternative outcomes. We present two approaches developed in the framework of the Chemical Reaction Network Theory for easy and efficient search of multiple steady state behavior in signaling networks (both with and without mass conservation), and apply them to search for sources of bistability at different levels of the interferon signaling pathway. Different type I interferon subtypes and/or doses are known to elicit differential bioactivities (ranging from antiviral, antiproliferative to immunomodulatory activities). How different signaling outcomes can be generated through the same receptor and activating the same JAK/STAT pathway is still an open question. Here, we detect bistability at the level of early STAT signaling, showing how two different cell outcomes are achieved under or above a threshold in ligand dose or ligand-receptor affinity. This finding could contribute to explain the differential signaling (antiviral vs apoptotic) depending on interferon dose and subtype (α vs ß) observed in type I interferons.

Stochastic modeling and numerical simulation of gene regulatory networks with protein bursting.

Gene expression is inherently stochastic. Advanced single-cell microscopy techniques together with mathematical models for single gene expression led to important insights in elucidating the sources of intrinsic noise in prokaryotic and eukaryotic cells. In addition to the finite size effects due to low copy numbers, translational bursting is a dominant source of stochasticity in cell scenarios involving few short lived mRNA transcripts with high translational efficiency (as is typically the case for prokaryotes), causing protein synthesis to occur in random bursts. In the context of gene regulation cascades, the Chemical Master Equation (CME) governing gene expression has in general no closed form solution, and the accurate stochastic simulation of the dynamics of complex gene regulatory networks is a major computational challenge. The CME associated to a single gene self regulatory motif has been previously approximated by a one dimensional time dependent partial integral differential equation (PIDE). However, to the best of our knowledge, multidimensional versions for such PIDE have not been developed yet. Here we propose a multidimensional PIDE model for regulatory networks involving multiple genes with self and cross regulations (in which genes can be regulated by different transcription factors) derived as the continuous counterpart of a CME with jump process. The model offers a reliable description of systems with translational bursting. In order to provide an efficient numerical solution, we develop a semilagrangian method to discretize the differential part of the PIDE, combined with a composed trapezoidal quadrature formula to approximate the integral term. We apply the model and numerical method to study sustained stochastic oscillations and the development of competence, a particular case of transient differentiation attained by certain bacterial cells under stress conditions. We found that the resulting probability distributions are distinguishable from those characteristic of other transient differentiation processes. In this way, they can be employed as markers or signatures that identify such phenomena from bacterial population experimental data, for instance. The computational efficiency of the semilagrangian method makes it suitable for purposes like model identification and parameter estimation from experimental data or, in combination with optimization routines, the design of gene regulatory networks under molecular noise.

Automated Design Framework for Synthetic Biology Exploiting Pareto Optimality.

In this work we consider Pareto optimality for automated design in synthetic biology. We present a generalized framework based on a mixed-integer dynamic optimization formulation that, given design specifications, allows the computation of Pareto optimal sets of designs, that is, the set of best trade-offs for the metrics of interest. We show how this framework can be used for (i) forward design, that is, finding the Pareto optimal set of synthetic designs for implementation, and (ii) reverse design, that is, analyzing and inferring motifs and/or design principles of gene regulatory networks from the Pareto set of optimal circuits. Finally, we illustrate the capabilities and performance of this framework considering four case studies. In the first problem we consider the forward design of an oscillator. In the remaining problems, we illustrate how to apply the reverse design approach to find motifs for stripe formation, rapid adaption, and fold-change detection, respectively.