BSTModelKit.jl: A Julia Package for Constructing, Solving, and Analyzing Biochemical Systems Theory Models

Abstract

We present BSTModelKit.jl, an open-source Julia package for constructing, solving, and analyzing Biochemical Systems Theory (BST) models of biochemical networks. The package implements S-system representations, a canonical power-law formalism for modeling metabolic and regulatory networks. BSTModelKit.jl provides a declarative model specification format, dynamic simulation via ordinary differential equation (ODE) integration, steady-state computation, and global sensitivity analysis using the Morris and Sobol methods. The package leverages the Julia scientific computing ecosystem, in particular the SciML suite of differential equation solvers, to provide efficient and flexible model analysis tools. We describe the mathematical formulation, software design, and demonstrate the package capabilities with illustrative examples.

Figures (3)

• Figure 1: Dynamic simulation of a feedback-inhibited linear pathway driven by a pulsed input. Top panel: reaction network schematic showing a linear pathway from $X_1$ to $X_4$ with a branch to byproduct $X_5$; solid arrows denote mass flow catalyzed by enzymes $E_1$, $E_2$, and $E_3$, while the dashed red line indicates feedback inhibition of $r_1$ by the product $X_4$ ($g_{X_4,r_1} = -0.5$). Middle panel: the time-varying input function $u_1(t)$ applied to the production of $X_1$, consisting of a high pulse ($u_1 = 10$, $t = 5$--$15$) and a medium pulse ($u_1 = 5$, $t = 25$--$35$) separated by returns to baseline ($u_1 = 1$). Bottom panel: simulated concentration trajectories for all five species over $t \in [0, 50]$; $X_4$ accumulation during each pulse suppresses $r_1$ via feedback, limiting upstream buildup, while the byproduct $X_5$ accumulates monotonically because it lacks a degradation pathway.
• Figure 2: Steady-state analysis of a branched pathway under varying enzyme levels. Top panel: reaction network schematic showing a linear pathway from source to $C$, where the pathway branches to $D$ (via $r_3$, catalyzed by $E_3$) and $E$ (via $r_4$, catalyzed by $E_4$); both $D$ and $E$ are degraded by first-order reactions $r_5$ and $r_6$, respectively; dashed red lines indicate feedback inhibition of $r_1$ by $E$ ($g_{E,r_1} = -0.5$) and $r_2$ by $D$ ($g_{D,r_2} = -0.5$). Bottom panel: steady-state concentrations of all five species as a function of the enzyme level $E_3$ (swept from 0.1 to 5.0, with $E_4 = 1.0$ held fixed), computed using the steadystate function. At low $E_3$, flux is directed predominantly toward $E$; as $E_3$ increases, the $D$ and $E$ curves cross near $E_3 \approx 1$ and the flux redistributes to favor the $D$ branch. The upstream species $A$ and $B$ decrease at high $E_3$ due to modulation by the feedback loops.
• Figure 3: Global sensitivity analysis of the time-integrated $X_4$ concentration ($J = \int_0^{20} X_4\,dt$) for the feedback-inhibited pathway shown in \ref{['fig:feedback']}. Seven parameters were varied: six rate constants ($\alpha_{r_1}, \alpha_{r_2}, \alpha_{r_3}, \alpha_{r_5}, \alpha_{r_0}, \alpha_{r_4}$; bounds at $\pm 50\%$ of nominal) and the feedback kinetic order $g_{X_4,r_1}$ (varied from $-2$ to $0$). Left panel: Morris screening (500 trajectories) showing mean absolute elementary effects ($\mu^{*}$) with standard deviation error bars for each parameter; $\alpha_{r_0}$ and $\alpha_{r_4}$ dominate. Right panel: Sobol variance decomposition (1000 quasi-random samples) showing first-order ($S_1$) and total-order ($S_T$) sensitivity indices with confidence intervals; $\alpha_{r_4}$ ($S_1 \approx 0.54$, $S_T \approx 0.59$) and $\alpha_{r_0}$ ($S_1 \approx 0.38$, $S_T \approx 0.45$) account for the majority of output variance, while the gap between $S_1$ and $S_T$ indicates mild parameter interactions.