toggleSwitch.py

# encoding: utf-8
from pylab import *
from sympy import *

def XDiff(x, y, a):
    return a / (1.0+ y **2) - x #Toggle switchモデル

def YDiff(x, y, a):
    return a / (1.0+ x **2) - y

def nullcline(range, a): #共通ヌルクライン
    return a / (1.0+ range **2)

def checkStablePoint(a, x, y): #det(J)が正かどうか確かめる
    if 4* x * y * (a **2) <= ((1.0+ x **2) **2) * ((1.0+ y **2) **2):
        return "ko" #安定
    else:
        return "wo" #不安定

x  = arange(0, 5.0, 0.1) #タンパク質濃度の定義域
a = [2.0, 4.0] #抑制速度
Nullcline = [] #共通ヌルクラインを作成する

for i in a:
    Nullcline.append(nullcline(x, i))
"""#1つ目の固定点解析
figure()
xlim(0, 4)
ylim(0, 4)
plot(Nullcline[0], x)
plot(x, Nullcline[0])
plot(1.0, nullcline(1.0, a[0]), checkStablePoint(a[0], 1.0, nullcline(1.0, a[0])), markersize =10)
#2つ目の固定点解析
figure()
xlim(0, 4)
ylim(0, 4)
plot(Nullcline[1], x)
plot(x, Nullcline[1])
xx = Symbol('xx') #交点を求めるため三次方程式を解く (sympy.solve)
vertex = solve(xx **5- a[1] * xx **4+2* xx **3-2* a[1] * xx **2  + (a[1] **2 + 1.0) * xx - a[1], xx)
print vertex
plot(vertex[0].evalf(), nullcline(vertex[0].evalf(), a[1]), checkStablePoint(a[1], vertex[0].evalf(), nullcline(vertex[0].evalf(), a[1])), markersize =10)
plot(vertex[2].evalf(), nullcline(vertex[2].evalf(), a[1]), checkStablePoint(a[1], vertex[2].evalf(), nullcline(vertex[2].evalf(), a[1])), markersize =10)
plot(vertex[3].evalf(), nullcline(vertex[3].evalf(), a[1]), checkStablePoint(a[1], vertex[3].evalf(), nullcline(vertex[3].evalf(), a[1])), markersize =10)

for i in range(2):
    figure()
    xlim(0, 4)
    ylim(0, 4)
    plot(Nullcline[i], x)
    plot(x, Nullcline[i])
    
    for j in arange(0, 4.1, 0.2):
        for k in arange(0, 4.1, 0.2):
            quiver(j, k, XDiff(j, k, a[i]), YDiff(j, k, a[i]))
"""
for i in range(2):
    figure()
    xlim(0, 4)
    ylim(0, 4)
    plot(Nullcline[i], x) #xヌルクライン
    plot(x, Nullcline[i]) #yヌルクライン

show() #全figure描画