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

`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)`

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

`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`Variable`

s, 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.

`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`

,`Variable`

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

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

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

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