• Added use_jump_ratelaws keyword argument to ode_model, sde_model, hybrid_model, ODEProblem, SDEProblem, and HybridProblem. When set to true, both drift and diffusion terms use the jump/stochastic rate law (binomial propensities) instead of the ODE rate law (power-based). This gives the mathematically correct CLE derived from the CME when species populations are integers. Defaults to false for backward compatibility.

Catalyst 16 is a major release that transitions from ModelingToolkit (MT) v9 to ModelingToolkitBase (the base for ModelingToolkit v11); introduces unified hybrid model support for mixed ODE/SDE/Jump systems; supports user-provided coupled ODEs, SDEs, jump processes, and jump-diffusions; and modernizes the conversion and problem-creation API.

Please also see the ModelingToolkit NEWS.md for all the changes that have occurred in ModelingToolkit as part of v10 and v11, and which are now relevant for Catalyst users.

• Catalyst now depends on and re-exports ModelingToolkitBase instead of ModelingToolkit. This ensures Catalyst remains fully MIT-licensed after the library split that occurred in ModelingToolkit v11. ModelingToolkitBase provides the core symbolic system infrastructure (types, accessors, problem construction) that Catalyst needs. With this update, Catalyst moves from ModelingToolkit v9 to ModelingToolkitBase 1.12+ (part of ModelingToolkit v11).

  Most commonly used functions (unknowns, parameters, equations, @mtkcompile, etc.) are available through ModelingToolkitBase. However, if you relied on structural_simplify, now part of mtkcompile, for reducing models with algebraic equations, you may need to explicitly load ModelingToolkit to obtain equivalent levels of model reduction and optimization:

  using Catalyst
  using ModelingToolkit  # only if you need MTK-specific features not in MTKBase

  The version of mtkcompile in ModelingToolkitBase is less feature filled than in ModelingToolkit, but please be aware that the latter version now loads AGPL-licensed libraries that may impose additional restrictions on your code.

• convert(ODESystem, rs) is replaced by ode_model(rs).

• convert(SDESystem, rs) is replaced by sde_model(rs).

• convert(JumpSystem, rs) is replaced by jump_model(rs).

• convert(NonlinearSystem, rs) is replaced by ss_ode_model(rs).

  The old Base.convert methods have been removed since ModelingToolkit now has one, unified, system type, System. The new functions accept the same keyword arguments as before.

  Before:

  osys = convert(ODESystem, rn; remove_conserved = true)

  After:

  osys = ode_model(rn; remove_conserved = true)

  Problem constructors (ODEProblem, SDEProblem, NonlinearProblem, etc.) still accept a ReactionSystem directly and continue to work as before.

• JumpInputs has been removed. Use JumpProblem(rs, u0, tspan, p; ...) directly instead.

• DiscreteProblem(rs::ReactionSystem, ...) has been removed. Use JumpProblem(rs, u0, tspan, p; ...) directly instead. DiscreteProblem is also no longer exported by Catalyst.

  Before (old API):

  jinputs = JumpInputs(rs, u0, tspan, p)
  jprob = JumpProblem(jinputs)

  After (new API):

  jprob = JumpProblem(rs, u0, tspan, p)

  For advanced usage where you need to reuse a converted system or customize the JumpSystem, use jump_model(rs; ...) to get the System, then call JumpProblem(sys, op, tspan; ...) from ModelingToolkitBase.

• JumpProblem(rs, ...) now produces pure jump systems only (symmetric with ODEProblem and SDEProblem). For hybrid ODE+Jump or SDE+Jump systems, use the new HybridProblem instead.

• jump_model no longer respects ODE/SDE metadata — it forces all reactions to Jump scale, only preserving VariableRateJump metadata if explicitly set. Use hybrid_model or HybridProblem for mixed-scale systems.

• The defaults field and keyword argument have been renamed to initial_conditions throughout Catalyst, following the same change in ModelingToolkitBase. This affects:

  • The ReactionSystem struct field (defaults → initial_conditions)
  • The keyword argument in conversion functions (ode_model, sde_model, jump_model, ss_ode_model, hybrid_model)
  • The accessor function (MT.defaults(rs) → MT.initial_conditions(rs))

• setdefaults! has been removed. Use default values in @species and @parameters declarations, or pass initial_conditions to the ReactionSystem constructor:

  # Setting initial conditions via the constructor
  rn = @reaction_network begin
      @species A(t) = 1.0
      @parameters k = 0.5
      k, A --> B
  end
  
  # Or via the keyword argument
  @named rn = ReactionSystem(rxs, t; initial_conditions = [A => 1.0, k => 0.5])

• default_u0 and default_p keyword arguments have been removed from conversion functions. Use initial_conditions instead:

  # Before:
  osys = convert(ODESystem, rn; default_u0 = [A => 1.0], default_p = [k => 0.5])
  
  # After:
  osys = ode_model(rn; initial_conditions = Dict(A => 1.0, k => 0.5))

• compose, extend, and flatten now only accept ReactionSystem subsystems. Composing a ReactionSystem with a generic MT.System is no longer supported.

  Migration: Convert your System to a ReactionSystem with algebraic/ODE equations:

  # Before (no longer works):
  @named constraints = System([x ~ a], t)
  extended = extend(constraints, rn)
  
  # After:
  @named constraints = ReactionSystem([x ~ a], t)
  extended = extend(constraints, rn)

• New validate_units/assert_valid_units functions defined on Reactions and ReactionSystems for checking units. For Reactions, consistency of terms within the rate expression are checked (i.e. if the rate is a sum of terms), and consistency of substrate and product species units are checked. All other checks are handled by the ReactionSystem.

• validate_units now correctly handles non-SI units (M, μM) in rate expressions. Previously, unit validation could produce spurious failures with non-SI units due to floating-point precision loss during unit expansion. Note that units should be specified using symbolic units (us"...") rather than concrete units (u"..."), as symbolic units use exact arithmetic and avoid precision issues. For example, you can now use μM concentrations:

  using DynamicQuantities
  @independent_variables t [unit=us"s"]
  @species X(t) [unit=us"μM"]
  @parameters k [unit=us"μM/s"]

• validate_units(rs::ReactionSystem) no longer multiplies substrate units into the rate for only_use_rate=true reactions. Previously, custom rate laws (e.g., Michaelis-Menten rates passed via the => arrow) would fail validation because substrate units were erroneously included.

• validate_units(rs::ReactionSystem) now also validates non-reaction equations (including those involving brownian and poissonian noise terms), checking that left- and right-hand side units match. Note, currently events and user-provided jumps are not unit checked.

• Unit checks now validate comparison/conditional expressions and exponent units in equations and rates. Comparisons are treated as unitless predicates with operand compatibility checks, ifelse branches must agree in units, and exponents are required to be unitless.

• Catalyst now provides explicit unit validation APIs validate_units (non-throwing) and assert_valid_units (throwing), rather than overloading ModelingToolkitBase.validate. These APIs and their report/error types are public APIs, but are no longer exported; call them as Catalyst.validate_units(...), Catalyst.assert_valid_units(...), etc. (on Julia 1.11+, this is marked via public).

• New DSL option: @unit_checks true|false. This controls constructor-time unit validation for @reaction_network and @network_component generated systems, matching the unit_checks keyword in programmatic ReactionSystem(...) and Reaction(...) construction (false by default).

• In all interfaces unit checking is opt-in: users with units should pass unit_checks = true to ReactionSystem and Reaction constructors, or use the @unit_checks true option within the DSL to enable validation at construction time.

• Note that ReactionSystems assume Reactions have already performed their associated checks, so for full unit checking pass unit_checks = true to both (or use the DSL @unit_checks true option).

• == no longer performs structural comparison on ReactionSystems. It now falls back to object identity (===), consistent with ModelingToolkitBase's System.

• All instances of structural_simplify (whenever as arguments to ODEProblem or as a separate function) have been replaced with mtkcompile, which serves an identical purpose (with exceptions of the move from ModelingToolkit to ModelingToolkitBase, as described previously).

The following are now considered internal, no longer exported, and could be removed at any time:

• isequivalent.
• symmap_to_varmap

Previously, Catalyst supported so-called LatticeReactionSystem spatial models. These have been renamed DiscreteSpaceReactionSystems. Several corresponding name changes are also affected by this. E.g.

• Catalyst.has_cartesian_lattice is now Catalyst.has_cartesian_dspace.
• Catalyst.lat_getp is not Catalyst.spat_getp.
• Catalyst.lattice_animation is now Catalyst.dspace_animation.

Furthermore, instead of creating spatial JumpProblems through spatial DiscreteProblems:

dprob = DiscreteProblem(lrs, u0, tspan, ps)
jprob = JumpProblem(lrs, dprob, NSM())

these are now created directly:

dprob = DiscreteProblem(lrs, u0, tspan, ps)
jprob = JumpProblem(lrs, u0, tspan, ps)

An aggregator can be designated through the aggregator key word argument (with NSM() being the default).

• New hybrid_model(rs::ReactionSystem; ...) function that converts a ReactionSystem to a unified ModelingToolkitBase.System containing ODE equations, Brownian noise terms (SDE), and/or jump processes, based on each reaction's PhysicalScale metadata.

  The existing ode_model, sde_model, and jump_model are now thin wrappers around hybrid_model. This enables true hybrid models mixing ODE, SDE, and Jump reactions in a single system.

  # Tag reactions with different physical scales
  rn = @reaction_network begin
      k1, S --> P, [physical_scale = PhysicalScale.ODE]
      k2, P --> S, [physical_scale = PhysicalScale.SDE]
      k3, S --> 0, [physical_scale = PhysicalScale.Jump]
  end
  sys = hybrid_model(rn)

  When building reactions programmatically, set the scale via the metadata keyword:

  @species S(t) P(t)
  @parameters k1 k2
  rx1 = Reaction(k1, [S], [P]; metadata = [:physical_scale => PhysicalScale.ODE])
  rx2 = Reaction(k2, [P], [S]; metadata = [:physical_scale => PhysicalScale.Jump])

  Users can also override scales via keyword arguments:

  sys = hybrid_model(rn; default_scale = PhysicalScale.ODE)

• New HybridProblem(rs, u0, tspan, p; ...) function creates problems for systems with per-reaction scale control. Uses hybrid_model internally and respects per-reaction PhysicalScale metadata.

  The return type depends on which reaction scales are present:

  • Pure ODE (only ODE-scale reactions) → ODEProblem
  • Pure SDE or ODE+SDE (no jumps) → SDEProblem
  • Any jumps present (ODE+Jump, SDE+Jump, ODE+SDE+Jump) → JumpProblem

  For SDE+Jump combinations, the returned JumpProblem wraps an SDEProblem internally, allowing simulation of hybrid systems with both diffusive noise and discrete jumps. For ODE+Jump combinations it wraps an ODEProblem.

  rn = @reaction_network begin
      k1, S --> P
      k2, P --> S
  end
  
  # Pure ODE via HybridProblem
  prob_ode = HybridProblem(rn, [:S => 100.0, :P => 0.0], (0.0, 10.0),
      [:k1 => 1.0, :k2 => 0.5]; default_scale = PhysicalScale.ODE)
  prob_ode isa ODEProblem  # true
  
  # Hybrid ODE+Jump
  rn2 = @reaction_network begin
      k1, S --> P, [physical_scale = PhysicalScale.ODE]
      k2, P --> S, [physical_scale = PhysicalScale.Jump]
  end
  prob = HybridProblem(rn2, [:S => 100.0, :P => 0.0], (0.0, 10.0),
      [:k1 => 1.0, :k2 => 0.5])
  prob isa JumpProblem  # true
  sol = solve(prob, Tsit5())
  
  # Hybrid SDE+Jump (requires SDE solver)
  rn3 = @reaction_network begin
      k1, S --> P, [physical_scale = PhysicalScale.SDE]
      k2, P --> S, [physical_scale = PhysicalScale.Jump]
  end
  prob = HybridProblem(rn3, [:S => 100.0, :P => 0.0], (0.0, 10.0),
      [:k1 => 1.0, :k2 => 0.5])
  prob isa JumpProblem  # true
  prob.prob isa SDEProblem  # true - underlying problem is SDE
  sol = solve(prob, SRIW1())  # use SDE solver from StochasticDiffEq

The affect part of events in the DSL (and only in the DSL) now uses => instead of ~. I.e.

@reaction_network begin
    @discrete_events begin
        [1.0] => [X => X + 10]
        [X ~ 1.0] => [X => X + 1]
    end
    (p,d), 0 <--> X
end

instead of

@reaction_network begin
    @discrete_events begin
        [1.0] => [X ~ X + 10]
        [X ~ 1.0] => [X ~=>~ X + 1]
    end
    (p,d), 0 <--> X
end

This is to align with a future expected ModelingToolkit update creating so-called "Assignment affects". Other aspects of events within the DSL remain unchanged (see next section regarding more general event updates).

Catalyst now supports several ModelingToolkit v11 features for events and callbacks. See the ModelingToolkit documentation for full details on these features.

• Pre(x) for pre-event values. Event affects that reference the current value of a variable or parameter must now use Pre(x) to access the value before the event fires. In the DSL (@continuous_events, @discrete_events), Pre is automatically wrapped — no user action needed. When constructing callbacks programmatically via SymbolicContinuousCallback or SymbolicDiscreteCallback, you must explicitly write Pre(x):

  # DSL (auto-wrapped, works as before):
  rn = @reaction_network begin
      @discrete_events 2.0 => [p => p + 0.1]   # becomes p ~ Pre(p + 0.1)
      (p, d), 0 <--> X
  end
  
  # Programmatic (explicit Pre required):
  continuous_events = [X ~ 1.0] => [X ~ Pre(X) - 0.5]

• @discretes for event-modified parameters. Parameters that are changed by events must now be declared as discrete variables using @discretes in the DSL (or @discretes from ModelingToolkitBase programmatically). In the DSL, parameters that appear on the LHS of event affects are automatically inferred as discretes — no explicit @discretes declaration is needed unless you want to set default values. Note that discretes should be created as time-dependent parameters. When constructing callbacks programmatically, you must pass the discrete_parameters keyword:

  # DSL (p auto-inferred as discrete):
  rn = @reaction_network begin
      @discrete_events 12 => [p => (p + 1) % 2]
      (p, d), 0 <--> X
  end
  
  # Programmatic (must declare and pass discrete_parameters):
  @discretes p(t) = 1.0
  discrete_events = SymbolicDiscreteCallback(
      12 => [p ~ Pre(p + 1) % 2];
      discrete_parameters = [p]
  )

• @discretes DSL option. A new @discretes option can be used in @reaction_network and @network_component to explicitly declare discrete variables with default values:

  rn = @reaction_network begin
      @discretes α(t)=5.0 β(t)=1.0
      @species V(t) = 0.0
      @continuous_events [V ~ 2.5] => [α => 0, β => 0]
      α, 0 --> V
      β, V --> 0
  end

• ReactionSystem now has brownians, poissonians, and jumps fields to support non-reaction differential equations containing Brownian motions and/or unit Poisson counting process-based jumps, and to also support explicitly passed coupled jumps (i.e. symbolic VariableRateJumps, ConstantRateJumps, or MassActionJumps). This mirrors the capability in ModelingToolkitBase's System struct and enables modeling hybrid stochastic systems where reactions couple with other equation types (ODEs, SDEs, PDMPs, jump-diffusions, etc).

  using Catalyst
  t = default_t()
  D = default_time_deriv()
  @brownians B
  @variables V(t)
  @species S(t) P(t)
  @parameters k λ
  
  # brownians is the 5th positional argument, jumps is a keyword arg
  @named sys = ReactionSystem([Reaction(k, [S], [P]), D(V) ~ -V + B],
      t, [V, S, P], [k, λ], [B])

  The two-argument constructor auto-discovers brownians and poissonians from equations:

  @brownians B
  @variables V(t)
  @parameters λ
  @named sys = ReactionSystem([D(V) ~ λ*V + B], t)

• User-provided brownians/jumps are merged with reaction-generated ones at conversion time in hybrid_model, sde_model, and HybridProblem.

• Conversion functions enforce scale constraints:

  • ode_model: Errors if brownians OR jumps are present (pure ODE only)
  • sde_model: Errors if jumps are present (ODE + SDE constraints allowed)
  • jump_model: Errors if brownians OR non-reaction equations are present
  • Use hybrid_model for mixed systems

• A new @brownians option can be used within the @reaction_network and @network_component macros to declare Brownian variables:

  rn = @reaction_network begin
      @species S(t) P(t)
      @parameters k
      @brownians B W
      k, S --> P
  end

  Brownian variables declared via @brownians are automatically included in the ReactionSystem's brownians field.

• A new @poissonians option can be used in @reaction_network and @network_component to declare Poissonian variables — symbolic Poisson counting processes with associated rates. Poissonians represent discrete jump noise and are the jump-process analogue of Brownians:

  rn = @reaction_network begin
      @parameters λ k d
      @variables X(t)
      @poissonians dN(λ)
      @equations D(X) ~ dN
      (k, d), 0 <--> S
  end

  Here dN(λ) represents a Poisson counting process with rate/intensity λ.

  Multiple poissonians and begin/end block syntax are supported:

  rn = @reaction_network begin
      @parameters λ₁ λ₂
      @poissonians begin
          dN₁(λ₁)
          dN₂(λ₂)
      end
      @equations begin
          D(V) ~ dN₁
          D(W) ~ dN₂
      end
      (p, d), 0 <--> X
  end

  Symbols appearing in poissonian rates must be pre-declared (via @parameters, @species, etc.) or interpolated — they are not auto-inferred from the rate expression.

• Poissonians are handled through composition (flatten, extend, compose) with union semantics and through equivalence checking (isequivalent), following the same patterns as brownians.

• HybridProblem automatically converts poissonians to ConstantRateJump (when the rate depends only on parameters) or VariableRateJump (when the rate depends on state variables) via mtkcompile:

  using Catalyst, JumpProcesses
  t = default_t()
  D = default_time_deriv()
  @parameters λ k d
  @species S(t)
  @variables X(t)
  @poissonians dN(λ)
  
  @named rn = ReactionSystem([Reaction(k, nothing, [S]), Reaction(d, [S], nothing), 
      D(X) ~ dN], t, [S, X], [λ, k, d]; poissonians = [dN])
  rn = complete(rn)
  
  # Simple Poisson counter: X increases by 1 at rate λ
  prob = HybridProblem(rn, [:S => 10.0, :X => 0.0], (0.0, 10.0),
      [:λ => 5.0, :k => 1.0, :d => 0.1])
  sol = solve(prob, Tsit5())
  
  # Jump-diffusion: combine Brownian noise and Poisson jumps
  @brownians dW
  @poissonians dN(λ)
  @parameters σ
  @named jd = ReactionSystem([D(X) ~ σ * dW + dN], t, [X], [λ, σ], [dW]; 
      poissonians = [dN])
  jd = complete(jd)
  prob = HybridProblem(jd, [:X => 0.0], (0.0, 10.0), [:λ => 3.0, :σ => 1.0])
  sol = solve(prob, SRIW1())  # SDE solver for the diffusion component
  
  # State-dependent rate → VariableRateJump (e.g. pure-birth / Yule process)
  @poissonians dN_birth(k * X)  # rate depends on state X

  We hope to add a more refined classification into VariableRateJumps, ConstantRateJumps, and MassActionJumps in the future.

• ode_model, sde_model, and jump_model error with informative messages when called on systems containing poissonians. Use HybridProblem instead, which handles the poissonian-to-jump conversion automatically.

• sde_model and SDEProblem now accept use_legacy_noise = true (default) which uses the traditional noise_eqs matrix approach for simple SDE systems without constraints, avoiding the overhead of mtkcompile. Set to false to use the new Brownian-variable-based approach where SDE noise terms are embedded directly in the equation RHS as stoich * sqrt(|ratelaw|) * B_j.

• SDEProblem(rs::ReactionSystem, ...) calls mtkcompile automatically when necessary — including when the system has user-declared brownians, when using Brownian-based noise (use_legacy_noise = false), or when the system has constraints/algebraic equations. Users do not need to pass structural_simplify = true in these cases.

• Parameters that appear in stoichiometric positions within the @reaction_network DSL are now automatically designated as Int64. There is no need to manually declare them with a type annotation:

  rn = @reaction_network begin
      k, n*A --> B    # n is inferred as a parameter and typed as Int64
  end

  If you need a non-integer type for a stoichiometric parameter (e.g., for fractional stoichiometry), explicitly declare it with the desired type — this overrides the automatic Int64 designation:

  rn = @reaction_network begin
      @parameters n::Float64
      k, n*A --> B    # n is Float64 because of the explicit declaration
  end

New metadata keys U0Map and ParameterMap with convenience accessors (has_u0_map, get_u0_map, set_u0_map and has_parameter_map, get_parameter_map, set_parameter_map) allow file parsers to store species/variable and parameter value mappings on a ReactionSystem in formats distinct from initial_conditions. These are set via the metadata keyword or set_* functions and are preserved through flatten, complete, and model conversions. Note that Catalyst does not use these at all, they, and other metadata, are intended to provide a way for users to cache additional information in a system.

Previously there were some limitations when applying remake on a NonlinearProblems where remove_conserved = true had been used. This is no longer the case. The associated warning has been removed, and the default is now remove_conserved = false for NonlinearProblems, consistent with all other problem types.

• ReactionSystem now supports a tstops keyword argument for specifying extra time points at which the integrator should stop. These can be numeric values or symbolic expressions of parameters. Tstops are automatically forwarded through model conversions (ode_model, sde_model, jump_model, hybrid_model) to the underlying System.

  @parameters t_switch
  @species X(t)
  @named rn = ReactionSystem([Reaction(k, nothing, [X])], t;
      tstops = [t_switch, 5.0])

• New @tstops DSL option for declaring tstops within @reaction_network and @network_component. Symbols appearing in tstop expressions that are not already declared are auto-discovered as parameters (consistent with how reaction rate parameters are inferred):

  rn = @reaction_network begin
      @tstops begin
          t_switch           # auto-discovered as a parameter
          0.5 * t_switch
      end
      (p, d), 0 <--> X
  end

• Tstops are validated at construction time — expressions containing unknowns (species/variables) or the independent variable are rejected with an informative error, since tstops must be computable before the solve starts.

• Tstops are handled through composition (flatten, extend) and equivalence checking (isequivalent), following the same patterns as events and brownians.

• Solver support: Symbolic tstops are supported for all problem types: ODEProblem, SDEProblem, JumpProblem, and HybridProblem.

• The Catalyst release process is changing; certain core dependencies of Catalyst will now be capped to ensure Catalyst releases are only installed with versions of dependencies for which Catalyst CI and doc build tests pass (at the time the release is made). If you need a dependency version increased, please open an issue and we can update it and make a new Catalyst release once testing against the newer dependency version is complete.
• BREAKING: New formula for inferring variables from equations (declared using the @equations options) in the DSL. The order of inference of species/variables/parameters is now:
  1. Every symbol explicitly declared using @species, @variables, and @parameters are assigned to the correct category.
  2. Every symbol used as a reaction reactant is inferred as a species.
  3. Every symbol not declared in (1) or (2) that occurs in an expression provided after @equations is inferred as a variable.
  4. Every symbol not declared in (1), (2), or (3) that occurs either as a reaction rate or stoichiometric coefficient is inferred to be a parameter. E.g. in

     @reaction_network begin
         @equations V1 + S ~ V2^2
         (p + S + V1), S --> 0
     end

     S is inferred as a species, V1 and V2 as variables, and p as a parameter. The previous special cases for the @observables, @compounds, and @differentials options still hold. Finally, the @require_declaration options (described in more detail below) can now be used to require everything to be explicitly declared.
• BREAKING: New formula for determining whether the default differentials have been used within an @equations option. Now, if any expression D(...) is encountered (where ... can be anything), this is inferred as usage of the default differential D. E.g. in the following equations D is inferred as a differential with respect to the default independent variable:

  @reaction_network begin
    @equations D(V) + V ~ 1
  end
  @reaction_network begin
    @equations D(D(V)) ~ 1
  end

  Please note that this cannot be used at the same time as D is used to represent a species, variable, or parameter (including if these are implicitly designated as such by e.g. appearing as a reaction reactant).
• BREAKING: Array symbolics support is more consistent with ModelingToolkit v9. Parameter arrays are no longer scalarized by Catalyst, while species and variables arrays still are (as in ModelingToolkit). As such, parameter arrays should now be specified as arrays in value mappings, i.e.

  @parameters k[1:4]
  pmap = [k => rand(4)]

  While one can still manually scalarize a parameter array, it is recommended not to do this as it has significant performance costs with ModelingToolkit v9. Note, scalarized parameter arrays passed to the two-argument ReactionSystem constructor may become unscalarized.
• BREAKING: We have introduced a restriction on bundling of reactions in the DSL. Now, bundling is not permitted if multiple rates are provided but only one set each of substrates/products. E.g. this model:

  @reaction_network begin
    (k1,k2), X --> Y
  end

  will now throw an error. The reason that users attempting to write bi-directional reactions but typing --> instead of <--> would get a wrong model. We decided that this kind of bundling was unlikely to be used, and throwing errors for people who made the typo was more important. If you use this type of bundling and it indeed is useful to you, please raise and issue and we will see if we can sort something out.
• BREAKING: Catalyst's network visualization capability has shifted from using Graphviz to GraphMakie.jl. To use this functionality, load the GraphMakie extension by installing Catalyst and GraphMakie, along with a Makie backend like CairoMakie or GLMakie. There are two new methods for visualizing graphs: plot_network and plot_complexes, which respectively display the species-reaction graph and complex graph.

  using Catalyst, GraphMakie, GLMakie
  brusselator = @reaction_network begin
     A, ∅ --> X
     1, 2X + Y --> 3X
     B, X --> Y
     1, X --> ∅
  end
  plot_network(brusselator)

• BREAKING: The letter Ø (used in Danish/Norwegian alphabet) is now considered the same as ∅ (empty set). It can no longer be used as a species/parameter.
• BREAKING: When converting a Catalyst ReactionSystem to a ModelingToolkit system, for example an ODESystem, Catalyst defined functions like hill(A,B,C,D) are now replaced with the explicit rational function they represent in the equations of the generated system. For example mm(X,v,K) will be replaced with v*X / (X + K). This can be disabled by passing the keyword argument expand_catalyst_funs = false. e.g.

  using Catalyst
  rn = @reaction_network begin
    hill(X,v,K,n), A --> 0
  end
  osys = convert(ODESystem, rn)

  generates an ODE system with D(A) ~ -((v*A(t)*(X^n)) / (K^n + X^n)), while

  osys = convert(ODESystem, rn; expand_catalyst_funs = false)

  generates an ODE system with D(A) ~ -A(t)*hill(X, v, K, n). This keyword argument can also be passed to problems defined over ReactionSystems, i.e. when calling ODEProblem(rn, u0, tspan, p; expand_catalyst_funs = false).
• It is no longer recommended to install and use the full OrdinaryDiffEq library to access specific ODE solvers. Instead, only install the specific OrdinaryDiffEq sub-libraries that contain the desired solver. This significantly reduces installation and package loading times. I.e. to use the default solver that auto-switches between explicit and implicit methods, install OrdinaryDiffEqDefault. To use Tsit5 install OrdinaryDiffEqTsit5, etc. The possible sub-libraries, each containing different solvers, can be viewed here.
• It should now be safe to use remake on problems which have had conservation laws removed with the exception of NonlinearProblems or NonlinearSystems. For NonlinearProblems it is safe to use remake if only updating u0 values, but it is not safe to update the value of the conserved constant, Γ. See the FAQ for details.
• Functional (e.g. time-dependent) parameters can now be used in Catalyst models. These can e.g. be used to incorporate arbitrary time-dependent functions (as a parameter) in a model. For more details on how to use these, please read: https://docs.sciml.ai/Catalyst/stable/model_creation/functional_parameters/.
• Scoped species/variables/parameters are now treated similar to the latest MTK releases (≥ 9.49).
• A tutorial on making interactive plot displays using Makie has been added.
• The BifurcationKit extension has been updated to v.4.
• There is a new DSL option @require_declaration that will turn off automatic inferring for species, parameters, and variables in the DSL. For example, the following will now error:

  rn = @reaction_network begin
    @require_declaration
    (k1, k2), A <--> B
  end

  When this flag is set, all symbolics must be explicitly declared.

  rn = @reaction_network begin
    @species A(t) B(t)
    @parameters k1 k2
    @require_declaration
    (k1, k2), A <--> B
  end

• Support for user-defined functions on the RHS when providing coupled equations for CRNs using the @equations macro. For example, the following now works:

  using Catalyst
  f(A, t) = 2*A*t
  rn = @reaction_network begin
    @equations D(A) ~ f(A,t)
  end

  Note that user-defined functions will not work on the LHS of equations.

• Symbolics 6 support.

• Support for simulating stochastic chemical kinetics models with explicitly time-dependent propensities (i.e. where the resulting JumpSystem contains VariableRateJumps). As such JumpProblems need to be defined over ODEProblems or SDEProblems instead of DiscreteProblems we have introduced a new input struct, JumpInputs, that handles selecting via analysis of the generated JumpSystem, i.e. one can now say

  using Catalyst, OrdinaryDiffEq, JumpProcesses, Plots
  rn = @reaction_network begin
      k*(1 + sin(t)), 0 --> A
  end
  jinput = JumpInputs(rn, [:A => 0], (0.0, 10.0), [:k => .5])
  # note that jinput.prob isa ODEProblem
  jprob = JumpProblem(jinput)
  sol = solve(jprob, Tsit5())
  plot(sol, idxs = :A)
  
  rn = @reaction_network begin
      k, 0 --> A
  end
  jinput = JumpInputs(rn, [:A => 0], (0.0, 10.0), [:k => .5])
  # note that jinput.prob isa DiscreteProblem
  jprob = JumpProblem(jinput)
  sol = solve(jprob)
  plot(sol, idxs = :A)

  When calling solve for problems with explicit time-dependent propensities, i.e. where jinput.prob isa ODEProblem, note that one must currently explicitly select an ODE solver to handle time-stepping and integrating the time-dependent propensities.
• Note that solutions to jump problems with explicit time-dependent propensities, i.e. a JumpProblem over an ODEProblem, require manual selection of the variables to plot. That is, currently plot(sol) will error in this case due to limitations in the SciMLBase plot recipe.

• Support for auto-algorithm selection in JumpProblems. For systems with only propensities that do not have an explicit time-dependence (i.e. that are not VariableRateJumps in JumpProcesses), one can now run model simulations via

  using Catalyst, JumpProcesses
  model = @reaction_network begin
    kB, S + E --> SE
    kD, SE --> S + E
    kP, SE --> P + E
  end
  u0 = [:S => 50, :E => 10, :SE => 0, :P => 0]
  tspan = (0., 200.)
  ps = [:kB => 0.01, :kD => 0.1, :kP => 0.1]
  dprob = DiscreteProblem(model, u0, tspan, ps)
  jprob = JumpProblem(model, dprob)
  sol = solve(jprob)

  For small systems this will just use Gillespie's Direct method, transitioning to using RSSA and RSSACR as system size increase. Once can still manually select a given SSA, but no longer needs to specify SSAStepper when calling solve, i.e.

  # use the SortingDirect method instead
  jprob = JumpProblem(model, dprob, SortingDirect())
  sol = solve(jprob)

• Latexify recipe improvements including display fixes for array symbolics.
• Deficiency one and concentration robustness checks.

The expansion of ReactionSystem models to spatial lattices has been enabled. Here follows a simple example where a Brusselator model is expanded to a 20x20 grid of compartments, with diffusion for species X, and then simulated using ODEs. Finally, an animation of the simulation is created.

using Catalyst, CairoMakie, OrdinaryDiffEq

# Create `LatticeReactionSystem` model.
brusselator = @reaction_network begin
    A, ∅ --> X
    1, 2X + Y --> 3X
    B, X --> Y
    1, X --> ∅
end
diffusion_rx = @transport_reaction D X
lattice = CartesianGrid((20,20))
lrs = LatticeReactionSystem(brusselator, [diffusion_rx], lattice)

# Create a spatial `ODEProblem`.
u0 = [:X => rand(20, 20), :Y => 10.0]
tspan = (0.0, 40.0)
ps = [:A => 1.0, :B => 4.0, :D => 0.2]
oprob = ODEProblem(lrs, u0, tspan, ps)

# Simulate the ODE and plot the results.
sol = solve(oprob, FBDF())
dspace_animation(sol, :X, lrs, "brusselator.mp4")

The addition of spatial modelling in Catalyst contains a large number of new features, all of which are described in the corresponding documentation.

Bug fix to address that independent variables, like time, should now be @parameters according to MTKv9. Converted internal time variables to consistently use default_t() to hopefully avoid such issues going forward.

Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which introduced several breaking changes to Catalyst. A summary of these (and how to handle them) can be found here. These are briefly summarised in the following bullet points:

• ReactionSystems must now be marked complete before they are exposed to most forms of simulation and analysis. With the exception of ReactionSystems created through the @reaction_network macro, all ReactionSystems are not marked complete upon construction. The complete function can be used to mark ReactionSystems as complete. To construct a ReactionSystem that is not marked complete via the DSL the new @network_component macro can be used.
• The states function has been replaced with unknowns. The get_states function has been replaced with get_unknowns.
• Support for most units (with the exception of s, m, kg, A, K, mol, and cd) has currently been dropped by ModelingToolkit, and hence they are unavailable via Catalyst too. Its is expected that eventually support for relevant chemical units such as molar will return to ModelingToolkit (and should then immediately work in Catalyst too).
• Problem parameter values are now accessed through prob.ps[p] (rather than prob[p]).
• ModelingToolkit currently does not support the safe application of the remake function, or safe direct mutation, for problems for which remove_conserved = true was used when updating the values of initial conditions. Instead, the values of each conserved constant must be directly specified.
• The reactionparams, numreactionparams, and reactionparamsmap functions have been deprecated and removed.
• To be more consistent with ModelingToolkit's immutability requirement for systems, we have removed API functions that mutate ReactionSystems such as addparam!, addreaction!, addspecies, @add_reactions, and merge!. Please use ModelingToolkit.extend and ModelingToolkit.compose to generate new merged and/or composed ReactionSystems from multiple component systems.

• default_t() and default_time_deriv() functions should be used for creating the default time independent variable and its differential. i.e.

  # do
  t = default_t()
  @species A(t)
  
  # avoid
  @variables t
  @species A(t)

• It is now possible to add metadata to individual reactions, e.g. using:

  rn = @reaction_network begin
      @parameters η
      k, 2X --> X2, [description="Dimerisation"]
  end
  getdescription(rn)

  a more detailed description can be found here.
• SDEProblem no longer takes the noise_scaling argument. Noise scaling is now handled through the noise_scaling metadata (described in more detail here)
• Fields of the internal Reaction structure have been changed. ReactionSystemss saved using serialize on previous Catalyst versions cannot be loaded using this (or later) versions.
• A new function, save_reactionsystem, which permits the writing of ReactionSystem models to files, has been created. A thorough description of this function can be found here
• Updated how compounds are created. E.g. use

  @variables t C(t) O(t)
  @compound CO2 ~ C + 2O

  to create a compound species CO2 that consists of C and two O.
• Added documentation for chemistry-related functionality (compound creation and reaction balancing).
• Added function isautonomous to check if a ReactionSystem is autonomous.
• Added function steady_state_stability to compute stability for steady states. Example:

  # Creates model.
  rn = @reaction_network begin
      (p,d), 0 <--> X
  end
  p = [:p => 1.0, :d => 0.5]
  
  # Finds (the trivial) steady state, and computes stability.
  steady_state = [2.0]
  steady_state_stability(steady_state, rn, p)

  Here, steady_state_stability takes an optional keyword argument tol = 10*sqrt(eps()), which is used to check that the real part of all eigenvalues are at least tol away from zero. Eigenvalues within tol of zero indicate that stability may not be reliably calculated.
• Added a DSL option, @combinatoric_ratelaws, which can be used to toggle whether to use combinatorial rate laws within the DSL (this feature was already supported for programmatic modelling). Example:

  # Creates model.
  rn = @reaction_network begin
      @combinatoric_ratelaws false
      (kB,kD), 2X <--> X2
  end

• Added a DSL option, @observables for creating observables (this feature was already supported for programmatic modelling).
• Added DSL options @continuous_events and @discrete_events to add events to a model as part of its creation (this feature was already supported for programmatic modelling). Example:

  rn = @reaction_network begin
      @continuous_events begin
          [X ~ 1.0] => [X ~ X + 1.0]
      end
      d, X --> 0
  end

• Added DSL option @equations to add (algebraic or differential) equations to a model as part of its creation (this feature was already supported for programmatic modelling). Example:

  rn = @reaction_network begin
      @equations begin
          D(V) ~ 1 - V
      end
      (p/V,d/V), 0 <--> X
  end

  couples the ODE $dV/dt = 1 - V$ to the reaction system.
• Coupled reaction networks and differential equation (or algebraic differential equation) systems can now be converted to SDESystems and NonlinearSystems.

• Added CatalystStructuralIdentifiabilityExtension, which permits StructuralIdentifiability.jl to be applied directly to Catalyst systems. E.g. use

  using Catalyst, StructuralIdentifiability
  goodwind_oscillator = @reaction_network begin
      (mmr(P,pₘ,1), dₘ), 0 <--> M
      (pₑ*M,dₑ), 0 <--> E
      (pₚ*E,dₚ), 0 <--> P
  end
  assess_identifiability(goodwind_oscillator; measured_quantities=[:M])

  to assess (global) structural identifiability for all parameters and variables of the goodwind_oscillator model (under the presumption that we can measure M only).
• Automatically handles conservation laws for structural identifiability problems (eliminates these internally to speed up computations).
• A more detailed of how this extension works can be found here.

• Add a CatalystBifurcationKitExtension, permitting BifurcationKit's BifurcationProblems to be created from Catalyst reaction networks. Example usage:

  using Catalyst
  wilhelm_2009_model = @reaction_network begin
      k1, Y --> 2X
      k2, 2X --> X + Y
      k3, X + Y --> Y
      k4, X --> 0
      k5, 0 --> X
  end
  
  using BifurcationKit
  bif_par = :k1
  u_guess = [:X => 5.0, :Y => 2.0]
  p_start = [:k1 => 4.0, :k2 => 1.0, :k3 => 1.0, :k4 => 1.5, :k5 => 1.25]
  plot_var = :X
  bprob = BifurcationProblem(wilhelm_2009_model, u_guess, p_start, bif_par; plot_var = plot_var)
  
  p_span = (2.0, 20.0)
  opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
  
  bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
  
  using Plots
  plot(bif_dia; xguide = "k1", guide = "X")

• Automatically handles elimination of conservation laws for computing bifurcation diagrams.
• Updated Bifurcation documentation with respect to this new feature.

• Added a CatalystHomotopyContinuationExtension extension, which exports the hc_steady_state function if HomotopyContinuation is exported. hc_steady_state finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example:

wilhelm_2009_model = @reaction_network begin
    k1, Y --> 2X
    k2, 2X --> X + Y
    k3, X + Y --> Y
    k4, X --> 0
end
ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
hc_steady_states(wilhelm_2009_model, ps)

• Added the ability to create species that represent chemical compounds and know their constituents. For example, water can be created and queried as

  @variables t
  @species H(t) O(t)
  @compound H2O(t) 2*H O
  iscompound(H2O) == true
  isspecies(H2O) == true   # compounds are also species, so can be used in reactions
  isequal(components(H2O), [H, O])
  coefficients(H2O) == [2, 1]

• Added reaction balancing via the balance_reaction command, which returns a vector of balanced reaction versions, i.e.

  @variables t
  @species H(t) O(t) C(t)
  @compound CH4(t) C 4H
  @compound O2(t) 2O
  @compound CO2(t) C 2O
  @compound H2O(t) 2H O
  
  # unbalanced reaction to balance
  rx = Reaction(1.0, [CH4, O2], [CO2, H2O])
  
  # calculate a balanced version, this returns a vector
  # storing a single balanced version of the reaction in this case
  brxs = balance_reaction(rx)
  
  # what one would calculate by hand
  balanced_rx = Reaction(1.0, [CH4, O2], [CO2, H2O], [1, 2], [1, 2])
  
  # testing equality
  @test isequal(balanced_rx, first(brxs))

• Note that balancing works via calculating the nullspace of an associated integer matrix that stores in entry (i,j) a signed integer representing the number of times the i'th atom appears within the jth compound. The entry is positive for a substrate and negative for a product. One cannot balance a reaction involving compounds of compounds currently. A non-empty solution vector is returned if the reaction can be balanced in exactly one way with minimal coefficients while preserving the set of substrates and products, i.e. if the dimension of the nullspace is one. If the dimension is greater than one we return a Reaction for each nullspace basis vector, but note that they may currently interchange substrates and products (i.e. we do not solve for if there is a linear combination of them that preserves the set of substrates and products). An empty Reaction vector indicates it is not possible to balance the reaction.

• Array parameters, species, and variables can be use in the DSL if explicitly declared with @parameters, @species, or @variables respectively, i.e.

  rn = @reaction_network begin
      @parameters k[1:2] a
      @variables (V(t))[1:2] W(t)
      @species (X(t))[1:2] Y(t)
      k[1]*a+k[2], X[1] + V[1]*X[2] --> V[2]*W*Y + B*C
  end

• Non-species states can be declared in the DSL using @variables, and custom independent variables (instead of just t) using @ivs. For the latter, the first independent variable is always interpreted as the time variable, and all discovered species are created to be functions of all the ivs. For example in

  rn = @reaction_network begin
      @ivs s x
      @variables A(s) B(x) C(s,x)
      @species D(s) E(x) F(s,x)
      k*C, A*D + B*E --> F + H
  end

  s will be the time variable, H = H(s,x) will be made a function of s and x, and A(s), B(x), and C(s,x) will be non-species state variables.
• Catalyst.isequal_ignore_names has been deprecated for isequivalent(rn1, rn2) to test equality of two networks and ignore their name. To include names in the equality check continue to use rn1 == rn2 or use isequivalent(rn1, rn2; ignorenames = false).

• BREAKING: Parameters should no longer be listed at the end of the DSL macro, but are instead inferred from their position in the reaction statements or via explicit declarations in the DSL macro. By default, any symbol that appears as a substrate or product is a species, while any other is a parameter. That is, parameters are those that only appear within a rate expression and/or as a stoichiometric coefficient. E.g. what previously was

  using Catalyst
  rn = @reaction_network begin
    p, 0 --> X
    d, X --> 0
  end p d

  is now

  using Catalyst
  rn = @reaction_network begin
    p, 0 --> X
    d, X --> 0
  end

  More generally, in the reaction system

  rn = @reaction_network begin
      k*k1*A, A --> B
      k2, k1 + k*A --> B
  end

  k and k2 are inferred as parameters by the preceding convention, while A, B and k1 are species.

• Explicit control over which symbols are treated as parameters vs. species is available through the new DSL macros, @species and @parameters. These can be used to designate when something should be a species or parameter, overriding the default DSL assignments. This allows setting that a symbol which would by default be interpreted as a parameter should actually be a species (or vice-versa). E.g. in:

  using Catalyst
  rn = @reaction_network begin
    @species X(t)
    k*X, 0 --> Y
  end

  X and Y will be considered species, while k will be considered a parameter. These options take the same arguments as standalone the @species (i.e. ModelingToolkit.@variables) and ModelingToolkit.@parameters macros, and support default values and setting metadata. E.g you can set default values using:

  using Catalyst
  rn = @reaction_network begin
    @species X(t)=1.0
    @parameters p=1.0 d=0.1
    p, 0 --> X
    d, X --> 0
  end

  or designate a parameter as representing a constant species using metadata:

  using Catalyst
  rn = @reaction_network begin
    @parameters Y [isconstantspecies=true]
    k, X + Y --> 0
  end

• BREAKING: A standalone @species macro was added and should be used in place of @variables when declaring symbolic chemical species, i.e.

  @parameters k
  @variables t
  @species A(t) B(t)
  rx = Reaction(k, [A], [B])
  @named rs = ReactionSystem([rx], t)

  This will no longer work as substrates and products must be species

  @parameters k
  @variables t A(t) B(t)
  rx = Reaction(k, [A], [B]) # errors as neither A or B are species
  rx = Reaction(k, [A], nothing) # errors as A is not a species
  rx = Reaction(k, nothing, [B]) # errors as B is not a species
  
  # this works as the rate or stoichiometry can be non-species
  @species C(t) D(t)
  rx = Reaction(k*A, [C], [D], [2], [B])
  @named rs = ReactionSystem([rx], t)

  @variables is now reserved for non-chemical species state variables (for example, arising from constraint equations). Internally, species are normal symbolic variables, but with added metadata to indicate they represent chemical species.

• To check if a symbolic variable is a species one can use isspecies:

  @variables t
  @species A(t)
  @variables B(t)
  isspecies(A) == true
  isspecies(B) == false

• BREAKING: Constraint subsystems and the associated keyword argument to ReactionSystem have been removed. Instead, one can simply add ODE or algebraic equations into the list of Reactions passed to a ReactionSystem. i.e. this should now work

  @parameters k α
  @variables t V(t)
  @species A(t)
  rx = Reaction(k*V, nothing, [A])
  D = Differential(t)
  eq = D(V) ~ α
  @named rs = ReactionSystem([rx, eq], t)
  osys = convert(ODESystem, rs)

  which gives the ODE model

  julia> equations(osys)
  2-element Vector{Equation}:
    Differential(t)(A(t)) ~ k*V(t)
    Differential(t)(V(t)) ~ α
  

  Mixing ODEs and algebraic equations is allowed and should work when converting to an ODESystem or NonlinearSystem (if only algebraic equations are included), but is not currently supported when converting to JumpSystems or SDESystems.

• API functions applied to a ReactionSystem, rs, now have:

  • species(rs) give the chemical species of a system.
  • states(rs) give all the variables, both chemical species and non-chemical species of a system.

  Catalyst now orders species before non-species in states(rs) such that states(rs)[1:length(species(rs))] and species(rs) should be the same. Similarly:

  • equations(rs) gives the set of Reactions and Equations of a system.
  • reactions(rs) gives the Reactions of a system.

  As with species, Reactions are always ordered before Equations so that equations(rs)[1:length(reactions(rs))] should be the same ordered list of Reactions as given by reactions(rs).

• Catalyst has been updated for Symbolics v5, and requires Symbolics v5.0.3 or greater and ModelingToolkit v8.47.0 or greater.

• The accessors for a given system, rs, that return the internal arrays at the top-level (i.e. ignoring sub-systems) now have

  • ModelingToolkit.get_states(rs) to get the list of all species and non-species variables.
  • Catalyst.get_species(rs) to get the list of all species variables. Note that get_states(rs)[1:length(get_species(rs))] should be the same ordered list of species as get_species(rs).
  • ModelingToolkit.get_eqs(rs) gives the list of all Reactions and then Equations in the system.
  • Catalyst.get_rxs(rs) gives the list of all Reactions, such that get_eqs(rs)[1:length(get_rx(rs))] is the same ordered list of Reactions as returned by get_rxs(rs).

• BREAKING: Chemical species specified or inferred via the DSL are now created via the same mechanism as @species, and therefore have the associated metadata that is missing from a normal symbolic variable.

• Deprecated functions params and merge have been removed.

• BREAKING: The old notation for the constants representing conserved quantities, _Conlaw, has been replaced with uppercase unicode gamma, "Γ". This can be entered in notebooks, the REPL, or many editors by typing the corresponding Latex command, "\Gamma", and hitting tab. This leads to much cleaner equations when Latexifying systems where conservation laws have been applied. The underlying symbol can also be accessed via Catalyst.CONSERVED_CONSTANT_SYMBOL.

• Modelingtoolkit symbolic continuous and discrete events are now supported when creating ReactionSystems via the continuous_events and discrete_events keyword arguments. As in ModelingToolkit, species, states, and parameters that appear only within events are not detected automatically, and hence the four-argument ReactionSystem constructor, where states and parameters are explicitly passed, must be used unless every variable, state, or parameter in the events appears within a Reaction or Equation too. See the ModelingToolkit docs for more information on using events. Note that JumpSystems only support discrete events at this time.

• Support for states/species that are functions of multiple variables. This enables (symbolically) building PDEs to solve with MethodOfLines. To use multiple independent variables one can say:

  using Catalyst
  using ModelingToolkit: scalarize
  @parameters k[1:7]
  @variables t x y U(x,y,t) V(x,y,t) W(x,y,t)
  rxs = [Reaction(k[1], [U, W], [V, W]),
        Reaction(k[2], [V], [W], [2], [1]),
        Reaction(k[3], [W], [V], [1], [2]),
        Reaction(k[4], [U], nothing),
        Reaction(k[5], nothing, [U]),
        Reaction(k[6], [V], nothing),
        Reaction(k[7], nothing, [V])]
  pars = scalarize(k)
  @named rn = ReactionSystem(rxs, t, [U, V, W], pars; spatial_ivs = [x, y])

  The spatial_ivs keyword lets Catalyst know which independent variables correspond to spatial variables. Note that rate expressions can depend on x and y too, i.e. k[1] * x + y*t would be valid. See the work in progress PDE tutorial to solve the resulting system and add spatial transport.

• API functions to generate substrate, product, and net stoichiometry matrices should now work with floating point stoichiometric coefficients. Note, symbolic coefficients are still not supported by such functions.

• BREAKING: Modified how constant and boundary condition species (in the SBML sense) work. Constant species should now be specified as ModelingToolkit @parameters with the isconstantspecies=true metadata, while non-constant boundary condition species should be specified as ModelingToolkit @variables with the isbcspecies=true metadata. As before, boundary condition species are treated as constant with respect to reactions, but since they are considered variables their dynamics should be defined in a constraint system. Moreover, it is required that BC species appear in a balanced manner (i.e. in each reaction for which a BC species is a reactant it must appear as a substrate and a product with the same stoichiometry). Right now only conversion of ReactionSystems to an ODESystem with a constraint ODESystem or NonlinearSystem, or conversion to a NonlinearSystem with a constraint NonlinearSystem, are supported. Constraints are not supported in SDESystem or JumpSystem conversion, and so boundary condition species are effectively constant when converting to those model types (but still left as states instead of parameters). Defining constant and boundary condition species is done by

  @parameters k A [isconstantspecies=true]
  @variables t  B(t) [isbcspecies=true] C(t)
  rx = Reaction(k, [A,B], [B,C], [1,2], [1,1])

  Here A is a constant species, B is a non-constant boundary condition species, and C is a normal species. Constant and boundary condition species can be used in creating Reactions like normal species as either substrates or products. Note that network API functions such as netstoichmat, conservationlaws, or reactioncomplexes ignore constant species. i.e. for A a constant species the reaction 2A + B --> C is treated as equivalent to B --> C with a modified rate constant, while B --> A would be identical to B --> 0. Boundary condition species are checked to be balanced by default when ReactionSystems are constructed, i.e.

  rx = Reaction(k, [A,B], [C], [1,2], [1])
  @named rs = ReactionSystem(rs, t)

  would error since B only appears as a substrate. This check can be disabled with

  @named rs = ReactionSystem(rs, t; balanced_bc_check=false)

  Note that network analysis functions assume BC species appear in a balanced manner, so may not work correctly if one appears in an unbalanced fashion. (Conversion to other system types should still work just fine.)

• BREAKING: Added the ability to eliminate conserved species when generating ODEs, nonlinear problems, SDEs, and steady state problems via the remove_conserved=true keyword that can be passed to convert or to ODEProblem, NonlinearProblem, SDEProblem, or SteadyStateProblem when called with a ReactionSystem. For example,

  rn = @reaction_network begin
     k, A + B --> C
     k2, C --> A + B
     end k k2
  osys = convert(ODESystem, rn; remove_conserved=true)
  equations(osys)

  gives

  Differential(t)(A(t)) ~ k2*(_ConLaw[2] - A(t)) - k*(A(t) + _ConLaw[1])*A(t)
  

  Initial conditions should still be specified for all the species in rn, and the conserved constants will then be calculated automatically. Eliminated species are stored as observables in osys and still accessible via solution objects. Breaking as this required modifications to the ReactionSystem type signature.
• BREAKING: Added an internal cache in ReactionSystems for network properties, and revamped many of the network analysis functions to use this cache (so just a ReactionSystem can be passed in). Most of these functions will now only calculate the chosen property the first time they are called, and in subsequent calls will simply returned that cached value. Call reset_networkproperties! to clear the cache and allow properties to be recalculated. The new signatures for rn a ReactionSystem are

  reactioncomplexmap(rn)
  reactioncomplexes(rn)
  complexstoichmat(rn)
  complexoutgoingmat(rn)
  incidencemat(rn)
  incidencematgraph(rn)
  linkageclasses(rn)
  deficiency(rn)
  sns = subnetworks(rn)
  linkagedeficiencies(rn)
  isreversible(rn)
  isweaklyreversible(rn, sns)

  Breaking as this required modifications to the ReactionSystem type signature.
• BREAKING ReactionSystems now store a default value for combinatoric_ratelaws=true. This default value can be set in the ReactionSystem constructor call as a keyword argument. Passing combinatoric_ratelaws as a keyword to convert or problem calls involving a ReactionSystem is still allowed, and will override the ReactionSystem's default.
• Fixed a bug where ODESystem constraint systems did not propagate continuous_events during calls to convert(ODESystem, rn::ReactionSystem).
• Added constant and boundary condition species (in the SBML sense). During conversion constant species are converted to parameters, while boundary condition species are kept as state variables. Note that boundary condition species are treated as constant with respect to reactions, so their dynamics must be defined in a constraint system. Right now only conversion of ReactionSystems to an ODESystem with a constraint ODESystem or NonlinearSystem, or conversion to a NonlinearSystem with a constraint NonlinearSystem, are supported. Constraints are not supported in SDESystem or JumpSystem conversion, and so boundary condition species are effectively constant when converting to those model types (but still left as states instead of parameters). Defining constant and boundary condition species is done by

  @variables t A(t) [isconstant=true] B(t) [isbc=true] C(t)

  Here A is a constant species, B is a boundary condition species, and C is a normal species. Note that network API functions do not make use of these labels, and treat all species as normal -- these properties are only made use of when converting to other system types.

• Added the ability to use symbolic stoichiometry expressions via the DSL. This should now work

  rn = @reaction_network rs begin
    t*k, (α+k+B)*A --> B
    1.0, α*A + 2*B --> k*C + α*D
  end k α

  Here Catalyst will try to preserve the order of symbols within an expression, taking the rightmost as the species and everything multiplying that species as stoichiometry. For example, we can interpret the above reaction as S1 A --> S2 b where S1 = (α+k+B) is the stoichiometry of the reactant A and 1 is the stoichiometry of the reactant B. For

  rn = @reaction_network rs begin
    1.0, 2X*(Y + Z) --> XYZ
  end

  all of X, Y and Z will be registered as species, with substrates (Y,Z) having associated stoichiometries of (2X,2X). As for rate expressions, any symbols that appear and are not defined as parameters will be declared to be species.

  In contrast, when declaring reactions

  rx = @reaction t*k, (k+α)*A --> B

  will work, with every symbol declared a parameter except the leftmost symbol in the reaction line. So

  rx = @reaction 1.0, 2X*(Y + Z) --> XYZ

  will make X a parameter and Y, Z and XYZ species.

• Symbolic stoichiometry supports interpolation of expressions in @reaction_network and @reaction.

• Added the ability to use symbolic variables, parameters and expressions for stoichiometric coefficients. See the new tutorial on Parametric Stoichiometry for details, and note the caveat about ModelingToolkit converting integer parameters to floating point types that must be worked around to avoid calls to factorial that involve floats.

• Added the ability to use floating point stoichiometry (currently only tested for generating ODE models). This should now work

  rn = @reaction_network begin
    k, 2.5*A --> 3*B
  end k

  or directly

  @parameters k b
  @variables t A(t) B(t) C(t) D(t)
  rx1 = Reaction(k,[B,C],[B,D], [2.5,1],[3.5, 2.5])
  rx2 = Reaction(2*k, [B], [D], [1], [2.5])
  rx3 = Reaction(2*k, [B], [D], [2.5], [2])
  @named mixedsys = ReactionSystem([rx1,rx2,rx3],t,[A,B,C,D],[k,b])
  osys = convert(ODESystem, mixedsys; combinatoric_ratelaws=false)

  Note, when using convert(ODESystem, mixedsys; combinatoric_ratelaws=false) the combinatoric_ratelaws=false parameter must be passed. This is also true when calling ODEProblem(mixedsys,...; combinatoric_ratelaws=false). This disables Catalyst's standard rescaling of reaction rates when generating reaction rate laws, see the docs. Leaving this out for systems with floating point stoichiometry will give an error message.

• Added @reaction macro

  rx = @reaction k*v, A + B --> C + D
  
  # is equivalent to
  @parameters k v
  @variables t A(t) B(t) C(t) D(t)
  rx == Reaction(k*v, [A,B], [C,D])

  Here k and v will be parameters and A, B, C and D will be variables. Interpolation of existing parameters/variables also works

  @parameters k b
  @variables t A(t)
  ex = k*A^2 + t
  rx = @reaction b*$ex*$A, $A --> C

  Any symbols arising in the rate expression that aren't interpolated are treated as parameters, while any in the reaction part (A + B --> C + D) are treated as species.

• Added symmap_to_varmap, setdefaults!, and updated all *Problem(rn,...) calls to allow setting initial conditions and parameter values using symbol maps. See the Catalyst API for details. These allow using regular Julia Symbols to specify parameter values and initial conditions. i.e. to set defaults we can do

  rn = @reaction_network begin
      α, S + I --> 2I
      β, I --> R
  end α β
  setdefaults!(rn, [:S => 999.0, :I => 1.0, :R => 0.0, :α => 1e-4, :β => .01])
  op    = ODEProblem(rn, [], (0.0,250.0), [])
  sol   = solve(op, Tsit5())

  To explicitly pass initial conditions and parameters using symbols we can do

  rn = @reaction_network begin
      α, S + I --> 2I
      β, I --> R
  end α β
  u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
  p  = (:α => 1e-4, :β => .01)
  op    = ODEProblem(rn, u0, (0.0,250.0), p)
  sol   = solve(op, Tsit5())

  In each case ModelingToolkit symbolic variables can be used instead of Symbols, e.g.

  @parameters α β
  @variables t S(t) I(t) R(t)
  setdefaults!(rn, [S => 999.0, I => 1.0, R => 0.0, α => 1e-4, β => .01])

• BREAKING: The order of the parameters in the ReactionSystem's .ps field has been changed (only when created through the @reaction_network macro). Previously they were ordered according to the order with which they appeared in the macro. Now they are ordered according the to order with which they appeared after the end part. E.g. in

  rn = @reaction_network begin
    (p,d), 0 <--> X
  end d p

  previously the order was [p,d], while now it is [d, p].

• Added support for @unpack observable_variable = rn and rn.observable_variable. This requires a new inner constructor definition for ReactionSystems, but is not considered breaking as the inner constructor is considered private.
• Support added for ModelingToolkit 7 and Symbolics 4.

• ReactionSystem(rxs::Vector{Reaction}, t) should now work and will infer the species and parameters.
• BREAKING: Any undeclared variables in the DSL are now inferred to be species. i.e. this no longer errors, and B is assumed to be a species

  rn = @reaction_network begin
    k*B, A --> C
  end k

• BREAKING: Internal changes mean the order of species or parameters in generated systems may have changed. Changes that induce different orders will not be considered breaking in the future.
• Added interpolation in the DSL for species, variables, and the network name. i.e. this is now valid

  @parameters k
  @variables t, A(t)
  spec = A
  rate = k*A
  name = :network
  rn = @reaction_network $name begin
    $rate*B, 2*$spec + B --> $spec + C
    end

• Added the ability to compose ReactionSystems via subsystems, and include either ODESystems or NonlinearSystems as subsystems. Note, if using non-ReactionSystem subsystems it is not currently possible to convert to a JumpSystem or SDESystem. It is also not possible to include either SDESystems or JumpSystems as subsystems.
• Added extend(sys, reactionnetwork, name=nameof(sys)) to extend ReactionSystems with constraint equations (algebraic equations or ODEs), or other ReactionSystems. Algebraic or differential constraints are stored as a NonlinearSystem or ODESystem within the ReactionSystem, and accessible via get_constraints(reactionnetwork).
• Added Catalyst.flatten(rn) to allow flattening of a ReactionSystem with sub-systems into one ReactionSystem. Non-ReactionSystem subsystems are merged into the constraints of the flattened ReactionSystem, and accessible via get_constraints.
• BREAKING: ReactionSystems are now always flattened when calling convert. This should only affect models that use subsystems.
• Added incidencematgraph, linkageclasses, deficiency, subnetworks, linkagedeficiency, isreversible and isweaklyreversible API functions.
• Deprecated merge, use ModelingToolkit.extend instead.
• Deprecated params and numparams (use ModelingToolkit.parameters to get all parameters of a system and all subsystems, or use reactionparams to get all parameters of a system and all ReactionSystem subsystems. The latter correspond to those parameters used within Reactions.)
• BREAKING: Added a custom hash for Reactions to ensure they work in Dicts and Sets properly, ensuring set-type comparisons between collections of Reactions work.
• Updated the docs and added a new tutorial on using compositional tooling.

1. BREAKING: netstoichmat, prodstoichmat and substoichmat are now transposed to be number of species by number of reactions. This is more consistent with the chemical reaction network literature for stoichiometry matrices.

2. reactioncomplexmap added to provide a mapping from reaction complexes to reactions they participate in.

3. Most API *mat functions now take an optional sparse keyword argument. If passed sparse=true a sparse matrix representation is generated, otherwise the default sparse=false value returns dense Matrix representations.

1. Network representations for the reaction complexes of a system along with associated graph functionality:

rn = @reaction_network begin
           k₁, 2A --> B
           k₂, A --> C
           k₃, C --> D
           k₄, B + D --> E
           k₅, B --> E
           k₆, D --> C
     end k₁ k₂ k₃ k₄ k₅ k₆
smap  = speciesmap(rn)
rcs,B = reactioncomplexes(rn; smap=smap)
Z     = complexstoichmat(rn; rcs=rcs)
Δ     = complexoutgoingmat(rn; B=B)
complexgraph(rn; complexdata=(rcs,B))

which gives

2. Support for units via ModelingToolkit and Unitful.jl in directly constructed ReactionSystems:

# ]add Unitful
using Unitful
@parameters α [unit=u"μM/s"] β [unit=u"s"^(-1)] γ [unit=u"μM*s"^(-1)]
@variables t [unit=u"s"] A(t) [unit=u"μM"] B(t) [unit=u"μM"] C(t) [unit=u"μM"]
rxs = [Reaction(α, nothing, [A]),
       Reaction(β, [A], [B]),
       Reaction(γ, [A,B], [B], [1,1], [2])]
@named rs = ReactionSystem(rxs, t, [A,B,C], [α,β,γ])

By default, during construction of rs Catalyst will call

validate(rs)

which will print warnings and return false if either

1. The species(rs) do not all have the same units.
2. The implicit (ODE) rate laws for each reaction do not have units of (species units) / (time units), where the time units are the units of t.

(Note, at this time the @reaction_network macro does not support units.)

3. Calculation of conservation laws

rn = @reaction_network begin
  (k₊,k₋), A + B <--> C
  end k₊ k₋
clawmat = conservationlaws(netstoichmat(rn))

giving

 1  -1  0
 0   1  1

and

cquants = conservedquantities(species(rn), clawmat)

giving

 A(t) - B(t)
 B(t) + C(t)

See the API docs for more details about each of these new features.

1. Basic unit validation has been added following its addition for all ModelingToolkit systems.

1. reactioncomplexes, ReactionComplex, reactionrates, complexstoichmat and complexoutgoingmat are added to allow the calculation of reaction complex-based network matrix representations.

BREAKING: This is a breaking release, with all ModelingToolkit ReactionSystem and Reaction functionality migrated to Catalyst.

1. Plain text arrows "<--" and "<-->" for backward and reversible reactions are available if using Julia 1.6 or higher:

rn = @reaction_network begin
  (k1,k2), A + B <--> C
  k3, 0 <-- C
end k1 k2 k3

2. BREAKING: Reaction networks can be named

rn = @reaction_network Reversible_Reaction begin
  k1, A --> B
  k2, B --> A
  end k1 k2
ModelingToolkit.nameof(rn) == :Reversible_Reaction

Note, empty networks can no longer be created with parameters, i.e. only

rn = @reaction_network          # uses a randomly generated name
rn = @reaction_network MyName   # is named MyName

are allowed.

3. Compositional modeling with generated ODESystems, see here for an example that composes three gene modules to make the repressilator.