Reakciókinetika
Sipos Szabó Eszter

Csomagok

AppendTo[$Path, "C:\Program Files\Wolfram \
Research\Mathematica\5.2\AddOns\Applications\ReactionKinetics\R\
eactionKinetics"];

Get["ReactionKinetics`ReactionKinetics`"]

AppendTo[$Path, "C:\Program Files\Wolfram \
Research\Mathematica\5.2\AddOns\Applications\ReplaceVariables"];

Get["ReplaceVariables`"]

Előkészületek

Az alábbi reakciót nézzük: S+M↔ES→P+M
A ReactionKinetics csomag segítségével felírjuk, majd megoldjuk a differenciálegyenleteket:

\(r1 = 
      ToStoichiometry[{"S" + "M" ↔ "ES" → "P" + "M"}]\)

\(rhs = RightHandSide[r1, Array[k\_# &, 3]]\)

lhs = variables /. r1 /. {(Subscript[c, d])_ -> (Subscript[c, d])'[t]}

r4 = rhs /. {(Subscript[c, d])_ -> (Subscript[c, d])[t]}

r5 = Thread[Equal[lhs, r4]]

r6 = 
      Thread[Equal[variables /. r1 /. {(Subscript[c, d])_ -> (Subscript[c, d])[0]}, {0, (Subscript[M, 0]), 
            0, (Subscript[S, 0])}]]

A diffegyenlet-rendszer megoldása:

nds = 
      NDSolve[Union[{r5, r6}] /. {(Subscript[M, 0]) -> 0.01, (Subscript[S, 0]) -> 0.1, 
            (Subscript[k, 1]) -> 100, (Subscript[k, 2]) -> 0.01, 
            (Subscript[k, 3]) -> 1}, variables /. 
            r1 /. {(Subscript[c, d])_ -> (Subscript[c, d])[t]}, {t, 0, 5}]

Plot[
      Evaluate[{(Subscript[c, "P"])[t], (Subscript[c, "S"])[t], (Subscript[c, "ES"])[t]} /. nds], {t, 0, 
        5}, PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1], 
          RGBColor[0, 1, 0]}]

Egy Module-ban összefoglalva az eddigieket:

ClearAll[kindemo]; 
    kindemo[mech_, rrates_, inival_, tmin_, tmax_] := 
      Module[{r1, de}, r1 = ToStoichiometry[mech]; Print[r1]; 
        de = Thread[
            Equal[variables /. r1 /. {(Subscript[c, d])_ -> (Subscript[c, d])'[t]}, 
              RightHandSide[r1, 
                  rrates] /. {(Subscript[c, d])_ -> 
                    (Subscript[c, d])[t]}]]; 
ini = 
          Thread[Equal[variables /. r1 /. {(Subscript[c, d])_ -> (Subscript[c, d])[0]}, 
              inival]]; 
NDSolve[Join[de, ini], 
            variables /. r1, {t, tmin, 
              tmax}][[1]]]

\(ered = 
      kindemo[{"S" + "M" ↔ "ES" → "P" + "M"}, {100, 0.01, 1}, {0, 0.01, 0, 0.1}, 0, 5]\)

Plot[Evaluate[((Subscript[c, "ES"]) /. ered)[t]], {t, 0, 5}]

A csomag Concentrations parancsával is megkaphatjuk a reakciókra felírt diffegyenletek megoldását:

vl = 
      ToStoichiometry[{"X" → 
            2 "X", "X" + "Y" → 
            2 "Y", "Y" → "A"}]

? Concentrations

ere = Concentrations[vl, {1, 1, 1}, {0.5, 2, 1}, 0, 15]

Plot[Evaluate[((Subscript[c, "X"])[t] /. ere)], {t, 0, 15}]

Egy bonyolult modell

Options[Concentrations]

Needs["Utilities`Notation`"];
Needs["Utilities`FilterOptions`"];
Needs["Graphics`"];
Needs["Algebra`InequalitySolve`"]

A reakcióegyenleteket a Micaelis-Menten kinetika alapján írjuk fel. Ezek különböznek az előző csomag által felírt egyenletektől, így kézzel írtam be a diffegyenleteket:

R5={((Subscript[c, "MPF"]))'[t]==(Subscript[k, 1])-((Subscript[k, 2]) + (Subscript[k, (Subscript[2, 1])]) (Subscript[c, "CDC20a"])[t] + 
                    (Subscript[k, (Subscript[2, 2])]) (Subscript[c, "CDH1a"])[t])(Subscript[c, "MPF"])[t],
((Subscript[c, "CDC20a"]))'[t]==-(Subscript[k, 6]) (Subscript[c, "CDC20a"])[
                    t]/((Subscript[J, 6]) +  (Subscript[c, "CDC20a"])[t]) + ((Subscript[k, 5])  + 
                    (Subscript[k, (Subscript[5, 1])]) (Subscript[c, "APCP"])[t])*( 
                    (Subscript[t, CDC20]) -  (Subscript[c, "CDC20a"])[t])/((Subscript[J, 5]) + (Subscript[t, CDC20]) -  
                      (Subscript[c, "CDC20a"])[t]),
((Subscript[c, "APCP"]))'[t]==-(Subscript[k, 3]) (Subscript[c, "APCP"])[
                    t]/((Subscript[J, 3]) +  (Subscript[c, "APCP"])[t]) + ((Subscript[k, 4])  + 
                    (Subscript[k, (Subscript[4, 1])]) (Subscript[c, "MPF"])[t])*( 
                    (Subscript[t, APC]) -  (Subscript[c, "APCP"])[t])/((Subscript[J, 4]) + (Subscript[t, APC]) -  
                      (Subscript[c, "APCP"])[t]),
((Subscript[c, "CDH1a"]))'[t]==-( 
                    (Subscript[k, 7]) + (Subscript[k, (Subscript[7, 1])]) (Subscript[c, "MPF"])[t]) (Subscript[c, "CDH1a"])[
                    t]/((Subscript[J, 7]) +  (Subscript[c, "CDH1a"])[t]) + 
              (Subscript[k, 8]) ( 
                    (Subscript[t, CDH1]) -  (Subscript[c, "CDH1a"])[t])/((Subscript[J, 8]) + (Subscript[t, CDH1]) -  
                      (Subscript[c, "CDH1a"])[t])
}
R6 = {(Subscript[c, "MPF"])[0] == 
          (Subscript[M, 0]), (Subscript[c, "CDC20a"])[0] == 0, (Subscript[c, "APCP"])[0] == 0, 
        (Subscript[c, "CDH1a"])[0] == (Subscript[C, 0])}

Az alábbi paraméterhalmazra oldjuk meg a diffegyenletrendszert!

NDSo = 
      NDSolve[Union[{R5, R6}] /. {(Subscript[M, 0]) -> 0, (Subscript[C, 0]) -> 0, 
            (Subscript[k, 1]) -> 1, (Subscript[k, 2]) -> 0.1, (Subscript[k, (Subscript[2, 1])]) -> 1, 
            (Subscript[k, (Subscript[2, 2])]) -> 3, (Subscript[k, 3]) -> 0.8, (Subscript[k, 4]) -> 0.1, 
            (Subscript[k, (Subscript[4, 1])]) -> 1, (Subscript[k, 5]) -> 0.1, (Subscript[k, (Subscript[5, 1])]) -> 1, 
            (Subscript[k, 6]) -> 1, (Subscript[k, 7]) -> 0.1, (Subscript[k, (Subscript[7, 1])]) -> 1, 
            (Subscript[k, 8]) -> 0.5,  (Subscript[J, 3]) -> 0.08,  (Subscript[J, 4]) -> 0.01,  
            (Subscript[J, 5]) -> 0.05, (Subscript[J, 6]) -> 0.5, (Subscript[J, 7]) -> 0.05, 
            (Subscript[J, 8]) -> 0.05,  (Subscript[t, APC]) -> 4,  (Subscript[t, CDH1]) -> 1,  
            (Subscript[t, CDC20]) -> 7}, {(Subscript[c, "MPF"])[t], (Subscript[c, "CDC20a"])[t], 
          (Subscript[c, "APCP"])[t], (Subscript[c, "CDH1a"])[t]}, {t, 0, 100}]

Plot[
      Evaluate[{(Subscript[c, "MPF"])[t], (Subscript[c, "CDC20a"])[t], (Subscript[c, "APCP"])[t], 
            (Subscript[c, "CDH1a"])[t]} /. NDSo], {t, 0, 100}, 
      PlotStyle -> {RGBColor[1, 0, 0], RGBColor[0, 0, 1], 
          RGBColor[0, 1, 0], RGBColor[0, 1, 1]}]

param = {(Subscript[M, 0]) -> 0, (Subscript[C, 0]) -> 0, (Subscript[k, 1]) -> 0, 
        (Subscript[k, 2]) -> 0.1, (Subscript[k, (Subscript[2, 1])]) -> 1, (Subscript[k, (Subscript[2, 2])]) -> 3, 
        (Subscript[k, 3]) -> 0.8, (Subscript[k, 4]) -> 0.1, (Subscript[k, (Subscript[4, 1])]) -> 1, 
        (Subscript[k, 5]) -> 0.1, (Subscript[k, (Subscript[5, 1])]) -> 1, (Subscript[k, 6]) -> 1, 
        (Subscript[k, 7]) -> 0.1, (Subscript[k, (Subscript[7, 1])]) -> 1, (Subscript[k, 8]) -> 0.5,  
        (Subscript[J, 3]) -> 0.08,  (Subscript[J, 4]) -> 0.01,  (Subscript[J, 5]) -> 0.05, 
        (Subscript[J, 6]) -> 0.5, (Subscript[J, 7]) -> 0.05, (Subscript[J, 8]) -> 0.05,  
        (Subscript[t, APC]) -> 4,  (Subscript[t, CDH1]) -> 1,  
        (Subscript[t, CDC20]) -> 7}

kisparam = {(Subscript[M, 0]) -> 0, (Subscript[C, 0]) -> 0, (Subscript[k, 2]) -> 0.1, 
        (Subscript[k, (Subscript[2, 1])]) -> 1, (Subscript[k, (Subscript[2, 2])]) -> 3, (Subscript[k, 3]) -> 0.8, 
        (Subscript[k, 4]) -> 0.1, (Subscript[k, (Subscript[4, 1])]) -> 1, (Subscript[k, 5]) -> 0.1, 
        (Subscript[k, (Subscript[5, 1])]) -> 1, (Subscript[k, 6]) -> 1, (Subscript[k, 7]) -> 0.1, 
        (Subscript[k, (Subscript[7, 1])]) -> 1, (Subscript[k, 8]) -> 0.5,  (Subscript[J, 3]) -> 0.08,  
        (Subscript[J, 4]) -> 0.01,  (Subscript[J, 5]) -> 0.05, (Subscript[J, 6]) -> 0.5, 
        (Subscript[J, 7]) -> 0.05, (Subscript[J, 8]) -> 0.05,  (Subscript[t, APC]) -> 4,  
        (Subscript[t, CDH1]) -> 1,  (Subscript[t, CDC20]) -> 7}

R5;
rhs2 = R5[[All, 2]];
rhsval = rhs2 /. kisparam;
tr = Together /@ rhsval

Solve[
      Thread[rhsval == 0], {(Subscript[c, "MPF"])[t], (Subscript[c, "CDC20a"])[t], 
        (Subscript[c, "APCP"])[t], (Subscript[c, "CDH1a"])[t]}]

ntr = Numerator /@ tr

NSolve[
      Thread[ntr == 0] /. {(Subscript[k, 1]) -> 2}, {(Subscript[c, "MPF"])[t], 
        (Subscript[c, "CDC20a"])[t], (Subscript[c, "APCP"])[t], (Subscript[c, "CDH1a"])[t]}]

Az egyensúlyi helyzeteket nem sikerült (Subscript[k, 1]) függvényében meghatározni. Nézzük pl.a (Subscript[k, 1])=2-t!

Egyensúlyi helyzetek:

\(TableForm[
      ns = NSolve[{k\_1 - 0.1`\ c\_"MPF"[t] - 
                  1.`\ c\_"CDC20a"[t]\ c\_"MPF"[t] - 
                  3.`\ c\_"CDH1a"[t]\ c\_"MPF"[t] == 
                0, \((\(\(0.35000000000000003`\)\(\)\) + 
                      3.5`\ c\_"APCP"[t] - 
                      6.3999999999999995`\ c\_"CDC20a"[t] + 
                      6.5`\ c\_"APCP"[t]\ c\_"CDC20a"[t] + 
                      0.9`\ c\_"CDC20a"[t]\^2 - 
                      1.`\ c\_"APCP"[
                          t]\ c\_"CDC20a"[t]\^2)\)/\((\((\(\(7.05`\)\(\)\) - 1.`\ c\_"CDC20a"[
                              t])\)\ \((\(\(0.5`\)\(\)\) + 
                          c\_"CDC20a"[t])\))\) == 
                0, \((\(\(0.032`\)\(\)\) - 
                      2.8160000000000003`\ c\_"APCP"[t] + 
                      0.7000000000000001`\ c\_"APCP"[t]\^2 + 
                      0.32`\ c\_"MPF"[t] + 
                      3.92`\ c\_"APCP"[t]\ c\_"MPF"[t] - 
                      1.`\ c\_"APCP"[t]\^2\ c\_"MPF"[
                          t])\)/\((\((\(\(4.01`\)\(\)\) - 
                          1.`\ c\_"APCP"[
                              t])\)\ \((\(\(0.08`\)\(\)\) + 
                          c\_"APCP"[t])\))\) == 
                0, \((\(\(0.025`\)\(\)\) + 
                      0.37`\ c\_"CDH1a"[t] - 0.4`\ c\_"CDH1a"[t]\^2 - 
                      1.05`\ c\_"CDH1a"[t]\ c\_"MPF"[t] + 
                      c\_"CDH1a"[t]\^2\ c\_"MPF"[
                          t])\)/\((\((\(\(1.05`\)\(\)\) - 
                          1.`\ c\_"CDH1a"[
                              t])\)\ \((\(\(0.05`\)\(\)\) + 
                          c\_"CDH1a"[t])\))\) == 
                0} /. {k\_1 -> 2}, {c\_"MPF"[t], c\_"CDC20a"[t], 
            c\_"APCP"[t], c\_"CDH1a"[t]}]]\)

A sajátértékek kiszámítása:

TableForm[jok = ns[[{1, 4, 7}]]]

ClearAll[jakab];
MatrixForm[
      jakab = {jakab1 = 
            D[rhs2 /. param, (Subscript[c, "MPF"])[t]] /. x_[t] -> x // 
              Simplify, 
jakab2 = 
            D[rhs2 /. param, (Subscript[c, "CDC20a"])[t]] /. x_[t] -> x // 
              Simplify, 
jakab3 = 
            D[rhs2 /. param, (Subscript[c, "APCP"])[t]] /. x_[t] -> x // 
              Simplify, 
jakab4 = 
            D[rhs2 /. param, (Subscript[c, "CDH1a"])[t]] /. x_[t] -> x // 
              Simplify}]

jobbak = jok /. {x_[t] -> x}

MatrixForm[jak1 = jakab /. jobbak[[1]]]
Eigenvalues[jak1]

Az összes sajátérték valós része negatív, így ez az egyesnúlyi állapot stabil.

MatrixForm[jak2 = jakab /. jobbak[[2]]]
Eigenvalues[jak2]

Az összes sajátérték valós része negatív, így ez az egyensúlyi állapot is stabil.

\(\(\(
\)\(Vizsgáljuk\ a\ trajektóriákat, és\ a\ \(fázismezőt!\) \(k\_1\) 
      paraméter\ értéke\ legyen\ \(α!\)\ Ez\ lesz\ a\ bifurkációs\ paraméterünk . 
        Az\ látszik\ az\ ábrákból, 
    hogy\ α - t\ 0.5 - nek\ választva, 
    a\ megoldások\ tartanak\ egy\ stabil\ pont\ felé, 
    ez\ a\ pont\ az\ alsó\ diagramról\ leolvasható . α = 
      2 - nél\ elveszti\ a\ stabilitását, és\ egy\ ciklusba\ megy\ át, 
    majd\ α = 
      4\ paraméterértéknél\ ismét\ csomópont\ jön\ \(\(létre\)\(.\)\)\)\)\)

NDS1[α_, {x0_, y0_}] := 
      Module[{}, 
        NDSolve[Union[{R5, (Subscript[c, "MPF"])[0] == (Subscript[M, 0]), 
                (Subscript[c, "CDC20a"])[0] == 0, (Subscript[c, "APCP"])[0] == 0, 
                (Subscript[c, "CDH1a"])[0] == (Subscript[C, 0])}] /. {(Subscript[M, 0]) -> x0, 
              (Subscript[C, 0]) -> y0, (Subscript[k, 1]) -> α, (Subscript[k, 2]) -> 0.1, 
              (Subscript[k, (Subscript[2, 1])]) -> 1, (Subscript[k, (Subscript[2, 2])]) -> 3, (Subscript[k, 3]) -> 0.8, 
              (Subscript[k, 4]) -> 0.1, (Subscript[k, (Subscript[4, 1])]) -> 1, (Subscript[k, 5]) -> 0.1, 
              (Subscript[k, (Subscript[5, 1])]) -> 1, (Subscript[k, 6]) -> 1, (Subscript[k, 7]) -> 0.1, 
              (Subscript[k, (Subscript[7, 1])]) -> 1, (Subscript[k, 8]) -> 0.5,  (Subscript[J, 3]) -> 0.08,  
              (Subscript[J, 4]) -> 0.01,  (Subscript[J, 5]) -> 0.05, (Subscript[J, 6]) -> 0.5, 
              (Subscript[J, 7]) -> 0.05, (Subscript[J, 8]) -> 0.05,  (Subscript[t, APC]) -> 4,  
              (Subscript[t, CDH1]) -> 1,  (Subscript[t, CDC20]) -> 7}, {(Subscript[c, "MPF"])[t], 
            (Subscript[c, "CDC20a"])[t], (Subscript[c, "APCP"])[t], (Subscript[c, "CDH1a"])[t]}, {t, 0, 
            200}]]

trajectory[α_, {x0_, y0_}] := 
      Module[{}, {(Subscript[c, "MPF"])[t], 
            (Subscript[c, "CDH1a"])[
              t]} /. NDS1[α, {x0, 
                y0}][[1]]]

PhaseAndField[α_, 
        ic : {{_, _} .. }, {MPFmin_, MPFmax_}, {CDH1amin_, 
          CDH1amax_}, {tmin_, tmax_}, opts___Rule]:=With[{optParametricPlot  =  
              FilterOptions[ParametricPlot, opts], 
            optPlotField = FilterOptions[PlotVectorField, opts], 
            optND  = FilterOptions[NDSolve, opts]}, 
DisplayTogether[
PlotVectorField[ {α-(0.1 + 3 + 3 (Subscript[c, "CDH1a"])[t])(Subscript[c, "MPF"])[t],-( 
                              0.1 + 1 (Subscript[c, "MPF"])[t]) (Subscript[c, "CDH1a"])[
                              t]/(0.05 +  (Subscript[c, "CDH1a"])[t]) + 
                        0.5 ( 
                              1 -  (Subscript[c, "CDH1a"])[t])/(0.05 + 1 -  
                                (Subscript[c, "CDH1a"])[t])},{(Subscript[c, "MPF"])[t], MPFmin, MPFmax},{(Subscript[c, "CDH1a"])[t], CDH1amin, CDH1amax},PlotPoints -> Automatic,optPlotField],
ParametricPlot[
                Evaluate[
                  Table[trajectory[α, 
                      ic[[i]]], {i, 
                      Length[ic]}]], {t, tmin, tmax}, optParametricPlot],Frame -> True,FrameLabel -> {"MPF", "CDH1a"},PlotLabel -> "α=" <> 
                  ToString[α]]]

PhaseAndField[2, 
        Table[{x0, 1}, {x0, 0.1, 2, 0.1}], {0, 1.5}, {0, 1}, {20, 
          80}
];

Plot[Evaluate[trajectory[2, {0.5, 1}]], {t, 0, 50}, 
        PlotRange -> {Automatic, {0, 3}}, 
        PlotStyle -> {{}, Dashing[{0.02, 0.02}]}, Frame -> True, 
        FrameLabel -> {"t", "concentration"}, 
        PlotLegend -> {"MPF", "CDH1a"}, 
        LegendPosition -> {-0.55, 0.35}, 
        LegendSize -> {0.7, 0.2}, 
        LegendShadow -> None];

PhaseAndField[4, 
        Table[{x0, 0.5}, {x0, 0.1, 2, 0.1}], {0.52, 0.7}, {0.07, 0.15}, {10, 
          100}
];

Plot[Evaluate[trajectory[4, {0.5, 0.5}]], {t, 0, 100}, 
        PlotRange -> {Automatic, {0, 3}}, 
        PlotStyle -> {{}, Dashing[{0.02, 0.02}]}, Frame -> True, 
        FrameLabel -> {"t", "concentration"}, 
        PlotLegend -> {"MPF", "CDH1a"}, 
        LegendPosition -> {-0.55, 0.35}, 
        LegendSize -> {0.7, 0.2}, 
        LegendShadow -> None];

PhaseAndField[0.5, 
        Table[{x0, 1}, {x0, 0.1, 2, 0.1}], {0, 0.4}, {0.4, 1}, {2, 
          200}
];

Plot[Evaluate[trajectory[0.5, {0.5, 1}]], {t, 0, 100}, 
        PlotRange -> {Automatic, {0, 3}}, 
        PlotStyle -> {{}, Dashing[{0.02, 0.02}]}, Frame -> True, 
        FrameLabel -> {"t", "concentration"}, 
        PlotLegend -> {"MPF", "CDH1a"}, 
        LegendPosition -> {-0.55, 0.35}, 
        LegendSize -> {0.7, 0.2}, 
        LegendShadow -> None];