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

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

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

今回、その[チュートリアル](https://pysb.readthedocs.io/en/stable/tutorial.html)を行ったので、
自分のために簡易的にまとめておく。

# 環境
```
$ python --version
Python 3.10.0
$ pip freeze | grep pysb
pysb==1.13.2
```
(なお、インストールは`pip install pysb`で大まかになんとかなる)

# 内容
## tutorial_a.py
```Python:tutorial_a.py
from pysb import *

Model()
Monomer('A')
Parameter('k', 3.0)
Rule('synthesize_A', None >> A(), k)
```
```Python: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行から始める。
```Python
from pysb import *
Model()
```
なお、これでここより下に記述されたものが、`from tutorial_a import model`でまとめてインポートできるように自動的になる。

ここより下に、
* 必須
    * `Monomer`
    * `Parameter`
    * `Rule`
* 非必須
    * `Observable`
    * `Compartment`

を定義することになる。定義時、第一引数が変数名としてglobal領域にオブジェクト作成されるのは、特徴的な動作である。(see. [Self-export](https://pysb.readthedocs.io/en/stable/tutorial.html#self-export))

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

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

例：MEK
```Python:Monomer 'MEK'の定義(tutorial_b.py)
Monomer('MEK', ['s218', 's222', 'k'], {'s218': ['u', 'p'], 's222': ['u', 'p']})
```

## tutorial_b.py
```Python: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`は三つの部位があるが、同様。

### モデルの内容のインタラクティブな確認
```Python
>>> 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
```Python: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)
```Python: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での記載と似た形式ではある。

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

例：
```Python: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)を解くためには初期条件が必要である。
初期条件は下記のように指定できる。
```Python: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
```Python: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
```

<details><summary><code>t</code>の中身</summary><div>

```Python
>>> 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.        ])

```
</div></details>


<details><summary><code>yout['obsBid']</code>の中身</summary><div>

```python
>>> 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])
```
</div></details>

### Plotting
```Python: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()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1166238/1ef33fd8-cdf0-bf3d-b97f-649ceaa9d346.png)

### Graph Visualization
`.../pysb/tools/`にある`render_reactions.py`や`render_species.py`を使って、dotファイルを作成できる。
dotファイルからは`dot`コマンドを使うなどして表示できる。
```shell:反応のグラフを出力する
$ 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/`であった。
```shell:出力されたdotファイルから画像を作成する。
$ dot -T png mymodel.dot -o mymodel.png
```
![mymodel.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1166238/92612749-146b-5adf-acfe-68eb0e6b94e7.png)

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

```Python: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使用)
```Python: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](https://pubmed.ncbi.nlm.nih.gov/29028896/)を用いることもできる。

## Modules
Moduleのリファレンスについては、[PySB Modules Reference](https://pysb.readthedocs.io/en/stable/modules/index.html)を参照して欲しい。

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

# 出典
https://pysb.readthedocs.io/en/stable/tutorial.html#first-steps