Self-oscillations in glycolysis. 1. A simple kinetic model

On the coupling of intracellular K + ${{\rm{K}}}^{+}$ to glycolytic oscillations in yeast

We have investigated the interplay between glycolytic oscillations and intracellularK + ${{\rm{K}}}^{+}$ concentration in the yeast Saccharomyces cerevisiae. IntracellularK + ${{\rm{K}}}^{+}$ concentration was measured using the fluorophore potassium-binding benzofuranisophthalate (PBFI). We found thatK + ${{\rm{K}}}^{+}$ is an essential ion for the occurrence of glycolytic oscillations and that intracellularK + ${{\rm{K}}}^{+}$ concentration oscillates synchronously with other variables such as nicotinamide adenine dinucleotide hydride (NADH), intracellular adenosine triphosphate (ATP), and mitochondrial membrane potential. We also investigated if glycolysis and intracellularK + ${{\rm{K}}}^{+}$ concentration oscillate in a number of yeast strains with mutations inK + ${{\rm{K}}}^{+}$ transporters in the plasma membrane, mitochondrial membrane and in the vacuolar membrane. Most of these strains are still capable of showing glycolytic oscillations, but two strains are not: (i) a strain with a deletion in the mitochondrial Mdm38pK + ∕ H + ${{\rm{K}}}^{+}\unicode{x02215}{{\rm{H}}}^{+}$ transporter and (ii) a strain with deletion of the late endosomal Nhx1pK + ∕ H + ${{\rm{K}}}^{+}\unicode{x02215}{{\rm{H}}}^{+}$ (Na + ∕ H + ${\text{Na}}^{+}\unicode{x02215}{{\rm{H}}}^{+}$ ) transporter. In these two mutant strains intracellularK + ${{\rm{K}}}^{+}$ concentration seems to be low, indicating that the two transporters may be involved in transport ofK + ${{\rm{K}}}^{+}$ into the cytosol. In the strain, Mdm38pΔ ${\rm{\Delta }}$ oscillations in glycolysis could be restored by addition of theK + ∕ H + ${{\rm{K}}}^{+}\unicode{x02215}{{\rm{H}}}^{+}$ exchange ionophore nigericin. Furthermore, in two nonoscillating mutant strains with a defective V-ATPase and deletion of the Arp1p protein the intracellularK + ${{\rm{K}}}^{+}$ is relatively high, suggesting that the V-ATPase is essential for transport ofK + ${{\rm{K}}}^{+}$ out of the cytosol and that the cytoskeleton may be involved in bindingK + ${{\rm{K}}}^{+}$ to reduce the concentration of free ion in the cytosol. Analyses of the time series of oscillations of NADH, ATP, mitochondrial membrane potential, and potassium concentration using data-driven modeling corroborate the conjecture thatK + ${{\rm{K}}}^{+}$ ion is essential for the emergence of oscillations and support the experimental findings using mutant strains.

Nonlinear dynamics in phosphoinositide metabolism

Phosphoinositides broadly impact membrane dynamics, signal transduction and cellular physiology. The orchestration of signaling complexity by this seemingly simple metabolic pathway remains an open question. It is increasingly evident that comprehending the complexity of the phosphoinositides metabolic network requires a systems view based on nonlinear dynamics, where the products of metabolism can either positively or negatively modulate enzymatic function. These feedback and feedforward loops may be paradoxical, leading to counterintuitive effects. In this review, we introduce the framework of nonlinear dynamics, emphasizing distinct dynamical regimes such as the excitable state, oscillations, and mixed-mode oscillations-all of which have been experimentally observed in phosphoinositide metabolisms. We delve into how these dynamical behaviors arise from one or multiple network motifs, including positive and negative feedback loops, coherent and incoherent feedforward loops. We explore the current understanding of the molecular circuits responsible for these behaviors. While mapping these circuits presents both conceptual and experimental challenges, redefining cellular behavior based on dynamical state, lipid fluxes, time delay, and network topology is likely essential for a comprehensive understanding of this fundamental metabolic network.

Network switches and their role in circadian clocks

Circadian rhythms are generated by complex interactions among genes and proteins. Self-sustained ∼24 h oscillations require negative feedback loops and sufficiently strong nonlinearities that are the product of molecular and network switches. Here, we review common mechanisms to obtain switch-like behavior, including cooperativity, antagonistic enzymes, multisite phosphorylation, positive feedback, and sequestration. We discuss how network switches play a crucial role as essential components in cellular circadian clocks, serving as integral parts of transcription-translation feedback loops that form the basis of circadian rhythm generation. The design principles of network switches and circadian clocks are illustrated by representative mathematical models that include bistable systems and negative feedback loops combined with Hill functions. This work underscores the importance of negative feedback loops and network switches as essential design principles for biological oscillations, emphasizing how an understanding of theoretical concepts can provide insights into the mechanisms generating biological rhythms.

Kinetic modelling of glycolytic oscillations

Glycolytic oscillations have been studied for well over 60 years, but aspects of their function, and mechanisms of regulation and synchronisation remain unclear. Glycolysis is amenable to mechanistic mathematical modelling, as its components have been well characterised, and the system can be studied at many organisational levels: in vitro reconstituted enzymes, cell free extracts, individual cells, and cell populations. In recent years, the emergence of individual cell analysis has opened new ways of studying this intriguing system.

From biological data to oscillator models using SINDy

Periodic changes in the concentration or activity of different molecules regulate vital cellular processes such as cell division and circadian rhythms. Developing mathematical models is essential to better understand the mechanisms underlying these oscillations. Recent data-driven methods like SINDy have fundamentally changed model identification, yet their application to experimental biological data remains limited. This study investigates SINDy's constraints by directly applying it to biological oscillatory data. We identify insufficient resolution, noise, dimensionality, and limited prior knowledge as primary limitations. Using various generic oscillator models of different complexity and/or dimensionality, we systematically analyze these factors. We then propose a comprehensive guide for inferring models from biological data, addressing these challenges step by step. Our approach is validated using glycolytic oscillation data from yeast.

Clustering of chemically propelled nanomotors in chemically active environments

Synthetic nanomotors powered by chemical reactions have been designed to act as vehicles for active cargo transport, drug delivery, and a variety of other uses. Collections of such motors, acting in consort, can self-assemble to form swarms or clusters, providing opportunities for applications on various length scales. While such collective behavior has been studied when the motors move in a chemically inactive fluid environment, when the medium in which they move is a chemical network that supports complex spatial and temporal patterns, through simulation and theoretical analysis we show that collective behavior changes. Spatial patterns in the environment can guide and control motor collective states, and interactions of the motors with their environment can give rise to distinctive spatiotemporal motor patterns. The results are illustrated by studies of the motor dynamics in systems that support Turing patterns and spiral waves. This work is relevant for potential applications that involve many active nanomotors moving in complex chemical or biological environments.

Tuning limit cycles with a noise: Survival and collapse

We consider a general class of limit cycle oscillators driven by an additive Gaussian white noise. Based on the separation of timescales, we construct the equation of motion for slow dynamics after appropriate averaging over the fast motion. The equation for slow motion whose coefficients are modified by noise characteristics is solved to obtain the analytic solution in the long time limit. We show that with increase of noise strength, the loop area of the limit cycle decreases until a critical value is reached, beyond which the limit cycle collapses. We determine the noise threshold from the condition for removal of secular divergence of the perturbation series and work out two explicit examples of Van der Pol and Duffing-Van der Pol oscillators for corroboration between the theory and numerics.

Self-organization of active colloids mediated by chemical interactions

Self-propelled colloidal particles exhibit rich non-equilibrium phenomena and have promising applications in fields such as drug delivery and self-assembled active materials. Previous experimental and theoretical studies have shown that chemically active colloids that consume or produce a chemical can self-organize into clusters with diverse characteristics depending on the effective phoretic interactions. In this paper, we investigate self-organization in systems with multiple chemical species that undergo a network of reactions and multiple colloidal species that participate in different reactions. Active colloids propelled by complex chemical reactions with potentially nonlinear kinetics can be realized using enzymatic reactions that occur on the surface of enzyme-coated particles. To demonstrate how the self-organizing behavior depends on the chemical reactions active colloids catalyze and their chemical environment, we consider first a single type of colloid undergoing a simple catalytic reaction, and compare this often-studied case with self-organization in binary mixtures of colloids with sequential reactions, and binary mixtures with nonlinear autocatalytic reactions. Our results show that in general active colloids at low particle densities can form localized clusters in the presence of bulk chemical reactions and phoretic attractions. The characteristics of the clusters, however, depend on the reaction kinetics in the bulk and on the particles and phoretic coefficients. With one or two chemical species that only undergo surface reactions, the space for possible self-organizations are limited. By considering the additional system parameters that enter the chemical reaction network involving reactions on the colloids and in the fluid, the design space of colloidal self-organization can be enlarged, leading to a variety of non-equilibrium structures.

Oscillations in three-reaction quadratic mass-action systems

It is known that rank-two bimolecular mass-action systems do not admit limit cycles. With a view to understanding which small mass-action systems admit oscillation, in this paper we study rank-two networks with bimolecular source complexes but allow target complexes with higher molecularities. As our goal is to find oscillatory networks of minimal size, we focus on networks with three reactions, the minimum number that is required for oscillation. However, some of our intermediate results are valid in greater generality. One key finding is that an isolated periodic orbit cannot occur in a three-reaction, trimolecular, mass-action system with bimolecular sources. In fact, we characterize all networks in this class that admit a periodic orbit; in every case, all nearby orbits are periodic too. Apart from the well-known Lotka and Ivanova reactions, we identify another network in this class that admits a center. This new network exhibits a vertical Andronov-Hopf bifurcation. Furthermore, we characterize all two-species, three-reaction, bimolecular-sourced networks that admit an Andronov-Hopf bifurcation with mass-action kinetics. These include two families of networks that admit a supercritical Andronov-Hopf bifurcation and hence a stable limit cycle. These networks necessarily have a target complex with a molecularity of at least four, and it turns out that there are exactly four such networks that are tetramolecular.

Biphasic amplitude oscillator characterized by distinct dynamics of trough and crest

Biphasic amplitude dynamics (BAD) of oscillation have been observed in many biological systems. However, the specific topology structure and regulatory mechanisms underlying these biphasic amplitude dynamics remain elusive. Here, we searched all possible two-node circuit topologies and identified the core oscillator that enables robust oscillation. This core oscillator consists of a negative feedback loop between two nodes and a self-positive feedback loop of the input node, which result in the fast and slow dynamics of the two nodes, thereby achieving relaxation oscillation. Landscape theory was employed to study the stochastic dynamics and global stability of the system, allowing us to quantitatively describe the diverse positions and sizes of the Mexican hat. With increasing input strength, the size of the Mexican hat exhibits a gradual increase followed by a subsequent decrease. The self-activation of input node and the negative feedback on input node, which dominate the fast dynamics of the input node, were observed to regulate BAD in a bell-shaped manner. Both deterministic and statistical analysis results reveal that BAD is characterized by the linear and nonlinear dependence of the oscillation trough and crest on the input strength. In addition, combining with computational and theoretical analysis, we addressed that the linear response of trough to input is predominantly governed by the negative feedback, while the nonlinear response of crest is jointly regulated by the negative feedback loop and the self-positive feedback loop within the oscillator. Overall, this study provides a natural and physical basis for comprehending the occurrence of BAD in oscillatory systems, yielding guidance for the design of BAD in synthetic biology applications.

Sloppy model analysis identifies bifurcation parameters without normal form analysis

Bifurcation phenomena are common in multidimensional multiparameter dynamical systems. Normal form theory suggests that bifurcations are driven by relatively few combinations of parameters. Models of complex systems, however, rarely appear in normal form, and bifurcations are controlled by nonlinear combinations of the bare parameters of differential equations. Discovering reparameterizations to transform complex equations into a normal form is often very difficult, and the reparameterization may not even exist in a closed form. Here we show that information geometry and sloppy model analysis using the Fisher information matrix can be used to identify the combination of parameters that control bifurcations. By considering observations on increasingly long timescales, we find those parameters that rapidly characterize the system's topological inhomogeneities, whether the system is in normal form or not. We anticipate that this novel analytical method, which we call time-widening information geometry (TWIG), will be useful in applied network analysis.

Boundary-bound reactions: Pattern formation with and without hydrodynamics

We study chemical pattern formation in a fluid between two flat plates and the effect of such patterns on the formation of convective cells. This patterning is made possible by assuming the plates are chemically reactive or release reagents into the fluid, both of which we model as chemical fluxes. We consider this as a specific example of boundary-bound reactions. In the absence of coupling with fluid flow, we show that the two-reagent system with nonlinear reactions admits chemical instabilities equivalent to diffusion-driven Turing instabilities. In the other extreme, when chemical fluxes at the two bounding plates are constant, diffusion-driven instabilities do not occur but hydrodynamic phenomena analogous to Rayleigh-Bénard convection are possible. Assuming we can influence the chemical fluxes along the domain and select suitable reaction systems, this presents a mechanism for the control of chemical and hydrodynamic instabilities and pattern formation. We study a generic class of models and find necessary conditions for a bifurcation to pattern formation. Afterwards, we present two examples derived from the Schnakenberg-Selkov reaction. Unlike the classical Rayleigh-Bénard instability, which requires a sufficiently large unstable density gradient, a chemohydrodynamic instability based on Turing-style pattern formation can emerge from a state that is uniform in density. We also find parameter combinations that result in the formation of convective cells whether gravity acts upwards or downwards relative to the reactive plate. The wave number of the cells and the direction of the flow at regions of high/low concentration depend on the orientation, hence, different patterns can be elicited by simply inverting the device. More generally, our results suggest methods for controlling pattern formation and convection by tuning reaction parameters. As a consequence, we can drive and alter fluid flow in a chamber without mechanical pumps by influencing the chemical instabilities.

Collective catalysis under spatial constraints: Phase separation and size-scaling effects on mass action kinetics

Chemical reactions are usually studied under the assumption that both substrates and catalysts are well-mixed (WM) throughout the system. Although this is often applicable to test-tube experimental conditions, it is not realistic in cellular environments, where biomolecules can undergo liquid-liquid phase separation (LLPS) and form condensates, leading to important functional outcomes, including the modulation of catalytic action. Similar processes may also play a role in protocellular systems, like primitive coacervates, or in membrane-assisted prebiotic pathways. Here we explore whether the demixing of catalysts could lead to the formation of microenvironments that influence the kinetics of a linear (multistep) reaction pathway, as compared to a WM system. We implemented a general lattice model to simulate LLPS of a collection of different catalysts and extended it to include diffusion and a sequence of reactions of small substrates. We carried out a quantitative analysis of how the phase separation of the catalysts affects reaction times depending on the affinity between substrates and catalysts, the length of the reaction pathway, the system size, and the degree of homogeneity of the condensate. A key aspect underlying the differences reported between the two scenarios is that the scale invariance observed in the WM system is broken by condensation processes. The main theoretical implications of our results for mean-field chemistry are drawn, extending the mass action kinetics scheme to include substrate initial "hitting times" to reach the catalysts condensate. We finally test this approach by considering open nonlinear conditions, where we successfully predict, through microscopic simulations, that phase separation inhibits chemical oscillatory behavior, providing a possible explanation for the marginal role that this complex dynamic behavior plays in real metabolisms.

Conditions for the origin of homochirality in primordial catalytic reaction networks

We study the generation of homochirality in a general chemical model (based on the homogeneous, fully connected Smoluchowski aggregation-fragmentation model) that obeys thermodynamics and can be easily mapped onto known origin of life models (e.g. autocatalytic sets, hypercycles, etc.), with essential aspects of origin of life modeling taken into consideration. Using a combination of theoretical modeling and numerical simulations, we look for minimal conditions for which our general chemical model exhibits spontaneous mirror symmetry breaking. We show that our model spontaneously breaks mirror symmetry in various catalytic configurations that only involve a small number of catalyzed reactions and nothing else. Of particular importance is that mirror symmetry breaking occurs in our model without the need for single-step autocatalytis or mutual inhibition, which may be of relevance for prebiotic chemistry.

Prediction of Oscillations in Glycolysis in Ethanol-Consuming Erythrocyte-Bioreactors

A mathematical model of energy metabolism in erythrocyte-bioreactors loaded with alcohol dehydrogenase and acetaldehyde dehydrogenase was constructed and analyzed. Such erythrocytes can convert ethanol to acetate using intracellular NAD and can therefore be used to treat alcohol intoxication. Analysis of the model revealed that the rate of ethanol consumption by the erythrocyte-bioreactors increases proportionally to the activity of incorporated ethanol-consuming enzymes until their activity reaches a specific threshold level. When the ethanol-consuming enzyme activity exceeds this threshold, the steady state in the model becomes unstable and the model switches to an oscillation mode caused by the competition between glyceraldehyde phosphate dehydrogenase and ethanol-consuming enzymes for NAD. The amplitude and period of metabolite oscillations first increase with the increase in the activity of the encapsulated enzymes. A further increase in these activities leads to a loss of the glycolysis steady state, and a permanent accumulation of glycolytic intermediates. The oscillation mode and the loss of the steady state can lead to the osmotic destruction of erythrocyte-bioreactors due to an accumulation of intracellular metabolites. Our results demonstrate that the interaction of enzymes encapsulated in erythrocyte-bioreactors with erythrocyte metabolism should be taken into account in order to achieve the optimal efficacy of these bioreactors.

Autoregulation of Transcription and Translation: A Qualitative Analysis

The regulation of both mRNA transcription and translation by down-stream gene products allows for a range of rich dynamical behaviours (e.g. homeostatic, oscillatory, excitability and intermittent solutions). Here, qualitative analysis is applied to an existing model of a gene regulatory network in which a protein dimer inhibits its own transcription and upregulates its own translation rate. It is demonstrated that the model possesses a unique steady state, conditions are derived under which limit cycle solutions arise and estimates are provided for the oscillator period in the limiting case of a relaxation oscillator. The analysis demonstrates that oscillations can arise only if mRNA is more stable than protein and the effect of nonlinear translation inhibition is sufficiently strong. Moreover, it is shown that the oscillation period can vary non-monotonically with transcription rate. Thus the proposed framework can provide an explanation for observed species-specific dependency of segmentation clock period on Notch signalling activity. Finally, this study facilitates the application of the proposed model to more general biological settings where post transcriptional regulation effects are likely important.

Optimal enzyme utilization suggests that concentrations and thermodynamics determine binding mechanisms and enzyme saturations

Deciphering the metabolic functions of organisms requires understanding the dynamic responses of living cells upon genetic and environmental perturbations, which in turn can be inferred from enzymatic activity. In this work, we investigate the optimal modes of operation for enzymes in terms of the evolutionary pressure driving them toward increased catalytic efficiency. We develop a framework using a mixed-integer formulation to assess the distribution of thermodynamic forces and enzyme states, providing detailed insights into the enzymatic mode of operation. We use this framework to explore Michaelis-Menten and random-ordered multi-substrate mechanisms. We show that optimal enzyme utilization is achieved by unique or alternative operating modes dependent on reactant concentrations. We find that in a bimolecular enzyme reaction, the random mechanism is optimal over any other ordered mechanism under physiological conditions. Our framework can investigate the optimal catalytic properties of complex enzyme mechanisms. It can further guide the directed evolution of enzymes and fill in the knowledge gaps in enzyme kinetics.

Glycolytic Wave Patterns in a Simple Reaction-diffusion System with Inhomogeneous Influx: Dynamic Transitions

An inhomogeneous profile of chemostatted species generates a rich variety of patterns in glycolytic waves depicted in a Selkov reaction-diffusion framework here. A key role played by diffusion amplitude and symmetry in the chemostatted species profile in dictating the fate of local spatial dynamics involving periodic, quasiperiodic, and chaotic patterns and transitions among them are investigated systematically. More importantly, various dynamic transitions, including wave propagation direction changes, are illustrated in interesting situations. Besides numerical results, our analytical formulation of the amplitude equation connecting complex Ginzburg-Landau and Lambda-omega representation shed light on the phase dynamics of the system. This systematic study of the glycolytic reaction-diffusion wave is in line with previous experimental results in open spatial reactor and will provide a knowledge about the dynamics that shape and control biological information processing and related phenomena.

Numerical and renormalization group analysis of the phase diagram of a stochastic cubic autocatalytic reaction-diffusion system

The renormalization group is a set of tools that can be used to incorporate the effect of fluctuations in a dynamical system as a rescaling of the system's parameters. Here, we apply the renormalization group to a pattern-forming stochastic cubic autocatalytic reaction-diffusion model and compare its predictions with numerical simulations. Our results demonstrate a good agreement within the range of validity of the theory and show that external noise can be used as a control parameter in such systems.

Action functional gradient descent algorithm for estimating escape paths in stochastic chemical reaction networks

We first derive the Hamilton-Jacobi theory underlying continuous-time Markov processes, and then we use the construction to develop a variational algorithm for estimating escape (least improbable or first passage) paths for a generic stochastic chemical reaction network that exhibits multiple fixed points. The design of our algorithm is such that it is independent of the underlying dimensionality of the system, the discretization control parameters are updated toward the continuum limit, and there is an easy-to-calculate measure for the correctness of its solution. We consider several applications of the algorithm and verify them against computationally expensive means such as the shooting method and stochastic simulation. While we employ theoretical techniques from mathematical physics, numerical optimization and chemical reaction network theory, we hope that our work finds practical applications with an inter-disciplinary audience including chemists, biologists, optimal control theorists and game theorists.

Systematic analysis of negative and positive feedback loops for robustness and temperature compensation in circadian rhythms

Temperature compensation and robustness to biological noise are two key characteristics of the circadian clock. These features allow the circadian pacemaker to maintain a steady oscillation in a wide range of environmental conditions. The presence of a time-delayed negative feedback loop in the regulatory network generates autonomous circadian oscillations in eukaryotic systems. In comparison, the circadian clock of cyanobacteria is controlled by a strong positive feedback loop. Positive feedback loops with substrate depletion can also generate oscillations, inspiring other circadian clock models. What makes a circadian oscillatory network robust to extrinsic noise is unclear. We investigated four basic circadian oscillators with negative, positive, and combinations of positive and negative feedback loops to explore network features necessary for circadian clock resilience. We discovered that the negative feedback loop system performs the best in compensating temperature changes. We also show that a positive feedback loop can reduce extrinsic noise in periods of circadian oscillators, while intrinsic noise is reduced by negative feedback loops.

Spontaneous emergence of enantioenriched chiral aldol reaction products from Achiral precursors in solution and origin of biological homochirality of sugars: a first-principles study

Experimental reports about observation of spontaneous mirror symmetry breaking and chiral amplification in stereoselective Mannich and aldol reactions, run under fully achiral initial conditions, have drawn a lot of attention, fuelled partly by the role these reactions could have played in chemical evolution as a cause for still puzzling observed homochirality of biomolecules, often considered a prerequisite for the origin of life. We have now revisited this still unresolved problem, using DFT computation of all combinatorially possible transition states and numerical solution of complete set of resulting coupled kinetic rate equations to model the aldol reaction rigorously "from the first principles" and without making any a priori assumptions. Spontaneous mirror symmetry breaking in this autocatalytic, reversible, closed and homogenous system is explained by a supercritical pitchfork bifurcation, occurring in concentrations of enantiomers due to time-delayed kinetic instability of racemic composition of reaction mixture, when reactants are initially provided in non-stoichiometric quantities. Same process, taking place under similar conditions in primordial "soup" of chemicals, might conceivably explain origin of biological homochirality of sugar molecules on early earth billions of years ago. Our results suggest that seemingly innocuous chemical reactions could exhibit unexpected and counter-intuitive emergent behaviour, when initial conditions are appropriately chosen. Chiral amplification in self-catalyzed aldol reaction occurs during approach of thermodynamic equilibrium in accord with principle of microscopic reversibility and second law of thermodynamics.

Foundations of static and dynamic absolute concentration robustness

Absolute Concentration Robustness (ACR) was introduced by Shinar and Feinberg (Science 327:1389-1391, 2010) as robustness of equilibrium species concentration in a mass action dynamical system. Their aim was to devise a mathematical condition that will ensure robustness in the function of the biological system being modeled. The robustness of function rests on what we refer to as empirical robustness-the concentration of a species remains unvarying, when measured in the long run, across arbitrary initial conditions. Even simple examples show that the ACR notion introduced in Shinar and Feinberg (Science 327:1389-1391, 2010) (here referred to as static ACR) is neither necessary nor sufficient for empirical robustness. To make a stronger connection with empirical robustness, we define dynamic ACR, a property related to long-term, global dynamics, rather than only to equilibrium behavior. We discuss general dynamical systems with dynamic ACR properties as well as parametrized families of dynamical systems related to reaction networks. We find necessary and sufficient conditions for dynamic ACR in complex balanced reaction networks, a class of networks that is central to the theory of reaction networks.

Signaling oscillations: Molecular mechanisms and functional roles

Mounting evidence shows that oscillatory activity is widespread in cell signaling. Here, we review some of this recent evidence, focusing on both the molecular mechanisms that potentially underlie such dynamical behavior, and the potential advantages that signaling oscillations might have in cell function. The biological processes considered include cellular differentiation and tissue maintenance, intermittent responses in pluripotent stem cells, and collective cell migration during wound healing. With the aid of mathematical modeling, we review recent examples in which delayed negative feedback has been seen to act as a unifying principle that underpins this wide variety of phenomena.

Noise-induced formation of heterogeneous patterns in the Turing stability zones of diffusion systems

We study a phenomenon of stochastic generation of waveform patterns for reaction-diffusion systems in the Turing stability zone where the homogeneous equilibrium is a single attractor. In this analysis, we use a distributed variant of the Selkov glycolytic model with diffusion and random forcing. It is shown that in the Turing stability zone, random disturbances can induce a diversity of metastable spatial patterns with different waveforms. We carry out the parametric analysis of statistical characteristics of evolution of these patterns, and reveal the dominant patterns in the stochastic flow of mixed spatial structures.

Dynamical systems defined on simplicial complexes: Symmetries, conjugacies, and invariant subspaces

We consider the general model for dynamical systems defined on a simplicial complex. We describe the conjugacy classes of these systems and show how symmetries in a given simplicial complex manifest in the dynamics defined thereon, especially with regard to invariant subspaces in the dynamics.

Provenance of life: Chemical autonomous agents surviving through associative learning

We present a benchmark study of autonomous, chemical agents exhibiting associative learning of an environmental feature. Associative learning systems have been widely studied in cognitive science and artificial intelligence but are most commonly implemented in highly complex or carefully engineered systems, such as animal brains, artificial neural networks, DNA computing systems, and gene regulatory networks, among others. The ability to encode environmental information and use it to make simple predictions is a benchmark of biological resilience and underpins a plethora of adaptive responses in the living hierarchy, spanning prey animal species anticipating the arrival of predators to epigenetic systems in microorganisms learning environmental correlations. Given the ubiquitous and essential presence of learning behaviors in the biosphere, we aimed to explore whether simple, nonliving dissipative structures could also exhibit associative learning. Inspired by previous modeling of associative learning in chemical networks, we simulated simple systems composed of long- and short-term memory chemical species that could encode the presence or absence of temporal correlations between two external species. The ability to learn this association was implemented in Gray-Scott reaction-diffusion spots, emergent chemical patterns that exhibit self-replication and homeostasis. With the novel ability of associative learning, we demonstrate that simple chemical patterns can exhibit a broad repertoire of lifelike behavior, paving the way for in vitro studies of autonomous chemical learning systems, with potential relevance to artificial life, origins of life, and systems chemistry. The experimental realization of these learning behaviors in protocell or coacervate systems could advance a new research direction in astrobiology, since our system significantly reduces the lower bound on the required complexity for autonomous chemical learning.

Turing Instabilities and Rotating Spiral Waves in Glycolytic Processes

We study single-frequency oscillations and pattern formation in the glycolytic process modeled by a reduction in the well-known Sel'kov's equations (Sel'kov in Eur J Biochem 4:79, 1968), which describe, in the whole cell, the phosphofructokinase enzyme reaction. By using averaging theory, we establish the existence conditions for limit cycles and their limiting average radius in the kinetic reaction equations. We analytically establish conditions on the model parameters for the appearance of unstable nonlinear modes seeding the formation of two-dimensional patterns in the form of classical spots and stripes. We also establish the existence of a Hopf bifurcation, which characterizes the reaction dynamics, producing glycolytic rotating spiral waves. We numerically establish parameter regions for the existence of these spiral waves and address their linear stability. We show that as the model tends toward a suppression of the relative source rate, the spiral wave solution loses stability. All our findings are validated by full numerical simulations of the model equations. Finally, we discuss in vitro evidence of spatiotemporal activity patterns found in glycolytic experiments, and propose plausible biological implications of our model results.

Non-linearity of Metabolic Pathways Critically Influences the Choice of Machine Learning Model

The use of machine learning (ML) in life sciences has gained wide interest over the past years, as it speeds up the development of high performing models. Important modeling tools in biology have proven their worth for pathway design, such as mechanistic models and metabolic networks, as they allow better understanding of mechanisms involved in the functioning of organisms. However, little has been done on the use of ML to model metabolic pathways, and the degree of non-linearity associated with them is not clear. Here, we report the construction of different metabolic pathways with several linear and non-linear ML models. Different types of data are used; they lead to the prediction of important biological data, such as pathway flux and final product concentration. A comparison reveals that the data features impact model performance and highlight the effectiveness of non-linear models (e.g., QRF: RMSE = 0.021 nmol·min-1 and R2 = 1 vs. Bayesian GLM: RMSE = 1.379 nmol·min-1 R2 = 0.823). It turns out that the greater the degree of non-linearity of the pathway, the better suited a non-linear model will be. Therefore, a decision-making support for pathway modeling is established. These findings generally support the hypothesis that non-linear aspects predominate within the metabolic pathways. This must be taken into account when devising possible applications of these pathways for the identification of biomarkers of diseases (e.g., infections, cancer, neurodegenerative diseases) or the optimization of industrial production processes.

Coherence of oscillations in the weak-noise limit

In a noisy environment, oscillations lose their coherence, which can be characterized by a quality factor. We determine this quality factor for oscillations arising from a driven Fokker-Planck dynamics along a periodic one-dimensional potential analytically in the weak-noise limit. With this expression, we can prove for this continuum model the analog of an upper bound that has been conjectured for the coherence of oscillations in discrete Markov network models. We show that our approach can also be adapted to motion along a noisy two-dimensional limit cycle. Specifically, we apply our scheme to the noisy Stuart-Landau oscillator and the thermodynamically consistent Brusselator as a simple model for a chemical clock. Our approach thus complements the fairly sophisticated extant general framework based on techniques from Hamilton-Jacobi theory with which we compare our results numerically.

The emergence of polyglot entrainment responses to periodic inputs in vicinities of Hopf bifurcations in slow-fast systems

Several distinct entrainment patterns can occur in the FitzHugh-Nagumo (FHN) model under external periodic forcing. Investigating the FHN model under different types of periodic forcing reveals the existence of multiple disconnected 1:1 entrainment segments for constant, low enough values of the input amplitude when the unforced system is in the vicinity of a Hopf bifurcation. This entrainment structure is termed polyglot to distinguish it from the single 1:1 entrainment region (monoglot) structure typically observed in Arnold tongue diagrams. The emergence of polyglot entrainment is then explained using phase-plane analysis and other dynamical system tools. Entrainment results are investigated for other slow-fast systems of neuronal, circadian, and glycolytic oscillations. Exploring these models, we found that polyglot entrainment structure (multiple 1:1 regions) is observed when the unforced system is in the vicinity of a Hopf bifurcation and the Hopf point is located near a knee of a cubic-like nullcline.

On the Dynamics of Higgins–Selkov, Selkov and Brusellator Oscillators

A complete algebraic characterization of the first integrals of the Higgins–Selkov, Selkov and Brusellator oscillators is given here. The existence of symmetries sometimes forces the existence of such first integrals. The nonexistence of centers for such oscillators is also proved. In order to determine the Puiseux integrability of such systems, the multiple Puiseux solutions are also studied.

A systems biology analysis of lipolysis and fatty acid release from adipocytes in vitro and from adipose tissue in vivo

Lipolysis and the release of fatty acids to supply energy fuel to other organs, such as between meals, during exercise, and starvation, are fundamental functions of the adipose tissue. The intracellular lipolytic pathway in adipocytes is activated by adrenaline and noradrenaline, and inhibited by insulin. Circulating fatty acids are elevated in type 2 diabetic individuals. The mechanisms behind this elevation are not fully known, and to increase the knowledge a link between the systemic circulation and intracellular lipolysis is key. However, data on lipolysis and knowledge from in vitro systems have not been linked to corresponding in vivo data and knowledge in vivo. Here, we use mathematical modelling to provide such a link. We examine mechanisms of insulin action by combining in vivo and in vitro data into an integrated mathematical model that can explain all data. Furthermore, the model can describe independent data not used for training the model. We show the usefulness of the model by simulating new and more challenging experimental setups in silico, e.g. the extracellular concentration of fatty acids during an insulin clamp, and the difference in such simulations between individuals with and without type 2 diabetes. Our work provides a new platform for model-based analysis of adipose tissue lipolysis, under both non-diabetic and type 2 diabetic conditions.

Pattern forming systems coupling linear bulk diffusion to dynamically active membranes or cells

Some analytical and numerical results are presented for pattern formation properties associated with novel types of reaction-diffusion (RD) systems that involve the coupling of bulk diffusion in the interior of a multi-dimensional spatial domain to nonlinear processes that occur either on the domain boundary or within localized compartments that are confined within the domain. The class of bulk-membrane system considered herein is derived from an asymptotic analysis in the limit of small thickness of a thin domain that surrounds the bulk medium. When the bulk domain is a two-dimensional disk, a weakly nonlinear analysis is used to characterize Turing and Hopf bifurcations that can arise from the linearization around a radially symmetric, but spatially non-uniform, steady-state of the bulk-membrane system. In a singularly perturbed limit, the existence and linear stability of localized membrane-bound spike patterns is analysed for a Gierer-Meinhardt activator-inhibitor model that includes bulk coupling. Finally, the emergence of collective intracellular oscillations is studied for a class of PDE-ODE bulk-cell model in a bounded two-dimensional domain that contains spatially localized, but dynamically active, circular cells that are coupled through a linear bulk diffusion field. Applications of such coupled bulk-membrane or bulk-cell systems to some biological systems are outlined, and some open problems in this area are discussed. This article is part of the theme issue 'Recent progress and open frontiers in Turing's theory of morphogenesis'.

Unified framework for localized patterns in reaction-diffusion systems; the Gray-Scott and Gierer-Meinhardt cases

A recent study of canonical activator-inhibitor Schnakenberg-like models posed on an infinite line is extended to include models, such as Gray-Scott, with bistability of homogeneous equilibria. A homotopy is studied that takes a Schnakenberg-like glycolysis model to the Gray-Scott model. Numerical continuation is used to understand the complete sequence of transitions to two-parameter bifurcation diagrams within the localized pattern parameter regime as the homotopy parameter varies. Several distinct codimension-two bifurcations are discovered including cusp and quadruple zero points for homogeneous steady states, a degenerate heteroclinic connection and a change in connectedness of the homoclinic snaking structure. The analysis is repeated for the Gierer-Meinhardt system, which lies outside the canonical framework. Similar transitions are found under homotopy between bifurcation diagrams for the case where there is a constant feed in the active field, to it being in the inactive field. Wider implications of the results are discussed for other pattern-formation systems arising as models of natural phenomena. This article is part of the theme issue 'Recent progress and open frontiers in Turing's theory of morphogenesis'.

The finding and researching algorithm for potentially oscillating enzymatic systems

Many processes in living organisms are subject to periodic oscillations at different hierarchical levels of their organization: from molecular-genetic to population and ecological. Oscillatory processes are responsible for cell cycles in both prokaryotes and eukaryotes, for circadian rhythms, for synchronous coupling of respiration with cardiac contractions, etc. Fluctuations in the numbers of organisms in natural populations can be caused by the populations' own properties, their age structure, and ecological relationships with other species. Along with experimental approaches, mathematical and computer modeling is widely used to study oscillating biological systems. This paper presents classical mathematical models that describe oscillatory behavior in biological systems. Methods for the search for oscillatory molecular-genetic systems are presented by the example of their special case - oscillatory enzymatic systems. Factors influencing the cyclic dynamics in living systems, typical not only of the molecular-genetic level, but of higher levels of organization as well, are considered. Application of different ways to describe gene networks for modeling oscillatory molecular-genetic systems is considered, where the most important factor for the emergence of cyclic behavior is the presence of feedback. Techniques for finding potentially oscillatory enzymatic systems are presented. Using the method described in the article, we present and analyze, in a step-by-step manner, first the structural models (graphs) of gene networks and then the reconstruction of the mathematical models and computational experiments with them. Structural models are ideally suited for the tasks of an automatic search for potential oscillating contours (linked subgraphs), whose structure can correspond to the mathematical model of the molecular-genetic system that demonstrates oscillatory behavior in dynamics. At the same time, it is the numerical study of mathematical models for the selected contours that makes it possible to confirm the presence of stable limit cycles in them. As an example of application of the technology, a network of 300 metabolic reactions of the bacterium Escherichia coli was analyzed using mathematical and computer modeling tools. In particular, oscillatory behavior was shown for a loop whose reactions are part of the tryptophan biosynthesis pathway.

Alteration of Relative Rates of Biodegradation and Regeneration of Cervical Spine Cartilage through the Restoration of Arterial Blood Flow Access to Rhomboid Fossa: A Hypothesis

We found the logical way to prove the existence of the mechanism that maintains the rates of biodegradation and regeneration of cervical spine cartilage. We demonstrate, that after we restore access to arterial blood flow through cervical vertebral arteries to rhomboid fossa it causes the prevalence of regeneration over biodegradation. This is in the frames of consideration of the human body as a dissipative structure. Then the recovery of the body should be considered as a reduction of the relative rates of decay below the regeneration ones. Then the recovery of cervical spine cartilage through redirecting of inner dissipative flow depends on the information about oxygen availability that is provided from oxygen detectors in the rhomboid fossa to the cerebellum. Our proposed approach explains already collected data, which satisfies all the scientific requirements. This allows us to draw conclusions that permit reconsidering the way of dealing with multiple chronic diseases.

Synchronization of the Glycolysis Reaction-Diffusion Model via Linear Control Law

The Selkov system, which is typically employed to model glycolysis phenomena, unveils some rich dynamics and some other complex formations in biochemical reactions. In the present work, the synchronization problem of the glycolysis reaction-diffusion model is handled and examined. In addition, a novel convenient control law is designed in a linear form and, on the other hand, the stability of the associated error system is demonstrated through utilizing a suitable Lyapunov function. To illustrate the applicability of the proposed schemes, several numerical simulations are performed in one- and two-spatial dimensions.

A trajectory-based loss function to learn missing terms in bifurcating dynamical systems

Missing terms in dynamical systems are a challenging problem for modeling. Recent developments in the combination of machine learning and dynamical system theory open possibilities for a solution. We show how physics-informed differential equations and machine learning—combined in the Universal Differential Equation (UDE) framework by Rackauckas et al.—can be modified to discover missing terms in systems that undergo sudden fundamental changes in their dynamical behavior called bifurcations. With this we enable the application of the UDE approach to a wider class of problems which are common in many real world applications. The choice of the loss function, which compares the training data trajectory in state space and the current estimated solution trajectory of the UDE to optimize the solution, plays a crucial role within this approach. The Mean Square Error as loss function contains the risk of a reconstruction which completely misses the dynamical behavior of the training data. By contrast, our suggested trajectory-based loss function which optimizes two largely independent components, the length and angle of state space vectors of the training data, performs reliable well in examples of systems from neuroscience, chemistry and biology showing Saddle-Node, Pitchfork, Hopf and Period-doubling bifurcations.

Synchronisation of glycolytic activity in yeast cells

Glycolysis is the central metabolic pathway of almost every cell and organism. Under appropriate conditions, glycolytic oscillations may occur in individual cells as well as in entire cell populations or tissues. In many biological systems, glycolytic oscillations drive coherent oscillations of other metabolites, for instance in cardiomyocytes near anorexia, or in pancreas where they lead to a pulsatile release of insulin. Oscillations at the population or tissue level require the cells to synchronize their metabolism. We review the progress achieved in studying a model organism for glycolytic oscillations, namely yeast. Oscillations may occur on the level of individual cells as well as on the level of the cell population. In yeast, the cell-to-cell interaction is realized by diffusion-mediated intercellular communication via a messenger molecule. The present mini-review focuses on the synchronisation of glycolytic oscillations in yeast. Synchronisation is a quorum-sensing phenomenon because the collective oscillatory behaviour of a yeast cell population ceases when the cell density falls below a threshold. We review the question, under which conditions individual cells in a sparse population continue or cease to oscillate. Furthermore, we provide an overview of the pathway leading to the onset of synchronized oscillations. We also address the effects of spatial inhomogeneities (e.g., the formation of spatial clusters) on the collective dynamics, and also review the emergence of travelling waves of glycolytic activity. Finally, we briefly review the approaches used in numerical modelling of synchronized cell populations.

Nonequilibrium thermodynamics of glycolytic traveling wave: Benjamin-Feir instability

Evolution of the nonequilibrium thermodynamic entities corresponding to dynamics of the Hopf instabilities and traveling waves at a nonequilibrium steady state of a spatially extended glycolysis model is assessed here by implementing an analytically tractable scheme incorporating a complex Ginzburg-Landau equation (CGLE). In the presence of self and cross diffusion, a more general amplitude equation exploiting the multiscale Krylov-Bogoliubov averaging method serves as an essential tool to reveal the various dynamical instability criteria, especially Benjamin-Feir (BF) instability, to estimate the corresponding nonlinear dispersion relation of the traveling wave pattern. The critical control parameter, wave-number selection criteria, and magnitude of the complex amplitude for traveling waves are modified by self- and cross-diffusion coefficients within the oscillatory regime, and their variabilities are exhibited against the amplitude equation. Unlike the traveling waves, a low-amplitude broad region appears for the Hopf instability in the concentration dynamics as the system phase passes through minima during its variation with the control parameter. The total entropy production rate of the uniform Hopf oscillation and glycolysis wave not only qualitatively reflects the global dynamics of concentrations of intermediate species but almost quantitatively. Despite the crucial role of diffusion in generating and shaping the traveling waves, the diffusive part of the entropy production rate has a negligible contribution to the system's total entropy production rate. The Hopf instability shows a more complex and colossal change in the energy profile of the open nonlinear system than in the traveling waves. A detailed analysis of BF instability shows a contrary nature of the semigrand Gibbs free energy with discrete and continuous wave numbers for the traveling wave. We hope the Hopf and traveling wave pattern around the BF instability in terms of energetics and dissipation will open up new applications of such dynamical phenomena.

Triggering avalanches by transverse perturbations in a rotating drum

We study the role of small-scale perturbations in the onset of avalanches in a rotating drum in the stick-slip regime. By vibrating the system along the axis of rotation with an amplitude orders of magnitude smaller than the particles’ diameter, we found that the order parameter that properly describes the system is the kinetic energy. We also show that, for high enough frequencies, the onset of the avalanche is determined by the amplitude of the oscillation, contrary to previous studies that showed that either acceleration or velocity was the governing parameter. Finally, we present a theoretical model that explains the transition between the continuous and discrete avalanche regimes as a supercritical Hopf bifurcation.

Autocatalytic and oscillatory reaction networks that form guanidines and products of their cyclization

Autocatalytic and oscillatory networks of organic reactions are important for designing life-inspired materials and for better understanding the emergence of life on Earth; however, the diversity of the chemistries of these reactions is limited. In this work, we present the thiol-assisted formation of guanidines, which has a mechanism analogous to that of native chemical ligation. Using this reaction, we designed autocatalytic and oscillatory reaction networks that form substituted guanidines from thiouronium salts. The thiouronium salt-based oscillator show good stability of oscillations within a broad range of experimental conditions. By using nitrile-containing starting materials, we constructed an oscillator where the concentration of a bicyclic derivative of dihydropyrimidine oscillates. Moreover, the mixed thioester and thiouronium salt-based oscillator show unique responsiveness to chemical cues. The reactions developed in this work expand our toolbox for designing out-of-equilibrium chemical systems and link autocatalytic and oscillatory chemistry to the synthesis of guanidinium derivatives and the products of their transformations including analogs of nucleobases.

So far, only a few chemical oscillators based on organic reactions have been developed. Here, the authors report both autocatalytic and oscillatory reaction networks that form substituted guanidines from thiouronium salts; when coupled to cascade cyclization, this reaction network produces oscillations in the production of pyrimidine-based heterocycles.

The Emergence of the Bilateral Symmetry in Animals: A Review and a New Hypothesis

Most biological organisms exhibit different kinds of symmetry; an Animal (Metazoa), which is our Darwinist ancestor, has bilateral symmetry, and many plants exhibit rotational symmetry. It raises some questions: I. How can the evolution from an undifferentiated cell without bilateral symmetry to a complex biological organism with symmetry, which is based on asymmetric DNA and enzymes, lead to the bilateral symmetry? II. Is this evolution to an organism with bilateral symmetry obtained by other factors than DNA and enzymatic reactions? The existing literature about the evolution of the bilateral symmetry has been reviewed, and a new hypothesis has been formulated based on these reviews. The hypothesis is that the morphogenesis of biosystems is connected with the metabolism and that the oscillating kinetics in the Glycolysis have played a role in the polarity of the biological cells and in the establishment of the bilateral symmetry in Animals.

Unbounded solutions of models for glycolysis

The Selkov oscillator, a simple description of glycolysis, is a system of two ordinary differential equations with mass action kinetics. In previous work the authors established several properties of the solutions of this system. In the present paper we extend this to prove that this system has solutions which diverge to infinity in an oscillatory manner at late times. This is done with the help of a Poincaré compactification of the system and a shooting argument. This system was originally derived from another system with Michaelis-Menten kinetics. A Poincaré compactification of the latter system is carried out and this is used to show that the Michaelis-Menten system, like that with mass action, has solutions which diverge to infinity in a monotone manner. It is also shown to admit subcritical Hopf bifurcations and thus unstable periodic solutions. We discuss to what extent the unbounded solutions cast doubt on the biological relevance of the Selkov oscillator and compare it with other models for the same biological system in the literature.

In silico Determination of Some Conditions Leading to Glycolytic Oscillations and Their Interference With Some Other Processes in E. coli Cells

Autonomous oscillations of species levels in the glycolysis express the self-control of this essential cellular pathway belonging to the central carbon metabolism (CCM), and this phenomenon takes place in a large number of bacteria. Oscillations of glycolytic intermediates in living cells occur according to the environmental conditions and to the cell characteristics, especially the adenosine triphosphate (ATP) recovery system. Determining the conditions that lead to the occurrence and maintenance of the glycolytic oscillations can present immediate practical applications. Such a model-based analysis allows in silico (model-based) design of genetically modified microorganisms (GMO) with certain characteristics of interest for the biosynthesis industry, medicine, etc. Based on our kinetic model validated in previous works, this paper aims to in silico identify operating parameters and cell factors leading to the occurrence of stable glycolytic oscillations in the Escherichia coli cells. As long as most of the glycolytic intermediates are involved in various cellular metabolic pathways belonging to the CCM, evaluation of the dynamics and average level of its intermediates is of high importance for further applicative analyses. As an example, by using a lumped kinetic model for tryptophan (TRP) synthesis from literature, and its own kinetic model for the oscillatory glycolysis, this paper highlights the influence of glycolytic oscillations on the oscillatory TRP synthesis through the PEP (phosphoenolpyruvate) glycolytic node shared by the two oscillatory processes. The numerical analysis allows further TRP production maximization in a fed-batch bioreactor (FBR).