PySBのチュートリアルを実行してみる

概要

PySBとはSystemBiologyを扱うためのPythonライブラリである。

document: https://pysb.readthedocs.io/en/stable/index.html
repository: https://github.com/pysb/pysb

今回、そのチュートリアルを行ったので、
自分のために簡易的にまとめておく。

環境

$ python --version
Python 3.10.0
$ pip freeze | grep pysb
pysb==1.13.2

(なお、インストールはpip install pysbで大まかになんとかなる)

内容

tutorial_a.py

tutorial_a.py

from pysb import *

Model()
Monomer('A')
Parameter('k', 3.0)
Rule('synthesize_A', None >> A(), k)

run_tutorial_a.py

from pysb.simulator import ScipyOdeSimulator
from tutorial_a import model

t = [0, 10, 20, 30, 40, 50, 60]
simulator = ScipyOdeSimulator(model, tspan=t)
simresult = simulator.run()
print(simresult.species)

• tutorial_a.pyでモデルを定義し、run_tutorial_a.pyでimportしシミュレーションしている。
• 元ページにあるfrom __future__ importはPython3では関係ないため無視した。
• 分子Aが1秒間に3コピー生成されるモデルを、0秒から60秒まで10秒ごとにAの個数をサンプルした。

Model作成

モデルを定義するには、下記の2行から始める。

from pysb import *
Model()

なお、これでここより下に記述されたものが、from tutorial_a import modelでまとめてインポートできるように自動的になる。

ここより下に、

• 必須
  • Monomer
  • Parameter
  • Rule
• 非必須
  • Observable
  • Compartment

を定義することになる。定義時、第一引数が変数名としてglobal領域にオブジェクト作成されるのは、特徴的な動作である。(see. Self-export)

Monomer(単量体)

Monomerは、モデル内では分割不可な要素（分子や化合物など）のこと。

• 第一引数：名前(name)
• 第二引数(オプション)：部位(sites)(結合したり状態変化したりするMonomer内の場所のList)
• 第三引数(オプション)：部位の状態(states)(部位ごとのとりうる状態のDict)

例：MEK

Monomer 'MEK'の定義(tutorial_b.py)

Monomer('MEK', ['s218', 's222', 'k'], {'s218': ['u', 'p'], 's222': ['u', 'p']})

tutorial_b.py

tutorial_b.py

from pysb import *

Model()
Monomer('Raf', ['s', 'k'], {'s': ['u', 'p']})
Monomer('MEK', ['s218', 's222', 'k'], {'s218': ['u', 'p'], 's222': ['u', 'p']})

• Rafには二つの部位('s','k')があり、's'については二つの状態('u', 'p')をとる。
• MEKは三つの部位があるが、同様。

モデルの内容のインタラクティブな確認

>>> from tutorial_b import model
>>> model.monomers
ComponentSet([
 Monomer('Raf', ['s', 'k'], {'s': ['u', 'p']}),
 Monomer('MEK', ['s218', 's222', 'k'], {'s218': ['u', 'p'], 's222': ['u', 'p']}),
 ])
>>> [m.name for m in model.monomers]
['Raf', 'MEK']
>>> model.monomers[0]
Monomer('Raf', ['s', 'k'], {'s': ['u', 'p']})
>>> model.monomers.keys()
['Raf', 'MEK']
>>> model.monomers['MEK']
Monomer('MEK', ['s218', 's222', 'k'], {'s218': ['u', 'p'], 's222': ['u', 'p']})
>>> model.monomers['MEK'].sites
['s218', 's222', 'k']

インタラクティブにも確認できる。モデルは自動的にmodelに収められている。

tutorial_c.py

tutorial_c.py

from pysb import *

Model()
Monomer('Raf', ['s', 'k'], {'s': ['u', 'p']})
Monomer('MEK', ['s218', 's222', 'k'], {'s218': ['u', 'p'], 's222': ['u', 'p']})
Parameter('kf', 1e-5)
Parameter('Raf_0', 7e4)
Parameter('MEK_0', 3e6)

• Parameter kfや初期条件Raf_0, MEK_0を定めている。

Parameter(定数)

Parameterは、生物学的な定数を示す定数。速度定数や初期値などに用いる。

• 第一引数：名前(name)
• 第二引数(オプション)：数値(value)(デフォルトは0)

注意点としては、単位が指定されないので、記述者が単位を揃えおく必要がある。（おすすめされているのは、分子数とS.I.単位を使うことである。）

mymodel.py(その1)

mymodel.py(その1)

from pysb import *

Model()
Monomer('C8', ['b'])
Monomer('Bid', ['b', 'S'], {'S': ['u', 't']})
Parameter('kf', 1.0e-07)
Parameter('kr', 1.0e-03)
Parameter('kc', 1.0)

Rule('C8_Bid_bind', C8(b=None) + Bid(b=None, S='u') | C8(b=1) % Bid(b=1, S='u'), kf, kr)
Rule('tBid_from_C8Bid', C8(b=1) % Bid(b=1, S='u') >> C8(b=None) + Bid(b=None, S='t'), kc)

Parameter('C8_0', 1000)
Parameter('Bid_0', 10000)
Initial(C8(b=None), C8_0)
Initial(Bid(b=None, S='u'), Bid_0)

Observable('obsC8', C8(b=None))
Observable('obsBid', Bid(b=None, S='u'))
Observable('obstBid', Bid(b=None, S='t'))

Ruleの記載が新しく追加されている。
C8とBidが一度結合し、その後解離してBidが活性化(S='t')されている。
Initialで初期条件を指定している。
Observableは後で見たい対象を指定している。

Rules(化学反応)

Rulesは分子や化合物の化学反応を定義する。Ruleの記載は、名前、反応パターン、生成のパターン、速度定数の指定、からなる。

BioNetGenやKappaでの記載と似た形式ではある。

• + は複合体形成を示す
• | は可逆反応を示す
• >>は不可逆反応を示す
• % は二つの物質が結合する反応を示す

例：

Ruleの定義例

Rule('C8_Bid_bind', C8(b=None) + Bid(b=None, S='u') | C8(b=1) % Bid(b=1, S='u'), kf, kr)

• 名前：'C8_Bid_bind'
• 反応パターン：C8 + Bid
• 生成パターン：（可逆反応で）C8 % Bid
• 速度定数：kfおよびkr

C8とBidが結合する反応。(bは結合を特定するindex。)

Observable(データ収集対象の指定)

Observableを使うことで、シミュレーション時のデータ収集対象を指定することができる。
（モデルが大きくなると対象は膨大になってしまうため）

Initial(初期条件の指定)

常微分方程式(ODE; Ordinary Differential Equation)を解くためには初期条件が必要である。
初期条件は下記のように指定できる。

mymodels.py(その1)抜粋 (初期条件の指定)

Parameter('C8_0', 1000)
Parameter('Bid_0', 10000)
Initial(C8(b=None), C8_0)
Initial(Bid(b=None, S='u'), Bid_0)

Simulation & Visualization

scipyとpylabも用いて、実際にシミュレーションしてみる

Simulation

Simulation (using interpreter)

>>> import mymodel as m
>>> from pysb.simulator import ScipyOdeSimulator
>>> import pylab as pl
>>> 
>>> t = pl.linspace(0, 20000)
>>> simres = ScipyOdeSimulator(m.model, tspan=t).run()
>>> yout = simres.all

tの中身

>>> t
array([     0.        ,    408.16326531,    816.32653061,   1224.48979592,
         1632.65306122,   2040.81632653,   2448.97959184,   2857.14285714,
         3265.30612245,   3673.46938776,   4081.63265306,   4489.79591837,
         4897.95918367,   5306.12244898,   5714.28571429,   6122.44897959,
         6530.6122449 ,   6938.7755102 ,   7346.93877551,   7755.10204082,
         8163.26530612,   8571.42857143,   8979.59183673,   9387.75510204,
         9795.91836735,  10204.08163265,  10612.24489796,  11020.40816327,
        11428.57142857,  11836.73469388,  12244.89795918,  12653.06122449,
        13061.2244898 ,  13469.3877551 ,  13877.55102041,  14285.71428571,
        14693.87755102,  15102.04081633,  15510.20408163,  15918.36734694,
        16326.53061224,  16734.69387755,  17142.85714286,  17551.02040816,
        17959.18367347,  18367.34693878,  18775.51020408,  19183.67346939,
        19591.83673469,  20000.        ])

yout['obsBid']の中身

>>> yout['obsBid']
array([10000.        ,  9600.82692793,  9217.57613337,  8849.61042582,
        8496.32045796,  8157.12260855,  7831.45589982,  7518.7808708 ,
        7218.58018014,  6930.35656027,  6653.63344844,  6387.95338333,
        6132.87596126,  5887.9786933 ,  5652.8553495 ,  5427.11687478,
        5210.38806188,  5002.31066362,  4802.53910592,  4610.74136092,
        4426.60062334,  4249.81001719,  4080.07733278,  3917.1205927 ,
        3760.66947203,  3610.46475238,  3466.25716389,  3327.80762075,
        3194.88629188,  3067.27263727,  2944.75491863,  2827.12948551,
        2714.20140557,  2605.78289392,  2501.69402243,  2401.76203172,
        2305.8208689 ,  2213.71139112,  2125.28052884,  2040.38151896,
        1958.87334783,  1880.62057855,  1805.49336521,  1733.36675338,
        1664.12107023,  1597.64120743,  1533.81668871,  1472.54158105,
        1413.71396601,  1357.23623273])

Plotting

Plotting (using simulation)

>>> pl.ion()
>>> pl.figure()
>>> pl.plot(t, yout['obsBid'], label="Bid")
>>> pl.plot(t, yout['obstBid'], label="tBid")
>>> pl.plot(t, yout['obsC8'], label="C8")
>>> pl.legend()
>>> pl.xlabel("Time (s)")
>>> pl.ylabel("Molecules/cell")
>>> pl.show()

Graph Visualization

.../pysb/tools/にあるrender_reactions.pyやrender_species.pyを使って、dotファイルを作成できる。
dotファイルからはdotコマンドを使うなどして表示できる。

反応のグラフを出力する

$ python [path-to-pysb]/pysb/tools/render_reactions.py [path-to-pysb-model-file]/mymodel.py > mymodel.dot

なお、[path-to-pysb]は私の環境では/Users/[myname]/.pyenv/versions/3.10.0/lib/python3.10/site-packages/であった。

出力されたdotファイルから画像を作成する。

$ dot -T png mymodel.dot -o mymodel.png

mymodel.py(その2)

Higher-order rulesを用いて記述をシンプルにする。
例えばcatalyze(Enz, Sub, "S", "i", "a", kf, kr, kc)のように。

mymodel.py(その2)

from pysb import *


def catalyze(enz, sub, site, state1, state2, kf, kr, kc):   # function call
    """2-step catalytic process"""                          # reaction name
    r1_name = '%s_assoc_%s' % (enz.name, sub.name)          # name of association reaction for rule
    r2_name = '%s_diss_%s' % (enz.name, sub.name)           # name of dissociation reaction for rule
    E = enz(b=None)                                         # define enzyme state in function
    S = sub({'b': None, site: state1})                      # define substrate state in function
    ES = enz(b=1) % sub({'b': 1, site: state1})             # define state of enzyme:substrate complex
    P = sub({'b': None, site: state2})                      # define state of product
    Rule(r1_name, E + S | ES, kf, kr)                       # rule for enzyme + substrate association (bidirectional)
    Rule(r2_name, ES >> E + P, kc)                          # rule for enzyme:substrate dissociation  (unidirectional)
   
# instantiate a model
Model()

# declare monomers
Monomer('C8', ['b'])
Monomer('Bid', ['b', 'S'], {'S':['u', 't']})

# input the parameter values
Parameter('kf', 1.0e-07)
Parameter('kr', 1.0e-03)
Parameter('kc', 1.0)

# OLD RULES
# Rule('C8_Bid_bind', C8(b=None) + Bid(b=None, S='u') | C8(b=1) % Bid(b=1, S='u'), *[kf, kr])
# Rule('tBid_from_C8Bid', C8(b=1) % Bid(b=1, S='u') >> C8(b=None) + Bid(b=None, S='t'), kc)
#
# NEW RULES
# Catalysis
catalyze(C8, Bid, 'S', 'u', 't', kf, kr, kc)


# initial conditions
Parameter('C8_0', 1000)
Parameter('Bid_0', 10000)
Initial(C8(b=None), C8_0)
Initial(Bid(b=None, S='u'), Bid_0)

# Observables
Observable('obsC8', C8(b=None))
Observable('obsBid', Bid(b=None, S='u'))
Observable('obstBid', Bid(b=None, S='t'))

Higher-order rules

上記ファイルのdef catalyze(を見るとニュアンスが分かる。
記述がシンプルに本質的になる。

提供されているmacroを使う。

.../pysb/macros.pyにその他の便利なmacroの記載がある。
下記に続く。

mymodels.py(その3 macro使用)

mymodels.py(その3 macro使用)

# import the pysb module and all its methods and functions
from pysb import *
from pysb.macros import *

# some functions to make life easy
site_name = 'b'
def catalyze_b(enz, sub, product, klist):
    """Alias for pysb.macros.catalyze with default binding site 'b'.
    """
    return catalyze(enz, site_name, sub, site_name, product, klist)

def bind_table_b(table):
    """Alias for pysb.macros.bind_table with default binding sites 'bf'.
    """
    return bind_table(table, site_name, site_name)

# Default forward, reverse, and catalytic rates
KF = 1e-6
KR = 1e-3
KC = 1

# Bid activation rates
bid_rates = [        1e-7, 1e-3, 1] #

# Bcl2 Inhibition Rates
bcl2_rates = [1.428571e-05, 1e-3] # 1.0e-6/v_mito

# instantiate a model
Model()

# declare monomers
Monomer('C8',    ['b'])
Monomer('Bid',   ['b', 'S'], {'S':['u', 't', 'm']})
Monomer('Bax',   ['b', 'S'], {'S':['i', 'a', 'm']})
Monomer('Bak',   ['b', 'S'], {'S':['i', 'a']})
Monomer('BclxL', ['b', 'S'], {'S':['c', 'm']})
Monomer('Bcl2', ['b'])
Monomer('Mcl1', ['b'])

# Activate Bid
catalyze_b(C8, Bid(S='u'), Bid(S='t'), [KF, KR, KC])

# Activate Bax/Bak
catalyze_b(Bid(S='m'), Bax(S='i'), Bax(S='m'), bid_rates)
catalyze_b(Bid(S='m'), Bak(S='i'), Bak(S='a'), bid_rates)

# Bid, Bax, BclxL "transport" to the membrane
equilibrate(Bid(b=None, S='t'),   Bid(b=None, S='m'), [1e-1, 1e-3])
equilibrate(Bax(b=None, S='m'),   Bax(b=None, S='a'), [1e-1, 1e-3])
equilibrate(BclxL(b=None, S='c'), BclxL(b=None, S='m'), [1e-1, 1e-3])


bind_table_b([[                  Bcl2,  BclxL(S='m'),       Mcl1],
              [Bid(S='m'), bcl2_rates,  bcl2_rates,   bcl2_rates],
              [Bax(S='a'), bcl2_rates,  bcl2_rates,         None],
              [Bak(S='a'),       None,  bcl2_rates,   bcl2_rates]])

# initial conditions
Parameter('C8_0',    1e4)
Parameter('Bid_0',   1e4)
Parameter('Bax_0',  .8e5)
Parameter('Bak_0',  .2e5)
Parameter('BclxL_0', 1e3)
Parameter('Bcl2_0',  1e3)
Parameter('Mcl1_0',  1e3)

Initial(C8(b=None), C8_0)
Initial(Bid(b=None, S='u'), Bid_0)
Initial(Bax(b=None, S='i'), Bax_0)
Initial(Bak(b=None, S='i'), Bak_0)
Initial(BclxL(b=None, S='c'), BclxL_0)
Initial(Bcl2(b=None), Bcl2_0)
Initial(Mcl1(b=None), Mcl1_0)

# Observables
Observable('obstBid', Bid(b=None, S='m'))
Observable('obsBax', Bax(b=None, S='a'))
Observable('obsBak', Bax(b=None, S='a'))
Observable('obsBaxBclxL', Bax(b=1, S='a')%BclxL(b=1, S='m'))

Compartments(区画)

Compartmentについては、｀mymodel_fxns.py`を見て欲しい。

Model calibration(較正)

Model calibrationにはPyDREAMを用いることもできる。

Modules

Moduleのリファレンスについては、PySB Modules Referenceを参照して欲しい。

感想

自分で見返すようにまとめたが、ひとまず速度重視でまとめてしまったので、後日訳語等の調整を行おうとは思う。

出典