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

Rはもっぱら統計解析にのみ使っていたのだけど，[[@T_Hash|http://twitter.com/T_Hash/statuses/687249842]]にRで数式解けるの? みたいな==挑戦状==質問をされたのでやってみた．Runge-Kutta-Gill 法で[[Michaelis-Menten kinetics|http://en.wikipedia.org/wiki/Michaelis-Menten_kinetics]]の連立微分方程式を解いてみる．

参考: [[R で微分方程式：odesolve|http://clinical-pk-pd.cocolog-nifty.com/blog/2007/11/r_odesolve_fe42.html]]

ode solverはRのライブラリ[[odesolve|http://cran.r-project.org/src/contrib/Descriptions/odesolve.html]]に実装されているのでまずはインストールする．

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である。この関数は初期値、時間、微分方程式のリスト、そしてパラメータが格納されたベクトルを引数にとる。結果としてそれぞれの時間での解が返される。

{{image 0, '画像の説明', nil, [400,400]}}