PySBのexamplesを確認してみる(第2回 bax_pore_sequential.py & run_bax_pore_sequential.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を確認してみる(第1回 bax_pore.py & run_bax_pore.py)](https://qiita.com/yo314159265/items/ddd16d80ce7f999f4fc9)。）

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

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

```Python:bax_pore_sequential.py
"""A version of BAX pore assembly where subunits are added to the growing
complex one at a time (contrast with bax_pore.py). Also implements cargo
transport (Smac).
"""

from __future__ import print_function
from pysb import *
from pysb.macros import assemble_pore_sequential, pore_transport, pore_species

Model()

# s1,s2 are ring-formation sites; t is the transport (cargo binding) site
Monomer('Bax', ['s1', 's2', 't'])
Annotation(Bax, 'http://identifiers.org/uniprot/Q07812')
# loc is the location, (m)itochondrion or (c)ytosol; t is the transport site
Monomer('Smac', ['loc', 't'], {'loc': ['m','c']})
Annotation(Smac, 'http://identifiers.org/uniprot/Q9NR28')

Parameter('Bax_0', 8e4)
Parameter('Smac_0', 1e5)
for p in Bax_0, Smac_0:
    Annotation(p, 'http://identifiers.org/doi/10.1371/journal.pcbi.1002482',
               'isDescribedBy')

# Start with monomeric Bax and Smac in the mitochondrion
Initial(Bax(s1=None, s2=None, t=None), Bax_0)
Initial(Smac(loc='m', t=None), Smac_0)

# Maximum number of subunits in a pore
max_size = 6
# Size at which pores are "competent" to transport cargo
min_transport_size = 4

# Build handy rate "sets" (in this formulation, rates don't vary with pore size)
assembly_rates = [[2e-4, 1e-3]] * (max_size - 1)
transport_rates = [[3e-5, 1e-3, 1e1]] * (max_size - min_transport_size + 1)

# Assemble the pore
# (specify t=None so the pore can't fall apart while it's bound to cargo)
assemble_pore_sequential(Bax(t=None), 's1', 's2', max_size, assembly_rates)
# Transport Smac
pore_transport(Bax, 's1', 's2', 't', min_transport_size, max_size,
               Smac(loc='m'), 't', Smac(loc='c'), transport_rates)

# Add an observable for each pore size
for size in range(1, max_size + 1):
    Observable('Bax%d' % size, pore_species(Bax, 's1', 's2', size))
# Observe unbound Smac in each compartment
Observable('mSmac', Smac(loc='m', t=None))
Observable('cSmac', Smac(loc='c'))


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_sequential.py for example
usage.""")
```
</div></details>

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

```Python:run_bax_pore_sequential.py
"""Simulate the bax_pore_sequential model and plot the results."""

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

from bax_pore_sequential import model, max_size


# System is very stiff, and using logspace instead of linspace to produce the
# vector of time points happens to help with the integration
t = logspace(-3, 5) # 1e-3 to 1e5
print("Simulating...")
x = ScipyOdeSimulator(model).run(tspan=t).all

# Plot trajectory of each pore
for i in range(1, max_size + 1):
    observable = 'Bax%d' % i
    # Map pore size to the central 50% of the YlOrBr color map
    color = plt.cm.YlOrBr(float(i) / max_size / 2 + 0.25)
    plt.plot(t, x[observable], c=color, label=observable)
# Plot Smac species
plt.plot(t, x['mSmac'], c='magenta', label='mSmac')
plt.plot(t, x['cSmac'], c='cyan', label='cSmac')

# Expand the limits a bit to show the min/max levels clearly
plt.ylim([-0.01e5, 1.01e5])
# Show time on a log scale 
plt.xscale('log')
plt.legend(loc='upper right')
plt.show()
```
</div></details>

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

## 生物学的解説(BAX, Smacとは)
出典：[BAX(タンパク質)](https://ja.wikipedia.org/wiki/BAX_(%E3%82%BF%E3%83%B3%E3%83%91%E3%82%AF%E8%B3%AA))、[Smac](https://www.uniprot.org/uniprot/Q9NR28)

* BAX
    * アポトーシスの調節因子（促進）。タンパク質。
    * 普段は細胞質に存在するが、アポトーシスシグナル伝達に伴って、ミトコンドリア膜などに移動。
    * ミトコンドリアの電位依存性アニオンチャネルと相互作用し、膜に孔を形成する。（つまり死ぬ）
    * Bcl-2-associated X protein (Bcl-2結合Xタンパク質)
    * 今回のモデルでは、単量体〜六量体までを形成し、四量体以上が運搬能を持つ。
* Smac
    * アポトーシスを促進する因子。タンパク質。
    * Second mitochondria-derived activator of caspase
    * なお、推奨名称はDiablo IAP-binding mitochondrial protein
    * 今回のモデルでは、BAXに運ばれるだけ。


## コード解説
### マクロのimport
```Python:bax_pore_sequential.py(抜粋)
from pysb import *
from pysb.macros import assemble_pore_sequential, pore_transport, pore_species

Model()
```
`from pysb import *`, `Model()`はお約束。
今回、3つのマクロ(`assemble_pore_sequential`, `pore_transport`, `pore_species`)を読み込んでいる。
それぞれのリンク、使用例、定義は下記の通り。

[**assemble_pore_sequential**](https://pysb.readthedocs.io/en/stable/modules/macros.htmll#pysb.macros.assemble_pore_sequential)
```Python:assemble_pore_sequential(使用例)
Model()
Monomer('Unit', ['p1', 'p2'])
assemble_pore_sequential(Unit, 'p1', 'p2', 3, [[1e-4, 1e-1]] * 2)
```
```Python:assemble_pore_sequential(定義抜粋)
def assemble_pore_sequential(subunit, site1, site2, max_size, ktable):
    """Generate rules to assemble a circular homomeric pore sequentially.
    ....
    """
    return assemble_polymer_sequential(subunit, site1, site2, max_size, ktable,
                                       closed=True)
```
[**pore_transport**](https://pysb.readthedocs.io/en/stable/modules/macros.html#pysb.macros.pore_transport)
```Python:pore_transport(使用例)
Model()
Monomer('Unit', ['p1', 'p2', 'sc_site'])
Monomer('Cargo', ['c_site', 'loc'], {'loc':['mito', 'cyto']})
pore_transport(Unit, 'p1', 'p2', 'sc_site', 3, 3,
               Cargo(loc='mito'), 'c_site', Cargo(loc='cyto'),
               [[1e-4, 1e-1, 1]])
```
```Python:pore_transport(定義抜粋)
def pore_transport(subunit, sp_site1, sp_site2, sc_site, min_size, max_size,
                   csource, c_site, cdest, ktable):
    """
    Generate rules to transport cargo through a circular homomeric pore.
    ....
    """

    def pore_transport_rule_name(rule_expression, size):
        ....
        return '%s_%d_%s' % (_monomer_pattern_label(subunit), size,
                             _monomer_pattern_label(cargo))

    components = ComponentSet()
    # Set up some aliases that are invariant with pore size
    subunit_free = subunit({sc_site: None})
    csource_free = csource({c_site: None})
    # If cdest is actually a variant of csource, we need to explicitly say that
    # it is no longer bound to the pore
    if cdest().monomer is csource().monomer:
        cdest = cdest({c_site: None})

    for size, klist in zip(range(min_size, max_size + 1), ktable):
        pore_free = pore_species(subunit_free, sp_site1, sp_site2, size)
        pore_bound = pore_free.copy()
        cargo_bond_num = size + 1
        pore_bound.monomer_patterns[0].site_conditions[sc_site] = cargo_bond_num
        csource_bound = csource({c_site: cargo_bond_num})
        
        pc_complex = pore_bound % csource_bound

        
        name_func = functools.partial(pore_transport_rule_name, size=size)
        components |= _macro_rule('pore_transport_complex',
                                  pore_free + csource_free | pc_complex,
                                  klist[0:2], ['kf', 'kr'],
                                  name_func=name_func)
        components |= _macro_rule('pore_transport_dissociate',
                                  pc_complex >> pore_free + cdest,
                                  [klist[2]], ['kc'],
                                  name_func=name_func)

    return components
```
ちょっと複雑めなコードだった。

[**pore_species**](https://pysb.readthedocs.io/en/stable/_modules/pysb/macros.html)
```Python:pore_species(使用例)
Model()
Monomer('Unit', ['p1', 'p2'])
pore_tetramer = pore_species(Unit, 'p1', 'p2', 4)
```
```Python:pore_species(定義抜粋)
def pore_species(subunit, site1, site2, size):
    """
    Return a ComplexPattern representing a circular homomeric pore.
    """
    return polymer_species(subunit, site1, site2, size, closed=True)
```


### 単量体の定義(BAX, Smac)
```Python:bax_pore_sequential.py(抜粋)
# s1,s2 are ring-formation sites; t is the transport (cargo binding) site
Monomer('Bax', ['s1', 's2', 't'])
Annotation(Bax, 'http://identifiers.org/uniprot/Q07812')
# loc is the location, (m)itochondrion or (c)ytosol; t is the transport site
Monomer('Smac', ['loc', 't'], {'loc': ['m','c']})
Annotation(Smac, 'http://identifiers.org/uniprot/Q9NR28')
```
単量体`BAX`は3つの部位(`'s1'`, `'s2'`, `'t'`)を持つ。
単量体`Smac`は2つの部位(`'loc'`, `'t'`)を持ち、部位`'loc'`は二つの状態(`'m'`, `'c'`)を取りうる。(`'loc'`は部位というよりも存在場所という属性に近い；(m)itochondrion, (c)ytosol)

### 初期条件等の指定
```Python:bax_pore_sequential.py(抜粋)
Parameter('Bax_0', 8e4)
Parameter('Smac_0', 1e5)
for p in Bax_0, Smac_0:
    Annotation(p, 'http://identifiers.org/doi/10.1371/journal.pcbi.1002482',
               'isDescribedBy')

# Start with monomeric Bax and Smac in the mitochondrion
Initial(Bax(s1=None, s2=None, t=None), Bax_0)
Initial(Smac(loc='m', t=None), Smac_0)
```
初期条件`Bax_0`, `Smac_0`を指定し、登録する。

### 引数用の定数の準備
```Python:bax_pore_sequential.py(抜粋)
# Maximum number of subunits in a pore
max_size = 6
# Size at which pores are "competent" to transport cargo
min_transport_size = 4

# Build handy rate "sets" (in this formulation, rates don't vary with pore size)
assembly_rates = [[2e-4, 1e-3]] * (max_size - 1)
transport_rates = [[3e-5, 1e-3, 1e1]] * (max_size - min_transport_size + 1)
```
* 孔を形成するサブユニットの最大数
    * `max_size`
    * 6
* 孔が荷物を運べるようになる大きさ
    * `min_transport_size`
    * 4
* 速度定数の集合
    * `assembly_rates`
    * `[[2e-4, 1e-3], [2e-4, 1e-3], ..., [2e-4, 1e-3]]` (max_size - 1 繰返し)
    * `transport_rates`
    * `[[3e-5, 1e-3, 1e1], [3e-5, 1e-3, 1e1], ..., [3e-5, 1e-3, 1e1]]` (max_size - min_transport_size + 1 繰返し)

### 孔形成を定義(assemble_pore_sequential)
```Python:bax_pore_sequential.py(抜粋)
# Assemble the pore
# (specify t=None so the pore can't fall apart while it's bound to cargo)
assemble_pore_sequential(Bax(t=None), 's1', 's2', max_size, assembly_rates)
```
* 第一引数(`Bax(t=None)`)：孔を形成する単量体
* 第二引数(`'s1'`)：単量体同士の結合部位その１
* 第三引数(`'s2'`)：単量体同士の結合部位その２
* 第四引数(`max_size`)：孔を形成する最大のsubunit(単量体)数
* 第五引数(`assembly_rates`)：二量体形成からmax量体形成まで、それぞれのforward, reverseの速度定数のリスト

### 孔による輸送を定義(pore_transport)
```Python:bax_pore_sequential.py(抜粋)
# Transport Smac
pore_transport(Bax, 's1', 's2', 't', min_transport_size, max_size,
               Smac(loc='m'), 't', Smac(loc='c'), transport_rates)
```

* 第一引数(`Bax`)：孔を形成する単量体
* 第二引数(`'s1'`)：単量体同士の結合部位その１
* 第三引数(`'s2'`)：単量体同士の結合部位その２
* 第四引数(`'t'`)：孔を通過させる(運搬する)荷物と結合する部位
* 第五引数(`min_transport_size`)：運搬の起こる最小のsubunit(単量体)数
* 第六引数(`max_size`)：運搬の起こる最大のsubunit(単量体)数
* 第七引数(`Smac(loc='m')`)：運搬される荷物
* 第八引数(`'t'`)：運搬される荷物の、単量体と結合する部位
* 第九引数(`Smac(loc='c')`)：運搬後の荷物の状態
* 第十引数(`transport_rates`)：最小subunit数〜最大subunit数までのそれぞれのforward, reverse, catalytic(transport reaction)の速度定数のリスト

### 観察対象の指定
```Python:bax_pore_sequential.py(抜粋)
# Add an observable for each pore size
for size in range(1, max_size + 1):
    Observable('Bax%d' % size, pore_species(Bax, 's1', 's2', size))
# Observe unbound Smac in each compartment
Observable('mSmac', Smac(loc='m', t=None))
Observable('cSmac', Smac(loc='c'))
```
観察対象を指定する。
* `BAX{d}`：`d`個のBAXサブユニットからなる孔
* `mSmac`：ミトコンドリアにある遊離Smac
* `cSmac`：細胞質にある遊離Smac


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

from bax_pore_sequential import model, max_size
```
今回は変化が遅いので、時間には`linspace`ではなく、`logspace`を用いる。

### シミュレート
```Python:run_bax_pore_sequential.py(抜粋)
# System is very stiff, and using logspace instead of linspace to produce the
# vector of time points happens to help with the integration
t = logspace(-3, 5) # 1e-3 to 1e5
print("Simulating...")
x = ScipyOdeSimulator(model).run(tspan=t).all
```
`ScipyOdeSimulator(model).run(tspan=t)`でシミュレート。
`logspace`をタイムスパンに使用。

### 描画
```Python:run_bax_pore_sequential.py(抜粋)
# Plot trajectory of each pore
for i in range(1, max_size + 1):
    observable = 'Bax%d' % i
    # Map pore size to the central 50% of the YlOrBr color map
    color = plt.cm.YlOrBr(float(i) / max_size / 2 + 0.25)
    plt.plot(t, x[observable], c=color, label=observable)
# Plot Smac species
plt.plot(t, x['mSmac'], c='magenta', label='mSmac')
plt.plot(t, x['cSmac'], c='cyan', label='cSmac')

# Expand the limits a bit to show the min/max levels clearly
plt.ylim([-0.01e5, 1.01e5])
# Show time on a log scale 
plt.xscale('log')
plt.legend(loc='upper right')
plt.show()
```
色の指定は`color = plt.cm.YlOrBr(float(i) / max_size / 2 + 0.25)`が使える。
対数軸の指定は、`plt.xscale('log')`のように行える。


出力
![bax_pore_sequential.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1166238/e7d16532-c337-7f76-dbf7-481c009046a9.png)
徐々に大きな孔が形成され、ある時期はとても速くSmacがミトコンドリアから細胞質に移動している。

# 感想
初見ではよくわからなかったものがある程度わかるようになって嬉しい。
理解しながらまとめたので文章がもたついてしまっており、またいつか書き直したいとは思う。
macroの使用例が見れてよかった。

# 出典等
文中参照