Program
& Abstracts

We discuss Feller's renewal theorem and present an elementary proof.

Structural identifiability answers a theoretical question about which parameter combinations of a mathematical model are uniquely recovered from ideal, noise-free input–output data. Identifiability has been long studied on linear compartment models, and is closely related to the algebraic structure of the input-output equation obtained from the model. In this work, we investigate how structural identifiability is affected under singular limits in which a single parameter of a model is sent to infinity. Such limits induce a collapse of the underlying linear compartment graph that resembles an edge contraction and results in a reduced model. Using the transfer function and differential algebra approaches, we show that the transfer function of the original model converges to the transfer function of the collapsed model. From this result, we can identify which parameter combinations survive the singular limit and characterize the resulting loss of identifiability. We illustrate these results on several families of linear compartment models, including mammillary, catenary, and cyclic networks, and discuss connections to model reduction and the geometry of parameter space under asymptotic limits.

Can a test tube compute? And if so, can it compute reliably when its parameters are unknown? Can it compute with the same efficiency over a wide range of inputs? Recent advances in synthetic biology have made it possible to deploy chemical reactions that implement computation inside a cell. Chemical reactions function as analog computers — they are ideally suited for solving differential equations and interface naturally with real-world analog signals. However, exact arithmetic on an analog computer presents novel challenges, and existing approaches fail on both counts: (i) computation time can depend on the input, in some cases resulting in impractically slow computation, and (ii) rate constants are assumed to be perfectly known and fixed — an assumption that is unrealistic in practice, as rate constants are subject to physical, biochemical, and physiological variability and are difficult to measure accurately.

We will discuss novel algorithms that address both of these issues. First, we present algorithms that perform arithmetic operations — identification, inversion, addition, multiplication, absolute difference, rectified subtraction, and nth roots — at input-independent speed, with the striking property that computation time does not scale with the number of elementary steps — a consequence of the inherently parallel nature of chemical computation. We then extend this program to transcendental functions, constructing reaction network modules for the exponential and logarithm that achieve arbitrary accuracy at input-independent speed without relying on truncated power series — the logarithm presenting significantly deeper mathematical challenges. Finally, we present an independent line of results showing that arithmetic operations can be implemented by reaction networks whose outputs converge to the correct values regardless of the specific rate constants, through a novel design principle that automatically cancels rate-constant effects.

In this talk I will introduce our recent progress on the numerical estimation of the speed of convergence of a Poisson-driven chemical reaction system towards its steady state or its quasi-stationary distribution (QSD). Our approach is based on the coupling method, which is a well-established approach for stochastic differential equations. In order to couple chemical reaction systems on high dimensional lattices, we developed a novel algorithm based on the strong approximation between a Wiener process and a centralized Poisson process. The diffusion approximation of a chemical reaction system, driven by Wiener processes, can “guide” the Poisson-driven chemical reaction system to couple. Then the speed of convergence is implied by the coupling inequality.

We address the problem of estimating the state of a chemical reaction network based on observed time trajectories of some components of the state. In the deterministic setting, our goal is the design of an observer with a tunable parameter which converges exponentially fast to the system state trajectory. The exponential convergence is guaranteed for all initial conditions inside a ball whose diameter increases indefinitely as the tunable parameter is increased. We provide an observer design methodology that provably works for a certain class of networks which includes oscillatory and even chaotic networks. We illustrate our approach via several numerical examples. We also present numerical examples where observations are corrupted by noise. 

A closed-form expression for the infinitesimal generator of first-order stochastic chemical reaction networks is derived. Based on this representation, scaling limits are studied via the associated martingale problem. 

A Law of Large Numbers is obtained, yielding convergence to the deterministic reaction rate equations. Furthermore, a Central Limit Theorem is established, describing fluctuations around the deterministic limit in terms of a diffusion process. 

The results extend the framework of Enger and Pfaffelhuber (2021, arXiv:2111.15396) by providing a closed-form expression for the infinitesimal generator for more general chemical reaction networks, i.e., allowing an arbitrary number of reactions and an arbitrary number of slow and fast species.

Reaction networks arising in systems biology are often modeled with mass-action kinetics.  These dynamical systems may exhibit biologically significant properties, such as multistationarity and absolute concentration robustness (ACR). Multistationarity refers to the capacity for two or more steady states (mathematically, this means two or more positive real solutions to a system of polynomials), while ACR pertains to when some species concentration is the same at all positive steady states. This talk presents results on the prevalence of multistationarity and ACR in randomly generated reaction networks; in particular, we show that these two properties rarely occur together (mirroring what is seen in applications).  Additionally, we consider the question of whether and when multistationarity and/or ACR are preserved under operations that enlarge the reaction network. For instance, although many operations have been proven to preserve multistationarity, little is known about how much the number of steady states can increase after such an operation.  We present results on operations that double the number of steady states (and in some cases preserve ACR) — and even double the number of steady states that are exponentially stable.

This is joint work with Badal Joshi, Nidhi Kaihnsa, and Tung D. Nguyen.

Chemical reaction networks (CRNs) have previously been engineered to implement artificial neural networks (NNs) by leveraging the convergence of specific chemical species to desired outputs. However, past approaches often relied on discrete or non-smooth activation functions, which cause discontinuities in the backpropagation so that the training is prone to error in noisy environments. In contrast, our CRNs employ smooth activation functions, enabling both training and feed-forward computations to be carried out more reliably in the presence of noise. This advancement positions our work as a significant step forward in the evolution of CRN-based NNs.

Chemical reaction networks admit descriptions at multiple levels — stoichiometric, thermodynamic, kinetic, mechanistic — yet no principled account explains why each is necessary. We construct a canonical bottom-up tower of categorical structures in which each level is the unique minimal extension whose automorphism group separates networks the previous level conflates, unifying existing results in CRNT as special cases of a single forced hierarchy.

Cells often make discrete decisions from complex chemical inputs, such as when lectin signaling pathways decode high-dimensional patterns of glycan-bound receptor activity. These decisions are implemented by biochemical reaction networks, yet the structural features that determine their computational power remain unclear. To address this, we present a mathematical framework for understanding the computations performed by steady states of chemical reaction networks. We show that computational expressivity is governed by nonlinear kinetic features, as well as cross-talk among conserved moiety groups which we quantify as the number of directed interaction paths in the reaction graph. Our results reveal a simple design principle: classification performance can be enhanced either by increasing the diversity of receptor inputs or by increasing cross-talk within the downstream network. We numerically illustrate this by optimizing the rates of signaling networks, demonstrating that highly cross-linked networks can accurately separate chemical classes even in low-dimensional activation spaces, and furthermore such networks can flexibly support multiple classification tasks without changing kinetic parameters. These findings identify chemical cross-talk as a key structural feature underlying computational flexibility in biochemical networks.

Multistationarity underlies biochemical switching and cellular decision-making. In this talk we study how multistationarity in the sequential n-site phosphorylation–dephosphorylation cycle is affected when only some species are open, meaning allowed to exchange with the environment (so-called semi-open networks). For n>1, opening any nonempty subset of substrate species preserves nondegenerate multistationarity, while opening the enzymes (kinase and phosphatase), even together with substrates, always destroys it. The latter follows from a reduction that combines detection of absolute concentration robustness (ACR) with projection onto a monostationary subsystem; biologically, enzyme exchange acts as a robust “off” switch while substrate exchange preserves switching and modular coupling. We further illustrate the general method on multi-layer cascade variants.

Chemical Reaction Networks provide a fundamental framework for modeling the stochastic dynamics of biochemical systems, where molecular species evolve through discrete and random noise reaction events. Parameter inference in Chemical Reaction Networks is a central problem in systems biology, but traditional methods such as maximum likelihood estimation are often intractable due to computational complexity and the lack of continuous-time data. In this study, we introduce a statistically grounded and computationally efficient estimator for reaction rate parameters using high-frequency discrete-time observations. Modeling the system as a Continuous-Time Markov Chain, our method handles general kinetics, including non-mass-action and higher-order reactions. Validation on synthetic and experimental datasets demonstrates its accuracy and robustness. This approach offers a simple and reliable framework for parameter inference in complex stochastic systems.

We discuss why the diffusion approximation of a stochastic CRN model needs to be corrected if we want to extend its validity beyond the first hitting time to the boundary. At least two different solutions have been proposed, either by adding a reflection term, or by spending a waiting time at the boundary and jumping back in the interior, like the original MC model does. Available results and new perspectives will be discussed.

Stochastic reaction networks are mathematical models with a wide range of applications in biochemistry, ecology, and epidemiology, and are often complex to analyze. Except for some special cases, it is generally difficult to predict how the abundances of all considered species evolve over time. A possible approach to address this issue is to develop tools to compare the model under study with a similar one whose behavior is better understood. The main contribution of our work is to provide direct and computable conditions that can be used to ensure the existence of an ordered coupling between two stochastic reaction networks and to identify which parameter changes in a given model lead to an increase or decrease in the count of certain species. We also make an algorithm available that implements our theory and we illustrate it with several applications.

Glycolytic oscillations are a well-studied example of biological rhythms, yet their underlying biochemical networks are sufficiently complex that analysis often relies on reduced ordinary differential equation (ODE) models. Most such reductions are derived from deterministic descriptions, even though the true dynamics arise from stochastic reaction events. In this talk, I present a continuous-time Markov chain formulation of a classical glycolytic pathway model and show how a rigorous law of large numbers limit yields a reduced ODE system through multiscale stochastic averaging. I will also discuss parameter inference when only slow variables are observable and outline results establishing statistical consistency of the resulting estimates. Together, these results provide a mathematically grounded link between stochastic reaction networks, reduced dynamical models, and data-driven inference.

This is joint work with Arnab Ganguly.

Biochemical and biomolecular networks underpin information processing and decision making in cells. There are many sources of complexity of these networks. In this talk we will discuss a structured systems dissection of two facets of cellular reaction networks, using systems approaches, computational and analytical work.

The first part of the talk focusses on substrate modification. The reversible modification of substrates is a basic ingredient of cellular biochemical networks. These modification systems are of interest because they are a way of establishing protein function and also because of their potential to act as complex information processors. We focus on multiple modifications of a substrate, exploring and providing an overview of the effect of key associated ingredients such as the commonality or difference in enzymes effecting multiple modifications, the chemical mechanism of multiple modifications, presence/absence of ordering of modifications and the behaviour of these modification systems as part of networks. The focus is on different qualitative behaviours of these modification systems.

The second part of the talk focusses on spatial organization. Cells are far from being well-mixed entities and exhibit different forms of spatial compartmentalization/location. Recent interest in spatial organization has been catalyzed by the development of imaging technologies, spatial proteomics, ways of implementing compartmentalization in synthetic biology and the interest in membraneless compartments We examine the effects of spatial localization and compartmentalization by combining studies at the basic biochemistry level, studies at the network level and the examination of concrete exemplars. We discuss the implications for both systems and synthetic biology.

Toric dynamical systems show simple dynamics: their unique positive equilibrium attracts every positive solution. Motivated by the fact that a number of non-toric mass-action systems behave similarly neatly, disguised toric systems were introduced a couple of years ago. By definition, the mass-action differential equation of such a system is dynamically identical to that of a toric mass-action system. In many instances, this concept allows us to conclude the global stability of reaction networks for a set of parameter values significantly larger than their toric variety.

In this talk, we take another step forward. We say that a mass-action system is disguised toric in the general sense if there exists a coordinate change that makes it disguised toric. Even the simplest coordinate transformation, a linear diagonal change of variables, leads to a powerful method. We illustrate the strength of this new concept via several examples.

Joint work with Gheorghe Craciun, Oskar Henriksson, Jiaxin Jin, and Diego Rojas La Luz.

Recent work has revealed that biochemical networks with bifunctional enzymes can display remarkably rich dynamics, including ultrasensitivity, switch-like responses, concentration robustness, and even exhaustion of species. This talk presents a dynamical systems analysis of such networks which sheds light on the subtle architectural differences which produces these vast differences in functional behavior. We employ the next-generation matrix method---only recently adapted to biochemical reaction networks---to characterize previously incomputable thresholds for the stability of boundary steady states. These thresholds are critical for determining when a mechanism will proceed and when it will shut down. Using bifurcation analysis, we further establish conditions for multistationarity, showing how multiple positive steady states can arise within a single stoichiometric compatibility class. This properties is theorized to underlie toggle-switch behavior in genetic networks. We also consider the capacity of such networks for absolute concentration robustness--an essential feature of metabolic regulation in which a species maintains a fixed steady-state concentration despite fluctuations in other species---and explore the connection between robustness and boundary stability.

The steady state equations of chemical reaction networks are vertically parametrized systems under common kinetic assumptions, which, in many ways, make them extremely well-behaved from an algebraic point of view. In this talk, I'll follow up on my Sunday tutorial on algebraic reaction network theory by giving an introduction to the vertically parametrized perspective, and our latest results on what it mean for multistationarity, absolute concentration robustness and steady-state invariants. This is joint work with Elisenda Feliu and Beatriz Pascual Escudero, and is primarily based on the preprints 2412.17798 and 2411.15134.

A response curve measures the output of a biological system at equilibrium against an input parameter, which could be a rate constant, the dose of a ligand, or the strength of an external signal.  While often monotonic, a response curve that is biphasic serves as a mechanism against over-activation. For example, it has been observed that a T-cell's response to antigen concentration is non-monotonic. It has been conjectured that an incoherent feedforward loop is necessary to explain biphasic response. In this talk, we establish necessary conditions for biphasic response for general non-linear systems. More precisely, we proved that either an incoherent feedforward loop, or a combination of positive and negative feedback loops, is necessary for biphasic response. 

Identifying network structure and reaction parameters from stochastic data remains a central challenge in reaction network theory. This talk presents a likelihood-based framework that uses multinomial logistic regression to infer both stoichiometry and connectivity from time-series trajectories. When full state observations are available, stoichiometric vectors can be recovered under mild conditions on reaction occurrence and uniqueness. The approach is illustrated on several catalytic network models, including the Togashi–Kaneko system and an SIR epidemic model. To address partial observability, Bayesian logistic regression is combined with differential-equation modeling, enabling parameter recovery from infection-only data in a COVID-like SIR setting. Overall, the results highlight that relatively simple regression-based methods can yield reliable and interpretable insight into stochastic network structure and dynamics.

Based on joint work with Boseung Choi and Hye-Won Kang.

Stochastic chemical reaction networks satisfying topological conditions, such as weak reversibility and deficiency zero, and kinetic conditions like propensity factorizability, admit product-form stationary distributions. In addition, the network translation solves some violations of topological conditions whenever translated networks hold the kinetic condition.

In this talk, I will introduce the framework of the "dummy species extension" to solve some violations of kinetic conditions. By adding an additional species, originally non-propensity factorizable network systems that satisfy certain requirements are transformed into propensity factorizable network systems within this framework, enabling the analytic derivation of their stationary distributions. Furthermore, this study categorizes network systems and establishes the necessary conditions for a network system to be dummy representable.

I will present a rigorous mathematical framework for analyzing a class of stochastic copolymerization processes, where finitely many types of monomers attach and detach at the tip of a polymer chain. These dynamics are modeled as a continuous-time Markov chain on an infinite tree-like state space. Sharp criteria for transience, null recurrence, and positive recurrence in terms of the attachment and detachment rates are established. In the transient regime, explicit formulas for the almost sure asymptotic composition of the polymer and its growth velocity are provided. The framework also naturally extends to models allowing finite blocks of monomers to attach or detach and the block-memory case, and is expected to facilitate the analysis of more complex polymerization mechanisms in future work.

Many biological systems operate far from thermodynamic equilibrium by sustained energy consumption – here, we focus on a system of dynamic microtubules, energy-utilizing biopolymers made of the protein tubulin. At the molecular level, microtubule dynamics is governed by a stochastic chemical reaction network consisting of polymerization, depolymerization, and nucleotide hydrolysis reactions at filament ends. Irreversible GTP hydrolysis breaks detailed balance in the network, generating sustained probability currents that propagate to system-level observables and give rise to multiple distinct nonequilibrium steady states (NESS). To characterize these nonequilibrium steady states, we model the molecular-scale reaction network as a continuous-time Markov process and simulate its dynamics using exact stochastic simulation (Gillespie) methods. This framework provides access to both trajectory-level fluctuations and steady-state probability distributions under broken detailed balance. We then analyze the resulting steady states using stochastic thermodynamics, quantifying probability currents, entropy production, and system-level observables such as time to reach polymer-mass NESS (tNESS), concentration of polymerized tubulin at NESS ([polymerized tubulin]), emergent from the underlying reaction network. We find that tNESS exhibits systematic and nontrivial dependence on both free tubulin concentration and the rate of GTP hydrolysis, the latter controlling the strength of irreversibility in the molecular-scale reaction network. Altogether we show how local irreversibility at the molecular reaction-network level governs both – the transient state and the nonequilibrium steady states at the system level.

An important research problem for chemical reaction networks is determining the maximum number of steady states for a network. Many researchers used the mixed volume bound to give an upper bound for the maximum number of steady states, however, the Newton-Okounkov bound (introduced by Obatake and Walker) is known to be able to give sharper bounds in some cases. In this poster we investigate a class of networks for which the Newton-Okounkov bound will always be sharper than the mixed volume bound; we call such networks Full Volume. We show that checking whether a network is Full Volume can be done easily using the newton polytope, and introduce reduced networks, a class of networks that requires significantly fewer checks to verify Full Volume.

In chemical reaction network theory, a common area of study is the analysis of steady state solutions. Steady states can be described by a variety in high dimensional parameter space, which is partitioned into connected components that give different numbers of steady state solutions. Computing this discriminant can be computationally difficult as our chemical reaction networks increase in complexity, but numerical methods allows us to analyze the connected components of the complement of the discriminant without computing the defining equation.

In population dynamics, the Allee effect refers to the phenomenon where a population has a higher growth rate at higher densities. It has been thoroughly studied and until recently the number of steady states was not well understood. I will present a case study of how numerical methods can be used to describe the number of positive steady states of the Allee effect.

[1] Breiding, Paul, John Cobb, Aviva K. Englander, Nayda Farnsworth, Jonathan D. Hauenstein, Oskar Henriksson, David K. Johnson, Jordy Lopez Garcia, and Deepak Mundayur. "Elimination Without Eliminating: Computing Complements of Real Hypersurfaces Using Pseudo-Witness Sets." arXiv preprint arXiv:2601.04383 (2026).

[2] Englander, Aviva K. and Jose Israel Rodriguez. "Towards Learning the Positive Real Discriminant of the Wnt Signaling Pathway Shuttle Model." ACM Communications in Computer Algebra 58, no. 3 (2025): 85-88.

[3] Song, Kuo, and Xiaoxian Tang. "Steady State Classification of Allee Effect System." arXiv preprint arXiv:2501.19062 (2025).

Targeted protein degradation (TPD) systems exhibit non-monotonic dose responses, in which total target protein first decreases and then increases at high drug doses. In the present study, we analyzed a kinetic model of the degrader-target-ligase system using structural sensitivity analysis to determine whether non-monotonic dose response is permissible from the network topology alone. By using Cramer's rule to determine the sign structure of an augmented Jacobian matrix, we expressed the dose-response sensitivity of total target concentration as a ratio of determinants whose signs can be taken from the network structure. Further, we showed that the numerator of this ratio contains monomials of both signs, establishing that non-monotonic dose response is structurally permissible. We also performed subnetwork analysis, which revealed that the mixed signs require the formation of the ternary degrader-E3 ligase-target complex through two different, competing pathways, consistent with an incoherent feedforward loop. To verify that the non-monotonic dose response is not only structurally permissible but genuinely realizable, we ran numerical simulations with biologically feasible parameters to confirm the non-monotonic response curve. To ensure that the observed response was stable, we performed an eigenvalue analysis that confirmed that the unique steady state was stable across the full dosage range. Additionally, we performed parameter sweeps across key kinetic ratios to understand when the effect of non-monotonicity was most pronounced. Finally, we introduced a cooperativity factor for ternary complex formation and showed that it creates a tradeoff between efficacy and dosing robustness.

We study methods for the estimation of parametric sensitivities (derivatives of system outputs with respect to system parameters)  in stochastic reaction networks when multiple simultaneous parameter perturbations are required.   Application areas include: (i) multiple, simultaneous first order derivatives, (ii) higher-order methods aimed at substantially reducing bias for first order derivatives, and (iii) the estimation of higher-order derivatives. 

We introduce a multi-path stacked coupling (MSC) framework, extending the method of Anderson and Yuan (2019), to efficiently couple multiple nearby parameterized processes within a single framework. MSC provides a general-purpose approach and we demonstrate both analytically and numerically that it substantially reduces the computation cost in natural settings.  The MSC framework utilizes Poisson point processes as opposed to the more common random time change representation popularized by Kurtz.

Chemical reaction networks (CRNs) provide a natural model for analog computation in which inputs and outputs are encoded by molecular abundances. Classical studies of CRN-based computation have mainly focused on deterministic mass-action systems. Yet in biomolecular implementations using DNA and protein-based circuits, low copy-number effects and intrinsic fluctuations can be non-negligible, so stochastic models are often more faithful than concentration-based deterministic models. In this research, we develop a framework for computation in stochastic CRNs at the level of the mean of stationary distributions. Specifically, we construct elementary arithmetic modules such as identification, addition, and multiplication, and analyze their ergodicity and mixing times. We then study how these modules can be interconnected to form composite circuits and examine how their computational behavior interacts. Our results provide a systematic framework for computation in stochastic chemical reaction networks.