From fcd64aa0bcd8899a2303a3ee653ea7f0e9800406 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 27 Jul 2024 15:01:25 -0400
Subject: [PATCH 1/4] init

---
 src/reactionsystem_conversions.jl | 549 ++++++++++++++++++------------
 1 file changed, 326 insertions(+), 223 deletions(-)

diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index f2c83507fc..ee4915da2d 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -1,7 +1,7 @@
 ### ODE & SDE Assembly ###
 
 """
-    oderatelaw(rx; combinatoric_ratelaw=true)
+ oderatelaw(rx; combinatoric_ratelaw = true)
 
 Given a [`Reaction`](@ref), return the symbolic reaction rate law used in
 generated ODEs for the reaction. Note, for a reaction defined by
@@ -17,21 +17,26 @@ the expression that is returned will be `k * (X(t)^2/2) * (Y(t)^3/6)`.
 
 Notes:
 - Allocates
-- `combinatoric_ratelaw=true` uses factorial scaling factors in calculating the
-    rate law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
-    `combinatoric_ratelaw=false` then the ratelaw is `k*S^2`, i.e. the scaling
-    factor is ignored.
+- `combinatoric_ratelaw = true` uses factorial scaling factors in calculating the
+ rate law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
+ `combinatoric_ratelaw = false` then the ratelaw is `k*S^2`, i.e. the scaling
+ factor is ignored.
 """
 function oderatelaw(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, only_use_rate = rx
     rl = rate
 
-    # if the stoichiometric coefficients are not integers error if asking to scale rates
+    # if the stoichiometric coefficients are not integers, error if asking to scale rates
     !all(s -> s isa Union{Integer, Symbolic}, substoich) &&
         (combinatoric_ratelaw == true) &&
-        error("Non-integer stoichiometric coefficients require the combinatoric_ratelaw=false keyword to oderatelaw, or passing combinatoric_ratelaws=false to convert or ODEProblem.")
+        error("Non-integer stoichiometric coefficients require the combinatoric_ratelaw = false keyword to oderatelaw, or passing combinatoric_ratelaws = false to convert or ODEProblem.")
 
+    # If we use mass action, update the rate law according to the reaction substrates.
     if !only_use_rate
+        # `coef` stores the number bit (i.e. division by factorial of stoichiometries).
+        # Symbolic bit (i.e. multiplication by substrate concentrations) are appended directly 
+        # to the rate law (`rl`). Cases with stoichiometries 1 are handled more efficiently
+        # using separate routines. 
         coef = eltype(substoich) <: Number ? one(eltype(substoich)) : 1
         for (i, stoich) in enumerate(substoich)
             combinatoric_ratelaw && (coef *= factorial(stoich))
@@ -44,21 +49,32 @@ end
 
 # Function returning `true` for species which shouldn't change from the reactions, 
 # including non-species variables.
-drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s))
+drop_dynamics(sp) = isconstant(sp) || isbc(sp) || (!isspecies(sp))
 
-function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserved = false)
+# For a `ReactionSystem`, create a vector with all the right-hand side expressions in the
+# corresponding ODE system.
+function assemble_oderhs(rs, isps; combinatoric_ratelaws = true, remove_conserved = false)
+    # Initiates the right-hand side equations as expressions with `0` only.
+    species_to_idx = Dict(x => i for (i, x) in enumerate(isps))
+    rhsvec = Any[0 for _ in isps]
+
+    # If we have eliminated conservation laws, create a dictionary taking each species eliminated
+    # in this process to the expression it is replaced by (if not, creates an empty dict).
     nps = get_networkproperties(rs)
-    species_to_idx = Dict(x => i for (i, x) in enumerate(ispcs))
-    rhsvec = Any[0 for _ in ispcs]
     depspec_submap = if remove_conserved
         Dict(eq.lhs => eq.rhs for eq in nps.conservedeqs)
     else
         Dict()
     end
 
+    # Loops through all reactions, for each, updates relevant right-hand side equations.
     for rx in get_rxs(rs)
+        # Computes the reaction's rate late. If conserved quantities have been removed, replaces
+        # these in the rate law.
         rl = oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws)
         remove_conserved && (rl = substitute(rl, depspec_submap))
+
+        # For each species which participates in the reaction, update its equation.
         for (spec, stoich) in rx.netstoich
             # dependent species don't get an ODE, so are skipped
             remove_conserved && (spec in nps.depspecs) && continue
@@ -67,10 +83,14 @@ function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserv
             drop_dynamics(spec) && continue
 
             i = species_to_idx[spec]
+            # If the relevant RHS equation is 0 (probably because it has not been
+            # modified yet), we use a different routine (compared to if it is non-zero).
             if _iszero(rhsvec[i])
                 if stoich isa Symbolic
                     rhsvec[i] = stoich * rl
                 else
+                    # If the stoichiometry is +/- 1, initiate using the (potentially negative)
+                    # rate law directly. If not, initiate with rate law * stoichiometry.
                     signedrl = (stoich > zero(stoich)) ? rl : -rl
                     rhsvec[i] = isone(abs(stoich)) ? signedrl : stoich * rl
                 end
@@ -78,6 +98,8 @@ function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserv
                 if stoich isa Symbolic
                     rhsvec[i] += stoich * rl
                 else
+                    # If the stoichiometry is negative, we subtract the rate law (multiplied by the
+                    # absolute value of the stoichiometry) from the RHS. If not, add it.
                     Δspec = isone(abs(stoich)) ? rl : abs(stoich) * rl
                     rhsvec[i] = (stoich > zero(stoich)) ? (rhsvec[i] + Δspec) :
                                 (rhsvec[i] - Δspec)
@@ -89,29 +111,36 @@ function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserv
     rhsvec
 end
 
-function assemble_drift(rs, ispcs; combinatoric_ratelaws = true, as_odes = true,
-        include_zero_odes = true, remove_conserved = false)
-    rhsvec = assemble_oderhs(rs, ispcs; combinatoric_ratelaws, remove_conserved)
+# For a `ReactionSystem`, creates a vector with all the differential equations corresponding to
+# the ODEs generated by the reaction rate equation. The `as_odes` determine whether the LHS should
+# be differential (when generating e.g. a `ODESystem`) or `0` (when generating e.g. a NonlinearSystem).
+function assemble_drift(rs, ind_sps; combinatoric_ratelaws = true, as_odes = true,
+                        include_zero_odes = true, remove_conserved = false)
+    rhsvec = assemble_oderhs(rs, ind_sps; combinatoric_ratelaws, remove_conserved)
     if as_odes
         D = Differential(get_iv(rs))
         eqs = [Equation(D(x), rhs)
-               for (x, rhs) in zip(ispcs, rhsvec) if (include_zero_odes || (!_iszero(rhs)))]
+               for (x, rhs) in zip(ind_sps, rhsvec) if (include_zero_odes || (!_iszero(rhs)))]
     else
         eqs = [Equation(0, rhs) for rhs in rhsvec if (include_zero_odes || (!_iszero(rhs)))]
     end
     eqs
 end
 
-# this doesn't work with constraint equations currently
-function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true,
-        remove_conserved = false)
+# this doesn't work with coupled equations currently
+# Assemblies the diffusion equations terms used in the chemical Langevin equation SDE.
+function assemble_diffusion(rs, us, ind_sps; combinatoric_ratelaws = true,
+                            remove_conserved = false)
     # as BC species should ultimately get an equation, we include them in the noise matrix
-    num_bcsts = count(isbc, get_unknowns(rs))
+    num_bcus = count(isbc, get_unknowns(rs))
 
     # we make a matrix sized by the number of reactions
     eqs = Matrix{Any}(undef, length(sts) + num_bcsts, length(get_rxs(rs)))
     eqs .= 0
     species_to_idx = Dict((x => i for (i, x) in enumerate(ispcs)))
+
+    # If we have eliminated conservation laws, create a dictionary taking each species eliminated
+    # in this process to the expression it is replaced by (if not, creates an empty dict).
     nps = get_networkproperties(rs)
     depspec_submap = if remove_conserved
         Dict(eq.lhs => eq.rhs for eq in nps.conservedeqs)
@@ -119,11 +148,15 @@ function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true,
         Dict()
     end
 
+    # Loops through all reactions, for each, updates relevant noise terms in the noise matrix.
     for (j, rx) in enumerate(get_rxs(rs))
+        # Computes the reactions (noise) rate late. This is the sqrt of the rate late, multiplied
+        # by any noise scaling term. If conserved quantities have been removed, these are substituted in.
         rlsqrt = sqrt(abs(oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws)))
         hasnoisescaling(rx) && (rlsqrt *= getnoisescaling(rx))
         remove_conserved && (rlsqrt = substitute(rlsqrt, depspec_submap))
 
+        # For each species which participate in the reaction, updates its corresponding noise terms.
         for (spec, stoich) in rx.netstoich
             # dependent species don't get an equation
             remove_conserved && (spec in nps.depspecs) && continue
@@ -132,6 +165,8 @@ function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true,
             drop_dynamics(spec) && continue
 
             i = species_to_idx[spec]
+            # Differentiates if the stoichiometry is a symbolic (e.g. a parameter). If not it might
+            # be +/- 1, in which case we initiate the term differently.
             if stoich isa Symbolic
                 eqs[i, j] = stoich * rlsqrt
             else
@@ -146,7 +181,7 @@ end
 ### Jumps Assembly ###
 
 """
-    jumpratelaw(rx; combinatoric_ratelaw=true)
+ jumpratelaw(rx; combinatoric_ratelaw = true)
 
 Given a [`Reaction`](@ref), return the symbolic reaction rate law used in
 generated stochastic chemical kinetics model SSAs for the reaction. Note,
@@ -164,42 +199,58 @@ binomial(Y,3)`.
 
 Notes:
 - Allocates
-- `combinatoric_ratelaw=true` uses binomials in calculating the rate law, i.e. for `2S ->
-  0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaw=false` then
-  the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling
-  factor.
+- `combinatoric_ratelaw = true` uses binomials in calculating the rate law, i.e. for `2S ->
+ 0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaw = false` then
+ the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling
+ factor.
 """
 function jumpratelaw(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, only_use_rate = rx
     rl = rate
+
+    # If the `only_use_rate` is unused, update the rate according to the stoichiometries.
     if !only_use_rate
+        # Initiates the coefficient as "1", using the same type as in `substoich`.
         coef = eltype(substoich) <: Number ? one(eltype(substoich)) : 1
+
+        # For each substrate, update the rate law (and coefficient) terms according to the
+        # substrate's stoichiometry.
         for (i, stoich) in enumerate(substoich)
-            s = substrates[i]
+            spec = substrates[i]
+
+            # For symbolic substrates (e.g. parameters), simply update the rate law with the
+            # appropriate expression. If not, we  can simplify the expression (and also treat
+            # the case with stoichiometry = 1 separately).
             if stoich isa Symbolic
-                rl *= combinatoric_ratelaw ? binomial(s, stoich) :
-                      factorial(s) / factorial(s - stoich)
+                rl *= combinatoric_ratelaw ? binomial(spec, stoich) :
+                      factorial(spec) / factorial(spec - stoich)
             else
-                rl *= s
+                rl *= spec
+                # Since `stoich` is a number, (if it is not 1) we add the binomial's denominator
+                # to the `coef` (which is a number) directly. For Symbolic stoichiometry everything
+                # is instead added to the rate law.
                 isone(stoich) && continue
                 for i in one(stoich):(stoich - one(stoich))
-                    rl *= (s - i)
+                    rl *= (spec - i)
                 end
                 combinatoric_ratelaw && (coef *= factorial(stoich))
             end
         end
+
+        # If we use a combinatoric rate law, we have generated a numeric coefficient which we
+        # will (if it is not 1) divde our rate law with (to generate the full rate law).
         combinatoric_ratelaw && !isequal(coef, one(coef)) && (rl /= coef)
     end
     rl
 end
 
-# if haveivdep=false then time dependent rates will still be classified as mass action
+# if haveivdep = false then time-dependent rates will still be classified as mass action
 """
 ```julia
 ismassaction(rx, rs; rxvars = get_variables(rx.rate),
-                              haveivdep = nothing,
-                              unknownset = Set(unknowns(rs)),
-                              ivset = nothing)
+ haveivdep = nothing,
+ unknownset = Set(unknowns(rs)),
+ ivset = nothing)
 ```
 
 True if a given reaction is of mass action form, i.e. `rx.rate` does not depend
@@ -210,25 +261,25 @@ depend explicitly on the independent variable (usually time).
 - `rx`, the [`Reaction`](@ref).
 - `rs`, a [`ReactionSystem`](@ref) containing the reaction.
 - Optional: `rxvars`, `Variable`s which are not in `rxvars` are ignored as
-  possible dependencies.
+ possible dependencies.
 - Optional: `haveivdep`, `true` if the [`Reaction`](@ref) `rate` field
-  explicitly depends on any independent variable (i.e. t or for spatial systems x,y,etc).
-  If not set, will be automatically calculated.
+ explicitly depends on any independent variable (i.e. t or for spatial systems x,y,etc).
+ If not set, will be automatically calculated.
 - Optional: `unknownset`, set of unknowns which if the rxvars are within mean rx is
-  non-mass action.
+ non-mass action.
 - Optional: `ivset`, a `Set` of the independent variables of the system. If not provided and
-  the system is spatial, i.e. `isspatial(rs) == true`, it will be created with all the
-  spatial variables and the time variable. If the rate expression contains any element of
-  `ivset`, then `ismassaction(rx,rs) == false`. Pass a custom set to control this behavior.
+ the system is spatial, i.e. `isspatial(rs) == true`, it will be created with all the
+ spatial variables and the time variable. If the rate expression contains any element of
+ `ivset`, then `ismassaction(rx,rs) == false`. Pass a custom set to control this behavior.
 
 Notes:
 - Non-integer stoichiometry is treated as non-mass action. This includes
-  symbolic variables/terms or floating point numbers for stoichiometric
-  coefficients.
+ symbolic variables/terms or floating point numbers for stoichiometric
+ coefficients.
 """
 function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
-        haveivdep::Union{Nothing, Bool} = nothing,
-        unknownset = Set(get_unknowns(rs)), ivset = nothing)
+    haveivdep::Union{Nothing, Bool} = nothing,
+    unknownset = Set(get_unknowns(rs)), ivset = nothing)
 
     # we define non-integer (i.e. float or symbolic) stoich to be non-mass action
     ((eltype(rx.substoich) <: Integer) && (eltype(rx.prodstoich) <: Integer)) ||
@@ -237,6 +288,8 @@ function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
     # if no dependencies must be zero order
     (length(rxvars) == 0) && return true
 
+    # Reactions that depend on independent variables are not mass action. If iv dependency
+    # was not marked when `ismassaction` was called, we check if the reaction depend on any ivs.
     if (haveivdep === nothing)
         if isspatial(rs)
             if (ivset === nothing)
@@ -249,22 +302,32 @@ function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
         else
             ivdep = any(var -> isequal(get_iv(rs), var), rxvars)
         end
-        ivdep && return false
+            ivdep && return false
     else
         haveivdep && return false
     end
+
+    # The `only_use_rate` option disables mass action for the reaction (=> not mass action reaction).
     rx.only_use_rate && return false
+
+    # not mass action if have a non-constant, variable species in the rate expression
     @inbounds for var in rxvars
-        # not mass action if have a non-constant, variable species in the rate expression
         (var in unknownset) && return false
     end
 
     return true
 end
 
+# From a reaction, create a `MassActionJump` structure (used by JumpPrcoesses) from it.
 @inline function makemajump(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, netstoich = rx
-    zeroorder = (length(substoich) == 0)
+
+    # If nothing changes in the reaction, throws an error. 
+    net_stoch = filter(p -> !drop_dynamics(p[1]), netstoich)
+    isempty(net_stoch) &&
+        error("$rx has no net stoichiometry change once accounting for constant and boundary condition species. This is not supported.")
+
+    # Creates the reactant stoichiometry vector (which mapping each substrate to its stoichiometry).
     reactant_stoch = Vector{Pair{Any, eltype(substoich)}}()
     @inbounds for (i, spec) in enumerate(substrates)
         # move constant species into the rate
@@ -279,23 +342,26 @@ end
         end
     end
 
+    # For non-zeroth order reaction (e.g. `0 --> X`), and if we use a combinatorial rate law,
+    # modifies the rate with its mass action coefficient (no need if coefficient is 1).
+    zeroorder = (length(substoich) == 0)
     if (!zeroorder) && combinatoric_ratelaw
         coef = prod(factorial, substoich)
         (!isone(coef)) && (rate /= coef)
     end
 
-    net_stoch = filter(p -> !drop_dynamics(p[1]), netstoich)
-    isempty(net_stoch) &&
-        error("$rx has no net stoichiometry change once accounting for constant and boundary condition species. This is not supported.")
-
     MassActionJump(Num(rate), reactant_stoch, net_stoch, scale_rates = false,
         useiszero = false)
 end
 
 # recursively visit each neighbor's rooted tree and mark everything in it as vrj
 function dfs_mark!(isvrjvec, visited, depgraph, i)
+    # Marks the current node as visited. Finds all neighbours.
     visited[i] = true
     nhbrs = depgraph[i]
+
+    # For each neighbour to the input node, if we have not visited it, mark it as visited and 
+    # recursively calls `dfs_mark!` on the neighbour.
     for nhbr in nhbrs
         if !visited[nhbr]
             isvrjvec[nhbr] = true
@@ -313,22 +379,32 @@ function get_depgraph(rs)
     eqeq_dependencies(jdeps, vdeps).fadjlist
 end
 
+# Creates the full vector of jumps (which is then used to generate a JumpProcesses JumpSystem).
 function assemble_jumps(rs; combinatoric_ratelaws = true)
+    # Initiates mass action, constant rate, and variable rate jump vectors. These are then filled
+    # and finally (concatenated and) return as the output.
     meqs = MassActionJump[]
     ceqs = ConstantRateJump[]
     veqs = VariableRateJump[]
+    
+    rxvars = []
     unknownset = Set(get_unknowns(rs))
+
+    # Error checks (all unknowns are species and the system has at least one reaction).
     all(isspecies, unknownset) ||
         error("Conversion to JumpSystem currently requires all unknowns to be species.")
-    rxvars = []
-
     isempty(get_rxs(rs)) &&
         error("Must give at least one reaction before constructing a JumpSystem.")
 
     # first we determine vrjs with an explicit time-dependent rate
+    # Next `isvrjvec` lists variable rate jumps, and `havevrjs` is 
     rxs = get_rxs(rs)
     isvrjvec = falses(length(rxs))
     havevrjs = false
+    
+    # For each reaction, check if any of its rate's variables is the time iv.
+    # This, and the next block, computes `isvrjvec` (which jumps are variable rate jumps) and 
+    # `havevrjs` (whether the system has any variable rate jumps).
     for (i, rx) in enumerate(rxs)
         empty!(rxvars)
         (rx.rate isa Symbolic) && get_variables!(rxvars, rx.rate)
@@ -353,20 +429,28 @@ function assemble_jumps(rs; combinatoric_ratelaws = true)
         end
     end
 
+    # Loops through all reactions. Determines its type, creates the corresponding jump, and adds 
+    # it to the correct vector.
     for (i, rx) in enumerate(rxs)
+        # Empties the `rxvars` cache and fills it with the variables in the reaction's rate.
+        # Does it in thsi way to avoid allocation.
         empty!(rxvars)
         (rx.rate isa Symbolic) && get_variables!(rxvars, rx.rate)
 
         isvrj = isvrjvec[i]
-        if (!isvrj) && ismassaction(rx, rs; rxvars, haveivdep = false, unknownset)
+        # If we are a mass action jump, create a mass action jump using `makemajump`.
+        if (!isvrj) && ismassaction(rx, rs; rxvars = rxvars, haveivdep = false,
+                        unknownset = unknownset)
             push!(meqs, makemajump(rx; combinatoric_ratelaw = combinatoric_ratelaws))
         else
+            # Computes the jump rate and jump affect functions (these are identical for vrjs and crjs).
             rl = jumpratelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws)
             affect = Vector{Equation}()
             for (spec, stoich) in rx.netstoich
                 # don't change species that are constant or BCs
                 (!drop_dynamics(spec)) && push!(affect, spec ~ spec + stoich)
             end
+
             if isvrj
                 push!(veqs, VariableRateJump(rl, affect))
             else
@@ -379,38 +463,49 @@ end
 
 ### Equation Coupling ###
 
-# merge constraint components with the ReactionSystem components
+# Adds any normal (differential or algebraic) equations to the equations (`eqs`) generated
+# Also adds any boundary value species to the unknowns vector.
+
+# Also removes any equations produced by constant rate species.
+# Also handles the addition of equations from boundary condition species.
+
 # also handles removing BC and constant species
-function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false)
+# Also handles the elimination of conserved quantities (if this was requested).
+function add_coupled_eqs!(eqs, rs::ReactionSystem, ind_us, ind_sps; remove_conserved = 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
+    rs_us = get_unknowns(rs)
+    us = any(isbc, rs_us) ? vcat(ind_us, filter(isbc, rs_us)) : ind_us
     ps = get_ps(rs)
 
-    # make dependent species observables and add conservation constants as parameters
+    # If conservation law removal is designated, adds conservation law constants as parameters
+    # and creates observables of the elimianted species.
     if remove_conserved
         nps = get_networkproperties(rs)
 
-        # add the conservation constants as parameters and set their values
+        # Adds the conservation constant as system parameters. For each, sets a default value
+        # according to its conservation equation.
         ps = vcat(ps, collect(eq.lhs for eq in nps.constantdefs))
         defs = copy(MT.defaults(rs))
         for eq in nps.constantdefs
             defs[eq.lhs] = eq.rhs
         end
 
-        # add the dependent species as observed
+        # Adds the (by conservation laws) eliminated species as observables.
         obs = copy(MT.observed(rs))
         append!(obs, nps.conservedeqs)
     else
+        # If we have not eliminated any conservation laws, the defaults and observables are unchanged.
         defs = MT.defaults(rs)
         obs = MT.observed(rs)
     end
 
-    ceqs = Equation[eq for eq in get_eqs(rs) if eq isa Equation]
+    # Appends any coupled (differential and algebraic) equations from the systems to the
+    # equations vector. Gives a note if coupled equations are combined with conservation law elimination.
+    ceqs = Equation[eq for eq in nonreactions(rs)]
     if !isempty(ceqs)
         if remove_conserved
             @info """
-                  Be careful mixing constraints and elimination of conservation laws.
+                  Be careful mixing coupled equations and elimination of conservation laws.
                   Catalyst does not check that the conserved equations still hold for the
                   final coupled system of equations. Consider using `remove_conserved =
                   false` and instead calling ModelingToolkit.structural_simplify to simplify
@@ -420,20 +515,20 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved
         append!(eqs, ceqs)
     end
 
-    eqs, sts, ps, obs, defs
+    eqs, us, ps, obs, defs
 end
 
-# used by flattened systems that don't support constraint equations currently
-function error_if_constraints(::Type{T}, sys::ReactionSystem) where {T <: MT.AbstractSystem}
+# used by flattened systems that don't support coupled equations currently
+function error_if_coupled(::Type{T}, sys::ReactionSystem) where {T <: MT.AbstractSystem}
     any(eq -> eq isa Equation, get_eqs(sys)) &&
         error("Can not convert to a system of type ", T,
-            " when there are constraint equations.")
+              " when there are coupled equations.")
     nothing
 end
 
 ### Utility ###
 
-# Throws an error when attempting to convert a spatial system to an unsupported type.
+# Throws an error when attempting to convert a spatial system to an unsuported type.
 function spatial_convert_err(rs::ReactionSystem, systype)
     isspatial(rs) && error("Conversion to $systype is not supported for spatial networks.")
 end
@@ -475,43 +570,44 @@ Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.ODESystem`.
 
 Keyword args and default values:
-- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate law,
-  i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
-  `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
-  ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
-  itself defaults to true).
-- `remove_conserved=false`, if set to `true` will calculate conservation laws of the
-  underlying set of reactions (ignoring constraint equations), and then apply them to reduce
-  the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
+- `combinatoric_ratelaws = true` uses factorial scaling factors in calculating the rate law,
+ i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
+ `combinatoric_ratelaws = false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
+ ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
+ itself defaults to true).
+- `remove_conserved = false`, if set to `true` will calculate conservation laws of the
+ underlying set of reactions (ignoring coupled equations), and then apply them to reduce
+ the number of equations.
 """
 function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs),
-        combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        include_zero_odes = true, remove_conserved = false, remove_conserved_warn = true,
-        checks = false, default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
-        kwargs...)
+                      combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+                      include_zero_odes = true, remove_conserved = false, checks = false,
+                      default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
+                      kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, ODESystem)
-    check_cons_warning(remove_conserved, remove_conserved_warn)
 
-    fullrs = Catalyst.flatten(rs)
-    remove_conserved && conservationlaws(fullrs)
-    ists, ispcs = get_indep_sts(fullrs, remove_conserved)
-    eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
-        include_zero_odes)
-    eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
+    # Flattens the system. Next (if `remove_conserved = true`) uses `conservationlaws` to compute
+    # these (these are then automatically cached in `rs`).
+    flatrs = Catalyst.flatten(rs)
+    remove_conserved && conservationlaws(flatrs)
 
-    ODESystem(eqs, get_iv(fullrs), us, ps;
-        observed = obs,
-        name,
-        defaults = _merge(defaults, defs),
-        checks,
-        continuous_events = MT.get_continuous_events(fullrs),
-        discrete_events = MT.get_discrete_events(fullrs),
-        kwargs...)
+    # Find independent unknowns and species. Assembles the drift ODE, and then the full set of equations.
+    # Here, `add_coupled_eqs!` performs some assembly in addition to adding coupled equations.
+    ind_us, ind_sps = get_indep_us(flatrs, remove_conserved)
+    eqs = assemble_drift(flatrs, ind_sps; combinatoric_ratelaws, remove_conserved,
+                         include_zero_odes)
+    eqs, us, ps, obs, defs = add_coupled_eqs!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
+
+    ODESystem(eqs, get_iv(flatrs), us, ps;
+              observed = obs,
+              name,
+              defaults = _merge(defaults, defs),
+              checks,
+              continuous_events = MT.get_continuous_events(flatrs),
+              discrete_events = MT.get_discrete_events(flatrs),
+              kwargs...)
 end
 
 """
@@ -522,16 +618,14 @@ Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.NonlinearSystem`.
 
 Keyword args and default values:
-- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate law,
-  i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
-  `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
-  ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
-  itself defaults to true).
-- `remove_conserved=false`, if set to `true` will calculate conservation laws of the
-  underlying set of reactions (ignoring constraint equations), and then apply them to reduce
-  the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
+- `combinatoric_ratelaws = true` uses factorial scaling factors in calculating the rate law,
+ i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
+ `combinatoric_ratelaws = false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
+ ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
+ itself defaults to true).
+- `remove_conserved = false`, if set to `true` will calculate conservation laws of the
+ underlying set of reactions (ignoring coupled equations), and then apply them to reduce
+ the number of equations.
 """
 function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
@@ -547,13 +641,17 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
         error("Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depend on $(get_iv(rs))) to a `NonlinearSystem`. This is not possible. if you are intending to compute system steady states, consider creating and solving a `SteadyStateProblem.")
     end
 
-    # Generates system equations.
-    fullrs = Catalyst.flatten(rs)
-    remove_conserved && conservationlaws(fullrs)
-    ists, ispcs = get_indep_sts(fullrs, remove_conserved)
-    eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
-        as_odes = false, include_zero_odes)
-    eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
+    # Flattens the system. Next (if `remove_conserved = true`) uses `conservationlaws` to compute
+    # these (and automatically cache these in `rs`).
+    flatrs = Catalyst.flatten(rs)
+    remove_conserved && conservationlaws(flatrs)
+
+    # Find independent unknowns and species. Assembles the drift ODE, and then the full set of equations.
+    # Here, `add_coupled_eqs!` performs some assembly in addition to adding coupled equations.
+    ind_us, ind_sps = get_indep_us(flatrs, remove_conserved)
+    eqs = assemble_drift(flatrs, ind_sps; combinatoric_ratelaws, remove_conserved,
+                         as_odes = false, include_zero_odes)
+    eqs, us, ps, obs, defs = add_coupled_eqs!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
 
     # Throws a warning if there are differential equations in non-standard format.
     # Next, sets all differential terms to `0`.
@@ -561,11 +659,11 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     eqs = [remove_diffs(eq.lhs) ~ remove_diffs(eq.rhs) for eq in eqs]
 
     NonlinearSystem(eqs, us, ps;
-        name,
-        observed = obs,
-        defaults = _merge(defaults, defs),
-        checks,
-        kwargs...)
+                    name,
+                    observed = obs,
+                    defaults = _merge(defaults, defs),
+                    checks,
+                    kwargs...)
 end
 
 # Ideally, when `ReactionSystem`s are converted to `NonlinearSystem`s, any coupled ODEs should be 
@@ -575,23 +673,23 @@ end
 function nonlinear_convert_differentials_check(rs::ReactionSystem)
     for eq in filter(is_diff_equation, equations(rs))
         # For each differential equation, checks in order: 
-        # If there is a differential on the right hand side.
+        # If there is a differential on the right-hand side.
         # If the lhs is not on the form D(...).
-        # If the lhs upper level function is not a differential w.r.t. time.
+        # If the lhs upper-level function is not a differential w.r.t. time.
         # If the content of the differential is not a variable (and nothing more).
-        # If either of this is a case, throws the warning.
-        if hasnode(Symbolics.is_derivative, eq.rhs) ||
-           !Symbolics.is_derivative(eq.lhs) ||
-           !isequal(Symbolics.operation(eq.lhs), Differential(get_iv(rs))) ||
-           (length(arguments(eq.lhs)) != 1) ||
-           !any(isequal(arguments(eq.lhs)[1]), nonspecies(rs))
+        # If either of this is a case, throw the warning.
+        if Symbolics._occursin(Symbolics.is_derivative, eq.rhs) ||
+                               !Symbolics.istree(eq.lhs) ||
+                               !isequal(Symbolics.operation(eq.lhs), Differential(get_iv(rs))) || 
+                               (length(arguments(eq.lhs)) != 1) ||
+                               !any(isequal(arguments(eq.lhs)[1]), nonspecies(rs))
             error("You are attempting to convert a `ReactionSystem` coupled with differential equations to a `NonlinearSystem`. However, some of these differentials are not of the form `D(x) ~ ...` where: 
-                    (1) The left-hand side is a differential of a single variable with respect to the time independent variable, and
-                    (2) The right-hand side does not contain any differentials.
-                This is generally not permitted.
-
-                If you still would like to perform this conversion, please use the `all_differentials_permitted = true` option. In this case, all differentials will be set to `0`. 
-                However, it is recommended to proceed with caution to ensure that the produced nonlinear equation makes sense for your intended application."
+                  (1) The left-hand side is a differential of a single variable with respect to the time-independent variable, and
+                  (2) The right-hand side does not contain any differentials.
+                  This is generally not permitted.
+                                  
+                  If you still would like to perform this conversion, please use the `all_differentials_permitted = true` option. In this case, all differential will be set to `0`. 
+                  However, it is recommended to proceed with caution to ensure that the produced nonlinear equation makes sense for your intended application."
             )
         end
     end
@@ -605,41 +703,40 @@ Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.SDESystem`.
 
 Notes:
-- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate law,
-  i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
-  `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
-  ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
-  itself defaults to true).
-- `remove_conserved=false`, if set to `true` will calculate conservation laws of the
-  underlying set of reactions (ignoring constraint equations), and then apply them to reduce
-  the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
+- `combinatoric_ratelaws = true` uses factorial scaling factors in calculating the rate law,
+ i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
+ `combinatoric_ratelaws = false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
+ ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
+ itself defaults to true).
+- `remove_conserved = false`, if set to `true` will calculate conservation laws of the
+ underlying set of reactions (ignoring coupled equations), and then apply them to reduce
+ the number of equations.
+- Does not currently support `ReactionSystem`s that includes coupled algebraic or
+ differential equations.
 """
 function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
-        name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        include_zero_odes = true, checks = false, remove_conserved = false,
-        remove_conserved_warn = true, default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
-        kwargs...)
+                      name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+                      include_zero_odes = true, checks = false, remove_conserved = false,
+                      default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
+                      kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, SDESystem)
     check_cons_warning(remove_conserved, remove_conserved_warn)
 
+    # Flattens the system. Next (if `remove_conserved = true`) uses `conservationlaws` to compute
+    # these (and automatically cache these in `rs`).
     flatrs = Catalyst.flatten(rs)
-
     remove_conserved && conservationlaws(flatrs)
-    ists, ispcs = get_indep_sts(flatrs, remove_conserved)
-    eqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws, include_zero_odes,
-        remove_conserved)
-    noiseeqs = assemble_diffusion(flatrs, ists, ispcs;
-        combinatoric_ratelaws, remove_conserved)
-    eqs, us, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved)
-
-    if any(isbc, get_unknowns(flatrs))
-        @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead."
-    end
+
+    # Find independent unknowns and species. Assembles the drift ODE, and then the full set of
+    # equations. Also assemblies the drift (stochasticity) terms.
+    # Here, `add_coupled_eqs!` performs some assembly in addition to adding coupled equations.
+    ind_us, ind_sps = get_indep_us(flatrs, remove_conserved)
+    eqs = assemble_drift(flatrs, ind_sps; combinatoric_ratelaws, include_zero_odes,
+                         remove_conserved)
+    noiseeqs = assemble_diffusion(flatrs, ind_us, ind_sps; combinatoric_ratelaws, remove_conserved)
+    eqs, us, ps, obs, defs = add_coupled_eqs!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
 
     SDESystem(eqs, noiseeqs, get_iv(flatrs), us, ps;
         observed = obs,
@@ -653,53 +750,52 @@ end
 
 """
 ```julia
-Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
+Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; combinatoric_ratelaws = true)
 ```
 
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.JumpSystem`.
 
 Notes:
-- `combinatoric_ratelaws=true` uses binomials in calculating the rate law, i.e. for `2S ->
-  0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaws=false` then
-  the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling factor.
-  Defaults to the value given when the `ReactionSystem` was constructed (which itself
-  defaults to true).
-- Does not currently support `ReactionSystem`s that include coupled algebraic or
-  differential equations.
+- `combinatoric_ratelaws = true` uses binomials in calculating the rate law, i.e. for `2S ->
+ 0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaws = false` then
+ the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling factor.
+ Defaults to the value given when the `ReactionSystem` was constructed (which itself
+ defaults to true).
+- Does not currently support `ReactionSystem`s that includes coupled algebraic or
+ differential equations.
 - Does not currently support continuous events as these are not supported by
-  `ModelingToolkit.JumpSystems`.
+ `ModelingToolkit.JumpSystems`.
 """
 function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs),
-        combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        remove_conserved = nothing, checks = false,
-        default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
-        kwargs...)
+                      combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+                      remove_conserved = nothing, checks = false,
+                      default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
+                      kwargs...)
+    # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, JumpSystem)
     (remove_conserved !== nothing) &&
         throw(ArgumentError("Catalyst does not support removing conserved species when converting to JumpSystems."))
 
+    # Flattens the system. Performs post-flattening errors and warning checks.
     flatrs = Catalyst.flatten(rs)
-    error_if_constraints(JumpSystem, flatrs)
-
+    error_if_coupled(JumpSystem, flatrs)
     (length(MT.continuous_events(flatrs)) > 0) &&
         (@warn "continuous_events will be dropped as they are not currently supported by JumpSystems.")
 
-    eqs = assemble_jumps(flatrs; combinatoric_ratelaws)
-
-    # handle BC species
-    sts, ispcs = get_indep_sts(flatrs)
-    any(isbc, get_unknowns(flatrs)) && (sts = vcat(sts, filter(isbc, get_unknowns(flatrs))))
+    # Computes species, parameters, and equations. Adds boundary condition species to species vector.
+    _, ind_sps = get_indep_us(flatrs, remove_conserved)
+    any(isbc, get_unknowns(flatrs)) && (ind_sps = vcat(ind_sps, filter(isbc, get_unknowns(flatrs))))
     ps = get_ps(flatrs)
-
-    JumpSystem(eqs, get_iv(flatrs), sts, ps;
-        observed = MT.observed(flatrs),
-        name,
-        defaults = _merge(defaults, MT.defaults(flatrs)),
-        checks,
-        discrete_events = MT.discrete_events(flatrs),
-        kwargs...)
+    eqs = assemble_jumps(flatrs; combinatoric_ratelaws)
+    
+    JumpSystem(eqs, get_iv(flatrs), ind_sps, ps;
+               observed = MT.observed(flatrs),
+               name,
+               defaults = _merge(defaults, MT.defaults(flatrs)),
+               checks,
+               discrete_events = MT.discrete_events(flatrs),
+               kwargs...)
 end
 
 ### Problems ###
@@ -717,10 +813,10 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
         remove_conserved, remove_conserved_warn)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
-    if structural_simplify
-        (osys = MT.structural_simplify(osys))
+    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.")
+        error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify = true` within `ODEProblem` call.")
     else
         osys = complete(osys)
     end
@@ -746,25 +842,25 @@ end
 
 # SDEProblem from AbstractReactionNetwork
 function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
-        p = DiffEqBase.NullParameters(), args...;
-        name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        include_zero_odes = true, checks = false, check_length = false, remove_conserved = false,
-        remove_conserved_warn = true, structural_simplify = false, kwargs...)
+                               p = DiffEqBase.NullParameters(), args...;
+                               name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+                               include_zero_odes = true, checks = false, check_length = false,
+                               remove_conserved = false, structural_simplify = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws,
-        include_zero_odes, checks, remove_conserved, remove_conserved_warn)
+                      include_zero_odes, checks, remove_conserved)
+    p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs))
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
-    if structural_simplify
-        (sde_sys = MT.structural_simplify(sde_sys))
+    if structural_simplify 
+        sde_sys = MT.structural_simplify(sde_sys)
     elseif has_alg_equations(rs)
-        error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.")
+        error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify = true` within `ODEProblem` call.")
     else
         sde_sys = complete(sde_sys)
     end
-
-    p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs))
+    
     return SDEProblem(sde_sys, u0map, tspan, pmap, args...; check_length,
         noise_rate_prototype = p_matrix, kwargs...)
 end
@@ -805,10 +901,10 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0,
         remove_conserved, remove_conserved_warn)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
-    if structural_simplify
-        (osys = MT.structural_simplify(osys))
+    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.")
+        error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify = true` within `ODEProblem` call.")
     else
         osys = complete(osys)
     end
@@ -820,10 +916,15 @@ end
 
 # convert symbol of the form :sys.a.b.c to a symbolic a.b.c
 function _symbol_to_var(sys, sym)
+    # The input is a normal Symbol, which exists in the input system (subsystems not considered).
     if hasproperty(sys, sym)
         var = getproperty(sys, sym, namespace = false)
     else
-        strs = split(String(sym), ModelingToolkit.NAMESPACE_SEPARATOR)   # need to check if this should be split of not!!!
+        # If sys does not contain sym it is possible that sys.subsystem contains it (or a 
+        # sub-sub system etc.). Here, "₊" splits the system and subsystem names in the symbol.
+        # We thus split the symbol to find, and go down the subsystem tree to attempt to find
+        # the symbol. 
+        strs = split(String(sym), "₊")   # need to check if this should be split of not!!!
         if length(strs) > 1
             var = getproperty(sys, Symbol(strs[1]), namespace = false)
             for str in view(strs, 2:length(strs))
@@ -837,7 +938,7 @@ function _symbol_to_var(sys, sym)
 end
 
 """
-    symmap_to_varmap(sys, symmap)
+ symmap_to_varmap(sys, symmap)
 
 Given a system and map of `Symbol`s to values, generates a map from
 corresponding symbolic variables/parameters to the values that can be used to
@@ -846,11 +947,11 @@ pass initial conditions and parameter mappings.
 For example,
 ```julia
 sir = @reaction_network sir begin
-    β, S + I --> 2I
-    ν, I --> R
+ β, S + I --> 2I
+ ν, I --> R
 end
 subsys = @reaction_network subsys begin
-    k, A --> B
+ k, A --> B
 end
 @named sys = compose(sir, [subsys])
 ```
@@ -858,15 +959,15 @@ gives
 ```
 Model sys with 3 equations
 Unknowns (5):
-  S(t)
-  I(t)
-  R(t)
-  subsys₊A(t)
-  subsys₊B(t)
+ S(t)
+ I(t)
+ R(t)
+ subsys₊A(t)
+ subsys₊B(t)
 Parameters (3):
-  β
-  ν
-  subsys₊k
+ β
+ ν
+ subsys₊k
 ```
 to specify initial condition and parameter mappings from *symbols* we can use
 ```julia
@@ -878,12 +979,14 @@ pmap   = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0])
 
 Notes:
 - Any `Symbol`, `sym`, within `symmap` must be a valid field of `sys`. i.e.
-  `sys.sym` must be defined.
+ `sys.sym` must be defined.
 """
 function symmap_to_varmap(sys, symmap::Tuple)
+    # If all inputs are of the form `:sym => value`, apply `_symbol_to_var` function. If not,
+    # pass input through unchanged.
     if all(p -> p isa Pair{Symbol}, symmap)
         return ((_symbol_to_var(sys, sym) => val for (sym, val) in symmap)...,)
-    else  # if not all entries map a symbol to value pass through
+    else
         return symmap
     end
 end
@@ -903,17 +1006,17 @@ symmap_to_varmap(sys, symmap) = symmap
 ### Other Conversion-related Functions ###
 
 # the following function is adapted from SymbolicUtils.jl v.19
-# later on (Spetember 2023) modified by Torkel and Shashi (now assumes input not on polynomial form, which is handled elsewhere, previous version failed in these cases anyway).
+# later on (Spetember 2023) modified by Torkel and Shashi (now assumes input not on polynomial form, which is handled elsewhere, the previous version failed in these cases anyway).
 # Copyright (c) 2020: Shashi Gowda, Yingbo Ma, Mason Protter, Julia Computing.
 # MIT license
 """
-    to_multivariate_poly(polyeqs::AbstractVector{BasicSymbolic{Real}})
+ to_multivariate_poly(polyeqs::AbstractVector{BasicSymbolic{Real}})
 
 Convert the given system of polynomial equations to multivariate polynomial representation.
 For example, this can be used in HomotopyContinuation.jl functions.
 """
 function to_multivariate_poly(polyeqs::AbstractVector{Symbolics.BasicSymbolic{Real}})
-    @assert length(polyeqs)>=1 "At least one expression must be passed to `multivariate_poly`."
+    @assert length(polyeqs) >= 1 "At least one expression must be passed to `multivariate_poly`."
 
     pvar2sym, sym2term = SymbolicUtils.get_pvar2sym(), SymbolicUtils.get_sym2term()
     ps = map(polyeqs) do x

From 92e60f313feec2b3b0eea6ed3f66fbc95ac73d78 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 27 Jul 2024 15:59:00 -0400
Subject: [PATCH 2/4] go through again

---
 src/reactionsystem_conversions.jl | 135 ++++++++++++++----------------
 1 file changed, 64 insertions(+), 71 deletions(-)

diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index ee4915da2d..65271a17f3 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -1,7 +1,7 @@
 ### ODE & SDE Assembly ###
 
 """
- oderatelaw(rx; combinatoric_ratelaw = true)
+    oderatelaw(rx; combinatoric_ratelaw = true)
 
 Given a [`Reaction`](@ref), return the symbolic reaction rate law used in
 generated ODEs for the reaction. Note, for a reaction defined by
@@ -26,17 +26,17 @@ function oderatelaw(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, only_use_rate = rx
     rl = rate
 
-    # if the stoichiometric coefficients are not integers, error if asking to scale rates
+    # if the stoichiometric coefficients are not integers error if asking to scale rates
     !all(s -> s isa Union{Integer, Symbolic}, substoich) &&
         (combinatoric_ratelaw == true) &&
         error("Non-integer stoichiometric coefficients require the combinatoric_ratelaw = false keyword to oderatelaw, or passing combinatoric_ratelaws = false to convert or ODEProblem.")
 
-    # If we use mass action, update the rate law according to the reaction substrates.
+    # if we use mass action, update the rate law according to the reaction substrates
     if !only_use_rate
-        # `coef` stores the number bit (i.e. division by factorial of stoichiometries).
-        # Symbolic bit (i.e. multiplication by substrate concentrations) are appended directly 
-        # to the rate law (`rl`). Cases with stoichiometries 1 are handled more efficiently
-        # using separate routines. 
+        # `coef` stores the numeric part (i.e. division by factorial of stoichiometries)
+        # the symbolic part (i.e. multiplication by substrate concentrations) is appended directly 
+        # to the rate law (`rl`)
+        # cases with stoichiometries 1 are handled more efficiently using separate routines
         coef = eltype(substoich) <: Number ? one(eltype(substoich)) : 1
         for (i, stoich) in enumerate(substoich)
             combinatoric_ratelaw && (coef *= factorial(stoich))
@@ -47,19 +47,19 @@ function oderatelaw(rx; combinatoric_ratelaw = true)
     rl
 end
 
-# Function returning `true` for species which shouldn't change from the reactions, 
-# including non-species variables.
+# function returning `true` for species which shouldn't change from the reactions, 
+# including non-species variables
 drop_dynamics(sp) = isconstant(sp) || isbc(sp) || (!isspecies(sp))
 
-# For a `ReactionSystem`, create a vector with all the right-hand side expressions in the
-# corresponding ODE system.
+# for a `ReactionSystem`, create a vector with all the right-hand side expressions in the
+# corresponding ODE system (used by the `assemble_drift` function)
 function assemble_oderhs(rs, isps; combinatoric_ratelaws = true, remove_conserved = false)
-    # Initiates the right-hand side equations as expressions with `0` only.
+    # initiates the right-hand side equations as expressions with `0` only.
     species_to_idx = Dict(x => i for (i, x) in enumerate(isps))
     rhsvec = Any[0 for _ in isps]
 
-    # If we have eliminated conservation laws, create a dictionary taking each species eliminated
-    # in this process to the expression it is replaced by (if not, creates an empty dict).
+    # if we have eliminated conservation laws, create a dictionary taking each species eliminated
+    # in this process to the expression it is replaced by (if not, creates an empty dict)
     nps = get_networkproperties(rs)
     depspec_submap = if remove_conserved
         Dict(eq.lhs => eq.rhs for eq in nps.conservedeqs)
@@ -67,14 +67,14 @@ function assemble_oderhs(rs, isps; combinatoric_ratelaws = true, remove_conserve
         Dict()
     end
 
-    # Loops through all reactions, for each, updates relevant right-hand side equations.
+    # loops through all reactions, for each, updates all relevant right-hand side equations
     for rx in get_rxs(rs)
-        # Computes the reaction's rate late. If conserved quantities have been removed, replaces
-        # these in the rate law.
+        # computes the reaction's rate late
+        # if conserved quantities have been removed, replaces these in the rate law
         rl = oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws)
         remove_conserved && (rl = substitute(rl, depspec_submap))
 
-        # For each species which participates in the reaction, update its equation.
+        # for each species which participates in the reaction, update its equation
         for (spec, stoich) in rx.netstoich
             # dependent species don't get an ODE, so are skipped
             remove_conserved && (spec in nps.depspecs) && continue
@@ -83,8 +83,8 @@ function assemble_oderhs(rs, isps; combinatoric_ratelaws = true, remove_conserve
             drop_dynamics(spec) && continue
 
             i = species_to_idx[spec]
-            # If the relevant RHS equation is 0 (probably because it has not been
-            # modified yet), we use a different routine (compared to if it is non-zero).
+            # if the relevant RHS equation is 0 (probably because it has not been
+            # modified yet), we use a different routine (compared to if it is non-zero)
             if _iszero(rhsvec[i])
                 if stoich isa Symbolic
                     rhsvec[i] = stoich * rl
@@ -112,7 +112,7 @@ function assemble_oderhs(rs, isps; combinatoric_ratelaws = true, remove_conserve
 end
 
 # For a `ReactionSystem`, creates a vector with all the differential equations corresponding to
-# the ODEs generated by the reaction rate equation. The `as_odes` determine whether the LHS should
+# the ODEs generated by the reaction rate equation. The `as_odes` argument determine whether the LHS should
 # be differential (when generating e.g. a `ODESystem`) or `0` (when generating e.g. a NonlinearSystem).
 function assemble_drift(rs, ind_sps; combinatoric_ratelaws = true, as_odes = true,
                         include_zero_odes = true, remove_conserved = false)
@@ -128,7 +128,7 @@ function assemble_drift(rs, ind_sps; combinatoric_ratelaws = true, as_odes = tru
 end
 
 # this doesn't work with coupled equations currently
-# Assemblies the diffusion equations terms used in the chemical Langevin equation SDE.
+# assemblies the diffusion equations terms used in the chemical Langevin equation SDE
 function assemble_diffusion(rs, us, ind_sps; combinatoric_ratelaws = true,
                             remove_conserved = false)
     # as BC species should ultimately get an equation, we include them in the noise matrix
@@ -139,8 +139,8 @@ function assemble_diffusion(rs, us, ind_sps; combinatoric_ratelaws = true,
     eqs .= 0
     species_to_idx = Dict((x => i for (i, x) in enumerate(ispcs)))
 
-    # If we have eliminated conservation laws, create a dictionary taking each species eliminated
-    # in this process to the expression it is replaced by (if not, creates an empty dict).
+    # if we have eliminated conservation laws, create a dictionary taking each species eliminated
+    # in this process to the expression it is replaced by (if not, creates an empty dict)
     nps = get_networkproperties(rs)
     depspec_submap = if remove_conserved
         Dict(eq.lhs => eq.rhs for eq in nps.conservedeqs)
@@ -148,15 +148,15 @@ function assemble_diffusion(rs, us, ind_sps; combinatoric_ratelaws = true,
         Dict()
     end
 
-    # Loops through all reactions, for each, updates relevant noise terms in the noise matrix.
+    # loops through all reactions, for each, updates relevant noise terms in the noise matrix
     for (j, rx) in enumerate(get_rxs(rs))
-        # Computes the reactions (noise) rate late. This is the sqrt of the rate late, multiplied
+        # Computes the reaction's (noise) rate law. This is the sqrt of the rate law, multiplied
         # by any noise scaling term. If conserved quantities have been removed, these are substituted in.
         rlsqrt = sqrt(abs(oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws)))
         hasnoisescaling(rx) && (rlsqrt *= getnoisescaling(rx))
         remove_conserved && (rlsqrt = substitute(rlsqrt, depspec_submap))
 
-        # For each species which participate in the reaction, updates its corresponding noise terms.
+        # for each species which participate in the reaction, updates its corresponding noise terms
         for (spec, stoich) in rx.netstoich
             # dependent species don't get an equation
             remove_conserved && (spec in nps.depspecs) && continue
@@ -165,8 +165,8 @@ function assemble_diffusion(rs, us, ind_sps; combinatoric_ratelaws = true,
             drop_dynamics(spec) && continue
 
             i = species_to_idx[spec]
-            # Differentiates if the stoichiometry is a symbolic (e.g. a parameter). If not it might
-            # be +/- 1, in which case we initiate the term differently.
+            # Different routine if the stoichiometry is a symbolic (e.g. a parameter)
+            # if not it might e +/- 1, in which case we initiate the term differently
             if stoich isa Symbolic
                 eqs[i, j] = stoich * rlsqrt
             else
@@ -208,18 +208,18 @@ function jumpratelaw(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, only_use_rate = rx
     rl = rate
 
-    # If the `only_use_rate` is unused, update the rate according to the stoichiometries.
+    # if the `only_use_rate` argument is unused, update the rate according to the stoichiometries
     if !only_use_rate
         # Initiates the coefficient as "1", using the same type as in `substoich`.
         coef = eltype(substoich) <: Number ? one(eltype(substoich)) : 1
 
-        # For each substrate, update the rate law (and coefficient) terms according to the
-        # substrate's stoichiometry.
+        # for each substrate, update the rate law (and coefficient) terms according to the
+        # substrate's stoichiometry
         for (i, stoich) in enumerate(substoich)
             spec = substrates[i]
 
             # For symbolic substrates (e.g. parameters), simply update the rate law with the
-            # appropriate expression. If not, we  can simplify the expression (and also treat
+            # appropriate expression. If not, we can simplify the expression (and also treat
             # the case with stoichiometry = 1 separately).
             if stoich isa Symbolic
                 rl *= combinatoric_ratelaw ? binomial(spec, stoich) :
@@ -237,8 +237,8 @@ function jumpratelaw(rx; combinatoric_ratelaw = true)
             end
         end
 
-        # If we use a combinatoric rate law, we have generated a numeric coefficient which we
-        # will (if it is not 1) divde our rate law with (to generate the full rate law).
+        # if we use a combinatoric rate law, we have generated a numeric coefficient which we
+        # will (if it is not 1) divide our rate law with (to generate the full rate law)
         combinatoric_ratelaw && !isequal(coef, one(coef)) && (rl /= coef)
     end
     rl
@@ -307,7 +307,8 @@ function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
         haveivdep && return false
     end
 
-    # The `only_use_rate` option disables mass action for the reaction (=> not mass action reaction).
+    # the `only_use_rate` option disables mass action for the reaction (e.g. if the  => arrow was
+    # used in the DSL)
     rx.only_use_rate && return false
 
     # not mass action if have a non-constant, variable species in the rate expression
@@ -318,16 +319,16 @@ function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
     return true
 end
 
-# From a reaction, create a `MassActionJump` structure (used by JumpPrcoesses) from it.
+# from a reaction, create a `MassActionJump` structure (as used by JumpPrcoesses)
 @inline function makemajump(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, netstoich = rx
 
-    # If nothing changes in the reaction, throws an error. 
+    # if nothing changes in the reaction (e.g. `X --> X`), throws an error
     net_stoch = filter(p -> !drop_dynamics(p[1]), netstoich)
     isempty(net_stoch) &&
         error("$rx has no net stoichiometry change once accounting for constant and boundary condition species. This is not supported.")
 
-    # Creates the reactant stoichiometry vector (which mapping each substrate to its stoichiometry).
+    # creates the reactant stoichiometry vector (which maps each substrate to its stoichiometry)
     reactant_stoch = Vector{Pair{Any, eltype(substoich)}}()
     @inbounds for (i, spec) in enumerate(substrates)
         # move constant species into the rate
@@ -342,8 +343,8 @@ end
         end
     end
 
-    # For non-zeroth order reaction (e.g. `0 --> X`), and if we use a combinatorial rate law,
-    # modifies the rate with its mass action coefficient (no need if coefficient is 1).
+    # for non-zeroth order reaction (e.g. `X1 --> X2`), and if we use a combinatorial rate law,
+    # modifies the rate with its mass action coefficient (no need if coefficient is 1)
     zeroorder = (length(substoich) == 0)
     if (!zeroorder) && combinatoric_ratelaw
         coef = prod(factorial, substoich)
@@ -360,7 +361,7 @@ function dfs_mark!(isvrjvec, visited, depgraph, i)
     visited[i] = true
     nhbrs = depgraph[i]
 
-    # For each neighbour to the input node, if we have not visited it, mark it as visited and 
+    # for each neighbour to the input node, if we have not visited it, mark it as visited and 
     # recursively calls `dfs_mark!` on the neighbour.
     for nhbr in nhbrs
         if !visited[nhbr]
@@ -372,39 +373,34 @@ function dfs_mark!(isvrjvec, visited, depgraph, i)
 end
 
 # get_depgraph(rs)[i] is the list of reactions with rates depending on species changed by
-# i'th reaction.
+# i'th reaction
 function get_depgraph(rs)
     jdeps = asgraph(rs)
     vdeps = variable_dependencies(rs)
     eqeq_dependencies(jdeps, vdeps).fadjlist
 end
 
-# Creates the full vector of jumps (which is then used to generate a JumpProcesses JumpSystem).
+# creates the full vector of jumps (which is then used to generate a JumpProcesses JumpSystem)
 function assemble_jumps(rs; combinatoric_ratelaws = true)
-    # Initiates mass action, constant rate, and variable rate jump vectors. These are then filled
-    # and finally (concatenated and) return as the output.
     meqs = MassActionJump[]
     ceqs = ConstantRateJump[]
     veqs = VariableRateJump[]
-    
     rxvars = []
     unknownset = Set(get_unknowns(rs))
 
-    # Error checks (all unknowns are species and the system has at least one reaction).
     all(isspecies, unknownset) ||
         error("Conversion to JumpSystem currently requires all unknowns to be species.")
     isempty(get_rxs(rs)) &&
         error("Must give at least one reaction before constructing a JumpSystem.")
 
     # first we determine vrjs with an explicit time-dependent rate
-    # Next `isvrjvec` lists variable rate jumps, and `havevrjs` is 
     rxs = get_rxs(rs)
     isvrjvec = falses(length(rxs))
     havevrjs = false
     
-    # For each reaction, check if any of its rate's variables is the time iv.
-    # This, and the next block, computes `isvrjvec` (which jumps are variable rate jumps) and 
-    # `havevrjs` (whether the system has any variable rate jumps).
+    # for each reaction, check if any of its rate's variables is the time iv
+    # this, and the next block, computes `isvrjvec` (which jumps are variable rate jumps) and 
+    # `havevrjs` (whether the system has any variable rate jumps)
     for (i, rx) in enumerate(rxs)
         empty!(rxvars)
         (rx.rate isa Symbolic) && get_variables!(rxvars, rx.rate)
@@ -429,21 +425,21 @@ function assemble_jumps(rs; combinatoric_ratelaws = true)
         end
     end
 
-    # Loops through all reactions. Determines its type, creates the corresponding jump, and adds 
-    # it to the correct vector.
+    # loops through all reactions.
+    # determines its type, creates the corresponding jump, and adds it to the correct vector.
     for (i, rx) in enumerate(rxs)
-        # Empties the `rxvars` cache and fills it with the variables in the reaction's rate.
-        # Does it in thsi way to avoid allocation.
+        # empties the `rxvars` cache and fills it with the variables in the reaction's rate
+        # does it in this way to avoid allocation (to make things run faster)
         empty!(rxvars)
         (rx.rate isa Symbolic) && get_variables!(rxvars, rx.rate)
 
         isvrj = isvrjvec[i]
-        # If we are a mass action jump, create a mass action jump using `makemajump`.
+        # if we are a mass action jump, create a mass action jump using `makemajump`
         if (!isvrj) && ismassaction(rx, rs; rxvars = rxvars, haveivdep = false,
                         unknownset = unknownset)
             push!(meqs, makemajump(rx; combinatoric_ratelaw = combinatoric_ratelaws))
         else
-            # Computes the jump rate and jump affect functions (these are identical for vrjs and crjs).
+            # computes the jump rate and jump affect functions (these are identical for vrjs and crjs)
             rl = jumpratelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws)
             affect = Vector{Equation}()
             for (spec, stoich) in rx.netstoich
@@ -463,38 +459,35 @@ end
 
 ### Equation Coupling ###
 
-# Adds any normal (differential or algebraic) equations to the equations (`eqs`) generated
-# Also adds any boundary value species to the unknowns vector.
-
-# Also removes any equations produced by constant rate species.
-# Also handles the addition of equations from boundary condition species.
-
+# adds any normal (differential or algebraic) equations to the equations (`eqs`) generated
+# also adds any boundary value species to the unknowns vector
+# also removes any equations produced by constant rate species
+# also handles the addition of equations from boundary condition species
 # also handles removing BC and constant species
-# Also handles the elimination of conserved quantities (if this was requested).
+# also handles the elimination of conserved quantities (if this was requested)
 function add_coupled_eqs!(eqs, rs::ReactionSystem, ind_us, ind_sps; remove_conserved = false)
     # if there are BC species, put them after the independent species
     rs_us = get_unknowns(rs)
     us = any(isbc, rs_us) ? vcat(ind_us, filter(isbc, rs_us)) : ind_us
     ps = get_ps(rs)
 
-    # If conservation law removal is designated, adds conservation law constants as parameters
+    # if conservation law removal is designated, adds conservation law constants as parameters
     # and creates observables of the elimianted species.
     if remove_conserved
         nps = get_networkproperties(rs)
-
-        # Adds the conservation constant as system parameters. For each, sets a default value
-        # according to its conservation equation.
+        
+        # make dependent species observables and add conservation constants as parameters
         ps = vcat(ps, collect(eq.lhs for eq in nps.constantdefs))
         defs = copy(MT.defaults(rs))
         for eq in nps.constantdefs
             defs[eq.lhs] = eq.rhs
         end
 
-        # Adds the (by conservation laws) eliminated species as observables.
+        # adds the (by conservation laws) eliminated species as observables
         obs = copy(MT.observed(rs))
         append!(obs, nps.conservedeqs)
     else
-        # If we have not eliminated any conservation laws, the defaults and observables are unchanged.
+        # if we have not eliminated any conservation laws, the defaults and observables are unchanged
         defs = MT.defaults(rs)
         obs = MT.observed(rs)
     end

From 32112822f555f39c5e9207f0220e2a31e5c0fcc5 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 27 Jul 2024 16:29:26 -0400
Subject: [PATCH 3/4] updates

---
 src/reactionsystem_conversions.jl | 307 +++++++++++++++---------------
 1 file changed, 154 insertions(+), 153 deletions(-)

diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index 65271a17f3..e9787646f4 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -18,9 +18,9 @@ the expression that is returned will be `k * (X(t)^2/2) * (Y(t)^3/6)`.
 Notes:
 - Allocates
 - `combinatoric_ratelaw = true` uses factorial scaling factors in calculating the
- rate law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
- `combinatoric_ratelaw = false` then the ratelaw is `k*S^2`, i.e. the scaling
- factor is ignored.
+    rate law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
+    `combinatoric_ratelaw = false` then the ratelaw is `k*S^2`, i.e. the scaling
+    factor is ignored.
 """
 function oderatelaw(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, only_use_rate = rx
@@ -47,9 +47,9 @@ function oderatelaw(rx; combinatoric_ratelaw = true)
     rl
 end
 
-# function returning `true` for species which shouldn't change from the reactions, 
+# Function returning `true` for species which shouldn't change from the reactions, 
 # including non-species variables
-drop_dynamics(sp) = isconstant(sp) || isbc(sp) || (!isspecies(sp))
+drop_dynamics(spec) = isconstant(spec) || isbc(spec) || (!isspecies(spec))
 
 # for a `ReactionSystem`, create a vector with all the right-hand side expressions in the
 # corresponding ODE system (used by the `assemble_drift` function)
@@ -165,7 +165,7 @@ function assemble_diffusion(rs, us, ind_sps; combinatoric_ratelaws = true,
             drop_dynamics(spec) && continue
 
             i = species_to_idx[spec]
-            # Different routine if the stoichiometry is a symbolic (e.g. a parameter)
+            # different routine if the stoichiometry is a symbolic (e.g. a parameter)
             # if not it might e +/- 1, in which case we initiate the term differently
             if stoich isa Symbolic
                 eqs[i, j] = stoich * rlsqrt
@@ -200,9 +200,9 @@ binomial(Y,3)`.
 Notes:
 - Allocates
 - `combinatoric_ratelaw = true` uses binomials in calculating the rate law, i.e. for `2S ->
- 0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaw = false` then
- the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling
- factor.
+  0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaw = false` then
+  the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling
+  factor.
 """
 function jumpratelaw(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, only_use_rate = rx
@@ -261,21 +261,21 @@ depend explicitly on the independent variable (usually time).
 - `rx`, the [`Reaction`](@ref).
 - `rs`, a [`ReactionSystem`](@ref) containing the reaction.
 - Optional: `rxvars`, `Variable`s which are not in `rxvars` are ignored as
- possible dependencies.
+  possible dependencies.
 - Optional: `haveivdep`, `true` if the [`Reaction`](@ref) `rate` field
- explicitly depends on any independent variable (i.e. t or for spatial systems x,y,etc).
- If not set, will be automatically calculated.
+  explicitly depends on any independent variable (i.e. t or for spatial systems x,y,etc).
+  If not set, will be automatically calculated.
 - Optional: `unknownset`, set of unknowns which if the rxvars are within mean rx is
- non-mass action.
+  non-mass action.
 - Optional: `ivset`, a `Set` of the independent variables of the system. If not provided and
- the system is spatial, i.e. `isspatial(rs) == true`, it will be created with all the
- spatial variables and the time variable. If the rate expression contains any element of
- `ivset`, then `ismassaction(rx,rs) == false`. Pass a custom set to control this behavior.
+  the system is spatial, i.e. `isspatial(rs) == true`, it will be created with all the
+  spatial variables and the time variable. If the rate expression contains any element of
+  `ivset`, then `ismassaction(rx,rs) == false`. Pass a custom set to control this behavior.
 
 Notes:
 - Non-integer stoichiometry is treated as non-mass action. This includes
- symbolic variables/terms or floating point numbers for stoichiometric
- coefficients.
+  symbolic variables/terms or floating point numbers for stoichiometric
+  coefficients.
 """
 function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
     haveivdep::Union{Nothing, Bool} = nothing,
@@ -302,7 +302,7 @@ function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
         else
             ivdep = any(var -> isequal(get_iv(rs), var), rxvars)
         end
-            ivdep && return false
+        ivdep && return false
     else
         haveivdep && return false
     end
@@ -322,6 +322,7 @@ end
 # from a reaction, create a `MassActionJump` structure (as used by JumpPrcoesses)
 @inline function makemajump(rx; combinatoric_ratelaw = true)
     @unpack rate, substrates, substoich, netstoich = rx
+    zeroorder = (length(substoich) == 0)
 
     # if nothing changes in the reaction (e.g. `X --> X`), throws an error
     net_stoch = filter(p -> !drop_dynamics(p[1]), netstoich)
@@ -345,7 +346,6 @@ end
 
     # for non-zeroth order reaction (e.g. `X1 --> X2`), and if we use a combinatorial rate law,
     # modifies the rate with its mass action coefficient (no need if coefficient is 1)
-    zeroorder = (length(substoich) == 0)
     if (!zeroorder) && combinatoric_ratelaw
         coef = prod(factorial, substoich)
         (!isone(coef)) && (rate /= coef)
@@ -357,7 +357,7 @@ end
 
 # recursively visit each neighbor's rooted tree and mark everything in it as vrj
 function dfs_mark!(isvrjvec, visited, depgraph, i)
-    # Marks the current node as visited. Finds all neighbours.
+    # marks the current node as visited. Finds all neighbours
     visited[i] = true
     nhbrs = depgraph[i]
 
@@ -373,7 +373,7 @@ function dfs_mark!(isvrjvec, visited, depgraph, i)
 end
 
 # get_depgraph(rs)[i] is the list of reactions with rates depending on species changed by
-# i'th reaction
+# i'th reaction.
 function get_depgraph(rs)
     jdeps = asgraph(rs)
     vdeps = variable_dependencies(rs)
@@ -563,44 +563,46 @@ Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.ODESystem`.
 
 Keyword args and default values:
-- `combinatoric_ratelaws = true` uses factorial scaling factors in calculating the rate law,
- i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
- `combinatoric_ratelaws = false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
- ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
- itself defaults to true).
-- `remove_conserved = false`, if set to `true` will calculate conservation laws of the
- underlying set of reactions (ignoring coupled equations), and then apply them to reduce
- the number of equations.
+- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate law,
+  i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
+  `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
+  ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
+  itself defaults to true).
+- `remove_conserved=false`, if set to `true` will calculate conservation laws of the
+  underlying set of reactions (ignoring constraint equations), and then apply them to reduce
+  the number of equations.
+- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
+  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs),
-                      combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-                      include_zero_odes = true, remove_conserved = false, checks = false,
-                      default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
-                      kwargs...)
+        combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+        include_zero_odes = true, remove_conserved = false, remove_conserved_warn = true,
+        checks = false, default_u0 = Dict(), default_p = Dict(),
+        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, ODESystem)
+    check_cons_warning(remove_conserved, remove_conserved_warn)
 
-    # Flattens the system. Next (if `remove_conserved = true`) uses `conservationlaws` to compute
-    # these (these are then automatically cached in `rs`).
     flatrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(flatrs)
 
     # Find independent unknowns and species. Assembles the drift ODE, and then the full set of equations.
     # Here, `add_coupled_eqs!` performs some assembly in addition to adding coupled equations.
     ind_us, ind_sps = get_indep_us(flatrs, remove_conserved)
-    eqs = assemble_drift(flatrs, ind_sps; combinatoric_ratelaws, remove_conserved,
-                         include_zero_odes)
-    eqs, us, ps, obs, defs = add_coupled_eqs!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
+    eqs = assemble_drift(fullrs, ind_sps; combinatoric_ratelaws, remove_conserved,
+        include_zero_odes)
+    eqs, us, ps, obs, defs = addconstraints!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
 
     ODESystem(eqs, get_iv(flatrs), us, ps;
-              observed = obs,
-              name,
-              defaults = _merge(defaults, defs),
-              checks,
-              continuous_events = MT.get_continuous_events(flatrs),
-              discrete_events = MT.get_discrete_events(flatrs),
-              kwargs...)
+        observed = obs,
+        name,
+        defaults = _merge(defaults, defs),
+        checks,
+        continuous_events = MT.get_continuous_events(flatrs),
+        discrete_events = MT.get_discrete_events(flatrs),
+        kwargs...)
 end
 
 """
@@ -611,14 +613,16 @@ Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.NonlinearSystem`.
 
 Keyword args and default values:
-- `combinatoric_ratelaws = true` uses factorial scaling factors in calculating the rate law,
- i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
- `combinatoric_ratelaws = false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
- ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
- itself defaults to true).
-- `remove_conserved = false`, if set to `true` will calculate conservation laws of the
- underlying set of reactions (ignoring coupled equations), and then apply them to reduce
- the number of equations.
+- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate law,
+  i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
+  `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
+  ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
+  itself defaults to true).
+- `remove_conserved=false`, if set to `true` will calculate conservation laws of the
+  underlying set of reactions (ignoring constraint equations), and then apply them to reduce
+  the number of equations.
+- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
+  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
@@ -634,8 +638,6 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
         error("Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depend on $(get_iv(rs))) to a `NonlinearSystem`. This is not possible. if you are intending to compute system steady states, consider creating and solving a `SteadyStateProblem.")
     end
 
-    # Flattens the system. Next (if `remove_conserved = true`) uses `conservationlaws` to compute
-    # these (and automatically cache these in `rs`).
     flatrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(flatrs)
 
@@ -643,8 +645,8 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     # Here, `add_coupled_eqs!` performs some assembly in addition to adding coupled equations.
     ind_us, ind_sps = get_indep_us(flatrs, remove_conserved)
     eqs = assemble_drift(flatrs, ind_sps; combinatoric_ratelaws, remove_conserved,
-                         as_odes = false, include_zero_odes)
-    eqs, us, ps, obs, defs = add_coupled_eqs!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
+        as_odes = false, include_zero_odes)
+    eqs, us, ps, obs, defs = addconstraints!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
 
     # Throws a warning if there are differential equations in non-standard format.
     # Next, sets all differential terms to `0`.
@@ -652,11 +654,11 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     eqs = [remove_diffs(eq.lhs) ~ remove_diffs(eq.rhs) for eq in eqs]
 
     NonlinearSystem(eqs, us, ps;
-                    name,
-                    observed = obs,
-                    defaults = _merge(defaults, defs),
-                    checks,
-                    kwargs...)
+        name,
+        observed = obs,
+        defaults = _merge(defaults, defs),
+        checks,
+        kwargs...)
 end
 
 # Ideally, when `ReactionSystem`s are converted to `NonlinearSystem`s, any coupled ODEs should be 
@@ -671,18 +673,18 @@ function nonlinear_convert_differentials_check(rs::ReactionSystem)
         # If the lhs upper-level function is not a differential w.r.t. time.
         # If the content of the differential is not a variable (and nothing more).
         # If either of this is a case, throw the warning.
-        if Symbolics._occursin(Symbolics.is_derivative, eq.rhs) ||
-                               !Symbolics.istree(eq.lhs) ||
-                               !isequal(Symbolics.operation(eq.lhs), Differential(get_iv(rs))) || 
-                               (length(arguments(eq.lhs)) != 1) ||
-                               !any(isequal(arguments(eq.lhs)[1]), nonspecies(rs))
+        if hasnode(Symbolics.is_derivative, eq.rhs) ||
+           !Symbolics.is_derivative(eq.lhs) ||
+           !isequal(Symbolics.operation(eq.lhs), Differential(get_iv(rs))) ||
+           (length(arguments(eq.lhs)) != 1) ||
+           !any(isequal(arguments(eq.lhs)[1]), nonspecies(rs))
             error("You are attempting to convert a `ReactionSystem` coupled with differential equations to a `NonlinearSystem`. However, some of these differentials are not of the form `D(x) ~ ...` where: 
-                  (1) The left-hand side is a differential of a single variable with respect to the time-independent variable, and
-                  (2) The right-hand side does not contain any differentials.
-                  This is generally not permitted.
-                                  
-                  If you still would like to perform this conversion, please use the `all_differentials_permitted = true` option. In this case, all differential will be set to `0`. 
-                  However, it is recommended to proceed with caution to ensure that the produced nonlinear equation makes sense for your intended application."
+                    (1) The left-hand side is a differential of a single variable with respect to the time independent variable, and
+                    (2) The right-hand side does not contain any differentials.
+                This is generally not permitted.
+
+                If you still would like to perform this conversion, please use the `all_differentials_permitted = true` option. In this case, all differentials will be set to `0`. 
+                However, it is recommended to proceed with caution to ensure that the produced nonlinear equation makes sense for your intended application."
             )
         end
     end
@@ -696,29 +698,28 @@ Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.SDESystem`.
 
 Notes:
-- `combinatoric_ratelaws = true` uses factorial scaling factors in calculating the rate law,
- i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
- `combinatoric_ratelaws = false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
- ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
- itself defaults to true).
-- `remove_conserved = false`, if set to `true` will calculate conservation laws of the
- underlying set of reactions (ignoring coupled equations), and then apply them to reduce
- the number of equations.
-- Does not currently support `ReactionSystem`s that includes coupled algebraic or
- differential equations.
+- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate law,
+  i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
+  `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
+  ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
+  itself defaults to true).
+- `remove_conserved=false`, if set to `true` will calculate conservation laws of the
+  underlying set of reactions (ignoring constraint equations), and then apply them to reduce
+  the number of equations.
+- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
+  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
-                      name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-                      include_zero_odes = true, checks = false, remove_conserved = false,
-                      default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
-                      kwargs...)
+        name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+        include_zero_odes = true, checks = false, remove_conserved = false,
+        remove_conserved_warn = true, default_u0 = Dict(), default_p = Dict(),
+        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, SDESystem)
     check_cons_warning(remove_conserved, remove_conserved_warn)
 
-    # Flattens the system. Next (if `remove_conserved = true`) uses `conservationlaws` to compute
-    # these (and automatically cache these in `rs`).
     flatrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(flatrs)
 
@@ -727,9 +728,14 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
     # Here, `add_coupled_eqs!` performs some assembly in addition to adding coupled equations.
     ind_us, ind_sps = get_indep_us(flatrs, remove_conserved)
     eqs = assemble_drift(flatrs, ind_sps; combinatoric_ratelaws, include_zero_odes,
-                         remove_conserved)
-    noiseeqs = assemble_diffusion(flatrs, ind_us, ind_sps; combinatoric_ratelaws, remove_conserved)
-    eqs, us, ps, obs, defs = add_coupled_eqs!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
+        remove_conserved)
+    noiseeqs = assemble_diffusion(flatrs, ind_us, ind_sps;
+        combinatoric_ratelaws, remove_conserved)
+    eqs, us, ps, obs, defs = addconstraints!(eqs, flatrs, ind_us, ind_sps; remove_conserved)
+
+    if any(isbc, get_unknowns(flatrs))
+        @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead."
+    end
 
     SDESystem(eqs, noiseeqs, get_iv(flatrs), us, ps;
         observed = obs,
@@ -743,52 +749,53 @@ end
 
 """
 ```julia
-Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; combinatoric_ratelaws = true)
+Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
 ```
 
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.JumpSystem`.
 
 Notes:
-- `combinatoric_ratelaws = true` uses binomials in calculating the rate law, i.e. for `2S ->
- 0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaws = false` then
- the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling factor.
- Defaults to the value given when the `ReactionSystem` was constructed (which itself
- defaults to true).
+- `combinatoric_ratelaws=true` uses binomials in calculating the rate law, i.e. for `2S ->
+  0` at rate `k` the ratelaw would be `k*S*(S-1)/2`. If `combinatoric_ratelaws=false` then
+  the ratelaw is `k*S*(S-1)`, i.e. the rate law is not normalized by the scaling factor.
+  Defaults to the value given when the `ReactionSystem` was constructed (which itself
+  defaults to true).
 - Does not currently support `ReactionSystem`s that includes coupled algebraic or
- differential equations.
+  differential equations.
 - Does not currently support continuous events as these are not supported by
- `ModelingToolkit.JumpSystems`.
+  `ModelingToolkit.JumpSystems`.
 """
 function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs),
-                      combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-                      remove_conserved = nothing, checks = false,
-                      default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
-                      kwargs...)
-    # Error checks.
+        combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+        remove_conserved = nothing, checks = false,
+        default_u0 = Dict(), default_p = Dict(),
+        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        kwargs...)
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, JumpSystem)
     (remove_conserved !== nothing) &&
         throw(ArgumentError("Catalyst does not support removing conserved species when converting to JumpSystems."))
 
-    # Flattens the system. Performs post-flattening errors and warning checks.
     flatrs = Catalyst.flatten(rs)
     error_if_coupled(JumpSystem, flatrs)
+
     (length(MT.continuous_events(flatrs)) > 0) &&
         (@warn "continuous_events will be dropped as they are not currently supported by JumpSystems.")
 
+    eqs = assemble_jumps(flatrs; combinatoric_ratelaws)
+
     # Computes species, parameters, and equations. Adds boundary condition species to species vector.
-    _, ind_sps = get_indep_us(flatrs, remove_conserved)
-    any(isbc, get_unknowns(flatrs)) && (ind_sps = vcat(ind_sps, filter(isbc, get_unknowns(flatrs))))
+    _, ind_sps = get_indep_us(flatrs)
+    any(isbc, get_unknowns(flatrs)) && (ind_sps = vcat(sts, filter(ind_sps, get_unknowns(flatrs))))
     ps = get_ps(flatrs)
-    eqs = assemble_jumps(flatrs; combinatoric_ratelaws)
-    
+
     JumpSystem(eqs, get_iv(flatrs), ind_sps, ps;
-               observed = MT.observed(flatrs),
-               name,
-               defaults = _merge(defaults, MT.defaults(flatrs)),
-               checks,
-               discrete_events = MT.discrete_events(flatrs),
-               kwargs...)
+        observed = MT.observed(flatrs),
+        name,
+        defaults = _merge(defaults, MT.defaults(flatrs)),
+        checks,
+        discrete_events = MT.discrete_events(flatrs),
+        kwargs...)
 end
 
 ### Problems ###
@@ -806,7 +813,7 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
         remove_conserved, remove_conserved_warn)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
-    if structural_simplify 
+    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.")
@@ -835,25 +842,25 @@ end
 
 # SDEProblem from AbstractReactionNetwork
 function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
-                               p = DiffEqBase.NullParameters(), args...;
-                               name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-                               include_zero_odes = true, checks = false, check_length = false,
-                               remove_conserved = false, structural_simplify = false, kwargs...)
+        p = DiffEqBase.NullParameters(), args...;
+        name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+        include_zero_odes = true, checks = false, check_length = false, remove_conserved = false,
+        remove_conserved_warn = true, structural_simplify = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws,
-                      include_zero_odes, checks, remove_conserved)
-    p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs))
+        include_zero_odes, checks, remove_conserved, remove_conserved_warn)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
-    if structural_simplify 
+    if structural_simplify
         sde_sys = MT.structural_simplify(sde_sys)
     elseif has_alg_equations(rs)
         error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify = true` within `ODEProblem` call.")
     else
         sde_sys = complete(sde_sys)
     end
-    
+
+    p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs))
     return SDEProblem(sde_sys, u0map, tspan, pmap, args...; check_length,
         noise_rate_prototype = p_matrix, kwargs...)
 end
@@ -894,10 +901,10 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0,
         remove_conserved, remove_conserved_warn)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
-    if structural_simplify 
-        osys = MT.structural_simplify(osys)
+    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.")
+        error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.")
     else
         osys = complete(osys)
     end
@@ -913,11 +920,7 @@ function _symbol_to_var(sys, sym)
     if hasproperty(sys, sym)
         var = getproperty(sys, sym, namespace = false)
     else
-        # If sys does not contain sym it is possible that sys.subsystem contains it (or a 
-        # sub-sub system etc.). Here, "₊" splits the system and subsystem names in the symbol.
-        # We thus split the symbol to find, and go down the subsystem tree to attempt to find
-        # the symbol. 
-        strs = split(String(sym), "₊")   # need to check if this should be split of not!!!
+        strs = split(String(sym), ModelingToolkit.NAMESPACE_SEPARATOR)   # need to check if this should be split of not!!!
         if length(strs) > 1
             var = getproperty(sys, Symbol(strs[1]), namespace = false)
             for str in view(strs, 2:length(strs))
@@ -931,7 +934,7 @@ function _symbol_to_var(sys, sym)
 end
 
 """
- symmap_to_varmap(sys, symmap)
+    symmap_to_varmap(sys, symmap)
 
 Given a system and map of `Symbol`s to values, generates a map from
 corresponding symbolic variables/parameters to the values that can be used to
@@ -940,11 +943,11 @@ pass initial conditions and parameter mappings.
 For example,
 ```julia
 sir = @reaction_network sir begin
- β, S + I --> 2I
- ν, I --> R
+    β, S + I --> 2I
+    ν, I --> R
 end
 subsys = @reaction_network subsys begin
- k, A --> B
+    k, A --> B
 end
 @named sys = compose(sir, [subsys])
 ```
@@ -952,15 +955,15 @@ gives
 ```
 Model sys with 3 equations
 Unknowns (5):
- S(t)
- I(t)
- R(t)
- subsys₊A(t)
- subsys₊B(t)
+  S(t)
+  I(t)
+  R(t)
+  subsys₊A(t)
+  subsys₊B(t)
 Parameters (3):
- β
- ν
- subsys₊k
+  β
+  ν
+  subsys₊k
 ```
 to specify initial condition and parameter mappings from *symbols* we can use
 ```julia
@@ -972,14 +975,12 @@ pmap   = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0])
 
 Notes:
 - Any `Symbol`, `sym`, within `symmap` must be a valid field of `sys`. i.e.
- `sys.sym` must be defined.
+  `sys.sym` must be defined.
 """
 function symmap_to_varmap(sys, symmap::Tuple)
-    # If all inputs are of the form `:sym => value`, apply `_symbol_to_var` function. If not,
-    # pass input through unchanged.
     if all(p -> p isa Pair{Symbol}, symmap)
         return ((_symbol_to_var(sys, sym) => val for (sym, val) in symmap)...,)
-    else
+    else  # if not all entries map a symbol to value pass through
         return symmap
     end
 end
@@ -999,11 +1000,11 @@ symmap_to_varmap(sys, symmap) = symmap
 ### Other Conversion-related Functions ###
 
 # the following function is adapted from SymbolicUtils.jl v.19
-# later on (Spetember 2023) modified by Torkel and Shashi (now assumes input not on polynomial form, which is handled elsewhere, the previous version failed in these cases anyway).
+# later on (Spetember 2023) modified by Torkel and Shashi (now assumes input not on polynomial form, which is handled elsewhere, previous version failed in these cases anyway).
 # Copyright (c) 2020: Shashi Gowda, Yingbo Ma, Mason Protter, Julia Computing.
 # MIT license
 """
- to_multivariate_poly(polyeqs::AbstractVector{BasicSymbolic{Real}})
+    to_multivariate_poly(polyeqs::AbstractVector{BasicSymbolic{Real}})
 
 Convert the given system of polynomial equations to multivariate polynomial representation.
 For example, this can be used in HomotopyContinuation.jl functions.
@@ -1021,4 +1022,4 @@ function to_multivariate_poly(polyeqs::AbstractVector{Symbolics.BasicSymbolic{Re
     end
 
     ps
-end
+end
\ No newline at end of file

From 24587b8c50dc27399118b6cb91bbae62addcce87 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 27 Jul 2024 19:03:46 -0400
Subject: [PATCH 4/4] up

---
 src/reactionsystem_conversions.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index e9787646f4..7393a95622 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -53,10 +53,10 @@ drop_dynamics(spec) = isconstant(spec) || isbc(spec) || (!isspecies(spec))
 
 # for a `ReactionSystem`, create a vector with all the right-hand side expressions in the
 # corresponding ODE system (used by the `assemble_drift` function)
-function assemble_oderhs(rs, isps; combinatoric_ratelaws = true, remove_conserved = false)
+function assemble_oderhs(rs, ind_sps; combinatoric_ratelaws = true, remove_conserved = false)
     # initiates the right-hand side equations as expressions with `0` only.
-    species_to_idx = Dict(x => i for (i, x) in enumerate(isps))
-    rhsvec = Any[0 for _ in isps]
+    species_to_idx = Dict(x => i for (i, x) in enumerate(ind_sps))
+    rhsvec = Any[0 for _ in ind_sps]
 
     # if we have eliminated conservation laws, create a dictionary taking each species eliminated
     # in this process to the expression it is replaced by (if not, creates an empty dict)