[WIP] Add new parts to coupled equations docs

docs/src/faqs.md

@@ -349,4 +349,39 @@ rn = @reaction_network begin
     @parameters k1 k2
     (k1, k2), A <--> B
 end
+```
+
+## [How do I designating custom differential operators when creating coupled differential equation models](@id faq_custom_differentials)
+When [coupling differential equations to reaction network models](@ref coupled_models_dsl), by default, `D` designates the differential operator:
+```@example faq_custom_differentials
+using Catalyst
+rn = @reaction_network begin
+    @equations D(V) ~ X - λ * V
+    (p*V,d), 0 <--> X
+end
+```
+However, it is possible to specify different symbols as differential operators. Potential reasons include that the symbol `D` is already used for a different purpose, a different differential symbol is preferred, or that multiple differentials with respect to multiple independent variables are needed.
+
+Non-default differentials can be declared using the `@differentials` option. E.g. here we declare the previous model, but call our differential `Δ` instead of `D`:
+```@example faq_custom_differentials
+rn = @reaction_network begin
+    @differentials Δ = Differential(t)
+    @equations Δ(V) ~ X - λ * V
+    (p*V,d), 0 <--> X
+end
+```
+Here, the left-hand side of the differential contains the differential's name only. The right-hand side is `Differential(t)`, where `t` is the default [time independent variable](@ref ref) used within Catalyst. If you wish to declare a differential with respect to another independent variable, replace `t` with its symbol. If you need to declare multiple differentials, this can be done by providing a `begin ... end` block (with one differential declaration on each line) to the `@differential` option.
+
+In [programmatic modelling](@ref programmatic_CRN_construction), we typically assign the default time differential to the symbol `D`:
+```@example faq_custom_differentials
+D = default_time_deriv()
+```
+However, we can assign it to any other symbol:
+```@example faq_custom_differentials
+Δ = default_time_deriv()
+```
+or declare it manually using a similar syntax as in the DSL:
+```@example faq_custom_differentials
+t = default_t()
+Dt = Differential(t)
 ```

docs/src/model_creation/coupled_non_crn_models.md

@@ -166,4 +166,136 @@ equations. To solve and plot the model we proceed like normal
 oprob = ODEProblem(growing_cell, [], (0.0, 1.0))
 sol = solve(oprob, Tsit5())
 plot(sol)
-```
+```
+
+## [Coupled reaction networks with differential equations](@id coupled_models_diffeqs)
+When declaring models using the DSL, differential equations can be inserted directly into the model through the `@equations` option. I.e. in the following model we have a protein in an inactive and active form ($Xᵢ$ and $Xₐ$, respectively). The conversion to the active form is promoted by the presence of a nutrient ($N$). If the nutrient is then degraded linearly with the concentration of $Xₐ$, then this can be described by the differential equation
+```math
+\mathbf{\frac{dN}{dt} = - d Xₐ(t) N(t)}.
+```
+This is declared with the following model:
+```@example coupled_diff_eqs
+using Catalyst
+rs = @reaction_network begin
+    @parameters d
+    @equations D(N) ~ -d*Xₐ*N
+    (kₐ*N, kᵢ), Xᵢ <--> Xₐ
+end
+```
+Here we note that:
+- Symbols not occurring as species within the `@equations` option are by default assumed to be [variables](@ref coupled_models_variables). Hence, `d` must [explicitly be declared as a parameter](@ref dsl_advanced_options_declaring_species_and_parameters).
+- Whenever the `D(X)` notation occurs within the `@equations` option, this is by default assumed to be the differential with respect to the default time independent variable (i.e. $dx/dt$). This means you have to *take extra care* if using `D` as e.g. a species name (typically by [declaring a custom symbol as a differential operator](@ref faq_custom_differentials)).
+
+The model can be simulated and plotted using normal syntax, providing an initial condition for `N` and a value for `d` through the normal initial condition and parameter value vectors.
+```@example coupled_diff_eqs
+using OrdinaryDiffEqDefault, Plots
+u0 = [:Xᵢ => 1.0, :Xₐ => 0.0, :N => 1.0]
+ps = [:kₐ => 2.0, :kᵢ => 0.5, :d => 1.0]
+oprob = ODEProblem(rs, u0, 10.0, ps)
+sol = solve(oprob)
+plot(sol)
+```
+
+!!! note
+    Species that occur within reactions cannot be subject to a differential operator within a `@equations` block.
+!!! note
+    When providing multiple equations to a model, write each equation on a separate line within a `@equations begin ... end` block.
+!!! note
+    Generally, the rules for declaring equations are the same as those used within the [ModelingToolkitBase.jl](https://github.com/SciML/ModelingToolkit.jl) package.
+
+## [Coupled reaction networks with algebraic equations](@id coupled_models_algeqs)
+Catalyst also permits the coupling of *algebraic equations* (equations not containing differentials). In practice, these are handled similarly to how differential equations are handled, but with the following differences:
+- The `mtkcompile = true` option must be provided to any `ODEProblem` (or other problem) generated from a model containing algebraic equations.
+- For variables involved in algebraic equations, a *guess* is provided, rather than an initial condition.
+
+A special case of algebraic equations (where a new variable is trivially produced by an algebraic expression) is *observables*. These are described in a separate section [here](@ref dsl_advanced_options_observables).
+
+We will demonstrate using a system where a species $X$ is produced and degraded, and can also bind/dissociate to form a dimer $X2$. Here, if the binding/unbinding dynamics is much faster than the production/degradation dynamics, $X2$ can be eliminated through an algebraic equation. This creates the following model:
+```@example coupled_eqs_alg_eq
+using Catalyst
+algebraic_crn = @reaction_network begin
+    @parameters kB kD
+    @equations kB*X^2 ~ kD*X2
+    (mmr(X2,v,K),d), 0 <--> X
+end
+```
+where we also made the production rate of $X$ depend on $X2$ through a Michaelis-Menten function.
+
+We can now simulate the model. Here, at each time step, $X2$'s value can be computed from $X$'s. Similarly, for the initial condition, if we know $X$, we can compute $X2$. Hence, we will not provide an initial condition value for $X2$. However, we will have to provide an initial *guess* for $X2$'s value, from which an internal [nonlinear solve-call](@ref steady_state_solving_nonlinear) will be used to compute the full ODE simulation initial condition.
+```@example coupled_eqs_alg_eq
+u0 = [:X => 1.0]
+ps = [:v => 2.0, :K => 1.0, :d => 0.2, :kB => 0.1, :kD => 0.4]
+guesses = [:X2 => 1.0]
+nothing # hide
+```
+Next, we provide `guesses` to our `ODEProblem` as an additional argument. Furthermore, we will use the `mtkcompile = true` argument, which is always required when simulating models containing algebraic equations. With these modifications, the model can be simulated using standard workflows.
+```@example coupled_eqs_alg_eq
+using OrdinaryDiffEqDefault, Plots
+oprob = ODEProblem(algebraic_crn, u0, 10.0, ps; mtkcompile = true, guesses)
+sol = solve(oprob)
+plot(sol)
+```
+
+There is no requirement on which values are provided as guesses and which as initial conditions. I.e. if we know the value of $X2$ we can provide this as the initial condition while instead providing a guess for $X$'s value.
+```@example coupled_eqs_alg_eq
+u0 = [:X2 => 1.0]
+guesses = [:X => 1.0]
+oprob = ODEProblem(algebraic_crn, u0, 10.0, ps; mtkcompile = true, guesses)
+nothing # hide
+```
+
+There is one exception to this, which is if the algebraic equation is formatted such that a variable is isolated on the equation left-hand side. I.e. if we were to declare the model using
+```@example coupled_eqs_alg_eq
+algebraic_crn_alt = @reaction_network begin
+    @parameters kB kD
+    @equations X2 ~ kB*X^2/kD
+    (mmr(X2,v,K),d), 0 <--> X
+end
+```
+Here, $X2$ can be fully eliminated from the system when the equations are generated, and we neither need to provide an initial condition nor a guess for it:
+```@example coupled_eqs_alg_eq
+u0 = [:X => 1.0]
+ps = [:v => 2.0, :K => 1.0, :d => 0.2, :kB => 0.1, :kD => 0.4]
+oprob = ODEProblem(algebraic_crn_alt, u0, 10.0, ps; mtkcompile = true)
+sol = solve(oprob)
+nothing # hide
+```
+A side effect of this is that $X2$ by default is not plotted with the solution, something which must be explicitly requested using the [`idxs` argument](@ref simulation_plotting_options):
+```@example coupled_eqs_alg_eq
+plot(sol; idxs = [:X, :X2])
+```
+In practice, we can check the equations that are generated by manually converting them to ODEs through `ode_model`:
+```@example coupled_eqs_alg_eq
+ode_model(algebraic_crn)
+```
+```@example coupled_eqs_alg_eq
+ode_model(algebraic_crn_alt)
+```
+
+## [Notes on *species* and *variables*](@id coupled_models_variables)
+Throughout this tutorial we have distinguished between *species* and *variables*. Functionally, both these are *unknown quantities* that we are simulating or solving for in `solve` calls. The differences are that
+- Only *species* can participate as reactants in reactions.
+- Only *variables* can be the subject of differentials in equations.
+
+Furthermore, species and variables are declared using different syntax. Either when declaring them [programmatically](@ref programmatic_CRN_construction), or through DSL options, we use `@species` and `@variables`, respectively. I.e. if we want to designate a [default initial condition value](@ref dsl_advanced_options_default_vals) for $N$ in the [nutrient model](@ref coupled_models_diffeqs) we use
+```@example coupled_eqs_variables
+using Catalyst # hide
+rs = @reaction_network begin
+    @parameters d
+    @variables N(t) = 1.0
+    @equations D(N) ~ -d*Xₐ*N
+    (kₐ*N, kᵢ), Xᵢ <--> Xₐ
+end
+nothing # hide
+```
+
+To retrieve either a model's species, variables, or all unknowns, use the `species`, `nonspecies`, and `unknowns` functions:
+```@example coupled_eqs_variables
+species(rs)
+```
+```@example coupled_eqs_variables
+nonspecies(rs)
+```
+```@example coupled_eqs_variables
+unknowns(rs)
+```