Theoretical and experimental techniques
Building Executable Models from Descriptions of Biochemistry—Application to Patient-Specific Modeling
2023 Volume 63 Issue 6 Pages 320-324
Details
2023 Volume 63 Issue 6 Pages 320-324
細胞内の反応ネットワークや反応速度，遺伝子発現情報などを取り込めるMechanisticかつDynamicな数理モデルは，細胞のシステムとしての振る舞いを予測・制御する上で最も有効なアプローチの一つである．本稿では，Pasmopy（Patient-Specific Modeling in Python）を用いた，テキストから数式を経由しない実行可能な数理モデルの生成および，がん公共オミクスデータを用いた患者固有モデルへの応用について紹介する．
近年のオミクスデータの観測技術の目覚ましい発展により，私たちは細胞状態を記述する大規模な情報を手にすることが可能となった．これらのビッグデータを解析するソフトウェアのユーザビリティも向上しており，今では誰もが細胞状態の分類を自由に行い可視化できる時代である．こうした背景を踏まえ，生物学における今後の課題は，こうした細胞の分類だけでなく，分類の背後にある分子メカニズムを理解し，細胞のシステムとしての振る舞いを予測・制御することであると筆者は捉えている．このためには，細胞内の反応ネットワークや反応速度，遺伝子発現情報などを取り込めるMechanisticかつDynamicな数理モデルが最も有効なアプローチの一つである1)．
生化学反応系の数理モデリングの魅力は携わってきた人によって千差万別であると思われるが，一つひとつの反応式やその背後にある数学はシンプルであっても，細胞内の複雑な反応ネットワークの情報を載せると，直感での予測が難しい振る舞いを私たちに見せてくれる．そして対象となる生化学反応系に対してより深い洞察を与えてくれるところに私たちは惹きつけられるのではないだろうか．この魅力に取り憑かれるにつれ，筆者のなかでこの手法を色々な人に幅広く活用してほしいという思いが強くなった．筆者が大阪大学大学院理学研究科生物科学専攻の大学院生だった頃，私の周りでもバイオインフォマティクス解析を自身の研究に導入する友人が増えてきたが，数理モデリングというアプローチはあまり浸透していないことに気がついた．理由を聞いてみると，どうやらなんとなく難しい・とっつきにくいものと思われているらしい．それならば私がなんとかするしかないと奮起し，話し言葉・書き言葉に近い表現から数理モデルを構築する「Text2Model」を作成した（図1）．先述の万人に活用してほしいという思いから，このソフトウェアはオープンソースで公開した（https://github.com/pasmopy/pasmopy）．これは従来の常識を打ち破る数式を経由しない数理モデリングの手法であり，数理モデリングをより身近に，気軽なものに感じさせてくれるはずである．本稿では，このText2Modelによるシグナル伝達系の数理モデルの作成から，がん公共オミクスデータを用いた患者固有モデルへの応用について紹介する．
Text2Modelの概要．ユーザーは会話や論文中で用いられる自然言語に近い表現から実行可能な数理モデルを生成することができる．
Text2Modelは生化学的なイベントを記述する一文を生化学反応式へ変換し，物質収支を考慮して常微分方程式の右辺にフラックスを付加する2)．ここでは文章からどのように反応式を生成するかを具体的に説明したい．まず，Text2Modelにはあらかじめサポートされる結合・解離，リン酸化，転写，分解反応のような10数個のイベントが用意されており，第一段階としてユーザーが記述した文章がどのイベントに分類されるか振り分けることから始まる．この振り分けで用いられているのは自然言語処理ではなく古典的なパターンマッチングである．よってText2Modelの辞書に登録されていない表現は基本的に受け入れられないことになる（が，タイポのような場合にはおそらくユーザーが意図する表現をエラーメッセージ上で提案してくれる）．生化学モデリングでよく使われる表現は辞書にあらかじめ格納されているが，個別のユーザーごとに「この生化学イベントはこう呼びたい」という希望があれば，後から辞書に登録することもできる．
第二段階は反応式の生成である．ここで行われる操作は，文章中の主語・述語から反応物・生成物を同定し，それを反応式に組み込むことである．どのイベントでどの反応キネティクス（Mass actionやMichaelis-Menten式など）が使われるかはイベントごとに事前に決まっており，例えば非線形性の強い転写反応では，Hillの式が使われる．このようなルールとキネティクスがかちっと決まっていては柔軟にモデルを作成できないという懸念もあるかもしれないが，実際には，Text2Modelがデフォルトでサポートしない複雑なキネティクスを表現することができる機能も備えており，フレキシブルなモデルの作り込みが可能である．
テキストから作成したモデルがユーザの意図しないような動作を招かないために，Text2Modelには様々な機能が備わっている．特に構築した生化学反応モデルが熱力学的に実現可能であることは重要である．生化学反応系の挙動は化学量論に基づく質量流量とエネルギー流量によって決定される．前者は物質収支式から容易に把握できるが，エネルギー流量を正しくモデル化するためには，パラメータに関する特定の制約を遵守する必要がある．具体的には，生化学反応ネットワーク内に初期状態と終状態が同一の閉じた経路がある場合，平衡解離定数はいわゆる「詳細釣り合い」を満たす．この条件下では，経路内全体の自由エネルギー変化はゼロであるため，平衡解離定数の積が1となる制約が課される3)．したがって，この詳細釣り合いの関係は，独立に振る舞う平衡解離定数の数を減らすことになる．この制約を無視すると，本来独立でないパラメータを独立に扱うこととなり，例えば感度解析などにおいて誤解を招く結果を生み出す危険性が報告されている4)．Text2Modelではテキスト情報からネットワーク内の閉路を自動で探し出し，熱力学的な制約を表示することが可能である．特にモデル化するネットワークが大きくなってくると，ネットワーク図を描いて閉路を発見することが困難であり，テキスト情報のみから見つけてもらえるありがたみを実感できるはずである．Text2Modelは複数の先行研究モデルでテストされ，その動作が保証されている．詳しくは，Pasmopy2)のオンラインドキュメントを参照されたい．
筆者はこのText2Modelを使い，乳がんで重要なErbBファミリーの活性化からc-Myc活性化までのプロセスをモデル化した．本項目では，構築した数理モデルをどのように患者固有に個別化したかを説明する．筆者は患者固有モデルを構築する以前の研究でがん細胞株固有モデルを作成する研究を行っており5)，そこでモデルの個別化の手法について思索していた．数理モデルのダイナミクスを形作る要因として，反応パラメータ（例：タンパク質の結合・解離定数）と初期条件（タンパク質の量）があるが，前者は各細胞株の置かれた遺伝的コンテクストというよりは反応に関わるタンパク質の生化学的なプロパティに依存するとし，細胞株間で共通であると仮定した．一方で後者の初期条件は，これこそまさに変異によるタンパク質の異常な発現量の増加を表現しうる箇所であり，この部分に遺伝子発現情報を取り込むことで数理モデルを細胞株ごとに個別化できると予想した．そこで，手元にあった乳がん4細胞株のデータのうち，3細胞株のデータをパラメータの訓練に回し，残りの1細胞株をCancer Cell Line Encyclopedia（CCLE）6)から取得した遺伝子発現情報のみからそのシグナル伝達の活性化動態を予測できるかテストした．その結果，3細胞株でトレーニングしたパラメータを用いることで残りの1細胞株のシグナル伝達動態を正確に予測することができ，さらに，乳がん由来ではないHeLa細胞株の遺伝子発現情報を入力としても，その活性化動態を予測することができた（図2）．このことから筆者は，がんの種類を超えて，パラメータの変動よりも遺伝子発現情報の取り込みこそががん細胞の数理モデルの個別化において鍵となると結論した．
数理モデルの細胞株ごとの個別化．乳がん由来の3細胞株の遺伝子発現データとシグナル伝達動態の情報からパラメータを訓練し（①），新しい遺伝子発現データを入力として（②），その細胞内シグナル伝達動態を予測する．点が実験データ，実線がモデルの予測結果を表す．文献5より改変．
特に患者からの情報としては，TCGAデータベース7)から得られるのはスナップショットの遺伝子発現情報のみであり，患者固有のシグナル伝達動態は遺伝子発現データから推定する必要がある．そこで，患者固有モデルの作成においても，遺伝子発現レベルと細胞内シグナル伝達系の活性化動態の関係が既知であるがん細胞株の情報を用いて反応パラメータを訓練し，実際に観測することが困難な患者固有のシグナル伝達動態を出力することを行った．反応パラメータの決定に関して，ある反応速度を司るパラメータはシグナル伝達ネットワークを介して他の制御にも影響を及ぼすことから，一つひとつのパラメータの値を独立に決定することはできない．したがって一般に，モデルの規模が大きくなると実験データを再現しうるパラメータの推定が難しくなる．Pasmopyでは残差平方和による実験データとシミュレーションの誤差の評価，差分進化（大域的最適化手法の一つ）8)による最適化をデフォルトで提供してくれるため，この最適化問題を解決する上で助けになるはずである．必要に応じてパラメータ推定に用いるオプティマイザを任意のものに取り替えることも可能である．
377人分の乳がん患者固有シミュレーション結果から，シグナル伝達動態の特徴量でクラスタリングすると一般に予後不良で知られるトリプルネガティブ乳がんにおいて，予後を分類できることがわかった（図3）．この分類能は，従来用いられてきた遺伝子パネルや数理モデルに使用した遺伝子のスナップショットの情報に基づく分類よりも精確であることを確認している．
患者固有モデルによる乳がんの層別化結果．ヒートマップの各行が一患者固有モデルのシミュレーション結果を表し，それぞれの患者に紐づいた予後（Prognosis）およびホルモン受容体とHER2受容体の発現量に基づくサブタイプ分類結果（Subtype）が可視化される．ヒートマップの下部は，層別化に用いたダイナミクスの特徴量であり，ここではEGF，HRG刺激時の最大値を使用した．クラスター1と2にBasal-likeの患者（赤色）が顕著に集中しており，予後に関して分類できる（1：不良；2：良好）．文献2より改変．
次に筆者らは，この明確な予後における差異の背後にある分子メカニズムに興味を持ち，例えばこれらの患者群では薬剤応答が異なるのか患者固有モデルを用いて検証した．先行研究において，薬剤標的の探索に生化学反応系の数理モデルに感度解析を行い，高感度なタンパク質を調べるという手法が用いられており9)，患者ごとに数理モデルを所持している現在では個別に薬剤標的の探索ができる．そう考えて，予後の不良な患者群と良好な患者群（図3においてクラスター1と2に分類される患者群）において感度解析を実施した．その結果，これらの群ではEGFRの感度に差があり，予後の悪い患者群はEGFR阻害剤に対して薬効が低いと予想された．しかし患者由来のデータとしてTCGAが提供してくれるのは遺伝子発現プロファイル，予後，乳がんのサブタイプといった情報のみであり，各患者がキナーゼ阻害剤にどのように応答するかを知ることはできない．そこで，患者ではなくがん細胞株の薬剤応答データを用いてこの予測の検証を試みた．具体的にはまず，患者群の予後を分けた指標，ここではErbBファミリーの発現量比を用いて細胞株を分類した（図4左）．ここではEGFRの絶対的な発現量ではなく，ERBB3，ERBB3，ERBB4に対するEGFRの相対的な量比を分類指標とするため，EGFR阻害剤への応答を直観的に予測することは難しい．そのため，それらの患者群に対するEGFR阻害剤への感受性の比較として数理モデルが非常に強力なツールとなる．薬剤応答データ解析の結果はモデルの予測と一致し，予後の悪い患者群に対応する細胞株群においてはEGFR阻害剤への感受性が有意に低いことがわかった（図4右）．以上の結果から，患者固有モデルは予後の分類のみならず，薬剤応答の予測にも到達できる可能性を示した．
細胞株データによる薬剤応答予測の検証．クラスター1と2を分類した指標と同様の指標で細胞株を分類（左）．予後の悪い患者群に対応する細胞株（紫）はEGFR阻害剤（erlotinib, lapatinib）への感受性が低い（右）．文献2より改変．
本稿で筆者は，患者固有モデルによるトリプルネガティブ乳がん患者の層別化事例について報告した．しかし，もしも予後の分類を最終的な研究の出口として捉えていたならば，遺伝子発現データに適切な重みづけを施すことで本研究と同等以上の層別化は達成されたと考えられる．では，手間のかかるMechanisticなモデルを経由するメリットはどこにあったのだろうか．本研究の層別化手法も，内部での操作はモデリングを通じた遺伝子発現データの重みづけと言えるが，出力がシグナル伝達系の活性化動態という目にみえる形で現れるため，個別の遺伝子を超えてシステムレベルでの結果の解釈ができる．実際にこの研究では，予後の異なるトリプルネガティブ乳がん患者群の活性化動態から，ErbBファミリーの発現量比が重要であるという予測に到達できた．また，患者固有モデルの内部は，生化学反応ネットワークを表現する常微分方程式と反応速度を表現するパラメータという生物学的に解釈可能な形をしているため，一患者レベルで反応ネットワーク内の生化学プロセス摂動に対する応答予測ができる．モデルによる層別化でどのダイナミクスの特徴量を使用すると分類能が向上するかはがんのタイプによって異なり，ここが面白いところである．今後，様々ながんの患者群において，予後の分類を達成する特徴量の抽出とそれを生み出モデル構造へ立ち返り，どの分子メカニズムが鍵となるかを探る研究が進む．
本研究においてはモデルの個別化に遺伝子発現量のみを使用しており，例えば変異によるタンパク質間相互作用の変化は取り込まれない．KRASやBRAFなど，いくつかの発がん突然変異による平衡解離定数への影響が研究されており，こうした情報もモデルの個別化に取り込むのが今後の課題である．
先ほど筆者は，数理モデル構築は「手間のかかる」作業と述べたが，Text2Modelがこの煩わしさから解放してくれる一助になると考えている．これまでのモデリングは，数理モデリングの経験があるか，あるいはプログラミングができる研究者が主に行っていた．数理モデリングに限らず，計算生物学の手法を試そうと思う時，計算環境の構築や使用するプログラミング言語の文法の理解，よくわからないエラーメッセージの解釈など，「生物の研究っぽくない」ところで時間を取られてしまい，なかなか生物の研究のところに到達できないということがしばしばある．一方で，Text2Modelでは数理モデルに数式やコードは登場せず，数理モデルは生化学イベントの一連の説明文として可視化される．したがって，実験研究者が主体となって，モデルの構築やデバックに積極的に関わることを可能にする．本稿が，患者固有モデリングという手法を知るきっかけに，また患者固有モデルという形でなくても，数理モデリングに興味を持ち，トライしてみたいと思うきっかけになれば，それは望外の喜びである．
本研究の遂行にあたり，トランスクリプトームデータの統合および患者固有シミュレーション結果の分類に大きな貢献をしてくださった山城紗和氏に深く感謝する．
井元宏明（いもと ひろあき）
アイルランド国立大学ダブリン校，日本学術振興会海外特別研究員
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/63-6ie2326p01.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
岡田眞里子（おかだ まりこ）
大阪大学蛋白質研究所教授，東京農工大学農学研究科博士（農学）
;%0A%09%09%09newWindow.document.open();%0A%09%09%09newWindow.document.write('<img src=%22./Graphics/63-6ie2326p02.jpg%22>');%0A%09%09%09newWindow.document.close();%0A%09%09)
