init

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

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](@id constraint_equations_coupling_constratins)
+## [Coupling ODE constraints via extending a system](@id constraint_equations_coupling_constraints)
 
 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

@@ -701,3 +701,65 @@ rn = @reaction_network begin
 end
 nothing # hide
 ```
+
+## 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

@@ -86,6 +86,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

src/expression_utils.jl

@@ -11,6 +11,13 @@ function get_tup_arg(ex::ExprValues, i::Int)
     return ex.args[i]
 end
 
+# Some options takes input on form that is either `@option ...` or `@option begin ... end`. 
+# This transforms input of the latter form to the former (with only one line in the `begin ... end` block)
+function option_block_form(expr)
+    (expr.head == :block) && return expr
+    return Expr(:block, expr)
+end
+
 # In variable/species/parameters on the forms like:
 # X
 # X = 1.0

src/networkapi.jl

@@ -97,6 +97,55 @@ function reactions(network)
     [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
 end
 
+"""
+    diff_equations(network)
+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).
+Notes:
+- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
+"""
+function diff_equations(network)
+    eqs = equations(network)
+    filter!(!isreaction, eqs)
+    systems = filter_nonrxsys(network)
+    isempty(systems) && (return rxs)
+    [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
+end
+
+"""
+    has_diff_equations(network)
+Given a [`ReactionSystem`](@ref), check whether it contain any differential equations (i.e. in addition to those generated through reactions).
+Notes:
+- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
+"""
+function has_diff_equations(network)
+    return !isempty(diff_equations(network))
+end
+
+"""
+    alg_equations(network)
+Given a [`ReactionSystem`](@ref), return a vector of all `Equations` in the system that are algebraic equations (does not contain any derivatives).
+Notes:
+- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
+"""
+function alg_equations(network)
+    eqs = equations(network)
+    filter!(!isreaction, eqs)
+    filter!(!isreaction, eqs)
+    systems = filter_nonrxsys(network)
+    isempty(systems) && (return rxs)
+    [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
+end
+
+"""
+    has_alg_equations(network)
+Given a [`ReactionSystem`](@ref), check whether it contain any algebraic equations.
+Notes:
+- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
+"""
+function has_alg_equations(network)
+    return !isempty(alg_equations(network))
+end
+
 """
     speciesmap(network)
 

src/reaction_network.jl

@@ -82,7 +82,8 @@ const forbidden_variables_error = let
 end
 
 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling,
+                     :differentials, :equations, :continuous_events, :discrete_events)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -353,20 +354,28 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     option_lines = Expr[x for x in ex.args if x.head == :macrocall]
 
     # Get macro options.
-    length(unique(arg.args[1] for arg in option_lines)) < length(option_lines) && error("Some options where given multiple times.")
+    if length(unique(arg.args[1] for arg in option_lines)) < length(option_lines)
+         error("Some options where given multiple times.")
+    end
     options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
                        option_lines))
 
     # Reads options.
     default_reaction_metadata = :([])
     check_default_noise_scaling!(default_reaction_metadata, options)
     compound_expr, compound_species = read_compound_options(options)
+    continuous_events_expr = read_events_option(options, :continuous_events)
+    discrete_events_expr = read_events_option(options, :discrete_events)
 
     # Parses reactions, species, and parameters.
     reactions = get_reactions(reaction_lines)
     species_declared = [extract_syms(options, :species); compound_species]
     parameters_declared = extract_syms(options, :parameters)
-    variables = extract_syms(options, :variables)
+    variables_declared = extract_syms(options, :variables)
+
+    # Reads more options.
+    vars_extracted, add_default_diff, equations = read_equations_options(options, variables_declared)
+    variables = vcat(variables_declared, vars_extracted)
 
     # handle independent variables
     if haskey(options, :ivs)
@@ -391,6 +400,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     species = vcat(species_declared, species_extracted)
     parameters = vcat(parameters_declared, parameters_extracted)
 
+    # Create differential expression.
+    diffexpr = create_differential_expr(options, add_default_diff, [species; parameters; variables])
+
     # Checks for input errors.
     (sum(length.([reaction_lines, option_lines])) != length(ex.args)) &&
         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)))
 
     # Creates expressions corresponding to actual code from the internal DSL representation.
     sexprs = get_sexpr(species_extracted, options; iv_symbols = ivs)
-    vexprs = haskey(options, :variables) ? options[:variables] : :()
+    vexprs = get_sexpr(vars_extracted, options, :variables)
     pexprs = get_pexpr(parameters_extracted, options)
     ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps")
     vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars")
@@ -412,6 +424,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     for reaction in reactions
         push!(rxexprs.args, get_rxexprs(reaction))
     end
+    for equation in equations
+        push!(rxexprs.args, equation)
+    end
 
     quote
         $ps
@@ -420,11 +435,13 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         $sps
         $observed_vars
         $comps
+        $diffexpr
 
         Catalyst.remake_ReactionSystem_internal(
             Catalyst.make_ReactionSystem_internal(
                 $rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms),
-                $pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs);
+                $pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs,
+                continuous_events = $continuous_events_expr, discrete_events = $discrete_events_expr);
             default_reaction_metadata = $default_reaction_metadata
         )
     end
@@ -810,6 +827,53 @@ function check_default_noise_scaling!(default_reaction_metadata, options)
     end
 end
 
+# Creates an expression declaring differentials.
+function create_differential_expr(options, add_default_diff, used_syms)
+    # Creates the differential expression.
+    # If differentials was provided as options, this is used as the initial expression.
+    # If the default differential (D(...)) was used in equations, this is added to the expression.
+    diffexpr = (haskey(options, :differentials) ? options[:differentials].args[3] : MacroTools.striplines(:(begin end)))
+    diffexpr = option_block_form(diffexpr)
+    add_default_diff && push!(diffexpr.args, :(D = Differential($(DEFAULT_IV_SYM))))
+
+    # Goes through all differentials, checking that they are correctly formatted and their symbol is not used elsewhere.
+    for dexpr in diffexpr.args
+        (dexpr.head != :(=)) && error("Differential declaration must have form like D = Differential(t), instead \"$(dexpr)\" was given.")
+        (dexpr.args[1] isa Symbol) || error("Differential left-hand side must be a single symbol, instead \"$(dexpr.args[1])\" was given.")
+        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.")
+        in(dexpr.args[1], forbidden_symbols_error) && error("A forbidden symbol ($(dexpr.args[1])) was used as a differential name.")
+    end
+
+    return diffexpr
+end
+
+# Read the events (continious or discrete) provided as options to the DSL. Returns an expression which evalutes to these.
+function read_events_option(options, event_type::Symbol)    
+    # Prepares the events, if required to, converts them to block form.
+    (event_type in [:continuous_events]) || error("Trying to read an unsupported event type.")
+    events_input = haskey(options, event_type) ? options[event_type].args[3] : :(begin end)
+    events_input = option_block_form(events_input)
+
+    # Goes throgh the events, checks for errors, and adds them to the output vector.
+    events_expr = :([])
+    for arg in events_input.args
+        # Formatting error checks.
+        # NOTE: Maybe we should move these deeper into the system (rather than the DSL), throwing errors more generally?
+        if (arg.head != :call) || (arg.args[1] != :(=>)) || length(arg.args) != 3
+            error("Events should be on form `condition => affect`, separated by a `=>`. This appears not to be the case for: $(arg).")
+        end
+        if (arg.args[2] != :call) && (event_type == :continuous_events)
+            error("The condition part of continious events (the left-hand side) must be a vector. This is not the case for: $(arg).")
+        end
+        if arg.args[3] != :call
+             error("The affect part of all events (the righ-hand side) must be a vector. This is not the case for: $(arg).")
+        end
+
+        push!(events_expr.args, arg)
+    end
+    return events_expr
+end
+
 ### Functionality for expanding function call to custom and specific functions ###
 
 #Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.

src/reactionsystem.jl

@@ -546,6 +546,11 @@ struct ReactionSystem{V <: NetworkProperties} <:
                             name, systems, defaults, connection_type, nps, cls, cevs, devs,
                             metadata, complete; checks::Bool = true)
 
+        # Filters away any potential obervables from `states` and `spcs`.
+        obs_vars = [obs_eq.lhs for obs_eq in observed]
+        states = filter(state -> !any(isequal(state, obs_var) for obs_var in obs_vars), states)
+        spcs = filter(spc -> !any(isequal(spc, obs_var) for obs_var in obs_vars), spcs)
+
         # unit checks are for ODEs and Reactions only currently
         nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
         if checks && isempty(sivs)
@@ -1409,7 +1414,7 @@ end
 
 # merge constraint components with the ReactionSystem components
 # also handles removing BC and constant species
-function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false)
+function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false, zero_derivatives = false)
     # if there are BC species, put them after the independent species
     rssts = get_unknowns(rs)
     sts = any(isbc, rssts) ? vcat(ists, filter(isbc, rssts)) : ists
@@ -1501,7 +1506,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
     ists, ispcs = get_indep_sts(fullrs, remove_conserved)
     eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
                          include_zero_odes)
-    eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
+    eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved, zero_derivatives=true)
 
     ODESystem(eqs, get_iv(fullrs), sts, ps;
               observed = obs,
@@ -1542,7 +1547,6 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     ists, ispcs = get_indep_sts(fullrs, remove_conserved)
     eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
                          as_odes = false, include_zero_odes)
-    error_if_constraint_odes(NonlinearSystem, fullrs)
     eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
 
     NonlinearSystem(eqs, sts, ps;
@@ -1663,12 +1667,20 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
                                check_length = false, name = nameof(rs),
                                combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                                include_zero_odes = true, remove_conserved = false,
-                               checks = false, kwargs...)
+                               checks = false, structural_simplify=false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
                    remove_conserved)
     osys = complete(osys)
+
+    # Handles potential Differential algebraic equations.
+    if structural_simplify 
+        (osys = MT.structural_simplify(osys))
+    # elseif has_alg_equations(rs)
+    #     error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.")
+    end
+
     return ODEProblem(osys, u0map, tspan, pmap, args...; check_length, kwargs...)
 end