曲线拟合耦合ODEs

Question

在曲线拟合/优化上有一个问题。我有三个偶联的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() 

Answer 1

编辑：我为了展示如何适合所有微分方程的实验数据做了一些修改的代码。

像这样：

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() 

如果你正在使用高峰邻垫（http://lorentz.sf.net），这是基于SciPy的互动式曲线拟合程序，您可以添加您的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选择合适的初始参数并点击'Fit'。