(* ::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. *) (* :Copyright: (c) 2020 by Robert B. Nachbar *) (* :Package Version: 1.1 *) (* :Mathematica Version: 12.1 *) (* :History: 1.0 14-Mar-2020 1.1 08-Oct-2020 *) (* :Caveats: *) (* :ToDo: o Transition with Missing[] rate o 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 o 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 *) Transition::usage = "Transition[c1, r, c2] specifies a Transition from compartment \ c1 to compartment c2 with transfer rate 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 specifiec \ with Plus and Times, respectively. Transition[\[EmptySet], r, c2] specifies a source (boundary or birth) to compartment c2 with transfer \ rate r, and \[EmptySet] is the \\[EmptySet] character. 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, \[EmptySet]] specifies a sink (boundary or death) from compartment c1 with transfer \ rate r, and \[EmptySet] is the \\[EmptySet] character. 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]."; 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 \ Graph of model comprised of Transitions after converting to a Petri net. \ Vertices representing compartments are colored blue, and those representing \ Transitions are colored red."; IncludeRates::usage ="IncludeRates is an input option to CompartmentalModelGraph \ that specifies whether or not rates are included in the output Graph."; $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, \"Time\" -> 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."; 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."; (* ::Section::Closed:: *) (*Other messages*) (* other messages *) KineticCompartmentalModel::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::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["`*"] (* ::Section::Closed:: *) (*Transition*) (* Transition *) Options[Transition] = { Catalysts -> {}, Inhibitors -> {}, KineticLaw -> Automatic }; Transition[Transition[c1_, rate_, c2_], opts:OptionsPattern[]] := Transition[c1, rate, c2, opts] Transition /: MakeBoxes[Transition[Global`\[EmptySet], rate_, c2_], form:StandardForm|TraditionalForm] := RowBox[{ Switch[rate, _Missing, "\[RightTeeArrow]", _, OverscriptBox["\[RightTeeArrow]", MakeBoxes[rate, form]] ], MakeBoxes[c2, form]} ] Transition /: MakeBoxes[Transition[c1_, rate_, Global`\[EmptySet]], 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] }] MakeExpression[RowBox[{OverscriptBox["\[RightTeeArrow]", r_], c2_}], StandardForm]:= MakeExpression[RowBox[{"Transition", "[", RowBox[{Global`\[EmptySet], ",", r, ",", c2}], "]"}], StandardForm] MakeExpression[RowBox[{"\[RightTeeArrow]", c2_}], StandardForm]:= MakeExpression[RowBox[{"Transition", "[", RowBox[{Global`\[EmptySet], ",", RowBox[{"Missing", "[", "\"\\"", "]"}], ",", c2}], "]"}], StandardForm] MakeExpression[RowBox[{c1_, OverscriptBox["\[RightArrowBar]", r_]}], StandardForm]:= MakeExpression[RowBox[{"Transition", "[", RowBox[{c1, ",", r, ",", Global`\[EmptySet]}], "]"}], StandardForm] MakeExpression[RowBox[{c1_, "\[RightArrowBar]"}], StandardForm]:= MakeExpression[RowBox[{"Transition", "[", RowBox[{c1, ",", RowBox[{"Missing", "[", "\"\\"", "]"}], ",", Global`\[EmptySet]}], "]"}], 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, ",", RowBox[{"Missing", "[", "\"\\"", "]"}], ",", 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, ",", RowBox[{"Missing", "[", "\"\\"", "]"}], ",", c1}], "]"}] }], "}"}], StandardForm] MakeExpression[RowBox[{c1_, UnderscriptBox["\[Equilibrium]" | "\[LeftRightArrow]", r2_], c2_}], StandardForm] := MakeExpression[RowBox[{"{", RowBox[{ RowBox[{"Transition", "[", RowBox[{c1, ",", RowBox[{"Missing", "[", "\"\\"", "]"}], ",", c2}], "]"}], ",", RowBox[{"Transition", "[", RowBox[{c2, ",", r2, ",", c1}], "]"}] }], "}"}], StandardForm] MakeExpression[RowBox[{c1_, "\[Equilibrium]" | "\[LeftRightArrow]", c2_}], StandardForm] := MakeExpression[RowBox[{"{", RowBox[{ RowBox[{"Transition", "[", RowBox[{c1, ",", RowBox[{"Missing", "[", "\"\\"", "]"}], ",", c2}], "]"}], ",", RowBox[{"Transition", "[", RowBox[{c2, ",", RowBox[{"Missing", "[", "\"\\"", "]"}], ",", 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] = Join[ {IncludeRates -> True}, Options[Graph] ]; CompartmentalModelGraph[model:{(_Transition | {___Transition})...}, options:OptionsPattern[CompartmentalModelGraph]] := With[{attemptEvaluation=iCompartmentalModelGraph[Flatten @ model, options]}, attemptEvaluation /; MatchQ[attemptEvaluation, Except[$Failed]] ] iCompartmentalModelGraph[model_, options:OptionsPattern[CompartmentalModelGraph]] := Module[{transitions, compartments, catalysts, inhibitors, kineticLaws, edges, catalystEdges, catalystVertices, inhibitorEdges, inhibitorVertices, vertices, i = 0, baseStyle, includeRatesQ}, includeRatesQ = OptionValue[IncludeRates]; If[MatchQ[includeRatesQ, Except[True | False]], Message[CompartmentalModelGraph::opttf, IncludeRates, includeRatesQ]; 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]}, r, b //. {Plus -> List, n_Integer c_ :> Table[c, n]}, opts ] /. l_List :> Flatten[l]) }(* // Echo*); i=0; transitions = transitions /. { Transition[a:Global`\[EmptySet], r_, b_, opts___] :> Transition[a, $R[++i, Global`\[EmptySet] (* needed for esf for source *)][r], b, opts], Transition[a_, r_, b_, opts___] :> Transition[a, $R[++i][r], b, opts] }(* // Echo*); edges = transitions /. { Transition[Global`\[EmptySet], r_, b_, opts___] :> If[includeRatesQ, Thread @ DirectedEdge[r, Thread @ $C[b]], (* else *) Nothing ], Transition[a_, r_, Global`\[EmptySet], opts___] :> If[includeRatesQ, Thread @ DirectedEdge[Thread @ $C[a], r], (* else *) Nothing ], Transition[a_, r_, b_, opts___] :> If[includeRatesQ, Thread /@ {DirectedEdge[Thread @ $C[a], r], DirectedEdge[r, Thread @ $C[b]]}, (* else *) Thread @ DirectedEdge[Thread @ $C[a], Thread @ $C[b]] ] } // Flatten(* // Echo*); If[includeRatesQ, catalystEdges = Cases[transitions, Transition[a_, r_, b_, ___, Catalysts -> m:{__}, ___] :> Map[DirectedEdge[$C[#], r]&, First /@ Flatten /@ List /@ m] ] // Flatten(* // Echo*); catalystVertices = {}; inhibitorEdges = Cases[transitions, Transition[a_, r_, b_, ___, Inhibitors -> m:{__}, ___] :> Map[DirectedEdge[$C[#], r]&, First /@ Flatten /@ List /@ m] ] // Flatten(* // Echo*); inhibitorVertices = {}, (* else *) catalystEdges = {}; catalystVertices = Cases[transitions, Transition[a_, r_, b_, ___, Catalysts -> m:{__}, ___] :> Map[$C, First /@ Flatten /@ List /@ m] ] // Flatten(* // Echo*); inhibitorEdges = {}; inhibitorVertices = Cases[transitions, Transition[a_, r_, b_, ___, Inhibitors -> m:{__}, ___] :> Map[$C, 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*); (* TODO: check that root vertex is present *) 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[_, Global`\[EmptySet]][_], 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]]) ] ] ] (* ::Section::Closed:: *) (*ResolveCompartmentalModel*) (* ResolveCompartmentalModel *) Options[ResolveCompartmentalModel] = { }; ResolveCompartmentalModel[OptionsPattern[]] /; Message[ResolveCompartmentalModel::argx, ResolveCompartmentalModel, 0, 1] := Null ResolveCompartmentalModel[t:(_Transition), options:OptionsPattern[]] := ResolveCompartmentalModel[{t}, options] ResolveCompartmentalModel[model:{(_Transition | {___Transition})...}, options:OptionsPattern[]] := Module[{transitions, transitionEdges, edges, graph, componentEdges}, transitions = Flatten @ model; 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[Global`\[EmptySet], _] | DirectedEdge[_, Global`\[EmptySet]]]&; 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[ D, N, O, _[], {___} ]; compartmentQ[c_] := MatchQ[s, _Symbol | _String | Subscript[__] | Superscript[__] | Subsuperscript[__] | CenterDot[__] | SuperStar[_] | SuperDagger[_] | SuperPlus[_] | SuperMinus[_] | Row[{__}, _] | _String[__?AtomQ] | (h_Symbol[___?AtomQ](* /; FreeQ[Context[h], "System`"]*)) ] Options[KineticCompartmentalModel] = { 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, compartmentExressions, badComp, 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] ]; If[checkComponents, nComp=Length[ResolveCompartmentalModel[transitionList]](* // Timing // Echo[#, "ResolveCompartmentalModel", First] & // Last*); If[nComp > 1, Message[KineticCompartmentalModel::subsys, nComp]] ]; transitions = (Union @ transitionList); compartmentExressions = Join @@ Cases[transitions, Transition[c1_, _, c2_] :> {c1, c2}] // Union(* // Echo[#, "compartment expressions"]&*); (* check for known problems *) badComp = Union @ Cases[compartmentExressions, forbiddenCompartmentsPattern, Infinity]; If[MatchQ[badComp, {__}], Message[KineticCompartmentalModel::fcomp, badComp]; Return[$Failed] ]; compartments = Variables[compartmentExressions] // DeleteCases[#, Global`\[EmptySet]]&(* // Echo[#, "compartments"]&*); (* check for potential problems *) Scan[If[!compartmentQ[#], Message[KineticCompartmentalModel::ifcomp, #]]&, compartments]; {reactants, rates, products} = Transpose @ Replace[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]]; stoichiometry = Outer[Coefficient, products - reactants, compartments](* // Timing // Echo[#, "stoichiometry", First] & // Last*); defaultOrders = (Thread[compartments -> #1]&) /@ Outer[Coefficient, reactants, compartments](* // Timing // Echo[#, "defaultOrders", First] & // Last*)(* // Echo[#, "defaultOrders"]&*); orders = compartments /. Dispatch[defaultOrders](* // Timing // Echo[#, "orders", First] & // Last*)(* // Echo[#, "orders"]&*); orderRules = ((Thread[compartments -> #1]&) /@ (vars^#1&) /@ orders)(* // Timing // Echo[#, "orderRules", First] & // Last*)(* // Echo[#, "orderRules"]&*); terms = rates (Times @@ compartments /. Dispatch[orderRules])(* // Timing // Echo[#, "terms", First] & // Last*); 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, "Time" -> 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 ] (* ::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[];