Gillespie算法的Python简单实现（实例）

文章目录

- Gillespie算法的Python简单实现（实例）
前言
一、Gillespie是什么？
二、题目
三、代码
- 1.引入库
- 2.类定义
- 3.设置参量
- 4.计算结果并输出拟合图像
部分结果
总结

前言

一位化院的朋友的期末作业，寻求帮助，因此了解学习了一下gillespie算法，并简单实现了该算法，实例化了mRNA翻译蛋白质的过程。

一、Gillespie是什么？

个人理解，Gillespie算法是一种将连续的化学反应转换为离散步骤模拟的算法。

二、题目

简单来说就是利用Gillespie算法拟合mRNA翻译为氨基酸链的过程，并带入多组参量计算氨基酸链生成速度。

三、代码

本文代码部分段来自python简单实现gillespie模拟

1.引入库

代码如下：

import numpy as np
from scipy.special import comb  # 科学计算
import matplotlib.pyplot as plt  # 绘图

2.类定义

封装了反应类和系统类，主要来源上述链接文章：

class Reaction:  # 封装的类，代表每一个化学反应def __init__(self, rate=0., num_lefts=None, num_rights=None):self.rate = rate  # 反应速率assert len(num_lefts) == len(num_rights)self.num_lefts = np.array(num_lefts)  # 反应前各个反应物的数目self.num_rights = np.array(num_rights)  # 反应后各个反应物的数目self.num_diff = self.num_rights - self.num_lefts  # 改变数def combine(self, n, s):  # 算组合数return np.prod(comb(n, s))def propensity(self, n):  # 算反应倾向函数return self.rate * self.combine(n, self.num_lefts)class System:  # 封装的类，代表多个化学反应构成的系统def __init__(self, num_elements):assert num_elements > 0self.num_elements = num_elements  # 系统内的反应物的类别数self.reactions = []  # 反应集合def add_reaction(self, rate=0., num_lefts=None, num_rights=None):assert len(num_lefts) == self.num_elementsassert len(num_rights) == self.num_elementsself.reactions.append(Reaction(rate, num_lefts, num_rights))def evolute(self, steps, inits=None):  # 模拟演化self.t = [0]  # 时间轨迹，t[0]是初始时间if inits is None:self.n = [np.ones(self.num_elements)]else:self.n = [np.array(inits)]  # 反应物数目，n[0]是初始数目for i in range(steps):A = np.array([rec.propensity(self.n[-1])for rec in self.reactions])  # 算每个反应的倾向函数A0 = A.sum()A /= A0  # 归一化得到概率分布t0 = -np.log(np.random.random()) / A0  # 按概率选择下一个反应发生的间隔# print(np.mean(t0))self.t.append(self.t[-1] + t0)d = np.random.choice(self.reactions, p=A)  # 按概率选择其中一个反应发生self.n.append(self.n[-1] + d.num_diff)

3.设置参量

L = 10
V0 = [0.2, 0.4, 0.6, 0.8, 1.0]
V1 = [0.2, 0.4, 0.6, 0.8, 1.0]
vd = 1

4.计算结果并输出拟合图像

for m in range(5):for n in range(5):v0 = V0[m]v1 = V1[n]num_elements = 3system = System(num_elements)system.add_reaction(v0, [0, 0, 0], [1, 0, 0])system.add_reaction(vd / L, [1, 0, 0], [0, 1, 0])system.add_reaction(v1, [0, 1, 0], [0, 0, 1])# 创建反应系统system.evolute(100000)  # 执行100000步x = system.ty1 = [i[1] for i in system.n]y2 = [i[2] for i in system.n]xp = x[1:]y2p = y2[1:]v2 = np.mean(np.divide(y2p, xp))print(v2)plt.clf()plt.plot(x, y1)  # 正在转录中的mRNA轨线图plt.xlim(0, x[-1] + 1)plt.xlabel('t/s')plt.ylabel('mRNA')plt.savefig('mRNA%d%d.png' % (m, n))# plt.show()plt.clf()plt.plot(x, y2)  # 蛋白质的轨线图plt.xlim(0, x[-1] + 1)plt.xlabel('t/s')plt.ylabel('protein')plt.text(xp[-1] / 2, 16000, 'v=%s' % v2, fontsize=12)plt.text(xp[-1] / 16, 30000, 'L=%s' % L, fontsize=12)plt.text(xp[-1] / 16, 27500, 'λ0=%s' % v0, fontsize=12)plt.text(xp[-1] / 16, 25000, 'λ1=%s' % v1, fontsize=12)plt.text(xp[-1] / 16, 22500, 'λd=%s' % vd, fontsize=12)plt.savefig('protein%d%d.png' % (m, n))# plt.show()

部分结果

可以看到，与λ0基本上线性相关，与λ1几乎无关，总结原因，主要是因为决速步在第一步，故与λ0正相关，而λ1前置步骤产率很慢，本文所取的λ1均已饱和，故几乎无关。拟合结果很贴合。

总结

学习了gillespie算法，总结了拟合结果及原因，练习了代码能力。

Gillespie算法的Python简单实现（实例）相关推荐

python简单程序实例-Python简单基础小程序的实例代码
1 九九乘法表 for i in range(9):#从0循环到8 i += 1#等价于 i = i+1 for j in range(i):#从0循环到i j += 1 print(j,'*',i, ...
python简单程序实例-python简单实例训练(21~30)
注意:我用的python2.7,大家如果用Python3.0以上的版本,请记得在print()函数哦!如果因为版本问题评论的,不做回复哦!! 21.题目:将一个正整数分解质因数.例如:输入90,打印出 ...
鲸鱼算法与Python简单可视化测试（1）
文章目录 1. 鲸鱼算法原理 2. 鲸鱼算法的python实现 2.1 问题描述 2.2 初始化 2.3 可视化 3. 可视化结果参考文献 1. 鲸鱼算法原理鲸鱼算法是一种根据鲸鱼围捕猎物行为抽象 ...
python简单程序实例-python实现的简单窗口倒计时界面实例
本文实例讲述了python实现的简单窗口倒计时界面.分享给大家供大家参考.具体分析如下: 下面的代码通过Tkinter制作windows窗口界面,然后时间了一个简单的倒计时功能,代码可以直接运行 # ...
python简单程序实例-python简单项目实例
语言多元化是PayPal编程文化中一个重要的组成部分.在C++和Java长期流行的同时,更多的团队选择了Jva和Scala.同时,Braintree的收购也引入了一个久经世故的Ruby社区.Pytho ...
python简单游戏实例_Python实现的简单算术游戏实例
本文实例讲述了Python实现的简单算术游戏.分享给大家供大家参考.具体实现方法如下: #!/usr/bin/env python from operator import add, sub from ...
EM算法及python简单实现
目录 1.摘要 2.EM算法简介 3.预备知识 3.1贝叶斯公式 3.2 极大似然估计 3.2.1似然函数 3.2.2问题描述 3.2.3极大似然函数估计值的求解步骤 3.3Jensen不等式(琴生不 ...
python简单程序实例-python下10个简单实例代码
注意:我用的python2.7,大家如果用Python3.0以上的版本,请记得在print()函数哦!如果因为版本问题评论的,不做回复哦!!! 1.题目:有1.2.3.4个数字,能组成多少个互不相同且 ...
OneR算法的Python简单实现
OneR算法就是,在已有数据中,根据具有相同特征值的个体最可能属于哪个类别进行分类.即取效果最好的那个特征进行分类. #-*- coding=utf-8 -*- # import numpy as n ...

Gillespie算法的Python简单实现（实例）

Gillespie算法的Python简单实现（实例）

文章目录

前言

一、Gillespie是什么？

二、题目

三、代码

1.引入库

2.类定义

3.设置参量

4.计算结果并输出拟合图像

部分结果

总结

Gillespie算法的Python简单实现（实例）相关推荐

最新文章

热门文章