- Added
@reaction
macroHererx = @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])
k
andv
will be parameters andA
,B
,C
andD
will be variables. Interpolation of existing parameters/variables also worksAny symbols arising in the rate expression that aren't interpolated are treated as parameters, while any in the reaction part (@parameters k b @variables t A(t) ex = k*A^2 + t rx = @reaction b*$ex*$A, $A --> C
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 JuliaSymbols
to specify parameter values and initial conditions. i.e. to set defaults we can doTo explicitly pass initial conditions and parameters using symbols we can dorn = @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())
In each case ModelingToolkit symbolic variables can be used instead ofrn = @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())
Symbol
s, 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 appeard after theend
part. E.g. inpreviously the order wasrn = @reaction_network begin (p,d), 0 <--> X end d p
[p,d]
, while now it is[d, p]
.
- Added support for
@unpack observable_variable = rn
andrn.observable_variable
. This requires a new inner constructor definition forReactionSystem
s, 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 speciesrn = @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
ReactionSystem
s via subsystems, and include eitherODESystem
s orNonlinearSystem
s as subsystems. Note, if using non-ReactionSystem
subsystems it is not currently possible to convert to aJumpSystem
orSDESystem
. It is also not possible to include eitherSDESystem
s orJumpSystems
as subsystems. - Added
extend(sys, reactionnetwork, name=nameof(sys))
to extendReactionSystem
s with constraint equations (algebraic equations or ODEs), or otherReactionSystem
s. Algebraic or differential constraints are stored as aNonlinearSystem
orODESystem
within theReactionSystem
, and accessible viaget_constraints(reactionnetwork)
. - Added
Catalyst.flatten(rn)
to allow flattening of aReactionSystem
with sub-systems into oneReactionSystem
. Non-ReactionSystem
subsystems are merged into the constraints of the flattenedReactionSystem
, and accessible viaget_constraints
. - BREAKING:
ReactionSystem
s are now always flattened when callingconvert
. This should only affect models that usesubsystem
s. - Added
incidencematgraph
,linkageclasses
,deficiency
,subnetworks
,linkagedeficiency
,isreversible
andisweaklyreversible
API functions. - Deprecated
merge
, useModelingToolkit.extend
instead. - Deprecated
params
andnumparams
(useModelingToolkit.parameters
to get all parameters of a system and all subsystems, or usereactionparams
to get all parameters of a system and allReactionSystem
subsystems. The latter correspond to those parameters used withinReaction
s.) - BREAKING: Added a custom
hash
forReaction
s to ensure they work inDict
s andSet
s properly, ensuring set-type comparisons between collections ofReaction
s 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
Uniftul.jl in directly constructed
ReactionSystem
s:
# ]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
- The
species(rs)
do not all have the same units. - 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
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 ODESystem
s, see
here
for an example that composes three gene modules to make the repressilator.