# 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;
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*20;
Ksvdweak = Ksvd*20;
Stable2 = 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 /. 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,
];
If[
Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0.0001 &&
Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0.0001,
];
]
]
]
(*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: ",
];
If[ (Stable2 == 1),
(*Print["MONOSTABLE ",  \[Alpha],"   ",\[Beta],",  Saddle: ",
];
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 = 3;
n2 = 2;
(*Ksrd=;*)
Knsrd = 10000;
Stable1 = 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*20;
Ksvdweak = Ksvd*20;
Stable2 = 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 /. 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,
];
If[
Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0.0001 &&
Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0.0001,
];
]
]
]
(*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: ",
];
If[ (Stable2 == 1),
(*Print["MONOSTABLE ",  \[Alpha],"   ",\[Beta],",  Saddle: ",
];
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 = 1;
n2 = 4;
(*Ksrd=;*)
Knsrd = 10000;
Stable1 = 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*20;
Ksvdweak = Ksvd*20;
Stable2 = 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 /. 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,
];
If[
Re[Eigenvalues[Jac[xtemp, ytemp]]1] > 0.0001 &&
Re[Eigenvalues[Jac[xtemp, ytemp]]2] < 0.0001,
];
]
]
]
(*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: ",
];
If[ (Stable2 == 1),
(*Print["MONOSTABLE ",  \[Alpha],"   ",\[Beta],",  Saddle: ",
];
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 Unbalanced Binding Strengths",
Directive[Large, Bold]], LabelStyle -> Directive[Medium, Bold]]
```