Enable equations in the DSL

HISTORY.md

@@ -24,6 +24,13 @@
   ```
   X's value will be `1.0` in the first vertex, but `0.0` in the remaining one (the system have 25 vertexes in total). SInce th parameters `p` and `d` are part of the non-spatial reaction network, their values are tied to vertexes. However, if the `D` parameter (which governs diffusion between vertexes) is given several values, these will instead correspond to the specific edges (and transportation along those edges.)
 
+- Update how compounds are created. E.g. use
+```julia
+@variables t C(t) O(t)
+@compound CO2 ~ C + 2O
+```
+to create a compound species `CO2` that consists of `C` and 2 `O`.
+- Added documentation for chemistry related functionality (compound creation and reaction balancing).
 - Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage:
 ```julia
 using Catalyst

Project.toml

@@ -61,6 +61,7 @@ HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -73,4 +74,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["BifurcationKit", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]
+test = ["BifurcationKit", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]

docs/pages.jl

@@ -7,6 +7,7 @@ pages = Any["Home" => "index.md",
                                             "catalyst_functionality/constraint_equations.md",
                                             "catalyst_functionality/parametric_stoichiometry.md",
                                             "catalyst_functionality/network_analysis.md",
+                                            "catalyst_functionality/chemistry_related_functionality.md",
                                             "Model creation examples" => Any["catalyst_functionality/example_networks/basic_CRN_examples.md",
                                                                              "catalyst_functionality/example_networks/hodgkin_huxley_equation.md",
                                                                              "catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md"]],

docs/src/api/catalyst_api.md

@@ -157,6 +157,10 @@ species
 nonspecies
 reactionparams
 reactions
+has_diff_equations
+diff_equations
+has_alg_equations
+alg_equations
 numspecies
 numparams
 numreactions
@@ -280,6 +284,17 @@ Base.convert
 ModelingToolkit.structural_simplify
 ```
 
+## Chemistry-related functionalities
+Various functionalities primarily relevant to modelling of chemical systems (but potentially also in biology). 
+```@docs
+@compound
+@compounds
+iscompound
+components
+coefficients
+component_coefficients
+```
+
 ## Unit validation
 ```@docs
 validate(rx::Reaction; info::String = "")

docs/src/catalyst_functionality/chemistry_related_functionality.md

@@ -0,0 +1,129 @@
+# [Chemistry-related functionality](@id chemistry_functionality)
+
+While Catalyst has primarily been designed around the modelling of biological systems, reaction network models are also common across chemistry. This section describes two types of functionality, that while of general interest, should be especially useful in the modelling of chemical systems.
+- The `@compound` option, which enables the user to designate that a specific species is composed of certain subspecies.
+- The `balance_reaction` function, which enables the user to balance a reaction so the same number of components occur on both sides.
+
+## Modelling with compound species
+
+#### Creating compound species programmatically
+We will first show how to create compound species through [programmatic model construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and their stoichiometries). In this example, we will create a CO₂ molecule, consisting of one C atom and two O atoms. First, we create species corresponding to the components:
+```@example chem1
+using Catalyst
+@variables t
+@species C(t) O(t) 
+```
+Next, we create the `CO2` compound species:
+```@example chem1
+@compound CO2 ~ C + 2O
+```
+Here, the compound is the first argument to the macro, followed by its component (with the left-hand and right-hand sides separated by a `~` sign). While non-compound species (such as `C` and `O`) have their independent variable (in this case `t`) designated, independent variables are not designated for compounds (these are instead directly inferred from their components). Components with non-unitary stoichiometries have this value written before the component (generally, the rules for designating the components of a compound are identical to those of designating the substrates or products of a reaction). The created compound, `CO2`, is a species in every sense, and can be used wherever e.g. `C` can be used:
+```@example chem1
+isspecies(CO2)
+```
+In its metadata, however, is stored information of its components, which can be retrieved using the `components` (returning a vector of its component species) and `coefficients` (returning a vector with each component's stoichiometry) functions:
+```@example chem1
+components(CO2)
+```
+```@example chem1
+coefficients(CO2)
+```
+Alternatively, we can retrieve the components and their stoichiometric coefficients as a single vector using the `component_coefficients` function:
+```@example chem1
+component_coefficients(CO2)
+```
+Finally, it is possible to check whether a species is a compound or not using the `iscompound` function:
+```@example chem1
+iscompound(CO2)
+```
+
+Compound components that are also compounds are allowed, e.g. we can create a carbonic acid compound (H₂CO₃) that consists of CO₂ and H₂O:
+```@example chem1
+@species H(t)
+@compound H2O ~ 2H + O
+@compound H2CO3 ~ CO2 + H2O
+```
+
+When multiple compounds are created, they can be created simultaneously using the `@compounds` macro, e.g. the previous code-block can be re-written as:
+```@example chem1
+@species H(t)
+@compounds begin
+    H2O ~ 2H + O
+    H2CO3 ~ CO2 + H2O
+end
+```
+
+#### Creating compound species within the DSL
+It is also possible to declare species as compound species within the `@reaction_network` DSL, using the `@compounds` options:
+```@example chem1
+rn = @reaction_network begin
+    @species C(t) H(t) O(t)
+    @compounds begin
+        C2O ~ C + 2O
+        H2O ~ 2H + O
+        H2CO3 ~ CO2 + H2O
+    end
+    (k1,k2), H2O+ CO2 <--> H2CO3
+end
+```
+When creating compound species using the DSL, it is important to note that *every component must be known to the system as a species, either by being declared using the `@species` option, or by appearing in a reaction*. E.g. the following is not valid
+```julia 
+rn = @reaction_network begin
+    @compounds begin
+        C2O ~ C + 2O
+        H2O ~ 2H + O
+        H2CO3 ~ CO2 + H2O
+    end
+    (k1,k2), H2O+ CO2 <--> H2CO3
+end
+```
+as the components `C`, `H`, and `O` are not declared as a species anywhere. Please also note that only `@compounds` can be used as an option in the DSL, not `@compound`.
+
+#### Designating metadata and default values for compounds
+Just like for normal species, it is possible to designate metadata and default values for compounds. Metadata is provided after the compound name, but separated from it by a `,`:
+```@example chem1
+@compound (CO2, [unit="mol"]) ~ C + 2O
+```
+Default values are designated using `=`, and provided directly after the compound name.:
+```@example chem1
+@compound (CO2 = 2.0) ~ C + 2O
+```
+If both default values and meta data are provided, the metadata is provided after the default value:
+```@example chem1
+@compound (CO2 = 2.0, [unit="mol"]) ~ C + 2O
+```
+In all of these cases, the side to the left of the `~` must be enclosed within `()`.
+
+## Balancing chemical reactions
+One use of defining a species as a compound is that they can be used to balance reactions to that the number of components are the same on both sides. Catalyst provides the `balance_reaction` function, which takes a reaction, and returns a balanced version. E.g. let us consider a reaction when carbon dioxide is formed from carbon and oxide `C + O --> CO2`. Here, `balance_reaction` enables us to find coefficients creating a balanced reaction (in this case, where the number of carbon and oxygen atoms are the same on both sides). To demonstrate, we first created the unbalanced reactions:
+```@example chem1
+rx = @reaction k, C + O --> $CO2
+```
+Here, the reaction rate (`k`) is not involved in the reaction balancing. We use interpolation for `CO2`, ensuring that the `CO2` used in the reaction is the same one we previously defined as a compound of `C` and `O`. Next, we call the `balance_reaction` function
+```@example chem1
+balance_reaction(rx)
+```
+which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that the reaction balancing problem may have several solutions. Typically, there is only a single solution (in which case this is the vector's only element). No, or an infinite number of, solutions is also possible.
+
+Let us consider a more elaborate example, the reaction between ammonia (NH₃) and oxygen (O₂) to form nitrogen monoxide (NO) and water (H₂O). Let us first create the components and the unbalanced reaction:
+```@example chem2
+using Catalyst # hide
+@variables t
+@species N(t) H(t) O(t) 
+@compounds begin
+    NH3 ~ N + 3H
+    O2 ~ 2O
+    NO ~ N + O
+    H2O ~ 2H + O
+end
+unbalanced_reaction = @reaction k, $NH3 + $O2 --> $NO + $H2O
+```
+We can now create a balanced version (where the amount of H, N, and O is the same on both sides):
+```@example chem2
+balanced_reaction = balance_reaction(unbalanced_reaction)[1]
+```
+
+Reactions declared as a part of a `ReactionSystem` (e.g. using the DSL) can be retrieved for balancing using the `reaction` function. Please note that balancing these will not mutate the `ReactionSystem`, but a new reaction system will need to be created using the balanced reactions.
+
+!!! note
+    Reaction balancing is currently not supported for reactions involving compounds of compounds.

docs/src/catalyst_functionality/constraint_equations.md

@@ -17,7 +17,7 @@ $\lambda$. For now we'll assume the cell can grow indefinitely. We'll also keep
 track of one protein $P(t)$, which is produced at a rate proportional $V$ and
 can be degraded.
 
-## Coupling ODE constraints via extending a system
+## [Coupling ODE constraints via extending a system](@id constraint_equations_coupling_constratins)
 
 There are several ways we can create our Catalyst model with the two reactions
 and ODE for $V(t)$. One approach is to use compositional modeling, create

docs/src/catalyst_functionality/dsl_description.md

@@ -567,3 +567,123 @@ species(rn)
 !!! note
     When using interpolation, expressions like `2$spec` won't work; the
     multiplication symbol must be explicitly included like `2*$spec`.
+
+## Including observables
+Sometimes, one might want to include observable variables. These are variables that can be computed directly from the other system variables (rather than having their values implicitly given through some differential equation). These can be introduced through the `@observables` option.
+
+Let us consider a simple example where two species ($X$ and $Y$) are produced and degraded at constant rates. They can also bind, forming a complex ($XY$). If we want to access the total amount of $X$ in the system we can create an observable that denotes this quantity ($Xtot = X + XY$). Here, we create observables for the total amount of $X$ and $Y$:
+```@example obs1
+using Catalyst # hide
+rn = @reaction_network begin
+  @observables begin
+    Xtot ~ X + XY
+    Ytot ~ Y + XY
+  end
+  (pX,dX), 0 <--> X
+  (pY,dY), 0 <--> Y
+  (kB,kD), X + Y <--> XY
+end
+```
+The `@observables` option is followed by one line for each observable formula (enclosed by a `begin ... end` block). The left-hand sides indicate the observables' names, and the right-hand sides how their values are computed. The two sides are separated by a `~`. 
+
+If we now simulate our model:
+```@example obs1
+using DifferentialEquations # hide
+u0 = [:X => 0.0, :Y => 0.0, :XY => 0.0]
+tspan = (0.0, 10.0)
+ps = [:pX => 1.0, :dX => 0.2, :pY => 1.0, :dY => 0.5, :kB => 1.0, :kD => 0.2]
+oprob = ODEProblem(rn, u0, tspan, ps)
+sol = solve(oprob)
+nothing # hide
+```
+we can index the solution using our observables (just like for [other variables](@ref simulation_structure_interfacing_solutions)). E.g. we can receive a vector with all $Xtot$ values using
+```@example obs1
+sol[:Xtot]
+```
+similarly, we can plot the values of $Xtot$ and $Ytot$ using
+```@example obs1
+plot(sol; idxs=[:Xtot, :Ytot], label=["Total X" "Total Y"])
+```
+
+If we only wish to provide a single observable, the `begin ... end` block is note required. E.g., to record only the total amount of $X$ we can use:
+```@example obs1
+using Catalyst # hide
+rn = @reaction_network begin
+  @observables Xtot ~ X + XY
+  (pX,dX), 0 <--> X
+  (pY,dY), 0 <--> Y
+  (kB,kD), X + Y <--> XY
+end
+```
+
+Finally, some general rules for creating observables:
+- Observables can depend on any species, parameters, or variables, but not on other observables.
+- All observables components must be declared somewhere (i.e., they cannot only appear as a part of the observables formula).
+- Only a single `@observables` option block can be used in each `@reaction_network` call.
+- The left-hand side of the observables expression must be a single symbol, indicating the observable's name.
+- Metadata can, however, be provided, e.g through `@observables (Xtot, [description="Total amount of X"]) ~ X + XY`.
+- The right-hand side of the observables expression can be any valid algebraic expression.
+- Observables are (by default, but this can be changed) considered `variables` (and not `species`). This can be changed by e.g. pre-declaring them using the `@species` option.
+
+
+## Incorporating (differential) equations into reaction network models
+Some models cannot be purely described as a reaction network. E.g. consider the growth of a cell, where the rate of change in cell's volume depends on some growth factor. Here, the cell's volume would be described by a normal equation. Such equations can be incorporated into a model using the `@equations` option. Here, we create a model where a growth factor ($G$) is produced and degraded at a linear rates, and the rate of change in cell volume ($V$) os linear to the amount of growth factor:
+```@example eqs1
+using Catalyst #hide
+rn = @reaction_network begin
+  @equations begin
+    D(V) ~ G
+  end
+  (p,d), 0 <--> G
+end
+```
+Here, `D(V)` indicates the (time) derivative with respect to `D`. The differential equation left and right hand sides are separated by a `~`. The left-hand side should contain differential only, the right hand side can contain any algebraic expression.
+
+We can check the differential equation corresponding to this reaction network using latexify:
+```@example eqs1
+using Latexify
+latexify(rn; form=:ode)
+```
+We can also simulate it using the normal syntax
+```@example eqs1
+using DifferentialEquations, Plots # hide
+u0 = [:G => 0.0, :V => 0.1]
+ps = [:p => 1.0, :d => 0.5]
+oprob = ODEProblem(rn, u0, (0.0, 1.0), ps)
+sol = solve(oprob)
+plot(sol)
+```
+Here, growth is indefinite. To improve the model, [a callback](@ref advanced_simulations_callbacks) can be used to half the volume (cell division) once some threshold is reached.
+
+When creating differential equations this way, the subject of the differential is automatically inferred to be a variable, however, any component on the right-hand side must be declare somewhere in the macro. E.g. to add a scaling parameter ($k$), we must declare that $k$ is a parmaeter using the `@paraemters` option:
+```@example eqs1
+rn = @reaction_network begin
+  @parameters k
+  @equations begin
+    D(V) ~ k*G
+  end
+  (p,d), 0 <--> G
+end
+nothing #hide
+```
+
+It is possible to add several equations to the model. In this case, each have a separate line. E.g. to keep track of a supply of nutrition ($N$) in the growth media, we can use:
+```@example eqs1
+rn = @reaction_network begin
+  @equations begin
+    D(V) ~ G
+    D(N) ~ -G
+  end
+  (p,d), 0 <--> G
+end
+nothing #hide
+```
+
+When only a single equation is added, the `begin ... end` statement can be omitted. E.g., the first model can be declared equivalently using:
+```@example eqs1
+rn = @reaction_network begin
+  @equations D(V) ~ G
+  (p,d), 0 <--> G
+end
+nothing # hide
+```

src/Catalyst.jl

@@ -79,6 +79,7 @@ export mm, mmr, hill, hillr, hillar
 # functions to query network properties
 include("networkapi.jl")
 export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap
+export has_diff_equations, diff_equations, has_alg_equations, alg_equations
 export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap
 export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
@@ -102,8 +103,8 @@ export Graph, savegraph, complexgraph
 
 # for creating compounds
 include("chemistry_functionality.jl")
-export @compound
-export components, iscompound, coefficients
+export @compound, @compounds
+export iscompound, components, coefficients, component_coefficients
 export balance_reaction
 
 ### Extensions ###