DiffEqBiological.jl API

# DiffEqBiological.jl API

## Reaction Network Generation Macros

DiffEqBiological has three macros for generating reaction networks. @reaction_network generates a complete network, including everything needed to construct ODE, SDE and Jump problems directly from the network. For small systems it is the recommended form to use.

@min_reaction_network constructs a network that just stores the basic information needed to represent the species, parameters and chemical reactions. This is sufficient for network analysis, such as calculating dependency graphs, but means the network must be extended to build mathematical models, see addodes!, addsdes!, and addjumps!.

@empty_reaction_network constructs an empty network. Both min_reaction_networks and empty_reaction_networks can be enlarged using addspecies!, addparam!, and addreaction!. Once the final chemistry for the network is set, addodes!, addsdes!, and addjumps! can be called to build corresponding mathematical models.

It is important to note for @reaction_network and @min_reaction_network that species which are used within the macro as part of a rate expression, but not as a substrate or product of some reaction, are not recognized as either a species or parameter. i.e. avoid

rn = @reaction_network begin
k*X, Y --> W
end k

as here X is never defined as either a species or parameter. This leads to internal problems in the representation of reactions that can not be corrected by subsequently calling addspecies!.

@reaction_network

Generates a subtype of an AbstractReactionNetwork that encodes a chemical reaction network, and complete ODE, SDE and jump representations of the system. See the Chemical Reaction Model docs for details on parameters to the macro.

@min_reaction_network

Generates a subtype of an AbstractReactionNetwork that only encodes a chemical reaction network. Use addodes!, addsdes! or addjumps! to complete the network for specific problem types. It accepts the same arguments as @reaction_network.

@empty_reaction_network networktype

Generates a subtype of an AbstractReactionNetwork that encodes an empty chemical reaction network. networktype is an optional parameter that specifies the type of the generated network. Use addspecies!, addparam! and addreaction! to extend the network. Use addodes!, addsdes! or addjumps! to complete the network for specific problem types.

## Basic properties

species(network)

Given an AbstractReactionNetwork, return a vector of species symbols.

speciesmap(network)

Given an AbstractReactionNetwork, return a Dictionary mapping from species symbol to species index.

params(network)

Given an AbstractReactionNetwork, return a vector of parameter symbols.

paramsmap(network)

Given an AbstractReactionNetwork, return a Dictionary mapping from parameter symbol to parameter index.

numspecies(network)

Return the number of species within the given AbstractReactionNetwork.

numparams(network)

Return the number of parameters within the given AbstractReactionNetwork.

numreactions(network)

Return the number of reactions within the given AbstractReactionNetwork.

## Reaction Properties

substrates(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a vector of symbols of species that correspond to substrates in the reaction. i.e. for

k*W, X + 3Y --> X + W

the returned vector would be [:X,:Y].

Allocates a new vector to store the symbols.

products(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a vector of symbols of species that correspond to products in the reaction. i.e. for

k*W, X + 3Y --> X + W

the returned vector would be [:X,:W].

Allocates a new vector to store the symbols.

dependents(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a vector of symbols of species the reaction rate law depends on. i.e. for

k*W, 2X + 3Y --> 5Z + W

the returned vector would be [:W,:X,:Y].

Non-allocating, returns underlying field within the reaction_network.

dependants(network, rxidx)

See documentation for dependents(network, rxidx).

ismassaction(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a boolean indicating whether the given reaction is of mass action form. For example, the reaction

2*k, 2X + 3Y --> 5Z + W

would return true, while reactions with state-dependent rates like

k*X, X + Y --> Z

would return false.

Non-allocating, returns underlying field within the reaction_network.

substratestoich(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a vector of pairs, mapping ids of species that serve as substrates in the reaction to the corresponding stoichiometric coefficient as a substrate.

Allocates a new vector to store the pairs.

substratesymstoich(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a Vector of ReactantStructs, mapping the symbols of species that serve as substrates in the reaction to the corresponding stoichiometric coefficient as a substrate.

Non-allocating, returns underlying field within the reaction_network.

productstoich(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a vector of pairs, mapping ids of species that are products in the reaction to the corresponding stoichiometric coefficient as a product.

Allocates a new vector to store the pairs.

productsymstoich(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a Vector of ReactantStructs, mapping the symbols of species that are products in the reaction to the corresponding stoichiometric coefficient as a product.

Non-allocating, returns underlying field within the reaction_network.

netstoich(network, rxidx)

Given an AbstractReactionNetwork and a reaction index, rxidx, return a vector of pairs, mapping ids of species that change numbers due to the reaction to the net stoichiometric coefficient of the species (i.e. net change in the species due to the reaction).

Allocates a new vector to store the pairs.

## Functions to Add Species, Parameters and Reactions to a Network

Both @min_reaction_network and @empty_reaction_network can be extended with additional species, parameters, and reactions.

Note, always add all species and parameter definitions before adding any reaction definitions. Other orderings may result in incorrect information stored within the generated network.

addspecies!(network, speciessym::Symbol)

Given an AbstractReaction network, add the species corresponding to the passed in symbol to the network (if it is not already defined).

addspecies!(network, speciesname::String)

Given an AbstractReaction network, add the species with name given by the passed in string to the network (if it is not already defined.

addparam!(network, param::Symbol)

Given an AbstractReaction network, add the parameter corresponding to the passed in symbol to the network (if it is not already defined).

addparam!(network, paramname::String)

Given an AbstractReaction network, add the parameter with name given by the passed in string to the network (if it is not already defined).

addreaction!(network, rateex::Union{Expr,Symbol,Int,Float64}, rxexpr::Expr)

Given an AbstractReaction network, add a reaction with the passed in rate and reaction expressions. i.e. a reaction of the form

k*X, 2X + Y --> 2W

would have rateex=:(k*X) and rxexpr=:(2X + Y --> W),

10.5, 0 --> X

would have rateex=10.5 and rxexpr=:(0 --> X), and

k, X+X --> Z

would have rateex=:k and rxexpr=:(X+X --> Z). All normal DSL reaction definition notation should be supported.

addreaction!(network, rateex::Union{Expr,Symbol,Int,Float64}, substrates, products)

Given an AbstractReaction network, add a reaction with the passed in rate, rateex, substrate stoichiometry, and product stoichiometry. Stoichiometries are represented as tuples of Pair{Symbol,Int}. i.e. a reaction of the form

k*X, 2X + Y --> 2W

would have rateex=:(k*X), substrates=(:X=>2, :Y=>2)andproducts=(W=>2,),

10.5, 0 --> X

would have rateex=10.5, substrates=() and products=(:X=>1,), and

k, X+X --> Z

would have rateex=:k, substrates=(:X=>2,) and products=(:Z=>2,). All normal DSL reaction definition notation should be supported for the rateex.

add_scale_noise_param!(network, scale_noise::Symbol)

Given an AbstractReaction network, add the parameter corresponding to the passed in symbol to the network (if it is not already defined), and register it as the noise scaling coefficient.

add_scale_noise_param!(network, scale_noise_name::String)

Given an AbstractReaction network, add the parameter with the passed in string as its name to the network (if it is not already defined), and register it as the noise scaling coefficient.

## Functions to Add ODEs, SDEs or Jumps to a Network

addodes!(network; build_jac=true, build_symfuncs=true)

Extend an AbstractReactionNetwork generated with the @min_reaction_network or @empty_reaction_network macros with everything needed to use ODE solvers.

Optional kwargs can be used to disable the construction of additional ODE solver components.

addsdes!(network; build_jac=true, build_symfuncs=true)

Extend an AbstractReactionNetwork generated with the @min_reaction_network or @empty_reaction_network macros with everything needed to use SDE solvers.

Optional kwargs can be used to disable the construction of additional SDE solver components.

addjumps!(network; build_jumps=true, build_regular_jumps=true, minimal_jumps=false)

Extend an AbstractReactionNetwork generated with the @min_reaction_network or @empty_reaction_network macros with everything needed to use jump SSA solvers.

Optional kwargs can be used to disable the construction of additional jump solver components.

Keyword arguments:

• build_jumps: if true jump rates and affects will be calculated for use in DiffEqJump SSAs.

• build_regular_jumps: if true a RegularJump representation of the stochastic chemical kinetics model will be calculated for use in τ-leaping methods.

• minimal_jumps: if true ConstantRate jumps are only constructed for non-mass action jumps. (Note, mass action jumps are still resolved within any jump simulation. This option simply speeds up the construction of the jump problem since it avoids building redundant ConstantRate jumps that encode MassActionJumps, which are subsequently ignored within jump simulations.)

## Generated Functions for Models

oderhsfun(network)

Given an AbstractReactionNetwork, return a function, f!(du,u,p,t), that evaluates the current value of the ODE model derivative functions, $du/dt = f(u,t)$, within du.

Note, for a network generated with the @min_reaction_network macro addodes! must be called first.

jacfun(network)

Given an AbstractReactionNetwork, return a function, jac!(J,u,p,t), that evaluates the current Jacobian matrix, J, of the ODE model, $du/dt = f(u,t)$. The Jacobian matrix has entries

$J_{i,j} = \partial f_i(u,t) / \partial u_j$.

Note, for a network generated with the @min_reaction_network macro addodes! must be called first.

paramjacfun(network)

Given an AbstractReactionNetwork, return a function, pjac(pJ,u,p,t), that evaluates the current parameter Jacobian matrix, pJ, of the ODE model, $du/dt = f(u,t)$. The parameter Jacobian matrix has entries

$pJ_{i,j} = \partial f_i(u,t) / \partial p_j$.

Note, for a network generated with the @min_reaction_network macro addodes! must be called first.

odefun(network)

Given an AbstractReactionNetwork, return a DiffEqBase.ODEFunction encoding an ODE model for the reaction network.

Note, for a network generated with the @min_reaction_network macro addodes! must be called first.

noisefun(network)

Given an AbstractReactionNetwork, return a function, g(η,u,p,t), that evaluates the current noise coefficients for each reaction in the Chemical Langevin Equation representation within η.

Note, for a network generated with the @min_reaction_network macro addsdes! must be called first.

sdefun(network)

Given an AbstractReactionNetwork, return a DiffEqBase.SDEFunction encoding a Chemical Langevin Equation SDE model for the reaction network.

Note, for a network generated with the @min_reaction_network macro addsdes! must be called first.

jumps(network)

Given an AbstractReactionNetwork, return a tuple of AbstractJumps encoding a stochastical chemical kinetics representation for the reaction network.

Note, for a network generated with the @min_reaction_network macro addjumps! must be called first.

regularjumps(network)

Given an AbstractReactionNetwork, return a RegularJump encoding a stochastical chemical kinetics representation of the reaction network for use in $\tau$-leaping approximations.

Note, for a network generated with the @min_reaction_network macro addjumps! must be called first.

## Generated Expressions

odeexprs(network)

Given an AbstractReactionNetwork, return a vector of the ODE expressions.

Note, for a network generated with the @min_reaction_network macro addodes! must be called first.

jacobianexprs(network)

Given an AbstractReactionNetwork, return a matrix with the ODE Jacobian expressions.

Note, for a network generated with the @min_reaction_network macro addodes! must be called first.

noiseexprs(network)

Given an AbstractReactionNetwork, return a vector of the SDE noise expressions for each reaction.

Note, for a network generated with the @min_reaction_network macro addsdes! must be called first.

jumpexprs(network)

Given an AbstractReactionNetwork, return a tuple of the jump rates and affects expressions.

Note, for a network generated with the @min_reaction_network macro addjumps! must be called first.

rateexpr(network, rxidx)

Given an AbstractReactionNetwork, return the reaction rate expression for the reaction with index rxidx. Note, for a reaction defined by

k*X*Y, X+Z --> 2X + Y

the expression that is returned will be :(k*X*Y), while the rate law used in ODEs and SDEs would be k*X^2*Y*Z.

oderatelawexpr(network, rxidx)

Given an AbstractReactionNetwork, return the reaction rate law expression used in generated ODEs for the reaction with index rxidx. 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 expression that is returned will be :(k * (X^2/2) * (Y^3/6)).

ssaratelawexpr(network, rxidx)

Given an AbstractReactionNetwork, return the reaction rate law expression used in generated stochastic chemical kinetic model SSAs for the reaction with index rxidx. 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 expression that is returned will be :(k * binomial(X,2) * binomial(Y,3)).

## Network Comparison Functions

==(rn1::DiffEqBase.AbstractReactionNetwork, rn2::DiffEqBase.AbstractReactionNetwork)

Tests whether the underlying species symbols, parameter symbols and reactions are the same in the two networks. Ignores order network components were defined, so the integer id of any individual species/parameters/reactions may be different between the two networks. Does not currently account for different reaction definitions, so "k, X+Y –> Y + Z" will not be the same as "k, Y+X –> Y + Z"

## Dependency Graphs

rxtospecies_depgraph(network)

Given an AbstractReactionNetwork, returns a Vector{Vector{Int}} mapping each reaction index to the indices of species that depend on it.

speciestorx_depgraph(network)

Given an AbstractReactionNetwork, returns a Vector{Vector{Int}} mapping each species index to the indices of reactions that depend on it.

rxtorx_depgraph(network)

Given an AbstractReactionNetwork, returns a Vector{Vector{Int}}` mapping each reaction index to the indices of reactions that depend on it.