From ba8189246316bc78ebd58e18d71f23de9566e42b Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 16 May 2024 16:42:05 -0400 Subject: [PATCH] rebase --- docs/src/02_model_creation/01_dsl_basics.md | 370 ----------- docs/src/02_model_creation/02_dsl_advanced.md | 626 ------------------ 2 files changed, 996 deletions(-) delete mode 100644 docs/src/02_model_creation/01_dsl_basics.md delete mode 100644 docs/src/02_model_creation/02_dsl_advanced.md diff --git a/docs/src/02_model_creation/01_dsl_basics.md b/docs/src/02_model_creation/01_dsl_basics.md deleted file mode 100644 index a230c2de00..0000000000 --- a/docs/src/02_model_creation/01_dsl_basics.md +++ /dev/null @@ -1,370 +0,0 @@ -# [The Catalyst DSL - Introduction](@id dsl_description) -In the [introduction to Catalyst](@ref introduction_to_catalyst) we described how the `@reaction_network` [macro](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros) can be used to create chemical reaction network (CRN) models. This macro enables a so-called [domain-specific language](https://en.wikipedia.org/wiki/Domain-specific_language) (DSL) for creating CRN models. This tutorial will give a basic introduction on how to create Catalyst models using this macro (from now onwards called the "*Catalyst DSL*"). A [follow-up tutorial](@ref ref) will describe some of the DSL's more advanced features. - -The Catalyst DSL generates a [`ReactionSystem`](@ref) (the [julia structure](https://docs.julialang.org/en/v1/manual/types/#Composite-Types) Catalyst uses to represent CRN models). These can be created through alternative methods (e.g. [programmatically](@ref programmatic_CRN_construction) or [compositionally](@ref compositional_modeling)). A summary of the various ways to create `ReactionSystems`s can be found [here](@ref ref). [Previous](@ref ref) and [following](@ref ref) tutorials describe how to simulate models once they have been created using the DSL. This tutorial will solely focus on model creation. - -Before we begin, we will first load the Catalyst package (which is required to run the code). -```@example dsl_1 -using Catalyst -``` - -### [Quick-start summary](@id dsl_description_quick_start) -The DSL is initiated through the `@reaction_network`, which is followed by one line for each reaction. Each reaction consists of a *rate*, followed lists first of the substrates and next of the products. E.g. a Michaelis-Menten enzyme kinetics system can be written as -```@example dsl_0 -rn = @reaction_network begin - (kB,kD), S + E <--> SE - kP, SE --> P + E -end -``` -Here, `<-->` is used to create a bi-directional reaction (with forward rate `kP` and backward rate `kD`). Next, the model (stored in the variable `rn`) can be used as input to various types of [simulations](@ref ref). - -## [Basic syntax](@id dsl_description_basic_syntax) -The basic syntax of the DSL is -```@example dsl_1 -rn = @reaction_network begin - 2.0, X --> Y - 1.0, Y --> X -end -``` -Here, you start with `@reaction_network begin`, next list all of the model's reaction, and finish with `end`. Each reaction consists of -- A *rate*. -- A (potentially empty) set of *substrates*. -- A (potentially empty) set of *products*. - -Each reaction line declares, in order, the rate, the substrate(s), and the product(s). The rate is separated from the substrate(s) by a `,`, and the substrate(s) from the production by a `-->` (other arrows, however, are [also possible](@ref ref)). In the above example, our model consists of two reactions. In the first one, `X` (the single substrate) becomes `Y` (the single product) at rate `2.0`. In the second reaction, `Y` becomes `X` at rate `1.0`. - -Finally, `rn = ` is used to store the model in the variable `rn` (a normal Julia variable, which does not need to be called `rn`). - -## [Defining parameters and species in the DSL](@id dsl_description_parameters_basics) -Typically, the rates are not constants, but rather parameters (which values can be set e.g. at [the beginning of each simulation](@ref ref)). To set parametric rates, simply use whichever symbol you wish to represent your parameter with. E.g. to set the above rates to `a` and `b`, we use: -```@example dsl_1 -rn1 = @reaction_network begin - a, X --> Y - b, Y --> X -end -``` - -Here we have used single-character symbols to designate all species and parameters. Multi-character symbols, however, are also permitted. E.g. we could call the rates `kX` and `kY`: -```@example dsl_1 -rn1 = @reaction_network begin - kX, X --> Y - kY, Y --> X -end -nothing # hide -``` -Generally, anything that is a [permitted Julia variable name](@id https://docs.julialang.org/en/v1/manual/variables/#man-allowed-variable-names) can be used to designate a species or parameter in Catalyst. - -## [Different types of reactions](@id dsl_description_reactions) - -### [Reactions with multiple substrates or products](@id dsl_description_reactions_multiples) -Previously, our reactions have had a single substrate and a single product. However, reactions with multiple substrates and/or products are possible. Here, all the substrates (or products) are listed and separated by a `+`. E.g. to create a model where `X` and `Y` bind (at rate `kB`) to form `XY` (which then can dissociate, at rate `kD`, to form `XY`) we use: -```@example dsl_1 -rn2 = @reaction_network begin - kB, X + Y --> XY - kD, XY --> X + Y -end -``` -Reactions can have any number of substrates and products, and their names do not need to have any relationship to each other, as demonstrated by the following mock model: -```@example dsl_1 -rn3 = @reaction_network begin - k, X + Y + Z --> A + B + C + D -end -``` - -### [Reactions with degradation or production](@id dsl_description_reactions_degradation_and_production) -Some reactions have no products, in which case the substrate(s) are degraded (i.e. removed from the system). To denote this, set the reaction's right-hand side to `0`. Similarly, some reactions have no substrates, in which case the product(s) are produced (i.e. added to the system). This is denoted by setting the left-hand side to `0`. E.g. to create a model where a single species `X` is both created (in the first reaction) and degraded (in a second reaction), we use: -```@example dsl_1 -rn4 = @reaction_network begin - p, 0 --> X - d, X --> 0 -end -``` - -### [Reactions with non-unitary stoichiometries](@id dsl_description_reactions_stoichiometries) -Reactions may include multiple copies of the same reactant (i.e. a substrate or a product). To specify this, the reactant is preceded by a number indicating its number of copies (also called the reactant's *stoichiometry*). E.g. to create a model where two copies of `X` dimerise to form `X2` (which then dissociate back to two `X` copies) we use: -```@example dsl_1 -rn5 = @reaction_network begin - kB, 2X --> X2 - kD, X2 --> 2X -end -``` -Reactants whose stoichiometries are not defined are assumed to have stoichiometry `1`. Any integer number can be used, furthermore, [decimal numbers and parameters can also be used as stoichiometries](@ref ref). A discussion of non-unitary (i.e. not equal to `1`) stoichiometries affecting the created model can be found [here](@ref ref). - -Stoichiometries can be combined with `()` to define them for multiple reactants. Here, the following (mock) model declares the same reaction twice, both with and without this notation: -```@example dsl_1 -rn6 = @reaction_network begin - k, 2X + 3(Y + 2Z) --> 5(V + W) - k, 2X + 3Y + 6Z --> 5V + 5W -end -nothing # hide -``` - -## [Bundling of similar reactions](@id dsl_description_reaction_bundling) - -### [Bi-directional (or reversible) reactions](@id dsl_description_reaction_bundling_reversible) -As is the case for the following two-state model: -```@example dsl_1 -rn7 = @reaction_network begin - k1, X1 --> X2 - k2, X2 --> X1 -end -``` -it is common that reactions occur in both directions (so-called *bi-directional* reactions). Here, it is possible to bundle the reactions into a single line by using the `<-->` arrow. When we do this, the rate term must include two separate rates (one for each direction, these are enclosed by a `()` and separated by a `,`). I.e. the two-state model can be declared using: -```@example dsl_1 -rn7 = @reaction_network begin - (k1,k2), X1 <--> X2 -end -``` -Here, the first rate (`k1`) denotes the *forward rate* and the second rate (`k2`) the *backwards rate*. - -Catalyst also permits writing pure backwards reactions. These use identical syntax to forward reactions, but with the `<--` arrow: -```@example dsl_1 -rn8 = @reaction_network begin - k, X <-- Y -end -``` -Here, the substrate(s) are on the right-hand side and the product(s) are on the left-hand side. Hence, the above model can be written identically using: -```@example dsl_1 -rn8 = @reaction_network begin - k, Y --> X -end -``` -Generally, using forward reactions is clearer than backwards ones, with the later generally not being used. - -### [Bundling similar reactions on a single line](@id dsl_description_reaction_bundling_similar) -There exist additional situations where models contain similar reactions (e.g. systems where all system components degrade at identical rates). Reactions which share either rates, substrates, or products can be bundled into a single line. Here, the parts which are different for the reactions are written using `(,)` (containing one separate expression for each reaction). E.g., let us consider the following model where species `X` and `Y` both degrade at the rate `d`: -```@example dsl_1 -rn8 = @reaction_network begin - d, X --> 0 - d, Y --> 0 -end -``` -These share both the rates (`d`) and product (`0`), however, the substrates are different (`X` and `Y`). Hence, the reactions can be bundled into a single line using the common rate and product expression while providing separate substrate expressions: -```@example dsl_1 -rn8 = @reaction_network begin - d, (X,Y) --> 0 -end -``` -This declaration of the model is identical to the previous one. Reactions can share any subset of the rate, substrate, and product expression (the cases where they share all, or none, however, do not make sense to use). I.e. if the two reactions also have different degradation rates: -```@example dsl_1 -rn9 = @reaction_network begin - dX, X --> 0 - dY, Y --> 0 -end -``` -This can be represented using: -```@example dsl_1 -rn9 = @reaction_network begin - (dX,dY), (X,Y) --> 0 -end -``` - -It is possible to use bundling for any number of reactions. E.g. in the following model we bundle the conversion of a species $X$ between its various forms (where all reactions use the same rate $k$): -```@example dsl_1 -rn10 = @reaction_network begin - k, (X0,X1,X2,X3) --> (X1,X2,X3,X4) -end -``` - -It is possible to combine bundling with bi-directional reactions. In this case, the rate is first split into the forward and backwards rates. These may then (or may not) indicate several rates. We exemplify this using the two following two (identical) networks, created with and without bundling. -```@example dsl_1 -rn11 = @reaction_network begin - kf, S --> P1 - kf, S --> P2 - kb_1, P1 --> S - kb_2, P2 --> S -end -``` -```@example dsl_1 -rn11 = @reaction_network begin - (kf, (kb_1, kb_2)), S <--> (P1,P2) -end -``` - -Like when we designated stoichiometries, reaction bundling can be applied very generally to create some truly complicated reactions: -```@example dsl_1 -rn12 = @reaction_network begin - ((pX, pY, pZ),d), (0, Y0, Z0) <--> (X, Y, Z1+Z2) -end -``` -However, like for the above model, bundling reactions too zealously can reduce (rather improve) a model's readability. - -## [Non-constant reaction rates](@id dsl_description_nonconstant_rates) -So far we have assumed that all reaction rates are constant (being either a number of a parameter). Non-constant rates that depend on one (or several) species are also possible. More generally, the rate can be any valid expression of parameters and species. - -Let us consider a model with an activator (`A`, which degraded at a constant rate) and a protein (`P`). The production rate of `P` depends both on `A` and a parameter (`kP`). We model this through: -```@example dsl_1 -rn_13 = @reaction_network begin - d, A --> 0 - kP*A, 0 --> P -end -``` -Here, `P`'s production rate will be reduced as `A` decays. We can [print the ODE this model produces with `Latexify`](@ref ref): -```@example dsl_1 -using Latexify -latexify(rn_13; form=:ode) -``` - -In this case, we can generate an equivalent model by instead adding `A` as both a substrate and a product to `P`'s production reaction: -```@example dsl_1 -rn_13_alt = @reaction_network begin - d, A --> 0 - kp, A --> A + P -end -``` -We can confirm that this generates the same ODE: -```@example dsl_1 -latexify(rn_13_alt; form=:ode) -``` -Here, while these models will generate identical ODE, SDE, and jump simulations, the chemical reaction network models themselves are not equivalent. Generally, as pointed out in the two notes below, using the second form is preferable. -!!! warn - While `rn_13` and `rn_13_alt` will generate equivalent simulations, for [jump simulations](@ref ref), the first model will have [reduced performance](@ref ref) (which generally are more performant when rates are constant). - -!!! danger - Catalyst automatically infers whether quantities appearing in the DSL are species or parameters (as described [here](@ref ref)). Generally, anything that does not appear as a reactant is inferred to be a parameter. This means that if you want to model a reaction activated by a species (e.g. `kp*A, 0 --> P`), but that species does not occur as a reactant, it will be interpreted as a parameter. This can be handled by [manually declaring the system species](@ref ref). A full example of how to do this for this example can be found [here](@ref ref). - -Above we used a simple example where the rate was the product of a species and a parameter. However, any valid Julia expression of parameters, species, and values can be used. E.g the following is a valid model: -```@example dsl_1 -rn_14 = @reaction_network begin - 2.0 + X^2, 0 --> X + Y - k1 + k2^k3, X --> ∅ - pi * X/(sqrt(2) + Y), Y → ∅ -end -``` - -### [Using functions in rates](@id dsl_description_nonconstant_rates_functions) -It is possible for the rate to contain Julia functions. These can either be functions from Julia's standard library: -```@example dsl_1 -rn_16 = @reaction_network begin - d, A --> 0 - kp*sqrt(A), 0 --> P -end -``` -or ones defined by the user: -```@example dsl_1 -custom_function(p1, p2, X) = (p1 + X) / (p2 + X) -rn_17 = @reaction_network begin - d, A --> 0 - custom_function(k1,k2,E), 0 --> P -end -``` - -### [Pre-defined functions](@id dsl_description_nonconstant_rates_available_functions) -Two functions frequently used within systems biology are the [*Michaelis-Menten*](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) and [*Hill*](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)) functions. These are pre-defined in Catalyst and can be called using `mm(X,v,K)` and `hill(X,v,K,n)`. E.g. a self-activation loop where `X` activates its own production through a Hill function can be created using: -```@example dsl_1 -custom_function(p1,p2,X) = (p1+E)/(p2+E) -rn_18 = @reaction_network begin - hill(X,v,K,n), 0 --> P - d, X --> 0 -end -``` - -Catalyst comes with the following predefined functions: -- The Michaelis-Menten function: $mm(X,v,K) = v * X/(X + K)$. -- The repressive Michaelis-Menten function: $mmr(X,v,K) = v * K/(X + K)$. -- The Hill function: $hill(X,v,K,n) = v * (X^n)/(X^n + K^n)$. -- The repressive Hill function: $hillr(X,v,K,n) = v * (K^n)/(X^n + K^n)$. -- The activating/repressive Hill function: $hillar(X,Y,v,K,n) = v * (X^n)/(X^n + Y^n + K^n)$. - -### [Time-dependant rates](@id dsl_description_nonconstant_rates_time) -Previously we have assumed that the rates are independent of the [time variable, $t$](@ref ref). However, time-dependent reactions are also possible. Here, simply use `t` to represent the time variable. E.g., to create a production/degradation model where the production rate decays as time progresses, we can use: -```@example dsl_1 -rn_14 = @reaction_network begin - kp/(1 + t), 0 --> P - d, P --> 0 -end -``` - -Like previously, `t` can be part of any valid expression. E.g. to create a reaction with a cyclic rate (e.g. to represent a [circadian system](https://en.wikipedia.org/wiki/Circadian_rhythm)) we can use: -```@example dsl_1 -rn_15 = @reaction_network begin - A*(sin(2π*f*t - ϕ)+1)/2, 0 --> P - d, P --> 0 -end -``` - -## [Non-standard stoichiometries](@id dsl_description_stoichiometries) - -### [Non-integer stoichiometries](@id dsl_description_stoichiometries_decimal) -Previously all stoichiometric constants have been integer numbers, however, decimal numbers are also permitted. Here we create a birth-death model where each production reaction produces 1.5 units of `X`: -```@example dsl_1 -rn_16 = @reaction_network begin - p, 0 --> 1.5X - d, X --> 0 -end -``` -It is also possible to have non-integer stoichiometric coefficients for substrates. However, in this case the [`combinatoric_ratelaw = false`](@ref ref) option must be used. We note that non-integer stoichiometric coefficients does not make sense in most fields, however, this features is available for use for models where it does make sense. - -### [Parametric stoichiometries](@id dsl_description_stoichiometries_parameters) -It is possible for stoichiometric coefficients to be parameters. E.g. here we create a generic polymerisation system where `n` copies of `X` binds to form `Xn`: -```@example dsl_1 -rn_17 = @reaction_network begin - (kB,kD), n*X <--> Xn -end -``` -Now we can designate the value of `n` through a parameter when we e.g. create an `ODEProblem`: -```@example dsl_1 -u0 = [:X => 5.0, :Xn => 1.0] -ps = [:kB => 1.0, :kD => 0.1, :n => 4] -oprob = ODEProblem(rn_17, u0, (0.0, 1.0), ps) -nothing # hide -``` - -## [Using special symbols](@id dsl_description_symbols) -Julia permits any Unicode characters to be used in variable names, thus Catalyst can use these as well. Below we describe some cases where this can be useful. No functionality is, however, tied to this. - -### [Using ∅ in degradation/production reactions](@id dsl_description_symbols_empty_set) -Previously, we described how `0` could be used to [create degradation or production reactions](@ref dsl_description_reactions_degradation_and_production). Catalyst permits the user to instead use the `∅` symbol. E.g. the production/degradation system can alternatively be written as: -```@example dsl_1 -rn4 = @reaction_network begin - p, ∅ --> X - d, X --> ∅ -end -``` - -### [Using special arrow symbols](@id dsl_description_symbols_arrows) -Catalyst uses `-->`, `<-->`, and `<--` to denote forward, bi-directional, and backwards reactions, respectively. Several unicode representations of these arrows are available. Here, -- `>`, `→`, `↣`, `↦`, `⇾`, `⟶`, `⟼`, `⥟`, `⥟`, `⇀`, and `⇁` can be used to represent forward reactions. -- `↔`, `⟷`, `⇄`, `⇆`, `⇌`, `⇋`, , and `⇔` can be used to represent bi-directional reactions. -- `<`, `←`, `↢`, `↤`, `⇽`, `⟵`, `⟻`, `⥚`, `⥞`, `↼`, , and `↽` can be used to represent backwards reactions. - -E.g. the production/degradation system can alternatively be written as: -```@example dsl_1 -rn4 = @reaction_network begin - p, ∅ → X - d, X → ∅ -end -``` - -### [Using special symbols to denote species or parameters](@id dsl_description_symbols_special) -A range of possible characters are available which can be incorporated into species and parameter names. This includes, but is not limited to: -- Greek letters (e.g `α`, `σ`, `τ`, and `Ω`). -- Superscript and subscript characters (to create e.g. `k₁`, `k₂`, `Xₐ`, and `Xᴾ`). -- Non-latin, non-greek, letters (e.g. `ä`, `Д`, `س`, and `א`). -- Other symbols (e.g. `£`, `ℂ`, `▲`, and `♠`). - -An example of how this can be used to create a neat-looking model can be found in [Schwall et al. (2021)](https://www.embopress.org/doi/full/10.15252/msb.20209832) where it was used to model a sigma factor V circuit in the bacteria *Bacillus subtilis*: -```@example dsl_1 -σᵛ_model = @reaction_network begin - v₀ + hill(σᵛ,v,K,n), ∅ → σᵛ + A - kdeg, (σᵛ, A, Aσᵛ) → ∅ - (kB,kD), A + σᵛ ↔ Aσᵛ - L, Aσᵛ → σᵛ -end -nothing # hide -``` - -This functionality can also be used to create less serious models: -rn_13 = @reaction_network begin - 🍦, 😢 --> 😃 -end - -It should be noted that the following symbols are *not permitted* to be used as species or parameter names: -- `pi` and `π` (used in Julia to denote [`3.1415926535897...`](https://en.wikipedia.org/wiki/Pi)). -- `ℯ` (used in Julia to denote [Euler's constant](https://en.wikipedia.org/wiki/Euler%27s_constant)). -- `t` (used to denote the [time variable](@ref dsl_description_nonconstant_rates_time)). -- `∅` ([used for production/degradation reactions](@ref dsl_description_symbols_empty_set)). -- `im` (used in Julia to represent [complex numbers](https://docs.julialang.org/en/v1/manual/complex-and-rational-numbers/#Complex-Numbers)). -- `nothing` (used in Julia to denote [nothing](https://docs.julialang.org/en/v1/base/constants/#Core.nothing)). -- `Γ` (used by Catalyst to represent [conserved quantities](@ref ref)). - diff --git a/docs/src/02_model_creation/02_dsl_advanced.md b/docs/src/02_model_creation/02_dsl_advanced.md deleted file mode 100644 index 7afb7a6427..0000000000 --- a/docs/src/02_model_creation/02_dsl_advanced.md +++ /dev/null @@ -1,626 +0,0 @@ -# [The Catalyst DSL - Advanced Features and Options](@id dsl_advanced_options) -Within the Catalyst DSL, each line can represent either *a reaction* or *an option*. The [previous DSL tutorial](@ref dsl_description) described how to create reactions. This one will focus on options. These are typically used to supplant a model with additional information. Examples include the declaration of initial condition/parameter default values, or the creation observables or events. - -All options designations begins with declaration starting with `@`, followed by its input. E.g. the `@observables` options allows for the generation of observables. Each option can only be used once within each use of `@reaction_network`. A full list of options can be found [here](@ref ref), with most (but not all) being described in more detail below. This tutorial will also describe some additional advanced DSL features that does not directly include using an option. - -As a first step, we import Catalyst (which is required to run the tutorial): -```@example dsl_advanced_1 -using Catalyst -``` - -## [Explicit specification of network species and parameters](@id dsl_advanced_options_declaring_species_and_parameters) -[Previously](@ref ref), we mentioned that the DSL automatically determines which symbols corresponds to species and which to parameters. This is done by designating everything that appear as either a substrate or a product as a species, and all remaining quantities as parameters (i.e. those only appearing within rates or [stoichiometric constants](@ref ref)). Sometimes, one might want to manually override this default behaviour for a given symbol. I.e. consider the following model, where the conversion of a protein `P` from its inactive form (`Pᵢ`) to its active form (`Pₐ`) is catalysed by an enzyme (`E`). Using the most natural description: -```@example dsl_advanced_1 -catalysis_sys = @reaction_network begin - k*E, Pᵢ --> Pₐ -end -``` -`E` (as well as `k`) will be considered a parameter, something we can confirm directly: -```@example dsl_advanced_1 -parameters(catalysis_sys) -``` -If we want `E` to be considered a species, we can designate this using the `@species` option: -```@example dsl_advanced_1 -catalysis_sys = @reaction_network begin - @species E(t) - k*E, Pᵢ --> Pₐ -end -parameters(catalysis_sys) -``` -!!! note - When declaring species using the `@species` option, the species symbol must be followed by `(t)`. The reason is that species are time-dependent variables, and this time-dependency must be explicitly specified ([designation of non-`t` dependant species is also possible](@ref ref)). - -Similarly, the `@parameters` option can be used to explicitly designate something as a parameter: -```@example dsl_advanced_1 -catalysis_sys = @reaction_network begin - @parameters k - k*E, Pᵢ --> Pₐ -end -``` -Here, while `k` is explicitly defined as a parameter, no information is provided about `E`. Hence, the default case will be used (setting `E` to a parameter). The `@species` and `@parameter` options can be used simultaneously (although a quantity cannot be declared *both* as a species and a parameter). They may be followed by a full list of all species/parameters, or just a subset. - -While designating something which would default to a parameter as a species is straightforward, the reverse (creating a parameter which occur as a substrate or product) is more involved. This is, however, possible, and described [here](@ref dsl_advanced_options_constant_species). - -Rather than listing all species/parameters on a single line after the options, a `begin ... end` block can be used (listing one species/parameter on each line). E.g. in the following example we use this notation to explicitly designate all species and parameters of the system: -```@example dsl_advanced_1 -catalysis_sys = @reaction_network begin - @species begin - E(t) - Pᵢ(t) - Pₐ(t) - end - @parameters begin - k - end - k*E, Pᵢ --> Pₐ -end -``` - -A side-effect of using the `@species` and `@parameter` options is that they specify *the order in which the species and parameters are stored*. I.e. lets check the order of the parameters in the parameters in the following dimerisation model: -```@example dsl_advanced_1 -dimerisation = @reaction_network begin - (p,d), 0 <--> X - (kB,kD), 2X <--> X2 -end -parameters(dimerisation) -``` -The default order is typically equal to the order with which the parameters (or species) are encountered in the DSL (this is, however, not guaranteed). If we specify the parameters using `@parameters`, the order used within the option is used instead: -```@example dsl_advanced_1 -dimerisation = @reaction_network begin - @parameters kB kD p d - (p,d), 0 <--> X - (kB,kD), 2X <--> X2 -end -parameters(dimerisation) -``` -!!! danger - Generally, Catalyst and the SciML ecosystem *does not* guarantee that parameter and species order is preserved throughout various operations on the model. Writing programs that depend on these orders is *strongly discouraged*. There are, however, some legacy packages which still depends on order (one example is provided [here](@ref ref)). In these situations, this might be useful. However, in these cases, it is recommended that the user is extra wary, and also checks the order manually. - -!!! note - The syntax of the `@species` and `@parameters` options is identical to that used by the `@species` and `@parameters` macros [used in programmatic modelling in Catalyst](@ref programmatic_CRN_construction) (for e.g. designating metadata or initial conditions). Hence, if one have learnt how to specify species/parameters using either approach, that knowledge can be transferred to the other one. - -Generally, there are four main reasons for specifying species/parameters using the `@species` and `@parameters` option: -1. To designate a quantity, that would otherwise have defaulted to a parameter, as a species. -2. To designate default values for parameters/species initial conditions (described [here](@ref dsl_advanced_options_default_vals)). -3. To designate metadata for species/parameters (described [here](@ref dsl_advanced_options_species_and_parameters_metadata)). -4. To designate a species or parameters that does not occur in reactions, but are still part of the model (e.g a [parametric initial condition](@ref dsl_advanced_options_parametric_initial_conditions)) - -!!!! warn - Catalyst's DSL automatically infer species and parameters from the input. However, it only does so for *quantities that appear in reactions*. Until now this has not been relevant. However, this tutorial will demonstrate cases where species/parameters that are not part of reactions are used. These *must* be designated using either the `@species` or `@parameters` options (or the `@variables` option, which is described [later](@ref dsl_advanced_options_variables)). - -### [Setting default values for species and parameters](@id dsl_advanced_options_default_vals) -When declaring species/parameters using the `@species` and `@parameters` options, one can also assign them default values (by appending them with `=` followed by the desired default value). E.g here we set `X`'s default initial condition value to $1.0$, and `p` and `d`'s default values to $1.0$ and $0.2$, respectively: -```@example dsl_advanced_defaults -using Catalyst # hide -rn = @reaction_network begin - @species X(t)=1.0 - @parameters p=1.0 d=0.1 - (p,d), 0 <--> X -end -``` -Next, if we simulate the model, we do not need to provide values for species or parameters that have default values. In this case all have default values, so both `u0` and `ps` can be empty vectors: -```@example dsl_advanced_defaults -using OrdinaryDiffEq, Plots -u0 = [] -tspan = (0.0, 10.0) -p = [] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob, Tsit5()) -plot(sol) -``` -It is still possible to provide values for some (or all) initial conditions/parameters in `u0`/`ps` (in which case these overrides the default values): -```@example dsl_advanced_defaults -u0 = [:X => 4.0] -p = [:d => 0.5] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob, Tsit5()) -plot(sol) -``` -It is also possible to declare a model with default values for only some initial conditions/parameters: -```@example dsl_advanced_defaults -using Catalyst # hide -rn = @reaction_network begin - @species X(t)=1.0 - (p,d), 0 <--> X -end - -tspan = (0.0, 10.0) -p = [:p => 1.0, :D => 0.2] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob, Tsit5()) -plot(sol) -``` -API for checking the default values of species and parameters can be found [here](@ref ref). - -### [Setting parametric initial conditions](@id dsl_advanced_options_parametric_initial_conditions) -In the previous section, we designated default values for species' initial conditions and parameters. However, the right-hand side of the designation accepts any valid expression (not only numeric values). While this can be used to set up some advanced default values, the most common use-case is to designate a species's initial condition as a parameter. E.e. in the following example we represent the initial condition of `X` using the parameter `X₀`. -```@example dsl_advanced_defaults -rn = @reaction_network begin - @species X(t)=X₀ - @parameters X₀ - (p,d), 0 <--> X -end -``` -Please note that as the parameter `X₀` does not occur as part of any reactions, Catalyst's DSL cannot infer whether it is a species or a parameter. This must hence be explicitly declared. We can now simulate our model while providing `X`'s value through the `X₀` parameter: -```@example dsl_advanced_defaults -u0 = [] -p = [:X₀ => 1.0, :p => 1.0, :d => 0.5] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob, Tsit5()) -plot(sol) -``` -It is still possible to designate $X$'s value in `u0`, in which case this overrides the default value. -```@example dsl_advanced_defaults -u0 = [:X => 0.5] -p = [:X₀ => 1.0, :p => 1.0, :d => 0.5] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob, Tsit5()) -plot(sol) -``` -Please note that `X₀` is still a parameter of the system, and as such its value must still be designated to simulate the model (even if it is not actually used). - -### [Designating metadata for species and parameters](@id dsl_advanced_options_species_and_parameters_metadata) -Catalyst permits the user to define *metadata* for species and parameters. This permits the user to assign additional information to these, which can be used for a variety of purposes. Some Catalyst features depend on using metadata (with each such case describing specifically how this is done). Here we will introduce how to set metadata, and describe some common metadata types. - -Whenever a species/parameter is declared using the `@species`/`@parameters` options, it can be followed by a `[]` within which the metadata is given. Each metadata entry consists of the metadata's name, followed by a `=`, followed by its value. E.g. the `description` metadata allows you to attach a [`String`](https://docs.julialang.org/en/v1/base/strings/) to a species/parameter. Here we create a simple model where we add descriptions to all species and parameters. -```@example dsl_advanced_metadata -using Catalyst # hide -two_state_system = @reaction_network begin - @species Xi(t) [description="The X's inactive form"] Xa(t) [description="The X's active form"] - @parameters kA [description="X's activation rate"] kD [description="X's deactivation rate"] - (ka,kD), Xi <--> Xa -end -``` -A metadata can be given to only a subset of a system's species/parameters, and a quantity can be given several metadata entries. To give several metadata, separate each by a `,`. Here we only provide a description for `kA`, for which we also provide a [bounds metadata](@ref https://docs.sciml.ai/ModelingToolkit/dev/basics/Variable_metadata/#Bounds), -```@example dsl_advanced_metadata -two_state_system = @reaction_network begin - @parameters kA [description="X's activation rate", bound=(0.01,10.0)] - (ka,kD), Xi <--> Xa -end -``` - -It is possible to add both default values and metadata to a parameter/species. In this case, first provide the default value, next the metadata. I.e. to in the above example set $kA$'s default value to $1.0$ we use -```@example dsl_advanced_metadata -two_state_system = @reaction_network begin - @parameters kA=1.0 [description="X's activation rate", bound=(0.01,10.0)] - (ka,kD), Xi <--> Xa -end -``` - -When designating metadata for species/parameters in `begin ... end` blocks the syntax changes slight. Here, a `,` must be inserted before the metadata (but after any potential default value). I.e. a version of the previous example can be written as -```@example dsl_advanced_metadata -two_state_system = @reaction_network begin - @parameters begin - kA, [description="X's activation rate", bound=(0.01,10.0)] - kD = 1.0, [description="X's deactivation rate"] - end - (ka,kD), Xi <--> Xa -end -``` - -Each metadata has its own getter functions. E.g. we can get the description of the parameter `kA` using `getdescription` (here we use [system indexing](@ref ref) to access the parameter): -```@example dsl_advanced_metadata -getdescription(two_state_system.kA) -``` - -It is not possible for the user to directly designate their own metadata. These have to first be added to Catalyst. Doing so is somewhat involved, and described in detail [here](@ref ref). A full list of metadata that can be used for species and/or parameters can be found [here](@ref ref). - -### [Designating constant-valued/fixed species parameters](@id dsl_advanced_options_constant_species) - -Catalyst enables the designation of parameters as `constantspecies`. These parameters can be used as species in reactions, however, their values are not changed by the reaction and remain constant throughout the simulation (unless changed by e.g. the [occurrence of an event]@ref ref). Practically, this is done by setting the parameter's `isconstantspecies` metadata to `true`. Here, we create a simple reaction where the species `X` is converted to `Xᴾ` at rate `k`. By designating `X` as a constant species parameter, we ensure that its quantity is unchanged by the occurrence of the reaction. -```@example dsl_advanced_constant_species -using Catalyst # hide -rn = @reaction_network begin - @parameters X [isconstantspecies=true] - k, X --> Xᴾ -end -``` -We can confirm that $Xᴾ$ is the only species of the system: -```@example dsl_advanced_constant_species -species(rn) -``` -Here, the produced model is actually identical to if $X$ had simply been a parameter in the reaction's rate: -```@example dsl_advanced_constant_species -rn = @reaction_network begin - k*X, 0 --> Xᴾ -end -``` - -A common use-case for constant species are when modelling systems where some species are present in such surplus that their amounts the reactions' effect on it is negligible. A system which is commonly modelled this way is the [Brusselator](https://en.wikipedia.org/wiki/Brusselator). - -### [Specifying non-species variables](@id dsl_advanced_options_variables) ---- MOVE THIS TO HYBRID SECTION --- -Chemical reaction network (CRN) models (which Catalyst creates) described how *species* are affected by the occurrence of *reaction events*. When they are converted to ODEs, the species become unknowns of the system. However, Catalyst permits the creation of [hybrid CRN models](@ref ref). These describe phenomenons which can only partially be modelled using CRNs. An example may be a bacterium. Here, we can use a CRN to model some internal system (e.g. controlling its growth rate). However, we might also want to model the bacterium's volume. Here, the volume cannot be considered a species (as it cannot plausibly be a reaction reactant). Instead, we should model it as a normal variable. Here, Catalyst provides the `@variables` option for adding non-species variables to the system. E.g. to create a model where a single growth factor ($G$) is produced and degraded, and where we also have a single volume variables ($V$) we can use: -```@example dsl_advanced_variables -using Catalyst # hide -rn = @reaction_network begin - @variables V(t) - (p,d), 0 <--> G -end -``` -Note that $V$ (like species) is time-dependant, and (like species) must be declared as such when the `@variables` option is used. We can now simulate our model (remembering to provide a value for $V$ as well as $G$): -```@example dsl_advanced_variables -using OrdinaryDiffEq, Plots -u0 = [:G => 0.1, :V => 1.0] -tspan = (0.0, 10.0) -p = [:p => 1.0, :d => 0.5] -oprob = ODEProblem(rn, u0, tspan, p) -sol = solve(oprob, Tsit5()) -plot(sol) -``` -Here, we have not actually described how $V$ interacting with our model, or how its value may change. Primarily variables are declared as part of hybrid CRN/equation modelling, which is described in more detail [here](@ref ref). - -You can set metadata and default initial condition values for variables using the same syntax as used for parameters and species. - -You can use the `variables` and `species` functions to retrieve a model's variables and species, respectively. The `unknowns` function retrieves both. - -## [Setting reaction metadata](@id dsl_advanced_options_reaction_metadata) ---- DISCUSS WHAT TO DO WITH THIS SECTION --- -Reactions can also have metadata. This is described in detail [here](@ref ref). - -## [Naming reaction networks](@id dsl_advanced_options_naming) -Each reaction network model has a name. It can be accessed using the `nameof` function. By default, some generic name is used: -```@example dsl_advanced_names -using Catalyst # hide -rn = @reaction_network begin - (p,d), 0 <--> X -end -nameof(rn) -``` -A specific name can be given as an argument between the `@reaction_network` and the `begin`. E.g. to name a network `my_network` we can use: -```@example dsl_advanced_names -rn = @reaction_network my_network begin - (p,d), 0 <--> X -end -nameof(rn) -``` - -A consequence of generic names being used by default is that networks, even if seemingly identical, by default are not. E.g. -```@example dsl_advanced_names -rn1 = @reaction_network begin - (p,d), 0 <--> X -end -rn2 = @reaction_network begin - (p,d), 0 <--> X -end -rn1 == rn2 -``` -The reason can be confirmed by checking that their respective (randomly generated) names are different: -```@example dsl_advanced_names -nameof(rn1) == nameof(rn2) -``` -By designating the networks to have the same name, however, identity is achieved. -```@example dsl_advanced_names -rn1 = @reaction_network my_network begin - (p,d), 0 <--> X -end -rn2 = @reaction_network my_network begin - (p,d), 0 <--> X -end -rn1 == rn2 -``` - -Setting model names is primarily useful for [hierarchical modelling](@ref ref), where network names are appended to the display name of subnetworks' species and parameters. - -## [Creating observables](@id dsl_advanced_options_observables) -Sometimes one might want to use observable variables. These are variables with values that can be computed directly from a system's state (rather than having their values implicitly given by reactions or equations). Observables can be designated using the `@observables` option. Here, the `@observables` option is followed by a `begin ... end` block with one line for each observable. Each line first gives the observable, followed by a `~` (*not* a `=`!), followed by an expression describing how to compute it. - -Let us consider a model where two species (`X` and `Y`) can bind to form a complex (`XY`, which also can dissociate back into `X` and `Y`). If we wish to create a representation for the total amount of `X` and `Y` in the system, we can do this by creating observables `Xtot` and `Ytot`: -```@example dsl_advanced_observables -using Catalyst # hide -rn = @reaction_network begin - @observables begin - Xtot ~ X + XY - Ytot ~ Y + XY - end - (kB,kD), X + Y <--> XY -end -``` -We can now simulate our model using normal syntax (initial condition values for observables should not, and can not, be provided): -```@example dsl_advanced_observables -using OrdinaryDiffEq -u0 = [:X => 1.0, :Y => 2.0, :XY => 0.0] -tspan = (0.0, 10.0) -ps = [:kB => 1.0, :kD => 1.5] -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(oprob, Tsit5()) -nothing # hide -``` - -Next, we can use [symbolic indexing](@ref simulation_structure_interfacing) of our solution object, but with the observable as input. E.g. we can use -```@example dsl_advanced_observables -sol[:Xtot] -``` -to get a vector with `Xtot`'s value throughout the simulation. We can also use -```@example dsl_advanced_observables -using Plots -plot(sol; idxs = [:Xtot, :Ytot]) -``` -to plot the observables (rather than the species). - -Observables can be defined using complicated expression containing species, parameters, and [variables](@ref ref) (but not other observables). In the following example (which uses a [parametric stoichiometry](@ref ref)) `X` polymerises to form a complex `Xn` containing `n` copies of `X`. Here, we create an observable describing the total number of `X` molecules in the system: -```@example dsl_advanced_observables -rn = @reaction_network begin - @observables Xtot ~ X + n*Xn - (kB,kD), n*X <--> Xn -end -nothing # hide -``` -!!! - If only a single observable is declared, the `begin .. end` block is not required and the observable can be declared directly after the `@observables` option. - -[Metadata](@ref dsl_advanced_options_species_and_parameters_metadata) can be supplied to an observable directly after the its declaration (but before its formula). If so, the metadata must be separated from the observable with a `,`, and the observable plus the metadata encapsulated by `()`. E.g. to add a [description metadata](@ref dsl_advanced_options_species_and_parameters_metadata) to our observable we can use -```@example dsl_advanced_observables -rn = @reaction_network begin - @observables (Xtot, [description="The total amount of X in the system."]) ~ X + n*Xn - (kB,kD), n*X <--> Xn -end -nothing # hide -``` - -Observables are by default considered [variables](@ref ref) (not species). To designate them as a species, they can be pre-declared using the `@species` option. I.e. Here `Xtot` becomes a species: -```@example dsl_advanced_observables -rn = @reaction_network begin - @species Xtot(t) - @observables (Xtot, [description="The total amount of X in the system."]) ~ X + n*XnXY - (kB,kD), n*X <--> Xn -end -nothing # hide -``` - -Some final notes regarding observables: -- The left-hand side of the observable declaration must contain a single symbol only (with the exception of metadata, which can also be supplied). -- All quantities appearing on the right-hand side must be declared elsewhere within the `@reaction_network` call (either by being part of a reaction, or through the `@species`, `@parameters`, or `@variables` options). -- Observables may not depend on other observables. -- Observables have their [dependent variable(s)](@ref ref) automatically assigned as the union of the dependent variables of the species and variables on which it depends. - -## [Creating events](@id dsl_advanced_options_events) ---- MOVE TO SEPARATE DOC PAGE --- -Sometimes one wishes to model events, describing things that can happen to it during a simulation. - - A chemical system where an amount of some species is added at a time point after the simulation's initiation. - - A simulation of a circadian rhythm, where light is turned on/off every 12 hours. - - A cell divides when some size variable reaches a certain threshold, halving the amount of each species in the system. - -Events are divided into *continuous* and *discrete* events, and these can be added directly to a system using the `continuous_events` and `discrete_events` options. Events can also be modelled through *callbacks*. These are different in that they are supplied in the simulation step (rather than on system creation), and generally provide more flexibility in how they may affect the system. Callbacks are described on a separate [page](@ref advanced_simulations_callbacks). - -The notation described below for creating continuous and discrete events is the same which is used in [ModelingToolkit to create events](https://docs.sciml.ai/ModelingToolkit/stable/basics/Events/), and which is used for [events for programmatic model creation](@ref ref). - -### [Continuous vs discrete events](@id dsl_advanced_options_events_continuous_vs_discrete) -Both continuous and discrete events combine some condition (for triggering the event) with some affect (describing their effects on the system). They differ in the following ways: -- They use slightly different notation. -- Discrete events' conditions are checked at *the end of* each simulation time step. For continuous events, the simulation instead finds the *exact time point* when the event is triggered at. -- Continuous events cannot be supplied to jump simulations. - -### [Continuous events](@id dsl_advanced_options_events_continuous) -Let us consider a simple system where species `X` degraded at a constant rate `d`. Next, we wish to add an event which adds `2.0` units of `X` whenever `X` reaches a critical threshold `1.0`. This can be done in the following manner: -```@example dsl_advanced_events -using Catalyst # hide -rn = @reaction_network begin - @continuous_events begin - X ~ 1.0 => [X ~ X + 2.0] - end - d, X --> 0 -end -nothing # hide -``` -Here, the `@continuous_events` option is followed by a `begin ... end` block. Next, each line corresponds to a separate event. Each event is created in the following manner: -- It combines a *condition* (denoting when the event will happen) with one (or several) *affects* (denoting what the event does to the system when it happens). -- The condition (left) and affect (right) are separated by a `=>`. -- The condition takes the form of an [equation](). Here, the event is triggered whenever the equation's two sides (separated by a `~`) are equal. -- The affect(s) are enclosed within `[]`. If there are multiple affects, these are separated by `,` (the example above contains a single affect). -- Each affect is a single equation that describes how a parameter, species, or [variable](@ref dsl_advanced_options_variables) is updated when the event is triggered. -- Each affect's equation's left-hand side must contain only the parameter/species/variable whose value should be updated. -- Each affect's equation's right-hand side is an expression describing its updated value. - -We can simulate the model we declared, just like any other model: -```@example dsl_advanced_events -using OrdinaryDiffEq, Plots -u0 = [:X => 2.0] -tspan = (0.0, 10.0) -ps = [:d => 1.0] - -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(ODEProblem, Tsit5()) -plot(sol) -``` -Inspecting the solution, we can confirm that whenever `X` reaches a value of `1.0`, `2.0` units of `X` is added to the system. - -In our example, we can also denote the critical quantities using parameters: -```@example dsl_advanced_events -rn = @reaction_network begin - @parameters X_thres X_add - @continuous_events begin - X ~ X_thres => [X ~ X + X_add] - end - d, X --> 0 -end -nothing # hide -``` -Here, since `X_thres` and `X_add` do not appear in any reactions, Catalyst cannot determine whether they are parameters or species. hence, they must be [explicitly designated as parameters by using the `@parameters` option](@ref dsl_advanced_options_declaring_species_and_parameters). Next, these can be designated as any value, and supplied to the `ODEProblem`: -```@example dsl_advanced_events -u0 = [:X => 2.0] -tspan = (0.0, 10.0) -ps = [:d => 1.0, :X_thres => 0.5, :X_add => 3.0] - -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(ODEProblem, Tsit5()) -plot(sol) -``` - -As previously noted, each continuous event can have multiple affects. The following system has two components (`X` and `Y`, one being produced and one being degraded). When their concentrations are equal, a continuous events reduce the concentration of `X` while increasing the concentration of `Y`: -```@example dsl_advanced_events -rn = @reaction_network begin - @continuous_events begin - X ~ Y => [X ~ X - 1.0, Y ~ Y + 1.0] - end - p, 0 --> X - d, Y --> 0 -end - -u0 = [:X => 1.0, :Y => 3.0] -tspan = (0.0, 10.0) -ps = [:p => 1.0, :d => 1.0] - -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(ODEProblem, Tsit5()) -plot(sol) -``` - -!!!warn - A single event (continuous or discrete) can update the value of either (one or several) species (and variables), or of (one or several) parameters. It is not possible for an event to update the values of both species/variables and parameters. - -In the above examples we have modelled a system with a single event. In these cases, the `begin end` block is not required, and the event can be provided on the same line as the `@continuous_events` option: -```@example dsl_advanced_events -rn = @reaction_network begin - @continuous_events X ~ Y => [X ~ X - 1.0, Y ~ Y + 1.0] - p, 0 --> X - d, Y --> 0 -end -nothing # hide -``` - -### [Discrete events](@id dsl_advanced_options_events_discrete) -Just like [continuous events](dsl_advanced_options_events_continuous), discrete events combine a condition with one or more affect statements. Here, discrete events' affects are created identically to those for continuous events. Discrete events' conditions are different. There exist 3 different types of discrete events, each with a different type of condition. All three types are created using the `@discrete_events` option, and a single system can contain a mix of all types. The three types are: -- Preset-time discrete events. -- Periodic discrete events. -- Conditional discrete events. - -#### [Preset-time discrete events](@id dsl_advanced_options_events_discrete_presettime) -*Present-time events* are events that happen at specific time points. Here, the condition is a vector with all the time points at which an event is triggered. E.g. here we create a production/degradation loop, where `2.0` units of `X` is added at time points `3.0` and `7.0` -```@example dsl_advanced_events -rn = @reaction_network begin - @discrete_events begin - [3.0, 7.0] => [X ~ X + 2.0] - end - (p,d), 0 <--> X -end - -u0 = [:X => 0.1, :Y => 3.0] -tspan = (0.0, 10.0) -ps = [:p => 1.0, :d => 0.5] - -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(ODEProblem, Tsit5()) -plot(sol) -``` - -The preset time points can also be parameters (in which case, they have to be [designated as such using the `@parameters` option](@ref dsl_advanced_options_declaring_species_and_parameters)): -```@example dsl_advanced_events -rn = @reaction_network begin - @parameters t1 t2 - @discrete_events begin - [t1, t2] => [X ~ X + 2.0] - end - (p,d), 0 <--> X -end -nothing -``` - -#### [Periodic discrete events](@id dsl_advanced_options_events_discrete_periodic) -When a discrete event's condition is a vector, a preset-time event is created. If it instead is a single value, a *periodic event* is created. These occur repeatedly throughout a simulation, with its period set by the affect term. E.g. here we create a system where `0.5` units of `X` is added every `1.0` time units. -```@example dsl_advanced_events -rn = @reaction_network begin - @discrete_events begin - 1.0 => [X ~ X + 0.5] - end - d, X --> 0 -end - -u0 = [:X => 1.0] -tspan = (0.0, 10.0) -ps = [:d => 1.0] - -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(ODEProblem, Tsit5()) -plot(sol) -``` -Like for preset-time events, periodic events' affects may contain parameters. - -#### [Conditional discrete events](@id dsl_advanced_options_events_discrete_conditional) -Finally, discrete events' condition may be a boolean expression (consisting of parameters, species, variables, and the time variable). Let's say that we want to create an event which, if the concentration of `X` is below a threshold `1.0`, adds `1.0` units of `X` to the system, then we can use the following discrete event: -```@example dsl_advanced_events -rn = @reaction_network begin - @discrete_events begin - X < 1.0 => [X ~ X + 2.0] - end - d, X --> 0 -end -``` -If we simulate the system using the same conditions as for our [similar, continuous, example](@ref dsl_advanced_options_events_continuous) we get the same result: -```@example dsl_advanced_events -u0 = [:X => 2.0] -tspan = (0.0, 10.0) -ps = [:d => 1.0] - -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(ODEProblem, Tsit5()) -plot(sol) -``` -So, how is modelling this event as a discrete or continuous event different? There are four differences: -1) For continuous events, the simulation method finds the exact time point when the condition triggers. Discrete events are triggered at the first time step when the condition holds. -2) This discrete event will be triggered whenever `X < 1.0` holds, not just when the concentration of `X` passes the threshold. E.g. it will be triggered if the initial concentration of `X` is less than `1.0`. -3) Only the discrete event event can be used with jump simulations. -4) The discrete event can be used to create more advanced conditions. - -E.g. using (4), we can modify our system so that the event is only triggered when time is less than `5.0` (after which `X` decays towards `0`): -```@example dsl_advanced_events -rn = @reaction_network begin - @discrete_events begin - (X < 1.0) & (t < 5.0) => [X ~ X + 2.0] - end - d, X --> 0 -end - -u0 = [:X => 2.0] -tspan = (0.0, 10.0) -ps = [:d => 1.0] - -oprob = ODEProblem(rn, u0, tspan, ps) -sol = solve(ODEProblem, Tsit5()) -plot(sol) -``` -!!!note - When we form composite boolean conditions for conditional discrete events, we use `&` to denote the AND operator (not `&&`, as this is currently not supported). - -!!!warn - Generally, discrete events including equality (`==`) will not be triggered. The reason is that the condition is only checked at the end of every time step. Hence, unless special precautions are taken to ensure that the simulator stops when the condition holds, it will unlikely be triggered. - - -## [Specifying non-time independent variables](@id dsl_advanced_options_ivs) - -As [described elsewhere](@ref ref), Catalyst's `ReactionSystem` models depends on a *time independent variable*, and potentially one or more *spatial independent variables*. By default, the independent variable `t` is used. We can declare another independent variable (which is automatically used as the default one) using the `@ivs` option. E.g. to use `τ` instead of `t` we can use -```@example dsl_advanced_ivs -using Catalyst # hide -rn = @reaction_network begin - @ivs τ - (ka,kD), Xi <--> Xa -end -nothing # hide -``` -We can confirm that `Xi` and `Xa` depend on `τ` (and not `t`): -```@example dsl_advanced_ivs -species(rn) -``` - -It is possible to designate several independent variables using `@ivs`. If so, the first one is considered the default (time) independent variable, while the following one(s) are considered spatial independent variable(s). If we want some species to depend on a non-default independent variable, this has to be explicitly declared: -```@example dsl_advanced_ivs -rn = @reaction_network begin - @ivs τ x - @species X(x) Y(x) - (p1,d1), 0 <--> X - (p2,d2), 0 <--> Y -end -species(rn) -``` -It is also possible to have species which depends on several independent variables: -```@example dsl_advanced_ivs -rn = @reaction_network begin - @ivs t x - @species Xi(t,x) Xa(t,x) - (ka,kD), Xi <--> Xa -end -species(rn) -``` - -!!! note - Setting spatial independent variables is primarily intended for modelling of spatial systems on continuous domains. Catalyst's support for this is currently under development. Hence, the utility of specifying spatial independent variables is currently limited. \ No newline at end of file