Catalyst.jl is a symbolic modeling package for analysis and high-performance simulation of chemical reaction networks and related dynamical systems. Models can be specified using an intuitive domain-specific language (DSL) or constructed programmatically. Catalyst supports ODE, steady-state ODE, SDE, stochastic chemical kinetics (jump), and hybrid simulations, including models that couple reactions with differential equations, events, and external noise (via Brownian Motions and/or Poisson Processes).

Built on ModelingToolkitBase.jl and Symbolics.jl, Catalyst leverages symbolic computation for sparsity exploitation, Jacobian construction, dependency graph analysis, and parallelism. Generated models integrate with the broader Julia and SciML ecosystems for sensitivity analysis, parameter estimation, bifurcation analysis, and more.

Catalyst can be installed as follows.

using Pkg

# (optional but recommended) create new environment in which to install Catalyst
Pkg.activate("catalyst_environment")

# install latest Catalyst release
Pkg.add("Catalyst")

Version 16 is a breaking release. Most breaking changes primarily affect libraries built on top of Catalyst. Please see HISTORY.md for the full list of breaking changes and migration guide.

Highlights:

• ModelingToolkitBase foundation — Without any reduction in core functionality, Catalyst now depends on ModelingToolkitBase instead of ModelingToolkit, to avoid adding new, non-MIT licensed dependencies.
• ModelingToolkit compatible - Catalyst is still compatible with ModelingToolkit for users who want to leverage more powerful, but non-MIT licensed, structural simplification libraries (e.g. for Catalyst-generated DAE models).
• Hybrid models — New HybridProblem and hybrid_model allow mixing ODE, SDE, and Jump reactions in a single system via per-reaction PhysicalScale metadata.
• Simplified jump API — Create jump problems directly with JumpProblem(rs, u0, tspan, ps), with no DiscreteProblem or JumpInputs intermediate.
• New DSL options — @brownians and @poissonians for coupling environmental noise, @discretes for event-modified parameters, @tstops for solver time stops, and => syntax for event affects.
• Modernized conversion API — ode_model, sde_model, jump_model, and ss_ode_model replace the old convert(ODESystem, rs) pattern, and all generate ModelingToolkitBase Systems.
• Unit validation — @unit_checks DSL option and validate_units/assert_valid_units functions with full support for non-SI units via DynamicQuantities.jl symbolic units.

The latest tutorials and information on using Catalyst are available in the stable documentation. The in-development documentation describes unreleased features in the current master branch.

An overview of the package, its features, and comparative benchmarking (as of version 13) can also be found in its corresponding research paper, Catalyst: Fast and flexible modeling of reaction networks.

• DSL for reaction networks — a readable, concise format for specifying models using chemical reaction notation.
• Multiple simulation types — generate and simulate ODE, steady-state ODE, SDE, jump, and hybrid models from a single ReactionSystem.
• Coupled models — combine reactions with differential equations, events, Brownian noise (@brownians), and Poisson jumps (@poissonians).
• Network analysis — compute linkage classes, deficiencies, reversibility, and other network properties.
• Compositional modeling — build models hierarchically using @network_component, compose, and extend.
• Spatial modeling — simulate reaction networks on discrete spatial domains.
• Steady state analysis — find and analyze steady states, stability, and bifurcation diagrams.
• Inverse problems — parameter estimation, sensitivity analysis, and structural identifiability.
• Model I/O — import from SBML and BioNetGen .net files, export to LaTeX and other formats.
• Visualization — reaction network graphs and LaTeX rendering.

Here we show a simple example where a model is created using the Catalyst DSL, and then simulated as an ordinary differential equation.

# Fetch required packages.
using Catalyst, OrdinaryDiffEqDefault, Plots

# Create model.
model = @reaction_network begin
    kB, S + E --> SE
    kD, SE --> S + E
    kP, SE --> P + E
end

# Create an ODE that can be simulated.
u0 = [:S => 50.0, :E => 10.0, :SE => 0.0, :P => 0.0]
tspan = (0., 200.)
ps = [:kB => 0.01, :kD => 0.1, :kP => 0.1]
ode = ODEProblem(model, u0, tspan, ps)

# Simulate ODE and plot results.
sol = solve(ode)
plot(sol; lw = 5)

The same model can be used as input to other types of simulations. E.g. here we instead generate and simulate a stochastic chemical kinetics jump process model for the reaction network. An exact realization of the jump process is sampled using an auto-selected stochastic simulation algorithm (SSA) (which for the small network in the current example ends up being Gillespie's Direct method):

# The initial conditions are now integers as we track exact populations for each species.
using JumpProcesses
u0_integers = [:S => 50, :E => 10, :SE => 0, :P => 0]
jprob = JumpProblem(model, u0_integers, tspan, ps)
jump_sol = solve(jprob)
plot(jump_sol; lw = 2)

This example demonstrates several Catalyst features composing together. We model a cell whose volume ($V$) grows proportionally to a phosphorylated growth factor ($G^P$), with environmental stochasticity ($\sigma,dW$) added via the @brownians DSL option. The phosphorylation of $G$ ($G \to G^P$) is driven by a cyclic sunlight signal $k_p(\sin(t)+1)$, and cell division occurs when the volume reaches a critical threshold $V_m$:

using Catalyst
cell_model = @reaction_network begin
    @parameters Vₘ g σ
    @brownians W
    @equations begin
        D(V) ~ g*Gᴾ + σ*W
    end
    @continuous_events begin
        [V ~ Vₘ] => [V => V/2]
    end
    kₚ*(sin(t)+1)/V, G --> Gᴾ
    kᵢ/V, Gᴾ --> G
end

We now study the system as a Chemical Langevin Dynamics SDE model:

u0 = [:V => 25.0, :G => 50.0, :Gᴾ => 0.0]
tspan = (0.0, 20.0)
ps = [:Vₘ => 50.0, :g => 0.3, :kₚ => 100.0, :kᵢ => 60.0, :σ => 0.5]
sprob = SDEProblem(cell_model, u0, tspan, ps)

This problem encodes the following stochastic differential equation model:

where $dW_1(t)$ and $dW_2(t)$ are the Chemical Langevin Equation noise terms from the reactions, and $dW(t)$ is an independent Brownian motion representing environmental stochasticity. Finally, we can simulate and plot the results:

using StochasticDiffEq, Plots
sol = solve(sprob, EM(); dt = 0.05)
plot(sol; xguide = "Time (au)", lw = 2)

Some features used here:

• Coupled differential equations modeled the cell volume alongside the reaction network.
• Events triggered cell division when the volume reached a threshold.
• @brownians added environmental noise to the volume equation.
• Specific solver and solver options were selected for the SDE simulation.
• The simulation was plotted using Plots.jl.

Catalyst integrates with a wide range of Julia packages:

Catalyst developers are active on the Julia Discourse and the Julia Slack channels #sciml-bridged and #sciml-sysbio. For bugs or feature requests, open an issue.

The software in this ecosystem was developed as part of academic research. If you would like to help support it, please star the repository as such metrics may help us secure funding in the future. If you use Catalyst as part of your research, teaching, or other activities, we would be grateful if you could cite our work:

@article{CatalystPLOSCompBio2023,
 doi = {10.1371/journal.pcbi.1011530},
 author = {Loman, Torkel E. AND Ma, Yingbo AND Ilin, Vasily AND Gowda, Shashi AND Korsbo, Niklas AND Yewale, Nikhil AND Rackauckas, Chris AND Isaacson, Samuel A.},
 journal = {PLOS Computational Biology},
 publisher = {Public Library of Science},
 title = {Catalyst: Fast and flexible modeling of reaction networks},
 year = {2023},
 month = {10},
 volume = {19},
 url = {https://doi.org/10.1371/journal.pcbi.1011530},
 pages = {1-19},
 number = {10},
}