医学と数理
医学と数理への招待
はじめに
「医学部での勉強は暗記ばかりで、つまらなさそう。」そのような印象をお持ちではないでしょうか?確かに、覚えるべき事項は多いですし、試験前は特に教科書と睨めっこをしながら知識を詰め込むことが多いかもしれません。しかし、我々の体はある日突然、神の気まぐれによって作り出された訳ではなく、そこには必ず何らかの因果関係が存在している筈でしょう。例えば、我々の体を構成している細胞は、その内側に膨大な種類と数の分子が相互作用をしながら、ある種の「生命らしさ」を保ち続けています。このような仕組みは、気が遠くなるほどの長い年月をかけて、生命が進化の道筋を辿り、その結果として現在の我々の巧妙に設計された「システム」としての生命の神秘とも言えるでしょう。
我々の体はこの世界に実体として存在している以上、その構成要素は物理化学的で普遍的な法則に支配されていると考えるのが自然です。そのような動機の下、人体をはじめとする、生物において見られる魅力的な現象を様々なスケールで解明しようとする研究が盛んに行われています。そこで、この記事ではそのような研究のいくつかを、前提知識の解説を交えながらご紹介します。
関連する研究分野
生命科学の研究手法は様々ですが、わかりやすく大別する概念としてウェットとドライという用語があります。前者は、実験室にて、細胞培養、遺伝子操作、PCRやウェスタンブロッティング、イメージング、行動実験など、実際の生物を使って行う研究手法を指します。一方、後者は主にコンピュータを利用して実験によって得られたデータを解析したり、シミュレーション実験を行ったりします。
一口に「生命現象を数学、物理学、化学や工学などを用いて研究する」と言っても、研究対象や手法によって幾つかの学問分野に大まかに分けることができます。下に列挙しているのは、その例です(ただし、厳密に定義が定まっていない分野も多く、示した記述は筆者なりの解釈です)。
その中でも今回は、「システム医学」に関する話題を中心に取り上げます。
常微分方程式の基礎
生命現象を定量的に扱う際に用いられる数理工学の手法は様々ですが、最も代表的なものは常微分方程式です。ゆえに、この節では簡単に常微分方程式の基礎事項に関する解説を行います。
生命科学に限らず、一般に、自然科学において、ある要素に対する作用とその時間変化を定量的に解析する手法として、微分方程式がよく用いられます。特に未知関数が1つの独立変数(この場合は時間と考えてもらって構いません)を持つ場合は、常微分方程式と呼ばれます。例えば高校の物理で学ぶ運動方程式などがその例です。力$${F}$$を定数とすれば、運動方程式はこのように書けます。
$$
m\frac{d^2 x(t)}{dt^2} =F
$$
この微分方程式は2回積分操作を施せば、
$$
x(t)=\frac{F}{2m}t^2+C_1 t+C_2
$$
と解くことができます($${C_1}$$ 、$${C_2}$$は定数)。このように、自然現象を微分方程式という形で適切にモデル化し、それを解析的に解くことができれば、ある粒子や物体の挙動を予測したり化学反応を定量的に表現したりすることができます。
以下に、よく登場する常微分方程式の例とその解を列挙しておきます[1]。
このように微分方程式を解析的に解くことは自然現象を解析・予測する上で非常に強力な武器のように思えますが、現実問題として全ての微分方程式が解けるわけではありません。そこで登場するのが、数値計算(いわゆるコンピュータ・シミュレーション)です。
シミュレーションの基礎
近年の機械学習ブームもあり、PythonやR、Juliaなどのプログラミング言語は比較的多くの人に馴染みのあるものとなりました。シミュレーションの実装も、これらの言語を用いて行われます。常微分方程式の数値解法は数多く存在しますが、一番代表的なものは以下に示すオイラー法というものです。
例えば、以下のような常微分方程式と初期条件が与えられたとします。
$$
\frac{dx(t)}{dt}=f(t,x)\\
x(t_0 )=x_0
$$
時間の刻み幅を $${h}$$とします。このときT時間後の数値解は、次の公式で定義されます。
$$
x_(n+1)=x_n+hf(t_n,x_n)\\
t_n≤T
$$
このように、微分方程式の右辺の1次近似を利用して、離散的に計算することで、微分方程式を計算する方法をオイラー法と呼びます。
例えば
$$
\frac{dx(t)}{dt}=-3x\\
x_0=2
$$
という微分方程式のオイラー法による数値解は
$$
x_{n+1}=x_n-3hx_n\\
x_0=2
$$
という漸化式を解くことにより
$$
x(t)=2(1-3h)^{\frac{t}{h-1}}
$$
と得られます。
実際に、この微分方程式の解析解
$$
x(t)=2e^{-3t}
$$
と比較すると、$${h→0}$$の極限で数値解が解析解に一致することがわかりますね。
このほかにも、オイラー法をより発展させた「ルンゲ=クッタ法」などの数値解法があります。また、多くのプログラミング言語ではこれらの数値解法を実装するためのライブラリが提供されていますので、初心者でも容易に実践することが可能です。
薬理学への応用
薬力学とのかかわり
医学における重要な分野の一つに薬理学という分野があります。この学問分野では、薬物が体内でどのような生理的・生化学的作用を及ぼすのか、またそれらは体内でどのように吸収・分布・代謝・排泄されるのか、ということが研究されています。
一般に、薬物を含む多くの化学物質は、細胞膜表面に存在する受容体という物質によって捉えられ、その後、リン酸化などの化学修飾を通して、細胞内に存在する下流のシグナル伝達分子群によって、情報が伝達されます。その信号は最終的には転写因子による遺伝子調節などの細胞内における様々な変化を引き起こします。
このような、薬物(化学物質)と受容体の相互作用は、先ほどのような微分方程式を用いてモデル化することが可能です。
以降、便宜的に薬物(Drug)を$${D}$$、受容体(Receptor)を$${R}$$と表します。薬物$${D}$$と受容体$${R}$$が結合して、複合体$${DR}$$を形成する可逆的な化学反応が次のように書けるとしましょう。
$$
D+R⇄DR
$$
この反応が質量作用の法則に従うと仮定し、結合速度定数および解離速度定数をそれぞれ$${k_{+1}}$$、$${k_{-1}}$$とすると
$$
\frac{d[DR]}{dt}=k_{+1} [D][R]-k_{-1} [DR]
$$
が成り立ちます。この反応が平衡状態に達すると、この値は0となるので
$$
\frac{[D][R]}{[DR]}=\frac{k_{-1}}{k_{+1}} =K_d
$$
が成り立ちます。このとき$${K_d}$$は解離定数と呼ばれます。また受容体の総量を$${[R]_{tot}=[R]+[DR]}$$とおくと
$$
[DR]=\frac{[R]_{tot}}{(1+\frac{K_d}{[D]})}
$$
が成り立ちます。ここで、薬理作用の強さ$${E}$$が$${[DR]}$$に正比例すると仮定すれば、この薬物の濃度反応関係は
$$
E=\frac{E_{max}}{(1+\frac{EC_{50}}{[D]})}
$$
と表すことができます。ここで、$${EC_{50}}$$は最大薬理作用$${E_{max}}$$の50%が得られる薬物濃度を表します。
同様に個体レベルでは薬物の投与量を用量と規定し、その用量と薬理作用の関係を用量反応関係といい、そのグラフを用量反応曲線と呼びます。用量反応曲線は、用量を対数表示することでS字型を示すことが多いです。
以上の例のように、受容体に結合して機能を促進するような薬物をアゴニストと呼びます。一方で、受容体に結合することでその機能を抑制する薬物も存在し、それらはアンタゴニストと呼ばれます。
そこで、次にアンタゴニストが与えるアゴニストと受容体の結合に対する影響を簡単な数理モデルも用いて評価してみましょう。
アンタゴニストには、受容体がアゴニストに結合する部位に直接結合することで、競合的にアゴニストの作用を阻害する、競合的アンタゴニストと呼ばれるものと、受容体のアロステリック部位と呼ばれる部分に結合することで受容体の構造変化を引き起こし、アゴニストとの結合を非可逆的に阻害する、非競合アンタゴニストと呼ばれるものが存在します。
まず、競合アンタゴニストについて考えてみましょう。競合アンタゴニストを$${I}$$、受容体と競合アンタゴニストの複合体を$${IR}$$とし、競合アンタゴニストと受容体の結合が次の式で表されるとします。
$$
I+R⇄IR
$$
この反応が質量作用の法則に従うと仮定し、競合アンタゴニストの解離定数を$${K_i}$$とします。
また、受容体の総量を$${[R]_{tot}=[R]+[DR]+[IR]}$$とおくと
$$
[DR]=\frac{[R]_{tot}}{1+\frac{K_d (1+\frac{[I]}{K_i })}{[D]}}
$$
が導かれます。つまり、競合アンタゴニストがない場合と同様の薬効を得るためには$${[D]}$$を$${1+\frac{[I]}{K_i}}$$ 倍にする必要があるということになりますね。
一方、非競合アンタゴニストの存在下では、一般的にアゴニストの$${EC_{50}}$$は変化させずに、$${E_{max}}$$を低下させることが多いです。
それぞれの条件下での用量反応曲線は図1のグラフのようになります。
薬物動態学とのかかわり
投与した薬物の血漿濃度の時間変化は、一般に吸収、分布、代謝、排泄のバランスによって規定されます。生体内における薬物の動態を定量的に評価することができれば、患者さんに薬物を投与する際の至適投与量や至適投与時間を予測できるでしょう。
そこで、この節では生体内薬物動態を簡単な微分方程式を用いてモデル化してみましょう。
薬物の吸収や排泄過程に比較して、生体内分布が急速に起こるのであれば、薬物は生体内に均一に分布すると考えることができるので、血漿濃度は吸収と排泄のバランスに依存します。
例えば、静脈内に1回薬物を投与する場合を考えてみましょう。ここで体内の薬物量を$${x}$$、薬物濃度を$${C}$$、薬物投与量を$${X_0}$$とします。
薬物の排泄速度が、体内の薬物濃度$${C}$$に比例すると仮定すると、薬物の減少速度は次の一次反応で書けます。
$$
\frac{dx}{dt}=-kx
$$
ここで、$${k}$$は排泄速度定数です。これは変数分離法で簡単に解くことができ
$$
x=X_0 e^{-kx}
$$
を得られます。薬物の分布容積を$${V}$$とすると、薬物濃度$${C}$$は
$$
C=C_0 e^{-kx}
$$
と表せます。ただし、$${C_0=\frac{X_0}{V}}$$です。この関数のグラフは片対数プロットを行うと、直線となるため、例えば血漿内薬物濃度を経時的に測定したデータをプロットし、外挿を施すことで$${C_0}$$や$${k}$$、その薬物の半減期などを推定することができます(ただし、その薬物が上記のような単一コンパートメントモデルで近似できる薬物動態を示す場合に限ります)。
また、別の状況として、持続点滴静脈注射を行う場合も考えられます。そこで、持続点滴静注速度を$${R}$$とすると、薬物の血漿濃度の時間変化は次の式で表されます。
$$
\frac{dx}{dt}=R-kx
$$
この解を、濃度$${C=\frac{x}{V}}$$に代入すると
$$
C=\frac{R}{kV}(1-e^{-kx})
$$
を得られます。投与開始から十分時間がたった時点での薬物濃度は
$$
C=\frac{R}{kV}
$$
と表せます。このことから、薬物の分布容積と排泄速度定数がわかれば、定常状態である血漿濃度を達成するにはどの程度の注入速度が必要かを求めることができます。
形態形成の数理
Turing Pattern
この項目に関しては以前、記事を書いたので、ご興味があればぜひご覧ください。
Clock and Wavefront モデル
ヒトをはじめとする脊椎動物の発生過程において一時的に形成される、非常に重要な構造として「体節」というものがあります。体節は特徴的な繰り返し構造を有し、それをもとに脊椎骨や骨格筋が形成されることが知られています。そのメカニズムについては古くからClock and Wavefrontモデルが提唱されています。
このモデルでは、まだ分節化していない中胚葉のそれぞれの細胞が持つ分子時計が作り出す時間的な振動が、空間的な周期性へ変換されると考えられています。さらに、後の実験でその分子時計の実体となるタンパク質が同定されました。しかし、この時間的な周期性というアナログな情報を、体節における空間的な分節というデジタルな情報に変換する機構についてはわかっていませんでした。
しかし、昨年ゼブラフィッシュにおいて、この仕組みを実現する分子機構の同定およびコンピュータシミュレーションによる体節形成の再現に成功したとする論文[7]が発表されました。紙面の都合上、ここでは詳細を説明することができないため、興味がおありの方は以下の論文を参照してください。
人体の動的特性
生命現象に見られる制御機構
酸素解離度の非線形性
ヘモグロビンは赤血球に含まれる鉄結合タンパク質であり、ヘムとグロビンという2種類の分子によって構成されています。この2つが鎖状に結合することでヘモグロビン鎖を形成します。同じ種類(グロビンには複数のタイプがあります)のヘモグロビン鎖が対をなし、さらにその2対が合わさった四量体としてヘモグロビン分子は存在しています。
ヘモグロビン鎖1本は酸素分子1つと結合するため、ヘモグロビン分子は合計4つの酸素分子と結合することで、血中において酸素を運搬する役割を果たしているというわけです。
ヘモグロビンに関する面白い性質として、環境の酸素分圧による酸素との結合が非線形に変化する、というものがあります。この性質のおかげで、ヘモグロビンは酸素分圧が高い肺において酸素とよく結合する一方で、酸素分圧の低い末梢組織においては酸素との結合度が低下するため、その飽和度の違いを利用して効率的に組織に酸素を供給することができます。
では、この酸素分圧に対する酸素飽和度の非線形な応答は、どのようにして生まれているのでしょうか。
ここで一般に、$${n}$$個の分子が同時に別の1分子に結合する高次反応を考えましょう。化学反応式は
$$
nA+B⇄A^n B
$$
と書けるので、
$$
\frac{d[A^n B]}{dt}=-k_b [A^n B]+k_f [A]^n [B]
$$
を得ます。ここで総和保存
$$
[B]+[A^n B]=1
$$
が成り立つとすると、平衡状態において
$$
[A^n B]=\frac{[A]^n}{K_d+[A]^n}
$$
を得ます。これはヒルの式と呼ばれます。横軸に$${[A]}$$、縦軸に$${[A^n B]}$$をとってグラフを描くとシグモイド曲線の形をとります。ここでnはヒル係数と呼ばれ、この値が大きくなればなるほど、よりシグモイド曲線の$${^{n}\sqrt{K_d}}$$あたりでの傾きが大きくなり閾値を持ったall-or-none反応に近づきます。
しかし、実際には1つの分子Bに分子Aが1つずつ結合して、最大でn個の分子Aが結合した状態をとる、という仮定の方が自然でしょう。ただし、このような場合でも、分子Aに協調性があれば、高次反応を示すことが知られています。つまり、1つ目の分子AがBに結合するときの$${K_d}$$値が、2つ目の$${K_d}$$よりもかなり大きいという状況は、分子Aの濃度を上げていったときに、1分子目はなかなか結合しないが、一旦1分子目が結合すれば2分子目以降はすんなり結合できる、ということを意味します。
実際に、ヘモグロビンにはこの性質があることが知られています。1つのヘムに酸素が結合すると、その変化が同じ分子内の別のヘムに伝わり(ヘム間相互作用)、他のヘムの酸素結合性が亢進するというわけです。これをホモトロピック作用と呼びます。
この仕組みによって、実際にヘモグロビンの酸素解離曲線のヒル係数は2.8という値になっています。
ちなみに、このように連続的な入力に対してスイッチ的に急激に変化して応答する性質のことを一般にウルトラ・センシティビティと呼びます[8,9,10,12]。
HPA系におけるホルモン調節
コルチゾールは抗ストレスホルモンとして知られています。生体の恒常性を乱すような心的・物理的ストレスにより、分泌が促進されます。コルチゾールの作用としては、糖新生、血圧の上昇、抗炎症作用などが挙げられます。
コルチゾールの分泌は、「視床下部」「下垂体前葉」「副腎皮質(束状帯)」という3つの組織によって調節されています。視床下部から分泌されたコルチコトロピンホルモン(CRH)が下垂体前葉に作用し、下垂体前葉から副腎皮質放出ホルモン(ACTH)が分泌されます。これが血流に乗って、副腎皮質へと作用し、その束状帯からコルチゾールが分泌されます。この一連の系はそれぞれの英語名(hypothalamus、pituitary、adrenal)の頭文字をとって、HPA系と呼ばれます。また、このHPA系において、コルチゾールは逆に視床下部や下垂体前葉に作用して、上位のホルモン産生・分泌を抑制します。このフィードバックによって、コルチゾールの血中濃度が過剰にならないように調整されているというわけです。
コルチゾールは脂溶性であるため細胞膜を通過して、細胞質中の受容体に結合します。これが核内に移動して、標的遺伝子のプロモーター領域に結合し、標的遺伝子転写速度を調節することで、ストレスへの応答作用を発揮します。
しかし、ストレスが長期間続く場合、コルチゾールの慢性的な過剰分泌が引き起こされ、肥満や高血圧、副腎肥大などの症状が現れることも知られています。
ではこのHPA系を、常微分方程式を用いてモデル化してみましょう。
ストレス刺激による脳から視床下部(H)の刺激を$${u(t)}$$、CRHの濃度を$${x_1 (t)}$$とすると、$${x_1 (t)}$$の時間変化は次の式でモデル化できるでしょう。
$$
\frac{dx_1}{dt}=\frac{q_1 Hu}{x_3} -α_1 x_1
$$
第1項の分母に含まれる$${x_3}$$は、$${x_3}$$(コルチゾール)によるCRHの産生の抑制を表しています。第1項がホルモンの産生を、第2項がホルモンの分解を表します。ここで、$${q_1}$$、$${α_1}$$、Hは定数としています。
同様に、視床下部(P)から分泌されるACTHの濃度を$${x_2 (t)}$$とすると、$${x_2 (t)}$$、$${x_3 (t)}$$の時間変化はそれぞれ
$$
\frac{dx_2}{dt}=\frac{q_2 Px_1}{x_3} -α_2 x_2\\
\frac{dx_3}{dt}=q_3 Ax_2-α_3 x_3
$$
と表せます。ここで、$${q_2}$$、$${α_2}$$、$${q_3}$$、$${α_3}$$、$${P}$$、$${A}$$は定数としています。
以上の3つの微分方程式を用いて、この系の定常状態におけるそれぞれのホルモンの濃度は
$$
x_3^{st}=(\frac{q_1 q_2 q_3}{α_1 α_2 α_3} uHPA)^{1/3} u^{1/3} P^{1/3} A^{1/3}\\
x_2^{st}≃u^{1/3} P^{-2/3} A^{1/3}\\
x_3^{st}≃u^{2/3} P^{-1/3} A^{-1/3}
$$
と計算できます。
しかし、このモデルには1つ問題点があります。このモデルは、分単位もしくは時間単位において、ストレス刺激を受けた場合にホルモンのレベルが上昇し、ストレス刺激提示が終わると元の値に戻るという現象をよく再現している一方で、妊娠、拒食症、薬物乱用、ステロイドの長期服用などの継続的なストレス刺激のあとではホルモンレベルが元の値に戻るのに月単位の時間がかかる、という現象を再現できていません。
そこで、器官のサイズの制御も考慮に入れたモデルが提唱されています。つまり、上の式のPとAを、それぞれ視床下部と副腎皮質の大きさを表す変数であると捉え直した上で、
$$
\frac{dP}{dt}=P(b_p x_1-a_p)\\
\frac{dx_3}{dt}=A(b_A x_2-a_A)
$$
という制御関係を仮定することで、週~月スケールでのストレス刺激に対するコルチゾール分泌量の変化が説明されました。
ファガーソン反射
下垂体後葉から分泌されるホルモンであるオキシトシンは、子宮平滑筋の収縮を促す効果があります。特に、分娩の際には、産道を胎児が通過する過程で子宮壁が伸展されると、反射的にオキシトシンの分泌が促進されることが知られており、これによってさらに子宮平滑筋が収縮されます。この反射はファガーソン反射と呼ばれ、人体におけるポジティブフィードバックの例であると捉えることができます[12]。
システム医学をはじめよう
システム分子医学
実は「生命を(もしくは疾患を)システムとして捉える」という理念は、1960年代には既に存在しており、ウィーナーが提唱した『サイバネティクス』という学際領域に代表されます。そこでは、通信理論や制御理論を導入することで、生命は個々の要素を統合するシステム的原理によって構成される、という捉え方が強調され、特にその後の生理学において、臓器や組織レベルでのホルモン調節・血圧調節系などのホメオスタシス機構に関する研究の発端となったと言えます(このような生理学の見方を「システム生理学」と呼ぶこともあります)。しかし、この時期のシステム的アプローチは、主に技術的な限界により、当時は臓器や組織のレベルの研究に留まりました。
そのような中、80~90年代ごろにはDNA、遺伝子、タンパク質についての分子生物学の研究が盛んとなり、対応する遺伝子の変異が存在するとほぼ確実に発症するという、単遺伝子性疾患の疾患原因遺伝子が同定されるなど、分子医学が急速に発展しました。
2000年代に突入すると徐々に網羅的な分子情報の獲得が技術的に可能となり始めました。その最たる例が「ヒトゲノム計画」でしょう。ゲノム情報以外にも、転写物やタンパク質、代謝物に加えエピゲノムに関するオミクス情報を取得する手法が次々と開発されました。これらの技術的な進歩と、それまでの分子生物学的研究において蓄積されてきた知見によって、システム医学的研究は、細胞レベルひいては分子レベルにおいて発展する基盤を獲得したと言えます。これは「システム分子医学」とも呼ぶことができ、その基本原理は「多くの疾患(特に多因子性)は、いくつかの遺伝子の変異やタンパク質の機能異常の効果は総合されて形成された分子パスウェイ/ネットワークの調節不全や歪みによって発症し進行する」ことを基本理念としてます[13]。
情報理論的アプローチ
前述の理念に基づき、システム分子医学では「細胞内の分子によって構成されるネットワーク」の同定が基本となります。一般に、この解析には各種オミクス情報が用いられ、遺伝子発現制御ネットワーク(GRN)やタンパク質相互作用ネットワーク(PPI)、代謝ネットワーク、シグナル伝達系など、細胞内のさまざまなスケールのネットワーク構造が解析対象となります。
そこで、このポスターでは情報理論において非常に重要な概念である相互情報量を応用したネットワーク同定手法について紹介します。
細胞における典型的な情報伝達機構は、「細胞外からの刺激が受容体を介して、細胞内のセカンドメッセンジャーに伝えられ、下流のシグナル伝達分子群に作用し、最終的には転写因子が特定の遺伝子群の発現を調節する」という流れに従います。
この時、細胞内のシグナル伝達機構はさまざまなノイズを受けることが考えられ、現実的には、外部からの刺激によって与えられる情報の全てが下流に存在する遺伝子群の発現調節にまで伝わるわけではないと考えられます。では、この情報伝達を定量的に評価するためにはどうしたら良いでしょうか。そこで有力な武器となるのが、情報理論という工学分野です。
情報理論の分野において非常に重要な概念の1つである相互情報量$${I}$$は以下のように定義されます。
$$
I(X;Y)=H(X)+H(Y)-H(X, Y)=H(X)-H(X | Y)
$$
ここで$${X}$$、$${Y}$$はそれぞれ離散型確率変数です。また、$${H(X)}$$、$${H(Y)}$$はそれぞれの確率変数のエントロピーと呼ばれる値であり、それぞれの確率分布を$${p(x)}$$、$${p(y)}$$とおくと、以下のように定義されます。
$$
H(X)=-Σ_x p(x)\log{p(x)}\\
H(Y)=-Σ_y p(y)\log{p(y)}
$$
$${H(X, Y)}$$は、$${X}$$、$${Y}$$の結合分布$${p(x,y)}$$に対するエントロピーであり
$$
H(X,Y)=-Σ_{x,y} p(x,y)\log{p(x,y)}
$$
とかけます。
また$${H(X | Y)}$$は$${X}$$が与えられた下での$${Y}$$の条件付きエントロピーであり
$$
H(X | Y)=-Σ_{x,y} p(x,y)\log{p(y | x)}
$$
と定式化されます。(第2式から第3式への以降はこれらの定義を駆使すれば導くことができます)
ここで、細胞外からの刺激を$${x}$$、標的分子(シグナル伝達経路の下流において情報を伝達される分子)の応答を$${y}$$とすると、$${H(X)}$$は外部からの刺激のもつ情報量、$${H(X | Y)}$$は通信路を介した情報伝達の際に失われる情報量(ノイズと捉えることもできる)であるため、相互情報量$${I(X;Y)}$$は実際に標的分子が受け取る情報量と解釈できるという訳です。
また、相互情報量$${I(X;Y)}$$は
$$
I(X;Y)=Σ_x p(x) Σ_y p(y|x)\log\frac{p(y│x)}{p(y)}
$$
と変形できることから、通信路$${p(y|x)}$$が与えられた下で、相互情報量は入力となる分布$${p(x)}$$に依存して決まることがわかり、$${p(x)}$$に関する相互情報量の最大値を通信路容量$${C}$$と呼びます。
$$
C=\underset{p(x)}{max}I(X;Y)
$$
このようにして、細胞内の情報伝達における情報の伝達を相互情報量などの定量的な指標を用いることができると、実験データから通信路である条件付き分布を推定することで情報伝達のネットワーク構造を推定することが可能になります。
実際に、Cheongら[15]は、培養細胞をTNFで刺激し、NF-κBの核移行とATF-2のリン酸化を測定し、TNFからNF-κBとATF-2への情報伝達の様式を調べました。
ここで、Cheongら[15]は情報を伝達するネットワーク構造の分岐に応じた二つの情報伝達のモデル、BushモデルとTreeモデルを提案し、入力分布に通信路容量を達成する分布を用いました。その結果、モデルを仮定しなかったときの相互情報量は、Treeモデルを仮定したときの相互情報量と非常に近く、Bushモデルの相互情報量はTNFと受容体複合体との相互情報量に近かったのです。ゆえに、TNF―NF-κBシグナル伝達の情報伝達は、ボトルネック状の構造をとっていると結論づけられました。
ネットワーク解析
ここでネットワーク理論(グラフ理論とも呼ばれます)に関する基本的な用語について説明します。ネットワークはノードと呼ばれる点とエッジと呼ばれる線の集合によって構成されます。ノード間をつなぐエッジが存在するとき、それら2つのノードは接続されているといいます。エッジに方向性がないネットワークを無向ネットワークといい、方向性があるものを有向ネットワークといいます。あるノードからエッジをたどって別のノードに移動し、そのノードからまたエッジをたどって別のノードに移動する、という操作を繰り返して得られるノードとエッジの列を歩道と呼び、始点と終点を除き、すべてのノードが異なるような歩道を経路と呼びます。経路におけるエッジの本数は経路長と呼ばれ、与えられた2つのノードを始点、終点とする経路のうち経路長が最短のものを最短経路と呼びます。
ノードに接続するエッジの本数は次数と呼ばれ、ネットワークにおける次数の算術平均は平均次数といいます。また、ノード間に張ることができるエッジ数の最大値に対する実際のエッジ数の割合はグラフ密度と呼ばれます。このほかにも、エッジの端点間の次数の相関を示す指標である次数相関や、任意のノードに隣接するノード同士がどれだけ密につながっているかを示すクラスタ係数などの指標も用いられます。
システム的アプローチをとる生命科学研究においては、実験によって同定されたネットワークやそれまでの知識を参照することにより構成されたネットワークを利用し、そのネットワークにおいて特に重要な役割を果たすノードを探索することで薬剤標的分子や病態関連分子などを同定する試みが行われています。その中でも特によく用いられる手法が「中心性解析」です。
中心性解析に用いられる、ネットワーク構造から算出される指標(中心性指標)にはさまざまなものが存在し、それぞれの指標の意義や特性は異なります。代表的ないくつかの中心性指標を以下に示します。
実際に、これらの中心性指標を導入することにより、酵素の活性部位の予測や脳ネットワークにおける病態と脳領域の中心性との関連についての研究、タンパク質相互作用ネットワークにおける必須タンパク質の同定などが行われています[16]。
また、ネットワーク構造を用いた解析手法として「ネットワーク可制御性」という概念も重要です。
ネットワーク可制御性解析は、与えられたネットワーク構造から、システムの振る舞いを制御するために重要なノード(このようなノードをドライバー・ノードという)を検出するために用いられ、必須遺伝子、疾患関連遺伝子、薬剤標的分子の推定などに用いられています。
そもそも「可制御性」という概念は制御理論に由来し、簡単に言えば「あるシステムにおいて、任意の初期状態から出発して、任意の時間に対して任意の目標状態に到達させることができるような、入力が存在すること」を指します。これは生命科学の文脈では、システムの状態が遺伝子発現パターンなどに対応すると捉えることができ、「可制御」であるとは、いくつかの遺伝子群の発現を制御することにより、そのネットワーク全体の発現パターンを制御することを通して、例えば、ある組織細胞から別の組織細胞に変化させたり、疾病状態から健康状態にしたりできることに対応します。
しかし、実はネットワーク構造のみではシステムの可制御性を判定することはできません。そこで役立つのが、構造可制御性という概念であり、この概念は、0でないほとんどのパラメータの値について、システムが可制御となることを意味します。また、あるネットワークが構造可制御性を満たすために、制御する必要のあるノード(ドライバー・ノード)の求め方も知られており、これによって必須遺伝子や疾患関連遺伝子の同定、脳構造ネットワークへの応用など、さまざまな試みが行われています[16,17]。
さらに学びたい方へ
近年では、機械学習ブームに伴い、基礎研究のみならず臨床現場においてもAIを活用した創薬・診断・治療などが盛んに行われつつあります。AIに限らず、このような数理工学の知見を基盤とした生命科学・医学分野への応用は、今後ますます発展を遂げることでしょう。
本記事にてご紹介した内容は、システム生物学・医学分野のごく一部にすぎません。これらの分野についてもっと勉強したい!という方のために、以下におすすめの教科書を列挙させていただきます。ぜひ参考にしていただければ幸いです。
鳥谷部祥一 (2022).『生物物理学』.日本評論社.
生物物理学に関する主要な話題に一通り触れることができ ます。金子邦彦, 澤井哲, 高木拓明, 古澤力 (2020).『細胞の理論生物学』.東京大学出版会.
東大の講義を元にして作られた理論生物学に関する代表 的な本です。畠山哲央, 姫岡優介 (2023).『システム生物学入門』.講談社.
システム生物学に「入門」するにはうってつけの本です。浜田道昭, 宇田新介 (2022).『システムバイオロジー』.コロナ社.
ある程度の数理工学の知識を必要としますが、確率過程や 情報理論に関する話題を取り扱っているという点で興味深 いです。近藤滋, 北野宏明, 金子邦彦, 黒田真也 (2010).『システムバイオロジー』.岩波書店.
こちらはシステム生物学に関するオムニバス的な書籍です。
参考文献一覧
医学と数理への招待
[1] 矢嶋信男 (2019).『常微分方程式』.岩波書店.
[2] 小林徹也先生の記事『システム生物学って何だったんですか?(前編)』(https://zenn.dev/crmind/articles/4f2fde14c808c2)
薬理学への応用
[3] 飯尾正光(監修) (2021).『標準薬理学 第8版』.医学書院.
形態形成の数理
[4] 金子邦彦, 澤井哲, 高木拓明, 古澤力 (2020).『細胞の理論生物学』.東京大学出版会.
[5] 近藤滋, 北野宏明, 金子邦彦, 黒田真也 (2010).『システムバイオロジー』.岩波書店.
[6] 佐藤純 (2021).『いますぐ始める数理生命科学』.コロナ社.
[7] Yabe, T., Uriu, K. & Takada, S. Ripply suppresses Tbx6 to induce dynamic-to-static conversion in somite segmentation. Nat Commun 14, 2115 (2023). https://doi.org/10.1038/s41467-023-37745-w
人体の動的特性
[8] 近藤滋, 北野宏明, 金子邦彦, 黒田真也 (2010).『システムバイオロジー』.岩波書店.
[9] 畠山哲央, 姫岡優介 (2023).『システム生物学入門』.講談社.
[10] 浜田道昭, 宇田新介 (2022).『システムバイオロジー』.コロナ社.
[11] Uri Alon, Systems Medicine (2024)
[12] 本間研一(監修) (2019).『標準生理学 第9版』.医学書院.
システム医学をはじめよう
[13] 田中博 (2012).『先制医療と創薬のための疾患システムバイオロジー』.培風館.
[14] 浜田道昭, 宇田新介 (2022).『システムバイオロジー』.コロナ社
[15] Cheong R, Rhee A, Wang CJ, Nemenman I, Levchenko A. Informa- tion transduction capacity of noisy biochemical signaling networks. Science. 2011;334(6054):354-358.
[16] 浜田道昭, 竹本和広 (2021).『生物ネットワーク解析』.コロナ社
[17] ホセ・ナチェル, 阿久津達也 (2023).『複雑ネットワークと制御理論』