diff --git a/Cargo.lock b/Cargo.lock
index c57089af0..736ab9f3c 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -380,6 +380,7 @@ dependencies = [
  "ustr",
  "uuid",
  "wasm-bindgen",
+ "web-sys",
 ]
 
 [[package]]
diff --git a/packages/catlog-wasm/src/analyses.rs b/packages/catlog-wasm/src/analyses.rs
index db5c1dfc8..e942e846a 100644
--- a/packages/catlog-wasm/src/analyses.rs
+++ b/packages/catlog-wasm/src/analyses.rs
@@ -20,3 +20,7 @@ pub struct LinearODEModelData(pub analyses::ode::LinearODEProblemData<Uuid>);
 #[derive(Serialize, Deserialize, Tsify)]
 #[tsify(into_wasm_abi, from_wasm_abi)]
 pub struct MassActionModelData(pub analyses::ode::MassActionProblemData<Uuid>);
+
+#[derive(Serialize, Deserialize, Tsify)]
+#[tsify(into_wasm_abi, from_wasm_abi)]
+pub struct SwitchingMassActionModelData(pub analyses::ode::SwitchingMassActionProblemData<Uuid>);
diff --git a/packages/catlog-wasm/src/theories.rs b/packages/catlog-wasm/src/theories.rs
index 12835c8ac..0ad893bb8 100644
--- a/packages/catlog-wasm/src/theories.rs
+++ b/packages/catlog-wasm/src/theories.rs
@@ -13,7 +13,7 @@ use catlog::dbl::theory;
 use catlog::one::Path;
 use catlog::stdlib::{analyses, models, theories, theory_morphisms};
 
-use super::model_morphism::{MotifsOptions, motifs};
+use super::model_morphism::{motifs, MotifsOptions};
 use super::{analyses::*, model::DblModel, theory::DblTheory};
 
 /// The empty or initial theory.
@@ -337,6 +337,39 @@ impl ThSymMonoidalCategory {
     }
 }
 
+/// The theory of state/aux interactions
+#[wasm_bindgen]
+pub struct ThModalStateAuxCategory(Rc<theory::UstrModalDblTheory>);
+
+#[wasm_bindgen]
+impl ThModalStateAuxCategory {
+    #[wasm_bindgen(constructor)]
+    pub fn new() -> Self {
+        Self(Rc::new(theories::th_modal_state_aux()))
+    }
+
+    #[wasm_bindgen]
+    pub fn theory(&self) -> DblTheory {
+        DblTheory(self.0.clone().into())
+    }
+
+    /// Simulates a mass-action ODE system with additional configurations.
+    #[wasm_bindgen(js_name = "massAction")]
+    pub fn state_aux_mass_action(
+        &self,
+        model: &DblModel,
+        data: SwitchingMassActionModelData,
+    ) -> Result<ODEResult, String> {
+        Ok(ODEResult(
+            analyses::ode::PetriNetMassActionFunctionAnalysis::default()
+                .build_numerical_system(model.modal()?, data.0)
+                .solve_with_defaults()
+                .map_err(|err| format!("{err:?}"))
+                .into(),
+        ))
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
diff --git a/packages/catlog/Cargo.toml b/packages/catlog/Cargo.toml
index 938839524..716afac63 100644
--- a/packages/catlog/Cargo.toml
+++ b/packages/catlog/Cargo.toml
@@ -28,6 +28,7 @@ tsify = { version = "0.5", features = ["js"], optional = true }
 ustr = "1"
 uuid = { version = "=1.11", optional = true }
 wasm-bindgen = { version = "0.2.100", optional = true }
+web-sys = { version = "0.3.77", features = ["console"] }
 
 [dev-dependencies]
 expect-test = "1.5"
diff --git a/packages/catlog/src/dbl/modal/model.rs b/packages/catlog/src/dbl/modal/model.rs
index d43cba7be..fbf7592eb 100644
--- a/packages/catlog/src/dbl/modal/model.rs
+++ b/packages/catlog/src/dbl/modal/model.rs
@@ -17,7 +17,7 @@ use crate::validate::{self, Validate};
 use crate::{one::computad::*, one::*, zero::*};
 
 /// Object in a model of a modal double theory.
-#[derive(Clone, Debug, PartialEq, Eq, From)]
+#[derive(Clone, Debug, PartialEq, Eq, From, Hash)]
 pub enum ModalOb<Id, ThId> {
     /// Generating object.
     #[from]
diff --git a/packages/catlog/src/dbl/modal/theory.rs b/packages/catlog/src/dbl/modal/theory.rs
index 6878c1575..113fdcaaa 100644
--- a/packages/catlog/src/dbl/modal/theory.rs
+++ b/packages/catlog/src/dbl/modal/theory.rs
@@ -54,7 +54,7 @@ double category of sets admits, besides the [plain](Self::Plain) list double
 monad, a number of variations decorating the spans of lists with extra
 combinatorial data.
  */
-#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)]
 pub enum List {
     /// Lists of objects and morphisms (of same length).
     Plain,
diff --git a/packages/catlog/src/one/graph.rs b/packages/catlog/src/one/graph.rs
index b162afde9..89d38061e 100644
--- a/packages/catlog/src/one/graph.rs
+++ b/packages/catlog/src/one/graph.rs
@@ -194,20 +194,20 @@ of [`FinGraph`].
  */
 pub trait ColumnarFinGraph:
     ColumnarGraph<
-        Vertices: FinSet<Elem = Self::V>,
-        Edges: FinSet<Elem = Self::E>,
-        Src: Column<Dom = Self::E, Cod = Self::V>,
-        Tgt: Column<Dom = Self::E, Cod = Self::V>,
-    >
+    Vertices: FinSet<Elem = Self::V>,
+    Edges: FinSet<Elem = Self::E>,
+    Src: Column<Dom = Self::E, Cod = Self::V>,
+    Tgt: Column<Dom = Self::E, Cod = Self::V>,
+>
 {
 }
 
 /// A columnar graph with mutable columns.
 pub trait MutColumnarGraph:
     ColumnarGraph<
-        Src: MutMapping<Dom = Self::E, Cod = Self::V>,
-        Tgt: MutMapping<Dom = Self::E, Cod = Self::V>,
-    >
+    Src: MutMapping<Dom = Self::E, Cod = Self::V>,
+    Tgt: MutMapping<Dom = Self::E, Cod = Self::V>,
+>
 {
     /// Variant of [`src_map`](ColumnarGraph::src_map) that returns a mutable
     /// reference.
diff --git a/packages/catlog/src/one/graph_algorithms.rs b/packages/catlog/src/one/graph_algorithms.rs
index 6395b16c1..7fd15fc64 100644
--- a/packages/catlog/src/one/graph_algorithms.rs
+++ b/packages/catlog/src/one/graph_algorithms.rs
@@ -1,6 +1,6 @@
 //! Algorithms on graphs.
 
-use std::collections::{HashSet, VecDeque};
+use std::collections::{HashMap, HashSet, VecDeque};
 use std::hash::Hash;
 
 use super::graph::*;
@@ -168,6 +168,136 @@ where
     result
 }
 
+fn out_neighbors<G>(graph: &G, v: &G::V) -> impl Iterator<Item = G::V>
+where
+    G: FinGraph,
+    G::V: Hash,
+{
+    graph.out_edges(v).map(|e| graph.tgt(&e))
+}
+
+fn in_neighbors<G>(graph: &G, v: &G::V) -> impl Iterator<Item = G::V>
+where
+    G: FinGraph,
+    G::V: Hash,
+{
+    graph.in_edges(v).map(|e| graph.src(&e))
+}
+
+#[derive(Clone, Debug)]
+struct VisitMap {
+    visited: Vec<bool>,
+}
+
+impl VisitMap {
+    fn visit(&mut self, idx: usize) -> bool {
+        let previous = self.visited[idx];
+        self.visited[idx] = true;
+        !previous
+    }
+
+    fn is_visited(&self, idx: usize) -> bool {
+        self.visited[idx]
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct Dfs<G>
+where
+    G: FinGraph,
+    G::V: Hash,
+{
+    //
+    pub stack: Vec<G::V>,
+    //
+    pub discovered: VisitMap,
+}
+
+// TODO
+// 1. replace String with Cycle error
+/** Computes a topological sorting for a given graph.
+
+This algorithm was borrowed from `petgraph`.
+ */
+pub fn toposort<G>(graph: &G) -> Result<Vec<G::V>, String>
+where
+    G: FinGraph,
+    G::V: Hash + std::fmt::Debug,
+{
+    // XXX dont clone
+    let n = graph.vertices().collect::<Vec<_>>().len();
+    let mut discovered = VisitMap {
+        visited: vec![false; n.clone()],
+    };
+    let mut finished = VisitMap {
+        visited: vec![false; n.clone()],
+    };
+    let mut finish_stack: Vec<G::V> = Vec::new();
+    let mut stack = Vec::new();
+
+    // we shouldn't need to do this
+    let gmap: HashMap<_, _> = HashMap::from_iter(graph.vertices().enumerate().map(|(k, v)| (v, k)));
+
+    for (idx, v) in graph.vertices().enumerate() {
+        if discovered.is_visited(idx) {
+            continue;
+        }
+        stack.push(v);
+        while let Some(nx) = stack.clone().last() {
+            if discovered.visit(gmap[&nx]) {
+                for succ in out_neighbors(graph, &nx) {
+                    if succ == *nx {
+                        return Err("self cycle".to_owned());
+                    }
+                    if !discovered.is_visited(gmap[&succ]) {
+                        stack.push(succ);
+                    }
+                }
+            } else {
+                stack.pop();
+                if finished.visit(gmap[&nx]) {
+                    finish_stack.push(nx.clone());
+                }
+            }
+        }
+    }
+    finish_stack.reverse();
+
+    // dfs.reset(g);
+    let mut discovered = VisitMap {
+        visited: vec![false; n.clone()],
+    };
+    for i in &finish_stack {
+        // dfs.move_to(i);
+        stack.clear();
+        stack.push(i.clone());
+        //
+        let mut cycle = false;
+        while let Some(j) = {
+            let mut out = None;
+            while let Some(node) = stack.pop() {
+                if discovered.visit(gmap[&node]) {
+                    for succ in in_neighbors(graph, &node) {
+                        if !discovered.is_visited(gmap[&succ]) {
+                            stack.push(succ);
+                        }
+                    }
+                    out = Some(node);
+                    break;
+                }
+            }
+            out
+        } {
+            if cycle {
+                return Err(format!("cycle detected involving node {:#?}", j).to_owned());
+            }
+            cycle = true;
+        }
+    }
+
+    Ok(finish_stack)
+}
+
 #[cfg(test)]
 mod tests {
     use super::GraphElem::*;
@@ -231,4 +361,46 @@ mod tests {
         let g = SkelGraph::cycle(1);
         assert_eq!(spec_order_all(&g), vec![Vertex(0), Edge(0)]);
     }
+
+    #[test]
+    fn toposorting() {
+        let mut g = SkelGraph::path(5);
+        assert_eq!(toposort(&g), Ok(vec![0, 1, 2, 3, 4]));
+
+        let mut g = SkelGraph::path(3);
+        g.add_vertices(1);
+        let _ = g.add_edge(2, 3);
+        let _ = g.add_edge(3, 0);
+        assert_eq!(toposort(&g), Err("cycle detected involving node 3".to_owned()));
+
+        let g = SkelGraph::triangle();
+        assert_eq!(toposort(&g), Ok(vec![0, 1, 2]));
+
+        let mut g = SkelGraph::path(4);
+        g.add_vertices(2);
+        let _ = g.add_edge(1, 4);
+        let _ = g.add_edge(4, 3);
+        let _ = g.add_edge(5, 2);
+        assert_eq!(toposort(&g), Ok(vec![5, 0, 1, 2, 4, 3]));
+
+        let mut g: HashGraph<u8, &str> = Default::default();
+        g.add_vertices(vec![0, 1, 2, 3, 4, 5]);
+        g.add_edge("0-1", 0, 1);
+        g.add_edge("1-2", 1, 2);
+        g.add_edge("2-3", 2, 3);
+        g.add_edge("1-4", 1, 4);
+        g.add_edge("4-3", 4, 3);
+        g.add_edge("5-2", 5, 2);
+        // TODO non-deterministic
+        // assert_eq!(toposort(&g), Ok(vec![5, 0, 1, 2, 4, 3]));
+    }
+
+    #[test]
+    #[test]
+    fn neighbors() {
+        let g = SkelGraph::triangle();
+        let out_neighbors = &g.out_edges(&1).map(|e| g.tgt(&e)).collect::<Vec<_>>();
+        // dbg!(&out_neighbors);
+        assert!(true);
+    }
 }
diff --git a/packages/catlog/src/simulate/ode/linear_ode.rs b/packages/catlog/src/simulate/ode/linear_ode.rs
index 6376e11ff..f786e53dc 100644
--- a/packages/catlog/src/simulate/ode/linear_ode.rs
+++ b/packages/catlog/src/simulate/ode/linear_ode.rs
@@ -1,6 +1,7 @@
 //! Constant-coefficient linear first-order differential equations.
 
 use nalgebra::{DMatrix, DVector};
+use std::collections::HashMap;
 
 #[cfg(test)]
 use super::ODEProblem;
@@ -30,6 +31,48 @@ impl ODESystem for LinearODESystem {
     }
 }
 
+// --------------------------------
+
+// TODO move into a context where Id is defined
+#[derive(Clone, Debug)]
+pub enum Condition {
+    Geq(usize, usize),
+}
+
+impl Condition {
+    pub fn call(&self, x: &DVector<f32>) -> bool {
+        match self {
+            Condition::Geq(left, right) => x.row(*left) >= x.row(*right),
+        }
+    }
+}
+
+/**
+ */
+#[derive(Clone)]
+pub struct HybridLinearODESystem {
+    /// Linear ODEs
+    pub subsystems: HashMap<Condition, LinearODESystem>,
+}
+
+impl HybridLinearODESystem {
+    // TODO I dont like checking for the first condition to be met.
+    pub fn get_subsystem(&self, x: &DVector<f32>) -> Option<&LinearODESystem> {
+        self.subsystems
+            .iter()
+            .find(|(cond, _)| cond.call(&x.clone()))
+            .map(|(_, sys)| sys)
+    }
+}
+
+impl ODESystem for HybridLinearODESystem {
+    fn vector_field(&self, dx: &mut DVector<f32>, x: &DVector<f32>, _t: f32) {
+        if let Some(system) = self.get_subsystem(x) {
+            system.vector_field(dx, x, _t)
+        }
+    }
+}
+
 #[cfg(test)]
 pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem<LinearODESystem> {
     use nalgebra::{dmatrix, dvector};
diff --git a/packages/catlog/src/simulate/ode/polynomial.rs b/packages/catlog/src/simulate/ode/polynomial.rs
index cd6bf65e1..931197340 100644
--- a/packages/catlog/src/simulate/ode/polynomial.rs
+++ b/packages/catlog/src/simulate/ode/polynomial.rs
@@ -1,13 +1,19 @@
 //! Polynomial differential equations.
 
-use std::collections::BTreeMap;
-use std::fmt::Display;
+use std::cmp::{Eq, Ord};
+use std::collections::{BTreeMap, HashMap};
+use std::fmt::{Debug, Display};
+use std::hash::Hash;
 use std::ops::Add;
 
+use itertools::{EitherOrBoth::*, Itertools};
+
 use derivative::Derivative;
 use nalgebra::DVector;
 use num_traits::{One, Pow, Zero};
 
+use crate::stdlib::analyses::ode::{ComputeGraph, EligibleFunctions, Output};
+
 #[cfg(test)]
 use super::ODEProblem;
 use super::ODESystem;
@@ -127,11 +133,32 @@ where
 Such a system is ready for use in numerical solvers: the coefficients are
 floating point numbers and the variables are consecutive integer indices.
  */
+#[derive(Clone)]
 pub struct NumericalPolynomialSystem<Exp> {
     /// Components of the vector field.
     pub components: Vec<Polynomial<usize, f32, Exp>>,
 }
 
+impl<Exp> Add for NumericalPolynomialSystem<Exp>
+where
+    Exp: Ord,
+{
+    type Output = Self;
+
+    fn add(self, rhs: Self) -> Self::Output {
+        let components = self
+            .components
+            .into_iter()
+            .zip_longest(rhs.components.into_iter())
+            .map(|pair| match pair {
+                Both(l, r) => l + r,
+                Left(x) | Right(x) => x,
+            })
+            .collect::<Vec<_>>();
+        Self { components }
+    }
+}
+
 impl<Exp> ODESystem for NumericalPolynomialSystem<Exp>
 where
     Exp: Clone + Ord,
@@ -144,6 +171,154 @@ where
     }
 }
 
+/**
+ */
+#[derive(Clone, Derivative)]
+#[derivative(Default(bound = ""))]
+pub struct NumericalPolynomialSwitchingSystem<Id, Exp>
+where
+    Id: Hash + Clone + Eq + Debug,
+    Exp: Clone,
+{
+    ///
+    pub ob_index: BTreeMap<Id, usize>,
+
+    // TODO public?
+    /// Null model
+    pub null_model: Option<NumericalPolynomialSystem<Exp>>,
+
+    /// Components of a switching system.
+    pub subsystems: Vec<(ComputeGraph<Id>, NumericalPolynomialSystem<Exp>)>,
+
+    /// Analysis data
+    pub functions: HashMap<Id, String>,
+}
+
+impl<Id, Exp> NumericalPolynomialSwitchingSystem<Id, Exp>
+where
+    Id: Hash + Clone + Eq + Debug + Ord,
+    Exp: Ord + Clone,
+{
+    // TODO do we need a `null_model`?
+    fn get_current(self, x: DVector<f32>) -> NumericalPolynomialSystem<Exp> {
+        self.subsystems
+            .into_iter()
+            .filter(|(graph, _)| graph.compute(&x, self.functions.clone(), self.ob_index.clone()))
+            .fold(self.null_model.unwrap(), |acc, (_, sys)| acc + sys)
+    }
+}
+
+// impl<Id: Hash + Clone + Eq + Debug, Exp> From<NumericalPolynomialSystem<Exp>>
+//     for NumericalPolynomialSwitchingSystem<Id, Exp>
+// where
+//     Exp: Clone,
+// {
+//     fn from(nps: NumericalPolynomialSystem<Exp>) -> Self {
+//         let subsystems = vec![(ComputeGraph::<Id>::new(), nps)];
+//         let functions = HashMap::from([(ustr::ustr("comparator"), EligibleFunctions::Geq())]);
+//         NumericalPolynomialSwitchingSystem {
+//             subsystems,
+//             functions,
+//         }
+//     }
+// }
+
+impl<Id> ComputeGraph<Id>
+where
+    Id: Eq + Hash + Debug + Clone + Ord,
+{
+    //
+    fn compute(
+        &self,
+        x: &DVector<f32>,
+        functions: HashMap<Id, String>,
+        ob_index: BTreeMap<Id, usize>,
+    ) -> bool {
+        // binding name
+        let mut bindings: Vec<(_, Output)> = Default::default();
+        for var in self.toposort.iter() {
+            if self.obs.contains(&var) {
+                let index = ob_index[var];
+                bindings.push((var, Output::Float(x[index])));
+            } else if self.borrows.keys().contains(&var) {
+                let index = ob_index[&self.borrows[var]];
+                bindings.push((var, Output::Float(x[index])));
+            } else if let Some(val) = self.auxes.get(&var) {
+                bindings.push((var, Output::Float(*val)));
+            } else if self.funcs.keys().contains(&var) {
+                let args = &self.funcs[&var]
+                    .clone()
+                    .into_iter()
+                    .map(|arg| self.fetch(x, arg, ob_index.clone()))
+                    .collect::<Vec<_>>();
+                let out = if let Some(function) = functions.get(&var) {
+                    let res = &args[0] >= &args[1];
+                    match function.as_str() {
+                        "Identity" => Output::Bool(true), // do nothing
+                        "Geq" => Output::Bool(args[0] >= args[1]),
+                        _ => Output::Float(args[0]),
+                    }
+                } else {
+                    Output::Bool(true)
+                };
+                bindings.push((var, out));
+            }
+        }
+        if let Some((_, value)) = bindings.last() {
+            match value {
+                Output::Bool(x) => *x,
+                _ => false,
+            }
+        } else {
+            true
+        }
+    }
+
+    // TODO parameterize
+    fn fetch(&self, x: &DVector<f32>, arg: Id, ob_index: BTreeMap<Id, usize>) -> f32 {
+        if self.obs.contains(&arg) {
+            x[ob_index[&arg]]
+        } else if self.borrows.keys().contains(&arg) {
+            x[ob_index[&self.borrows[&arg]]]
+        } else {
+            7f32 // TODO remove
+        }
+    }
+}
+
+// impl<Id, Var, Exp> From<Vec<(ComputeGraph<Id>, PolynomialSystem<Var, f32, Exp>)>>
+//     for NumericalPolynomialSwitchingSystem<Id, Exp>
+// where
+//     Var: Ord,
+//     Exp: Ord,
+//     Id: Hash + Clone + Eq + Debug,
+// {
+//     fn from(subsystems: Vec<(ComputeGraph<Id>, PolynomialSystem<Var, Coef, Exp>)>) -> Self {
+//         let subsystems = subsystems
+//             .into_iter()
+//             .map(|(id, poly)| (id, poly.to_numerical()))
+//             .collect::<Vec<_>>();
+//         NumericalPolynomialSwitchingSystem { subsystems }
+//     }
+// }
+
+impl<Id, Exp> ODESystem for NumericalPolynomialSwitchingSystem<Id, Exp>
+where
+    Exp: Clone + Ord + Display + num_traits::One,
+    f32: Pow<Exp, Output = f32>,
+    Id: Eq + Hash + Ord + Clone + Debug,
+{
+    fn vector_field(&self, dx: &mut DVector<f32>, x: &DVector<f32>, _t: f32) {
+        let subsystem = self.clone().get_current(x.clone());
+        for i in 0..dx.len() {
+            dx[i] = match subsystem.components.get(i) {
+                Some(p) => p.eval(|var| x[*var]),
+                None => 0f32,
+            }
+        }
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use expect_test::expect;
@@ -208,4 +383,87 @@ mod tests {
         "#]];
         expected.assert_eq(&textplot_ode_result(&problem, &result));
     }
+
+    #[test]
+    fn water_simulation() {
+        let th = std::rc::Rc::new(crate::stdlib::th_modal_state_aux());
+        let model = crate::stdlib::water(th);
+        let system = crate::stdlib::analyses::ode::PetriNetMassActionFunctionAnalysis::default()
+            .build_switching_system(&model);
+
+        let functions = HashMap::from([(ustr::ustr("comparator"), String::from("Geq"))]);
+        let sys = system.to_numerical(functions);
+
+        // container, lake, sediment, watershed
+        let initial = DVector::from_column_slice(&[6.0, 4.0, 1.0, 4.0]);
+        // let initial = DVector::from_column_slice(&[1.0, 2.0, 4.0]);
+        let problem = ODEProblem::new(sys, initial).end_time(10.0);
+        let result = problem.solve_rk4(0.1).unwrap();
+        let expected = expect![[r#"
+            ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠒⢣⠀⣀⠤⠤⠒⠒⠉⠒⠒⠒⠉⠉⠒⠒⠒⠒⠉⠉⠉⠉⠉⠒⠒⠒⠒⠒⠒⠂ 4.5
+            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠉⠀⠀⠀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⣱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⢇⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠂⠘⡄⠀⠀⠀⠀⠀⡔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡁⠀⠱⡀⠀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⠀⠀⢣⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠂⠀⠀⠀⢇⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⡀
+            ⡁⠀⠀⠀⡎⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⣀⠤⠤⠤⠤⠤⠒⠒⠒⠒⠒⠒⠒⠒⠊⠀⠀⠀⠀⠀⠀⠀
+            ⠄⠀⠀⡜⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠒⠒⢲⠓⠒⠒⠚⢖⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡁⢀⠇⠀⠀⠀⠀⠀⠣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⡜⠀⠀⠀⠀⠀⠀⠀⠑⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⢲⠁⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠒⠤⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠒⠒⠒⠤⠤⠤⠤⣀⣀⣀⣀⣀⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠉⠉⠉⠁ 0.0
+            0.0                                            5.0
+        "#]];
+        expected.assert_eq(&textplot_ode_result(&problem, &result));
+    }
+
+    #[test]
+    fn water_no_logic_simulation() {
+        let th = std::rc::Rc::new(crate::stdlib::th_modal_state_aux());
+        let model = crate::stdlib::water_no_logic(th);
+        let system = crate::stdlib::analyses::ode::PetriNetMassActionFunctionAnalysis::default()
+            .build_switching_system(&model);
+
+        let functions = HashMap::from([(ustr::ustr("comparator"), String::from("Geq"))]);
+        let sys = system.to_numerical(functions);
+
+        // sediment, watershed, lake
+        let initial = DVector::from_column_slice(&[1.0, 2.0, 4.0]);
+        let problem = ODEProblem::new(sys, initial).end_time(10.0);
+        let result = problem.solve_rk4(0.1).unwrap();
+        let expected = expect![[r#"
+            ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠒⢣⠀⣀⠤⠤⠒⠒⠉⠒⠒⠒⠉⠉⠒⠒⠒⠒⠉⠉⠉⠉⠉⠒⠒⠒⠒⠒⠒⠂ 4.5
+            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠉⠀⠀⠀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⣱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⢇⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠂⠘⡄⠀⠀⠀⠀⠀⡔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡁⠀⠱⡀⠀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⠀⠀⢣⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠂⠀⠀⠀⢇⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⡀
+            ⡁⠀⠀⠀⡎⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⣀⠤⠤⠤⠤⠤⠒⠒⠒⠒⠒⠒⠒⠒⠊⠀⠀⠀⠀⠀⠀⠀
+            ⠄⠀⠀⡜⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠒⠒⢲⠓⠒⠒⠚⢖⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡁⢀⠇⠀⠀⠀⠀⠀⠣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⡜⠀⠀⠀⠀⠀⠀⠀⠑⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⢲⠁⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠒⠤⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠒⠒⠒⠤⠤⠤⠤⣀⣀⣀⣀⣀⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀
+            ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠉⠉⠉⠁ 0.0
+            0.0                                            5.0
+        "#]];
+        expected.assert_eq(&textplot_ode_result(&problem, &result));
+    }
 }
diff --git a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs
index bc58f28d1..a95869773 100644
--- a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs
+++ b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs
@@ -4,13 +4,15 @@ Such ODEs are based on the *law of mass action* familiar from chemistry and
 mathematical epidemiology.
  */
 
-use std::collections::{BTreeMap, HashMap};
-use std::fmt::Debug;
+use std::collections::{BTreeMap, HashMap, HashSet};
+use std::fmt::{Debug, Display};
 use std::hash::Hash;
 
+use itertools::Itertools;
+
 use nalgebra::DVector;
 use num_traits::Zero;
-use ustr::{Ustr, ustr};
+use ustr::{ustr, Ustr};
 
 #[cfg(feature = "serde")]
 use serde::{Deserialize, Serialize};
@@ -19,11 +21,13 @@ use tsify::Tsify;
 
 use super::ODEAnalysis;
 use crate::dbl::{
-    model::{DiscreteTabModel, FgDblModel, ModalDblModel, MutDblModel, TabEdge},
-    theory::{ModalMorType, ModalObType, TabMorType, TabObType},
+    model::{DiscreteTabModel, FgDblModel, ModalDblModel, ModalOb, MutDblModel, TabEdge},
+    theory::{ModalMorType, ModalObType, ModeApp, TabMorType, TabObType},
+};
+use crate::one::{FgCategory, FinGraph, HashGraph};
+use crate::simulate::ode::{
+    NumericalPolynomialSwitchingSystem, NumericalPolynomialSystem, ODEProblem, PolynomialSystem,
 };
-use crate::one::FgCategory;
-use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem};
 use crate::zero::{alg::Polynomial, rig::Monomial};
 
 /// Data defining a mass-action ODE problem for a model.
@@ -148,6 +152,7 @@ impl StockFlowMassActionAnalysis {
         &self,
         model: &DiscreteTabModel<Id, Ustr>,
     ) -> PolynomialSystem<Id, Parameter<Id>, u8> {
+        // associate each flow to its domain.
         let mut terms: HashMap<Id, Monomial<Id, u8>> = model
             .mor_generators_with_type(&self.flow_mor_type)
             .map(|flow| {
@@ -223,6 +228,505 @@ fn into_numerical_system<Id: Eq + Clone + Hash + Ord>(
     ODEAnalysis::new(problem, ob_index)
 }
 
+// ------------------------------------------------------------------------- //
+
+// trait Evaluable<T> {
+//     type Output;
+//     fn eval(&self) -> Self::Output;
+// }
+
+// pub struct GeqFunction<T>(T, T);
+// pub struct IdentityFunction<T>(T);
+
+// impl<T: PartialOrd> Evaluable<T> for GeqFunction<T> {
+//     type Output = bool;
+
+//     fn eval(&self) -> Self::Output {
+//         self.0 >= self.1
+//     }
+// }
+
+// impl<T> Evaluable<T> for IdentityFunction<T>
+// where
+//     T: Clone,
+// {
+//     type Output = T;
+
+//     fn eval(&self) -> Self::Output {
+//         self.0.clone()
+//     }
+// }
+
+// pub enum EligibleFunctions<Y> {
+//     Geq(Y),
+//     Identity(Y),
+// }
+
+// impl<T> EligibleFunctions<Y> where T: PartialOrd {
+//     pub fn eval(&self, args: Vec<T>) -> Y {
+
+#[derive(Debug)]
+pub enum Output {
+    Bool(bool),
+    Float(f32),
+}
+
+#[derive(Debug, Clone)]
+#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
+#[cfg_attr(
+    feature = "serde-wasm",
+    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
+)]
+pub enum EligibleFunctions {
+    Geq(),
+    Identity(),
+}
+
+impl EligibleFunctions {
+    pub fn eval(self, args: &Vec<f32>) -> Output {
+        match self {
+            EligibleFunctions::Geq() => Output::Bool(args[0] >= args[1]),
+            EligibleFunctions::Identity() => Output::Float(args[0]),
+        }
+    }
+}
+
+/// Data defining a mass-action ODE problem for a model with function parameters
+#[derive(Clone)]
+#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
+#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
+#[cfg_attr(
+    feature = "serde-wasm",
+    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
+)]
+pub struct SwitchingMassActionProblemData<Id>
+where
+    Id: Eq + Hash + Clone + Debug,
+{
+    /// mass action problem data
+    pub mass: MassActionProblemData<Id>,
+
+    /// functions associated to T(aux) --> aux
+    pub functions: HashMap<Id, String>,
+}
+
+#[derive(Debug, Clone)]
+pub struct ComputeGraph<Id>
+where
+    Id: Clone + Eq + Hash + Debug,
+{
+    /// compute graph
+    pub graph: HashGraph<Id, String>, // TODO parameterize edge type
+
+    pub obs: Vec<Id>,
+
+    pub funcs: HashMap<Id, Vec<Id>>,
+
+    pub auxes: HashMap<Id, f32>,
+
+    pub borrows: HashMap<Id, Id>,
+
+    // /// in-ports
+
+    // /// out-ports
+    // inp
+    ///
+    pub toposort: Vec<Id>,
+}
+
+impl<Id: Eq + Hash + Clone + Debug> ComputeGraph<Id> {
+    pub fn new() -> Self {
+        Self {
+            graph: HashGraph::<Id, String>::default(),
+            obs: Vec::<Id>::new(),
+            auxes: HashMap::<Id, f32>::new(),
+            funcs: HashMap::<Id, Vec<Id>>::new(), // store inputs
+            borrows: HashMap::<Id, Id>::new(),
+            toposort: vec![],
+        }
+    }
+
+    // TODO need to verify that arrow is an arrow. Otherwise, return an identity function
+    // XXX update this
+    /// This extracts a compute graph associated to an arrow in a model.
+    pub fn complete(model: &ModalDblModel<Id, Ustr>, arrow: Id) -> Self {
+        let mut cg = ComputeGraph::<Id>::new();
+        match model.mor_generator_dom(&arrow.clone()) {
+            ModalOb::Generator(id) => {
+                cg.add_vertex(arrow.clone());
+                cg.add_vertex(id.clone());
+                cg.connect(id, arrow.clone());
+            }
+            ModalOb::List(_, xs) => {
+                cg.add_vertex(arrow.clone());
+                for x in xs {
+                    // TODO `x` may not be Generator
+                    cg.add_vertex(x.clone().unwrap_generator());
+                    cg.connect(x.clone().unwrap_generator(), arrow.clone());
+                }
+            }
+            _ => todo!(),
+        }
+        match model.mor_generator_cod(&arrow) {
+            ModalOb::Generator(id) => {
+                cg.add_vertex(arrow.clone());
+                cg.add_vertex(id.clone());
+                cg.connect(arrow.clone(), id.clone());
+            }
+            ModalOb::List(_, xs) => {
+                cg.add_vertex(arrow.clone());
+                for x in xs {
+                    // TODO `x` may not be Generator
+                    cg.add_vertex(x.clone().unwrap_generator());
+                    cg.connect(arrow.clone(), x.unwrap_generator());
+                }
+            }
+            _ => todo!(),
+        }
+
+        cg.toposort = crate::one::graph_algorithms::toposort(&cg.graph).expect("!");
+        cg
+    }
+
+    pub fn add_vertex(&mut self, v: Id) -> bool {
+        self.graph.add_vertex(v)
+    }
+
+    pub fn add_edge(&mut self, e: String, dom: Id, cod: Id) -> bool {
+        self.graph.add_edge(e, dom, cod)
+    }
+
+    // TODO check that the name exists?
+    /// Connects a src and tgt by an edge whose name is generated automatically.
+    pub fn connect(&mut self, dom: Id, cod: Id) -> bool {
+        self.graph.add_edge(format!("{:?}=>{:?}", dom.clone(), cod.clone()), dom, cod)
+    }
+}
+
+#[derive(Clone)]
+pub struct SwitchingSystem<Id>
+where
+    Id: Eq + Clone + Hash + Debug,
+{
+    pub null_model: PolynomialSystem<Id, Parameter<Id>, u8>,
+    pub subsystems: Vec<(ComputeGraph<Id>, PolynomialSystem<Id, Parameter<Id>, u8>)>,
+    pub ob_index: BTreeMap<Id, usize>,
+}
+
+impl<Id> SwitchingSystem<Id>
+where
+    Id: Hash + Clone + Eq + Ord + Debug,
+{
+    pub fn to_numerical(
+        &self,
+        functions: HashMap<Id, String>,
+    ) -> NumericalPolynomialSwitchingSystem<Id, u8> {
+        let subsystems: Vec<(ComputeGraph<Id>, _)> = self
+            .subsystems
+            .clone()
+            .into_iter()
+            .map(|(id, poly)| (id, poly.extend_scalars(|p| p.eval(|_| 1.0)).to_numerical()))
+            .collect();
+        // TODO option is used here to avoid implementing Default.
+        let null_model =
+            Some(self.null_model.clone().extend_scalars(|p| p.eval(|_| 1.0)).to_numerical());
+        NumericalPolynomialSwitchingSystem {
+            ob_index: self.ob_index.clone(),
+            functions,
+            null_model,
+            subsystems,
+        }
+    }
+}
+
+/// Convenience struct
+#[derive(Clone, Debug)]
+pub struct PetriNetSystemData<Id>
+where
+    Id: Eq + Clone + Hash + Debug,
+{
+    pub flowneg: HashMap<Id, Id>,
+    pub flowpos: HashMap<Id, Id>,
+    pub outpos_dom: HashSet<Id>,
+    pub outneg_dom: HashSet<Id>,
+    pub mediators: Vec<Id>,
+}
+
+impl<Id> PetriNetSystemData<Id>
+where
+    Id: Hash + Eq + Clone + Debug,
+{
+    fn make(model: &ModalDblModel<Id, Ustr>, pna: &PetriNetMassActionFunctionAnalysis) -> Self {
+        let mut flowneg = HashMap::new();
+        let mut flowpos = HashMap::new();
+        let outpos_dom: HashSet<_> = HashSet::from_iter::<_>(
+            model
+                .mor_generators_with_type(&pna.outpos_mor_type)
+                .map(|p| {
+                    let dom = model.mor_generator_dom(&p).unwrap_generator();
+                    flowpos.insert(dom.clone(), model.mor_generator_cod(&p).unwrap_generator());
+                    dom
+                })
+                .collect::<Vec<_>>(),
+        );
+        let outneg_dom: HashSet<_> = HashSet::from_iter::<_>(
+            model
+                .mor_generators_with_type(&pna.outneg_mor_type)
+                .map(|p| {
+                    let dom = model.mor_generator_dom(&p).unwrap_generator();
+                    flowneg.insert(dom.clone(), model.mor_generator_cod(&p).unwrap_generator());
+                    dom
+                })
+                .collect::<Vec<_>>(),
+        );
+        let mediators: Vec<_> = outpos_dom.intersection(&outneg_dom).cloned().collect();
+
+        Self {
+            flowneg,
+            flowpos,
+            outpos_dom,
+            outneg_dom,
+            mediators,
+        }
+    }
+}
+
+// TODO rename
+/**
+ */
+pub struct PetriNetMassActionFunctionAnalysis {
+    /// Object type for states
+    pub state_ob_type: ModalObType<Ustr>,
+    /// Object type for auxiliary variables
+    pub aux_ob_type: ModalObType<Ustr>,
+    /// Morphism type for functions between auxiliary variables
+    pub fun_mor_type: ModalMorType<Ustr>,
+    /// Morphism type for "borrowing," or maps from states to auxiliary variables
+    pub borrow_mor_type: ModalMorType<Ustr>,
+    /// Morphism type for positively-signed flows
+    pub outpos_mor_type: ModalMorType<Ustr>,
+    /// Morphism type for negatively-signed flows
+    pub outneg_mor_type: ModalMorType<Ustr>,
+}
+
+impl Default for PetriNetMassActionFunctionAnalysis {
+    fn default() -> Self {
+        Self {
+            state_ob_type: ModalObType::new(ustr("State")),
+            aux_ob_type: ModalObType::new(ustr("Auxiliary")),
+            fun_mor_type: ModalMorType::One(ModeApp::new(ustr("function"))),
+            borrow_mor_type: ModalMorType::One(ModeApp::new(ustr("borrow"))),
+            outpos_mor_type: ModalMorType::One(ModeApp::new(ustr("out-pos"))),
+            outneg_mor_type: ModalMorType::One(ModeApp::new(ustr("out-neg"))),
+        }
+    }
+}
+
+impl PetriNetMassActionFunctionAnalysis {
+    fn build_graph<Id: Clone + Eq + Hash + Debug>(
+        &self,
+        model: &ModalDblModel<Id, Ustr>,
+    ) -> ComputeGraph<Id> {
+        let mut cg = ComputeGraph::<Id>::new();
+        cg.obs = model.ob_generators_with_type(&self.state_ob_type).collect::<Vec<Id>>();
+        // TODO these should just be auxes entering the graph
+        // cg.auxes = HashMap::from_iter(
+        //     model.ob_generators_with_type(&self.aux_ob_type).map(|aux| (aux, 4.5f32)),
+        // );
+        cg.borrows =
+            HashMap::from_iter(model.mor_generators_with_type(&self.borrow_mor_type).map(|f| {
+                (
+                    model.mor_generator_cod(&f.clone()).unwrap_generator(),
+                    model.mor_generator_dom(&f.clone()).unwrap_generator(),
+                )
+            }));
+        let arrows = model.mor_generators_with_type(&self.fun_mor_type).collect::<Vec<_>>();
+        for arrow in arrows {
+            let mut inputs: Vec<Id> = vec![];
+            match model.mor_generator_dom(&arrow) {
+                ModalOb::Generator(id) => {
+                    inputs.push(id.clone());
+                    cg.add_vertex(arrow.clone());
+                    cg.add_vertex(id.clone());
+                    cg.connect(id, arrow.clone());
+                }
+                ModalOb::List(_, xs) => {
+                    cg.add_vertex(arrow.clone());
+                    for x in xs {
+                        inputs.push(x.clone().unwrap_generator());
+                        // TODO `x` may not be Generator
+                        cg.add_vertex(x.clone().unwrap_generator());
+                        cg.connect(x.clone().unwrap_generator(), arrow.clone());
+                    }
+                }
+                _ => todo!(),
+            }
+            cg.funcs.insert(arrow.clone(), inputs);
+            //
+            match model.mor_generator_cod(&arrow) {
+                ModalOb::Generator(id) => {
+                    cg.add_vertex(arrow.clone());
+                    cg.add_vertex(id.clone());
+                    cg.connect(arrow.clone(), id.clone());
+                }
+                ModalOb::List(_, xs) => {
+                    cg.add_vertex(arrow.clone());
+                    for x in xs {
+                        // TODO `x` may not be Generator
+                        cg.add_vertex(x.clone().unwrap_generator());
+                        cg.connect(arrow.clone(), x.unwrap_generator());
+                    }
+                }
+                _ => todo!(),
+            }
+        }
+
+        cg.toposort = crate::one::graph_algorithms::toposort(&cg.graph).expect("!");
+        cg
+    }
+
+    pub fn build_system<Id>(
+        &self,
+        model: &ModalDblModel<Id, Ustr>,
+        model_data: PetriNetSystemData<Id>, // TODO rename
+        filter_ids: Vec<Id>,
+    ) -> PolynomialSystem<Id, Parameter<Id>, u8>
+    where
+        Id: Eq + Clone + Hash + Debug + Ord + Display,
+    {
+        // TODO what are these
+        let terms: Vec<(Id, Polynomial<Id, Parameter<Id>, u8>)> = model_data
+            .mediators
+            .iter()
+            // keep the mediator if it is in the list
+            .filter(|m| filter_ids.clone().contains(m))
+            .map(|mediator| {
+                let param = Parameter::generator(mediator.clone());
+                let term = Monomial::generator(model_data.flowneg[mediator].clone());
+                (mediator.clone(), [(param, term)].into_iter().collect())
+            })
+            .collect();
+
+        let mut sys = PolynomialSystem::new();
+        for ob in model.ob_generators_with_type(&self.state_ob_type) {
+            sys.add_term(ob, Polynomial::zero());
+        }
+        // TODO
+        for (mediator, term) in terms.iter() {
+            sys.add_term(model_data.flowneg[&mediator.clone()].clone(), -term.clone());
+            sys.add_term(model_data.flowpos[&mediator].clone(), term.clone());
+        }
+
+        for (var, component) in sys.clone().components.iter() {
+            println!("d{var} = {component}")
+        }
+        println!("--------------------------");
+
+        sys
+    }
+
+    pub fn build_switching_system<Id>(&self, model: &ModalDblModel<Id, Ustr>) -> SwitchingSystem<Id>
+    where
+        Id: Clone + Eq + Hash + Debug + Ord + Display,
+    {
+        let modeldata = PetriNetSystemData::make(model.into(), &self);
+
+        // TODO ensure they are connected components
+        let programs = self.build_graph(model);
+
+        // build subsystem with no programs.
+        // This means we exclude all flows are are mediated by
+        // any program
+        let affected_mediators = modeldata
+            .clone()
+            .mediators
+            .into_iter()
+            .filter(|m| programs.graph.vertices().contains(m));
+        println!("INIT MODEL");
+        let init_model = self.build_system(
+            &model,
+            modeldata.clone(),
+            modeldata
+                .clone()
+                .mediators
+                .into_iter()
+                .filter(|m| !programs.graph.vertices().contains(m))
+                .collect::<Vec<_>>(),
+        );
+        // println!("DISCRETE MODEL");
+        let null_model = self.build_system(&model, modeldata.clone(), vec![]);
+        let mut subsystems = affected_mediators
+            .map(|m| {
+                println!("MEDIATOR: {:?}", &m);
+                let sys = self.build_system(&model, modeldata.clone(), vec![m]);
+                (programs.clone(), sys)
+            })
+            .collect::<Vec<(ComputeGraph<Id>, PolynomialSystem<Id, Parameter<Id>, u8>)>>();
+        // subsystems.push((ComputeGraph::<Id>::new(), init_model));
+        let ob_index: BTreeMap<Id, usize> =
+            null_model.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect();
+        let null_model = init_model.clone();
+        SwitchingSystem {
+            null_model,
+            subsystems,
+            ob_index,
+        }
+    }
+
+    pub fn build_numerical_system<Id: Eq + Clone + Hash + Ord + Debug + std::fmt::Display>(
+        &self,
+        model: &ModalDblModel<Id, Ustr>,
+        data: SwitchingMassActionProblemData<Id>,
+    ) -> ODEAnalysis<Id, NumericalPolynomialSwitchingSystem<Id, u8>> {
+        into_numerical_switching_system(self.build_switching_system(model), data)
+    }
+}
+
+fn into_numerical_switching_system<Id: Eq + Clone + Hash + Ord + Debug>(
+    system: SwitchingSystem<Id>,
+    data: SwitchingMassActionProblemData<Id>,
+) -> ODEAnalysis<Id, NumericalPolynomialSwitchingSystem<Id, u8>> {
+    let ob_index: BTreeMap<_, _> = system.ob_index.clone();
+    let n = ob_index.len();
+    let null_model = Some(
+        system
+            .null_model
+            .clone()
+            .extend_scalars(|p| {
+                p.eval(|flow| data.mass.rates.get(flow).copied().unwrap_or_default())
+            })
+            .to_numerical(),
+    );
+    let initial_values = ob_index
+        .keys()
+        .map(|ob| data.mass.initial_values.get(ob).copied().unwrap_or_default());
+    let x0 = DVector::from_iterator(n, initial_values);
+    let subsystems = system
+        .subsystems
+        .into_iter()
+        .map(|(graph, sys)| {
+            let sys = sys
+                .extend_scalars(|poly| {
+                    poly.eval(|flow| data.mass.rates.get(flow).copied().unwrap_or_default())
+                })
+                .to_numerical();
+            (graph, sys)
+        })
+        .collect::<Vec<_>>();
+    let sys = NumericalPolynomialSwitchingSystem {
+        ob_index: ob_index.clone(),
+        subsystems,
+        null_model,
+        functions: data.functions,
+    };
+    let problem = ODEProblem::new(sys, x0).end_time(data.mass.duration);
+    ODEAnalysis {
+        problem,
+        variable_index: ob_index,
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use expect_test::expect;
@@ -255,4 +759,13 @@ mod tests {
         "#]);
         expected.assert_eq(&sys.to_string());
     }
+
+    #[test]
+    fn water_bath() {
+        let th = Rc::new(th_modal_state_aux());
+        let model = water(th);
+        let sys = PetriNetMassActionFunctionAnalysis::default().build_switching_system(&model);
+        // TODO
+        assert!(true);
+    }
 }
diff --git a/packages/catlog/src/stdlib/analyses/ode/mod.rs b/packages/catlog/src/stdlib/analyses/ode/mod.rs
index 63adba85a..e1ef038bd 100644
--- a/packages/catlog/src/stdlib/analyses/ode/mod.rs
+++ b/packages/catlog/src/stdlib/analyses/ode/mod.rs
@@ -7,6 +7,8 @@ use derivative::Derivative;
 use derive_more::Constructor;
 use ode_solvers::dop_shared::IntegrationError;
 
+use nalgebra::DVector;
+
 #[cfg(feature = "serde")]
 use serde::{Deserialize, Serialize};
 #[cfg(feature = "serde-wasm")]
diff --git a/packages/catlog/src/stdlib/models.rs b/packages/catlog/src/stdlib/models.rs
index e0fc91bac..3daa4abed 100644
--- a/packages/catlog/src/stdlib/models.rs
+++ b/packages/catlog/src/stdlib/models.rs
@@ -1,7 +1,7 @@
 //! Standard library of models of double theories.
 
 use std::rc::Rc;
-use ustr::{Ustr, ustr};
+use ustr::{ustr, Ustr};
 
 use crate::dbl::{model::*, theory::*};
 use crate::one::Path;
@@ -141,6 +141,117 @@ pub fn catalyzed_reaction(th: Rc<UstrModalDblTheory>) -> UstrModalDblModel {
     model
 }
 
+/**
+ */
+pub fn water(th: Rc<UstrModalDblTheory>) -> UstrModalDblModel {
+    let (state_type, aux_type) =
+        (ModalObType::new(ustr("State")), ModalObType::new(ustr("Auxiliary")));
+    let mut model = UstrModalDblModel::new(th);
+    let (watershed, lake, sediment) = (ustr("watershed"), ustr("lake"), ustr("sediment"));
+    model.add_ob(watershed, state_type.clone());
+    model.add_ob(lake, state_type.clone());
+    model.add_ob(sediment, state_type.clone());
+    let blake = ustr("&lake");
+    model.add_ob(blake, aux_type.clone());
+    let borrow = ModalMorType::One(ModeApp::new(ustr("borrowing")));
+    model.add_mor(
+        ustr("borrow-lake"),
+        ModalOb::Generator(lake.clone()),
+        ModalOb::Generator(blake.clone()),
+        borrow.clone(),
+    );
+    let container = ustr("container");
+    model.add_ob(container, state_type.clone());
+    let bcontainer = ustr("&container");
+    model.add_ob(bcontainer, aux_type.clone());
+    model.add_mor(
+        ustr("borrow-container"),
+        ModalOb::Generator(container.clone()),
+        ModalOb::Generator(bcontainer.clone()),
+        borrow.clone(),
+    );
+    let (comparator, lake_sediment) = (ustr("comparator"), ustr("lake_sediment"));
+    model.add_mor(
+        comparator,
+        ModalOb::List(List::Plain, vec![blake.into(), bcontainer.into()]).into(),
+        ModalOb::Generator(lake_sediment.clone()),
+        ModalMorType::One(ModeApp::new(ustr("function"))),
+    );
+    // watershed outflow
+    let watershed_lake = ustr("watershed_lake");
+    model.add_ob(watershed_lake.clone(), aux_type.clone());
+    model.add_mor(
+        ustr("watershed-lake-outflow"),
+        ModalOb::Generator(watershed_lake.clone()),
+        ModalOb::Generator(watershed.clone()),
+        ModalMorType::One(ModeApp::new(ustr("out-neg"))),
+    );
+    model.add_mor(
+        ustr("watershed-lake-inflow"),
+        ModalOb::Generator(watershed_lake.clone()),
+        ModalOb::Generator(lake.clone()),
+        ModalMorType::One(ModeApp::new(ustr("out-pos"))),
+    );
+    //
+    model.add_mor(
+        ustr("lake-sediment-outflow"),
+        ModalOb::Generator(lake_sediment.clone()),
+        ModalOb::Generator(lake.clone()),
+        ModalMorType::One(ModeApp::new(ustr("out-neg"))),
+    );
+    model.add_mor(
+        ustr("lake-sediment-infow"),
+        ModalOb::Generator(lake_sediment.clone()),
+        ModalOb::Generator(sediment.clone()),
+        ModalMorType::One(ModeApp::new(ustr("out-pos"))),
+    );
+    model
+}
+
+/**
+ */
+pub fn water_no_logic(th: Rc<UstrModalDblTheory>) -> UstrModalDblModel {
+    let (state_type, aux_type) =
+        (ModalObType::new(ustr("State")), ModalObType::new(ustr("Auxiliary")));
+    let mut model = UstrModalDblModel::new(th);
+    let (watershed, lake, sediment) = (ustr("watershed"), ustr("lake"), ustr("sediment"));
+    model.add_ob(watershed, state_type.clone());
+    model.add_ob(lake, state_type.clone());
+    model.add_ob(sediment, state_type.clone());
+    let container = ustr("container");
+    model.add_ob(container, aux_type.clone());
+    let lake_sediment = ustr("lake_sediment");
+    // watershed outflow
+    let watershed_lake = ustr("watershed_lake");
+    model.add_ob(watershed_lake.clone(), aux_type.clone());
+    model.add_mor(
+        ustr("watershed-lake-outflow"),
+        ModalOb::Generator(watershed_lake.clone()),
+        ModalOb::Generator(watershed.clone()),
+        ModalMorType::One(ModeApp::new(ustr("out-neg"))),
+    );
+    model.add_mor(
+        ustr("watershed-lake-inflow"),
+        ModalOb::Generator(watershed_lake.clone()),
+        ModalOb::Generator(lake.clone()),
+        ModalMorType::One(ModeApp::new(ustr("out-pos"))),
+    );
+    //
+    model.add_mor(
+        ustr("lake-sediment-outflow"),
+        ModalOb::Generator(lake_sediment.clone()),
+        ModalOb::Generator(lake.clone()),
+        ModalMorType::One(ModeApp::new(ustr("out-neg"))),
+    );
+    model.add_mor(
+        ustr("lake-sediment-infow"),
+        ModalOb::Generator(lake_sediment.clone()),
+        ModalOb::Generator(sediment.clone()),
+        ModalMorType::One(ModeApp::new(ustr("out-pos"))),
+    );
+    model
+}
+
 #[cfg(test)]
 mod tests {
     use super::super::theories::*;
@@ -182,4 +293,11 @@ mod tests {
         let th = Rc::new(th_sym_monoidal_category());
         assert!(catalyzed_reaction(th).validate().is_ok());
     }
+
+    #[test]
+    fn modal_state_aux() {
+        let th = Rc::new(th_modal_state_aux());
+        // dbg!(water(th));
+        assert!(true)
+    }
 }
diff --git a/packages/catlog/src/stdlib/theories.rs b/packages/catlog/src/stdlib/theories.rs
index 948281e1b..061e29d34 100644
--- a/packages/catlog/src/stdlib/theories.rs
+++ b/packages/catlog/src/stdlib/theories.rs
@@ -3,7 +3,7 @@
 use ustr::ustr;
 
 use crate::dbl::theory::*;
-use crate::one::{Path, fp_category::UstrFpCategory};
+use crate::one::{fp_category::UstrFpCategory, Path};
 
 /** The empty theory, which has a single model, the empty model.
 
@@ -151,6 +151,11 @@ pub fn th_sym_monoidal_category() -> UstrModalDblTheory {
     th_list_algebra(List::Symmetric)
 }
 
+///
+pub fn th_modal_state_aux() -> UstrModalDblTheory {
+    th_modal_sf(List::Plain)
+}
+
 /** The theory of a strict algebra of a list monad.
 
 This is a modal double theory, parametric over which variant of the double list
@@ -209,6 +214,22 @@ fn th_generalized_multicategory(list: List) -> UstrModalDblTheory {
     th
 }
 
+///
+fn th_modal_sf(list: List) -> UstrModalDblTheory {
+    let m = Modality::List(list);
+    let mut th: UstrModalDblTheory = Default::default();
+    let (state, aux) = (ustr("State"), ustr("Auxiliary"));
+    th.add_ob_type(state);
+    th.add_ob_type(aux);
+    let function = ustr("function");
+    th.add_mor_type(function, ModeApp::new(aux).apply(m), ModeApp::new(aux));
+    let (borrow, outpos, outneg) = (ustr("borrow"), ustr("out-pos"), ustr("out-neg"));
+    th.add_mor_type(borrow, ModeApp::new(state), ModeApp::new(aux));
+    th.add_mor_type(outpos, ModeApp::new(aux), ModeApp::new(state));
+    th.add_mor_type(outneg, ModeApp::new(aux), ModeApp::new(state));
+    th
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -237,6 +258,7 @@ mod tests {
         assert!(th_monoidal_category().validate().is_ok());
         assert!(th_lax_monoidal_category().validate().is_ok());
         assert!(th_multicategory().validate().is_ok());
+        assert!(th_modal_state_aux().validate().is_ok());
     }
 
     #[test]
diff --git a/packages/frontend/src/notebook/notebook_editor.tsx b/packages/frontend/src/notebook/notebook_editor.tsx
index 25bc7577f..c278e5007 100644
--- a/packages/frontend/src/notebook/notebook_editor.tsx
+++ b/packages/frontend/src/notebook/notebook_editor.tsx
@@ -301,6 +301,7 @@ export function NotebookEditor<T>(props: {
                         };
 
                         const cell = props.notebook.cellContents[cellId];
+						console.log(props.notebook.cellContents);
                         invariant(cell, `Failed to find contents for cell '${cellId}'`);
 
                         if (cell.tag !== "rich-text") {
diff --git a/packages/frontend/src/stdlib/analyses/index.ts b/packages/frontend/src/stdlib/analyses/index.ts
index f1f066d95..03e7729fa 100644
--- a/packages/frontend/src/stdlib/analyses/index.ts
+++ b/packages/frontend/src/stdlib/analyses/index.ts
@@ -7,3 +7,4 @@ export * from "./model_graph";
 export * from "./petri_net_visualization";
 export * from "./stock_flow_diagram";
 export * from "./submodel_graphs";
+export * from "./switching_mass_action";
diff --git a/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx b/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx
new file mode 100644
index 000000000..1c6279311
--- /dev/null
+++ b/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx
@@ -0,0 +1,189 @@
+import { createMemo } from "solid-js";
+
+import type {
+    DblModel,
+    ODEResult,
+    SwitchingMassActionModelData,
+    SwitchingMassActionProblemData,
+} from "catlog-wasm";
+import type { ModelAnalysisProps } from "../../analysis";
+import {
+    type ColumnSchema,
+    FixedTableEditor,
+    Foldable,
+    createNumericalColumn,
+} from "../../components";
+import type { MorphismDecl, ObjectDecl } from "../../model";
+import type { ModelAnalysisMeta } from "../../theory";
+import { ODEResultPlot } from "../../visualization";
+import { createModelODEPlot } from "./simulation";
+
+import "./simulation.css";
+
+/** Configuration for a mass-action ODE analysis of a model. */
+export type SwitchingMassActionContent = SwitchingMassActionProblemData<string>;
+
+type Simulator = (model: DblModel, data: SwitchingMassActionModelData) => ODEResult;
+
+/** Configure a mass-action ODE analysis for use with models of a theory. */
+export function configureSwitchingMassAction(options: {
+    id?: string;
+    name?: string;
+    description?: string;
+    help?: string;
+    simulate: Simulator;
+    isState?: (ob: ObjectDecl) => boolean;
+    isTransition?: (mor: MorphismDecl) => boolean;
+}): ModelAnalysisMeta<SwitchingMassActionContent> {
+    const {
+        id = "switching-mass-action",
+        name = "Switching Mass-action dynamics",
+        description = "Simulate a switching system using the law of mass action",
+        help = "mass-action",
+        ...otherOptions
+    } = options;
+    return {
+        id,
+        name,
+        description,
+        help,
+        component: (props) => <SwitchingMassAction title={name} {...otherOptions} {...props} />,
+        initialContent: () => ({
+            mass: { rates: {}, initialValues: {}, duration: 10 },
+            functions: {},
+        }),
+    };
+}
+
+/** Analyze a model using mass-action dynamics. */
+export function SwitchingMassAction(
+    props: ModelAnalysisProps<SwitchingMassActionContent> & {
+        simulate: Simulator;
+        isState?: (ob: ObjectDecl) => boolean;
+        isTransition?: (mor: MorphismDecl) => boolean;
+        title?: string;
+    },
+) {
+	console.log(props.liveModel.formalJudgments());
+    const obDecls = createMemo<ObjectDecl[]>(() => {
+        return props.liveModel
+            .formalJudgments()
+            .filter((jgmt) => jgmt.tag === "object")
+			.filter((ob) => ob.obType.content === "State")
+            .filter((ob) => props.isState?.(ob) ?? true);
+    }, []);
+
+	const auxDecls = createMemo<ObjectDecl[]>(() => {
+        return props.liveModel
+            .formalJudgments()
+            .filter((jgmt) => jgmt.tag === "object")
+			.filter((ob) => ob.obType.content === "Auxiliary");
+    }, []);
+
+    // const morDecls = createMemo<MorphismDecl[]>(() => {
+    //     return props.liveModel
+    //         .formalJudgments()
+    //         .filter((jgmt) => jgmt.tag === "morphism")
+			// // .filter((mor) => mor.morType.content === "out-pos")
+    //         .filter((mor) => props.isTransition?.(mor) ?? true);
+    // }, []);
+
+	const functionDecls = createMemo<MorphismDecl[]>(() => {
+		return props.liveModel.formalJudgments().filter((jgmt) => jgmt.tag === "morphism").filter((mor) => mor.morType.content === "function").filter((mor) => props.isTransition?.(mor) ?? true);
+	}, []);
+
+    const obSchema: ColumnSchema<ObjectDecl>[] = [
+        {
+            contentType: "string",
+            header: true,
+            content: (ob) => ob.name,
+        },
+        createNumericalColumn({
+            name: "Initial value",
+            data: (ob) => props.content.mass.initialValues[ob.id],
+            validate: (_, data) => data >= 0,
+            setData: (ob, data) =>
+                props.changeContent((content) => {
+                    content.mass.initialValues[ob.id] = data;
+                }),
+        }),
+    ];
+
+    const auxSchema: ColumnSchema<ObjectDecl>[] = [
+        {
+            contentType: "string",
+            header: true,
+            content: (mor) => mor.name,
+        },
+        createNumericalColumn({
+            name: "Rate",
+            data: (mor) => props.content.mass.rates[mor.id],
+            default: 1,
+            validate: (_, data) => data >= 0,
+            setData: (mor, data) =>
+                props.changeContent((content) => {
+					console.log(content.mass.rates);
+                    content.mass.rates[mor.id] = data;
+                }),
+        }),
+    ];
+
+    // TODO Duration is a simulation-level quantity and should be moved outside of "mass"
+    const toplevelSchema: ColumnSchema<null>[] = [
+        createNumericalColumn({
+            name: "Duration",
+            data: (_) => props.content.mass.duration,
+            validate: (_, data) => data >= 0,
+            setData: (_, data) =>
+                props.changeContent((content) => {
+                    content.mass.duration = data;
+                }),
+        }),
+    ];
+
+
+
+	const functionSchema: ColumnSchema<MorphismDecl>[] = [
+        {
+            contentType: "string",
+            header: true,
+            content: (ob) => ob.name,
+        },
+        {
+            contentType: "enum",
+            name: "Function",
+            variants() {
+                return ["Identity", "Geq"];
+            },
+            content: () => "Geq",
+            setContent: (ob, value) =>
+                props.changeContent((content) => {
+                    if (value === null) {
+                        delete content.functions[ob.id];
+                    } else {
+                        content.functions[ob.id] = value;
+                    }
+                }),
+        },
+        // TODO createEnumColumn
+    ];
+
+    const plotResult = createModelODEPlot(
+        () => props.liveModel,
+        (model: DblModel) => props.simulate(model, props.content),
+    );
+
+    return (
+        <div class="simulation">
+            <Foldable title={props.title}>
+                <div class="parameters">
+                    <FixedTableEditor rows={obDecls()} schema={obSchema} />
+                    <FixedTableEditor rows={auxDecls()} schema={auxSchema} />
+					<FixedTableEditor rows={functionDecls()} schema={functionSchema} />
+                    <FixedTableEditor rows={[null]} schema={toplevelSchema} />
+                </div>
+            </Foldable>
+            <ODEResultPlot result={plotResult()} />
+        </div>
+    );
+}
diff --git a/packages/frontend/src/stdlib/theories.ts b/packages/frontend/src/stdlib/theories.ts
index 216246196..a495d7de2 100644
--- a/packages/frontend/src/stdlib/theories.ts
+++ b/packages/frontend/src/stdlib/theories.ts
@@ -84,6 +84,17 @@ stdTheories.add(
     async () => (await import("./theories/primitive-stock-flow")).default,
 );
 
+stdTheories.add(
+    {
+        id: "stock-flow-with-auxes",
+        name: "Stock and flow with logic",
+        description: "Model accumulation (stocks) and change (flows)",
+        group: "System Dynamics",
+    },
+    async () => (await import("./theories/stock-flow-aux")).default,
+);
+
+
 stdTheories.add(
     {
         id: "reg-net",
diff --git a/packages/frontend/src/stdlib/theories/stock-flow-aux.ts b/packages/frontend/src/stdlib/theories/stock-flow-aux.ts
new file mode 100644
index 000000000..89846ba1d
--- /dev/null
+++ b/packages/frontend/src/stdlib/theories/stock-flow-aux.ts
@@ -0,0 +1,86 @@
+import { ThModalStateAuxCategory } from "catlog-wasm";
+
+import { Theory } from "../../theory";
+import * as analyses from "../analyses";
+import type { TheoryMeta } from "../types";
+
+import styles from "../styles.module.css";
+import svgStyles from "../svg_styles.module.css";
+
+export default function createStockFlowAuxTheory(theoryMeta: TheoryMeta): Theory {
+    const thCategoryStockFlowAux = new ThModalStateAuxCategory();
+
+    return new Theory({
+        ...theoryMeta,
+        theory: thCategoryStockFlowAux.theory(),
+        onlyFreeModels: true,
+        modelTypes: [
+            {
+                tag: "ObType",
+                obType: { tag: "Basic", content: "State" },
+                name: "Stock",
+                description: "Thing with an amount",
+                shortcut: ["S"],
+                cssClasses: [styles.box],
+                svgClasses: [svgStyles.box],
+            },
+            {
+                tag: "MorType",
+				morType: { tag: "Basic", content: "out-neg" },
+                name: "Outflow",
+                description: "Flow from one stock to another",
+                shortcut: ["+"],
+                arrowStyle: "double",
+            },
+			{
+                tag: "MorType",
+				morType: { tag: "Basic", content: "out-pos" },
+                name: "Inflow",
+                description: "Flow from one stock to another",
+                shortcut: ["-"],
+                arrowStyle: "double",
+            },
+			{
+                tag: "ObType",
+                obType: { tag: "Basic", content: "Auxiliary" },
+                name: "Auxiliary",
+                description: "Variable used in computation",
+                shortcut: ["A"],
+                // cssClasses: [styles.box],
+                // svgClasses: [svgStyles.box],
+            },
+            {
+                tag: "MorType",
+                morType: { tag: "Basic", content: "borrow" },
+                name: "Borrow",
+                description: "Influence of a stock on a flow",
+                preferUnnamed: true,
+                shortcut: ["&"],
+            },
+			{
+				tag: "MorType",
+				morType: { tag: "Basic", content: "function" },
+				name: "Function",
+				description: "Function consuming auxiliary variables",
+				preferUnnamed: true,
+				shortcut: ["f"],
+			}
+        ],
+        modelAnalyses: [
+            analyses.configureStockFlowDiagram({
+                id: "diagram",
+                name: "Visualization",
+                description: "Visualize the stock and flow diagram",
+                help: "visualization",
+            }),
+            analyses.configureSwitchingMassAction({
+                simulate(model, data) {
+					return thCategoryStockFlowAux.massAction(model, data);
+                },
+                isTransition(mor) {
+                    return mor.morType.tag === "Basic";
+                },
+            }),
+        ],
+    });
+}