HIT 系统生物学实验一：定量模拟生化反应参考代码

最新推荐文章于 2024-07-23 14:36:35 发布

huiyiiiiii

最新推荐文章于 2024-07-23 14:36:35 发布

阅读量285

点赞数 6

文章标签： python 算法

本文链接：https://blog.csdn.net/huiyiiiiii/article/details/138274781

版权

本文介绍了如何利用Python和Matlab中的龙格-库塔算法模拟生化反应，通过米氏方程探讨底物S、酶E、中间产物ES和产物P随时间的变化，以及产物生成速率与底物浓度的关系。

摘要由CSDN通过智能技术生成

实验目的：

深入理解米氏方程

掌握使用龙格－库塔算法编程定量模拟生化反应

实验内容：

1：应用质量作用定律和原子守恒定律推导米氏方程

2：模拟仿真底物S、酶E、中间产物ES、产物P的与时间关系曲线

3：模拟仿真产物P生成速率与底物浓度之间关系曲线

定量模拟生化反应初始参数为：

k1=1000（1/秒.摩尔）=1（1/秒.毫摩尔）,

k_1=1（1/秒）;k2=0.05（1/秒）;

E0=0.5mM(毫摩尔）=0.0005M(摩尔）;

提示仿真算法采用Matlab中的龙格-库塔法ode23函数：

原理及算法：

由质量作用定律得：

由原子守恒定律得：

故：

其中：

实验环境：

python3.9 matplotlib3.6.2 pycharm64

参考代码：

## funcX：计算物质X的导数
## y：物质X	t：步长	l_:传入的一系列参数
def funcP(y, t, l_):
    k2 = l_[0]
    ES = l_[1]
    return k2 * ES
def funcS(y, t, l_):
    k1 = l_[0]
    E = l_[1]
    k_1 = l_[2]
    ES = l_[3]
    return -k1 * E * y + k_1 * ES
def funcES(y, t, l_):
    k1 = l_[0]
    E = l_[1]
    S = l_[2]
    k_1 = l_[3]
    k2 = l_[4]
    return k1 * E * S - (k_1 + k2) * y
def funcE(y, t, l_):
    k1 = l_[0]
    S = l_[1]
    k_1 = l_[2]
    k2 = l_[3]
    ES = l_[4]
return -k1 * y * S + (k_1 + k2) * ES
## 龙格－库塔算法
## y:物质	x：时间 dt：
def runge_kutta(y, x, dt, f, l_):
    k1 = dt * f(y, x, l_)
    k2 = dt * f(y + 0.5 * k1, x + 0.5 * dt, l_)
    k3 = dt * f(y + 0.5 * k2, x + 0.5 * dt, l_)
    k4 = dt * f(y + k3, x + dt, l_)
    return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.
## 产物生成速率与底物浓度关系函数：
def fun(k1, k_1, k2, E, S):
return (k2 * E * S) / (S + ((k_1 + k2) / k1))
## main 函数
if __name__ == '__main__':
## 初始化各物质浓度与反应参数
    t = 0.
    dt = 0.001
    k1 = 1e3
    k_1 = 1
    k2 = 0.05
    E = 0.5e-3
    E_total = E
    P = 0.
    ES = 0.
    S = 1e-3
    Es, Ss, ESs, Ps, ts = [], [], [], [], []
##  联立迭代计算
    while t <= 300:
        P0 = runge_kutta(P, t, dt, funcP, [k2, ES])
        S0 = runge_kutta(S, t, dt, funcS, [k1, E, k_1, ES])
        ES0 = runge_kutta(ES, t, dt, funcES, [k1, E, S, k_1, k2])
        E0 = runge_kutta(E, t, dt, funcE, [k1, S, k_1, k2, ES])
        P = P0
        S = S0
        ES = ES0
        E = E0
        t += dt
        Ps.append(P)
        Ss.append(S)
        ESs.append(ES)
        Es.append(E)
        ts.append(t)
##  输出底物S、酶E、中间产物ES、产物P的与时间关系曲线
    plt.plot(ts, Ps, label='P')
    plt.plot(ts, Ss, label='S')
    plt.plot(ts, ESs, label='ES')
    plt.plot(ts, Es, label='E')
    plt.legend()
    plt.show()
##  计算输出产物P生成速率与底物浓度关系曲线
    Vs, Ss = [], []
    for i in range(5000):
        Vs.append(fun(k1, k_1, k2, E_total, i / 100000))
        Ss.append(i / 100000)
    plt.plot(Ss, Vs, label='v')
    plt.axhline(2.5 * pow(10, -5), color="red", linestyle='--')
    plt.plot([0, 0.00105], [1.25 * pow(10, -5), 1.25 * pow(10, -5)], linestyle='--', color="red")
    plt.plot([0.00105, 0.00105], [0, 1.25 * pow(10, -5)], linestyle='--', color="red")
    plt.xlim(0, 0.05)
    plt.ylim(0, 2.6e-5)
    plt.show()