Circadiane Tagesrhytmik
Um mal mit komplexen Systemen zu experimentieren habe ich ein System aus Differentialgleichungen über den molekularen circadianen Rhythmus der Fruchtfliege implementiert und daraus ein Video gemacht:
Hier ist der Sourcecode für die numerische Simulation im Video. Dieses Skript erzeugt eine Reihe von Bilddateien, die später mit ffmpeg zu Videos konvertiert werden. Diese Videos wiederum wurden dann mit Synfig Studio zu einer Animation zusammengestellt, aber das geht natürlich auch prima in Blender 3D.
Dieser Quellcode ist allerdings nichts für zartbesaitete und erfordert ein gewisses Grundverständnis in numpy und matplotlib.
from numpy import *
from numpy.random import uniform
from pylab import *
from copy import copy
"""
How do cells tell the time on a molecular level,
without depending on light and dark?
The circadian rhythm of the fruit fly is a good example for a moderately complicated system in biology.
It uses proteins that eventually downregulate their own transcription to achieve oscillations that can both be
entrained by light-and-dark cycles and persist in continuous darkness.
Here is a link to the publication where I took the model and parameters from:
http://www.ncbi.nlm.nih.gov/pubmed/9486845
This script will create a set of images in a subdirectory.
"""
class System:
def __init__(self, number=1):
# Setup the system for a number of "number" cells
self.number = number
self.variables = [uniform(0.5,0.6, self.number) for i in range(10)]
self.variables[0] = repeat(0.289, self.number)
self.variables[1] = repeat(1.279, self.number)
self.variables[2] = repeat(0.068, self.number)
self.variables[3] = repeat(0.076, self.number)
self.variables[4] = repeat(0.026, self.number)
self.variables[5] = repeat(0.430, self.number)
self.variables[6] = repeat(0.703, self.number)
self.variables[7] = repeat(23.042, self.number)
self.variables[8] = repeat(0.11, self.number)
self.variables[9] = repeat(0.6, self.number)
self.t = 0
def gradient(self):
# Compute gradient from the differential equations
n = 1
M_p, M_t, P_0, P_1, P_2,\
T_0, T_1, T_2, \
C, C_n= self.variables
v_sp = 0.8
v_st = 1.0
k_sp = k_st = 0.9
v_mt = 0.7
v_dp = 2.0
if self.t > 24*10:
if self.t % 24 < 12:
v_dt = 2.0
else:
v_dt = 5.0
else:
v_dt = 2.0
v_dt = 2.0
vmp = 0.8
K_mp = 0.2
K_mt = 0.2
K_d = k_dc = K_dn = 0.01
k_d = 0.01
K_dp = 0.20
K_dt = K_dp = 0.2
K_ip = K_it = 1.0
V_1p = V_1t = 8.0
V_2p = V_2t = 1.0
V_3P = V_3t = 8.0
K_1p = K_2p = K_3p = K_4p = K_1t = K_2t = K_3t = K_4t = 2
V_1p = V_1t = 8.0
V_2p = V_2t = 1.0
V_3p = V_3t = 8.0
V_4p = V_4t = 1.0
k_1 = 1.2
k_2 = 0.2
k_3 = 1.2
k_4 = 0.6
dM_p = v_sp * (K_ip / (K_ip + C_n)) \
- vmp * (M_p/(K_mp+M_p)) \
- K_d*M_p
dP_0 = k_sp*M_p \
- V_1p * (P_0/ (K_1p + P_0)) \
+ V_2p * (P_1 / (K_2p+P_1)) \
- K_d * P_0
dP_1 = V_2p * (P_0/(K_1p+P_0)) \
- V_2p* (P_1 / (K_2p + P_1)) \
- V_3p* (P_1 / (K_3p + P_1)) \
+ V_4p* (P_2 / (K_4p + P_2)) \
- K_d * P_1
dP_2 = V_3p * (P_1/(K_3p+P_1)) \
- V_4p* (P_2 / (K_4p + P_2)) \
- k_3*P_2*T_2 \
+ k_4 * C \
- v_dp * (P_2 /(K_dp + P_2)) \
- k_d * P_2
dM_t = v_st * (K_it / (K_it+ C_n)) \
- v_mt * (M_t / (K_mt+ M_t)) \
- k_d * M_t
dT_0 = k_st * M_t \
- V_1t * ( T_0 / (K_1t+T_0)) \
+ V_2t * ( T_1 / (K_2t+T_1)) \
- k_d * T_0
dT_1 = V_1t * (T_0 / (K_1t + T_0)) \
- V_2t * (T_1 / (K_2t + T_1)) \
- V_3t * (T_1 / (K_3t + T_1)) \
+ V_4t * (T_2 / (K_4t + T_2)) \
- k_d * T_1
dT_2 = V_3t * (T_1 / (K_3t + T_1)) \
- V_4t * (T_2 / (K_4t + T_2)) \
- k_3 * P_2 * T_2 \
+ k_4 * C \
- v_dt * (T_2 / (K_dt + T_2)) \
- k_d * T_2
dC = k_3 * P_2 * T_2 - k_4 * C - k_1*C + k_2 * C_n - k_dc*C
dC_n = k_1 * C - k_2 * C_n - K_dn * C_n
return dM_p, dM_t, dP_0, dP_1, dP_2,\
dT_0, dT_1, dT_2, \
dC, dC_n
def iterate(self):
# Perform Runge-Kutta integration with the gradient above
h = 0.02
M_p, M_t, P_0, P_1, P_2,\
T_0, T_1, T_2, \
C, C_n= self.variables
dM_p, dM_t, dP_0, dP_1, dP_2,\
dT_0, dT_1, dT_2, \
dC, dC_n = self.gradient()
M_p += h*dM_p
P_0 += h*dP_0
P_1 += h*dP_1
P_2 += h*dP_2
M_t += h*dM_t
T_0 += h*dT_0
T_1 += h*dT_1
T_2 += h*dT_2
C += h*dC
C_n += h*dC_n
self.t += h / 4.0
self.variables = M_p, M_t, P_0, P_1, P_2,\
T_0, T_1, T_2, \
C, C_n
system = System()
values = list([list() for i in range(10)])
time_axis = []
fps = 24.0 # Frames per second
hours_per_second = 24.0 / 4.0 # Hours per second in video
start_time = 24*2 # Start video after two simulated days
frame = 0
for i in range(1, int(24*100*4/ 0.02)):
system.iterate()
if i % 240==0:
print system.t
if (i % 50 ==0) and (system.t> 0):
for j in range(10):
values[j].append(system.variables[j][0])
time_axis.append(system.t)
if system.t> start_time:
if system.t - start_time> frame / fps * hours_per_second:
frame+=1
print frame, system.t
M_p, M_t, P_0, P_1, P_2,\
T_0, T_1, T_2, \
C, C_n= values
start_index = 0 #len(time_axis)- int(24 / 0.02 / 4.0)
cla()
subplot(1,1,1)
plot(time_axis[start_index:], P_0[start_index:], "r", linewidth=5)
plot(time_axis[start_index:], P_1[start_index:], "g", linewidth=5)
plot(time_axis[start_index:], P_2[start_index:], "b", linewidth=5)
ylim (0,0.2)
xlim(system.t - 24*1, system.t)
savefig("anim/P_i%05i.png" % (frame),dpi=25)
cla()
subplot(1,1,1)
plot(time_axis[start_index:], M_p[start_index:], "r", linewidth=5)
plot(time_axis[start_index:], M_t[start_index:], "g", linewidth=5)
xlim(system.t - 24*1, system.t)
savefig("anim/rna%05i.png" % (frame),dpi=40)
cla()
subplot(1,1,1)
plot(time_axis[start_index:], T_0[start_index:], "r", linewidth=5)
plot(time_axis[start_index:], T_1[start_index:], "g", linewidth=5)
plot(time_axis[start_index:], array(T_2[start_index:])/ 4.0, "b", linewidth=5)
xlim(system.t - 24*1, system.t)
savefig("anim/T_i%05i.png" % (frame),dpi=40)
cla()
subplot(1,1,1)
plot(time_axis[start_index:], C[start_index:], "r", linewidth=5)
plot(time_axis[start_index:], array(C_n[start_index:]), "b", linewidth=5)
xlim(system.t - 24*1, system.t)
savefig("anim/C%05i.png" % (frame),dpi=40)