Skip to content

Commit a299c84

Browse files
committed
init
1 parent afced2c commit a299c84

File tree

9 files changed

+474
-10
lines changed

9 files changed

+474
-10
lines changed

docs/src/api/catalyst_api.md

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -157,6 +157,10 @@ species
157157
nonspecies
158158
reactionparams
159159
reactions
160+
has_diff_equations
161+
diff_equations
162+
has_alg_equations
163+
alg_equations
160164
numspecies
161165
numparams
162166
numreactions

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: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -701,3 +701,65 @@ 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 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:
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 parmaeter using the `@paraemters` 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+
746+
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:
747+
```@example eqs1
748+
rn = @reaction_network begin
749+
@equations begin
750+
D(V) ~ G
751+
D(N) ~ -G
752+
end
753+
(p,d), 0 <--> G
754+
end
755+
nothing #hide
756+
```
757+
758+
When only a single equation is added, the `begin ... end` statement can be omitted. E.g., the first model can be declared equivalently using:
759+
```@example eqs1
760+
rn = @reaction_network begin
761+
@equations D(V) ~ G
762+
(p,d), 0 <--> G
763+
end
764+
nothing # hide
765+
```

src/Catalyst.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,7 @@ export mm, mmr, hill, hillr, hillar
8888
# functions to query network properties
8989
include("networkapi.jl")
9090
export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap
91+
export has_diff_equations, diff_equations, has_alg_equations, alg_equations
9192
export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap
9293
export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap
9394
export dependants, dependents, substoichmat, prodstoichmat, netstoichmat

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/networkapi.jl

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -97,6 +97,55 @@ function reactions(network)
9797
[rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
9898
end
9999

100+
"""
101+
diff_equations(network)
102+
Given a [`ReactionSystem`](@ref), return a vector of all `Equations` in the system that are differential equations (contains a derivative with respect to any variable).
103+
Notes:
104+
- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
105+
"""
106+
function diff_equations(network)
107+
eqs = equations(network)
108+
filter!(!isreaction, eqs)
109+
systems = filter_nonrxsys(network)
110+
isempty(systems) && (return rxs)
111+
[rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
112+
end
113+
114+
"""
115+
has_diff_equations(network)
116+
Given a [`ReactionSystem`](@ref), check whether it contain any differential equations (i.e. in addition to those generated through reactions).
117+
Notes:
118+
- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
119+
"""
120+
function has_diff_equations(network)
121+
return !isempty(diff_equations(network))
122+
end
123+
124+
"""
125+
alg_equations(network)
126+
Given a [`ReactionSystem`](@ref), return a vector of all `Equations` in the system that are algebraic equations (does not contain any derivatives).
127+
Notes:
128+
- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
129+
"""
130+
function alg_equations(network)
131+
eqs = equations(network)
132+
filter!(!isreaction, eqs)
133+
filter!(!isreaction, eqs)
134+
systems = filter_nonrxsys(network)
135+
isempty(systems) && (return rxs)
136+
[rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
137+
end
138+
139+
"""
140+
has_alg_equations(network)
141+
Given a [`ReactionSystem`](@ref), check whether it contain any algebraic equations.
142+
Notes:
143+
- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
144+
"""
145+
function has_alg_equations(network)
146+
return !isempty(alg_equations(network))
147+
end
148+
100149
"""
101150
speciesmap(network)
102151

src/reaction_network.jl

Lines changed: 69 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])
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)
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
@@ -810,6 +827,53 @@ function check_default_noise_scaling!(default_reaction_metadata, options)
810827
end
811828
end
812829

830+
# Creates an expression declaring differentials.
831+
function create_differential_expr(options, add_default_diff, used_syms)
832+
# Creates the differential expression.
833+
# If differentials was provided as options, this is used as the initial expression.
834+
# If the default differential (D(...)) was used in equations, this is added to the expression.
835+
diffexpr = (haskey(options, :differentials) ? options[:differentials].args[3] : MacroTools.striplines(:(begin end)))
836+
diffexpr = option_block_form(diffexpr)
837+
add_default_diff && push!(diffexpr.args, :(D = Differential($(DEFAULT_IV_SYM))))
838+
839+
# Goes through all differentials, checking that they are correctly formatted and their symbol is not used elsewhere.
840+
for dexpr in diffexpr.args
841+
(dexpr.head != :(=)) && error("Differential declaration must have form like D = Differential(t), instead \"$(dexpr)\" was given.")
842+
(dexpr.args[1] isa Symbol) || error("Differential left-hand side must be a single symbol, instead \"$(dexpr.args[1])\" was given.")
843+
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.")
844+
in(dexpr.args[1], forbidden_symbols_error) && error("A forbidden symbol ($(dexpr.args[1])) was used as a differential name.")
845+
end
846+
847+
return diffexpr
848+
end
849+
850+
# Read the events (continious or discrete) provided as options to the DSL. Returns an expression which evalutes to these.
851+
function read_events_option(options, event_type::Symbol)
852+
# Prepares the events, if required to, converts them to block form.
853+
(event_type in [:continuous_events]) || error("Trying to read an unsupported event type.")
854+
events_input = haskey(options, event_type) ? options[event_type].args[3] : :(begin end)
855+
events_input = option_block_form(events_input)
856+
857+
# Goes throgh the events, checks for errors, and adds them to the output vector.
858+
events_expr = :([])
859+
for arg in events_input.args
860+
# Formatting error checks.
861+
# NOTE: Maybe we should move these deeper into the system (rather than the DSL), throwing errors more generally?
862+
if (arg.head != :call) || (arg.args[1] != :(=>)) || length(arg.args) != 3
863+
error("Events should be on form `condition => affect`, separated by a `=>`. This appears not to be the case for: $(arg).")
864+
end
865+
if (arg.args[2] != :call) && (event_type == :continuous_events)
866+
error("The condition part of continious events (the left-hand side) must be a vector. This is not the case for: $(arg).")
867+
end
868+
if arg.args[3] != :call
869+
error("The affect part of all events (the righ-hand side) must be a vector. This is not the case for: $(arg).")
870+
end
871+
872+
push!(events_expr.args, arg)
873+
end
874+
return events_expr
875+
end
876+
813877
### Functionality for expanding function call to custom and specific functions ###
814878

815879
#Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.

src/reactionsystem.jl

Lines changed: 16 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -551,6 +551,11 @@ struct ReactionSystem{V <: NetworkProperties} <:
551551
(p isa Symbolics.BasicSymbolic) || error("Parameter $p is not a `BasicSymbolic`. This is required.")
552552
end
553553

554+
# Filters away any potential obervables from `states` and `spcs`.
555+
obs_vars = [obs_eq.lhs for obs_eq in observed]
556+
states = filter(state -> !any(isequal(state, obs_var) for obs_var in obs_vars), states)
557+
spcs = filter(spc -> !any(isequal(spc, obs_var) for obs_var in obs_vars), spcs)
558+
554559
# unit checks are for ODEs and Reactions only currently
555560
nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
556561
if checks && isempty(sivs)
@@ -1415,7 +1420,7 @@ end
14151420

14161421
# merge constraint components with the ReactionSystem components
14171422
# also handles removing BC and constant species
1418-
function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false)
1423+
function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false, zero_derivatives = false)
14191424
# if there are BC species, put them after the independent species
14201425
rssts = get_unknowns(rs)
14211426
sts = any(isbc, rssts) ? vcat(ists, filter(isbc, rssts)) : ists
@@ -1507,7 +1512,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
15071512
ists, ispcs = get_indep_sts(fullrs, remove_conserved)
15081513
eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
15091514
include_zero_odes)
1510-
eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
1515+
eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved, zero_derivatives=true)
15111516

15121517
ODESystem(eqs, get_iv(fullrs), sts, ps;
15131518
observed = obs,
@@ -1548,7 +1553,6 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
15481553
ists, ispcs = get_indep_sts(fullrs, remove_conserved)
15491554
eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
15501555
as_odes = false, include_zero_odes)
1551-
error_if_constraint_odes(NonlinearSystem, fullrs)
15521556
eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
15531557

15541558
NonlinearSystem(eqs, sts, ps;
@@ -1669,12 +1673,20 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
16691673
check_length = false, name = nameof(rs),
16701674
combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
16711675
include_zero_odes = true, remove_conserved = false,
1672-
checks = false, kwargs...)
1676+
checks = false, structural_simplify=false, kwargs...)
16731677
u0map = symmap_to_varmap(rs, u0)
16741678
pmap = symmap_to_varmap(rs, p)
16751679
osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
16761680
remove_conserved)
16771681
osys = complete(osys)
1682+
1683+
# Handles potential Differential algebraic equations.
1684+
if structural_simplify
1685+
(osys = MT.structural_simplify(osys))
1686+
# elseif has_alg_equations(rs)
1687+
# error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.")
1688+
end
1689+
16781690
return ODEProblem(osys, u0map, tspan, pmap, args...; check_length, kwargs...)
16791691
end
16801692

0 commit comments

Comments
 (0)