You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
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.
4
+
- The `@compound` option, allowing the user to designate that a specific species is composed of certain subspecies.
5
+
- The `balance_reaction` function enabling the user to balance a reactions so the same number of compounds occur on both sides.
6
+
7
+
## Modelling with compound species
8
+
Defining compound species is currently only supported for [programmatic construction](@ref programmatic_CRN_construction) of reactions and reaction network models. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ compound species, consisting of 1 C species and 2 O species. First we create species corresponding to the components:
9
+
```@example chem1
10
+
@variables t
11
+
@species C(t) O(t)
12
+
```
13
+
Next, we create the `CO2` compound:
14
+
```@example chem1
15
+
@compound CO2(t) = C + 2O
16
+
```
17
+
Here, the compound is the first argument to the macro, followed by its component. The `(t)` indicates that `CO2` is a time-dependant variable. Components with non-unitary stoichiometry have this value written before the component (generally, the rules for designating the components of a compounds are identical to hose 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` could be used:
18
+
```@example chem1
19
+
isspecies(CO2)
20
+
```
21
+
However, in its metadata is stored the 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:
22
+
```@example chem1
23
+
components(CO2)
24
+
```
25
+
```@example chem1
26
+
coefficients(CO2)
27
+
```
28
+
Alternatively, we can retrieve the components and their stoichiometric coefficients as a single vector using the `component_coefficients` function:
29
+
```@example chem1
30
+
component_coefficients(CO2)
31
+
```
32
+
Finally, it is possible to check whether a species is a compound or not using the `iscompound` function:
33
+
```@example chem1
34
+
iscompound(CO2)
35
+
```
36
+
37
+
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:
38
+
```@example chem1
39
+
@species H(t)
40
+
@compound H2O(t) = 2H + O
41
+
@compound H2CO3(t) = CO2 + H2O
42
+
```
43
+
44
+
When multiple compounds are created, they can be created simultaneously using the `@compounds` macro, e.g. the previous code-block could have been executed using:
45
+
```@example chem1
46
+
@species H(t)
47
+
@compounds begin
48
+
H2O(t) = 2H + O
49
+
H2CO3(t) = CO2 + H2O
50
+
end
51
+
```
52
+
53
+
One use defining a species as a compound is that they can be used to balance reactions to that the number of compounds are the same on both side.
54
+
55
+
## Balancing chemical reactions
56
+
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` enable us to find coefficients creating a balanced reaction (in this case, where the number of carbon and oxygen atoms are teh same on both sides). To demonstrate, we first created the unbalanced reactions:
57
+
```@example chem1
58
+
rx = @reaction k, C + O --> $CO2
59
+
```
60
+
Here, we set a reaction rate `k` (which is not involved in the reaction balancing). We also use interpolation of `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
61
+
```@example chem1
62
+
balance_reaction(rx)
63
+
```
64
+
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).
65
+
66
+
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:
Macro that creates a compound species, which is composed of smaller constituent species.
16
+
17
+
Example:
18
+
```julia
19
+
@variables t
20
+
@species C(t) O(t)
21
+
@compound CO2(t) = C + 2O
22
+
```
23
+
24
+
Notes:
25
+
- The constituent species must be defined before using the `@compound` macro.
26
+
"""
27
+
macrocompound(expr)
28
+
make_compound(MacroTools.striplines(expr))
29
+
end
46
30
47
-
coeffs_expr =Expr(:vect, coeffs...)
48
-
species_expr =Expr(:vect, species...)
31
+
# Function managing the @compound macro.
32
+
functionmake_compound(expr)
33
+
# Error checks.
34
+
!(expr isa Expr) || (expr.head != :(=)) &&error("Malformed expression. Compounds should be declared using a \"=\".")
35
+
(length(expr.args) !=2) &&error("Malformed expression. Compounds should be consists of two expression, separated by a \"=\".")
36
+
((expr.args[1] isa Symbol) || (expr.args[1].head !=:call)) &&error("Malformed expression. There must be a single compound which depend on an independent variable, e.g. \"CO2(t)\".")
37
+
38
+
# Extracts the composite species name, and a Vector{ReactantStruct} of its components.
components = :([]) # Extra step required here to get something like :([C, O]), rather than :([:C, :O])
43
+
foreach(comp ->push!(components.args, comp.reactant), composition) # E.g. [C, O]
44
+
coefficients =getfield.(composition, :stoichiometry) # E.g. [1, 2]
45
+
46
+
# Creates the found expressions that will create the compound species.
47
+
# The `Expr(:escape, :(...))` is required so that teh expressions are evaluated in the scope the users use the macro in (to e.g. detect already exiting species).
48
+
species_declaration_expr =Expr(:escape, :(@species$species_expr)) # E.g. `@species CO2(t)`
Macro that creates several compound species, which each is composed of smaller constituent species. Uses the same syntax as `@compound`, but with one compound species one each line.
57
66
58
-
# Construct the expression to set the coefficients metadata
# Creates a separate compound call for each compound.
89
+
compound_calls = [make_compound(line) for line in expr.args] # Crates a vector containing the quote for each compound.
90
+
returnExpr(:block, vcat([c_call.args for c_call in compound_calls]...)...) # Combines the quotes to a single one. Don't want the double "...". But had problems getting past them for various metaprogramming reasons :(.)
68
91
end
69
92
70
-
# Check if a species is a compound
93
+
## Compound Getters ###
94
+
95
+
"""
96
+
iscompound(s)
97
+
98
+
Returns `true` if the input is a compound species (else false).
99
+
"""
71
100
iscompound(s::Num) =iscompound(MT.value(s))
72
101
functioniscompound(s)
73
102
MT.getmetadata(s, CompoundSpecies, false)
74
103
end
75
104
76
-
coefficients(s::Num) =coefficients(MT.value(s))
77
-
functioncoefficients(s)
78
-
MT.getmetadata(s, CompoundCoefficients)
79
-
end
105
+
"""
106
+
components(s)
80
107
108
+
Returns a vector with a list of all the components of a compound species (created using e.g. the @compound macro).
109
+
"""
81
110
components(s::Num) =components(MT.value(s))
82
111
functioncomponents(s)
83
112
MT.getmetadata(s, CompoundComponents)
84
113
end
85
114
115
+
"""
116
+
coefficients(s)
117
+
118
+
Returns a vector with a list of all the stoichiometric coefficients of the components of a compound species (created using e.g. the @compound macro).
119
+
"""
120
+
coefficients(s::Num) =coefficients(MT.value(s))
121
+
functioncoefficients(s)
122
+
MT.getmetadata(s, CompoundCoefficients)
123
+
end
124
+
125
+
"""
126
+
component_coefficients(s)
127
+
128
+
Returns a Vector{Pari{Symbol,Int64}}, listing a compounds species (created using e.g. the @compound macro) all the coefficients and their stoichiometric coefficients.
@@ -143,6 +189,7 @@ function create_matrix(reaction::Catalyst.Reaction)
143
189
return A
144
190
end
145
191
192
+
# Internal function used by "balance_reaction".
146
193
functionget_balanced_stoich(reaction::Reaction)
147
194
# Create the reaction matrix A that is m atoms by n compounds
148
195
A =create_matrix(reaction)
@@ -171,6 +218,52 @@ function get_balanced_stoich(reaction::Reaction)
171
218
return stoichvecs
172
219
end
173
220
221
+
"""
222
+
balance_reaction(reaction::Reaction)
223
+
224
+
Returns a vector of all possible stoichiometrically balanced `Reaction` objects for the given `Reaction`.
225
+
226
+
Example:
227
+
```julia
228
+
@variables t
229
+
@species Si(t) Cl(t) H(t) O(t)
230
+
@compound SiCl4(t) = Si + 4Cl
231
+
@compound H2O(t) = 2H + O
232
+
@compound H4SiO4(t) = 4H + Si + 4O
233
+
@compound HCl(t) = H + Cl
234
+
rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl])
235
+
balance_reaction(rx) # Exactly one solution.
236
+
```
237
+
238
+
```julia
239
+
@variables t
240
+
@species C(t) H(t) O(t)
241
+
@compound CO(t) = C + O
242
+
@compound CO2(t) = C + 2O
243
+
@compound H2(t) = 2H
244
+
@compound CH4(t) = C + 4H
245
+
@compound H2O(t) = 2H + O
246
+
rx = Reaction(1.0, [CO, CO2, H2], [CH4, H2O])
247
+
balance_reaction(rx) # Multiple solutions.
248
+
```
249
+
250
+
```julia
251
+
@variables t
252
+
@species Fe(t) S(t) O(t) H(t) N(t)
253
+
@compound FeS2(t) = Fe + 2S
254
+
@compound HNO3(t) = H + N + 3O
255
+
@compound Fe2S3O12(t) = 2Fe + 3S + 12O
256
+
@compound NO(t) = N + O
257
+
@compound H2SO4(t) = 2H + S + 4O
258
+
rx = Reaction(1.0, [FeS2, HNO3], [Fe2S3O12, NO, H2SO4])
259
+
brxs = balance_reaction(rx) # No solution.
260
+
```
261
+
262
+
Notes:
263
+
- Balancing reactions that contain compounds of compounds is currently not supported.
264
+
- A reaction may not always yield a single solution; it could have an infinite number of solutions or none at all. When there are multiple solutions, a vector of all possible `Reaction` objects is returned. However, substrates and products may be interchanged as we currently do not solve for a linear combination that maintains the set of substrates and products.
265
+
- If the reaction cannot be balanced, an empty `Reaction` vector is returned.
266
+
"""
174
267
functionbalance_reaction(reaction::Reaction)
175
268
# Calculate the stoichiometric coefficients for the balanced reaction.
0 commit comments