1

对曲线拟合/优化有疑问。我有三个耦合的 ODE，它们描述了与消失的底物和正在形成的两种产物的生化反应。我找到了帮助我创建代码来解决 ODE 的示例（如下）。现在我想优化未知速率常数（k、k3 和 k4）以适应实验数据 P，它是来自产品 y[1] 的信号。这样做最简单的方法是什么？谢谢。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Experimental data 
P = [29.976,193.96,362.64,454.78,498.42,517.14,515.76,496.38,472.14,432.81,386.95,
352.93,318.93,279.47,260.19,230.92,202.67,180.3,159.09,137.31,120.47,104.51,99.371,
89.606,75.431,67.137,58.561,55.721]

# Three coupled ODEs
def conc (y, t) : 
    a1 = k * y[0] 
    a2 = k2 * y[0]
    a3 = k3 * y[1]
    a4 = k4 * y[1]
    a5 = k5 * y[2]
    f1 = -a1 -a2
    f2 = a1 -a3 -a4
    f3 = a4 -a5
    f = np.array([f1, f2, f3])
    return f


# Initial conditions for y[0], y[1] and y[2]
y0 = np.array([50000, 0.0, 0.0])

# Times at which the solution is to be computed.
t = np.linspace(0.5, 54.5, 28)

# Experimentally determined parameters.
k2 = 0.071
k5 = 0.029

# Parameters which would have to be fitted
k = 0.002
k3 = 0.1
k4 = 0.018

# Solve the equation
y = odeint(conc, y0, t)

# Plot data and the solution.
plt.plot(t, P, "bo")
#plt.plot(t, y[:,0]) # Substrate
plt.plot(t, y[:,1]) # Product 1
plt.plot(t, y[:,2]) # Product 2
plt.xlabel('t')
plt.ylabel('y')
plt.show()
4

1 回答 1

1

编辑：我对代码进行了一些更改，以展示如何适应所有 ODE 的实验数据。

像这样：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from scipy.odr import Model, Data, ODR
# Experimental data 
P = [29.976,193.96,362.64,454.78,498.42,517.14,515.76,496.38,472.14,432.81,386.95,
352.93,318.93,279.47,260.19,230.92,202.67,180.3,159.09,137.31,120.47,104.51,99.371,
89.606,75.431,67.137,58.561,55.721]

# Times at which the solution is to be computed.
t = np.linspace(0.5, 54.5, 28)


def coupledODE(beta, x):
    k, k3, k4 = beta

    # Three coupled ODEs
    def conc (y, t) : 
        a1 = k * y[0] 
        a2 = k2 * y[0]
        a3 = k3 * y[1]
        a4 = k4 * y[1]
        a5 = k5 * y[2]
        f1 = -a1 -a2
        f2 = a1 -a3 -a4
        f3 = a4 -a5
        f = np.array([f1, f2, f3])
        return f


    # Initial conditions for y[0], y[1] and y[2]
    y0 = np.array([50000, 0.0, 0.0])

    # Experimentally determined parameters.
    k2 = 0.071
    k5 = 0.029

    # Parameters which would have to be fitted
    #k = 0.002
    #k3 = 0.1
    #k4 = 0.018

    # Solve the equation
    y = odeint(conc, y0, x)

    return y[:,1]
    # in case you are only fitting to experimental findings of ODE #1

    # return y.ravel()
    # in case you have experimental findings of all three ODEs

data = Data(t, P)
# with P being experimental findings of ODE #1

# data = Data(np.repeat(t, 3), P.ravel())
# with P being a (3,N) array of experimental findings of all ODEs

model = Model(coupledODE)
guess = [0.1,0.1,0.1]
odr = ODR(data, model, guess)
odr.set_job(2)
out = odr.run()
print out.beta
print out.sd_beta

f = plt.figure()
p = f.add_subplot(111)
p.plot(t, P, 'ro')
p.plot(t, coupledODE(out.beta, t))
plt.show()

如果您使用的是基于 scipy 的交互式曲线拟合程序peak-o-mat ( http://lorentz.sf.net )，您可以添加您的 ODE 模型并将其保存到 userfunc.py（请参阅自定义部分在文档中）：

import numpy as np
from scipy.integrate import odeint
from peak_o_mat import peaksupport as ps

def coupODE(x, k, k3, k4):

    # Three coupled ODEs
    def conc (y, t) : 
        a1 = k * y[0] 
        a2 = k2 * y[0]
        a3 = k3 * y[1]
        a4 = k4 * y[1]
        a5 = k5 * y[2]
        f1 = -a1 -a2
        f2 = a1 -a3 -a4
        f3 = a4 -a5
        f = np.array([f1, f2, f3])
        return f


    # Initial conditions for y[0], y[1] and y[2]
    y0 = np.array([50000, 0.0, 0.0])

    # Times at which the solution is to be computed.
    #t = np.linspace(0.5, 54.5, 28)

    # Experimentally determined parameters.
    k2 = 0.071
    k5 = 0.029

    # Parameters which would have to be fitted
    #k = 0.002
    #k3 = 0.1
    #k4 = 0.018

    # Solve the equation
    y = odeint(conc, y0, x)
    print y
    return y[:,1]

ps.add('ODE',
          func='coupODE(x,k,k3,k4)',
          info='thre coupled ODEs',
          ptype='MISC')

您需要将数据准备为包含时间和实验数据两列的文本文件。将数据导入 peak-o-mat，输入“ODE”作为拟合模型，为 k、k3、k4 选择适当的初始参数并点击“拟合”。

于 2013-07-16T16:20:28.987 回答