PySBのexamplesを確認してみる(第1回 bax_pore.py & run_bax_pore.py)

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

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

今回、その[examples](https://github.com/pysb/pysb/tree/master/pysb/examples)の一部の内容を確認したので、
自分のために簡易的にまとめておく。

（次回は[PySBのexamplesを確認してみる(第2回 bax_pore_sequential.py & run_bax_pore_sequential.py)](https://qiita.com/yo314159265/items/178d93cb33958e2269e5)。）

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

# 内容
## コード内容（bax_pore.py）
<details><summary>コード内容（bax_pore.py）</summary><div>

```Python:bax_pore.py
"""A version of BAX pore assembly where the subunits dimerize and then
tetramerize instead of assembling sequentially (contrast with
bax_pore_sequential.py). Inhibition of pore formation by Mcl-1 is also
implemented.
"""

from __future__ import print_function
from pysb import *

Model()

# Each BAX-BAX bond must always involve a t1 site on one monomer and a
# t2 site on the other.
Monomer('BAX', ['t1', 't2', 'inh'])
Annotation(BAX, 'http://identifiers.org/uniprot/Q07812')
Monomer('MCL1', ['b'])
Annotation(MCL1, 'http://identifiers.org/uniprot/Q07820')

# Two lone monomers form a dimer.
Parameter('kdimf', 1e-6)
Parameter('kdimr', 1e-7)
Rule('bax_dim',
     BAX(t1=None, t2=None) + BAX(t1=None, t2=None) |
     BAX(t1=1, t2=None) % BAX(t1=None, t2=1),
     kdimf, kdimr)

# Two lone dimers form a tetramer, with a higher rate than the dimerization.
Parameter('ktetf', 1e-3)
Parameter('ktetr', 1e-4)
Rule('bax_tet',
     BAX(t1=1, t2=None) % BAX(t1=None, t2=1) + BAX(t1=2, t2=None) % BAX(t1=None, t2=2) |
     BAX(t1=1, t2=3) % BAX(t1=4, t2=1) % BAX(t1=2, t2=4) % BAX(t1=3, t2=2),
     ktetf, ktetr)

# An inhibitory protein can bind to a BAX subunit at any time.
Parameter('kbaxmcl1f', 1e-5)
Parameter('kbaxmcl1r', 1e-6)
Rule('bax_inh_mcl1',
     BAX(inh=None) + MCL1(b=None) |
     BAX(inh=1)    % MCL1(b=1),
     kbaxmcl1f, kbaxmcl1r)

# Initial conditions
Parameter('BAX_0', 8e4)
Initial(BAX(t1=None, t2=None, inh=None), BAX_0)
Parameter('MCL1_0', 2e4)
Initial(MCL1(b=None), MCL1_0)
for p in BAX_0, MCL1_0:
    Annotation(p, 'http://identifiers.org/doi/10.1371/journal.pcbi.1002482',
               'isDescribedBy')

# We must fully specify all four BAX-BAX bonds, otherwise the pattern
# is too loose, match a given species multiple times (beyond the
# factor of four expected due to the rotational symmetry of the
# tetramer), resulting in erroneously high values.
Observable('BAX4', BAX(t1=1, t2=3) % BAX(t1=4, t2=1) % BAX(t1=2, t2=4) % BAX(t1=3, t2=2))
# Same all-bonds requirement here.  However since the BAX tetramer is
# considered inhibited when even one subunit has an inhibitor bound,
# we only need to explicitly write inh=ANY on one of the monomer
# patterns.
Observable('BAX4_inh', BAX(inh=ANY, t1=1, t2=3) % BAX(t1=4, t2=1) % BAX(t1=2, t2=4) % BAX(t1=3, t2=2))


if __name__ == '__main__':
    print(__doc__, "\n", model)
    print("""
NOTE: This model code is designed to be imported and programatically
manipulated, not executed directly. The above output is merely a
diagnostic aid. Please see run_bax_pore.py for example usage.""")
```
</div></details>

## コード内容（run_bax_pore.py）
<details><summary>コード内容（run_bax_pore.py）</summary><div>

```Python:run_bax_pore.py
"""Simulate the bax_pore model and plot the results."""

from __future__ import print_function
import matplotlib.pyplot as plt
from numpy import linspace
from pysb.simulator import ScipyOdeSimulator

from bax_pore import model


t = linspace(0, 100)
print("Simulating...")
x = ScipyOdeSimulator(model).run(tspan=t).all

plt.plot(t, x['BAX4'])
plt.plot(t, x['BAX4_inh'])
plt.legend(['BAX4', 'BAX4_inh'], loc='upper left')
plt.show()
```
</div></details>

## コード構成
**bax_pore.py**でモデルを定義し、
**run_bax_pore.py**でそのモデルをimportし、シミュレート・出力する。

## 生物学的解説(BAX, MCL1, Bax poreとは)
出典：[Wikipedia: BAX(タンパク質)](https://ja.wikipedia.org/wiki/BAX_(%E3%82%BF%E3%83%B3%E3%83%91%E3%82%AF%E8%B3%AA))、[MCL1](https://en.wikipedia.org/wiki/MCL1)、[Bax Pore関連](https://www.nature.com/articles/s41598-017-02825-7)

* BAX
    * アポトーシスの調節因子（促進）。タンパク質。
    * 普段は細胞質に存在するが、アポトーシスシグナル伝達に伴って、ミトコンドリア膜などに移動。
    * ミトコンドリアの電位依存性アニオンチャネルと相互作用し、膜に孔を形成する。（つまり死ぬ）
    * Bcl-2-associated X protein (Bcl-2結合Xタンパク質)
    * 今回のモデルでは、単体、二量体、四量体を形成する。
* MCL1
    * アポトーシスの抑制因子。タンパク質。
    * Induced myeloid leukemia cell differentiation protein
    * 今回のモデルでは、Baxに結合することで、穴の形成を阻害する。
* BAX pore
    * poreは孔
    * 参考画像はこちら(https://www.nature.com/articles/s41598-017-02825-7/figures/1)


## コード解説
### 単量体の定義(BAX, MCL1)
```Python:bax_pore.py(抜粋)
from pysb import *

Model()

# Each BAX-BAX bond must always involve a t1 site on one monomer and a
# t2 site on the other.
Monomer('BAX', ['t1', 't2', 'inh'])
Annotation(BAX, 'http://identifiers.org/uniprot/Q07812')
Monomer('MCL1', ['b'])
Annotation(MCL1, 'http://identifiers.org/uniprot/Q07820')
```
`from pysb import *`, `Model()`はお約束。
Monomerの`BAX`は3つの部位(`'t1'`, `'t2'`, `'inh'`)を持つ。
Monomerの`MCL1`は1つの部位(`'b'`)を持つ。
`Annotation()`でアノテーションを入れられる。

### 二量体形成(単量体＋単量体＝二量体)
```Python:bax_pore.py(抜粋)
# Two lone monomers form a dimer.
Parameter('kdimf', 1e-6)
Parameter('kdimr', 1e-7)
Rule('bax_dim',
     BAX(t1=None, t2=None) + BAX(t1=None, t2=None) |
     BAX(t1=1, t2=None) % BAX(t1=None, t2=1),
     kdimf, kdimr)
```
片方の`'t1'`と片方の`'t2'`が結合する。
`kdimf`, `kdimr`の`k`は確かドイツ語由来だった気がする。`dim`はdimer(二量体)から、
`f`, `r`はおそらく**forward**, **backward**から。

### 四量体形成(二量体＋二量体＝四量体)
```Python:bax_pore.py(抜粋)
# Two lone dimers form a tetramer, with a higher rate than the dimerization.
Parameter('ktetf', 1e-3)
Parameter('ktetr', 1e-4)
Rule('bax_tet',
     BAX(t1=1, t2=None) % BAX(t1=None, t2=1) + BAX(t1=2, t2=None) % BAX(t1=None, t2=2) |
     BAX(t1=1, t2=3) % BAX(t1=4, t2=1) % BAX(t1=2, t2=4) % BAX(t1=3, t2=2),
     ktetf, ktetr)
```
空いていた`'t1'`,`'t2'`も結合に繋がり塞がることになる。
四量体は英語では**tetramer**。
なお、二量体形成反応より大分早い反応となっている。

### 抑制蛋白質の結合(BAX + MCL1 = BAX・MCL1)
```Python:bax_pore.py(抜粋)
# An inhibitory protein can bind to a BAX subunit at any time.
Parameter('kbaxmcl1f', 1e-5)
Parameter('kbaxmcl1r', 1e-6)
Rule('bax_inh_mcl1',
     BAX(inh=None) + MCL1(b=None) |
     BAX(inh=1)    % MCL1(b=1),
     kbaxmcl1f, kbaxmcl1r)
```
抑制蛋白`MCL1`は、`BAX`の`'inh'`にいつでも結合できる。

### 初期条件の指定
```Python:bax_pore.py(抜粋)
# Initial conditions
Parameter('BAX_0', 8e4)
Initial(BAX(t1=None, t2=None, inh=None), BAX_0)
Parameter('MCL1_0', 2e4)
Initial(MCL1(b=None), MCL1_0)
for p in BAX_0, MCL1_0:
    Annotation(p, 'http://identifiers.org/doi/10.1371/journal.pcbi.1002482',
               'isDescribedBy')
```
初期条件`BAX_0`と`MCL1_0`を定義し登録している。Annotationもつけている。

### 観察対象の指定
```Python:bax_pore.py(抜粋)
Observable('BAX4', BAX(t1=1, t2=3) % BAX(t1=4, t2=1) % BAX(t1=2, t2=4) % BAX(t1=3, t2=2))
Observable('BAX4_inh', BAX(inh=ANY, t1=1, t2=3) % BAX(t1=4, t2=1) % BAX(t1=2, t2=4) % BAX(t1=3, t2=2))
```
観察対象を指定する。
* `BAX4`：四量体
* `BAX4_inh`：4つのうち一つでも`'inh'`にMCL1が結合している四量体

`inh=ANY`は、相手が誰でも結合はしている状態を示す。(cf. [pysb.core.ANY](https://pysb.readthedocs.io/en/stable/modules/core.html?highlight=ANY#pysb.core.ANY))

### モデルの読み込み
```Python:run_bax_pore.py(抜粋)
import matplotlib.pyplot as plt
from numpy import linspace
from pysb.simulator import ScipyOdeSimulator

from bax_pore import model
```
`from bax_pore import model`で、
bax_pore.pyでの定義が自動的に格納された`model`をimportできる。

### シミュレート
```Python:run_bax_pore.py(抜粋)
t = linspace(0, 100)
print("Simulating...")
x = ScipyOdeSimulator(model).run(tspan=t).all
```
`ScipyOdeSimulator(model).run(tspan=t)`でシミュレート。

### 結果の描画
```Python:run_bax_pore.py(抜粋)
plt.plot(t, x['BAX4'])
plt.plot(t, x['BAX4_inh'])
plt.legend(['BAX4', 'BAX4_inh'], loc='upper left')
plt.show()
```
出力
![BaxPore.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1166238/15501c18-d1e8-b0d4-36b1-449611a0e0b2.png)

四量体形成が速いので四量体が蓄積される。
阻害部位への結合は、四量体形成を阻害するものではなく、`BAX`の状態に関係なく作用するので、大体一定の割合を占める。（ネガティブフィードバックはないという程度の意味）

# 感想
書き写しながらまとめることで、目があまり滑らなくてよかった。
訳語や解説等は調整が必要かとは思うが、気が向いたら書き直そうと思う。
`Annotation()`, `pysb.core.ANY`が今回初めて知ったものだった。

# 出典等
文中参照