Equations and events reupload

Project.toml

@@ -6,6 +6,7 @@ version = "13.5.1"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
@@ -46,20 +47,21 @@ JumpProcesses = "9.3.2"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "9.7.1"
+ModelingToolkit = "9.10.0"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
 Setfield = "1"
 StructuralIdentifiability = "0.5.1"
 SymbolicUtils = "1.0.3"
-Symbolics = "5.14"
+Symbolics = "5.27"
 Unitful = "1.12.4"
 julia = "1.9"
 
 [extras]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
@@ -80,4 +82,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["BifurcationKit", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"]
+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

@@ -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,66 @@ rn = @reaction_network begin
 end
 nothing # hide
 ```
+
+## Incorporating (differential) equations into reaction network models
+Some models cannot be purely described as reaction networks. 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 parameter using the `@parameters` option:
+```@example eqs1
+rn = @reaction_network begin
+  @parameters k
+  @equations begin
+    D(V) ~ k*G
+  end
+  (p,d), 0 <--> G
+end
+nothing #hide
+```
+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).
+
+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
+```

ext/CatalystHomotopyContinuationExtension.jl

@@ -2,6 +2,7 @@ module CatalystHomotopyContinuationExtension
 
 # Fetch packages.
 using Catalyst
+import DynamicPolynomials
 import ModelingToolkit as MT
 import HomotopyContinuation as HC
 import Setfield: @set

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -53,7 +53,8 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     eqs_pars_funcs = vcat(equations(ns), conservedequations(rs))
     eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs)
     eqs_intexp = make_int_exps.(eqs)
-    return Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
+    ss_poly = Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp))
+    return poly_type_convert(ss_poly)
 end
 
 # If u0s are not given while conservation laws are present, throws an error.
@@ -94,7 +95,7 @@ end
 function reorder_sols!(sols, ss_poly, rs::ReactionSystem)
     var_names_extended = String.(Symbol.(HC.variables(ss_poly)))
     var_names = [Symbol(s[1:prevind(s, findlast('_', s))]) for s in var_names_extended]
-    sort_pattern = indexin(MT.getname.(species(rs)), var_names)
+    sort_pattern = indexin(MT.getname.(unknowns(rs)), var_names)
     foreach(sol -> permute!(sol, sort_pattern), sols)
 end
 
@@ -104,4 +105,13 @@ function filter_negative_f(sols; neg_thres=-1e-20)
         (neg_thres < sol[idx] < 0) && (sol[idx] = 0)
     end
     return filter(sol -> all(>=(0), sol), sols)
+end
+
+# Sometimes (when polynomials are created from hybrid CRN/DAEs), the steady state polynomial have the wrong type.
+# This converts it to the correct type, which homotopy continuation can handle.
+const WRONG_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}}}
+const CORRECT_POLY_TYPE = Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, DynamicPolynomials.Graded{DynamicPolynomials.LexOrder}, Float64}}
+function poly_type_convert(ss_poly)
+    (typeof(ss_poly) == WRONG_POLY_TYPE) && return convert(CORRECT_POLY_TYPE, ss_poly)
+    return ss_poly
 end

src/Catalyst.jl

@@ -37,11 +37,15 @@ import ModelingToolkit: check_variables,
                         check_parameters, _iszero, _merge, check_units,
                         get_unit, check_equations, iscomplete
 
+# Event-related MTK stuff.
+import ModelingToolkit: PresetTimeCallback
+
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
 import MacroTools, Graphs
 import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne
 import DataStructures: OrderedDict, OrderedSet
 import Parameters: @with_kw_noshow
+import Symbolics: occursin, wrap
 
 # globals for the modulate
 function default_time_deriv()

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/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
@@ -799,6 +816,100 @@ function make_observed_eqs(observables_expr)
     end 
 end
 
+# Reads the variables options. Outputs:
+# `vars_extracted`: A vector with extracted variables (lhs in pure differential equations only).
+# `dtexpr`: If a differentialequation is defined, the default derrivative (D ~ Differential(t)) must be defined.
+# `equations`: a vector with the equations provided.
+function read_equations_options(options, variables_declared)
+    # Prepares the equations. First, extracts equations from provided option (converting to block form if requried).
+    # Next, uses MTK's `parse_equations!` function to split input into a vector with the equations.
+    eqs_input = haskey(options, :equations) ? options[:equations].args[3] : :(begin end)
+    eqs_input = option_block_form(eqs_input)
+    equations = Expr[]
+    ModelingToolkit.parse_equations!(Expr(:block), equations, Dict{Symbol, Any}(), eqs_input)
+
+    # Loops through all equations, checks for lhs of the form `D(X) ~ ...`.
+    # When this is the case, the variable X and differential D are extracted (for automatic declaration).
+    # Also performs simple error checks.
+    vars_extracted = []
+    add_default_diff = false
+    for eq in equations
+        if (eq.head != :call) || (eq.args[1] != :~)
+             error("Malformed equation: \"$eq\". Equation's left hand and right hand sides should be separated by a \"~\".")
+        end
+
+        # Checks if the equation have the format D(X) ~ ... (where X is a symbol). This means that the 
+        # default differential has been used. X is added as a declared variable to the system, and
+        # we make a note that a differential D = Differential(iv) should be made as well.
+        lhs = eq.args[2]
+        # if lhs: is an expression. Is a function call. The function's name is D. Calls a single symbol.
+        if (lhs isa Expr) && (lhs.head == :call) && (lhs.args[1] == :D) && (lhs.args[2] isa Symbol)
+            diff_var = lhs.args[2]
+            if in(diff_var, forbidden_symbols_error)
+                error("A forbidden symbol ($(diff_var)) was used as an variable in this differential equation: $eq")
+            end
+            add_default_diff = true
+            in(diff_var, variables_declared) || push!(vars_extracted, diff_var)
+        end
+    end
+
+    return vars_extracted, add_default_diff, equations
+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)
+
+    # 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
+
+    # If the default differential D has been used, but not pre-declared using the @differenitals
+    # options, add this declaration to the list of declared differentials.
+    if add_default_diff && !any(diff_dec.args[1] == :D for diff_dec in diffexpr.args) 
+        push!(diffexpr.args, :(D = Differential($(DEFAULT_IV_SYM))))
+    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, :discrete_events]) || error("Trying to read an unsupported event type.")
+    events_input = haskey(options, event_type) ? options[event_type].args[3] : MacroTools.striplines(:(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 isa Expr) && (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 isa Expr) && (arg.args[2] isa Expr) && (arg.args[2].head != :vect) && (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 isa Expr) && (arg.args[3] isa Expr) && (arg.args[3].head != :vect)
+             error("The affect part of all events (the righ-hand side) must be a vector. This is not the case for: $(arg).")
+        end
+
+        # Adds the correctly formatted event to the event creation expression.
+        push!(events_expr.args, arg)
+    end
+
+    return events_expr
+end
+
 # Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a
 # default metadata value to the `default_reaction_metadata` vector.
 function check_default_noise_scaling!(default_reaction_metadata, options)