update constraint eq tutorial

Author: isaacsas

docs/src/model_creation/constraint_equations.md

@@ -21,7 +21,9 @@ The easiest way to include ODEs and algebraic equations is to just include them
 when using the DSL to specify a model. Here we include an ODE for $V(t)$ along
 with degradation and production reactions for $P(t)$:
 ```@example ceq1
-rn = @reaction_network dynamic_cell begin
+using Catalyst, OrdinaryDiffEqTsit5, Plots
+
+rn = @reaction_network growing_cell begin
     # the growth rate
     @parameters λ = 1.0
 
@@ -44,7 +46,7 @@ end
 We can now create an `ODEProblem` from our model and solve it to see how $V(t)$
 and $P(t)$ evolve in time:
 ```@example ceq1
-oprob = ODEProblem(growing_cell, [], (0.0, 1.0))
+oprob = ODEProblem(rn, [], (0.0, 1.0))
 sol = solve(oprob, Tsit5())
 plot(sol)
 ```
@@ -54,6 +56,7 @@ As an alternative to the previous approach, we could have also constructed our
 `ReactionSystem` all at once using the symbolic interface:
 ```@example ceq2
 using Catalyst, OrdinaryDiffEqTsit5, Plots
+
 t = default_t()
 D = default_time_deriv()
 
@@ -84,8 +87,9 @@ equations), we must import the time differential operator from Catalyst. We do
 this through `D = default_time_deriv()`. Here, `D(V)` denotes the differential
 of the variable `V` with respect to time.
 
-```@example ceq1
+```@example ceq2b
 using Catalyst, OrdinaryDiffEqTsit5, Plots
+
 t = default_t()
 D = default_time_deriv()
 
@@ -116,14 +120,14 @@ systems together Catalyst requires that the systems have not been marked as
 
 We can now merge the two systems into one complete `ReactionSystem` model using
 [`ModelingToolkit.extend`](@ref):
-```@example ceq1
+```@example ceq2b
 @named growing_cell = extend(osys, rn)
 growing_cell = complete(growing_cell)
 ```
 
 We see that the combined model now has both the reactions and ODEs as its
 equations. To solve and plot the model we proceed like normal
-```@example ceq1
+```@example ceq2b
 oprob = ODEProblem(growing_cell, [], (0.0, 1.0))
 sol = solve(oprob, Tsit5())
 plot(sol)
@@ -143,7 +147,82 @@ details on the types of events that can be represented symbolically. A
 lower-level approach for creating events via the DifferentialEquations.jl
 callback interface is illustrated [here](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) tutorial.
 
-Let's first create our equations and unknowns/species again
+```@example ceq3a
+using Catalyst, OrdinaryDiffEqTsit5, Plots
+
+rn = @reaction_network growing_cell begin
+    # the growth rate
+    @parameters λ = 1.0
+
+    # assume there is no protein initially
+    @species P(t) = 0.0
+
+    # set the initial volume to 1.0
+    @variables V(t) = 1.0
+
+    # the reactions
+    V, 0 --> P
+    1.0, P --> 0
+
+    # the coupled ODE for V(t)
+    @equations begin
+        D(V) ~ λ * V
+    end
+
+    # every 1.0 time unit we half the volume of the cell and the number of proteins
+    @continuous_events begin 
+        [V ~ 2.0] => [V ~ V/2, P ~ P/2]
+    end
+end
+```
+We can now create and simulate our model
+```@example ceq3a
+oprob = ODEProblem(rn, [], (0.0, 10.0))
+sol = solve(oprob, Tsit5())
+plot(sol)
+Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide
+```
+We can also model discrete events. Here at a time `switch_time` we will set the parameter `k_on` to be
+zero:
+```@example ceq3a
+rn = @reaction_network param_off_ex begin
+    @parameters switch_time
+    k_on, A --> B
+    k_off, B --> A
+    
+    @discrete_events begin
+        (t == switch_time) => [k_on ~ 0.0]
+    end
+end
+
+u0 = [:A => 10.0, :B => 0.0]
+tspan = (0.0, 4.0)
+p = [:k_on => 100.0, :switch_time => 2.0, :k_off => 10.0]
+oprob = ODEProblem(rn, u0, tspan, p)
+sol = solve(oprob, Tsit5(); tstops = 2.0)
+plot(sol)
+```
+Note that for discrete events we need to set a stop time via `tstops` so that
+the ODE solver can step exactly to the specific time of our event. In the
+previous example we just manually set the numeric value of the parameter in the
+`tstops` kwarg to `solve`, however, it can often be convenient to instead get
+the value of the parameter from `oprob` and pass this numeric value. This helps
+ensure consistency between the value passed via `p` and/or symbolic defaults and
+what we pass as a `tstop` to `solve`. We can do this as
+```@example ceq3a
+oprob = ODEProblem(rn, u0, tspan, p)
+switch_time_val = oprob.ps[:switch_time]
+sol = solve(oprob, Tsit5(); tstops = switch_time_val)
+plot(sol)
+```
+For a detailed discussion on how to directly use the lower-level but more
+flexible DifferentialEquations.jl event/callback interface, see the
+[tutorial](https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Event-handling-using-callbacks)
+on event handling using callbacks.
+
+## [Adding events via the symbolic interface](@id constraint_equations_events_symbolic)
+Let's repeat the previous two models using the symbolic interface. We first
+create our equations and unknowns/species again
 ```@example ceq3
 using Catalyst, OrdinaryDiffEqTsit5, Plots
 t = default_t()
@@ -171,7 +250,9 @@ sol = solve(oprob, Tsit5())
 plot(sol)
 Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide
 ```
-We can also model discrete events. Similar to our example with continuous events, we start by creating reaction equations, parameters, variables, and unknowns.
+We can again also model discrete events. Similar to our example with continuous
+events, we start by creating reaction equations, parameters, variables, and
+unknowns.
 ```@example ceq3
 t = default_t()
 @parameters k_on switch_time k_off
@@ -193,22 +274,7 @@ Simulating our model,
 rs2 = complete(rs2)
 
 oprob = ODEProblem(rs2, u0, tspan, p)
-sol = solve(oprob, Tsit5(); tstops = 2.0)
-plot(sol)
-```
-Note that for discrete events we need to set a stop time via `tstops` so that
-the ODE solver can step exactly to the specific time of our event. In the
-previous example we just manually set the numeric value of the parameter in the
-`tstops` kwarg to `solve`, however, it can often be convenient to instead get
-the value of the parameter from `oprob` and pass this numeric value. This helps
-ensure consistency between the value passed via `p` and/or symbolic defaults and
-what we pass as a `tstop` to `solve`. We can do this as
-```julia
 switch_time_val = oprob.ps[:switch_time]
 sol = solve(oprob, Tsit5(); tstops = switch_time_val)
 plot(sol)
 ```
-For a detailed discussion on how to directly use the lower-level but more
-flexible DifferentialEquations.jl event/callback interface, see the
-[tutorial](https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Event-handling-using-callbacks)
-on event handling using callbacks.