(* ::Package:: *) (* :Title: CompartmentalModeling.wl *) (* :Author: Robert B. Nachbar *) (* :Summary: Functions for constructing, visualizing, and simulating compartmental models. Application areas supported: epidemiology, population ecology, pharmacokinetics, chemical kinetics, systems biology. *) (* :Copyright: (c) 2020, 2021 by Robert B. Nachbar, Jr. *) (* The MIT License Copyright 2021 Robert B. Nachbar, Jr. Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. *) (* :Package Version: 1.7 *) (* :Mathematica Version: 12.1, 12.2, 12.3 *) (* :History: 1.0 14-Mar-2020 1.1 08-Oct-2020 1.2 28-Dec-2020 1.3 15-Feb-2021 1.4 15-Apr-2021 1.5 26-Apr-2021 1.6 05-Jun-2021 1.7 08-Jun-2021 *) (* :Caveats: *) (* :ToDo: o Transition with Missing[] rate x input/output o kinetics o boundary conditions; e.g., solvent o Verbatim wrapper on rate means it is used explicitly & coupled with stoichiometry for ODEs, otherwise mass action kinetics assumed o ComprtmentalModelGraph x option to show compartments only o additional compartments and arrows for catalysts and inhibitors *) (* ::Section:: *) (*BeginPackage*) BeginPackage["CompartmentalModeling`"] ClearAll["`*"] (* ::Section::Closed:: *) (*Usage messages*) (* usage messages *) NullCompartment::usage = "NullCompartment or \[EmptySet] is a symbol that \ represents a source or a sink in a Transition."; (*\[EmptySet]; *) Transition::usage = "Transition[c1, r, c2] specifies a Transition from compartment \ c1 to compartment c2 with transfer rate constant r. In the notebook FrontEnd, a \ Transition is rendered as c1 \!\(\*OverscriptBox[\(\[RightArrow]\), \(r\)]\) c2. The same \ form can be used as input where the arrow is \\[RightArrow]. Multiple \ compartments on either side of the arrow and/or stoichiometry are specified \ with Plus and Times, respectively. Transition[NullCompartment, r, c2] specifies a source (boundary or birth) to compartment c2 with transfer \ rate constant r. In the notebook FrontEnd, a source \ Transition is rendered as \!\(\*OverscriptBox[\(\[RightTeeArrow]\), \(r\)]\) c2. The same \ form can be used as input where the arrow is \\[RightTeeArrow]. Transition[c1, r, NullCompartment] specifies a sink (boundary or death) from compartment c1 with transfer \ rate constant r. In the notebook FrontEnd, a sink \ Transition is rendered as c1 \!\(\*OverscriptBox[\(\[RightArrowBar]\), \(r\)]\). The same \ form can be used as input where the arrow is \\[RightArrowBar]. Transition[c1, c2] is equivalent to Transition[c1, Missing[\"NotSpecified\"], c2]."; (*Catalysts::usage = "Catalysts is an option to Transition and specifies \ the list of catalysts, which are of the form \!\(\* StyleBox[\"catalyst\",\nFontSlant->\"Italic\"]\) or {\!\(\* StyleBox[\"catalyst\",\nFontSlant->\"Italic\"]\), \!\(\* StyleBox[\"rate\",\nFontSlant->\"Italic\"]\)}."; *) (*Inhibitors::usage = "Inhibitors is an option to Transition and specifies \ the list of inhibitors, which are of the form \!\(\* StyleBox[\"inhibitor\",\nFontSlant->\"Italic\"]\) or {\!\(\* StyleBox[\"inhibitor\",\nFontSlant->\"Italic\"]\), rate}."; *) (*KineticLaw::usage = "KineticLaw is an option to Transition and specifies \ the kinetic law for the velocity of the transition. Automatic means that \ mass action will be used."; *) CompartmentalModelGraph::usage = "CompartmentalModelGraph[model] gives a \ directed Graph of model comprised of Transitions. Compartment symbols label the \ vertices and can be colored; rate constants label the edges. CompartmentalModelGraph[model, GraphTheme -> \"PetriNet\"] gives a \ directed Graph of model after converting it to a Petri net; \ vertices representing compartments are colored blue, and those representing \ Transitions are colored red."; (*GraphTheme::usage = "GraphTheme is an option for CompartmentalModelGraph that \ specifies an overall theme for visualization elements and styles. Valid values \ are \"PetriNet\", \"CompartmentsOnly\", and Automatic."; *) (*CompartmentColors::usage = "CompartmentColors is an option for CompartmentalModelGraph that \ specifies the rules for coloring vertices."; *) $C::usage = "Vertex wrapper for a compartment node in the output of CompartmentalModelGraph."; $R::usage = "Vertex wrapper for a rate node in the output of CompartmentalModelGraph."; ResolveCompartmentalModel::usage = "ResolveCompartmentalModel[model] returns a \ list of the independent components of the input compartmental model that is \ expressed as a list of transitions, which are considered coupled when they have \ compartments and/or symbolic rates in common."; KineticCompartmentalModel::usage = "KineticCompartmentalModel[model, t] returns an \ Association <|\"Variables\" -> v, \"ComponentTransitions\" -> c, \ \"StoichiometricMatrix\" -> S, \"Rates\" -> r, \"Equations\" -> e, \"TimeSymbol\" -> t|> \ for the set of transitions in model; the argument t is the time \ variable. The length of the component transition vector c may be longer than the \ input model because equilibria are decomposed into their forward and reverse reactions."; (*CheckForIndependentComponents::usage = "CheckForIndependentComponents is an input option to \ KineticCompartmentalModel and specifies whether or not the model is checked for independent \ components, which can be time consuming."; *) StoichiometryTable::usage = "StoichiometryTable[modelData] gives a formatted table \ of the stoichiometric information about the transitions in the modelData \ Association."; DeriveTransitions::usage = "DeriveTransitions[equations] returns a list of Transition \ objects derived from the input first order ordinary differential equations."; DefinePropensityFunction::usage = "DefinePropensityFunction[rates, sMatrix, vars, t] \ returns a Function or CompiledFunction that computes the propensity for the given \ state vector of length k and time. The argument rates is a list of n expressions for \ the rates of the events, sMatrix is the n\[Times]k stoichiometric matrix, vars is a list \ of k state variables used in the rate expressions, and t is the symbol for time used \ in the rate expressions."; StochasticSolve::usage = "StochasticSolve[rates, sMatrix, initalState, vars, \ {t, tmin, tmax}] returns a list of rules with InterpolatingFunction objects on the rhs \ for the stochastic simulation generated with the chemical master equation method given \ the list of n expressions for the rates of the events, the n\[Times]k stoichiometric matrix, \ the list of k state variables used in the rate expressions, and the list of time \ symbol, start time, and end time."; $CompartmentalModelingVersion::usage = "$CompartmentalModelingVersion is a string that gives \ the version of the CompartmentalModeling package you have loaded."; (* ::Section::Closed:: *) (*Other messages*) (* other messages *) General::fcomp = "The expressions `1` were detected and cannot be used \ as compartment variables."; KineticCompartmentalModel::ifcomp = "Warning: the compartment expression `1` appears to be \ ill-formed."; KineticCompartmentalModel::xtrm = "Warning: unexpected products or powers of compartments `1`."; KineticCompartmentalModel::ntzr = "Warning: compartment `1` has a net contribution \ of zero."; KineticCompartmentalModel::subsys = "Warning: the input consists of `1` independent \ subsystems."; StoichiometryTable::argt = "Argument `1` at position `2` is expected to be `3`."; General::unbal = "Warning: The transitions `` are not balanced."; General::npnuml = "The propensity `1` is not a vector of numbers when the arguments are `2`."; StochasticSolve::mxst = "Maximum number of `1` steps reached at the point `2` == `3`."; (* ::Section:: *) (*Begin*) Begin["`Private`"] ClearAll["`*"] $version = "1.7"; $releaseDate = DateObject[{2021, 6, 8}]; $CompartmentalModelingVersion = StringTemplate["`ver` (`date`)"][<| "ver" -> $version, "date" -> DateString[$releaseDate, {"MonthName", " ", "DayShort", ", ", "Year"}] |>]; (* ::Section::Closed:: *) (*Transition*) (* NullCompartment *) (*\[EmptySet] = NullCompartment*) (* typsetting rules *) NullCompartment /: MakeBoxes[NullCompartment, form:StandardForm | TraditionalForm] := "\[EmptySet]" (* Transition *) Options[Transition] = { "Catalysts" -> {}, "Inhibitors" -> {}, "KineticLaw" -> Automatic }; Transition[Transition[c1_, rate_, c2_, opts1:OptionsPattern[]], opts2:OptionsPattern[]] := Module[{opts}, opts = Sequence @@ Normal @ Association[opts1, Association[opts2]]; Transition[c1, rate, c2, opts] ] Transition[c1_, c2_, opts:OptionsPattern[]] := Transition[c1, Missing["NotSpecified"], c2, opts] (* rules for post-processing Transitions generated from equations *) Transition[c1_, rate_, Times[n_Integer?Negative, c2_]] := Transition[c1, -rate, -n c2] Transition[{(lhs_)..}, rate:{__}, {(rhs_)..}] := Transition[lhs, Total@rate, rhs] (* typsetting rules *) (* TODO: handle options *) Transition /: MakeBoxes[Transition[NullCompartment, rate_, c2:Except[NullCompartment]], form:StandardForm|TraditionalForm] := RowBox[{ Switch[rate, _Missing, "\[RightTeeArrow]", _, OverscriptBox["\[RightTeeArrow]", MakeBoxes[rate, form]] ], MakeBoxes[c2, form]} ] Transition /: MakeBoxes[Transition[c1:Except[NullCompartment], rate_, NullCompartment], form:StandardForm|TraditionalForm] := RowBox[{ MakeBoxes[c1, form], Switch[rate, _Missing, "\[RightArrowBar]", _, OverscriptBox["\[RightArrowBar]", MakeBoxes[rate, form]] ] }] Transition /: MakeBoxes[Transition[c1_, rate_, c2_], form:StandardForm|TraditionalForm] := RowBox[{ MakeBoxes[c1, form], Switch[rate, _Missing, "\[RightArrow]", _, OverscriptBox["\[RightArrow]", MakeBoxes[rate, form]] ], MakeBoxes[c2, form] }] (* not supported in InputForm input *) (*MakeExpression[RowBox[{"Transition", "[", RowBox[{"\[EmptySet]", ",", r_, ",", c2_}], "]"}], StandardForm] := MakeExpression[RowBox[{"Transition", "[", RowBox[{NullCompartment, ",", r, ",", c2}], "]"}], StandardForm]*) (* not supported in InputForm input *) (*MakeExpression[RowBox[{"Transition", "[", RowBox[{c1_, ",", r_, ",", "\[EmptySet]"}], "]"}], StandardForm] := MakeExpression[RowBox[{"Transition", "[", RowBox[{c1, ",", r, ",", NullCompartment}], "]"}], StandardForm]*) MakeExpression[RowBox[{OverscriptBox["\[RightTeeArrow]", r_], c2_}], StandardForm] := MakeExpression[RowBox[{"Transition", "[", RowBox[{NullCompartment, ",", r, ",", c2}], "]"}], StandardForm] MakeExpression[RowBox[{"\[RightTeeArrow]", c2_}], StandardForm] := MakeExpression[RowBox[{"Transition", "[", RowBox[{NullCompartment, ",", c2}], "]"}], StandardForm] MakeExpression[RowBox[{c1_, OverscriptBox["\[RightArrowBar]", r_]}], StandardForm] := MakeExpression[RowBox[{"Transition", "[", RowBox[{c1, ",", r, ",", NullCompartment}], "]"}], StandardForm] MakeExpression[RowBox[{c1_, "\[RightArrowBar]"}], StandardForm] := MakeExpression[RowBox[{"Transition", "[", RowBox[{c1, ",", NullCompartment}], "]"}], StandardForm] MakeExpression[RowBox[{c1_, OverscriptBox["\[RightArrow]"(* | "\[Rule]"*), r_], c2_}], StandardForm] := MakeExpression[RowBox[{"Transition", "[", RowBox[{c1, ",", r, ",", c2}], "]"}], StandardForm] MakeExpression[RowBox[{c1_, "\[RightArrow]" (*| "\[Rule]"*), c2_}], StandardForm] := MakeExpression[RowBox[{"Transition", "[", RowBox[{c1, ",", c2}], "]"}], StandardForm] MakeExpression[RowBox[{c1_, UnderoverscriptBox["\[Equilibrium]"(* | "\[LeftRightArrow]"*), r2_, r1_], c2_}], StandardForm] := MakeExpression[RowBox[{"{", RowBox[{ RowBox[{"Transition", "[", RowBox[{c1, ",", r1, ",", c2}], "]"}], ",", RowBox[{"Transition", "[", RowBox[{c2, ",", r2, ",", c1}], "]"}] }], "}"}], StandardForm] MakeExpression[RowBox[{c1_, OverscriptBox["\[Equilibrium]"(* | "\[LeftRightArrow]"*), r1_], c2_}], StandardForm] := MakeExpression[RowBox[{"{", RowBox[{ RowBox[{"Transition", "[", RowBox[{c1, ",", r1, ",", c2}], "]"}], ",", RowBox[{"Transition", "[", RowBox[{c2, ",", c1}], "]"}] }], "}"}], StandardForm] MakeExpression[RowBox[{c1_, UnderscriptBox["\[Equilibrium]"(* | "\[LeftRightArrow]"*), r2_], c2_}], StandardForm] := MakeExpression[RowBox[{"{", RowBox[{ RowBox[{"Transition", "[", RowBox[{c1, ",", c2}], "]"}], ",", RowBox[{"Transition", "[", RowBox[{c2, ",", r2, ",", c1}], "]"}] }], "}"}], StandardForm] MakeExpression[RowBox[{c1_, "\[Equilibrium]"(* | "\[LeftRightArrow]"*), c2_}], StandardForm] := MakeExpression[RowBox[{"{", RowBox[{ RowBox[{"Transition", "[", RowBox[{c1, ",", c2}], "]"}], ",", RowBox[{"Transition", "[", RowBox[{c2, ",", c1}], "]"}] }], "}"}], StandardForm] (* ::Section:: *) (*CompartmentalModelGraph*) (* CompartmentalModelGraph *) $C=.; $R=.; (* esf[{p1_, ___, p2_}, e:(u_\[DirectedEdge]v_)] /; MemberQ[$catalystEdges, e] := {Arrowheads[{{Small,1.0,{Graphics[Circle[{0,0},0.5]],0.5}}}], ColorData["Crayola"]["Green"], Arrow[{p1, p2}]} esf[{p1_, ___, p2_}, e:(u_\[DirectedEdge]v_)] := {Arrowheads[{{Small, 1}}], Arrow[{p1, p2}]} *) Options[CompartmentalModelGraph] = Sort @ Join[ {"GraphTheme" -> Automatic, "CompartmentColors" -> Automatic}, Options[Graph] ]; (* model may contain lists of Transitions due to use of equilibrium *) CompartmentalModelGraph[model:{(_Transition | {___Transition})...}, options:OptionsPattern[{CompartmentalModelGraph, "VertexColor" -> Automatic}]] := With[{attemptEvaluation=iCompartmentalModelGraph[Flatten @ model, options]}, attemptEvaluation /; MatchQ[attemptEvaluation, Except[$Failed]] ] iCompartmentalModelGraph[model_, options:OptionsPattern[{CompartmentalModelGraph, "VertexColor" -> Automatic}]] := Catch @ Module[{transitions, compartments, rootVertex, edges, catalystEdges, catalystVertices, inhibitorEdges, inhibitorVertices, vertices, i, baseStyle, theme, edgeLabels, vertexColor, compartmentColors, colorPattern, textColor}, theme = OptionValue["GraphTheme"]; If[MatchQ[theme, Except[Automatic | "CompartmentsOnly" | "PetriNet"]], Message[CompartmentalModelGraph::optvg, "GraphTheme", theme, "Automatic, \"CompartmentsOnly\", or \"PetriNet\""]; Return[$Failed] ]; vertexColor = OptionValue["VertexColor"]; If[MatchQ[vertexColor, Except[Automatic]], Message[CompartmentalModelGraph::obsalt, "VertexColor" -> vertexColor, "CompartmentColors" -> vertexColor] ]; compartmentColors = Switch[vertexColor, Automatic, OptionValue["CompartmentColors"], _, vertexColor]; colorPattern = _RGBColor | _CMYKColor | _Hue | _GrayLevel ; If[MatchQ[compartmentColors, Except[Automatic | Rule[_, colorPattern] | RuleDelayed[_, colorPattern] | {(Rule[_, colorPattern] | RuleDelayed[_, colorPattern])...}]], Message[CompartmentalModelGraph::optvg, "CompartmentColors", compartmentColors, "Automatic, a rule, or a list of rules"]; Return[$Failed] ]; baseStyle = {OptionValue[BaseStyle]} /. Automatic -> {} // Flatten; transitions = Union @ model /. { Transition[a_, r_, b_, opts___] :> (Transition[ a //. {Plus -> List, n_Integer c_ :> Table[c, n]} /. l_List :> Flatten[l], r, b //. {Plus -> List, n_Integer c_ :> Table[c, n]} /. l_List :> Flatten[l], opts ]) }(* // Echo*); compartments = Union @@ Join[ Cases[transitions, Transition[a_, r_, b_, opts___] :> Flatten[{a, b}]], Cases[transitions, HoldPattern["Catalysts" | "Inhibitors" -> m:{__}] :> First /@ Flatten /@ List /@ m, {2}] ] // DeleteCases[#, NullCompartment]&; rootVertex = FirstCase[OptionValue[GraphLayout], ("RootVertex" -> r_) | ("RootVertex" :> r_) :> r, Missing[], Infinity]; If[MatchQ[rootVertex, Except[_Missing]] && FreeQ[compartments, Verbatim[rootVertex], {1}], Message[CompartmentalModelGraph::optvg, "RootVertex", rootVertex, StringForm["one of ``", compartments]]; Throw[$Failed] ]; Switch[theme, "PetriNet", i = 0; transitions = transitions /. { Transition[a:NullCompartment, r_, b_, opts___] :> Transition[a, $R[++i, NullCompartment (* needed for esf for source *)][r], b, opts], Transition[a_, r_, b_, opts___] :> Transition[a, $R[++i][r], b, opts] }(* // Echo*); edges = transitions /. { Transition[NullCompartment, r_, b_, opts___] :> Thread @ DirectedEdge[r, Thread @ $C[b]], Transition[a_, r_, NullCompartment, opts___] :> Thread @ DirectedEdge[Thread @ $C[a], r], Transition[a_, r_, b_, opts___] :> Thread /@ {DirectedEdge[Thread @ $C[a], r], DirectedEdge[r, Thread @ $C[b]]} } // Flatten(* // Echo*); catalystEdges = Cases[transitions, Transition[a_, r_, b_, ___, "Catalysts" -> m:{__}, ___] :> Map[DirectedEdge[$C[#], r]&, First /@ Replace[m, mi:Except[{_, _}] :> {mi}, {1}]] ] // Flatten(* // Echo*); catalystVertices = {}; inhibitorEdges = Cases[transitions, Transition[a_, r_, b_, ___, "Inhibitors" -> m:{__}, ___] :> Map[DirectedEdge[$C[#], r]&, First /@ Replace[m, mi:Except[{_, _}] :> {mi}, {1}]] ] // Flatten(* // Echo*); inhibitorVertices = {}; edges = Join[ Style[#, ColorData["Crayola"]["MidnightBlue"], Thick] & /@ edges, Style[#, ColorData["Crayola"]["MountainMeadow"], Thick] & /@ catalystEdges, Style[#, ColorData["Crayola"]["Scarlet"], Thick] & /@ inhibitorEdges ](* // Echo*); vertices = Union @ Join[Join @@ List @@@ Cases[edges, _DirectedEdge, Infinity], catalystVertices, inhibitorVertices](* // Echo*); With[{vsf=With[{style=baseStyle}, Function[{position, vertex, box}, Switch[vertex, $C[_], Inset[Framed[Style[First@vertex, Black, Sequence @@ style], FrameStyle -> Blue, Background -> LightBlue], position ], $R[_, NullCompartment][_], Inset[Framed[Style[First@vertex, Black, Sequence @@ style], RoundingRadius -> 25, FrameStyle -> Red, Background -> LightRed], position ], $R[_][_], Inset[Framed[Style[First@vertex, Black, Sequence @@ style], FrameStyle -> Red, Background -> LightRed], position ] ] ] ]}, Graph[vertices, edges, VertexShapeFunction -> vsf(*, EdgeShapeFunction -> esf*), PerformanceGoal -> "Quality", FilterRules[{options}, Options[Graph]] /. ("RootVertex" -> root_) :> ("RootVertex" -> FirstCase[vertices, $C[root] | $R[__][root]]) ] ], _(* Automatic | "CompartmentsOnly" *), i = 0; transitions = transitions /. a:NullCompartment :> $C[++i][a](* // Echo*); edges = transitions /. Transition[a_, r_, b_, opts___] :> Thread @ DirectedEdge[a, b] // Flatten(* // Echo*); edgeLabels = DeleteCases[transitions, Transition[_, _Missing, __]] /. Transition[a_, r_, b_, opts___] :> Thread[(Thread @ DirectedEdge[a, b]) -> Framed[Style[r, Black, Sequence @@ baseStyle], FrameStyle -> White, Background -> White]] // Flatten(* // Echo*); catalystEdges = {}; catalystVertices = Cases[transitions, Transition[a_, r_, b_, ___, "Catalysts" -> m:{__}, ___] :> (First /@ Flatten /@ List /@ m) ] // Flatten(* // Echo*); inhibitorEdges = {}; inhibitorVertices = Cases[transitions, Transition[a_, r_, b_, ___, "Inhibitors" -> m:{__}, ___] :> (First /@ Flatten /@ List /@ m) ] // Flatten(* // Echo*); edges = Join[ Style[#, ColorData["Crayola"]["MidnightBlue"], Thick] & /@ edges, Style[#, ColorData["Crayola"]["MountainMeadow"], Thick] & /@ catalystEdges, Style[#, ColorData["Crayola"]["Scarlet"], Thick] & /@ inhibitorEdges ](* // Echo*); vertices = Union @ Join[Join @@ List @@@ Cases[edges, _DirectedEdge, Infinity], catalystVertices, inhibitorVertices](* // Echo*); compartmentColors = Replace[compartmentColors, Automatic -> {}] // Append[#, _ -> LightBlue]&; textColor = MapAt[If[First@ColorConvert[#, "Grayscale"] > 0.5, Black, White]&, compartmentColors, {All, 2}]; With[{vsf=With[{style=Sequence @@ baseStyle}, Function[{position, vertex, box}, If[MatchQ[vertex, $C[_][_]], (* NullCompartment *) Inset[First@vertex, position ], (* else *) (* all other compartemnts *) Inset[Framed[Style[vertex, vertex /. textColor, style], FrameStyle -> Darker@(vertex /. compartmentColors), Background -> (vertex /. compartmentColors)], position ] ] ] ]}, Graph[vertices, edges, VertexShapeFunction -> vsf, EdgeLabels -> If[MatchQ[theme, Automatic], edgeLabels, None], PerformanceGoal -> "Quality", FilterRules[{options}, Options[Graph]] /. ("RootVertex" -> root_) :> ("RootVertex" -> FirstCase[vertices, root]) ] ] ] ] (* ::Section::Closed:: *) (*ResolveCompartmentalModel*) (* ResolveCompartmentalModel *) Options[ResolveCompartmentalModel] = { }; ResolveCompartmentalModel[{}, options:OptionsPattern[]] := {} ResolveCompartmentalModel[OptionsPattern[]] /; Message[ResolveCompartmentalModel::argx, ResolveCompartmentalModel, 0, 1] := Null ResolveCompartmentalModel[t:(_Transition), options:OptionsPattern[]] := ResolveCompartmentalModel[{t}, options] ResolveCompartmentalModel[model:{(_Transition | {___Transition})...}, options:OptionsPattern[]] := With[{attemptEvaluation=iResolveCompartmentalModel[Flatten @ model, options]}, attemptEvaluation /; MatchQ[attemptEvaluation, Except[$Failed]] ] iResolveCompartmentalModel[transitions:{__Transition}, options:OptionsPattern[]] := Module[{transitionEdges, edges, graph, componentEdges}, edges = transitions /. { Transition[a_, r_, b_] :> (Transition[a //. {Plus -> List, n_Integer c_ :> Table[c, n]}, r, b //. {Plus -> List, n_Integer c_ :> Table[c, n]}] /. l_List :> Flatten[l]) } /. { Transition[c1_, r_, c2_] :> {Thread @ DirectedEdge[c1, r], Thread @ DirectedEdge[r, c2]} }; transitionEdges = Thread[{transitions, Flatten /@ edges}](* // Echo*); edges = Flatten @ edges // DeleteCases[#, DirectedEdge[NullCompartment, _] | DirectedEdge[_, NullCompartment]]&; graph = Graph[edges]; componentEdges = EdgeList /@ WeaklyConnectedGraphComponents[graph](* // Echo*); Cases[transitionEdges, {trns_, edg_} /; (MatchQ[Intersection[edg, #], {__}]) :> trns]& /@ componentEdges ] (* ::Section::Closed:: *) (*KineticCompartmentalModel*) (* KineticCompartmentalModel *) Attributes[debugTiming] = {HoldFirst}; debugTiming[expr_, label_] := expr // Timing // Echo[#, label, First] & // Last $forbiddenCompartmentsPattern = Alternatives[ C, D, E, N, O, _[], {___} ]; compartmentQ[c_] := (MatchQ[c, _Symbol] && !MatchQ[c, C | D | E | N | O]) || MatchQ[c, _String | Subscript[__] | Superscript[__] | Subsuperscript[__] | CenterDot[__] | SuperStar[_] | SuperDagger[_] | SuperPlus[_] | SuperMinus[_] | _String[__?AtomQ] | (h_Symbol[__?AtomQ] /; Context[h] != "System`") ] parameterQ[c_] := (MatchQ[c, _Symbol] && !MatchQ[c, C | D | E | N | O]) || MatchQ[c, _?NumericQ | _Times | _Plus | Power[_, -1] | Subscript[__] | Superscript[__] | Subsuperscript[__] | CenterDot[__] | SuperStar[_] | SuperDagger[_] | SuperPlus[_] | SuperMinus[_] | _String[__?AtomQ] | (h_Symbol[__?AtomQ] /; Context[h] != "System`") ] Options[getCompartments] = { "Compartments" -> Automatic }; getCompartments[transitionList:{___Transition}, callingFunction_Symbol, opts:OptionsPattern[]] := Catch @ Module[{transitions, compartmentExressions, badComp, compartments}, transitions = Union @ transitionList; compartments = OptionValue["Compartments"] // Replace[#, c:Except[Automatic] :> Flatten[{c}]]&; compartmentExressions = Switch[{transitionList, compartments}, {{___}, Automatic}, Join @@ Cases[transitions, Transition[c1_, _, c2_] :> {c1, c2}] // Union(* // Echo[#, "compartment expressions"]&*), {{}, {__}}, compartments, {{__}, _}, Message[callingFunction::optvg, "Compartments", compartments, "Automatic when the transition list is not empty"]; Throw[$Failed], {{}, {}}, Message[callingFunction::optvg, "Compartments", compartments, "a nonempty list of compartments or Automatic when the transition list is empty"]; Throw[$Failed] ](* // Echo[#, "compartmentExressions"]&*); compartmentExressions = Join @@ Replace[Expand[compartmentExressions], {HoldPattern[Plus[c__]] :> {c}, c_ :> {c}}, {1}] /. Times[n_Integer, c_] :> c // Union(* // Echo[#, "compartmentExressions"]&*); (* check for known problems *) badComp = Union @ Cases[compartmentExressions, $forbiddenCompartmentsPattern(*, Infinity*)]; If[MatchQ[badComp, {__}], Message[callingFunction::fcomp, badComp]; Throw[$Failed] ]; compartments = Variables[compartmentExressions] // DeleteCases[#, NullCompartment]&(* // Echo[#, "compartments"]&*); (* check for potential problems *) Scan[If[!compartmentQ[#], Message[callingFunction::ifcomp, #]]&, compartments]; compartments ] balancedTransitionQ[Transition[lhs_, _, rhs_]] := Module[{vars}, vars = Union @@ Variables /@ {lhs, rhs}; 0 == Total@Coefficient[rhs - lhs, vars] ] Options[KineticCompartmentalModel] = Sort @ Join[ Options[getCompartments], { "CheckForIndependentComponents" -> True } ]; KineticCompartmentalModel[_, OptionsPattern[]] /; Message[KineticCompartmentalModel::argr, KineticCompartmentalModel, 2] := Null KineticCompartmentalModel[trans:(_Transition), t_, options:OptionsPattern[]] := KineticCompartmentalModel[{trans}, t, options] KineticCompartmentalModel[transitionList:{(_Transition | {___Transition})...}, t_Symbol, options:OptionsPattern[]] := With[{attemptEvaluation=iKineticCompartmentalModel[Flatten @ transitionList, t, options]}, attemptEvaluation /; MatchQ[attemptEvaluation, Except[$Failed]] ] iKineticCompartmentalModel[transitionList:{(_Transition)...}, t_Symbol, options:OptionsPattern[KineticCompartmentalModel]] := Module[{checkComponents, nComp, transitions, reactants, rates, products, stoichiometry, crossTerms, defaultOrders, orders, orderRules, terms, lhs, rhs, optlist, ok, compartments, vars, eqns, flux}, optlist = {options} // Flatten; ok = True; Scan[If[FreeQ[Keys @ Options[KineticCompartmentalModel], Keys @ #], Message[KineticCompartmentalModel::optx, #, HoldForm[KineticCompartmentalModel[transitionList, t, options]]]; ok=False]&, optlist]; If[!ok, Return[$Failed]]; checkComponents = OptionValue["CheckForIndependentComponents"]; If[MatchQ[checkComponents, Except[True | False]], Message[KineticCompartmentalModel::opttf, "CheckForIndependentComponents", checkComponents]; Return[$Failed] ]; transitions = (Union @ transitionList); compartments = getCompartments[transitions, KineticCompartmentalModel, FilterRules[{options}, Options[getCompartments]]]; If[compartments === $Failed, Return[$Failed]]; If[checkComponents, nComp=Length[ResolveCompartmentalModel[transitionList, FilterRules[{options}, Options[ResolveCompartmentalModel]]]](* // Echo[#, "ResolveCompartmentalModel"]&*); If[nComp > 1, Message[KineticCompartmentalModel::subsys, nComp]] ]; {reactants, rates, products} = Replace[Transpose @ (List @@@ transitions), {} -> {{}, {}, {}}]; crossTerms = Union @ Cases[ Replace[{reactants, products}, Times[_?NumericQ, x__] :> Times[x], {2}], _Times | _Power, {2} ](* // Timing // Echo[#, "crossTerms", First] & // Last*); If[MatchQ[crossTerms, {__}], Message[KineticCompartmentalModel::xtrm, crossTerms] ]; vars = Through[compartments[t]](* // Echo[#, "vars"]&*); stoichiometry = Outer[Coefficient, products - reactants, compartments](* // Echo[#, "stoichiometry"]&*); defaultOrders = (Thread[compartments -> #1]&) /@ Outer[Coefficient, reactants, compartments](* // Echo[#, "defaultOrders"]&*); orders = compartments /. #& /@ Dispatch /@ defaultOrders(* // Echo[#, "orders"]&*); orderRules = ((Thread[compartments -> #1]&) /@ (vars^#1&) /@ orders)(* // Echo[#, "orderRules"]&*); terms = rates (Times @@ compartments /. Dispatch[orderRules])(* // Echo[#, "terms"]&*); flux = terms; lhs = D[vars, t]; rhs = flux . stoichiometry(* // Timing // Echo[#, "rhs", First] & // Last*); eqns = lhs == rhs // Thread; (* check for possible misspecification of model *) Scan[If[MatchQ[#, Derivative[1][_][t] == 0], Message[KineticCompartmentalModel::ntzr, #[[1, 0, 1]]] ]&, eqns]; <| "Variables" -> Head /@ vars, "ComponentTransitions" -> transitions, "StoichiometricMatrix" -> stoichiometry, "Rates" -> flux, "Equations" -> eqns, "TimeSymbol" -> t |> ] (* ::Section::Closed:: *) (*StoichiometryTable*) (* StoichiometryTable *) StoichiometryTable[modelData_Association, opts:OptionsPattern[Grid]] := Module[{keys, type, vars, trans, transHeading, sMatrix, rates, k}, keys = Keys @ modelData; type = Which[ Intersection[keys, {"Variables", "ComponentReactions", "StoichiometricMatrix", "Rates"}] == {"ComponentReactions", "Rates", "StoichiometricMatrix", "Variables"}, "SystemsBiologyModel", Intersection[keys, {"Variables", "ComponentTransitions", "StoichiometricMatrix", "Rates"}] == {"ComponentTransitions", "Rates", "StoichiometricMatrix", "Variables"}, "CompartmentalModel", True, $Failed ]; If[MatchQ[type, $Failed], Message[StoichiometryTable::argt, modelData, 1, "an Association with keys \"Variables\", \"ComponentTransitions\", \"StoichiometricMatrix\", \"Rates\""]; Return[$Failed] ]; vars = modelData["Variables"]; k = Length @ vars; {trans, transHeading} = Switch[type, "CompartmentalModel", {modelData["ComponentTransitions"], "Transition"}, "SystemsBiologyModel", {modelData["ComponentReactions"], "Reaction"} ]; sMatrix = modelData["StoichiometricMatrix"]; rates = modelData["Rates"]; MapThread[Join, {List /@ trans, sMatrix, List /@ rates}] // Join[ {{"", "Stoichiometry", Sequence @@ ConstantArray[SpanFromLeft, k - 1], ""}, {transHeading, Sequence @@ vars, (*"Propensity"*)"Rate"}}, #]& // Grid[#, Dividers -> {{False, True, {False}, True, False}, {False, True, True, False}}, Alignment -> {{Center, {Right}, Center}, Automatic, {{1, 2} -> Center, Sequence@@Map[{2, # + 1} -> Center&, Range@k]}}, Sequence @@ FilterRules[{opts}, Options@Grid], Spacings -> {{{1.5}}, {{0.5}}} ]& // Framed ] (* DeriveTransitions *) canonicalizeTerm[term_] := Replace[Abs[term], Times[_Integer, r__] :> Times[r]] extractTransitions[{terms:{__}, eqns:{(Derivative[1][_][t_] == _) ..}}, vars:{_[t_] ..}] := Module[{termVars, lhs, rateAndRhs, rate, rhs}, (*Echo["---------------------------"]; Echo[terms, "terms"]; Echo[eqns, "eqns"]; Echo[vars, "vars"]; *) termVars = Cases[terms, Alternatives @@ vars, \[Infinity]] // Union // Apply[Alternatives, #] &(* // Echo[#, "termVars"]&*); lhs = ((Exponent[First@terms, #] & /@ vars).Head /@ vars) /. (0 -> NullCompartment)(* // Echo[#, "lhs"]&*); rateAndRhs = Cases[eqns, (dv : Except[D[#, t] & /@ termVars]) == (Plus[term : Alternatives @@ terms, __] | term : Alternatives @@ terms) :> {term, {Replace[dv, Derivative[1][v_][t] :> v], Replace[term, {Times[n_Integer, __] :> n, _ -> 1}]}} ](* // Echo[#, "rateAndRhs"]&*); If[MatchQ[rateAndRhs, {}], rateAndRhs = Cases[eqns, (dv : D[#, t] & /@ termVars) == (Plus[term : Alternatives @@ terms, __] | term : Alternatives @@ terms) :> {term, {Replace[dv, Derivative[1][v_][t] :> v], Replace[term, {Times[n_Integer, __] :> Max[0, n], _ -> 1}]}} ](* // Echo[#, "rateAndRhs"]&*) ]; {rate, rhs} = rateAndRhs // {First@First@#, Last /@ #} &(* // Echo[#, "{rate, rhs}"]&*); rate = Replace[rate, Times[_Integer, r__] :> Times[r]] /. Thread[vars -> 1]; rhs = (Plus @@ Times @@@ rhs) /. (0 -> NullCompartment); Transition[lhs, rate, rhs] ] Options[DeriveTransitions] = { "CheckBalance" -> False }; DeriveTransitions[eqn:(Derivative[1][_][_] == _), opts:OptionsPattern[]] := DeriveTransitions[{eqn}, opts] DeriveTransitions[{}, opts:OptionsPattern[]] := {} DeriveTransitions[eqns : {(Derivative[1][_][t_] == _) ...}, opts:OptionsPattern[]] := With[{attemptEvaluation=iDeriveTransitions[eqns, opts]}, attemptEvaluation /; MatchQ[attemptEvaluation, Except[$Failed]] ] (* TODO: need to ensure that equations are first order ODEs *) iDeriveTransitions[eqnsIn:{__Equal}, opts:OptionsPattern[DeriveTransitions]] := Catch @ Module[{checkQ, eqns=Expand@eqnsIn, vars, terms, rateTerms, termsAndEqns, transitions, unbalanced}, checkQ = OptionValue["CheckBalance"]; If[MatchQ[checkQ, Except[True | False]], Message[DeriveTransitions::opttf, "checkBalance", checkQ]; Throw[$Failed] ]; vars = Cases[eqns, v_'[t_] :> v[t], Infinity] // Union; terms = MapAt[Replace[#, {p_Plus :> List @@ p, x_ :> {x}}] &, List @@@ eqns, {All, 2}]; rateTerms = Union @@ Last /@ terms; rateTerms = Values@GroupBy[rateTerms, canonicalizeTerm]; termsAndEqns = {#, Cases[eqns, _ == rhs_ /; MemberQ[rhs, Alternatives @@ #, {0, 1}]]} & /@ rateTerms; transitions = extractTransitions[#, vars] & /@ termsAndEqns(* // Echo*); transitions = Replace[ transitions, {Transition[lhs : NullCompartment, r_, p_Plus] :> Splice[List @@ Thread[Transition[lhs, r, p], Plus]], Transition[lhs_, r_, p_Plus] /; MemberQ[p, Times[_?Negative, __]] :> Splice[List @@ Thread[Transition[lhs, r, p], Plus]]}, {1}](* // Echo*); transitions = Sort @ Values @ GroupBy[transitions, {First@#, Last@#} &, Thread[#, Transition] &]; If[checkQ && MatchQ[unbalanced=Select[transitions, balancedTransitionQ /* Not], {__}], Message[DeriveTransitions::unbal, unbalanced] ]; transitions ] (* ::Section::Closed:: *) (*DefinePropensityFunction*) (* DefinePropensityFunction *) Options[DefinePropensityFunction] = { Compiled -> True }; DefinePropensityFunction[rates : {__}, sMatrix_?(MatrixQ[#1, IntegerQ] &), vars : {__}, time_Symbol, opts:OptionsPattern[]] := With[{attemptEvaluation=iDefinePropensityFunction[rates, sMatrix, vars, time, opts]}, attemptEvaluation /; MatchQ[attemptEvaluation, Except[$Failed]] ] iDefinePropensityFunction[rates_, sMatrix_, varsIn_, time_, opts:OptionsPattern[DefinePropensityFunction]] := Module[{vars, compileQ}, compileQ = OptionValue[Compiled]; If[MatchQ[compileQ, Except[True | False]], Message[DefinePropensityFunction::opttf, Compiled, compileQ]; Return[$Failed] ]; If[Abs@Max[sMatrix] == 1, dpf[rates, sMatrix, varsIn, time, compileQ], (* else *) dpfWithMask[rates, sMatrix, varsIn, time, compileQ] ] ] dpf[rates_, sMatrix_, varsIn_, time_, compileQ_] := If[compileQ, With[{vars = Through[varsIn[time]], states = Quiet[Table[s[[i]], {i, Length[varsIn]}]]}, With[{prop = rates /. Thread[vars -> states] /. time -> t}, Compile[{{s, _Integer, 1}, {t, _Real}}, prop, CompilationTarget -> "WVM"(*"C"*) ] ] ], (* else *) (* s$ and t$ needed for localization in Function *) With[{vars = Through[varsIn[time]], states = Quiet[Table[s$[[i]], {i, Length[varsIn]}]]}, With[{prop = rates /. Thread[vars -> states] /. time -> t$}, Function[{s, t}, prop] ] ] ] dpfWithMask[rates_, sMatrix_, varsIn_, time_, compileQ_] := If[compileQ, With[{vars = Through[varsIn[time]], states = Quiet[Table[s[[i]], {i, Length[varsIn]}]]}, With[{prop = rates /. Thread[vars -> states] /. time -> t, mask = Function[row, Boole[AllTrue[states + row, NonNegative]]] /@ sMatrix}, Compile[{{s, _Integer, 1}, {t, _Real}}, prop mask, CompilationTarget -> "WVM"(*"C"*) ] ] ], (* else *) (* s$ and t$ needed for localization in Function *) With[{vars = Through[varsIn[time]], states = Quiet[Table[s$[[i]], {i, Length[varsIn]}]]}, With[{prop = rates /. Thread[vars -> states] /. time -> t$, mask = Function[row, Boole[AllTrue[states + row, NonNegative]]] /@ sMatrix}, Function[{s, t}, prop mask] ] ] ] (* ::Section::Closed:: *) (*StochasticSolve*) (* StochasticSolve *) $propensity = Null; Options[StochasticSolve] = { Compiled -> True, MaxSteps -> Automatic }; StochasticSolve[rates:{__}, (sMatrix_)?(MatrixQ[#1, IntegerQ] & ), (stateIn_)?(VectorQ[#1, IntegerQ] & ), vars:{__}, {t_Symbol, (tmin_)?NumberQ, (tmax_)?NumberQ}, opts:OptionsPattern[]] /; Dimensions[sMatrix] == Length /@ {rates, stateIn} && Length[vars] == Length[First[sMatrix]] := With[{attemptEvaluation=iStochasticSolve[rates, sMatrix, stateIn, vars, {t, tmin, tmax}, opts]}, attemptEvaluation /; MatchQ[attemptEvaluation, Except[$Failed]] ] iStochasticSolve[rates_, sMatrix_, stateIn_, vars_, {t_, tmin_, tmax_}, opts:OptionsPattern[StochasticSolve]] := Module[{compiledQ, maxSteps, p, time, state, times, data}, compiledQ = OptionValue[Compiled]; If[MatchQ[compiledQ, Except[True | False]], Message[StochasticSolve::opttf, Compiled, compiledQ]; Return[$Failed] ]; maxSteps = OptionValue[MaxSteps]; If[MatchQ[maxSteps, Except[_Integer?NonNegative | Infinity | Automatic]], Message[StochasticSolve::ioppfa, MaxSteps, maxSteps]; Return[$Failed] ]; If[MatchQ[maxSteps, Automatic | Infinity], maxSteps = 500000 ]; $propensity = DefinePropensityFunction[rates, sMatrix, vars, t, Compiled -> compiledQ(*False*)]; (* test for fully numeric result *) If[MatchQ[p=$propensity[stateIn, tmin], Except[{__Real}]], Message[StochasticSolve::npnuml, p, <|"state" -> stateIn, t -> tmin|>]; Return[$Failed] ]; data = If[compiledQ(*False*), stochasticSolveLoopCompiled[stateIn, tmin, tmax, sMatrix, maxSteps], (* else *) stochasticSolveLoop[stateIn, tmin, tmax, sMatrix, maxSteps] ]; If[Length@data == maxSteps && data[[-1, 1]] < tmax, Message[StochasticSolve::mxst, maxSteps, t, data[[-1, 1]]] ]; times = First /@ data; data = Transpose[Rest /@ data]; data = (Join @@ Replace[SplitBy[Thread[{times, #1}], Last], {a_, __} :> {a}, {1}])& /@ data; (* the following is a hack to force the interpolation the behave the desired way *) data = Module[{x, y}, {x, y} = Thread @ #; x = Join[{First@x}, Rest@x, {Last@x + Sqrt[$MachineEpsilon]}]; y = Join[{First@y}, Most@y, {Last@y}]; Thread[{x, y}] ]& /@ data; Thread[vars -> (Interpolation[#, InterpolationOrder -> 0] & ) /@ data] ] stochasticSolveLoop = Function[{stateIn, tmin, tmax, sMatrix, maxSteps}, Module[{time = tmin, state = stateIn, stateList, cum, rand, result, step = 0}, stateList = CreateDataStructure["DynamicArray"]; step++; stateList["Append", (*Echo[#, "time & state"]&@*)Prepend[N@state, time]]; While[step < maxSteps && time < tmax && Last[cum = (*Echo[#, "cum"]&@*)Accumulate[(*Echo[#, "propensity"]&@*)$propensity[state, time]]] > 0, time += Log[1./RandomReal[]]/Last[cum]; rand = (*Echo[#, "rand"]&@*)RandomReal[Last[cum]]; state += sMatrix[[(*Echo[#, "idx"]&@*)First[FirstPosition[cum, c_ /; rand <= c]]]]; step++; stateList["Append", (*Echo[#, "time & state"]&@*)Prepend[N@state, time]] ]; result = Normal[stateList]; stateList["DropAll"]; result ] ] stochasticSolveLoopCompiled = Compile[{{stateIn, _Integer, 1}, {tmin, _Real}, {tmax, _Real}, {sMatrixIn, _Integer, 2}, {maxSteps, _Integer}}, Module[{time = tmin, state = N@stateIn, sMatrix = N@sMatrixIn, cum = ConstantArray[0., Length@sMatrix], rand, s = 1, result = ConstantArray[0., {maxSteps, Length@stateIn + 1}], step = 0}, (*result[[++step]] = Prepend[N@state, time];*) result[[++step, 1]] = time; Do[result[[step, i + 1]] = state[[i]],{i, Length[state]}]; While[step < maxSteps && time < tmax && Last[cum = Accumulate[$propensity[state, time]]] > 0, time += Log[1./RandomReal[]]/Last[cum]; rand = RandomReal[Last[cum]]; For[s = 1, s <= Length@cum && rand > cum[[s]], s++, Null]; state += sMatrix[[s]]; (*result[[++step]] = Prepend[N@state, time];*) result[[++step, 1]] = time; Do[result[[step, i + 1]] = state[[i]],{i, Length[state]}]; ]; Take[result, step] ], {{propensity, _Real, 1}}, CompilationOptions -> {"InlineExternalDefinitions" -> True}, "CompilationTarget" -> "WVM"(*, CompilationTarget -> "C"*) ] (* ::Section:: *) (*End, EndPackage*) End[]; EndPackage[];