Skip to content

Commit

Permalink
Merge b97ee07 into bc6d339
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed Jun 8, 2024
2 parents bc6d339 + b97ee07 commit 33a1e4d
Show file tree
Hide file tree
Showing 11 changed files with 28 additions and 35 deletions.
13 changes: 4 additions & 9 deletions docs/src/model_creation/dsl_advanced.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# [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 supply a model with additional information. Examples include the declaration of initial condition/parameter default values, or the creation of observables.

All option designations begin with a declaration starting with `@`, followed by its input. E.g. the `@observables` option 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 do not involve using an option.
All option designations begin with a declaration starting with `@`, followed by its input. E.g. the `@observables` option allows for the generation of observables. Each option can only be used once within each use of `@reaction_network`. This tutorial will also describe some additional advanced DSL features that do not involve using an option.

As a first step, we import Catalyst (which is required to run the tutorial):
```@example dsl_advanced_explicit_definitions
Expand Down Expand Up @@ -130,7 +130,6 @@ oprob = ODEProblem(rn, u0, tspan, p)
sol = solve(oprob)
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 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.g. in the following example we represent the initial condition of `X` using the parameter `X₀`.
Expand Down Expand Up @@ -198,13 +197,11 @@ two_state_system = @reaction_network begin
end
```

Each metadata has its own getter functions. E.g. we can get the description of the parameter `kA` using `ModelingToolkit.getdescription` (here we use [system indexing](@ref ref) to access the parameter):
Each metadata has its own getter functions. E.g. we can get the description of the parameter `kA` using `ModelingToolkit.getdescription`:
```@example dsl_advanced_metadata
ModelingToolkit.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 constraint_equations_events). 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.
Expand Down Expand Up @@ -391,11 +388,11 @@ 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.
- Observables have their dependent variable(s) automatically assigned as the union of the dependent variables of the species and variables on which it depends.

## [Specifying non-time independent variables](@id dsl_advanced_options_ivs)

As [described elsewhere](@ref ref), Catalyst's `ReactionSystem` models depend 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
Catalyst's `ReactionSystem` models depend 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
Expand Down Expand Up @@ -466,5 +463,3 @@ A reaction's metadata can be accessed using specific functions, e.g. `Catalyst.h
rx = @reaction p, 0 --> X, [description="A production reaction"]
Catalyst.getdescription(rx)
```

A list of all available reaction metadata can be found [here](@ref ref).
10 changes: 5 additions & 5 deletions docs/src/model_creation/dsl_basics.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# [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 dsl_advanced_options) 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 introduction_to_catalyst) and [following](@ref simulation_intro) tutorials describe how to simulate models once they have been created using the DSL. This tutorial will solely focus on model creation.
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)). [Previous](@ref introduction_to_catalyst) and [following](@ref simulation_intro) 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_basics_intro
Expand Down Expand Up @@ -217,7 +217,7 @@ 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, the first model will have [reduced performance](@ref ref) (which generally are more performant when rates are constant).
While `rn_13` and `rn_13_alt` will generate equivalent simulations, for jump simulations, the first model will have reduced performance (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 dsl_advanced_options_declaring_species_and_parameters)). 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 dsl_advanced_options_declaring_species_and_parameters).
Expand Down Expand Up @@ -282,7 +282,7 @@ end
```

!!! warn
Jump simulations cannot be performed for models with time-dependent rates without additional considerations, which are discussed [here](@ref ref).
Jump simulations cannot be performed for models with time-dependent rates without additional considerations.

## [Non-standard stoichiometries](@id dsl_description_stoichiometries)

Expand All @@ -294,7 +294,7 @@ rn_16 = @reaction_network begin
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 do not make sense in most fields, however, this feature is available for use for models where it does make sense.
It is also possible to have non-integer stoichiometric coefficients for substrates. However, in this case the `combinatoric_ratelaw = false` option must be used. We note that non-integer stoichiometric coefficients do not make sense in most fields, however, this feature 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` bind to form `Xn`:
Expand Down Expand Up @@ -368,4 +368,4 @@ It should be noted that the following symbols are *not permitted* to be used as
- `` ([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)).
- `Γ` (used by Catalyst to represent conserved quantities).
2 changes: 1 addition & 1 deletion docs/src/model_creation/examples/basic_CRN_library.md
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ oplt = plot(osol; idxs = X₁ + X₂, title = "Reaction rate equation (ODE)")
splt = plot(ssol; idxs = X₁ + X₂, title = "Chemical Langevin equation (SDE)")
plot(oplt, splt; lw = 3, ylimit = (99,101), size = (800,450), layout = (2,1))
```
Catalyst has special methods for working with conserved quantities, which are described [here](@ref ref).
Catalyst has special methods for working with conserved quantities, which are described here.

## [Michaelis-Menten enzyme kinetics](@id basic_CRN_library_mm)
[Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions, it can be simplified to a single function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ plot!(sol_n10; idxs = :X10, label = "n = 10")
## [Modelling linear pathways using programmatic, generative, modelling](@id programmatic_generative_linear_pathway_generative)
Above, we investigated the impact of linear pathways' lengths on their behaviours. Since the models were implemented using the DSL, we had to implement a new model for each pathway (in each case writing out all reactions). Here, we will instead show how [programmatic modelling](@ref programmatic_CRN_construction) can be used to generate pathways of arbitrary lengths.

First, we create a function, `generate_lp`, which creates a linear pathway model of length `n`. It utilises [*vector variables*](@ref ref) to create an arbitrary number of species, and also creates an [observable](@ref ref) for the final species of the chain.
First, we create a function, `generate_lp`, which creates a linear pathway model of length `n`. It utilises *vector variables* to create an arbitrary number of species, and also creates an observable for the final species of the chain.
```@example programmatic_generative_linear_pathway_generative
using Catalyst # hide
t = default_t()
Expand Down
2 changes: 1 addition & 1 deletion docs/src/model_creation/model_visualisation.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ Here, we note that the output of `latexify(brusselator)` is identical to how a m
latexify(brusselator; form = :ode)
```
!!! note
Internally, `latexify(brusselator; form = :ode)` calls `latexify(convert(ODESystem, brusselator))`. Hence, if you have already [generated the `ODESystem` corresponding to your model](@ref ref), it can be used directly as input to `latexify`.
Internally, `latexify(brusselator; form = :ode)` calls `latexify(convert(ODESystem, brusselator))`. Hence, if you have already generated the `ODESystem` corresponding to your model, it can be used directly as input to `latexify`.

!!! note
It should be possible to also generate SDEs through the `form = :sde` input. This feature is, however, currently broken.
Expand Down
4 changes: 2 additions & 2 deletions docs/src/model_simulation/ode_simulation_performance.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# [Advice for performant ODE simulations](@id ode_simulation_performance)
We have previously described how to perform ODE simulations of *chemical reaction network* (CRN) models. These simulations are typically fast and require little additional consideration. However, when a model is simulated many times (e.g. as a part of solving an [inverse problem](@ref ref)), or is very large, simulation run
We have previously described how to perform ODE simulations of *chemical reaction network* (CRN) models. These simulations are typically fast and require little additional consideration. However, when a model is simulated many times (e.g. as a part of solving an inverse problem), or is very large, simulation run
times may become noticeable. Here we will give some advice on how to improve performance for these cases.

Generally, there are few good ways to, before a simulation, determine the best options. Hence, while we below provide several options, if you face an application for which reducing run time is critical (e.g. if you need to simulate the same ODE many times), it might be required to manually trial these various options to see which yields the best performance ([BenchmarkTools.jl's](https://github.com/JuliaCI/BenchmarkTools.jl) `@btime` macro is useful for this purpose). It should be noted that the default options typically perform well, and it is primarily for large models where investigating alternative options is worthwhile. All ODE simulations of Catalyst models are performed using the OrdinaryDiffEq.jl package, [which documentation](https://docs.sciml.ai/DiffEqDocs/stable/) provides additional advice on performance.
Expand Down Expand Up @@ -304,7 +304,7 @@ nothing # hide
```
Here have we increased the number of simulations to 10,000, since this is a more appropriate number for GPU parallelisation (as compared to the 100 simulations we performed in our CPU example).
!!! note
Currently, declaration of static vectors requires [symbolic, rather than symbol, form](@ref ref) for species and parameters. Hence, we here first [`@unpack` these](@ref ref) before constructing `u0` and `ps` using `@SVector`.
Currently, declaration of static vectors requires symbolic, rather than symbol, form for species and parameters. Hence, we here first `@unpack` these before constructing `u0` and `ps` using `@SVector`.

We can now simulate our model using a GPU-based ensemble algorithm. Currently, two such algorithms are available, `EnsembleGPUArray` and `EnsembleGPUKernel`. Their differences are that
- Only `EnsembleGPUKernel` requires arrays to be static arrays (although it is still advantageous for `EnsembleGPUArray`).
Expand Down
Loading

0 comments on commit 33a1e4d

Please sign in to comment.