ReactionSystem
A ReactionSystem
represents a system of chemical reactions. Conversions are provided to generate corresponding chemical reaction ODE models, chemical Langevin equation SDE models, and stochastic chemical kinetics jump process models. As a simple example, the code below creates a SIR model, and solves the corresponding ODE, SDE, and jump process models.
using ModelingToolkit, OrdinaryDiffEq, StochasticDiffEq, DiffEqJump
@parameters β γ t
@variables S(t) I(t) R(t)
rxs = [Reaction(β, [S,I], [I], [1,1], [2])
Reaction(γ, [I], [R])]
rs = ReactionSystem(rxs, t, [S,I,R], [β,γ])
u₀map = [S => 999.0, I => 1.0, R => 0.0]
parammap = [β => 1/10000, γ => 0.01]
tspan = (0.0, 250.0)
# solve as ODEs
odesys = convert(ODESystem, rs)
oprob = ODEProblem(odesys, u₀map, tspan, parammap)
sol = solve(oprob, Tsit5())
# solve as SDEs
sdesys = convert(SDESystem, rs)
sprob = SDEProblem(sdesys, u₀map, tspan, parammap)
sol = solve(sprob, EM(), dt=.01)
# solve as jump process
jumpsys = convert(JumpSystem, rs)
u₀map = [S => 999, I => 1, R => 0]
dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)
jprob = JumpProblem(jumpsys, dprob, Direct())
sol = solve(jprob, SSAStepper())
System Constructors
ModelingToolkit.Reaction
— Typestruct Reaction{S, T<:Number}
One chemical reaction.
Fields
rate
The rate function (excluding mass action terms).
substrates
Reaction substrates.
products
Reaction products.
substoich
The stoichiometric coefficients of the reactants.
prodstoich
The stoichiometric coefficients of the products.
netstoich
The net stoichiometric coefficients of all species changed by the reaction.
only_use_rate
false
(default) ifrate
should be multiplied by mass action terms to give the rate law.true
ifrate
represents the full reaction rate law.
Examples
using ModelingToolkit
@parameters t k[1:20]
@variables A(t) B(t) C(t) D(t)
rxs = [Reaction(k[1], nothing, [A]), # 0 -> A
Reaction(k[2], [B], nothing), # B -> 0
Reaction(k[3],[A],[C]), # A -> C
Reaction(k[4], [C], [A,B]), # C -> A + B
Reaction(k[5], [C], [A], [1], [2]), # C -> A + A
Reaction(k[6], [A,B], [C]), # A + B -> C
Reaction(k[7], [B], [A], [2], [1]), # 2B -> A
Reaction(k[8], [A,B], [A,C]), # A + B -> A + C
Reaction(k[9], [A,B], [C,D]), # A + B -> C + D
Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]), # A+3B+4C -> 2C + 3D
Reaction(k[13], [A,B], nothing, [3,1], nothing), # 3A+B -> 0
Reaction(k[14], nothing, [A], nothing, [2]), # 0 -> 2A
Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true), # A -> 0 with custom rate
Reaction(k[16], [A], [B]; only_use_rate=true), # A -> B with custom rate.
Reaction(k[17]*A*exp(B), [C], [D], [2], [1]), # 2C -> D with non constant rate.
Reaction(k[18]*B, nothing, [B], nothing, [2]), # 0 -> 2B with non constant rate.
Reaction(k[19]*t, [A], [B]), # A -> B with non constant rate.
Reaction(k[20]*t*A, [B,C], [D],[2,1],[2]) # 2A +B -> 2C with non constant rate.
]
Notes:
nothing
can be used to indicate a reaction that has no reactants or no products. In this case the corresponding stoichiometry vector should also be set tonothing
.- The three-argument form assumes all reactant and product stoichiometric coefficients are one.
ModelingToolkit.ReactionSystem
— Typestruct ReactionSystem <: ModelingToolkit.AbstractSystem
A system of chemical reactions.
Fields
eqs
The reactions defining the system.
iv
Independent variable (usually time).
states
Dependent (state) variables representing amount of each species.
ps
Parameter variables.
observed
name
The name of the system
systems
systems: The internal systems
Example
Continuing from the example in the Reaction
definition:
rs = ReactionSystem(rxs, t, [A,B,C,D], k)
Composition and Accessor Functions
sys.eqs
orequations(sys)
: The reactions that define the system.sys.states
orstates(sys)
: The set of chemical species in the system.sys.parameters
orparameters(sys)
: The parameters of the system.sys.iv
orindependent_variable(sys)
: The independent variable of the reaction system, usually time.
Query Functions
ModelingToolkit.oderatelaw
— Functionoderatelaw(rx; combinatoric_ratelaw=true)
Given a Reaction
, return the reaction rate law Operation
used in generated ODEs for the reaction. Note, for a reaction defined by
k*X*Y, X+Z --> 2X + Y
the expression that is returned will be k*X(t)^2*Y(t)*Z(t)
. For a reaction of the form
k, 2X+3Y --> Z
the Operation
that is returned will be k * (X(t)^2/2) * (Y(t)^3/6)
.
Notes:
- Allocates
combinatoric_ratelaw=true
uses factorial scaling factors in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S^2/2!
. Ifcombinatoric_ratelaw=false
then the ratelaw isk*S^2
, i.e. the scaling factor is ignored.
ModelingToolkit.jumpratelaw
— Functionjumpratelaw(rx; rxvars=get_variables(rx.rate), combinatoric_ratelaw=true)
Given a Reaction
, return the reaction rate law Operation
used in generated stochastic chemical kinetics model SSAs for the reaction. Note, for a reaction defined by
k*X*Y, X+Z --> 2X + Y
the expression that is returned will be k*X^2*Y*Z
. For a reaction of the form
k, 2X+3Y --> Z
the Operation
that is returned will be k * binomial(X,2) * binomial(Y,3)
.
Notes:
rxvars
should give theVariable
s, i.e. species and parameters, the rate depends on.- Allocates
combinatoric_ratelaw=true
uses binomials in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S*(S-1)/2
. Ifcombinatoric_ratelaw=false
then the ratelaw isk*S*(S-1)
, i.e. the rate law is not normalized by the scaling factor.
ModelingToolkit.ismassaction
— Functionismassaction(rx, rs; rxvars = get_variables(rx.rate),
haveivdep = any(var -> isequal(get_iv(rs),var), rxvars),
stateset = Set(states(rs)))
True if a given reaction is of mass action form, i.e. rx.rate
does not depend on any chemical species that correspond to states of the system, and does not depend explicitly on the independent variable (usually time).
Arguments
rx
, theReaction
.rs
, aReactionSystem
containing the reaction.- Optional:
rxvars
,Variable
s which are not inrxvars
are ignored as possible dependencies. - Optional:
haveivdep
,true
if theReaction
rate
field explicitly depends on the independent variable. - Optional:
stateset
, set of states which if the rxvars are within mean rx is non-mass action.
Transformations
Base.convert
— FunctionBase.convert(::Type{<:ODESystem},rs::ReactionSystem)
Convert a ReactionSystem
to an ODESystem
.
Notes:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate
law, i.e. for 2S -> 0
at rate k
the ratelaw would be k*S^2/2!
. If combinatoric_ratelaws=false
then the ratelaw is k*S^2
, i.e. the scaling factor is ignored.
Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
Convert a ReactionSystem
to an NonlinearSystem
.
Notes:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate
law, i.e. for 2S -> 0
at rate k
the ratelaw would be k*S^2/2!
. If combinatoric_ratelaws=false
then the ratelaw is k*S^2
, i.e. the scaling factor is ignored.
Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
Convert a ReactionSystem
to an SDESystem
.
Notes:
combinatoric_ratelaws=true
uses factorial scaling factors in calculating the rate
law, i.e. for 2S -> 0
at rate k
the ratelaw would be k*S^2/2!
. If combinatoric_ratelaws=false
then the ratelaw is k*S^2
, i.e. the scaling factor is ignored.
noise_scaling=nothing::Union{Vector{Operation},Operation,Nothing}
allows for linear
scaling of the noise in the chemical Langevin equations. If nothing
is given, the default value as in Gillespie 2000 is used. Alternatively, an Operation
can be given, this is added as a parameter to the system (at the end of the parameter array). All noise terms are linearly scaled with this value. The parameter may be one already declared in the ReactionSystem
. Finally, a Vector{Operation}
can be provided (the length must be equal to the number of reactions). Here the noise for each reaction is scaled by the corresponding parameter in the input vector. This input may contain repeat parameters.
Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
Convert a ReactionSystem
to an JumpSystem
.
Notes:
combinatoric_ratelaws=true
uses binomials in calculating the rate law, i.e. for2S -> 0
at ratek
the ratelaw would bek*S*(S-1)/2
. Ifcombinatoric_ratelaws=false
then the ratelaw isk*S*(S-1)
, i.e. the rate law is not normalized by the scaling factor.