Team:Duke/Modeling/Codes/Kinetic3

{{Template:Team:Duke}}

<div class="right_box">
= '''Mathematical Modeling of Bistable Toggle Switch''' =
 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;
 n1 = 5;
 n2 = 0;
 (*Ksrd=;*)
 Knsrd = 10000;
 Stable1 = 0;
 Saddle1 = 0;
 Unstable1 = 0;
 matUp = {};
 matDown = {};
 RangeUp = {};
 RangeDown = {};
 For[Ksud = 0.1*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
 Print[Ksud]
  For[Ksvd = 0.1*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
   Ksudweak = Ksud*10;
   Ksvdweak = Ksvd*10;
   Stable2 = 0;
   Saddle2 = 0;
   Unstable2 = 0;
   delEud = kb*T*Log[10, Ksud/Knsrd];
   delEvd = kb*T*Log[10, Ksvd/Knsrd];
   delEudweak = kb*T*Log[10, Ksudweak/Knsrd];
   delEvdweak = kb*T*Log[10, Ksvdweak/Knsrd];
   xdot[x_, y_] := 
    K1/(1 + \[Alpha] ((1 + (y/Nns)*E^(-delEvd/(kb*T)))^
           n1)*((1 + (y/Nns)*E^(-delEvdweak/(kb*T)))^n2)) - x; 
   ydot[x_, y_] := 
    K2/(1 + \[Alpha] ((1 + (x/Nns)*E^(-delEud/(kb*T)))^
           n1)*((1 + (x/Nns)*E^(-delEudweak/(kb*T)))^n2)) - 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 /. fps[[i]];
     ytemp = y /. fps[[i]];
     (*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 = 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;
 n1 = 5;
 n2 = 0;
 (*Ksrd=;*)
 Knsrd = 10000;
 Stable1 = 0;
 Saddle1 = 0;
 Unstable1 = 0;
 matUp1 = {};
 matDown1 = {};
 RangeUp1 = {};
 RangeDown1 = {};
 For[Ksud = 0.1*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
  Print[Ksud]
  For[Ksvd = 0.1*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
   Ksudweak = Ksud*10;
   Ksvdweak = Ksvd*10;
   Stable2 = 0;
   Saddle2 = 0;
   Unstable2 = 0;
   delEud = kb*T*Log[10, Ksud/Knsrd];
   delEvd = kb*T*Log[10, Ksvd/Knsrd];
   delEudweak = kb*T*Log[10, Ksudweak/Knsrd];
   delEvdweak = kb*T*Log[10, Ksvdweak/Knsrd];
   xdot[x_, y_] := 
    K1/(1 + \[Alpha] ((1 + (y/Nns)*E^(-delEvd/(kb*T)))^
           n1)*((1 + (y/Nns)*E^(-delEvdweak/(kb*T)))^n2)) - x + 300; 
   ydot[x_, y_] := 
    K2/(1 + \[Alpha] ((1 + (x/Nns)*E^(-delEud/(kb*T)))^
           n1)*((1 + (x/Nns)*E^(-delEudweak/(kb*T)))^n2)) - y + 300; 
   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 /. fps[[i]];
     ytemp = y /. fps[[i]];
     (*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]}];*)
     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 = 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;
 n1 = 5;
 n2 = 0;
 (*Ksrd=;*)
 Knsrd = 10000;
 Stable1 = 0;
 Saddle1 = 0;
 Unstable1 = 0;
 matUp2 = {};
 matDown2 = {};
 RangeUp2 = {};
 RangeDown2 = {};
 For[Ksud = 0.1*10^-9, Ksud <= 7000*10^-9, Ksud = Ksud*1.1,
 Print[Ksud]
  For[Ksvd = 0.1*10^-9, Ksvd <= 7000*10^-9, Ksvd = Ksvd*1.1,
   Ksudweak = Ksud*10;
   Ksvdweak = Ksvd*10;
   Stable2 = 0;
   Saddle2 = 0;
   Unstable2 = 0;
   delEud = kb*T*Log[10, Ksud/Knsrd];
   delEvd = kb*T*Log[10, Ksvd/Knsrd];
   delEudweak = kb*T*Log[10, Ksudweak/Knsrd];
   delEvdweak = kb*T*Log[10, Ksvdweak/Knsrd];
   xdot[x_, y_] := 
    K1/(1 + \[Alpha] ((1 + (y/Nns)*E^(-delEvd/(kb*T)))^
           n1)*((1 + (y/Nns)*E^(-delEvdweak/(kb*T)))^n2)) - x + 600; 
   ydot[x_, y_] := 
    K2/(1 + \[Alpha] ((1 + (x/Nns)*E^(-delEud/(kb*T)))^
           n1)*((1 + (x/Nns)*E^(-delEudweak/(kb*T)))^n2)) - y + 600; 
   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 /. fps[[i]];
     ytemp = y /. fps[[i]];
     (*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]}];*)
     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;
   ]
  ]
 
 RangeUp;
 RangeDown;
 ListPlot[{RangeUp, RangeDown}];
 matUp;
 matDown;
 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 Leaky Promoters", 
   Directive[Large, Bold]], LabelStyle -> Directive[Medium, Bold]]

</div>