To run the notebook, choose Kernel->Evaluation->Evaluate Notebook.
        To run individual lines of code, click the line of code and type Shift+Enter.
        (Be sure to evaluate all the lines of code hidden in the Modules section before attempting to run the examples below. To un-hide the Modules section, double-click the blue downward arrow at the far right.)

WARNING! This notebook runs *incredibly* slowly with version 5. It takes about 15 minutes to execute the entire notebook with version 4.2; I recommend you use this previous version (over 2 hours with 5.2).

Turns off irrelevant warnings:

Off[General::spell]
Off[General::spell1]
Off[General::unfl]
Off[ReplaceAll::reps]
Off[InterpolatingFunction::dmval]

Loads auxiliary packages required below:

<< LinearAlgebra`MatrixManipulation`
<< Graphics`Graphics`

Modules

Modeules are little programs - The variable_ parameters appearing in the square brackets are the program inputs. The variables appearing in the {} braces are specific to that particular Module. Return[output] is the program output.

General Modules

These are two short-cuts I use below to integrate exponentials.

Matt[aVar_, bVar_] := 
        Exp[Log[aVar]*bVar];
Matt2[Exp[1], bVar_] := 
        1/(s - D[bVar, tp]);
Matt2[aVar_, bVar_] := Power[aVar, bVar];

Here are two different ways to calculate the matrix exponential - one using the standard routine supplied by Mathematica, and one using a digonalization. The standard routine works much faster.

myExpMat[matrix_] := MatrixExp[matrix*tp];

myExpMatDiag[matrix_] := 
        Module[{LVar, VVar, pVar, p1Var}, 
LVar = 
            Eigenvalues[matrix]; 
VVar = 
            Eigenvectors[matrix]; 
pVar = 
            Transpose[VVar]; 
p1Var = 
            Inverse[pVar]; 

Return[
            pVar . DiagonalMatrix[Exp[LVar*tp]] . 
              p1Var]
];

A set of modules to extract the solutions of interest from a long list of numerical roots of a high-degree polynomial.

myRootGet[{aVar_ -> bVar_}] = bVar;

myRootCheck[aVar_ -> bVar_] := 
        If[Abs[Im[bVar]] < 10^(-8) && 
            Re[bVar] ≥ 0, {aVar -> Re[Chop[bVar]]}, 
          Continue[]];

myRootExtraction[ssSoln_] := 
        Module[{lenS, rapSoln, kInd, iInd, jInd, rap, 
            newTab}, 
lenS = 
            Length[ssSoln]; 
For[kInd = 0; rapSoln = {}, 
            kInd < lenS, kInd++; 
rap[kInd] = 
              Flatten[Map[myRootCheck, 
                  ssSoln[[kInd]]]]; 
newTab = 
              Table[rap[kInd][[k, 2]], {k, 1, 
                  Length[rap[kInd]]}]; 
If[
              VectorNorm[N[newTab]] ≥ 10^(-6), 
              rapSoln = 
                Append[rapSoln, rap[kInd]]];]; 
Return[
            rapSoln];
];

A routine taken from the MathGroup archive:  http://forums.wolfram.com/mathgroup/archive/1995/May/msg00228.html.
          It uses the Plot function to adaptively sample values within a given domain. It is used below to compute the phase-plots.

\(\(\(SetAttributes[AdaptiveSample, HoldAll];\)\(
\)
    \)\)
AdaptiveSample[f_, {var_, min_, max_}, opts___] := 
      Block[{plot, cmp, pts, bend, div, sampler}, 
        sampler = Plot; 
If[MatchQ[Hold[f], Hold[{_, _}]], 
          sampler = ParametricPlot]; 
{cmp, pts, bend, 
            div} = {Compiled, PlotPoints, MaxBend, 
                PlotDivision} /. {opts} /. 
            Options[sampler]; 
plot = 
          sampler[f, {var, min, max}, Compiled -> cmp, 
            PlotPoints -> pts, MaxBend -> bend, 
            PlotDivision -> div, Axes -> None, Frame -> False, 
            DisplayFunction -> Identity, Prolog -> {}, 
            Epilog -> {}]; 
First /@ 
          Cases[First[plot], _Line, Infinity]]

LNA Modules

These modules are the core of the effective stability approximation. 

myRenorm takes as inputs the matrix J0, J1(t), J1(t') and the covariance matrix Ξ as inputs.
It returns the Laplace transform of Jc.

myRenorm[Alist_, A1alist_, A1blist_, Clist_] := 
        Module[{lenA, lenC, sMat, dMat, invSm, G1, Corr1, MatExp, repVar, 
            repMean, intg, 
            laps}, 

lenA = 
            Length[Alist]; 
lenC = 
            Length[Clist]; 
G1 = 
            myExpMat[Alist]; 

Corr1 = 
            G1 . Clist; 
MatExp = 
            Expand[A1alist . G1 . 
                A1blist]; 

repVar = 
            Flatten[Join[
                Table[{nF[Min[i, j]][t1]*nF[Max[i, j]][t2] -> 
                      Corr1[[Min[i, j], Max[i, j]]], nF[Min[i, j]][
                          t2]*nF[Max[i, j]][t1] -> 
                      Corr1[[Max[i, j], Min[i, j]]]}, {i, 1, lenA}, {j, 
                    i + 1, lenA}], 
                Table[nF[i][t1]*nF[i][t2] -> 
                    Corr1[[i, i]], {i, 1, 
                    lenC}]]]; 
repMean = 
            Flatten[Table[{nF[i][t1] -> 0, nF[i][t2] -> 
                    0}, {i, 1, 
                  lenC}]]; 

intg = 
            ReleaseHold[TrigToExp[Expand[MatExp /. repVar /. repMean]] /. 
                  Power -> Matt /. 
                Power -> Matt2];  
Return[
            intg];
];

RenormalizedJacobian takes the propensity vector ν, the stoichiometry matrix S and the steady-state around which the stability is calculated.
The output is the two matrices J0 and the Laplace transform of Jc that will be used in the resolvent equation to calculate the effective eigenvalues.

RenormalizedJacobian[VList_, SMatrix_, SolnList_] := 
        Module[{lenV, varvec, perList, menList, V2, S2, A0, A1, B, Cee, eX, 
            reX, var1, vareq, vsoln, matB, vecC, J1, iInd, cofs, tabby, 
            omSt}, 

lenV = 
            Length[SolnList]; 
varVec = 
            Table[SolnList[[i, 1]], {i, 1, 
                lenV}]; 
perList = 
            Table[varVec[[i]] -> 
                varVec[[i]] + L*nF[i][t], {i, 1, 
                lenV}]; 
menList = 
            Table[nF[i][t] -> 0, {i, 1, 
                lenV}]; 

S2 = 
            ExpandAll[SMatrix /. perList]; 
V2 = 
            ExpandAll[
              VList /. 
                perList]; 

A0 = Table[
                  Flatten[
                    Table[D[(S2 . V2)[[i]], varVec[[j]]], {j, 1, 
                        lenV}]], {i, 1, lenV}] /. menList /. 
              SolnList; 
A1[x_] = 
            Chop[Normal[
                      Series[Expand[
                          Table[Flatten[
                                Table[
                                  D[(S2 . V2)[[i]], 
                                    varVec[[j]]], {j, 1, lenV}]], {i, 1, 
                                lenV}] - A0], {L, 0, 1}]] /. L -> 1 /. 
                  SolnList] /. 
              t -> 
                x; 
B = (S2 . 
                    DiagonalMatrix[Flatten[V2]] . Transpose[S2]) /. 
                menList /. 
              SolnList; 

eX = 
            Table[X[Min[i, j], Max[i, j]], {j, 1, lenV}, {i, 1, 
                lenV}]; 
reX = 
            Flatten[Table[
                X[Min[i, j], Max[i, j]], {j, 1, lenV}, {i, j, 
                  lenV}]]; 

var1 = 
            Expand[((A0 . eX) + (eX . Transpose[A0]) + 
                  B)]; 
vareq = 
            Flatten[Table[
                0 == var1[[i, j]], {i, 1, lenV}, {j, i, 
                  lenV}]]; 
{matB, vecC} = 
            LinearEquationsToMatrices[vareq, 
              reX]; 
vsolna = 
            Inverse[matB] . vecC; 
vsoln = 
            Table[reX[[i]] -> vsolna[[i]], {i, 1, 
                Length[reX]}]; 
Cee = 
            Expand[eX /. vsoln]; 

J1 = 
            myRenorm[A0, A1[t1], A1[t2], Cee]; 
Return[{A0, 
              J1}];
];

Hurwitz Stability Modules

These are modules used below to invoke the Routh-Hurwitz criterion to establish the stability of the effective eigenvalues.

FindLastNonzero[a_] := 
        Module[{i, n = Length[a]}, 
For[i = n, 
            i > 0, i--, 
            If[a[[i]] == 0, Continue[], 
              Return[i], 
Return[i]]]; 
          Return[0]
];

HurwitzDeterminant[b_, k_] := 
        Module[{n = Length[b] - 1, i, j, m, A}, 
If[
            k < 1 || k > n || n ≤ 0,  Return[0]]; 
          If[k == n, Return[b[[n + 1]]]]; 
If[
            k == 1, Return[b[[2]]]];  
          A = Table[0, {k}, {k}]; 
For[i = 1, 
            i ≤ k, i++, 
For[j = 1, 
              j ≤ k, j++, 
              m = 2*j - i + 1; 
If[
                m > 0 && m ≤ n + 1, 
                A[[i, j]] = b[[m]]  ]]]; 
Return[
            Det[A]]
];

RouthHurwitzStability[b_] := 
        Module[{n = Length[b] - 1, i, c}, 
n = 
            FindLastNonzero[b] - 1;  
          If[n < 3, Return[b]]; 
c = Table[0, {n + 1}];  
          c[[1]] = b[[1]]; 
For[i = 2, 
            i ≤ n + 1, i++, 
            c[[i]] = 
              HurwitzDeterminant[b, i - 1]]; 
Return[
            c]
];

LienardChipartStability[b_] := 
        Module[{n = Length[b] - 1, i, c}, 
If[n < 3, 
            Return[b]]; 
If[EvenQ[n], 
            c = Table[{b[[i]], HurwitzDeterminant[b, i]}, {i, 1, n, 
                  2}]]; 
If[OddQ[n], 
            c = Table[{b[[i + 1]], HurwitzDeterminant[b, i + 1]}, {i, 1, 
                  n - 1, 2}]; 
c = {b[[1]], c}]; 
          AppendTo[c, b[[n + 1]]]; 
          c = Flatten[c]; 
Return[
            c]
];

Stability Modules

Because the resolvent equation is such a devil to solve, I have constructed a variety of routines to compute the effective eigenvalues. 
    
    pertEig uses the expansion of the resolvent equation discussed in the text.
    pertEigExact simply calculates the roots of the resolvent equation directly.
    pertEigExplicit attempts an analytic solution for the roots.
    pertEigHurwitz uses the Routh-Hurwitz criterion to determine whether the system is stable or not, without calculating the roots directly.
    
    myStable___ take as inputs the propensity vector ν, the stoichiometry matrix S, the steady-state around which the stability is calculated and the value of the system size (Vcell in the write-up).
    The output is a list of the effective eignvalues found by whatever method is specified.

pertEig[J0_, J1_, omParm_] := 
        Module[{nJ = Length[J0], LVar, VVar, pVar, p1Var, corrTab, res, 
            corrS}, 
{LVar, VVar} = 
            Eigensystem[J0]; 
pVar = 
            Transpose[VVar]; 
p1Var = 
            Inverse[pVar]; 

corrTab = 
            Expand[p1Var . J1 . pVar]; 
corrS = 
            Chop[Table[
                Limit[corrTab[[i, i]], s -> LVar[[i]]], {i, 1, 
                  nJ}]]; 

res = 
            Table[LVar[[i]] + corrS[[i]]/omParm, {i, 1, 
                nJ}]; 
Return[
            res];
];

pertEigExact[J0_, J1_, eps_] := 
        Module[{nJ = Length[J0], res, eq}, 
eq = 
            Numerator[
              Together[
                Det[s*IdentityMatrix[nJ] - J0 - 
                    eps^2*J1]]]; 

res = 
            Map[myRootGet, 
              NSolve[eq == 0, s]]; 
Return[
            res];
];

pertEigExplicit[J0_, J1_, eps_] := 
        Module[{nJ = Length[J0], res, eq}, 
eq = 
            Numerator[
              Together[
                Det[s*IdentityMatrix[nJ] - J0 - 
                    eps^2*J1]]]; 

res = 
            Solve[eq == 0, s]; 
Return[
            res];
];

myStable[VList_, SMatrix_, SolnList_, omParm_] := 
        Module[{J0Var, J1Var, eigRes}, 
{J0Var, J1Var} = 
            RenormalizedJacobian[VList, SMatrix, 
              SolnList]; 
eigRes = 
            Max[Re[pertEig[J0Var, J1Var, omParm]]]; 
Return[
            eigRes];
];

myStableExplicit[VList_, SMatrix_, SolnList_, omParm_] := 
        Module[{J0Var, J1Var, eigRes}, 
{J0Var, J1Var} = 
            RenormalizedJacobian[VList, SMatrix, 
              SolnList]; 
Return[{J0Var, 
              J1Var}];
];

myStableOmega[VList_, SMatrix_, SolnList_, omParm_] := 
        Module[{J0Var, J1Var, eigRes}, 
{J0Var, J1Var} = 
            RenormalizedJacobian[VList, SMatrix, 
              SolnList]; 
eigRes = 
            pertEig[J0Var, J1Var, omParm]; 
Return[
            eigRes];
];

myStableOmegaExact[VList_, SMatrix_, SolnList_, omParm_] := 
        Module[{J0Var, J1Var, eigRes}, 
{J0Var, J1Var} = 
            RenormalizedJacobian[VList, SMatrix, 
              SolnList]; 
eigRes = 
            pertEigExact[J0Var, J1Var, 
              1/Sqrt[omParm]]; 
Return[
            eigRes];
];

myStableExact[VList_, SMatrix_, SolnList_, omParm_] := 
        Module[{J0Var, J1Var, eigRes}, 
{J0Var, J1Var} = 
            RenormalizedJacobian[VList, SMatrix, 
              SolnList]; 
eigRes = 
            Max[Re[Evaluate[
                  pertEigExact[J0Var, J1Var, 
                    1/Sqrt[omParm]]]]]; 
Return[
            eigRes];
];

pertEigHurwitz[J0VarH_, J1VarH_, omParmH_] := 
        Module[{nJSize = Length[J0VarH], lapdenom, cofs, 
            tabby}, 
lapdenom = 
            Det[Expand[
                s*IdentityMatrix[nJSize] - J0VarH - 
                  J1VarH/omParmH]]; 
cofs = 
            CoefficientList[
              Numerator[Together[Chop[lapdenom, 10^(-6)]]], 
              s]; 
tabby = 
            Min[Re[LienardChipartStability[
                  Re[Chop[cofs/cofs[[1]], 
                      10^(-6)]]]]]; 
Return[-tabby];
];

myStableHurwitz[VList_, SMatrix_, SolnList_, omParm_] := 
        Module[{J0Var, J1Var, eigRes}, 
{J0Var, J1Var} = 
            RenormalizedJacobian[VList, SMatrix, 
              SolnList]; 
eigRes = 
            Max[Re[Evaluate[
                  pertEigHurwitz[J0Var, J1Var, 
                    omParm]]]]; 
Return[
            eigRes];
];

Phase Diagram Modules

This module puts all the preceding ingredients together to compute a phase-plot with the phase boundaries corrected by the effective stability approximation. 
      The input is the propensity vector and stoichiometry matrix {ν,S}, a list of the state variables, a list of the two phase plot axes {xPar,yPar}, the domain and range of the phase plot. (It is a good idea not to set either xMin or yMin to 0 since the order of the system may change and the numerics get messed up; instead set xMin or yMin to some small number, eg. 0.001). 
      Finally, there is a flag to allow options to be set. Options are:
      
      reolutioin in the vertical direction: delY->50 (Default; can be any integer)
      system size (Vcell): Omega->1 (Default; can be any real number)
      whether the phase-plot countours are shaded: ContourShading->True (Default; can be True or False)
      countour threshold for the number of roots:Contours->{1/2,3/2,5/2,7/2}(Default; can be any list of real numbers)
      whether or not to use the effective stability analysis:Noise->False (Default; set to True to use the ESA)
      whether or not to use the exact numerical solution of the resolvent equation, or simply the perturbation approximation:ExactPolynomial->False(Default; set to True to use exact solution)

\(\(phasePlotNoise[{sMatrix_, vVector_}, 
          listVars_, {xPar_, yPar_}, {{xMin_, xMax_}, {yMin_, yMax_}}, 
          opts_] := Module[{nVariables = Length[listVars], listRates, 
            nReactions, numberY, omegaPar, shadeCont, testChoice, 
            noiseChoice, hurwitzChoice, exactChoice, EquilibEqs, contList, 
            Jac, deltaY, stabEps, xInd, yInd, resInter, totalData, 
            inter}, 

listRates = 
            Flatten[sMatrix . vVector]; 
          nReactions = Length[listRates]; 
If[
            nVariables ≠ nReactions, 
            Print["Rate vector and state vector must have the same dimension!"]; 
            Abort[]]; 

numberY = \(delY /. 
                opts\) /. {delY -> 
                  50}; 
omegaPar = \(Omega /. 
                opts\) /. {Omega -> 
                  1}; 
shadeCont = \(ContourShading /. 
                opts\) /. {ContourShading -> 
                  True}; 
dipFun = \(DisplayFunction /. 
                opts\) /. {DisplayFunction -> $DisplayFunction}; 
contList = \(Contours /. 
                opts\) /. {Contours -> {1/2, 3/2, 5/2, 
                    7/2}}; 
testChoice = \(Bistable /. 
                opts\) /. {Bistable -> 
                  False}; 
noiseChoice = \(Noise /. 
                opts\) /. {Noise -> 
                  False}; 
hurwitzChoice = \(Hurwitz /. 
                opts\) /. {Hurwitz -> 
                  False}; 
exactChoice = \(ExactPolynomial /. 
                opts\) /. {ExactPolynomial -> 
                  False}; 

If[
            testChoice || hurwitzChoice || exactChoice, 
            noiseChoice = 
              True]; 

EquilibEqs = 
            Table[listRates[\([i]\)] == 0, {i, 1, 
                nReactions}]; 
Jac = 
            Table[D[listRates[\([i]\)], listVars[\([j]\)]], {i, 1, 
                nReactions}, {j, 1, 
                nReactions}]; 
deltaY = \((yMax - yMin)\)/
              numberY; 
stabEps = 
            0; 

Which[
            Not[noiseChoice], 
testFunc[EquilibEqsVar_, 
                JacVar_, 
                listVarsVar_, {vVectorVar_, sMatrixVar_}, {xParameter_, 
                  yParameter_}, omegaParVar_] := 
              Module[{soln, iInd}, 
soln = 
                  myRootExtraction[
                    NSolve[EquilibEqsVar /. {xPar -> xParameter, 
                          yPar -> yParameter}, 
                      listVarsVar]]; 
Return[
                  Apply[Plus, 
                    Table[
If[
                        Max[Re[
                              Eigenvalues[\(JacVar /. {xPar -> 
                                        xParameter, 
                                      yPar -> yParameter}\) /. 
                                  soln[\([iInd]\)]]]] < 0, 1, 0], {iInd, 1, 
                        Length[
                          soln]}]]]
], 
testChoice, 
testFunc[EquilibEqsVar_, 
                JacVar_, 
                listVarsVar_, {vVectorVar_, sMatrixVar_}, {xParameter_, 
                  yParameter_}, omegaParVar_] := 
              Module[{soln, iInd, res1}, 
soln = 
                  myRootExtraction[
                    NSolve[EquilibEqsVar /. {xPar -> xParameter, 
                          yPar -> yParameter}, 
                      listVarsVar]]; 
res1 = 
                  Apply[Plus, 
                    Table[
If[
                        Max[Re[
                              Eigenvalues[\(JacVar /. {xPar -> 
                                        xParameter, 
                                      yPar -> yParameter}\) /. 
                                  soln[\([iInd]\)]]]] < 0, 1, 0], {iInd, 1, 
                        Length[soln]}]]; 
If[
                  res1 > 1, 
Return[
                    Apply[Plus, 
                      Table[
If[
                          myStable[
                              vVectorVar /. {xPar -> xParameter, 
                                  yPar -> yParameter}, 
                              sMatrixVar /. {xPar -> xParameter, 
                                  yPar -> yParameter}, soln[\([iInd]\)], 
                              omegaParVar] < 0, 1, 0], {iInd, 1, 
                          Length[soln]}]]], 
Return[
                    res1]
]
], 
hurwitzChoice, 
testFunc[EquilibEqsVar_, 
                JacVar_, 
                listVarsVar_, {vVectorVar_, sMatrixVar_}, {xParameter_, 
                  yParameter_}, omegaParVar_] := 
              Module[{soln, iInd}, 
soln = 
                  myRootExtraction[
                    NSolve[EquilibEqsVar /. {xPar -> xParameter, 
                          yPar -> yParameter}, 
                      listVarsVar]]; 
Return[
Apply[
                    Plus, 
Table[
If[
                        Max[Re[
                              Eigenvalues[\(JacVar /. {xPar -> 
                                        xParameter, 
                                      yPar -> yParameter}\) /. 
                                  soln[\([iInd]\)]]]] < 
                          0, 
If[
myStableHurwitz[vVectorVar /. {xPar -> xParameter, 
                                  yPar -> yParameter}, 
                              sMatrixVar /. {xPar -> xParameter, 
                                  yPar -> yParameter}, soln[\([iInd]\)], 
                              omegaParVar] ≤ 0, 1, 0], 0], {iInd, 
                        1, Length[
                          soln]}]]]
], 
exactChoice, 
testFunc[EquilibEqsVar_, 
                JacVar_, 
                listVarsVar_, {vVectorVar_, sMatrixVar_}, {xParameter_, 
                  yParameter_}, omegaParVar_] := 
              Module[{soln, iInd}, 
soln = 
                  myRootExtraction[
                    NSolve[EquilibEqsVar /. {xPar -> xParameter, 
                          yPar -> yParameter}, 
                      listVarsVar]]; 
Return[
                  Apply[Plus, 
                    Table[
If[
                        Max[Re[
                              Eigenvalues[\(JacVar /. {xPar -> 
                                        xParameter, 
                                      yPar -> yParameter}\) /. 
                                  soln[\([iInd]\)]]]] < 0, 
                        If[Chop[
                              myStableExact[
                                vVectorVar /. {xPar -> xParameter, 
                                    yPar -> yParameter}, 
                                sMatrixVar /. {xPar -> xParameter, 
                                    yPar -> yParameter}, 
                                soln[\([iInd]\)], omegaParVar]] ≤ 
                            0, 1, 0], 0], {iInd, 1, 
                        Length[
                          soln]}]]]
], 
True, 
testFunc[EquilibEqsVar_, JacVar_, 
                listVarsVar_, {vVectorVar_, sMatrixVar_}, {xParameter_, 
                  yParameter_}, omegaParVar_] := 
              Module[{soln, iInd}, 
soln = 
                  myRootExtraction[
                    NSolve[EquilibEqsVar /. {xPar -> xParameter, 
                          yPar -> yParameter}, 
                      listVarsVar]]; 
Return[
                  Apply[Plus, 
                    Table[
If[
                        Max[Re[
                              Eigenvalues[\(JacVar /. {xPar -> 
                                        xParameter, 
                                      yPar -> yParameter}\) /. 
                                  soln[\([iInd]\)]]]] < 0, 
                        If[Chop[
                              myStable[
                                vVectorVar /. {xPar -> xParameter, 
                                    yPar -> yParameter}, 
                                sMatrixVar /. {xPar -> xParameter, 
                                    yPar -> yParameter}, 
                                soln[\([iInd]\)], omegaParVar]] ≤ 
                            0, 1, 0], 0], {iInd, 1, 
                        Length[
                          
                          soln]}]]]
]]; 

For[yInd = yMin; dX = \((xMax - xMin)\)/\((10*numberY)\); 
            totalData = {}; hiData = {}; loData = {}, yInd ≤ yMax, 
            yInd += deltaY, 
resInter = 
              Interpolation[
                Flatten[
                  AdaptiveSample[
                    testFunc[EquilibEqs, Jac, 
                      listVars, {vVector, sMatrix}, {xInd, yInd}, 
                      omegaPar], \ {xInd, xMin, xMax}, \ 
                    PlotPoints\  -> 30\ ], 1], 
                InterpolationOrder -> 1]; 
totalData = 
              Append[totalData, 
                Table[resInter[xInd], {xInd, xMin, 
                    xMax, \((xMax - xMin)\)/\((5*
                          numberY)\)}]]; 
If[
              testChoice, 
c = xMin; 
              While[\((resInter[c] < 2)\) && \((c < xMax)\), c += dX]; 
              If[c < xMax && c ≠ xMin, 
                loData = Append[loData, {c, yInd}]]; 
c = 
                xMax; While[\((resInter[c] < 2)\) && \((c > xMin)\), 
                c -= dX]; 
              If[c > xMin && c ≠ xMax, 
                hiData = 
                  Append[hiData, {c, 
                      yInd}]];];
]; 
If[
            testChoice, 
            Return[Show[Graphics[{Line[loData], Line[hiData]}], 
                PlotRange -> {{xMin, xMax}, {yMin, yMax}}, 
                Axes -> True, AspectRatio -> 1, 
                DisplayFunction -> dipFun]], 
inter = 
              ListInterpolation[
                Transpose[totalData], {{xMin, xMax}, {yMin, yMax}}, 
                InterpolationOrder -> 1]; 
Return[
              ContourPlot[inter[x, y], {x, xMin, xMax}, {y, yMin, yMax}, 
                Contours -> contList, 
                ColorFunction -> \((GrayLevel[1 - #] &)\), 
                PlotPoints -> numberY, PlotRange -> All, 
                ContourShading -> shadeCont, 
                DisplayFunction -> 
                  dipFun]]];

];\)\)

\(\(phasePlotOmega[{sMatrix_, vVector_}, listVars_, 
          yPar_, {yMin_, yMax_}, {OmMin_, OmMax_}, opts_] := 
        Module[{nVariables = Length[listVars], listRates, nReactions, 
            numberY, omegaPar, shadeCont, testChoice, noiseChoice, 
            hurwitzChoice, exactChoice, EquilibEqs, contList, Jac, deltaY, 
            stabEps, xInd, yInd, resInter, totalData, inter, 
            logOmega}, 

listRates = 
            Flatten[sMatrix . vVector]; 
          nReactions = Length[listRates]; 
If[
            nVariables ≠ nReactions, 
            Print["Rate vector and state vector must have the same dimension!"]; 
            Abort[]]; 

numberY = \(delY /. 
                opts\) /. {delY -> 
                  50}; 
shadeCont = \(ContourShading /. 
                opts\) /. {ContourShading -> 
                  True}; 
dipFun = \(DisplayFunction /. 
                opts\) /. {DisplayFunction -> $DisplayFunction}; 
contList = \(Contours /. 
                opts\) /. {Contours -> {1/2, 3/2, 5/2, 
                    7/2}}; 
testChoice = \(Bistable /. 
                opts\) /. {Bistable -> 
                  False}; 
noiseChoice = \(Noise /. 
                opts\) /. {Noise -> 
                  False}; 
hurwitzChoice = \(Hurwitz /. 
                opts\) /. {Hurwitz -> 
                  False}; 
exactChoice = \(ExactPolynomial /. 
                opts\) /. {ExactPolynomial -> 
                  False}; 
logOmega = \(LogOmega /. 
                opts\) /. {LogOmega -> 
                  True}; 

EquilibEqs = 
            Table[listRates[\([i]\)] == 0, {i, 1, 
                nReactions}]; 
Jac = 
            Table[D[listRates[\([i]\)], listVars[\([j]\)]], {i, 1, 
                nReactions}, {j, 1, 
                nReactions}]; 
deltaY = \((yMax - yMin)\)/
              numberY; 
stabEps = 
            0; 

For[yInd = yMin; 
            dX = \((OmMax - OmMin)\)/\((10*numberY)\); totalData = {}; 
            hiData = {}; loData = {}, yInd ≤ yMax, 
            yInd += deltaY, 
soln = 
              myRootExtraction[
                NSolve[EquilibEqs /. {yPar -> yInd}, 
                  listVars]]; 
If[
              exactChoice, 
resInter = 
                Interpolation[
                  Flatten[
                    AdaptiveSample[
                      Max[Re[myStableOmegaExact[
                            vVector /. {yPar -> yInd}, 
                            sMatrix /. {yPar -> yInd}, soln[\([1]\)], 
                            If[logOmega, 10^xInd, xInd]]]], \ {xInd, OmMin, 
                        OmMax}, \ PlotPoints\  -> 30\ ], 1], 
                  InterpolationOrder -> 
                    1], 

stabList = 
                myStableOmega[vVector /. {yPar -> yInd}, 
                  sMatrix /. {yPar -> yInd}, soln[\([1]\)], 
                  If[logOmega, 10^xInd, xInd]]; 
resInter = 
                Interpolation[
                  Flatten[
                    AdaptiveSample[
                      Max[Re[stabList]], \ {xInd, OmMin, OmMax}, \ 
                      PlotPoints\  -> 30\ ], 1], 
                  InterpolationOrder -> 
                    1]]; 

totalData = 
              Append[totalData, 
                Table[resInter[xInd], {xInd, OmMin, 
                    OmMax, \((OmMax - OmMin)\)/\((5*
                          numberY)\)}]];
]; 

inter = 
            ListInterpolation[
              Transpose[totalData], {{OmMin, OmMax}, {yMin, yMax}}, 
              InterpolationOrder -> 1]; 
Return[
            ContourPlot[inter[x, y], {x, OmMin, OmMax}, {y, yMin, yMax}, 
              Contours -> {0}, 
              ColorFunction -> \((GrayLevel[1 - #] &)\), 
              PlotPoints -> numberY, PlotRange -> All, 
              ContourShading -> shadeCont, 
              DisplayFunction -> 
                dipFun]];

];\)\)

Transcriptional Autoactivator

Clear[Ka, bC, bBurst];

solnB is the discreteness parameter in activator synthesis/degradation

solnB = 
        Flatten[Solve[(bBurst + 1)/2/Ka == bC, 
            bBurst]];

gA is the promoter activity function.

gA[
          a_] = (1/w + (a/Ka)^2)/(1 + (a/Ka)^2);

vVec is the propensity vector and sMat is the stoichiometry matrix

vVec = 
        Transpose[{{a0/bBurst*gA[a], a}}];
sMat = {{bBurst, -1}};

The nominal rate constants

Ka = 25;

ssvar = 
        Join[{a0 -> Ka*10^wa, w -> 10^wf}, solnB];

Clear[bC]

bC = 0;
Timing[
        pAct1N0 = 
          phasePlotNoise[{sMat, vVec} /. ssvar, {a}, {wf, 
              wa}, {{1, 2}, {0.2, 0.7}}, {delY -> 100, Omega -> 1, 
              ContourShading -> False, Noise -> True, 
              DisplayFunction -> Identity, 
              Contours -> {3/2}}]];

bC = 0.1;
Timing[
        pAct1N01 = 
          phasePlotNoise[{sMat, vVec} /. ssvar, {a}, {wf, 
              wa}, {{1, 2}, {0.2, 0.7}}, {delY -> 100, Omega -> 1, 
              ContourShading -> False, Noise -> True, 
              DisplayFunction -> Identity, 
              Contours -> {3/2}}]];

bC = 0.2;
Timing[
        pAct1N02 = 
          phasePlotNoise[{sMat, vVec} /. ssvar, {a}, {wf, 
              wa}, {{1, 2}, {0.2, 0.7}}, {delY -> 100, Omega -> 1, 
              ContourShading -> False, Noise -> True, 
              DisplayFunction -> Identity, 
              Contours -> {3/2}}]];

bC = 0.3;
Timing[
        pAct1N03 = 
          phasePlotNoise[{sMat, vVec} /. ssvar, {a}, {wf, 
              wa}, {{1, 2}, {0.2, 0.7}}, {delY -> 100, Omega -> 1, 
              ContourShading -> False, Noise -> True, 
              DisplayFunction -> Identity, 
              Contours -> {3/2}}]];

Show[pAct1N0, pAct1N01, pAct1N02, pAct1N03, 
      DisplayFunction -> $DisplayFunction, 
      PlotRange -> {{1, 2}, {0, 1}}]

Compare with figure 2a in the PNAS paper.

Clear[gA]

vVec = 
        Transpose[{{a0/bBurst*gA[a], a}}];
sMat = {{bBurst, -1}};

To compute the effective eigenvalue as shown in Eq. 27 of the Supplementary Material:

FullSimplify[
      Evaluate[myStable[vVec, sMat, {a -> As}, Omega]]]

Toggle Switch

Clear[Ka, bC, bBurst, gA];

The toggle switch example corresponds to the example detailed in the tutorial notes write-up. The Bistable->True option can be used in systems that are strictly bistable/monostable in order to smooth the phase boundary.

gr[wConst_, kConst_, 
        xConst_] := (((1 + ((xConst (2 + (xConst/(2 kConst))))/(2 kConst wConst))))/((((1 + (xConst/(2 kConst))))^2) ))

vVec is the propensity vector and sMat is the stoichiometry matrix

vVec = 
        Transpose[{{r10/bBurst*gr[w2, kr2, r2], 
              r1, 
r20/bBurst*gr[w1, kr1, r1], 
              r2}}];
sMat = {
{bBurst, -1, 0, 
            0}, 
{0, 0, bBurst, -1}};

The nominal rate constants

ssvar = {w1 -> 100, w2 -> 100, kr1 -> 5, 
          kr2 -> 5, bBurst -> 10};

Deterministic phase-plot for the symmetric toggle switch.

Timing[
      pTog100N = 
        phasePlotNoise[{sMat, vVec} /. ssvar, {r1, r2}, {r10, 
            r20}, {{1, 250}, {1, 250}}, {delY -> 50, Omega -> 100, 
            ContourShading -> False, Contours -> {3/2}, 
            DisplayFunction -> Identity, Noise -> True, 
            Bistable -> True}]]

Phase-plot for the symmetric toggle switch computed using the ESA.

Timing[
      pTog1N = phasePlotNoise[{sMat, vVec} /. ssvar, {r1, r2}, {r10, 
            r20}, {{100, 250}, {100, 250}}, {delY -> 50, 
            Omega -> 1, ContourShading -> False, 
            Noise -> True, Bistable -> True, 
            Contours -> {3/2}, 
            DisplayFunction -> Identity}]]

Show[Show[pTog1N, DisplayFunction -> Identity], pTog100N, 
      DisplayFunction -> $DisplayFunction, 
      PlotRange -> {{0, 250}, {0, 250}}]

Clear[gr]

vVec = 
        Transpose[{{r10/bBurst*gr[w2, kr2, r2], 
              r1, 
r20/bBurst*gr[w1, kr1, r1], 
              r2}}];
sMat = {
{bBurst, -1, 0, 
            0}, 
{0, 0, bBurst, -1}};

In this case, an analytic expression for the effective eigenvalues is not very useful.

Simplify[
      Evaluate[myStable[vVec, sMat, {r1 -> R1s, r2 -> R2s}, 
          Omega]]]

Vilar oscillator - reduced model

In this case, the approximation to the solution of the resolvent equation breaks-down since there are regions of the parameter space where the eigenvalues become close to degenerate. Consequently, the ExactPolynomial->True option is used to ensure the effective eigenvalues are computed numerically.

Clear[Ka, bC, bBurst];
Clear[a, r, c, gFA, gFR];

Vcell = 100;
solnB = 
        Flatten[Solve[(burstA + 1)/2/Ka/Vcell == bC, 
            burstA]];

vVec is the propensity vector and sMat is the stoichiometry matrix

vVec = 
        Transpose[{{gamA/burstA*gFA[a], dA*a, gamC*a*r, gamR/burstR*gFR[a], 
              dR*r, dA*c}}];

sMat = {
{burstA, -1, -1, 0, 0, 
            0}, 
{0, 0, -1, burstR, -1, 
            1}, 
{0, 0, 1, 0, 0, -1}};

gFA and GFR are the promoter activity functions.

gFA[
          A_] = (fA^(-1) + (A/Ka))/(1 + (A/
                  Ka));
gFR[A_] = ((A/Kr))/(1 + (A/Kr));

The nominal rate parameters

ssVars = 
        Join[{gamA -> 25, gamR -> 5, Ka -> 0.5, Kr -> 1, 
            dA -> 1, gamC -> 200, fA -> 10, 
            a0 -> gamA/dA}, solnB];

Determinisitic phase boundary.

Timing[
      pVilarRed0 = 
        phasePlotNoise[{sMat, vVec} /. burstR -> 10 //. ssVars, {a, 
            r, c}, {bC, dR}, {{0, 0.1}, {0.01, 0.12}}, {delY -> 100, 
            ContourShading -> False, Noise -> False, 
            DisplayFunction -> Identity}]]

Phase boundary determined by the ESA.

Timing[
      pVilarRed10 = 
        phasePlotNoise[{sMat, vVec} /. burstR -> 10 //. ssVars, {a, 
            r, c}, {bC, dR}, {{0, 0.1}, {0.01, 0.12}}, {delY -> 100, 
            Omega -> Vcell, ContourShading -> False, 
            ExactPolynomial -> True, Noise -> True, 
            DisplayFunction -> Identity}]]

datPlot is a plot of the stochastic simulation results corresponding to a noise-to-signal ratio of 0.2.

datPlot = 
        ListPlot[{{0.02, 0.065}, {0.03, 0.062}, {0.04, 0.054}, {0.05, 
              0.05}, {0.06, 0.048}, {0.07, 0.047}, {0.08, 0.046}, {0.09, 
              0.046}, {0.1, 0.045}}, PlotStyle -> {PointSize[0.02]}, 
          DisplayFunction -> Identity];

Show[pVilarRed10, pVilarRed0, datPlot, 
      PlotRange -> {{0, 0.1}, {0, 0.15}}, 
      DisplayFunction -> $DisplayFunction]

Compare with figure 3a in the PNAS paper. (Note that the x-axis in the PNAS article is mislabelled - the final entry should read '10', not '1').

TimeUsed[]/60