読者です 読者をやめる 読者になる 読者になる

kivantium活動日記

プログラムを使っていろいろやります

COBRApyを使ったFlux Balance Analysis

生物の代謝を数学的にシミュレーションする手法の一つにFlux balance analysis(FBA, 日本語だと「代謝流束均衡解析」)というものがあります。一般的な代謝のシミュレーションでは速度定数を実験的に求めてから微分方程式を解く必要があり、実験の手間も計算量も膨大になります。一方、FBAでは化学反応をStoichiometric coefficient(日本語では「化学量論係数」)で表し、定常状態で最適な(例えば成長速度が最も早くなる)パラメータを線形計画法で求めるという大胆な簡略化を行うことで高速な解析を可能にしています。原理の説明はこのスライドが分かりやすかったです。

この記事ではWhat is flux balance analysis? : Nature Biotechnology : Nature Researchという解説記事のSupplementary Tutorialで説明されているMATLABCOBRA Toolbox 2.0を用いたプログラムを、PythonのCOBRApyという似たライブラリ上で動かしてみることでFBAのだいたいの感じをつかんでみることを目標とします。

インストール

まずはCOBRApyをインストールします。基本的にINSTALL.mdに従えば大丈夫です。
pipでインストールできるはずなのですが、Ubuntu 14.04では動かなかったのでGitHubからクローンしてインストールしました。

sudo apt-get install python-pip python-dev
sudo pip install --upgrade pip
sudo pip install cython
git clone https://github.com/opencobra/cobrapy.git
cd cobrapy
python setup.py develop --user

以下はOptionalですが、入れないとWarningが出てうるさいので入れたほうがいいです。scipyはpip一発では入らないのでいくつか依存パッケージを先にインストールしておく必要がありました。

sudo apt-get install build-essential gfortran libatlas-base-dev
sudo pip install numpy scipy python-libsbml 

インストールが終わったらpythonを起動し

from cobra.test import test_all
test_all()

を実行してFalseが返ってきたらインストール成功です。

あとはGetting Startedから順番にDocumentを読めばだいたいのことが分かると思います。

動かしてみる

まず、好気条件下で乾燥重量1 gのE. coliが1時間あたりグルコースを最大18.5 mmol吸収できるとき(原文の18.5 mmol gDW^{-1} h^{-1}はこういう意味だと思うのですが、間違っているかもしれません)の成長率を求めてみます。

# -*- coding: utf-8 -*-
import cobra.test

# E. coliの「教科書的な」コア代謝モデルを読み込む(「教科書的な」の意味が私にはよく分からない)
model = cobra.test.create_test_model("textbook")
# グルコースの最大吸収量を指定
model.reactions.get_by_id("EX_glc__D_e").lower_bound = -18.5
# 酸素の最大吸収量を指定(1000はどんなに吸っても枯渇しない十分大きな数くらいの意味)
model.reactions.get_by_id("EX_o2_e").lower_bound = -1000
# 最適化する目的関数をE. coliの成長速度に指定
model.objective = "Biomass_Ecoli_core"

# 最適化した結果を出力
print(model.optimize().f)

出力は1.65307185116となり、論文の1.6531/hrという値に一致します。

ちなみに、このシミュレーションの妥当性はStoichiometric Flux Balance Models Quantitatively Predict Growth and Metabolic By-Product Secretion in Wild-Type Escherichia coli W3110で検証されていますが、

the predictions they yield have not been directly experimentally verified

と書いてありました……。

同じ実験を嫌気環境下で行ったシミュレーションは

import cobra.test
model = cobra.test.create_test_model("textbook")
model.reactions.get_by_id("EX_o2_e").lower_bound = 0
model.reactions.get_by_id("EX_glc__D_e").lower_bound = -18.5

model.objective = "Biomass_Ecoli_core"

print(model.optimize().f)

となります。出力値は0.470565171089で、論文の0.4706/hrと一致します。


大腸菌に炭素源としてグルコースを与えず、代わりにコハク酸を与えた場合をシミュレーションします。

import cobra.test
model = cobra.test.create_test_model("textbook")
model.reactions.get_by_id("EX_succ_e").lower_bound = -20 
model.reactions.get_by_id("EX_glc__D_e").lower_bound = 0

model.objective = "Biomass_Ecoli_core"

print(model.optimize().f)

出力は0.840134176847でこれも論文の0.8401/hrと一致します。

以上のようにCOBRApyを用いることでFlux balance analysisを簡単に行うことができます。