_ Rで常微分方程式を解いてみた．Michaelis-Menten kineticsを例として

Rはもっぱら統計解析にのみ使っていたのだけど，@T_HashにRで数式解けるの? みたいな挑戦状質問をされたのでやってみた．Runge-Kutta-Gill 法でMichaelis-Menten kineticsの連立微分方程式を解いてみる．

参考: R で微分方程式：odesolve

ode solverはRのライブラリodesolveに実装されているのでまずはインストールする．

install.packages("odesolve")

コードは以下の通り．

library(odesolve)

dydt <- function(t, y, p){
  k10 <- p['k10']
  k01 <- p['k01']
  k20 <- p['k20']
  E0  <- p['E0']
  
  y1 <- -k10*E0*y[1]+(k10*y[1]+k01)    *y[2]
  y2 <-  k10*E0*y[1]-(k10*y[1]+k01+k20)*y[2] 
  y3 <- k20*y[2]

  list(c(y1,y2,y3))
}

params <- c(k10 = 1e3,
                  k01 = 1,
                  k20 = 0.05,
                  E0  = 0.5e-3)
times <- c(0, 0.1*(1:250))
y <- lsoda(c(1e-3,0,0), times, dydt, params)

S  <- y[,2]
ES <- y[,3] 
E  <- params['E0'] - ES
P  <- y[,4]

matplot(times, matrix(c(S, ES, E, P), ncol=4),
        type="l",
        ylab = "concentration (M)",
        xlab = "time (s)"
)
legend(20,1e-3, c("[S]", "[ES]", "[E]", "[P]"), pch="-", col=c(1:4))

すると以下ような図ができる．簡単に解説すると、lsoda関数が ode solverである。この関数は初期値、時間、微分方程式のリスト、そしてパラメータが格納されたベクトルを引数にとる。結果としてそれぞれの時間での解が返される。