曲線擬合耦合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'。