はじめに¶
このページの目的は、Python/Jupyterを用いて遺伝子制御ネットワークやその進化過程のシミュレーションを行うための基礎を学ぶことです。そのため、Pythonの使い方全般を説明するのではなく、シミュレーションのコードを書けるようになるというゴールに向けて、最低限の知識を記載しています。ここに書かれていない使い方も多くあるので、適宜googleなどで検索をしながら必要な部分は補って下さい。また、質問やコメントはchikara.furusawa@riken.jpまでお願いします。
簡単な計算¶
Jypeterを立ち上げ、「File」→「New Notebook」→「Python3」とすると、In [ ]: と表示された灰色の行が出力されるはずです。これをセルと呼びます。
このセルに、式やコマンドを入力してコードを実行することが出来ます。
あと、右上にある「Code」とあるプルダウンメニューを「Makedown」とすると、Notebookにテキストを書き込むことが出来ます。 （このweb pageも、Jupyther Notebookにコードやテキストを入力して作成しています）
では、コードを走らせてみましょう。下のようにセルに数式を入力して、ShiftキーとReturnキーを同時に押すと、そのコードが実行されます。 （自分で入力しても良いですし、このページからコピー＆ペーストで入力してもかまいません。後半の長いコードはコピペが楽です）
3+2
上のコードは、電卓のようですが、「2+3」という式を評価して値を出力するということを行っています。以下、四則演算などの方法を示します。
（「#」以降はコメントとして無視されます。以降、コードの中に解説を入れるときに用います）
3-2 # 差
3*2 # 積
3/2 # 商
3**2 # べき乗
3 % 2 # 剰余
次に、変数というものを考えます。これは、値が入っている箱のようなものをイメージしてもらうと良いかと思います。下のコードでは、aとbという箱に値を代入して、それを演算しています。このとき、記号「=」は等号ではなく、「右辺の値を左辺の変数に代入する」という操作を意味します。
a = 3
b = 2
a + b
aという箱（変数）に3、bという箱に2を入れて、「a+b」のところで、その箱の中身を足すという計算をやっています。
あと、すでに使われている変数をセルに入力すると、その変数がどんな値を持っているかが表示されます。
a
配列（ベクトル）¶
次に、複数の変数が1次元に並んでいる配列と呼ばれる構造を扱います。この配列によって、ベクトルや行列などを表現します。Pythonでは、多くのデータをこの配列を用いて扱います。
Pythonにはリストと呼ばれるデータ構造など、配列を扱うための複数の方法がありますが、この解説では数値シミュレーションを実行することを主眼としているので、その中でもPythonで数値計算を行うためのライブラリである「numpy」を用いて配列を扱うことにります。
まず、今ある環境でnumpyが使えるようにする作業（インポートと呼ばれる）を以下のように行います。
import numpy as np
これは、numpyというライブラリをインポートして、以降は「np」という略号で用いるということを意味しています。これ以降に、「np.何とか」というコードが出てきたら、それはnumpyを用いたコードであるという意味になります。この操作は、1回やっておけば後で繰り返す必要はありません。
numpyでは、例えば以下のようにして1次元の配列（ベクトル）を作ります。要素の値が1, 2, 3である長さ3のベクトルです。
a = np.array([1,2,3])
配列の中身を知りたいときには、単にその配列の名前を入力して実行します。
a
配列もいろいろと演算をできます。
a*2 # 積
np.sum(a) # 要素の和（1+2+3=6）、平均値はnp.mean
np.max(a) # 最大値（最小値はnp.min）
len(a) # ベクトルの長さ
ベクトルを生成するときに、いちいち値を手で入力することは滅多にありません。その代わり、以下のコードをよく使います。
b = np.zeros(10)
b
c = np.arange(0,9)
c
np.zeros(n)は長さnで全ての要素がゼロのベクトルを出力します。上のOut[22]で「0.,」とゼロのあとにピリオドがついているのは、その配列のゼロが少数として扱われているという意味です。整数が良ければ、np.zeros(10, dtype='int')とします。
np.arange(x,y)はxからyまでの連続した値を要素としたベクトルを出力します。なにも指定をしなければ、値の間隔は1です。以下のように間隔を変えることも可能です（3つめの入力が間隔になります）。また、最初の入力xを省略して、np.arrange(y)とすると、0からy-1まで間隔1の値を出力します。
c = np.arange(0,3,0.1)
c
c = np.arange(10)
c
配列の値には以下のようにしてアクセスすることが出来ます。例えば、配列aのn番目の値を表示するためには、a[n-1]とします。なぜマイナス1が必要かというと、Pythonでは配列のインデックス（座標）がゼロから始まるからです。
a = np.arange(0,8,2)
a
a[2]
上の長さ4のベクトルの値はそれぞれ、a[0]=0、a[1]=2、a[2]=4、a[3]=6 となります。a[3]まででベクトルは終わりであることに注意してください。
また、特定のベクトルの要素に値を代入することもできます。
a = np.arange(0,8,2)
a
a[2]=10 #3番目の値を10に変える
a
配列（行列）¶
次に、2次元の行列を扱います。例えば3x3の行列は以下のように作成できます。
a = np.array([[0,1,2],[3,4,5],[6,7,8]])
a
以下の行列の要素にアクセスすることが出来ます。ここでも、3x3の行列の座標は行と列がそれぞれ0,1,2の範囲になることに気をつけてください。
a[1,1] #どの値が出るか、上の行列を見て考えてください。
また、行列も特定の要素の値を以下のように変えることが出来ます。
a = np.array([[0,1,2],[3,4,5],[6,7,8]])
a
a[1,1]=99
a
行列も自分で値を入力して作成することは少なく、例えば以下のコマンドで作成できます。
a = np.zeros((3,3)) #カッコが二重であることに注意してください。これには意味があるのですが、ちょっとここでは説明を省略します。
a
b = np.eye(3) #これは3x3の単位行列を出力します。
b
行列の積は、np.dot()を用います。行列とベクトルの積も計算できます。
a = np.array([[0,1,2],[3,4,5],[6,7,8]])
b = np.array([[1,2,3],[1,2,3],[1,2,3]])
c = np.dot(a,b)
c
行列から特定の行ベクトルや列ベクトルを取り出したい場合があります。例えば、3x3の行列から1行目を取り出したいときには、以下のようにします。
a = np.array([[0,1,2],[3,4,5],[6,7,8]])
a
a[0,:]
行列の範囲は縦横ともに0,1,2であることに注意してください。上のa[0,:]とすることにより、1行目（座標は0）の3つの値が取り出せています。このコロン「:」は、この場合はその行にある値を全部という意味になります。
繰り返し¶
同じ処理を複数回繰り返すことが簡単にできるのが、計算機プログラムの良いところです。Pythonでは以下のように繰り返しの処理を書きます。
for i in [0,1,2,3,4]: #変数iを0,1,2,3,4と変化させて、「:」以下の命令を実行する
print(i) # 変数iを表示する
上のコードは、以下のように書くこともできます。こちらの方が楽なので、普通はこちらも使います。
for i in np.arange(5): #これは、np.arange(5)が[0,1,2,3,4]というベクトルを出力することを用いてる
print(i)
繰り返しを用いると、たくさんの計算を行うことも簡単です。
total = 0
for i in np.arange(10):
total = total + i
print(total)
このとき、最後のprint()文をどこに置くかによって、出力は変わります。for文の繰り返しループの中にあるか外にあるかは、インデント（字下げ）で表現されます。上のコードと以下のコードの振る舞いの違いを考えてみてください。
total = 0
for i in np.arange(10):
total = total + i
print(total)
この場合は、print文はループの中にあります。 Jupyterは自動でインデントを行ってくれますが、入れ子のfor文を作ることも多いので、混乱をしないようにしてください。
条件分岐¶
プログラム中で、特定の条件が成り立ったときのみに命令を実行したい場合が良くあります。そうした場合にはif文を用います。
a=3
if a>2:
print('aは2より大きい')
else:
print('aは2より小さいか等しい')
もしif以下の条件文（この場合はa>2）が成り立てばif以下の命令を実行、成り立たなければelse以下の命令を実行します。
if文などで用いる条件文（比較演算子と呼ぶこともあります）は、以下のようなものがあります。
a == 10 # aが10に等しい
a != 10 #aが10に等しくない
a > 10 #aが10よりも大きい
a >= 10 #aが10よりも大きいか等しい
a < 10 #aが10よりも小さい
a <= 10 # aが10よりも小さいか等しい
複数の条件を論理演算子を用いて組み合わせて用いることも出来ます。論理演算子は、and、or、notがあります。
a=10
b=3
if a == 10 and b < 5:
print('aが10に等しく、かつbは5より小さい')
for文による繰り返しと、if文による条件分岐を組み合わせて用いることにより、いろいろなことが出来ます。 例えば以下は0から100までの値のうち、3でも5でも割り切れるものを出力するコードです。
for i in np.arange(100):
if i % 3 ==0 and i % 5==0: #記号%は余りを計算する。それが0であるということは、割り切れるということ
print(i)
自分で関数を定義する¶
一般にプログラミング言語では、さまざまな処理をまとめて関数を定義して用います。
関数とは大雑把にいえば、入力を関数という「ブラックボックス」に入れると、それに応じて何かが出力されるという構造のことです。
Pythonでは、以下のようにして自分で関数を定義することが出来ます。
def my_function( x ): #def という命令の後に、関数の名前を入れ、カッコの中に入力となる変数を指定する
y=x*2
return y #return の後に、出力となる変数を指定する
上のコードをセルに入れて、Shift+Returnキーを押すと、そのときは何も出力はされませんが、my_functionという関数が使えるようになります。
この関数my_functionは、iという入力を受けたときに、それを2倍して出力するという働きをします。伝統的に、入力のことを引数（ひきすう）、出力のことを返り値と呼んだりもします。
上のようにmy_functionを定義したら、それに入力をいれると、それに応じた出力を返します。
my_function(3)
といった感じです。この関数の場合、入力をベクトルにすると、出力もベクトルになります。
x=np.array([1,2,3])
my_function(x)
入力ベクトルの要素が2倍されて出力されています。また、入力をカンマで区切って複数にすることもできます。
def my_function2(a,b):
y=a**b
return y
my_function2(2,3)
関数の定義の中で、for文による繰り返しや、if文による条件分岐を使うことにより、複雑な機能を持つ関数を作れます。
例えば、以下は入力ベクトルの値の和を出力する関数です。
def my_sum(x):
sum=0 #最初に変数sumにゼロを入れておかないと、値が確定しない（と思う）
for j in np.arange(len(x)): #len(x)で入力ベクトルの長さとなるので、例えば長さ4のベクトルxを入力すると、forループは4回繰り返されて、そこでjは0,1,2,3と変化する。
sum=sum+x[j] #入力ベクトルxのj番目の値を、合計数sumに足す。
return sum
a=np.array([1,2,3,4])
my_sum(a)
入力ベクトルの標準偏差を求める関数であれば、例えば以下のようになります。データ数$N$のときの標準偏差$\sigma$は以下のように求められます。
$\sigma$=$\sqrt{\frac{1}{N} \sum_{i=1}^{N} (x_i-\bar{x})^2}$
ここで、$\bar{x}$は$x$の平均値です。なので関数の方でもまず平均値を求めた上で、上の和をとるようにします。
def my_sd(x):
av=0
for j in np.arange(len(x)): #このfor文は入力xの平均値を計算している
av=av+x[j]
av=av/len(x) #合計を要素数で割ることにより、変数avに平均値が入る
var=0
for j in np.arange(len(x)): #このfor文は入力xの分散を計算している
var=var+(x[j]-av)*(x[j]-av)
var=var/len(x) #これで変数varに分散が入る
return np.sqrt(var) #np.sqrt()という関数は、入力の平方根を返す関数なのでそれを利用している
x=np.array([1,2,3,4])
my_sd(x) #答えがあっているかの確認は、次節で説明しているnp.std()の出力と比較すればよい
(追記：12月12日)
関数を用いるときに、用いる変数が「関数の内側にあるか」「関数の外側にあるか」に気をつけて下さい。関数の内側で用いる変数（例えば、上の標準偏差を求める関数my_sd()の例では、avとかvarとか）は関数内部のローカル変数なので、関数の外から参照することは出来ません。逆に、関数の外で定義した変数は、関数の中から参照はできますが、その値を（基本的には）変えることができません。
(追記：12月15日)
一つの関数で、複数の値や配列を出力したい場合があります。そんな時は、関数の"return"の所に出力したい変数をカンマでつなげて書くだけで、複数の値や配列を出力することが出来ます。例えば、上のmy_sd()関数で、標準偏差と平均値を両方出力したいときには、次のような関数にします。
def my_sd(x):
av=0
for j in np.arange(len(x)): #このfor文は入力xの平均値を計算している
av=av+x[j]
av=av/len(x) #合計を要素数で割ることにより、変数avに平均値が入る
var=0
for j in np.arange(len(x)): #このfor文は入力xの分散を計算している
var=var+(x[j]-av)*(x[j]-av)
var=var/len(x) #これで変数varに分散が入る
return np.sqrt(var) , av #これで標準偏差と平均値が両方出力される
out=my_sd(x)
out
標準偏差と平均値（この場合2.5）がカンマでつながって出力されています。この場合、out[0]で標準偏差を、out[1]で平均値を取り出すことが出来ます。このようにカンマでつながった変数の組のことを、Pythonでは「タプル」と呼びます。タプルによる出力は変数の型が異なっていても良く、片方がベクトル、片方が配列とかでも使えます。
練習問題：以下のコードからコメントアウト#を取り除き、関数の空白部分に適当なコードを入れて、入力ベクトルxに含まれる値の最大値を返す関数を作成してください。ただし入力の値は常に正であるとします。
#def my_max(x):
#max=0
#for j in np.arange(len(x)):
###
### ここに何か入れてください
###
#return max
解答例は以下になります。
用意されている様々な関数¶
上では最大値や標準偏差を返す関数を作成しましたが、もちろんPythonやnumpyにはいろいろな関数が用意されています。例えば数学関数系では、以下のようなものがあります。
| 機能 | 関数 |
|---|---|
| 平方根 | np.sqrt(x) |
| べき乗 | np.power(x,y) |
| 絶対値 | np.fabs(x) |
| 自然対数, 常用対数 | np.log(x) , np.log10(x) |
| 三角関数 | np.sin(x), np.cos(x), np.tan(x) |
統計関係では、例えば以下のものがあります。
| 機能 | 関数 |
|---|---|
| 合計 | np.sum(x) |
| 平均値 | np.average(x) |
| 標準偏差 | np.std(x) |
入力xがベクトルの場合は、そのまま要素の平均値などを出力します。入力として行列を入れると、行あるいは列方向の平均値などをベクトルとして出力します。オプションとして、axis=0と指定すると行方向、axis=1と指定すると列方向の計算をします。例えば、
x = np.array([[0,1,2],[3,4,5]])
np.average(x,axis=0)
np.average(x,axis=1)
他にも、いろいろと便利な関数が用意されています。乱数はシミュレーションで使うので覚えておいてください。
| 機能 | 関数 |
|---|---|
| 1次元配列（ベクトル）の長さ | len(x) |
| 配列の次元数 | np.shape(x) |
| 0から1の範囲の一様分布に従う乱数（長さn） | np.random.rand(n) |
| 平均0、標準偏差1の正規分布に従う乱数（長さn） | np.random.randn(n) |
| 乱数の初期化（整数nを用いて初期化） | np.random.seed(n) |
グラフのプロット¶
Pythonでグラフを作成する方法は様々なものがありますが、ここは広く持ちられているライブラリであるmatplotlibを用いた方法を紹介します。
まず、以下のようにしてライブラリをインポートすると、以降はmatplotlibの関数をpltとして参照することになります。
import matplotlib.pyplot as plt
では、試しにグラフをプロットしてみましょう。xとして0から2$\pi$までを0.01刻みで増加させた（長さが629となる）ベクトルを用意し、yに関数sin(x)を代入します。その後に、横軸をx、縦軸をyとしてプロットしてみます。
x=np.arange(0,2*np.pi,0.01)#np.piはnumpyで用意されている定数（円周率）。この行で、xに0から2piまでの範囲を0.01刻みで分割した長さ629のベクトルが代入される。
y=np.sin(x)
plt.plot(x,y)
plt.plotは入力を2つとすると、最初の入力を横軸、次の入力を縦軸とします。入力を一つとすると、下のように横軸をベクトルのインデックス（何番目の要素か）、縦軸をその要素の値としてプロットします。
x=np.arange(0,2*np.pi,0.01)
y=np.sin(x)
plt.plot(y) #上のグラフと縦軸の値は同一だけど、横軸が異なっていることに注意
複数のグラフを書くときには以下のようにします。加えて、範囲やラベルの指定方法を説明しています。コメントを見ると理解できると思います。
x=np.arange(0,2*np.pi,0.01)
y1=np.sin(x)
y2=np.cos(x) #ここまでデータの用意
plt.plot(x,y1, label = "Data 1", color="r") #sin(x)の描画。"label"でデータのラベルを指定。"color"で色を指定
plt.plot(x,y2, label = "Data 2", color="b") #cos(x)の描画
#色の指定方法はいろいろとあるが、赤/青/緑/黄/シアン/マゼンダ/黒/白をr/b/g/y/c/m/k/wで指定可能
plt.legend() # 凡例を表示。場所は自動的に決められる。指定したければ、ple.legend(loc='upper right')などを用いる
plt.xlim(1.0,6.0) #x軸の範囲を最小値、最大値で指定
plt.ylim(-1.2,1.2) #y軸の範囲を最大値、最小値で指定
plt.title("Graph Title") #グラフのタイトル
plt.xlabel("X-axis") #x軸のタイトル
plt.ylabel("Y-axis") #y軸のタイトル
plt.show() #上の指定やオプションなどを全部ひっくるめてグラフを表示
線ではなくて点でデータをプロットしたいときには、以下のようにします。
x=np.arange(0,2*np.pi,0.01)
y=np.sin(x)+0.5*np.random.randn(len(x)) #sin(x)に、正規分布に従う乱数を足しています。xの要素数はlen(x)なので、その数の乱数を用意しています。
plt.plot(x,y,"o",color="r") #この"o"がプロットの形を指定しています。そのほか、"."（小さい点）、"+"、"x"、"*"などがあります。
その他、思いついたplt.plotのためのオプションです（適宜追加をします）。
| 機能 | 関数 |
|---|---|
| 対数プロット（x軸） | plt.xscale("log") |
| 対数プロット（y軸） | plt.yscale("log") |
| x軸の表示範囲の設定（xminからxmaxまで） | plt.xlim(xmin, xmax) |
| x軸の目盛りの配分（xminからxmaxまでをn-1分割） | plt.xticks(xmin, xmax, n) |
オイラー法による微分方程式の数値解法¶
微分方程式の振る舞いを数値的に解析する手法は様々なものがありますが、まずここでは最も簡単なものであるオイラー法を解説します。
$dx(t)/dt=f(x)$
の形の微分方程式を考えます。ここで$x$は状態を表すベクトルとし、$f(x)$は$x$のみで決まり時間には依存しないとします。
この微分方程式について、一つの初期条件からの時間発展を調べたいとします。しかし、微分扱うために必要な無限小の時間変化$dt$を計算機上で扱うのは難しいため、その代わりに有限の小さい時間変化$\Delta t$を使います。微分の定義を考えると、時間間隔が十分に短いという極限操作が必要ですが、$\Delta t$が十分に小さければ以下のように近似的に扱うことができるはずです。
$dx/dt = \lim_{\Delta t \rightarrow 0}\{x(t+\Delta t)-x(t)\}/\Delta t \fallingdotseq \{x(t+\Delta t)-x(t)\}/\Delta t = f(x)$
この近似後の式を変形すると、
$x(t+\Delta t) = x(t)+f(x)\Delta t$
となります。つまり、小さい時間間隔$\Delta t$後の状態$x(t+\Delta t)$は、今の状態$x(t)$に$f(x)\Delta t$を加えたものとなります。この近似的な時間発展方程式に従い、$x(t)$の時間発展を求める手法のことをオイラー法と呼びます。
注意点 :オイラー法は基本的に誤差が大きい（$\Delta t$のオーダーを持つ）手法です。ルンゲ・クッタ法など、より誤差の少ない方法も提案されており、広く用いられています。本入門では、最も簡単な手法としてオイラー法を採用しますが、本格的な研究を進める場合には、数値的な誤差に関する適切な扱いが必要になります。
では簡単な例として、以下の1変数の微分方程式の時間発展を考えます。
$dx/dt = a - bx$
これは、一定のレートで生成され、その濃度に応じて分解されていく化学反応のモデルと解釈することもできます。ここで、$a$と$b$は定数とし、$t=0$での初期条件を$x(0)=0$、また$a=1$、$b=0.1$とします。この微分方程式の時間発展を、$t=0$から$100$の範囲で、時間刻み$\Delta t=0.01$としてオイラー法によって求めるコード（の例）は以下になります。
def get_dx(x,a,b): #関数get_dxが入力xに対してその微分を返す
dx=a-b*x
return dx
a=1.0 #パラメータa
b=0.1 #パラメータb
dt=0.01 #時間刻み
t_max=100 #tの最大値
times=np.arange(0,t_max,dt) #timesは出力される時刻を要素としたベクトル
total_steps=len(times) #total_stepsは総ステップ数
x=np.zeros(total_steps) #この行で、長さtotal_stepsの出力用ベクトルxを用意しておく
x[0]=0 #xの初期値を最初の要素に代入
for i in np.arange(total_steps-1): #オイラー法のメインループ。繰り返し回数がtotal_steps-1であることに注意。
x[i+1]=x[i]+get_dx(x[i],a,b)*dt #オイラー法の計算
plt.plot(times,x) #横軸を時刻t、縦軸をxとしてグラフをプロット
定常値である$x=a/b$に近づいていく振る舞いが見れます。
もちろん、変数が複数ある場合にも、同じように数値的に振る舞いを解析することが出来ます。例えば以下の簡単な微分方程式
$dx_0/dt = a-bx_0$
$dx_1/dt = (x_0)^2-cx_1$
の時間発展をオイラー法で求めるコード例は以下になります（$t=0\sim40$、$\Delta t=0.01$としています）。
def get_dx(x,a,b,c):#関数get_dxが入力xに対してその微分を返す。ただし入力x、出力dxは要素数2のベクトルとしている
dx=np.zeros(2) #複数の値をベクトルで返すためには、まずベクトルを用意する必要がある
dx[0]=a-b*x[0]
dx[1]=x[0]*x[0]-c*x[1]
return dx
a=1.0 #パラメータa
b=0.1 #パラメータb
c=10.0 #パラメータc
dt=0.01 #時間刻み
t_max=40 #tの最大値
times=np.arange(0,t_max,dt) #timesは出力される時刻を要素としたベクトル
total_steps=len(times) #total_stepsは総ステップ数
x=np.zeros((total_steps,2)) #（total_steps行、2列）の出力用の2次元行列であるxを用意
x[0,0]=0 #xの初期値を最初の要素に代入。x[0,0]とx[0,1]の最初の要素は時刻、次の要素は変数のインデックス
x[0,1]=0
for i in np.arange(total_steps-1): #オイラー法のメインループ
x[i+1,:]=x[i,:]+get_dx(x[i,:],a,b,c)*dt #オイラー法の計算。ここでx[i,:]は、i番目の時刻の状態xを表す要素数2のベクトルになる
plt.plot(times,x[:,0],label="x_0")
plt.plot(times,x[:,1],label="x_1")
plt.legend()
plt.xlabel("time")
plt.ylabel("x")
plt.show()
$x_1$は$x_0$が上昇してきた後に立ち上がることが見て取れます。
練習問題1
2つの種（例えばウサギとオオカミ）が「被食者ー捕食者」の関係にある生態系の簡単なモデルとして、以下のLotka-Volterra方程式を考えます。
$dx_0/dt = ax_0-bx_0x_1$
$dx_1/dt = cx_0x_1-dx_1$
ここで、$x_0$と$x_1$はそれぞれ被食者と捕食者の個体数、$a,b,c,d$は正の定数パラメータとします。$bx_0x_1$の項は食べられることによる被食者の減少、$cx_0x_1$の項は食べることによる捕食者の増加をそれぞれ表しています。$a=b=c=d=1$、$\Delta t=0.001$、$t=0$での初期条件を$x_0=2$、 $x_1=1$としたとき、$t=0 \sim 30$の範囲の時間発展をプロットしてください。また、パラメータやモデルを適当に変えて、その振る舞いの変化を観察してください。
解答例は以下になります。
練習問題2（optional: これは当日に説明をする予定ですので、この問題に解答しておく必要はありません。ただし、どんなモデルを考えるかを理解しておいてください）
n個の遺伝子のタンパク質発現量を$\vec{x}=(x_0, x_1, \cdots, x_{n-1})$として、それらの間の制御ネットワークを以下のようにモデル化します。
$i$番目のタンパク質発現量の時間発展を、
$dx_i/dt = f(\sum_j W_{ji} x_j + \theta)- \gamma x_i$
で表せると仮定します。ただし、
$f(z)=1/\{1+\exp(\beta z)\}$
とします。ここで、$f(...)$の項はタンパク質の生成を、次項の$-\gamma x_i$は$\gamma$を定数としてタンパク質の分解・希釈を表すとします。$\beta=-10$とすると、$f=1/\{1+\exp(\beta z)\}$のグラフは$z$を横軸とすると以下のようになります。
def f(z):
return 1.0/(1.0+np.exp(-10*z))
z=np.arange(-2,2,0.01)
plt.plot(z,f(z))
図に示されるように、入力zが正の値をとれば1、負の値を取れば0となる階段のような関数になります（こうした形の関数はシグモイド関数と呼ばれます）。
さて、タンパク質の生成は他の遺伝子からの制御$W_{ji}$によって決まります。$W_{ji}$は発現制御を表す行列で、$i$番目のタンパク質が$j$番目のタンパク質の発現を活性化するときに$W_{ji}$は正の値をとり、抑制するときには負、相互作用がないときには0の値をとるとします。
上のモデル式では、$\sum_j W_{ji} x_j + \theta$が遺伝子$i$に対する発現制御の入力で、これが正になれば発現がオン、負になれば発現はオフになるというオンーオフ型の発現制御のモデルになっています。
このタンパク質間の制御ネットワークモデルについて、$n=10$とし、$W_{ji}$に$-1 \sim 1$のランダムな値を入れ、その振る舞いを出力するコードを作成してください。
数値計算ライブラリscipyを用いた微分方程式の数値解法¶
前節でのオイラー法のプログラムは構造が単純で理解がしやすいのですが、弱点が2つあります。
- 数値誤差が比較的大きい
- 実行時間が長い
このうち2については、前節のプログラムがforループを持ちていることが大きな要因になっています。基本的に、pythonのようにプログラムを逐次解釈しながら実行するインタープリタと呼ばれる言語は、forループを用いると実行時間が大きくなってしまう傾向があります。なので、プログラムを構築するときには、出来るだけ実行時間に影響を与えてしまうようなforループを避けることが重要になります。
上記の2つの問題を回避する方法として、すでに構築されている数値解析用のライブラリを用います。ここでは、scipyと呼ばれるライブラリを用います。
それを解説する準備として、まず前節のプログラムを少し変更して、scipyによるプログラムとの対比が簡単になるようにします。変更点は2点で、1つはget_dxの引数に時間刻みを表したベクトルtimesを加えたことです。これは、scipyの書式に合わせるためで、結果に影響は与えません。もう一つは、オイラー法でxの時間発展を求める関数my_odeを定義したことです。この関数の最初の引数は、関数get_dxであることに注意してください。このように、関数の引数（入力）に別の関数を用いることもできます。この関数my_odeは、微分方程式の右辺に対応する関数、xの初期条件、時間刻みを表したベクトル、微分方程式のパラメータを引数としています。扱っている微分方程式は前節の2遺伝子モデルと変わりません。
def get_dx(x,times,a,b,c):#関数getが入力xに対してその微分を返す。ただし入力は出力dxは任意の長さのベクトルとしている
dx=np.zeros(len(x)) #複数の値をベクトルで返すためには、まずベクトルを用意する必要がある。その長さをxの次元（len(x)）として求めている。
dx[0]=a-b*x[0]
dx[1]=x[0]*x[0]-c*x[1]
return dx
a=1.0 #パラメータa
b=0.1 #パラメータb
c=10.0 #パラメータc
dt=0.01 #時間刻み
t_max=40 #tの最大値
times=np.arange(0,t_max,dt) #timesは出力される時刻を要素としたベクトル
x_ini=np.array([0,0]) #x_iniは初期条件を与えるベクトル。今は要素数が2で、ともに0を代入している
def my_ode(get_dx,x_ini,times,a,b,c): #これが時間発展を求める関数
xx=np.zeros((len(times),len(x_ini))) #xxは、関数内で用いるために用意した出力用行列。要素数は、(時間刻みの数, xの次元)としている
xx[0,:]=x_ini #初期条件の代入
for i in np.arange(len(times)-1):
xx[i+1,:]=xx[i,:]+get_dx(xx[i,:],times,a,b,c)*dt #オイラー法の計算
return xx
上のセルを実行すると、関数my_odeを定義することができます。get_dxの中を変えることにより、微分方程式を数値的に解析するために、そこそこ汎用的に用いることが出来る関数となっています。
この関数より、例えば下のようにしてget_dxで定義している2遺伝子モデルの振る舞いを調べることが出来ます。
x=my_ode(get_dx,x_ini,times,a,b,c)
plt.plot(times,x[:,0],label="x_0")
plt.plot(times,x[:,1],label="x_1")
plt.legend()
plt.xlabel("time")
plt.ylabel("x")
plt.show()
ここまでは前節の結果と同じである。では、これをscipyを用いて書き換えてみます。まず、以下のようにしてscipyの中のodeintという関数を使えるようにします。
from scipy.integrate import odeint
次に、その関数odeintを用いて上と同じモデルの振る舞いを解析してみます。
x=odeint(get_dx,x_ini,times,args=(a,b,c)) #ここでmy_odeの代わりにodeintを用いてる。関数get_dxへの引数をargs=(...)で与える
plt.plot(times,x[:,0],label="x_0")
plt.plot(times,x[:,1],label="x_1")
plt.legend()
plt.xlabel("time")
plt.ylabel("x")
plt.show()
同じ結果が得られています。自作した関数my_ode()とscipyライブラリのodeint()の書式の違いは、後者で関数get_dx()への引数をargs=(...)という形で与えている部分だけです。このargs=(...)のカッコの内部は、変数の組を表すタプルという型になっています。odeint()のargs=(...)の部分は、常にタプルで渡す必要があり、関数get_dx()への引数が一つであっても、args=(a,)と最後にカンマをいれてタプルであることを明示する必要があります。
さて、同じ結果が得られるのであれば、どちらの関数でも良いと思うかもしれませんが、この2つの関数では実行に要する時間が全く違います。 それを定量的に調べるために、時間の最大値を表すパラメータt_maxを1000倍にして実行時間が長くして差が解りやすくなるようにします。
また、jupyterで実行時間を求めるためには、以下のようにtimeライブラリをimportした後に、"%time 関数"という書式を用います。
まず、t_maxを1000倍にしたときの関数my_ode()での実行時間を求めてみます。
t_max=40000 #tの最大値を1000倍にした
times=np.arange(0,t_max,dt) #timesは出力される時刻を要素としたベクトルをもう一度作っておく
import time #timeライブラリを使えるようにする
%time x=my_ode(get_dx,x_ini,times,a,b,c)
%time x=odeint(get_dx,x_ini,times,args=(a,b,c))
環境に依るとは思いますが、私のマシン（Thinkpad X1 Carbon 2017）では100倍くらい実行速度が違います。これは、関数odeint()の内部でforループのような繰り返し演算を用いるのではなく、行列計算をうまく使っているからです。
関数odeint()では、主にadams法と呼ばれる手法を用いて微分方程式を数値的に積分しています。これは、大雑把に言うと、オイラー法のように直前の時間ステップの情報のみから時間発展を求めるのではなく、それ以前の時刻における情報も用いて解を計算する方法となっています。
cythonやnumbaと呼ばれる部分的にC言語に変換したり、機械語にコンパイルする手法により、より高速化を図ることも出来ますが、今回の演習では用いないつもりです。興味のある人は調べてください。