Edit on GitHub

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 — Type
struct Reaction{S<:Variable, 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) if rate should be multiplied by mass action terms to give the rate law. true if rate 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 to nothing.
  • The three-argument form assumes all reactant and product stoichiometric coefficients are one.
source
ModelingToolkit.ReactionSystem — Type
struct 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.

  • pins

  • 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)
source

Composition and Accessor Functions

  • sys.eqs or equations(sys): The reactions that define the system.
  • sys.states or states(sys): The set of chemical species in the system.
  • sys.parameters or parameters(sys): The parameters of the system.
  • sys.iv or independent_variable(sys): The independent variable of the reaction system, usually time.

Query Functions

ModelingToolkit.oderatelaw — Function
oderatelaw(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. for 2S -> 0 at rate k the ratelaw would be k*S^2/2!. If combinatoric_ratelaw=false then the ratelaw is k*S^2, i.e. the scaling factor is ignored.
source
ModelingToolkit.jumpratelaw — Function
jumpratelaw(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 the Variables, i.e. species and parameters, the rate depends on.
  • Allocates
  • combinatoric_ratelaw=true uses binomials in calculating the rate law, i.e. for 2S -> 0 at rate k the ratelaw would be k*S*(S-1)/2. If combinatoric_ratelaw=false then the ratelaw is k*S*(S-1), i.e. the rate law is not normalized by the scaling factor.
source
ModelingToolkit.ismassaction — Function
ismassaction(rx, rs; rxvars = get_variables(rx.rate),
                              haveivdep = any(var -> isequal(rs.iv,convert(Variable,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, the Reaction.
  • rs, a ReactionSystem containing the reaction.
  • Optional: rxvars, Variables which are not in rxvars are ignored as possible dependencies.
  • Optional: haveivdep, true if the Reaction 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.
source

Transformations

Base.convert — Function
Base.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.

source
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.

source
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. for 2S -> 0 at rate k the ratelaw would be k*S*(S-1)/2. If combinatoric_ratelaws=false then the ratelaw is k*S*(S-1), i.e. the rate law is not normalized by the scaling factor.
source
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.

source