Repository files navigation

Catalyst.jl is a symbolic modeling package for analysis and high-performancesimulation of chemical reaction networks. Catalyst defines symbolicReactionSystems,which can be created programmatically or easilyspecified using Catalyst's domain-specific language (DSL). LeveragingModelingToolkit.jl andSymbolics.jl, Catalyst enableslarge-scale simulations through auto-vectorization and parallelism. SymbolicReactionSystems can be used to generate ModelingToolkit-based models, allowingthe easy simulation and parameter estimation of mass action ODE models, ChemicalLangevin SDE models, stochastic chemical kinetics jump process models, and more.Generated models can be used with solvers throughout the broader Julia andSciML ecosystems, including higher-level SciML packages (e.g.for sensitivity analysis, parameter estimation, machine learning applications,etc).

Catalyst can be installed as follows.

using Pkg# (optional but recommended) create new environment in which to install CatalystPkg.activate("catalyst_environment")# install latest Catalyst releasePkg.add("Catalyst")

NOTE: Version 15 is a breaking release, though most breaking changes arelikely to only impact libraries being developed on top of Catalyst. Pleasesee theHISTORY.md file for a summary of breaking changes and newfunctionality.

The latest tutorials and information on using Catalyst are available in thestabledocumentation. Thein-developmentdocumentation describes unreleased features inthe current master branch.

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

• The Catalyst DSL provides a simple and readable format for manually specifying reaction network models using chemical reaction notation.
• CatalystReactionSystems provides a symbolic representation of reaction networks, built onModelingToolkit.jl andSymbolics.jl.
• TheCatalyst.jl API provides functionality for building networks programmatically and for composing multiple networks together.
• Leveraging ModelingToolkit, generated models can be converted to symbolic reaction rate equation ODE models, symbolic Chemical Langevin Equation models, and symbolic stochastic chemical kinetics (jump process) models. These can be simulated using anyDifferentialEquations.jlODE/SDE/jump solver, and can be used withinEnsembleProblems for carrying outparallelized parameter sweeps and statistical sampling. Plot recipes are available forvisualization of all solutions.
• Non-integer (e.g.Float64) stoichiometric coefficientsare supported for generating ODE models, and symbolic expressions for stoichiometric coefficientsare supported for all system types.
• Anetwork analysis suite permits the computation of linkage classes, deficiencies, reversibility, and other network properties.
• Conservation laws can be detected and utilized to reduce system sizes, and to generate non-singular Jacobians (e.g. during conversion to ODEs, SDEs, and steady state equations).
• Catalyst reaction network models can becoupled with differential and algebraic equations (which are then incorporated during conversion to ODEs, SDEs, and steady state equations).
• Models can becoupled with events that affect the system and its state during simulations.
• By leveraging ModelingToolkit, users have a variety of options for generating optimized system representations to use in solvers. These include construction ofdense or sparse Jacobians,multithreading or parallelization of generated derivative functions,automatic classification of reactions into optimized jump types for Gillespie type simulations,automatic construction of dependency graphs for jump systems, and more.
• Symbolics.jl symbolic expressions and JuliaExprs can be obtained for all rate laws and functions determining the deterministic and stochastic terms within resulting ODE, SDE, or jump models.
• Steady states (and theirstabilities) can be computed for model ODE representations.

• OrdinaryDiffEq.jl Can be used to numerically solver generated reaction rate equation ODE models.
• StochasticDiffEq.jl can be used to numerically solve generated Chemical Langevin Equation SDE models.
• JumpProcesses.jl can be used to numerically sample generated Stochastic Chemical Kinetics Jump Process models.
• Support forparallelization of all simulations, including parallelization ofODE andSDE simulations on GPUs usingDiffEqGPU.jl.
• Latexify can be used togenerate LaTeX expressions corresponding to generated mathematical models or the underlying set of reactions.
• GraphMakie can be used to generate andvisualize reaction network graphs.
• Model steady states can becomputed through homotopy continuation usingHomotopyContinuation.jl (which can findall steady states of systems with multiple ones), byforward ODE simulations usingSteadyStateDiffEq.jl, or bynumerically solving steady-state nonlinear equations usingNonlinearSolve.jl.
• BifurcationKit.jl can be used tocompute bifurcation diagrams of model steady states (including finding periodic orbits).
• DynamicalSystems.jl can be used to compute modelbasins of attraction,Lyapunov spectrums, and other dynamical system properties.
• Optimization.jl andPEtab.jl can all be used tofit model parameters to data.
• GlobalSensitivity.jl can be used to performglobal sensitivity analysis of model behaviors.
• SciMLSensitivity.jl can be used to compute local sensitivities of functions containing forward model simulations.
• StructuralIdentifiability.jl can be used toperform structural identifiability analysis.

• CatalystReactionSystems can beimported from SBML files viaSBMLImporter.jl andSBMLToolkit.jl, andfrom BioNetGen .net files and various stoichiometric matrix network representations usingReactionNetworkImporters.jl.
• MomentClosure.jl allows generation of symbolic ModelingToolkitODESystems that represent moment closure approximations to moments of the Chemical Master Equation, from reaction networks defined in Catalyst.
• FiniteStateProjection.jl allows the construction and numerical solution of Chemical Master Equation models from reaction networks defined in Catalyst.
• DelaySSAToolkit.jl can augment Catalyst reaction network models with delays, and can simulate the resulting stochastic chemical kinetics with delays models.

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

# Fetch required packages.using Catalyst, OrdinaryDiffEqDefault, Plots# Create model.model=@reaction_networkbegin    kB, S+ E--> SE    kD, SE--> S+ E    kP, SE--> P+ Eend# 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 weinstead generate and simulate a stochastic chemical kinetics jump process modelfor the reaction network. An exact realization of the jump process is sampledusing an auto-selected stochastic simulation algorithm (SSA) (which for thesmall 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 JumpProcessesu0_integers= [:S=>50,:E=>10,:SE=>0,:P=>0]jinput=JumpInputs(model, u0_integers, tspan, ps)jprob=JumpProblem(jinput)jump_sol=solve(jprob)plot(jump_sol; lw=2)

In the above example, we used basic Catalyst workflows to simulate a simplemodel. Here we instead show how various Catalyst features can compose to createa much more advanced model. Our model describes how the volume of a cell ($V$)is affected by a growth factor ($G$). The growth factor only promotes growthwhile in its phosphorylated form ($G^P$). The phosphorylation of$G$ ($G \to G^P$)is promoted by sunlight, which is modeled as the cyclic sinusoid$k_a (\sin(t) + 1)$.When the cell reaches a critical volume ($V_m$) it undergoes cell division. First, wedeclare our model:

using Catalystcell_model=@reaction_networkbegin@parameters Vₘ g@equationsbeginD(V)~ g*Gᴾend@continuous_eventsbegin        [V~ Vₘ]=> [V~ V/2]end    kₚ*(sin(t)+1)/V, G--> Gᴾ    kᵢ/V, Gᴾ--> Gend

We now study the system as a Chemical Langevin Dynamics SDE model, which can be generated as follows

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]sprob=SDEProblem(cell_model, u0, tspan, ps)

This problem encodes the following stochastic differential equation model:

where the$dW_1(t)$ and$dW_2(t)$ terms represent independent Brownian Motions, encoding the noise added by the Chemical Langevin Equation. Finally, we can simulate and plot the results.

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

Some features we used here:

• The cell volume wasmodeled as a differential equation, which was coupled to the reaction network model.
• The cell divisions were created byincorporating events into the model.
• We designated a specific numericsolver and corresponding solver options.
• The model simulation wasplotted using Plots.jl.

Catalyst developers are active on theJulia Discourse andtheJulia 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 helpsupport it, please star the repository as such metrics may help us secure funding in the future. Ifyou use Catalyst as part of your research, teaching, or other activities, we would be grateful if youcould 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},}