Python/生物分子物理-尝试编写一个简单的随机模拟系统以展示条件行为！
*编辑于2010年6月17日
我正在努力理解如何改进我的代码（让它更符合Python的风格）。此外，我也想写出更直观的“条件语句”，来描述生物化学中常见的场景。下面程序中的条件标准我在答案#2中解释过，但我对代码不太满意——它运行得很好，但不够明显，也不容易处理更复杂的条件场景。欢迎提出想法，欢迎评论和批评。这是我第一次在stackoverflow发帖——如果需要，请评论一下礼仪。
这段代码生成了一系列值，作为以下练习的解决方案：
“在你选择的编程语言中，实现吉莱斯皮的第一反应算法，以研究反应A--->B的时间行为，其中A转变为B只有在另一种化合物C存在时才能发生，而C与D之间动态相互转化，如下图所示的Petri网所建模。假设在反应开始时有100个A分子，1个C分子，且没有B或D。将kAB设置为0.1 s-1，kCD和kDC都设置为1.0 s-1。模拟系统在100秒内的行为。”
def sim():
# Set the rate constants for all transitions
kAB = 0.1
kCD = 1.0
kDC = 1.0
# Set up the initial state
A = 100
B = 0
C = 1
D = 0
# Set the start and end times
t = 0.0
tEnd = 100.0
print "Time\t", "Transition\t", "A\t", "B\t", "C\t", "D"
# Compute the first interval
transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC)
# Loop until the end time is exceded or no transition can fire any more
while t <= tEnd and transition >= 0:
print t, '\t', transition, '\t', A, '\t', B, '\t', C, '\t', D
t += interval
if transition == 0:
A -= 1
B += 1
if transition == 1:
C -= 1
D += 1
if transition == 2:
C += 1
D -= 1
transition, interval = transitionData(A, B, C, D, kAB, kCD, kDC)
def transitionData(A, B, C, D, kAB, kCD, kDC):
""" Returns nTransition, the number of the firing transition (0: A->B,
1: C->D, 2: D->C), and interval, the interval between the time of
the previous transition and that of the current one. """
RAB = kAB * A * C
RCD = kCD * C
RDC = kDC * D
dt = [-1.0, -1.0, -1.0]
if RAB > 0.0:
dt[0] = -math.log(1.0 - random.random())/RAB
if RCD > 0.0:
dt[1] = -math.log(1.0 - random.random())/RCD
if RDC > 0.0:
dt[2] = -math.log(1.0 - random.random())/RDC
interval = 1e36
transition = -1
for n in range(len(dt)):
if dt[n] > 0.0 and dt[n] < interval:
interval = dt[n]
transition = n
return transition, interval
if __name__ == '__main__':
sim()
4 个回答
我不太了解Gillespie算法，但我猜你已经确认程序能正确收敛到平衡状态。所以我理解你的问题是：
“这是一个能运行的物理程序，我该如何让它更符合Python的风格？”
可能更符合Python风格的做法是像下面这样：
R = [ kAB * A * C, kCD * C, kAB * A * C]
dt = [(-math.log(1-random.random())/x,i) for i,x in enumerate(R) if x > 0]
if dt:
interval,transition = min(dt)
else:
transition = None
如果你想在物理学中使用Python，我建议你学习numpy。因为numpy在处理很多问题时速度更快。所以这里有一些未经测试的numpy解决方案的部分。把以下内容添加到你程序的开头：
from numpy import log, array, isinf, zeros
from numpy.random import rand
然后你可以把TransitionData里面的内容替换成类似下面的东西：
R = array([ kAB * A * C, kCD * C, kAB * A * C])
dt = -log(1-rand(3))/R
transition = dt.argmin()
interval = dt[transition]
if isinf(interval):
transition = None
我不确定用抛出StopIteration异常代替返回None是否更符合Python风格，但这只是个小细节。
你还应该把浓度存储在一个数据结构里。如果你使用numpy，我建议你用数组。同样，你可以用一个数组dABCD来存储浓度的变化（你可能能想出更好的变量名）。把以下代码添加到你的循环外：
ABCD = array([A,B,C,D])
dABCD = zeros(3,4)
dABCD[0,0] = -1#A decreases when transition == 0
dABCD[0,1] = 1 #B increases when transition == 0
dABCD[1,2] = -1#C decreases when transition == 1
dABCD[1,3] = 1 #D increases when transition == 1
..... etc
现在你可以把主循环替换成类似下面的内容：
while t <= tEnd:
print t, '\t', transition, '\t', ABCD
transition, interval = transitionData(ABCD, kAB, kCD, kDC)
if transition != None:
t += interval
ABCD += dABCD[transition,:]
else:
break;
else:
print "Warning: Stopping due to timeout. The system did not equilibrate"
可能还有更多需要做的事情。比如，dABCD可能应该是一个稀疏数组，但我希望这些想法能给你一个好的开始。
关于简单随机模拟化学反应的数学信息：
通常，这种过程会被模拟成离散事件，每个事件发生的概率用'P'表示，这个概率是基于一个特定的速率常数'k'和在时间间隔'dt'内可能发生的事件数量'n'计算出来的，公式是P=1-e**(-kdtn)。在这里，我们忽略了每个事件的实际发生时间（大约为0），而是关注事件发生的时间间隔。对于那些熟悉组合问题或伯努利试验的人来说，1/e的出现是很容易理解的。例如，当N=K且N趋近于无穷大时，从N中不选择特定元素的概率接近于1/e。因此，在一个随机化学反应（一级反应）中，分子不发生反应（不被选择）的概率是1/e的某个幂，这个幂取决于时间间隔、速率常数以及相关的分子数量和速率常数。相反，1-(1/e)^xyz则表示任何特定分子发生反应（被选择）的概率。
在模拟方面，合理的做法是将我们的总时间间隔分成越来越小的间隔，并使用随机数生成器来预测在某个时间间隔内事件是否发生。例如，如果我们将一个事件的dt分成10个更小的间隔，那么0到0.1之间的数字表示事件发生，而0.1到1.0之间的数字则表示事件没有发生。然而，关于事件确切发生的时间仍然存在不确定性，因此我们必须将间隔做得更小，但这样做很快就会变得无效，因为这种方法的误差依然存在。
解决这个问题的方法是对上述公式的两边取自然对数（这里用‘ln’表示，在Python中默认是log()），然后解出dt，得到dt= (-ln(1-P))/(k*n)。然后随机生成概率P，从而为每个事件提供一个明确的dt。