非線形関数フィッティング¶
酵素反応は一般に基質濃度を[S]として、 $$-\frac{d[S]}{dt}=k[S]$$ として表されます。この解を解析的に求めることは簡単で、時間と基質濃度の関係は以下の式で与えられます。 $$[S] = Ae^{-kt}\qquad (\text{Aは積分定数})$$ この式は非線形関数ですが、誤差関数を指定して（今回は二乗誤差）パラメータを最尤推定すれば曲線フィッティングが可能です。
一応今回の実験のプロトコルや原理をかいつまんで書いておくので、非線形関数フィッティングだけ見たい場合はここへどうぞ。
原理¶
リゾチームは、真正細菌の細胞壁を構成する多糖類を加水分解する酵素である。 Wikipedia
この時の反応は酵素反応であり、細胞溶液の濃度は、式(2)の曲線に乗ると考えられます。
なお、直接濃度を測ることは難しいので、$\mathrm{OD}_{600}$ の値を用います。なお、$\mathrm{OD}_{600}$ と濃度の関係は $\mathrm{OD}_{600}=a[S]+b$ という関係で表せることがわかっているので、以下の式で表されます。
$$\begin{aligned}
\mathrm{OD}_{600}
&= aAe^{-kt} + b\\
&= Ce^{-kt} + b
\end{aligned}$$
プロトコル¶
- 0.3 mg/mL菌液(Mycrococcus lysodeikticus)25 mLを分注した。
- リゾチームを500 mM酢酸ナトリウムpH 4.5に溶かし、0.1 mg/mLを作製した。
- 菌液をよく撹拌してから1.4 mLを光学セルに入れ、600 nmの吸光度を測定した。これを、t=0 秒後の吸光度の値とした。
- 0.1 mg/mLのリゾチーム溶液を0.1 mL入れて反応を開始した。リゾチームを添加して120 秒後までの吸光度を10 秒毎に測定した。なお、この時毎回の測定の間にピペッティングをした。
結果¶
| t[s] | 0 | 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | 100 | 110 | 120 | 130 | 140 | 150 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| $\mathrm{OD}_{600}$ | 0.939 | 0.677 | 0.576 | 0.508 | 0.451 | 0.419 | 0.381 | 0.362 | 0.339 | 0.321 | 0.299 | 0.282 | 0.275 | 0.276 | 0.267 | 0.252 |
import numpy as np
import matplotlib.pyplot as plt
time = [0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150]
OD = [0.939, 0.677, 0.576, 0.508, 0.451, 0.419, 0.381, 0.362, 0.339, 0.321, 0.299, 0.282, 0.275, 0.276, 0.267, 0.252]
plt.plot(time, OD)
plt.scatter(time, OD, label="data")
plt.title("The relationship between time and $\mathrm{OD}_{600}$")
plt.xlabel("time [s]")
plt.ylabel("$\mathrm{OD}_{600}$")
plt.legend()
plt.show()
from scipy.optimize import curve_fit
以下のパラメトリック関数を定義します。 $$\mathrm{OD}_{600} = Ce^{-kt} + b$$
def func(t,C,k,b): # 順番に注意。t(xdata)が最初。
return C * np.exp(-k*t) + b
# フィッティング
popt, _ = curve_fit(func, time, OD, method='dogbox')
trf（Trust Region Reflective algorithm, 信頼領域法）¶
反復法の１つで、基本的な概念は、シンプルな関数 $q$ で 点 $\mathbf{x}$ の近傍 $N$ （信頼領域）での目的関数 $f$ の挙動を近似するというもの。 $$\underset{\mathbf{s}}{\min}\{q(\mathbf{s}), \mathbf{s}\in N\}$$ ここで、$f(\mathbf{x}+\mathbf{s})<f(\mathbf{x})$ の場合、現在の点が $\mathbf{x}+\mathbf{s}$ に更新される。そうでない場合は現在の点は変更されず、信頼領域 $N$ は縮小され、テストステップの計算が繰り返される。なお、標準的な信頼領域法では、
- 二次近似 $q$ は、$\mathbf{x}$ における $f$ のテイラー近似のはじめの $2$ 項によって決められる。
- 近傍 $N$ は球体または楕円体。
dogbox（dogleg algorithm with rectangular trust regions）¶
信頼領域法の一種。経路に沿って、信頼領域（rectangula）内で目的関数 $f$ を最小化する $\mathbf{x}$ を見つける。
lm（Levenberg-Marquardt algorithm）※default¶
ニュートン法でヘッセ行列・またその逆行列が必要となるが、それらを用いるのは大変。
→ ヤコビ行列の積でヘッセ行列を近似しよう（ガウス・ニュートン法）
→ 最急降下法と組み合わせよう（初期は最急降下法、近傍でガウス・ニュートン法に切り替えて収束速度を上げる？）
optimal_func = lambda t: func(t,popt[0],popt[1],popt[2])
X = np.linspace(min(time), max(time), 1000)
Y = np.vectorize(optimal_func)(X)
plt.plot(X,Y, label="$\mathrm{OD}_{600} = Ce^{-kt}+b$")
plt.scatter(time, OD, label="data")
plt.title("The relationship between time and $\mathrm{OD}_{600}$")
plt.xlabel("time [s]")
plt.ylabel("$\mathrm{OD}_{600}$")
plt.legend()
plt.show()
popt[0],popt[1],popt[2]
各種パラメータを代入すると、以下の式になった。 $$\mathrm{OD}_{600} = 0.62e^{-0.32t}+0.27$$
print("RMSE: {:.3f}".format(np.sqrt(np.sum((OD-np.vectorize(optimal_func)(time))**2)/len(OD))))
補正¶
OD1 = [od*(14/15) if i==0 else od for i,od in enumerate(OD)]
# フィッティング
popt1, _ = curve_fit(func, time, OD1, method='dogbox')
optimal_func1 = lambda t: func(t,popt1[0],popt1[1],popt1[2])
Y1 = np.vectorize(optimal_func1)(X)
plt.plot(X,Y1, label="$\mathrm{OD}_{600} = Ce^{-kt}+b$")
plt.scatter(time, OD1, label="data")
plt.title("The relationship between time and $\mathrm{OD}_{600}$")
plt.xlabel("time [s]")
plt.ylabel("$\mathrm{OD}_{600}$")
plt.legend()
plt.show()
popt1[0],popt1[1],popt1[2]
各種パラメータを代入すると、以下の式になった。 $$\mathrm{OD}_{600} = 0.59e^{-0.028t}+0.26$$
print("RMSE: {:.3f}".format(np.sqrt(np.sum((OD1-np.vectorize(optimal_func1)(time))**2)/len(OD1))))
補正②¶
続いて考えられるものとして、分光光度計の測定時間の誤差が挙げられます。というのも、分光光度計を押してから計測に数秒かかっているので、どのタイミングで測られたのかわかりません。（TAの方に聞いてもわかりませんでした笑）
誤差があるのはt=0 秒後とt=10 秒後の間だけなので、初期値を取り除いてみます。
time2 = [t for i,t in enumerate(time) if i!=0]
OD2 = [od for i,od in enumerate(OD) if i!=0]
# フィッティング
popt2, _ = curve_fit(func, time2, OD2, method="dogbox")
optimal_func2 = lambda t: func(t,popt2[0],popt2[1],popt[2])
X2 = np.linspace(min(time2), max(time2), 1000)
Y2 = np.vectorize(optimal_func2)(X2)
plt.plot(X2,Y2, label="$\mathrm{OD}_{600} = Ce^{-kt}+b$")
plt.scatter(time2, OD2, label="data")
plt.title("The relationship between time and $\mathrm{OD}_{600}$")
plt.xlabel("time [s]")
plt.ylabel("$\mathrm{OD}_{600}$")
plt.legend()
plt.show()
print("RMSE: {:.3f}".format(np.sqrt(np.sum((OD2-np.vectorize(optimal_func2)(time2))**2)/len(OD2))))
各種パラメータを代入すると、以下の式になった。 $$\mathrm{OD}_{600} = 0.53e^{-0.022t}+0.24$$
ここで、初期値以外の点を利用してフィッティングさせた関数を利用して、何秒誤差があったのかを調べます。解析的に逆関数を求めることができるので、上のフィッティング式の逆関数を求めます。
$$t = \frac{1}{k}\ln\frac{C}{\mathrm{OD}_{600}-b}$$def inv_func2(od,C,k,b):
return 1/k * np.log(C/(od-b))
inv_optimal_func2 = lambda od: inv_func2(od,popt2[0],popt2[1],popt2[2])
# この値をとるtを求めたい。
OD[0]
optimal_func2(inv_optimal_func2(OD[0]))
# 誤差がある…値が小さいから？？関数ミスった？？
# パラメータ。確かに値は小さい。
popt2[0],popt2[1],popt2[2]
inv_optimal_func2(OD[0])
# -12秒は流石に嘘。