Skip to content

Commit 2276e1c

Browse files
authored
Merge pull request #801 from SciML/equations_and_events_remake
Equations and events reupload
2 parents 63ffa92 + 8601ed1 commit 2276e1c

File tree

19 files changed

+2019
-96
lines changed

19 files changed

+2019
-96
lines changed

.github/workflows/CI.yml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ jobs:
1616
group:
1717
- Core
1818
version:
19-
- '1'
19+
- '1.10.2'
2020
steps:
2121
- uses: actions/checkout@v4
2222
- uses: julia-actions/setup-julia@v1

Project.toml

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ version = "13.5.1"
66
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
77
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
88
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
9+
DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
910
DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
1011
Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
1112
JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
@@ -46,20 +47,21 @@ JumpProcesses = "9.3.2"
4647
LaTeXStrings = "1.3.0"
4748
Latexify = "0.14, 0.15, 0.16"
4849
MacroTools = "0.5.5"
49-
ModelingToolkit = "9.7.1"
50+
ModelingToolkit = "9.11.0"
5051
Parameters = "0.12"
5152
Reexport = "0.2, 1.0"
5253
Requires = "1.0"
5354
RuntimeGeneratedFunctions = "0.5.12"
5455
Setfield = "1"
5556
StructuralIdentifiability = "0.5.1"
5657
SymbolicUtils = "1.0.3"
57-
Symbolics = "5.14"
58+
Symbolics = "5.27"
5859
Unitful = "1.12.4"
5960
julia = "1.9"
6061

6162
[extras]
6263
BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
64+
DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
6365
DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
6466
Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
6567
HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
@@ -80,4 +82,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
8082
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
8183

8284
[targets]
83-
test = ["BifurcationKit", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]
85+
test = ["BifurcationKit", "DiffEqCallbacks", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]

docs/src/catalyst_functionality/constraint_equations.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ $\lambda$. For now we'll assume the cell can grow indefinitely. We'll also keep
1717
track of one protein $P(t)$, which is produced at a rate proportional $V$ and
1818
can be degraded.
1919

20-
## [Coupling ODE constraints via extending a system](@id constraint_equations_coupling_constratins)
20+
## [Coupling ODE constraints via extending a system](@id constraint_equations_coupling_constraints)
2121

2222
There are several ways we can create our Catalyst model with the two reactions
2323
and ODE for $V(t)$. One approach is to use compositional modeling, create

docs/src/catalyst_functionality/dsl_description.md

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -701,3 +701,66 @@ rn = @reaction_network begin
701701
end
702702
nothing # hide
703703
```
704+
705+
## Incorporating (differential) equations into reaction network models
706+
Some models cannot be purely described as reaction networks. E.g. consider the growth of a cell, where the rate of change in the cell's volume depends on some growth factor. Here, the cell's volume would be described by a normal ODE. 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$) is linear in the amount of growth factor:
707+
```@example eqs1
708+
using Catalyst #hide
709+
rn = @reaction_network begin
710+
@equations begin
711+
D(V) ~ G
712+
end
713+
(p,d), 0 <--> G
714+
end
715+
```
716+
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.
717+
718+
We can check the differential equation corresponding to this reaction network using latexify:
719+
```@example eqs1
720+
using Latexify
721+
latexify(rn; form=:ode)
722+
```
723+
We can also simulate it using the normal syntax
724+
```@example eqs1
725+
using DifferentialEquations, Plots # hide
726+
u0 = [:G => 0.0, :V => 0.1]
727+
ps = [:p => 1.0, :d => 0.5]
728+
oprob = ODEProblem(rn, u0, (0.0, 1.0), ps)
729+
sol = solve(oprob)
730+
plot(sol)
731+
```
732+
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.
733+
734+
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 parameter using the `@parameters` option:
735+
```@example eqs1
736+
rn = @reaction_network begin
737+
@parameters k
738+
@equations begin
739+
D(V) ~ k*G
740+
end
741+
(p,d), 0 <--> G
742+
end
743+
nothing #hide
744+
```
745+
If the differential does not appear isolated on the lhs, its subject variable must also be explicitly declared (as it is not inferred for these cases).
746+
747+
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:
748+
```@example eqs1
749+
rn = @reaction_network begin
750+
@equations begin
751+
D(V) ~ G
752+
D(N) ~ -G
753+
end
754+
(p,d), 0 <--> G
755+
end
756+
nothing #hide
757+
```
758+
759+
When only a single equation is added, the `begin ... end` statement can be omitted. E.g., the first model can be declared equivalently using:
760+
```@example eqs1
761+
rn = @reaction_network begin
762+
@equations D(V) ~ G
763+
(p,d), 0 <--> G
764+
end
765+
nothing # hide
766+
```

ext/CatalystHomotopyContinuationExtension.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@ module CatalystHomotopyContinuationExtension
22

33
# Fetch packages.
44
using Catalyst
5+
import DynamicPolynomials
56
import ModelingToolkit as MT
67
import HomotopyContinuation as HC
78
import Setfield: @set

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -53,7 +53,8 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
5353
eqs_pars_funcs = vcat(equations(ns), conservedequations(rs))
5454
eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs)
5555
eqs_intexp = make_int_exps.(eqs)
56-
return Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
56+
ss_poly = Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
57+
return poly_type_convert(ss_poly)
5758
end
5859

5960
# If u0s are not given while conservation laws are present, throws an error.
@@ -94,7 +95,7 @@ end
9495
function reorder_sols!(sols, ss_poly, rs::ReactionSystem)
9596
var_names_extended = String.(Symbol.(HC.variables(ss_poly)))
9697
var_names = [Symbol(s[1:prevind(s, findlast('_', s))]) for s in var_names_extended]
97-
sort_pattern = indexin(MT.getname.(species(rs)), var_names)
98+
sort_pattern = indexin(MT.getname.(unknowns(rs)), var_names)
9899
foreach(sol -> permute!(sol, sort_pattern), sols)
99100
end
100101

@@ -104,4 +105,13 @@ function filter_negative_f(sols; neg_thres=-1e-20)
104105
(neg_thres < sol[idx] < 0) && (sol[idx] = 0)
105106
end
106107
return filter(sol -> all(>=(0), sol), sols)
108+
end
109+
110+
# Sometimes (when polynomials are created from coupled CRN/DAEs), the steady state polynomial have the wrong type.
111+
# This converts it to the correct type, which homotopy continuation can handle.
112+
const WRONG_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}}}
113+
const CORRECT_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}, Float64}}
114+
function poly_type_convert(ss_poly)
115+
(typeof(ss_poly) == WRONG_POLY_TYPE) && return convert(CORRECT_POLY_TYPE, ss_poly)
116+
return ss_poly
107117
end

src/Catalyst.jl

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -30,7 +30,8 @@ using ModelingToolkit: Symbolic, value, istree, get_unknowns, get_ps, get_iv, ge
3030

3131
import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_variables!,
3232
modified_unknowns!, validate, namespace_variables,
33-
namespace_parameters, rename, renamespace, getname, flatten
33+
namespace_parameters, rename, renamespace, getname, flatten,
34+
is_alg_equation, is_diff_equation
3435

3536
# internal but needed ModelingToolkit functions
3637
import ModelingToolkit: check_variables,
@@ -42,6 +43,7 @@ import MacroTools, Graphs
4243
import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne
4344
import DataStructures: OrderedDict, OrderedSet
4445
import Parameters: @with_kw_noshow
46+
import Symbolics: occursin, wrap
4547

4648
# globals for the modulate
4749
function default_time_deriv()

src/expression_utils.jl

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,13 @@ function get_tup_arg(ex::ExprValues, i::Int)
1111
return ex.args[i]
1212
end
1313

14+
# Some options takes input on form that is either `@option ...` or `@option begin ... end`.
15+
# This transforms input of the latter form to the former (with only one line in the `begin ... end` block)
16+
function option_block_form(expr)
17+
(expr.head == :block) && return expr
18+
return Expr(:block, expr)
19+
end
20+
1421
# In variable/species/parameters on the forms like:
1522
# X
1623
# X = 1.0

src/reaction_network.jl

Lines changed: 117 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -82,7 +82,8 @@ const forbidden_variables_error = let
8282
end
8383

8484
# Declares the keys used for various options.
85-
const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling)
85+
const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling,
86+
:differentials, :equations, :continuous_events, :discrete_events)
8687

8788
### The @species macro, basically a copy of the @variables macro. ###
8889
macro species(ex...)
@@ -353,20 +354,28 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
353354
option_lines = Expr[x for x in ex.args if x.head == :macrocall]
354355

355356
# Get macro options.
356-
length(unique(arg.args[1] for arg in option_lines)) < length(option_lines) && error("Some options where given multiple times.")
357+
if length(unique(arg.args[1] for arg in option_lines)) < length(option_lines)
358+
error("Some options where given multiple times.")
359+
end
357360
options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
358361
option_lines))
359362

360363
# Reads options.
361364
default_reaction_metadata = :([])
362365
check_default_noise_scaling!(default_reaction_metadata, options)
363366
compound_expr, compound_species = read_compound_options(options)
367+
continuous_events_expr = read_events_option(options, :continuous_events)
368+
discrete_events_expr = read_events_option(options, :discrete_events)
364369

365370
# Parses reactions, species, and parameters.
366371
reactions = get_reactions(reaction_lines)
367372
species_declared = [extract_syms(options, :species); compound_species]
368373
parameters_declared = extract_syms(options, :parameters)
369-
variables = extract_syms(options, :variables)
374+
variables_declared = extract_syms(options, :variables)
375+
376+
# Reads more options.
377+
vars_extracted, add_default_diff, equations = read_equations_options(options, variables_declared)
378+
variables = vcat(variables_declared, vars_extracted)
370379

371380
# handle independent variables
372381
if haskey(options, :ivs)
@@ -391,6 +400,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
391400
species = vcat(species_declared, species_extracted)
392401
parameters = vcat(parameters_declared, parameters_extracted)
393402

403+
# Create differential expression.
404+
diffexpr = create_differential_expr(options, add_default_diff, [species; parameters; variables], tiv)
405+
394406
# Checks for input errors.
395407
(sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
396408
error("@reaction_network input contain $(length(ex.args) - sum(length.([reaction_lines,option_lines]))) malformed lines.")
@@ -402,7 +414,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
402414

403415
# Creates expressions corresponding to actual code from the internal DSL representation.
404416
sexprs = get_sexpr(species_extracted, options; iv_symbols = ivs)
405-
vexprs = haskey(options, :variables) ? options[:variables] : :()
417+
vexprs = get_sexpr(vars_extracted, options, :variables; iv_symbols = ivs)
406418
pexprs = get_pexpr(parameters_extracted, options)
407419
ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps")
408420
vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars")
@@ -412,6 +424,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
412424
for reaction in reactions
413425
push!(rxexprs.args, get_rxexprs(reaction))
414426
end
427+
for equation in equations
428+
push!(rxexprs.args, equation)
429+
end
415430

416431
quote
417432
$ps
@@ -420,11 +435,13 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
420435
$sps
421436
$observed_vars
422437
$comps
438+
$diffexpr
423439

424440
Catalyst.remake_ReactionSystem_internal(
425441
Catalyst.make_ReactionSystem_internal(
426442
$rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
427-
$pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs);
443+
$pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs,
444+
continuous_events = $continuous_events_expr, discrete_events = $discrete_events_expr);
428445
default_reaction_metadata = $default_reaction_metadata
429446
)
430447
end
@@ -799,6 +816,101 @@ function make_observed_eqs(observables_expr)
799816
end
800817
end
801818

819+
# Reads the variables options. Outputs:
820+
# `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
821+
# `dtexpr`: If a differentialequation is defined, the default derrivative (D ~ Differential(t)) must be defined.
822+
# `equations`: a vector with the equations provided.
823+
function read_equations_options(options, variables_declared)
824+
# Prepares the equations. First, extracts equations from provided option (converting to block form if requried).
825+
# Next, uses MTK's `parse_equations!` function to split input into a vector with the equations.
826+
eqs_input = haskey(options, :equations) ? options[:equations].args[3] : :(begin end)
827+
eqs_input = option_block_form(eqs_input)
828+
equations = Expr[]
829+
ModelingToolkit.parse_equations!(Expr(:block), equations, Dict{Symbol, Any}(), eqs_input)
830+
831+
# Loops through all equations, checks for lhs of the form `D(X) ~ ...`.
832+
# When this is the case, the variable X and differential D are extracted (for automatic declaration).
833+
# Also performs simple error checks.
834+
vars_extracted = Vector{Symbol}()
835+
add_default_diff = false
836+
for eq in equations
837+
if (eq.head != :call) || (eq.args[1] != :~)
838+
error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")
839+
end
840+
841+
# Checks if the equation have the format D(X) ~ ... (where X is a symbol). This means that the
842+
# default differential has been used. X is added as a declared variable to the system, and
843+
# we make a note that a differential D = Differential(iv) should be made as well.
844+
lhs = eq.args[2]
845+
# if lhs: is an expression. Is a function call. The function's name is D. Calls a single symbol.
846+
if (lhs isa Expr) && (lhs.head == :call) && (lhs.args[1] == :D) && (lhs.args[2] isa Symbol)
847+
diff_var = lhs.args[2]
848+
if in(diff_var, forbidden_symbols_error)
849+
error("A forbidden symbol ($(diff_var)) was used as an variable in this differential equation: $eq")
850+
end
851+
add_default_diff = true
852+
in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
853+
end
854+
end
855+
856+
return vars_extracted, add_default_diff, equations
857+
end
858+
859+
# Creates an expression declaring differentials. Here, `tiv` is the time independent variables,
860+
# which is used by the default differential (if it is used).
861+
function create_differential_expr(options, add_default_diff, used_syms, tiv)
862+
# Creates the differential expression.
863+
# If differentials was provided as options, this is used as the initial expression.
864+
# If the default differential (D(...)) was used in equations, this is added to the expression.
865+
diffexpr = (haskey(options, :differentials) ? options[:differentials].args[3] : MacroTools.striplines(:(begin end)))
866+
diffexpr = option_block_form(diffexpr)
867+
868+
# Goes through all differentials, checking that they are correctly formatted and their symbol is not used elsewhere.
869+
for dexpr in diffexpr.args
870+
(dexpr.head != :(=)) && error("Differential declaration must have form like D = Differential(t), instead \"$(dexpr)\" was given.")
871+
(dexpr.args[1] isa Symbol) || error("Differential left-hand side must be a single symbol, instead \"$(dexpr.args[1])\" was given.")
872+
in(dexpr.args[1], used_syms) && error("Differential name ($(dexpr.args[1])) is also a species, variable, or parameter. This is ambigious and not allowed.")
873+
in(dexpr.args[1], forbidden_symbols_error) && error("A forbidden symbol ($(dexpr.args[1])) was used as a differential name.")
874+
end
875+
876+
# If the default differential D has been used, but not pre-declared using the @differenitals
877+
# options, add this declaration to the list of declared differentials.
878+
if add_default_diff && !any(diff_dec.args[1] == :D for diff_dec in diffexpr.args)
879+
push!(diffexpr.args, :(D = Differential($(tiv))))
880+
end
881+
882+
return diffexpr
883+
end
884+
885+
# Read the events (continious or discrete) provided as options to the DSL. Returns an expression which evalutes to these.
886+
function read_events_option(options, event_type::Symbol)
887+
# Prepares the events, if required to, converts them to block form.
888+
(event_type in [:continuous_events, :discrete_events]) || error("Trying to read an unsupported event type.")
889+
events_input = haskey(options, event_type) ? options[event_type].args[3] : MacroTools.striplines(:(begin end))
890+
events_input = option_block_form(events_input)
891+
892+
# Goes throgh the events, checks for errors, and adds them to the output vector.
893+
events_expr = :([])
894+
for arg in events_input.args
895+
# Formatting error checks.
896+
# NOTE: Maybe we should move these deeper into the system (rather than the DSL), throwing errors more generally?
897+
if (arg isa Expr) && (arg.head != :call) || (arg.args[1] != :(=>)) || length(arg.args) != 3
898+
error("Events should be on form `condition => affect`, separated by a `=>`. This appears not to be the case for: $(arg).")
899+
end
900+
if (arg isa Expr) && (arg.args[2] isa Expr) && (arg.args[2].head != :vect) && (event_type == :continuous_events)
901+
error("The condition part of continious events (the left-hand side) must be a vector. This is not the case for: $(arg).")
902+
end
903+
if (arg isa Expr) && (arg.args[3] isa Expr) && (arg.args[3].head != :vect)
904+
error("The affect part of all events (the righ-hand side) must be a vector. This is not the case for: $(arg).")
905+
end
906+
907+
# Adds the correctly formatted event to the event creation expression.
908+
push!(events_expr.args, arg)
909+
end
910+
911+
return events_expr
912+
end
913+
802914
# Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a
803915
# default metadata value to the `default_reaction_metadata` vector.
804916
function check_default_noise_scaling!(default_reaction_metadata, options)

0 commit comments

Comments
 (0)