Michaelis-Menten酵素動力學一般解：Lambert W function

Michaelis-Menten速率積分式

在描述一個酵素催化反應的速率時最典型的方法就是使用Michaelis-Menten表達式

$$\frac{d[A]}{dt} =- \frac{V_{max}[A]}{[A]+K_{M}} $$

其中在基質濃度$[A]$很大的時候A的消耗速率能達到最大值$V_{max}$
而$K_{M}$則表達了要多高的濃度$[A]$才可以達到最大消耗速率的一半$\frac{V_{max}}{2}$

在學習基本的化學動力學時我們都聽過零級，一級，以及二級反應
利用積分的方式我們可以得到$[A]$隨著時間t消耗的關係式
這樣我們便能藉由動力學實驗配合積分速率式得到相應的速率常數
然而我們能夠用同樣的方式處理Michaelis-Menten速率式嗎？
答案是肯定的

首先我們需要將$[A]$和時間項t分開至等號兩側並積分

$$\int_{[A]_{0}}^{[A]} 1+ \frac{K_{M}}{[A]} d[A] =  -\int_0^t V_{max}dt $$

積分後整理可以得到

 $$ln([A])+ \frac{[A]}{K_{M}} =ln([A]_{0})+ \frac{[A]_{0}-V_{max}t}{K_{M}} $$

如果想得到$[A](t)$的解析解乍看之下難以辦到

我們這時需要引入Lambert W function的概念

Lambert W function (朗伯W函數）

Lambert W function (或稱為Omega function)有以下的定義：

$$y e^{y}=x \Longleftrightarrow y=W(x) $$

注意在一般實數應用裡$x \geq  -\frac{1}{e} $函數才會有定義
這個特殊函數可以被來在解許多指數方程式
只要能夠化成 $ye^{y}=x$ 的形式

在原先的問題中為了得到指數項
只需要將非對數項用$z=ln( e^{z}) $表示即可

$$ln([A] e^{\frac{[A]}{K_{M}}}) =ln([A]_{0} e^{\frac{[A]_{0}-V_{max}t}{K_{M}}} ) $$

將兩邊自然對數拿掉可得

$$[A] e^{\frac{[A]}{K_{M}}} =[A]_{0} e^{\frac{[A]_{0}-V_{max}t}{K_{M}}}  $$

為了湊出Lambert W function的形式
指數和指數前的係數必須一樣
所以可以考慮兩邊同除$K_{M}$

$$\frac{[A]}{K_{M}} e^{\frac{[A]}{K_{M}}} =\frac{[A]_{0} }{K_{M}} e^{\frac{[A]_{0}-V_{max}t}{K_{M}}}$$

如此以來最後的解就是

$$\frac{[A]}{K_{M}}  =W(\frac{[A]_{0} }{K_{M}} e^{\frac{[A]_{0}-V_{max}t}{K_{M}}})$$

利用無因次化的$[A]$對無因次化的t作圖，可得到在不同初始濃度下基質對時間的變化

利用積分速率式求取未知參數

如果我們現在有一組基質濃度對時間的動力學數據
可以用MATLAB的fitnlm函數來得到未知數據$[A]_{0}, K_{M}, V_{max}$
我先假設真實$K_{M}=0.5 mM, V_{max} = 0.1 mM/min, [A]_{0} = 0.5 mM$
經過非線性迴歸後可得$K_{M}=0.48 mM, V_{max} = 0.09 mM/min, [A]_{0} = 0.68 mM$

% Real parameters
KM_real = 0.5; % mM
Vmax_real = 0.1; % mM/min
A0_real = 0.5; % mM
% Corrupt A with some random noise
A = y(1,:)*KM_real; 
A = A + 0.03*randn(1,length(A));
t = x*KM_real/Vmax_real;
% Specify the model function for regression
modelfun = @(b,t) b(1).*lambertw(b(3)./b(1).*exp(b(3)-b(2).*t./b(1)));
% Initial guess of parameters
beta0 = [0.5 0.1 0.5]; 
% Non-linear fitting
mdl = fitnlm(t,A,modelfun,beta0);
% b(1) = KM, b(2) = Vmax, b(3) = A0
KM_est = table2array(mdl.Coefficients(1,1))
Vmax_est = table2array(mdl.Coefficients(2,1))
A0_est = table2array(mdl.Coefficients(3,1))
A_est = modelfun([KM_est Vmax_est A0_est],t);
scatter(t,A,'r','filled','MarkerFaceAlpha',0.4);hold on
plot(t,A_est,'k','LineWidth',2);
ylabel('[A](mM)');
xlabel('t (min)');
legend('Raw data','Estimation')