Butterfly wing patterns is one of most typical biological pattern formation. Explaining the development of these patterns of gene G switching from inactive to active status by a biochemical signal substance S can refer to the paper of Lewis et al. (1977). Gene may normally be inactive but can be switched on to produce a pigment or other gene product when the concentration of S exceeds a certain threshold. Let \(g(t)\) denote the concentration of the gene product, and assume that the concentration \(s_{0}\) of S is fixed. The model is \[\frac{dg}{dt} = k_{1}s_{0} -k_{2}g + \frac{k_{3}g^{2}}{k^{2}_{4} +g^{2}}\], which can be transformed to \[\frac{dx}{d\tau} = s-rx + \frac{x^{2}}{1+x^{2}}\], where \[\begin{align} x &=\frac{g}{k_{4}} \\ s &=\frac{k_{1}s_{0}}{k_{3}} \\ r &=\frac{k_{2}k_{4}}{k_{3}} \\ \tau &=\frac{k_{3}}{k_{4}} t \end{align}\]
There are three real fixed points: one is at zero, the other two real fixed points exit when \(1-4r^{2} \geq 0\) or \(r\leq r_{c}=\frac{1}{2}\) for all positive constants of k.
\[\begin{equation} x = \left\{ \begin{array}{cc} 0 :& (un)stable/semistable \\ \\ \frac{1+ \sqrt{1-4r^{2}}}{2r} :& unstable \\ \frac{1- \sqrt{1-4r^{2}}}{2r} :& stable \end{array} \right. \end{equation}\]
Interpretation: if increasing s from 0, and then returning back s = 0, what happen to the gene concentration or x :
- for \(r> 1/2\), the semistable fixed point \(x=0\) increases very slowly with the increment of s; and decreases very slowly with the decrement of s when back to s= 0. Therefore, when s increases from 0 and then decreases back to s= 0, the gene concentration still remains zero and inactivated or still at \(x=0\).
- for \(0\leq r< 1/2\), the stable fixed point of \(x=\frac{1- \sqrt{1-4r^{2}}}{2r}\) shifts to left(decreasing) with increment of s ; then shifts back to right(increasing) with decrement of s. Therefore, when s increases from 0 and then decreases back to s= 0, the gene concentration start decreasing from s = 0 and then increasing when back to s=0.
- for \(r = 1/2\) the system will jump between the stable fixed point of \(x=0\) and the semi-stable fixed point at \(x=\frac{1}{2r} = 1\) and back to stable fixed point of \(x=0\) finally. The change rate of jump process depends on the change rate of s. Therefore, when s increases from 0 and then decreases back to s= 0, the gene concentration start to jump from 0(inactive) to 1(active) and back to 0(inactive). The switching rate of gene concentration depends on the change rate of s.
geneswitching = function(s, r, x) {
den = (1 + x * x)
norm_r = x * x
dxdtau = s - (r * x) + (norm_r/den)
return(dxdtau)
}library(latex2exp)
x = -250:250/100
s = 0/40
r = 0
dxdtau = geneswitching(s = s, r = r, x = x)
plot(x, dxdtau, type = "l", col = "blue", lwd = 0.3, ylim = c(-0.1, 0.1), xlim = c(-1,
2), xlab = "x or gene concentration")
####
s = 1/40
r = 0
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "green", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
grid()
##
s = 2/40
r = 0
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "red", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
title(TeX("$\\frac{dx}{d\\tau} = s-rx + \\frac{x^{2}}{1+x^{2}}$"))
legend("bottomleft", border = "black", legend = c("s=0/40, r = 0", "s=1/40, r = 0",
"s=2/40, r = 0"), col = c("blue", "green", "red"), lty = 1:2, cex = 0.8)
grid()library(latex2exp)
x = -250:250/100
s = 0/40
r = 0.25
dxdtau = geneswitching(s = s, r = r, x = x)
plot(x, dxdtau, type = "l", col = "blue", lwd = 0.3, ylim = c(-0.1, 0.1), xlim = c(-1,
2), xlab = "x or gene concentration")
####
s = 1/40
r = 0.25
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "green", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
grid()
##
s = 2/40
r = 0.25
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "red", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
title(TeX("$\\frac{dx}{d\\tau} = s-rx + \\frac{x^{2}}{1+x^{2}}$"))
legend("bottomleft", border = "black", legend = c("s=0/40, r = 0.25", "s=1/40, r = 0.25",
"s=2/40, r = 0.25"), col = c("blue", "green", "red"), lty = 1:2, cex = 0.8)
grid()library(latex2exp)
x = -250:250/100
s = 0/40
r = 0.45
dxdtau = geneswitching(s = s, r = r, x = x)
plot(x, dxdtau, type = "l", col = "blue", lwd = 0.3, ylim = c(-0.1, 0.1), xlim = c(-1,
2), xlab = "x or gene concentration")
####
s = 1/40
r = 0.45
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "green", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
grid()
##
s = 2/40
r = 0.45
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "red", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
title(TeX("$\\frac{dx}{d\\tau} = s-rx + \\frac{x^{2}}{1+x^{2}}$"))
legend("bottomleft", border = "black", legend = c("s=0/40, r = 0.45", "s=1/40, r = 0.45",
"s=2/40, r = 0.45"), col = c("blue", "green", "red"), lty = 1:2, cex = 0.8)
grid()library(latex2exp)
x = -250:250/100
s = 0/40
r = 0.5
dxdtau = geneswitching(s = s, r = r, x = x)
plot(x, dxdtau, type = "l", col = "blue", lwd = 0.3, ylim = c(-0.1, 0.1), xlim = c(-1,
2), xlab = "x or gene concentration")
####
s = 1/40
r = 0.5
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "green", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
grid()
##
s = 2/40
r = 0.5
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "red", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
title(TeX("$\\frac{dx}{d\\tau} = s-rx + \\frac{x^{2}}{1+x^{2}}$"))
legend("bottomleft", border = "black", legend = c("s=0/40, r = 0.5", "s=1/40, r = 0.5",
"s=2/40, r = 0.5"), col = c("blue", "green", "red"), lty = 1:2, cex = 0.8)
grid()library(latex2exp)
x = -250:250/100
s = 0/40
r = 0.7
dxdtau = geneswitching(s = s, r = r, x = x)
plot(x, dxdtau, type = "l", col = "blue", lwd = 0.3, ylim = c(-0.1, 0.1), xlim = c(-1,
2), xlab = "x or gene concentration")
####
s = 1/40
r = 0.7
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "green", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
grid()
##
s = 2/40
r = 0.7
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "red", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
title(TeX("$\\frac{dx}{d\\tau} = s-rx + \\frac{x^{2}}{1+x^{2}}$"))
legend("bottomleft", border = "black", legend = c("s=0/40, r = 0.7", "s=1/40, r = 0.7",
"s=2/40, r = 0.7"), col = c("blue", "green", "red"), lty = 1:2, cex = 0.8)
grid()library(latex2exp)
x = -250:250/100
s = 0/40
r = 1
dxdtau = geneswitching(s = s, r = r, x = x)
plot(x, dxdtau, type = "l", col = "blue", lwd = 0.3, ylim = c(-0.1, 0.1), xlim = c(-1,
2), xlab = "x or gene concentration")
####
s = 1/40
r = 1
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "green", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
grid()
##
s = 2/40
r = 1
dxdtau = geneswitching(s = s, r = r, x = x)
lines(x, dxdtau, type = "l", col = "red", lwd = 0.3)
abline(v = 0, col = "black", lwd = 3, lty = 2)
abline(h = 0, col = "black", lwd = 3, lty = 2)
title(TeX("$\\frac{dx}{d\\tau} = s-rx + \\frac{x^{2}}{1+x^{2}}$"))
legend("bottomleft", border = "black", legend = c("s=0/40, r = 1", "s=1/40, r = 1",
"s=2/40, r = 1"), col = c("blue", "green", "red"), lty = 1:2, cex = 0.8)
grid()x = -100:100/100
den = (1 + x * x) * (1 + x * x)
norm_r = 2 * x
r = norm_r/den
norm_s = x * x * (1 - x * x)
s = norm_s/den
plot(r, s, type = "o", col = "red", lwd = 1)
abline(h = 0, col = "blue")
grid()The saddle-node bifurcation curve in (r, s) space is as follow:
\[ \begin{equation} \left\{ \begin{array}{rl} s &= \frac{x^{2}(1-x^{2})}{(1+x^{2})^{2}} \\ \\ r &= \frac{2x}{(1+x^{2})^{2}} \end{array} \right. \end{equation}\]
x = 1:100/100
den = (1 + x * x) * (1 + x * x)
norm_r = 2 * x
r = norm_r/den
norm_s = x * x * (1 - x * x)
s = norm_s/den
plot(r, s, type = "o", col = "red", lwd = 1)
abline(h = 0, col = "blue", lwd = 3, lty = 2)
abline(v = 0.5, col = "green", lwd = 3, lty = 2)
title("Saddle-node bifurcation Curve in (r,s) space for 3 real fixed points")
legend("topleft", border = "black", legend = c("Saddle-node Bifurcation Boundary",
"Bifurcation at r = 0.5", "Bifurcation at r = 0.5"), col = c("red", "blue",
"green"), lty = 1:2, cex = 0.8)
grid()Bifurcation Diagram
r = 100:500/1000
x_fixed_point_01 = (1 + sqrt(1 - 4 * r * r))/(2 * r)
x_fixed_point_02 = (1 - sqrt(1 - 4 * r * r))/(2 * r)
x_fixed_point_03 = 0
plot(r, x_fixed_point_02, type = "l", col = "blue", lwd = 0.7, ylim = c(0, 10),
ylab = "x_fixed_points")
lines(r, x_fixed_point_01, type = "l", col = "red", lwd = 0.7)
abline(h = 0, col = "green", lwd = 3, lty = 2)
title(" Bifurcation Diagram: 3 Real fixed points ")
legend("topright", border = "black", legend = c("unstable", "semistable at x = 0",
"stable"), col = c("red", "green", "blue"), lty = 1:2, cex = 0.8)
grid()