|
| 1 | +# [Coupling chemical reaction networks with differential or algebraic equations](@id coupled_equations) |
| 2 | +Catalyst is primarily designed to model [chemical reaction networks (CRNs)](@ref ref), which are defined by a set of species and reaction events. These can then be used to generate e.g. ODEs which can be simulated. However, there are a large number phenomena, many relevant to the fields utilising CRNs, that cannot be described by this formalism. To model these, more general equations (either [differential](https://en.wikipedia.org/wiki/Differential_equation) or [algebraic](https://en.wikipedia.org/wiki/Algebraic_equation)) are needed. Here, Catalyst provides the functionality of coupling CRN models with algebraic or differential equations models. For creating pure non-CRN models, we recommend using [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl). |
| 3 | + |
| 4 | +## [Basic example](@id coupled_equations_basic_example) |
| 5 | +Let us consider a simple model of a growth factor ($G$), which controls a cell's volume ($V$). Here, we can model the production and degradation of the growth factor using a simple CRN in Catalyst: |
| 6 | +```@example coupled_eqs_1 |
| 7 | +using Catalyst |
| 8 | +crn = @reaction_network begin |
| 9 | + (p,d), 0 <--> G |
| 10 | +end |
| 11 | +``` |
| 12 | +CRNs, however, model quantities that exists in (potentially very large) discrete quantities (like number of protein molecules) and change in discrete reaction events (like production and degradation events). This does not hold for a cell's volume. Here we will instead model the cell volume using the following differential equation: |
| 13 | +```math |
| 14 | +\frac{d\mathbf{V}}{dt} = \mathbf{G(t)} - \mathbf{G(V)}, |
| 15 | +``` |
| 16 | +In this case, we can insert this differential equation directly into our model using the `@equations` option, and `D(V)` to describe the differential with respect to $V$: |
| 17 | +```@example coupled_eqs_1 |
| 18 | +coupled_crn = @reaction_network begin |
| 19 | + @equations begin |
| 20 | + D(V) ~ G - V |
| 21 | + end |
| 22 | + (p,d), 0 <--> G |
| 23 | +end |
| 24 | +nothing # hide |
| 25 | +``` |
| 26 | +We can check the ODE this model would generate [using Latexify](@ref ref): |
| 27 | +```@example coupled_eqs_1 |
| 28 | +using Latexify |
| 29 | +latexify(coupled_crn; form=:ode) |
| 30 | +``` |
| 31 | +Finally, the model can be simualted using standard syntax (while providing initial conditions for $V$ as well as for $G$): |
| 32 | +```@example coupled_eqs_1 |
| 33 | +using OrdinaryDiffEq, Plots |
| 34 | +u0 = [:G => 0.1, :V => 1.0] |
| 35 | +tspan = (0.0, 10.0) |
| 36 | +ps = [:p => 1.0, :D => 0.2] |
| 37 | +oprob = ODEProblem(coupled_crn, u0, tspan, ps) |
| 38 | +sol = solve(oprob) |
| 39 | +plot(sol) |
| 40 | +``` |
| 41 | + |
| 42 | +## Variables and species |
| 43 | +The unknowns of CRN models are called *species*. These are unknown variables of a system which can partake in reaction events. Catalyst support another form of unknowns, *variables*. While variables (like species) may have values that change continuously throughout a simulation, they cannot be substrates or products of a reaction. They can, however, be the subject of differential equations. Species, *cannot* be part of differential equations. |
| 44 | + |
| 45 | +## [Coupling CRNs with differential equation](@id coupled_equations_differential_equations) |
| 46 | +Like [demonstrated in the previous example](@ref coupled_equations_basic_example), differential equations can be added to CRN models through the `@equations` option (in the DSL, non-DSL approaches for coupled CRN/equations [are described later](@ref ref)). These are then saved as part of the generated `ReactionSystem` model, and added to the system of ODEs generated from the CRN. The `@equations` option is followed by a `begin ... end` block, where each line correspond to a single equation. Here we create a simple growth model which describes the growth of a cell's volume ($V$), which depends on a growth factor ($G$), and also depletes a nutrient supply ($N$): |
| 47 | +```@example coupled_eqs_2 |
| 48 | +growth_model = @reaction_network begin |
| 49 | + @equations begin |
| 50 | + D(V) ~ G - V |
| 51 | + D(N) ~ -V |
| 52 | + end |
| 53 | + (N*p,d), 0 <--> G |
| 54 | +end |
| 55 | +latexify(growth_model; form=:ode) |
| 56 | +``` |
| 57 | +Here, while $G$ is inferred to be a species, $V$ and $N$ are inferred to be variables. We can check this using the `species` and `variables` functions: |
| 58 | +```@example coupled_eqs_2 |
| 59 | +species(growth_model) |
| 60 | +``` |
| 61 | +```@example coupled_eqs_2 |
| 62 | +variables(growth_model) |
| 63 | +``` |
| 64 | +Previously, we have described how Catalyst automatically infers what is a species and what is a parameter depending on where in the model it occurs. Catalyst can also infer variables, however, the only situation where this is done is *terms that occur as the subject of the default differential on the left-hand side of a differential equation* E.g. above both $D$ and $V$ are inferred to be variables.Catalyst does not infer any other components to be species or variables. This means that if we want to include parameters in our differential equations, these must explicitly be declared as such: |
| 65 | +```@example coupled_eqs_2 |
| 66 | +growth_model = @reaction_network begin |
| 67 | + @parameters g d |
| 68 | + @equations begin |
| 69 | + D(V) ~ g*G - V |
| 70 | + D(N) ~ -d*V |
| 71 | + end |
| 72 | + (N*p,d), 0 <--> G |
| 73 | +end |
| 74 | +parameters(growth_model) |
| 75 | +``` |
| 76 | +Here, with the following exceptions: |
| 77 | +- Variables can neither be substrates or products of reactions. |
| 78 | +- Species cannot be subject of differentials in differential equations. |
| 79 | +- Neither species nor variables can occur as stoichiometric coefficients within reactions |
| 80 | + |
| 81 | +species and variables may occur throughout the model. E.g. in our example above, the production rate of $G$ depends on $N$, while the growth grate of $V$ depends on $G$. |
| 82 | + |
| 83 | +Generally, differential equations are formed similarly as in the [ModelingToolkit.jl package](https://github.com/SciML/ModelingToolkit.jl) (on which [Catalyst is built](@ref ref)). In practise, they follow exactly the same rules as when forming expressions for [rates in Catalyst](@ref ref). The only noteworthy exceptions are: |
| 84 | +- The `~` (not `=`!) is used to separate the left-hand and right-hand sides of the equations. |
| 85 | +- By default, `D(...)` is used to represent the differential operator with respect to time (or the [default independent variable](@ref ref)). |
| 86 | + |
| 87 | +Differential equations does not need to be written in a pure `D(...) ~ ...` form (with the differential isolated on the left-hand side). However, if this is not the case, `D` is not automatically inferred to be the differential operator. Here, a differential must be [manually declared](@ref ref). Higher-order differentials are also permitted, how to model these are described [here](@ref ref). In these cases, what is a variable is not inferred at all, and this must hence be declared explicitly. |
| 88 | + |
| 89 | +## [Coupling CRNs with algebraic equation](@id coupled_equations_algebraic_equations) |
| 90 | +Catalyst also permits the coupling of CRN models with *algebraic equations* (equations not containing differentials). In practise, these are handled similarly to how differential equations are handled, but with the following differences: |
| 91 | +- The equations does not contain any differentials. |
| 92 | +- The `structural_simplify = true` option must be provided to any `ODEProblem` (or other problem) generated from a problem containing an algebraic equation. |
| 93 | +- Catalyst will not infer what is a variable |
| 94 | + |
| 95 | +A special case of algebraic equations (where a new variable is trivially produced by an algebraic expression) are *observables*. These are described in a separate section [here](@ref ref). |
| 96 | + |
| 97 | +## [Converting coupled crn/equations to non-ODE systems](@id coupled_equations_conversions) |
| 98 | +As [previously described](@ref ref), Catalyst permits the conversion of `ReactionSystem`s to a range of model types. WIth the exception of `JumpProblem`s, all these can be generated from CRNs coupled to equation models. Below, we briefly describe each case. |
| 99 | + |
| 100 | +### [Converting coupled models to SteadyState problems](@id coupled_equations_conversions_steady_state) |
| 101 | + |
| 102 | +### [Converting coupled models to NonlinearProblems](@id coupled_equations_conversions_nonlinear) |
| 103 | + |
| 104 | +### [Converting coupled models to SDEs](@id coupled_equations_conversions_SDEs) |
| 105 | +Generally, Catalyst models containing variables (whenever these values' are given through differential/algebraic equations, or remain constant) can be converted to SDEs. Here, while the chemical Langevin equation is sued to generate noise for all species, equations for variables will remain deterministic. However, if the variables depend on species values, noise from these may still filter through. |
| 106 | + |
| 107 | +!!! note |
| 108 | + Adding support for directly adding noise to variables during SDE conversion is a work in progress. If this is a feature you are interested in, please raise an issue. |
| 109 | + |
| 110 | +### [Other applications of coupled CRN/equation models](@id coupled_equations_conversions_extensions) |
| 111 | +In other sections, we describe how to perform [bifurcation analysis](@ref ref), [identifiability analysis](@ref ref), or [compute steady states using homotopy continuation](@ref ref) using chemical reaction network models. Generally, most types of analysis supported by Catalyst is supported for coupled models. below follows some exceptions and notes: |
| 112 | +- Methods related to steady state analysis internally converts all systems to `NonlinearSystem`s. This means that all differential terms are set to `0`. |
| 113 | +- Homotopy continuation can only be applied if all equations also corresponds to (rational) polynomials. |
| 114 | + |
| 115 | + |
| 116 | +## [Special cases](@id coupled_equations_other) |
| 117 | + |
| 118 | +### [Defining custom differentials](@id coupled_equations_other_custom_differentials) |
| 119 | + |
| 120 | +### [Using higher-order differentials](@id coupled_equations_other_higher_order_differentials) |
| 121 | + |
| 122 | + |
| 123 | +## [Coupled equations for programmatic and compositional modelling](@id coupled_equations_non_dsl) |
| 124 | + |
| 125 | +### [Coupled equations for programmatic modelling](@id coupled_equations_non_dsl_programmatic) |
| 126 | + |
| 127 | +### [Coupled equations for compositional modelling](@id coupled_equations_non_dsl_compositional) |
0 commit comments