廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

メモ:逐次反応の微分方程式

逐次反応の微分方程式

反応物 A が反応物 B を生成し, 反応物 B が C を生成する状況を考える.

\displaystyle A\stackrel{k_1}{\longrightarrow }B\stackrel{k_2}{\longrightarrow }C

  • A は A の量に比例する速度で減少する.
  • B は B の量に比例する速度で増加し, C の量に比例する速度で減少する.
  • B は B の量に比例する速度で増加する.

とすると微分方程式系,

\displaystyle \frac{dA}{dt} = -k_1 A

\displaystyle \frac{dB}{dt} =  k_1 A - k_2 B

\displaystyle \frac{dC}{dt} = k_2 B

が導かれる.

任意の時点で

\displaystyle A+B+C = A_0 + B_0 +C_0

であることに注意する.

ここで初期濃度を A(0) = A_0, B(0) = B_0, C(0) = C_0 と表記した.

上記の微分方程式系は解析的に解くことができ,

\displaystyle A(t) = A_0 \exp(-k_1 t)

\displaystyle B(t) = B_0 \exp(-k_2 t)+\frac{A0}{1-k_2/k_1} ( \exp(-k_2 t)-\exp(-k_1 t) ) )

\displaystyle C(t) = C_0 + A_0 (1-\exp(-k_1 t)) +\\ ~~~~B_0 \left[ (1-\exp(-k_2 t)) - \frac{A_0/B_0}{1-k_2/k_1}(\exp(-k_2 t)- \exp(-k_1 t)) \right]

である.

R による計算

f:id:abrahamcow:20150301023521p:plain

deSolve パッケージによって求めた数値解を点で, 解析解を曲線で表している.

両者が一致していることがわかる.

コードは以下の通り.

library(deSolve)

rxnrate=function(t,c,parms){
  
  # rate constant passed through a list called parms
  k1=parms$k1
  k2=parms$k2
  
  # c is the concentration of species
  
  # derivatives dc/dt are computed below
  r=rep(0,length(c))
  r[1]=-k1*c["A"] #dcA/dt
  r[2]=k1*c["A"]-k2*c["B"] #dcB/dt
  r[3]=k2*c["B"] #dcC/dt
  
  # the computed derivatives are returned as a list
  # order of derivatives needs to be the same as the order of species in c
  return(list(r))
  
}
cinit=c(A=10,B=0.1,C=0.1)
parms=list(k1=3,k2=1)

times <- seq(0, 20, by=0.1)

out=ode(y=cinit,times=times,func=rxnrate,parms=parms)

matplot(out[,1],out[,-1],pch=16,lty=1,lwd=2,xlab="time",ylab="concentration",
        col=c("royalblue","orange2","forestgreen"))

Afn <- function(t,A0,k1){
  A0*exp(-k1*t)
}

Bfn <- function(t,A0,B0,k1,k2){
  B0*exp(-k2*t)+((A0/(1-k2/k1)) *( exp(-k2*t)-exp(-k1*t) ) )
}

Cfn <- function(t,A0,B0,C0,k1,k2){
  C0 + A0*(1-exp(-k1*t)) + 
    B0*( (1-exp(-k2*t)) - ((A0/B0)/(1-k2/k1))*(exp(-k2*t)- exp(-k1*t)) )
}


curve(Afn(t=x,A0=cinit["A"],k1=parms$k1),add=TRUE,col="royalblue",lwd=3)
curve(Bfn(t=x,A0=cinit["A"],B0=cinit["B"],k1=parms$k1, k2=parms$k2),add=TRUE,col="orange2",lwd=3)
curve(Cfn(t=x,A0=cinit["A"],B0=cinit["B"],C0=cinit["C"],k1=parms$k1, k2=parms$k2),add=TRUE,col="forestgreen",lwd=3)

legend(locator(1),legend=c("A","B","C"),lty=1,lwd=3,pch=16,
       col=c("royalblue","orange2","forestgreen"),bty = "n")

追記

Rate equation - Wikipedia, the free encyclopedia にふつうに解が書いてあった. これなら  B_0=0 でも大丈夫.

\displaystyle A(t)=A_0 e^{-k_1 t}

 \displaystyle B(t) =\left\{ \begin{array}{l}   A _{0}\frac{k_{1}}{k_{2}-k_{1}}\left( e^{-k_{1}t}-e^{-k_{2}t} \right)+B_{0}e^{-k_{2}t}& k_{1}\ne k_{2}  \\    A_{0}k_{1}te^{-k_{1}t}+B_{0}e^{-k_{1}t} & \text{otherwise}  \\ \end{array} \right.

\displaystyle C(t)  =\left\{ \begin{array}{l}   A  _{0}\left( 1+\frac{k_{1}e^{-k_{2}t}-k_{2}e^{-k_{1}t}}{k_{2}-k_{1}} \right)+ B  _{0}\left( 1-e^{-k_{2}t} \right)+C  _{0} & k_{1}\ne k_{2}  \\    A  _{0}\left( 1-e^{-k_{1}t}-k_{1}te^{-k_{1}t} \right)+  B  _{0}\left( 1-e^{-k_{1}t} \right)+ C  _{0} & \text{otherwise}  \\ \end{array} \right.

数値解との残差平方和は, ほぼ0になってる.

library(deSolve)

rxnrate=function(t,c,parms){
  
  # rate constant passed through a list called parms
  k1=parms$k1
  k2=parms$k2
  
  # c is the concentration of species
  
  # derivatives dc/dt are computed below
  r=rep(0,length(c))
  r[1]=-k1*c["A"] #dcA/dt
  r[2]=k1*c["A"]-k2*c["B"] #dcB/dt
  r[3]=k2*c["B"] #dcC/dt
  
  # the computed derivatives are returned as a list
  # order of derivatives needs to be the same as the order of species in c
  return(list(r))
  
}
cinit=c(A=1,B=0,C=0)
parms=list(k1=3,k2=1)

times <- seq(0, 20, by=0.1)

out=ode(y=cinit,times=times,func=rxnrate,parms=parms)

Afn <- function(t,A0,k1){
  A0*exp(-k1*t)
}

Bfn <- function(t,A0,B0,k1,k2){
   A0 * (k1/(k2 - k1)) * ( exp(-k1*t)-exp(-k2*t) ) + B0 *  exp(-k2*t) 
}

Cfn <- function(t,A0,B0,C0,k1,k2){
  A0*( (1 + ((k1*exp(-k2*t) - k2* exp(-k1*t) )/(k2-k1))) ) + B0*(1 - exp(-k2*t)) + C0
}
> sum((out[,2]-Afn(t=times,A0=cinit["A"],k1=parms$k1))^2)
[1] 4.344727e-12
> sum((out[,3]-Bfn(t=times,A0=cinit["A"],B0=cinit["B"],k1=parms$k1, k2=parms$k2))^2)
[1] 9.169868e-12
> sum((out[,4]-Cfn(t=times,A0=cinit["A"],B0=cinit["B"],C0=cinit["C"],k1=parms$k1, k2=parms$k2))^2)
[1] 9.510325e-13