Team:Duke/Modeling/Codes/Kinetic1

From 2013.igem.org

(Difference between revisions)
(→Mathematical Modeling of Bistable Toggle Switch)
(→Mathematical Modeling of Bistable Toggle Switch)
Line 2: Line 2:
<div class="right_box">
<div class="right_box">
-
= '''Mathematical Modeling of Bistable Toggle Switch''' =
+
  Clear["Global`*"]
-
 
+
  k1 = 100;
-
[[[
+
  d1 = 0.007;
-
88
+
  k2 = 100;
-
]]]
+
  d2 = 0.007;
-
Clear["Global`*"]
+
  K1 = k1/d1;
-
k1 = 100;
+
  K2 = k2/d2;
-
d1 = 0.007;
+
  Nns = 5*10^6;
-
k2 = 100;
+
  delEpd = -1.0421*10^-19;
-
d2 = 0.007;
+
  kb = 1.3806*10^-23;
-
K1 = k1/d1;
+
  T = 298;
-
K2 = k2/d2;
+
  P = 3000;
-
Nns = 5*10^6;
+
  \[Alpha] = Nns*E^(delEpd/(kb*T))/P;
-
delEpd = -1.0421*10^-19;
+
  n = 5;
-
kb = 1.3806*10^-23;
+
  (*Ksrd=;*)
-
T = 298;
+
  Knsrd = 10000;
-
P = 3000;
+
  Stable1 = 0;
-
\[Alpha] = Nns*E^(delEpd/(kb*T))/P;
+
  Saddle1 = 0;
-
n = 5;
+
  Unstable1 = 0;
-
(*Ksrd=;*)
+
  matUp = {};
-
Knsrd = 10000;
+
  matDown = {};
-
Stable1 = 0;
+
  RangeUp = {};
-
Saddle1 = 0;
+
  RangeDown = {};
-
Unstable1 = 0;
+
  For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
-
matUp = {};
+
  Print[Ksud]
-
matDown = {};
+
    For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
-
RangeUp = {};
+
    Stable2 = 0;
-
RangeDown = {};
+
    Saddle2 = 0;
-
For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
+
    Unstable2 = 0;
-
Print[Ksud]
+
    delEud = kb*T*Log[10, Ksud/Knsrd];
-
  For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
+
    delEvd = kb*T*Log[10, Ksvd/Knsrd];
-
  Stable2 = 0;
+
    xdot[x_, y_] :=  
-
  Saddle2 = 0;
+
      K1/(1 + \[Alpha] (1 + (y/Nns)*E^(-delEvd/(kb*T)))^n) - x;  
-
  Unstable2 = 0;
+
    ydot[x_, y_] :=  
-
  delEud = kb*T*Log[10, Ksud/Knsrd];
+
      K2/(1 + \[Alpha] (1 + (x/Nns)*E^(-delEud/(kb*T)))^n) - y;  
-
  delEvd = kb*T*Log[10, Ksvd/Knsrd];
+
    Jac[x_,  
-
  xdot[x_, y_] :=  
+
      y_] := {{Derivative[1, 0][xdot][x, y],  
-
    K1/(1 + \[Alpha] (1 + (y/Nns)*E^(-delEvd/(kb*T)))^n) - x;  
+
        Derivative[0, 1][xdot][x, y]}, {Derivative[1, 0][ydot][x, y],  
-
  ydot[x_, y_] :=  
+
        Derivative[0, 1][ydot][x, y]}};
-
    K2/(1 + \[Alpha] (1 + (x/Nns)*E^(-delEud/(kb*T)))^n) - y;  
+
    fps = NSolve[xdot[x, y] == 0 && ydot[x, y] == 0, {x, y}, Reals];
-
  Jac[x_,  
+
    (*Print[N[fps,3]]*)
-
    y_] := {{Derivative[1, 0][xdot][x, y],  
+
    (*Print[Length[fps]-1]*)
-
      Derivative[0, 1][xdot][x, y]}, {Derivative[1, 0][ydot][x, y],  
+
    For[i = 1, i < Length[fps] + 1, i++,
-
      Derivative[0, 1][ydot][x, y]}};
+
      xtemp = x /. fps[[i]];
-
  fps = NSolve[xdot[x, y] == 0 && ydot[x, y] == 0, {x, y}, Reals];
+
      ytemp = y /. fps[[i]];
-
  (*Print[N[fps,3]]*)
+
      (*Print[N[xtemp,3]]*)
-
  (*Print[Length[fps]-1]*)
+
      (*Print[N[ytemp,3]]*)
-
  For[i = 1, i < Length[fps] + 1, i++,
+
      If[Im[xtemp] == 0 && Im[ytemp] == 0,
-
    xtemp = x /. fps[[i]];
+
        If[xtemp > 0 && ytemp > 0,
-
    ytemp = y /. fps[[i]];
+
        (*Print["              ",i-1,"  R1=",xtemp,"  R2=",ytemp,  
-
    (*Print[N[xtemp,3]]*)
+
        "  Ksud=",Ksud,"  Ksvd=", Ksvd]*)
-
    (*Print[N[ytemp,3]]*)
+
        If[
-
    If[Im[xtemp] == 0 && Im[ytemp] == 0,
+
          Re[Eigenvalues[Jac[xtemp, ytemp]][[1]]] > 0 &&  
-
      If[xtemp > 0 && ytemp > 0,
+
          Re[Eigenvalues[Jac[xtemp, ytemp]][[2]]] > 0,
-
      (*Print["              ",i-1,"  R1=",xtemp,"  R2=",ytemp,  
+
          Unstable2 += 1;
-
      "  Ksud=",Ksud,"  Ksvd=", Ksvd]*)
+
          ];
-
      If[
+
        If[
-
        Re[Eigenvalues[Jac[xtemp, ytemp]][[1]]] > 0 &&  
+
          Re[Eigenvalues[Jac[xtemp, ytemp]][[1]]] < 0 &&  
-
        Re[Eigenvalues[Jac[xtemp, ytemp]][[2]]] > 0,
+
          Re[Eigenvalues[Jac[xtemp, ytemp]][[2]]] < 0,
-
        Unstable2 += 1;
+
          Stable2 += 1;
-
        ];
+
          ];
-
      If[
+
        If[
-
        Re[Eigenvalues[Jac[xtemp, ytemp]][[1]]] < 0 &&  
+
          Re[Eigenvalues[Jac[xtemp, ytemp]][[1]]] < 0.0001 &&  
-
        Re[Eigenvalues[Jac[xtemp, ytemp]][[2]]] < 0,
+
            Re[Eigenvalues[Jac[xtemp, ytemp]][[2]]] > 0.0001,
-
        Stable2 += 1;
+
          Saddle2 += 1;
-
        ];
+
          ];
-
      If[
+
          If[
-
        Re[Eigenvalues[Jac[xtemp, ytemp]][[1]]] < 0.0001 &&  
+
-
        Re[Eigenvalues[Jac[xtemp, ytemp]][[2]]] > 0.0001,
+
-
        Saddle2 += 1;
+
-
        ];
+
-
      If[
+
         Re[Eigenvalues[Jac[xtemp, ytemp]][[1]]] > 0.0001 &&  
         Re[Eigenvalues[Jac[xtemp, ytemp]][[1]]] > 0.0001 &&  
         Re[Eigenvalues[Jac[xtemp, ytemp]][[2]]] < 0.0001,
         Re[Eigenvalues[Jac[xtemp, ytemp]][[2]]] < 0.0001,
Line 110: Line 105:
-
k1 = 100;
+
  k1 = 100;
-
d1 = 0.007;
+
  d1 = 0.007;
-
k2 = 50;
+
  k2 = 50;
-
d2 = 0.007;
+
  d2 = 0.007;
-
K1 = k1/d1;
+
  K1 = k1/d1;
-
K2 = k2/d2;
+
    K2 = k2/d2;
-
Nns = 5*10^6;
+
    Nns = 5*10^6;
-
delEpd = -1.0421*10^-19;
+
  delEpd = -1.0421*10^-19;
-
kb = 1.3806*10^-23;
+
  kb = 1.3806*10^-23;
-
T = 298;
+
  T = 298;
-
P = 3000;
+
  P = 3000;
-
\[Alpha] = Nns*E^(delEpd/(kb*T))/P;
+
  \[Alpha] = Nns*E^(delEpd/(kb*T))/P;
-
n = 5;
+
  n = 5;
-
(*Ksrd=;*)
+
  (*Ksrd=;*)
-
Knsrd = 10000;
+
  Knsrd = 10000;
-
Stable1 = 0;
+
  Stable1 = 0;
-
Saddle1 = 0;
+
  Saddle1 = 0;
-
Unstable1 = 0;
+
  Unstable1 = 0;
-
matUp1 = {};
+
  matUp1 = {};
-
matDown1 = {};
+
  matDown1 = {};
-
RangeUp1 = {};
+
  RangeUp1 = {};
-
RangeDown1 = {};
+
  RangeDown1 = {};
-
For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
+
  For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
  Print[Ksud]
  Print[Ksud]
   For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
   For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
Line 215: Line 210:
-
k1 = 100;
+
  k1 = 100;
-
d1 = 0.007;
+
  d1 = 0.007;
-
k2 = 10;
+
  k2 = 10;
-
d2 = 0.007;
+
  d2 = 0.007;
-
K1 = k1/d1;
+
  K1 = k1/d1;
-
K2 = k2/d2;
+
  K2 = k2/d2;
-
Nns = 5*10^6;
+
  Nns = 5*10^6;
-
delEpd = -1.0421*10^-19;
+
  delEpd = -1.0421*10^-19;
-
kb = 1.3806*10^-23;
+
  kb = 1.3806*10^-23;
-
T = 298;
+
  T = 298;
-
P = 3000;
+
  P = 3000;
-
\[Alpha] = Nns*E^(delEpd/(kb*T))/P;
+
  \[Alpha] = Nns*E^(delEpd/(kb*T))/P;
-
n = 5;
+
  n = 5;
-
(*Ksrd=;*)
+
  (*Ksrd=;*)
-
Knsrd = 10000;
+
  Knsrd = 10000;
-
Stable1 = 0;
+
  Stable1 = 0;
-
Saddle1 = 0;
+
  Saddle1 = 0;
-
Unstable1 = 0;
+
  Unstable1 = 0;
-
matUp2 = {};
+
  matUp2 = {};
-
matDown2 = {};
+
  matDown2 = {};
-
RangeUp2 = {};
+
  RangeUp2 = {};
-
RangeDown2 = {};
+
  RangeDown2 = {};
-
For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
+
  For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
  Print[Ksud]
  Print[Ksud]
   For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
   For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
Line 335: Line 330:
   Style["Bistable Regions for Unbalanced Promoter Strengths",  
   Style["Bistable Regions for Unbalanced Promoter Strengths",  
   Directive[Large, Bold]], LabelStyle -> Directive[Medium, Bold]]
   Directive[Large, Bold]], LabelStyle -> Directive[Medium, Bold]]
-
 
-
 
-
 
-
 
-
 
-
</div>
 

Revision as of 09:27, 22 September 2013

  • 1 1
  • 2 2
  • 3 3
  • 4 4
  • 5 5
  • 6 6
  • 7 7
  • 8 8
  • 9 9
  • 10 10 


 Clear["Global`*"]
 k1 = 100;
 d1 = 0.007;
 k2 = 100;
 d2 = 0.007;
 K1 = k1/d1;
 K2 = k2/d2;
 Nns = 5*10^6;
 delEpd = -1.0421*10^-19;
 kb = 1.3806*10^-23;
 T = 298;
 P = 3000;
 \[Alpha] = Nns*E^(delEpd/(kb*T))/P;
 n = 5;
 (*Ksrd=;*)
 Knsrd = 10000;
 Stable1 = 0;
 Saddle1 = 0;
 Unstable1 = 0;
 matUp = {};
 matDown = {};
 RangeUp = {};
 RangeDown = {};
 For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
  Print[Ksud]
   For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
    Stable2 = 0;
    Saddle2 = 0;
    Unstable2 = 0;
    delEud = kb*T*Log[10, Ksud/Knsrd];
    delEvd = kb*T*Log[10, Ksvd/Knsrd];
    xdot[x_, y_] := 
     K1/(1 + \[Alpha] (1 + (y/Nns)*E^(-delEvd/(kb*T)))^n) - x; 
    ydot[x_, y_] := 
     K2/(1 + \[Alpha] (1 + (x/Nns)*E^(-delEud/(kb*T)))^n) - y; 
    Jac[x_, 
      y_] := {{Derivative[1, 0][xdot][x, y], 
       Derivative[0, 1][xdot][x, y]}, {Derivative[1, 0][ydot][x, y], 
       Derivative[0, 1][ydot][x, y]}};
    fps = NSolve[xdot[x, y] == 0 && ydot[x, y] == 0, {x, y}, Reals];
    (*Print[N[fps,3]]*)
    (*Print[Length[fps]-1]*)
    For[i = 1, i < Length[fps] + 1, i++,
      xtemp = x /. fpsi;
      ytemp = y /. fpsi;
      (*Print[N[xtemp,3]]*)
      (*Print[N[ytemp,3]]*)
      If[Im[xtemp] == 0 && Im[ytemp] == 0,
       If[xtemp > 0 && ytemp > 0,
        (*Print["              ",i-1,"  R1=",xtemp,"  R2=",ytemp, 
        "  Ksud=",Ksud,"   Ksvd=", Ksvd]*)
        If[
         Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0 && 
          Re[Eigenvalues[Jac[xtemp, ytemp]]2] > 0,
         Unstable2 += 1;
         ];
        If[
         Re[Eigenvalues[Jac[xtemp, ytemp]]1] < 0 && 
          Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0,
         Stable2 += 1;
         ];
        If[
         Re[Eigenvalues[Jac[xtemp, ytemp]]1] < 0.0001 && 
            Re[Eigenvalues[Jac[xtemp, ytemp]]2] > 0.0001,
         Saddle2 += 1;
         ];
          If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0.0001 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0.0001,
       Saddle2 += 1;
       ];
      ]
     ]
    ]
   (*Print[Stable2, "  ", Stable1]*)
   If[ (Stable2 == 2 && Stable1 == 1),
    Print["Switch  (UP)"];
    (*Print[{\[Alpha],\[Beta]}];*)
    matUp = Insert[matUp, {Log[10, Ksud], Log[10, Ksvd]}, 1];
    RangeUp = Insert[RangeUp, {Ksud, Ksvd}, 1];
    ];
  If[ ((Stable2 == 1 && Stable1 == 2)),
   Print["Switch (DOWN)"];
   (*Print[{\[Alpha],\[Beta]}];*)
   matDown = Insert[matDown, {Log[10, Ksud], Log[10, Ksvd]}, 1];
   RangeDown = Insert[RangeDown, {Ksud, Ksvd}, 1];
   ];
  If[(Stable2 == 2),
   (*Print["BISTABLE ", \[Alpha], "   ",\[Beta], ",  Saddle: ",
   Saddle2 ];*)
   ];
  If[ (Stable2 == 1),
   (*Print["MONOSTABLE ",  \[Alpha],"   ",\[Beta],",  Saddle: ",
   Saddle2 ];*)
   ];
  Saddle1 = Saddle2;
  Unstable1 = Unstable2;
  Stable1 = Stable2;
  ]
]



 k1 = 100;
 d1 = 0.007;
 k2 = 50;
 d2 = 0.007;
 K1 = k1/d1;
   K2 = k2/d2;
   Nns = 5*10^6;
 delEpd = -1.0421*10^-19;
 kb = 1.3806*10^-23;
 T = 298;
 P = 3000;
 \[Alpha] = Nns*E^(delEpd/(kb*T))/P;
 n = 5;
 (*Ksrd=;*)
 Knsrd = 10000;
 Stable1 = 0;
 Saddle1 = 0;
 Unstable1 = 0;
 matUp1 = {};
 matDown1 = {};
 RangeUp1 = {};
 RangeDown1 = {};
 For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
Print[Ksud]
 For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
  Stable2 = 0;
  Saddle2 = 0;
  Unstable2 = 0;
  delEud = kb*T*Log[10, Ksud/Knsrd];
  delEvd = kb*T*Log[10, Ksvd/Knsrd];
  xdot[x_, y_] := 
   K1/(1 + \[Alpha] (1 + (y/Nns)*E^(-delEvd/(kb*T)))^n) - x; 
  ydot[x_, y_] := 
   K2/(1 + \[Alpha] (1 + (x/Nns)*E^(-delEud/(kb*T)))^n) - y; 
  Jac[x_, 
    y_] := {{Derivative[1, 0][xdot][x, y], 
     Derivative[0, 1][xdot][x, y]}, {Derivative[1, 0][ydot][x, y], 
     Derivative[0, 1][ydot][x, y]}};
  fps = NSolve[xdot[x, y] == 0 && ydot[x, y] == 0, {x, y}, Reals];
  (*Print[N[fps,3]]*)
  (*Print[Length[fps]-1]*)
  For[i = 1, i < Length[fps] + 1, i++,
    xtemp = x /. fpsi;
    ytemp = y /. fpsi;
    (*Print[N[xtemp,3]]*)
    (*Print[N[ytemp,3]]*)
    If[Im[xtemp] == 0 && Im[ytemp] == 0,
     If[xtemp > 0 && ytemp > 0,
      (*Print["              ",i-1,"  R1=",xtemp,"  R2=",ytemp, 
      "  Ksud=",Ksud,"   Ksvd=", Ksvd]*)
      If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] > 0,
       Unstable2 += 1;
       ];
      If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] < 0 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0,
       Stable2 += 1;
       ];
      If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] < 0.0001 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] > 0.0001,
       Saddle2 += 1;
       ];
      If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0.0001 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0.0001,
       Saddle2 += 1;
       (*saddlePoint2=Insert[saddlePoint2,{\[Alpha],\[Beta],Re[
       Eigenvalues[Jac[xtemp,ytemp]]2]},1];*)
       ];
      ]
     ]
    ]
   (*Print[Stable2, "  ", Stable1]*)
   If[ (Stable2 == 2 && Stable1 == 1),
    Print["Switch  (UP)"];
    (*Print[{\[Alpha],\[Beta]}];*)
    matUp1 = Insert[matUp1, {Log[10, Ksud], Log[10, Ksvd]}, 1];
    RangeUp1 = Insert[RangeUp1, {Ksud, Ksvd}, 1];
    ];
  If[ ((Stable2 == 1 && Stable1 == 2)),
   Print["Switch (DOWN)"];
   (*Print[{\[Alpha],\[Beta]}];*)
   matDown1 = Insert[matDown1, {Log[10, Ksud], Log[10, Ksvd]}, 1];
   RangeDown1 = Insert[RangeDown1, {Ksud, Ksvd}, 1];
   ];
  If[(Stable2 == 2),
   (*Print["BISTABLE ", \[Alpha], "   ",\[Beta], ",  Saddle: ",
   Saddle2 ];*)
   ];
  If[ (Stable2 == 1),
   (*Print["MONOSTABLE ",  \[Alpha],"   ",\[Beta],",  Saddle: ",
   Saddle2 ];*)
   ];
  Saddle1 = Saddle2;
  Unstable1 = Unstable2;
  Stable1 = Stable2;
  ]
]




 k1 = 100;
 d1 = 0.007;
 k2 = 10;
 d2 = 0.007;
 K1 = k1/d1;
 K2 = k2/d2;
 Nns = 5*10^6;
 delEpd = -1.0421*10^-19;
 kb = 1.3806*10^-23;
 T = 298;
 P = 3000;
 \[Alpha] = Nns*E^(delEpd/(kb*T))/P;
 n = 5;
 (*Ksrd=;*)
 Knsrd = 10000;
 Stable1 = 0;
 Saddle1 = 0;
 Unstable1 = 0;
 matUp2 = {};
 matDown2 = {};
 RangeUp2 = {};
 RangeDown2 = {};
 For[Ksud = 0.01*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
Print[Ksud]
 For[Ksvd = 0.01*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
  Stable2 = 0;
  Saddle2 = 0;
  Unstable2 = 0;
  delEud = kb*T*Log[10, Ksud/Knsrd];
  delEvd = kb*T*Log[10, Ksvd/Knsrd];
  xdot[x_, y_] := 
   K1/(1 + \[Alpha] (1 + (y/Nns)*E^(-delEvd/(kb*T)))^n) - x; 
  ydot[x_, y_] := 
   K2/(1 + \[Alpha] (1 + (x/Nns)*E^(-delEud/(kb*T)))^n) - y; 
  Jac[x_, 
    y_] := {{Derivative[1, 0][xdot][x, y], 
     Derivative[0, 1][xdot][x, y]}, {Derivative[1, 0][ydot][x, y], 
     Derivative[0, 1][ydot][x, y]}};
  fps = NSolve[xdot[x, y] == 0 && ydot[x, y] == 0, {x, y}, Reals];
  (*Print[N[fps,3]]*)
  (*Print[Length[fps]-1]*)
  For[i = 1, i < Length[fps] + 1, i++,
    xtemp = x /. fpsi;
    ytemp = y /. fpsi;
    (*Print[N[xtemp,3]]*)
    (*Print[N[ytemp,3]]*)
    If[Im[xtemp] == 0 && Im[ytemp] == 0,
     If[xtemp > 0 && ytemp > 0,
      (*Print["              ",i-1,"  R1=",xtemp,"  R2=",ytemp, 
      "  Ksud=",Ksud,"   Ksvd=", Ksvd]*)
      If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] > 0,
       Unstable2 += 1;
       ];
      If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] < 0 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0,
       Stable2 += 1;
       ];
      If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] < 0.0001 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] > 0.0001,
       Saddle2 += 1;
       ];
      If[
       Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0.0001 && 
        Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0.0001,
       Saddle2 += 1;
       (*saddlePoint2=Insert[saddlePoint2,{\[Alpha],\[Beta],Re[
       Eigenvalues[Jac[xtemp,ytemp]]2]},1];*)
       ];
      ]
     ]
    ]
   (*Print[Stable2, "  ", Stable1]*)
   If[ (Stable2 == 2 && Stable1 == 1),
    Print["Switch  (UP)"];
    (*Print[{\[Alpha],\[Beta]}];*)
    matUp2 = Insert[matUp2, {Log[10, Ksud], Log[10, Ksvd]}, 1];
    RangeUp2 = Insert[RangeUp2, {Ksud, Ksvd}, 1];
    ];
  If[ ((Stable2 == 1 && Stable1 == 2)),
   Print["Switch (DOWN)"];
   (*Print[{\[Alpha],\[Beta]}];*)
   matDown2 = Insert[matDown2, {Log[10, Ksud], Log[10, Ksvd]}, 1];
   RangeDown2 = Insert[RangeDown2, {Ksud, Ksvd}, 1];
   ];
  If[(Stable2 == 2),
   (*Print["BISTABLE ", \[Alpha], "   ",\[Beta], ",  Saddle: ",
   Saddle2 ];*)
   ];
  If[ (Stable2 == 1),
   (*Print["MONOSTABLE ",  \[Alpha],"   ",\[Beta],",  Saddle: ",
   Saddle2 ];*)
   ];
  Saddle1 = Saddle2;
  Unstable1 = Unstable2;
  Stable1 = Stable2;
  ]
]



RangeUp1; RangeDown1; ListPlot[{RangeUp1, RangeDown1}]; matUp1; matDown1; ListLinePlot[{matUp, matDown, matUp1, matDown1, matUp2, matDown2}, 

PlotStyle -> {Blue, Blue, Red, Red, Green, Green}, 
PlotRange -> {{-10, -5}, {-10, -5}}, InterpolationOrder -> 1, 
Mesh -> Full, MeshStyle -> Directive[Gray], 
AxesLabel -> {"Log(\!\(\*SubscriptBox[\(K\), \(d, Repressor\\\ \

1\)]\))", "Log(\!\(\*SubscriptBox[\(K\), \(d, Repressor\\\ 2\)]\))"}, 

PlotLabel -> 
 Style["Bistable Regions for Unbalanced Promoter Strengths",
Directive[Large, Bold]], LabelStyle -> Directive[Medium, Bold]]
Retrieved from "http://2013.igem.org/Team:Duke/Modeling/Codes/Kinetic1"