数理生物学早稲田大学理工学部電気・情報生命工学科高松敦子 2015 年度版 1 1 最新版は http : //www.f.waseda.jp/atsuko ta/M athBio.html 参照のこと目次第 1 章現象をモデル化するということ ( [1, 11]) 3 第 2 章細胞内反応と酵素反応 2.1 細胞内の化学反応と酵素：ミカエリス-メンテンの式 . . . . 2.1.1 化学反応と活性化エネルギー . . . . . . . . . . . . 2.1.2 質量作用の法則による化学反応の記述 . . . . . . . 2.1.3 可逆反応と平衡定数 . . . . . . . . . . . . . . . . . 2.1.4 酵素反応速度論 (Enzyme Kinetics) . . . . . . . . . 2.1.5 準定常状態近似 (quasi-steady state approximation) 2.1.6 平衡近似 (equilibrium approxmation) . . . . . . . . 2.2 酵素反応の阻害と協同的反応 . . . . . . . . . . . . . . . . . 2.2.1 競合阻害 (competitive inhibition) . . . . . . . . . . 2.2.2 非競合阻害 (noncompetitive inhibition) . . . . . . . 2.2.3 基質阻害 (substrate inhibition) . . . . . . . . . . . 2.2.4 不競合阻害 (uncompetitive inhibition) . . . . . . . 2.2.5 協同的反応 . . . . . . . . . . . . . . . . . . . . . . 2.2.6 ヒルの式 (Hill equation) . . . . . . . . . . . . . . . 6 6 7 7 8 10 12 14 15 16 18 20 21 23 26 第 3 章生体反応の制御とリズム 3.1 生体反応の制御 [6, 7, 11–13] . . . . . . . . . . . . . . . . . 3.1.1 1 成分系の例-細胞分化のモデルの例から- [11] . . . 3.1.2 2 成分系の例-ウイルス感染モデルの例から- [8] . . . 3.1.3 2 成分系の解析方法 . . . . . . . . . . . . . . . . . . 3.1.4 安定性解析 . . . . . . . . . . . . . . . . . . . . . . 3.1.5 線形システムの相平面における解の挙動 . . . . . . 3.1.6 生体反応の制御 . . . . . . . . . . . . . . . . . . . . 3.1.7 相平面解析による振動反応の一般論 . . . . . . . . . 3.1.8 Hopf 分岐 . . . . . . . . . . . . . . . . . . . . . . . 3.2 2 成分系における振動反応の実例 . . . . . . . . . . . . . . 3.2.1 Ca2+ 振動 [9] . . . . . . . . . . . . . . . . . . . . . 3.2.2 基質消費型振動 (Substrate-depletion oscillation) [6] 3.2.3 活性-抑制型振動 (Activetor-inhibiotr oscillation) [6] 29 31 34 38 41 41 42 45 46 51 54 54 57 63 1 3.3 3.4 3 成分系の振動反応 . . . . . . . . . . . . . . . . . . . . . . 66 時間遅れフィードバック系の振動反応 . . . . . . . . . . . . 69 第 4 章膜と興奮系 4.1 振動反応から興奮系の反応へ . . . . . . . . . . . . . . . . . 4.2 興奮性膜のモデル-Hodgkin-Huxley モデル- . . . . . . . . . 4.2.1 膜 . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 膜のモデル . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Hodgkin-Huxley モデルの導入 . . . . . . . . . . . . 4.2.4 膜電位固定法 (Voltage-clamp method) . . . . . . . 4.2.5 イオンゲート . . . . . . . . . . . . . . . . . . . . . 4.2.6 Hodgkin-Huxley 方程式 . . . . . . . . . . . . . . . 4.3 4 変数 Hodgkin-Huxley 方程式から 2 変数モデルへ-FitzHughNagumo 方程式- . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Hodgkin-Huxley モデルの定性的解析 . . . . . . . . 4.3.2 FitzHugh-Nagumo 方程式 . . . . . . . . . . . . . . 2 74 74 74 74 76 77 79 81 84 88 88 90 第1章現象をモデル化するということ ( [1, 11]) 生命現象のみならず一般に何かしらの「現象」を数理モデル化するには、まず注目した変数の時間発展を考えなければならない。ここでは、ポピュレーションダイナミクス (個体数変動論) を例に取って考えてみよう。例えばいわゆる「ねずみ算」、つまり、ネズミの個体数が増えていく過程を考える。１匹のネズミが 4 匹ずつ子供を産むとする。ただし、ここでは簡単のためにネズミの雄雌の区別はなく単体で子供を産めるものとする。また、観測時間内にネズミが死ぬことはなく、継続的に子供も生めるとする。第 0 世代に 1 匹のネズミがいたとすると、第 1 世代は 1 + 4 · 1 = 5 匹、第 2 世代は 5 + 4 · 5 匹となり、第 k 世代の個体数 nk は、 nk = nk−1 + 4nk−1 となることがわかる。ここで、話を一般に拡張するために、1 匹のネズミが 1 回に子供を生む数を成長率 α とおくと、 nk − nk−1 = αnk−1 とかける。更に時間も取り入れて、1 世代の時間を ∆t とすると、単位時間当たりの個体数の増殖速度は ∆n = nk − nk−1 = αnk−1 ∆t とかける。∆t → 0 の極限を取ると、以下の微分方程式が導出できる。 dn = αn dt (1.1) この式が個体数変動のダイナミクスを表す式となる。ここで、当然であるが個体数は正の値をとるのでここでは、n ≥ 0 としておく。この式は t = 0 のときの個体数を n(0) = n0 とした初期条件の下で解くことができて、 n = n0 eαt 3 (1.2) となり、「ねずみ算」とは、指数関数的に個体数が増えていく過程であることがわかる。この式に従った成長の仕方をマルサス 1 の法則という。しかし、ネズミとはいえ無限に個体数が増えていくことはなくやがてある値に収束していくことは予想できるだろう。つまり、今対象としている生物には生息場所や食料の量などに限界があり個体数の成長にも限界があるということである。このことを式 (1.1) に取り入れてみよう。つまり、個体数が増えるほど餌や生息場の取り合いが生じて成長率が個体数に比例して減少するという形にしてみる。 dn = (a − bn)n dt (1.3) ここで、α = a − bn とおいており、a は個体数が非常に少ない場合の成長率を、b は限界に達する率を表している。この式をロジスティック方程式 (logistic equation) という。式 (1.1) の場合と同様に初期条件 n(0) = n0 の下に解くことができて、a − bn0 > 0 の場合には、 n(t) = nmax −at 1 + ( nmax n0 − 1)e (1.4) と求められる。但し、nmax = a/b とおいている。式 (1.4) で t → ∞ とすると n = nmax となることが直ぐわかる。これは、初期値 n0 に関係なくある一定値 nmax に収束することを意味する。同じことが、微分方程式 (1.3) からも導出可能である。時間がどれだけ経っても状態 n がそれ以上変化しない条件を求めればよい。すなわち、 dn =0 dt となる条件を求めればよく、 (a − bn)n = 0 より、n = 0 または n = nmax が求まる。これら２つの解を平衡解または固定点、あるいは定常解という。各自、初期値 n0 をいろいろと変えて式 (1.4) のグラフを描いてみよう。どの初期値から出発しても無限に増殖することはなく最終的に n = nmax に収束し、モデルの狙いであった成長限界の存在が確認できるだろう。成長限界の値である nmax のことを環境収容力というが、環境収容力を超えた初期値から出発しても減衰してやはり環境収容力値へと収束する。では、初期値 n0 = 0 の場合はどうなるだろうか。式 (1.4) から n = 0 となり、これは、先ほど求めた平衡解のうちの１つと一致することが確認できる。つまり、元の個体数が 0 であれば当然 1 Malthus: イギリスの経済学者で人口の増え方について数学的に考察した。 4 個体数は増加せず、ずっと 0 のままである。 (練習問題 1-1) 各自、式 (1.1), (1.3) からそれぞれ解 (1.2), (1.4) を導出せよ。このときラプラス変換は用いないこと。次に、様々な初期条件 n(0) = a a n0 (1) 0 < n0 < 2b , (2) 2b < n0 < ab , (3) n0 > ab のとき、t − n プロットはどのようになるか、描いて確認せよ。環境収容力の存在も確かめておこう。さらに、n が小さい場合式 (1.3) のグラフは式 (1.1) のグラフと比較してどうか。 5 第2章 2.1 細胞内反応と酵素反応細胞内の化学反応と酵素：ミカエリス-メンテンの式エネルギー (AB)* Ea A+B C 分子の状態図 2.1: 化学反応における各分子のエネルギー状態細胞内では様々な生体分子が相互作用し化学反応が起こっている。通常の化学反応と大きく異なる特徴は酵素の存在である。酵素が果たす役割は体温程度の温度で化学反応を効率的に触媒することにある。このセクションでは細胞内化学反応を記述する上で最も基礎となるミカエリス-メンテンの (Michaelis-Menten equation) の式を導出する 1 。 1 [2] や [3] を参照。前者は生化学の入門書、第 4 章を参照、酵素反応の基礎的取り扱いを押さえたい人はこちら。後者は応用数学者が書いたもの、chapter 1 を参照。数学的取り扱いについて学びたい人はこちら。 6 2.1.1 化学反応と活性化エネルギーある化学反応が起こるかどうかは、反応分子同士が衝突すること、かつ、その分子が反応に十分なエネルギーを持っているかに依存する。このことを次の化学反応式で表せる反応を例に考えてみよう。 k A + B⇀C (2.1) 各分子のエネルギー状態を図示したものが、図 2.1 であるが、A + B から C への反応は直接進むのではなく、中間状態 AB∗ を経る必要がある。これらの状態エネルギー差が活性化エネルギー (activation energy)Ea であり、反応分子はこれ以上のエネルギーを持っている必要がある。その確率は、分子のエネルギー分布がボルツマン分布に従うとすれば 2 、 k = A exp(−Ea /RT ) (2.2) と表せる。ここで k は反応速度定数 (rate constant) と言い、この式をアレニウスの式 (arrhenius equation) という。T は絶対温度、A は頻度因子 (frequency factor) である。3 酵素はこの Ea を低くする作用がある。式 (2.2) の両辺の対数を取ると、 ln k = ln A − Ea 1 · R T (2.3) と表される。反応速度の温度依存性を実験により計測し、この式に従ってプロットを作成すると (アレニウスプロットと言う)、A および Ea を見積もることができる。 2.1.2 質量作用の法則による化学反応の記述細胞内のダイナミクス（動的挙動）を数式で記述するための準備段階として、化学反応を微分方程式で記述する方法について解説する。化学反応式 (2.1) において、反応物 A, B の濃度を [A], [B]、生成物 C の濃度を [C] と表すとき、生成物の生成速度は、つまり、反応物が反応できる確率は、先ほど示した反応速度定数 k と反応物の濃度 [A], [B] に依存するので 2 化学反応の専門書を参考のこと。たとえば [4]。化学としての反応速度論を扱っている。酵素反応についても言及。第 1 章∼第 3 章を参照。 √ 3 同じく [4] などを参照。頻度因子は分子の大きさ、形状、質量に依存し、 T に比例する。生体内で化学反応を考える場合には、温度変化があったとしても非常に小さい範囲で、式 (2.2) で T について 1/2 乗で変化する項は、指数的に変化する項に対してほどんど一定と見なせるので、A は一定として取り扱う場合が多い。 7 (質量作用の法則 (low of mass action) より)、次の微分方程式の形に書ける 4。 d[C] = k[A][B] (2.4) dt (練習問題 2-1) 次の化学反応式について、生成物 P についての生成速度の式を立てよ。 (1) k A + B + C⇀P (2) k A + B⇀P + C (3) k A⇀P + C (4) k 2A ⇀ P 2.1.3 可逆反応と平衡定数次のような順反応と逆反応からなる反応を可逆反応 (reversible reaction) という。 k+ A + B!C k− (2.5) このとき A の濃度は順反応で減少し、逆反応で増加するので、その反応速度方程式は d[A] = (自分で考えて下さい) (2.6) dt とかける。この反応は、時間が十分経過すれば、順反応速度と逆反応速度が釣り合うようになり、やがて平衡状態 (equilibrium state) に達する。このとき、見かけ上、反応物質 A の濃度は時間変化しなくなるので 5 、 d[A] =0 dt 4 この式は、あくまで、巨視的に見た場合の近似式なので念のため注意。適切な濃度で十分に攪拌された系を想定している。 5 平衡状態と定常状態 (steady state) は、見かけ上の変化が見られなくなる点で非常によく似ているが、厳密には異なるので注意が必要である。定常状態は、注目した物理量が動的であるが時間的に不変な状態である。例えば、後で学ぶ酵素反応で反応中間体の反応速度に着目した場合、反応中間体の濃度に変化がなくても、基質の消費も生成物の生成も続いているので、反応系全体で見ればそれは平衡状態ではなく、定常状態である。一方、反応系全体で各成分の濃度に変化が無くなれば、平衡状態となる。 8 とおける。式 (2.6) にこれを代入し、平衡状態での A,B,C の濃度をそれぞれ [A]eq , [B]eq , [C]eq とおくと、 K≡ [C]eq k+ = k− [A]eq [B]eq (2.7) となる。K を平衡定数といい、アレニウスの式 (2.2) からわかるように、温度の関数である。注意したいことは、酵素が存在していても平衡定数 K の値が変わることはなく、酵素は単に活性化エネルギーを小さくするという点である。なぜならば、例えば酵素が順反応の反応速度を大きくするように働く場合には、同時に逆反応の反応速度も速くなるからだ。 (a) 酵素がない場合 + (b) 酵素がある場合 + - 図 2.2: 化学反応における各分子のエネルギー状態. 酵素の存在は活性化エネルギーの大きさに影響するだけで分子のエネルギー状態には影響しない. (練習問題 2-2) 次の化学反応式について、反応物 A および生成物 P についての生成速度の式を立てよ。 (1) k+ A + B + C!P k− (2) k+ A + B!P + C k− (3) k+ A!P + C k− (4) k+ 2A ! P k− 9 (練習問題 2-3) 酵素の存在で活性化エネルギーが小さくなったときでも、平衡定数の値が変わらないことを式 (2.2) および式（2.7）を用いて証明せよ。図 2.2 を参考にすると良い。 2.1.4 酵素反応速度論 (Enzyme Kinetics) 基質 (Substrate)S が酵素 (Enzyme)E を触媒 (catalyst) として生成物 (Product)P となる反応を考えよう。酵素反応には、図 2.3 に示したように、基質濃度 [S] がどんなに高くても (酵素濃度 [E] は一定という条件で)、反応速度 V(生成物 P の生成速度) は頭打ちになるという特徴がある。つまり、基質濃度に対して常に線形で増加するわけではない。酵素反応の特徴を纏めると次のようになる。 V V∼const. [S] [S] 図 2.3: 酵素反応速度の基質濃度依存性. 全酵素濃度は一定であるとする. 反応速度は基質濃度が低いとき基質濃度にほぼ比例するが、基質濃度が増加するに従い、反応速度は飽和し基質濃度がどんなに増加しても一定となる. (1) 基質濃度 [S] が低濃度のとき、一定の酵素濃度では [S] に比例する。 (2) 基質濃度 [S] が高濃度のとき、一定になる。従って、単純に質量作用の法則を適用することはできない。その理由を考えてみよう。いま、酵素反応について、反応物を S, E、生成物を P と 10 して、 k S + E⇀P + E とのように反応すると考えると、 V = d[P] = k[S][E] dt のように書けると思うかもしれない。しかし、これでは基質濃度 [S] に対して常に [E] に比例して増加してしまうことになってしまう。そこで Michaelis と Menten(1913) は、反応中間体 (複合体 (Complex))C の存在を仮定し次の反応式を提案した。 k+1 k S + E ! C ⇀2 P + E k−1 (2.8) [S] が小さいとき、式 (2.8) の最初の反応に平衡状態が近似的に成立すると仮定した 6 。さらに、[S] 大きいとき遊離 E はすぐさま S と結合するので遊離 E は不足して反応は飽和する。これらの仮定は酵素反応に見られる現象の特徴をうまく説明している。彼らの功績から酵素反応速度に関する関係式をミカエリス-メンテンの式 (Michaelis-Menten equation)(2.18) という。彼らによる導出方法は平衡近似として §2.1.6 で説明する。この式の導出法に関して後に Briggs と Haldane (1925) によって修正が加えられ一般化された。ここでは彼らの方法に従って解説する。彼らは、 [S] ≫ [E] を仮定し、複合体 C が分解して酵素 E が解放されるとすぐさま別の基質 S と結合することに着目した。また、反応 (2.8) の後半では生成物 P はすぐさま除去され 7 、その逆反応は無視できるものとしている。これらのことから、k2 > k−1 の場合を考えると、S と E から C が生成する速度と P へと崩壊する速度が均しい状態、準定常状態 (quasi-steady state) 8 を仮定した。§2.1.5 では、この仮定を用いてミカエリス-メンテン式を導出する。その前に、化学種の生成速度を整理しておこう。簡単のためにそれぞれの化学種の濃度を s = [S], e = [E], c = [C], p = [P] とおくことにする。質量作用の法則から、それぞれの化学物質の生成速度に関する方程式を次のように書くことができる。 6 厳密には k2 ≪ k−1 のときに成立する。実際に細胞内では生成物は直ぐに別の反応に用いられることになるので、この仮定は妥当である。 8 準定常状態では、ある物理量が、長い時間スケールで見たら変化しているが、短い時間スケールで見たらほとんど変化していないので、変化がないものと近似できるような場合である。どのような時間スケールでも変化しない定常状態とは異なる。 7 11 ds dt de dt dc dt dp dt = k−1 c − k+1 se (2.9) = k−1 c + k2 c − k+1 se (2.10) = −k−1 c − k2 c + k+1 se (2.11) = k2 c (2.12) いま、ここで行いたいのは、ミカエリス–メンテンの式の導出である。つまり、生成物の反応速度 V と基質の濃度 s の関係式の導出を行う。V = dp dt なので、式 (2.12) を最終的に用いることがわかる。従って、残りの式から c を求めればよい。ところで、式 (2.10) と式 (2.11) の右辺は符号が異なるだけなので、次の関係式が成立する。 de dc + =0 dt dt e + c = const ≡ e0 (2.13) (2.14) ここで、e0 は遊離酵素、結合酵素を含めた全酵素の濃度である。この式から変数 e は消去されて、本質的には c と s だけの 2 変数だけ考えれば良いことがわかる。 2.1.5 準定常状態近似 (quasi-steady state approximation) 前のセクションで触れたように Briggs と Haldane は C の生成速度が崩壊速度と均しいとし、C の準定常状態を仮定した。つまり式 (2.11) で、 dc =0 dt (2.15) とおける。この条件と、式 (2.14) より、e を消去すると、 c= e0 s s + Km (2.16) が求まる。ここで、 Km ≡ k−1 + k2 k+1 (2.17) とおいている。Km をミカエリス定数 (Michaelis constant) という。生成物の反応速度 V は式 (2.12) より次式で与えられるので、式 (2.16) より、 dp = k2 c dt Vmax s = s + Km V = 12 (2.18) が導出できる。ここで、 Vmax = k2 e0 とおいた。式 (2.18) をミカエリス-メンテンの式 (Michaelis-Menten equation) という。これを図示したものが、図 2.4(a) であり、節 2.1.4 の冒頭で挙げた酵素反応の特徴をよく表している。反応速度がちょうど反応速度の最大値 Vmax の半分となったときの基質濃度 s が Km と等しくなる。(各自確かめておくこと。) つまり、Km は反応がその最大速度に達するのに必要な基質濃度の目安となる。この値が小さいほど、基質と親和性の高い酵素ということになる。これは次の Lineweaver-Burk プロットという形式に変形でき、実験ではこのプロット (図 2.4(b)) から、Vmax , Km を求める 9。 1/V = s + Km Vmax s = 1/Vmax + (2.19) Km 1/s Vmax (2.20) (練習問題 2-4) ミカエリス-メンテンの式 (2.18) で、本節の最初に整理したミカエリス-メンテン則の性質が表現できているか以下の方法で各自確認してみよう。各条件を s = const., s が非常に小さい場合, s が非常に大きい場合として、式 (2.18) がどう近似されるか示して議論しなさい。 (練習問題 2-5) ミカエリス-メンテンの式 (2.18) を線形にプロットする手法として、Lineweaver-Burk プロット (1/s–1/V ) の他に、Eadie-Hofstee プロット (V /s–V ) および Hanes-woolf プロット (s–s/V ) と呼ばれるものがある。これらのプロットを作成するために式 (2.20) を参考に、それぞれ対応する関係式を導出なさい。図 2.4(b) に習ってプロットを作成せよ。さらに、3 つの手法の利点、欠点について議論せよ。 9 コンピュータによる解析が発達する以前は、このような手法でデータを線形の関係に変換してから線形回帰によって各係数を求めていた。しかし線形の関係に変換する際に、値によって誤差が異なり係数の見積もりを誤ることがある。現在では加工前のデータから非線形回帰によって各係数を求める。 13 (a) (b) V 1/V Vmax Vmax /2 1/Vmax s Km slope = Km /Vmax 1/s -1/K m 図 2.4: ミカエリス-メンテンの式 (2.18) による s−V プロット (a), 1/s−1/V プロット (b). 2.1.6 平衡近似 (equilibrium approxmation) 前セクションでは複合体の準定常状態を仮定して、ミカエリス-メンテン式を導出したが、遊離酵素 E と基質 S から複合体 C 形成する反応 (反応式 (2.8) の１番目の反応) が平衡状態にあると仮定して導出する方法もあり、この方法を平衡近似という。歴史的にはこちらの方が先であり、複雑な系を概観するには適しているので紹介する。式 (2.7) より平衡定数は、その反応の反応物、生成物の濃度比で決まることを思い出しておこう。このとき、Ks = k−1 /k+1 とおけば 10 、 se − Ks c = 0 (2.21) となることがわかり、これより c が求まる。従って、式 (2.14) と併せて、 dp = k2 c dt Vmax s = s + Ks V = (2.22) が求まる。式 (2.18) と式の形は同じであるが、Ks は Km とは定義が異なることに注意しよう。 10 解離定数という。節 2.1.3 で定義した平衡定数の逆数となっていることに注意。 14 関連図書 [1] 関村利郎, 竹内康博, 梯正之, 山村則男, 理論生物学入門, 現代図書, 2007 [2] 猪飼篤, 生化学, 化学入門コース 8, 岩波書店, 1996 [3] J. Keener and J. Sneyd, Mathematical Physiology, Springer, 1991 [4] 斉藤勝裕, 反応速度論：化学を新しく理解するためのエッセンス, 三共出版, 1997. [5] A. Bruce et al., 中村佳子, 松原謙一監訳, Essential 細胞生物学原書第 2 版, 江南堂, 2006. [6] C. P. Fall, E. S. Marland, J. M. Wagner, J.J. Tyson, Computational Cell Biology, Springer, 2000. [7] J. D. Murray, Mathematical Biology I, 3rd Ed., Springer, 2001 [8] 稲葉寿編著, 感染症の数理モデル, 培風館, 2008. [9] A. Goldbeter, G. Dupont, and M.J.Berridge (1990) Mimimal model for signal-induced Ca2+ Oscillations and for their frequency encoding through protein phosphorylation. Proc. Natl. Acad. Sci. USA, 87, 1461-1465. [10] N. MacDonald (1989), Biological delay systems: linear stability theory, Cambridge University Press. [11] L. Edelsterin-Keshet, Mathematical Models in Biology, SIAM, 1988 [12] G. D. Vries, T. Hillen, M. Lewis, J. Muller, B. Schonfish, A Course in Mathematical Biology, SIAM, 2006 [13] D. Kaplan and L. Glass, Understanding Nonlinear Dynamics, Springer, 1995 94 [14] A. L. Hodgkin and A. F. Huxley (1952) Quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117, 500-544. 95