1. Introduction
- そもそも、なぜモデルを作るのか？
- 現象を定量的に理解するため。（観測データの再現には、どのパラメータが重要？）
- 現象を予測するため。（面白い現象を見るには、どんな実験条件が適切？）
- 実験で確かめるのが困難なことを知りたいから。（コスト・技術・倫理などの問題で、そもそも実験できないことを知りたい）
- 構造は既知とする。パラメータはどうやって決める？
- 文献・データベースを参照（→実験条件がバラバラ。使って大丈夫？）
- 実験的に測る（→手間がかかる。そもそも測定できる？）
- 再現したい現象とモデルに依存して決める。
今日のテーマ
シグナル分子活性の経時変化のデータから、生化学反応の速度パラメータを推定する。
2. 細胞運命決定機構とRas, Rap1
- システム生物学の講義参照。
- ERKの時間波形が異なるだけで、異なる現象（増殖/分化）を制御
- 複雑なモデルをシンプルに記述する。（ERKの一過性活性化はRas（のみ）に、持続性活性化はRap1（のみ）に依存している。）
3. パラメータ推定
- パラメータの決め方
- 事象をうまく説明するモデルを作りたい。
- 客観的かつ定量的な指標（評価関数）を定義（ex.残差二乗和）
- 最適化
4. 最適化
- 勾配（微分）を用いた最適化：局所的な最適化に向いている。収束性や収束速度なども理論が進んでいる。
- 最急降下法
- 共役勾配法
- ニュートン法
- 信頼領域法
- など
- 確率を用いた最適化：大域的な最適化に向いている。ただし、得られた解の最適性は保証されない。
- Nelder-Mead法
- シミュレーテッドアニーリング（SA）
- タブー探索
- 遺伝的アルゴリズム（GA）
- 進化的プログラミング（EP）
- 粒子群最適化（PSO）
- など
5. 課題
In [1]:
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from scipy.io import loadmat
課題1：Ras, Rap1のODEモデルを作成¶
- Rap1モデルの $\mathrm{S},\mathrm{GEF_a}, \mathrm{GAP_a}, \mathrm{Rap1_a}$ の常微分方程式 (ODE) を記述せよ。
- Rasモデルの $\mathrm{S}, \mathrm{GEF_a}, \mathrm{GAP_a}$ RasaのODEを記述せよ。
- ※ 上の図では触媒を $\mathrm{pR}$ と記述しているが、以下では $\mathrm{S}$ と記述する。
- また、適当に反応速度定数 $k_i$ を定義している。（見れば対応関係はすぐに掴めるはず。）
- 不活性型 $\mathrm{hoge_i}$ と活性型 $\mathrm{hoge_a}$ の総和は保存されていると考える（簡単のため、$\mathrm{hoge_i} + \mathrm{hoge_a} = 1$ とおいている。）
Rap1モデル¶
$$ \begin{aligned} \frac{d\mathrm{S}}{dt} &= 0\\ \frac{d\mathrm{GEF_a}}{dt} &= k_1\left[\mathrm{S}\right]\left[\mathrm{GEF_i}\right] - k_2\left[\mathrm{GEF_a}\right]\\ &=k_1\left[\mathrm{S}\right]\left(1-\left[\mathrm{GEF_a}\right]\right) - k_2\left[\mathrm{GEF_a}\right]\\ \frac{d\mathrm{GAP_a}}{dt} &= 0\\ \frac{d\mathrm{Rap1_a}}{dt} &= k_3\left[\mathrm{GEF_a}\right]\left[\mathrm{Rap1_i}\right] - k_4\left[\mathrm{GAP_a}\right]\left[\mathrm{Rap1_a}\right]\\ &= k_3\left[\mathrm{GEF_a}\right]\left(1-\left[\mathrm{Rap1_a}\right]\right) - k_4\left[\mathrm{GAP_a}\right]\left[\mathrm{Rap1_a}\right] \end{aligned} $$モデルを関数化¶
In [2]:
# Rap model の微分方程式
def Rap1ODE(S,GEFa,GAPa,Rap1a,k1,k2,k3,k4,dt):
dS = 0 * dt
dGEFa = (k1*S*(1-GEFa) - k2*GEFa) * dt
dGAPa = 0 * dt
dRap1a = (k3*(1-Rap1a)*GEFa - k4*Rap1a*GAPa) * dt
return S+dS,GEFa+dGEFa,GAPa+dGAPa,Rap1a+dRap1a
In [3]:
def Rap1model(params,dt=0.01,min_t=0,max_t=100,
S=1,GEFa=0,GAPa=0.005,Rap1a=0,both=False):
max_t+=dt
k1,k2,k3,k4=params
time = np.arange(min_t,max_t,dt)
Rap1a_vals = np.zeros(len(time))
Gefa_vals = np.zeros(len(time))
for i,t in enumerate(time):
S,GEFa,GAPa,Rap1a = Rap1ODE(S,GEFa,GAPa,Rap1a,k1,k2,k3,k4,dt)
Rap1a_vals[i] = Rap1a
Gefa_vals[i] = GEFa
if both: return time, Rap1a_vals, Gefa_vals
return time, Rap1a_vals
適当なパラメタでシミュレーション¶
In [4]:
rap_params=[0.5,0.5,0.2,10]
In [5]:
time,Rap1_vals = Rap1model(params=rap_params)
In [6]:
plt.plot(time,Rap1_vals)
plt.title("Simple Rap1 model simulation")
plt.xlabel("time")
plt.ylabel("conc.")
plt.show()
Rasモデル¶
$$ \begin{aligned} \frac{d\mathrm{S}}{dt} &= 0\\ \frac{d\mathrm{GEF_a}}{dt} &= k_1\left[\mathrm{S}\right]\left[\mathrm{GEF_i}\right] - k_2\left[\mathrm{GEF_a}\right]\\ &=k_1\left[\mathrm{S}\right]\left(1-\left[\mathrm{GEF_a}\right]\right) - k_2\left[\mathrm{GEF_a}\right]\\ \frac{d\mathrm{GAP_a}}{dt} &= k_3\left[\mathrm{S}\right]\left[\mathrm{GAP_i}\right] - k_4\left[\mathrm{GAP_a}\right]\\ &= k_3\left[\mathrm{S}\right]\left(1-\left[\mathrm{GAP_a}\right]\right) - k_4\left[\mathrm{GAP_a}\right]\\ \frac{d\mathrm{Ras_a}}{dt} &= k_5\left[\mathrm{GEF_a}\right]\left[\mathrm{Ras_i}\right] - k_6\left[\mathrm{GAP_a}\right]\left[\mathrm{Ras_a}\right]\\ &= k_5\left[\mathrm{GEF_a}\right]\left(1-\left[\mathrm{Ras_a}\right]\right) - k_6\left[\mathrm{GAP_a}\right]\left[\mathrm{Ras_a}\right] \end{aligned} $$モデルを関数化¶
In [7]:
# Ras model の微分方程式
def RasODE(S,GEFa,GAPa,Rasa,k1,k2,k3,k4,k5,k6,dt):
dS = 0 * dt
dGEFa = (k1*S*(1-GEFa) - k2*GEFa) * dt
dGAPa = (k3*S*(1-GAPa) - k4*GAPa) * dt
dRasa = (k5*GEFa*(1-Rasa) - k6*GAPa*Rasa) * dt
return S+dS,GEFa+dGEFa,GAPa+dGAPa,Rasa+dRasa
In [8]:
def Rasmodel(params,dt=0.01,min_t=0,max_t=100,
S=1,GEFa=0,GAPa=0,Rasa=0,both=False):
max_t+=dt
k1,k2,k3,k4,k5,k6=params
time = np.arange(min_t,max_t,dt)
Rasa_vals = np.zeros(len(time))
GEFa_vals = np.zeros(len(time))
for i,t in enumerate(time):
S,GEFa,GAPa,Rasa = RasODE(S,GEFa,GAPa,Rasa,k1,k2,k3,k4,k5,k6,dt)
Rasa_vals[i] = Rasa
GEFa_vals[i] = GEFa
if both: return time, Rasa_vals, GEFa_vals
return time, Rasa_vals
適当なパラメタでシミュレーション¶
In [9]:
ras_params = [0.5,5,0.0005,0.005,0.05,100]
In [10]:
time,Rasa_vals = Rasmodel(ras_params)
In [11]:
plt.plot(time, Rasa_vals)
plt.title("Simple Ras model simulation")
plt.xlabel("time")
plt.ylabel("conc.")
plt.show()
課題2：1つ or 2つの未知パラメータを推定¶
- 最小二乗法に基づき、$\mathrm{Rap1}$ モデルのパラメータ $k_4$ を手で推定せよ。（省略）
- 最小二乗法に基づき、$\mathrm{Rap1}$ モデルのパラメータ $k_4$ を自動で推定せよ。
- 最小二乗法に基づき、$\mathrm{Rap1}$ モデルのパラメータ $k_2,k_4$ を自動で推定せよ。
課題2.2¶
$k4$ の候補を総なめし、最も値が良かったものを採用する、という方法をとる。
データの読み込み¶
In [12]:
train_data_kadai2_2 = loadmat("data_Rap1_1para.mat")
In [13]:
train_t_kadai2_2 = train_data_kadai2_2["t"].ravel()
train_x_kadai2_2 = train_data_kadai2_2["x"].ravel()
評価関数（RSS）の定義¶
In [14]:
def CalRap1RSS(train_t, train_x, params,
dt=0.1, min_t=0, max_t=100,
S=1,GEFa=0,GAPa=0.005,Rap1a=0):
max_t+=dt
k1,k2,k3,k4=params
time = np.arange(min_t,max_t,dt)
timing = train_t/dt
rss=0
for i,t in enumerate(time):
S,GEFa,GAPa,Rap1a = Rap1ODE(S,GEFa,GAPa,Rap1a,k1,k2,k3,k4,dt)
if GEFa>100: return 1e10
if i in timing:
rss += (train_x[np.argmax(timing==i)]-Rap1a)**2
return rss
パラメタを変化させて評価関数を最適化¶
In [15]:
k1=0.5
k2=0.5
k3=0.2
In [16]:
k4range = np.arange(0.1,10.05,0.05)
In [17]:
params_kadai2_2 = [[k1,k2,k3,k4] for k4 in k4range]
In [18]:
RSS_kadai2_2 = [
CalRap1RSS(
train_t_kadai2_2,
train_x_kadai2_2,
param_kadai2_2,
dt=1
) for param_kadai2_2 in params_kadai2_2
]
In [19]:
opt_id_kadai2_2 = np.argmin(RSS_kadai2_2)
opt_k4_kadai2_2 = k4range[opt_id_kadai2_2]
opt_RSS_kadai2_2 = RSS_kadai2_2[opt_id_kadai2_2]
In [20]:
plt.plot(k4range, RSS_kadai2_2)
plt.scatter(opt_k4_kadai2_2, opt_RSS_kadai2_2, label="optimized k4", s=300, marker="*", color="red")
plt.text(opt_k4_kadai2_2+1, opt_RSS_kadai2_2, "optimized k4={:.2f}".format(opt_k4_kadai2_2))
plt.title("The relationship between $k_4$ and RSS")
plt.xlabel("$k_4$")
plt.ylabel("RSS")
plt.legend()
plt.show()
最適なパラメタでシミュレーション¶
In [21]:
min_param_kadai2_2 = params_kadai2_2[opt_id_kadai2_2]
In [22]:
time,Rap1_vals = Rap1model(params=min_param_kadai2_2)
In [23]:
plt.plot(time, Rap1_vals)
plt.scatter(train_t_kadai2_2, train_x_kadai2_2, label="data")
plt.title("The best $k_4$ with sample data")
plt.xlabel("time")
plt.ylabel("conc.")
plt.legend()
plt.show()
課題2.3¶
$k2, k4$ の候補を総なめし、最も値が良かったものを採用する、という方法をとる。なお、値が収束するまで繰り返すことはしない。
データの読み込み¶
In [24]:
train_data_kadai2_3 = loadmat("data_Rap1_2para.mat")
In [25]:
train_t_kadai2_3 = train_data_kadai2_3["t"].ravel()
train_x_kadai2_3 = train_data_kadai2_3["x"].ravel()
パラメタを変化させて評価関数を最適化¶
In [26]:
k1 = 0.5
k3 = 0.2
In [27]:
k2range = np.arange(0.1,3.05,0.05)
k4range = np.arange(0.1,3.05,0.05)
In [28]:
params_kadai2_3 = [[k1,k2,k3,k4] for k4 in k4range for k2 in k2range]
In [29]:
RSS_kadai2_3 = [
CalRap1RSS(
train_t_kadai2_3,
train_x_kadai2_3,
param_kadai2_3,
dt=0.1
) for param_kadai2_3 in params_kadai2_3
]
In [30]:
RSS_kadai2_3 = np.array(RSS_kadai2_3).reshape(len(k2range),len(k4range))
In [31]:
opt_k2_index, opt_k4_index = np.unravel_index(np.argmin(RSS_kadai2_3), RSS_kadai2_3.shape)
optk2 = k2range[opt_k2_index]
optk4 = k4range[opt_k4_index]
optRSS = RSS_kadai2_3[opt_k2_index][opt_k4_index]
In [32]:
X, Y = np.meshgrid(k2range, k4range)
In [33]:
fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("$k_2$")
ax.set_ylabel("$k_4$")
ax.set_zlabel("RSS")
ax.plot_wireframe(X, Y, RSS_kadai2_3)
ax.scatter3D(optk2,optk4,optRSS, s=300, marker="*", color="red")
ax.set_title("The relationship between RSS and $k_2$ and $k_4$")
ax.text3D(optk2,optk4-1,optRSS,"optimized k2={:.2f}, k4={:.2f}".format(optk2,optk4))
plt.show()
最適なパラメタでシミュレーション¶
In [34]:
min_param_kadai2_3 = params_kadai2_3[np.argmin(RSS_kadai2_3)]
In [35]:
time,Rap1_vals = Rap1model(params=min_param_kadai2_3, dt=0.1)
In [36]:
plt.plot(time, Rap1_vals)
plt.scatter(train_t_kadai2_2, train_x_kadai2_2, label="data")
plt.title("The best k4 with sample data")
plt.xlabel("time")
plt.ylabel("conc.")
plt.legend()
plt.show()
課題3：4つの未知パラメータをEPで推定¶
- 最小二乗法に基づき、$\mathrm{Rap1}$ モデルのパラメータ $k_1,k_2,k_3,k_4$ を手で推定せよ。（省略）
- 最小二乗法に基づき、$\mathrm{Rap1}$ モデルのパラメータ $k_4$ をEPで推定せよ。
EP（Evolutionary programming）¶
- 乱数を初期値として，$N$ 個の個体（解の候補）を生成する。
- 個体のそれぞれのコピーをつくる。
- 2.で作った各コピーに正規乱数を加える。（変異）
- 3.の操作で作られた新しい $N$ 個の個体と、元の個体を混ぜた $2N$ 個の各個体に対してスコアを求める。スコアは以下のように決める。
- a）スコアを決める対象の個体を除いて、2N個体の集団から無作為に $q$ 個体を選ぶ。
- b）選んだ $q$ 個体のうち，スコアを決める対象の個体より評価関数の値（RSSなど）が悪い個体（最小二乗法の場合は，対象の個体よりRSSの値が大きい個体）の数をその個体のスコアとする。
- スコアの下位 $N$ 個の個体を削除する。（淘汰）
- 終了条件を満たすまで 2~5（1世代に相当）の操作を繰り返す。
- 残った中で最も評価関数の値が良い個体を"解"とする。
データの読み込み¶
In [37]:
data = loadmat("data_Rap1_4para.mat")
ts = data["t"]
xs = data["x"]
_, train_GEFa, _, train_Rap1a = xs.T
In [38]:
train_data_kadai3 = loadmat("data_Rap1_4para.mat")
In [39]:
train_t_kadai3 = train_data_kadai3["t"].ravel()
train_x_kadai3 = train_data_kadai3["x"]
In [40]:
print("訓練データの中身\n{}".format(train_x_kadai3))
In [41]:
_, train_GEFa, _, train_Rap1a = train_x_kadai3.T
In [42]:
print("GEFa:\n{}".format(train_GEFa))
print("Rap1a:\n{}".format(train_Rap1a))
諸々の関数定義¶
In [43]:
# 摂動(dt)が大きいと、パラメータの値によっては急激な変化によってうまくシミュレーションができないことがある。
# そこで、濃度があまりにも大きい（or負の数になった）時には大きなペナルティを与える。
# そうすることでdtを比較的大きな値に設定することができるので、短時間で学習が行える。
def AvoidOverFlow(values,minv=0,maxv=1e3):
return sum([v<minv or maxv<v for v in values])
In [44]:
def CalRap1RSS_both(train_t, train_GEFa, train_Rap1a, params,
dt=0.01, min_t=0, max_t=100, maxrss=1e5,
S=1,GEFa=0,GAPa=0.005,Rap1a=0):
max_t+=dt
k1,k2,k3,k4=params
time = np.arange(min_t,max_t,dt)
timing = train_t/dt
rss=0
for i,t in enumerate(time):
S,GEFa,GAPa,Rap1a = Rap1ODE(S,GEFa,GAPa,Rap1a,k1,k2,k3,k4,dt)
if AvoidOverFlow([S,GEFa,GAPa,Rap1a]): return maxrss
if i in timing:
index = np.argmax(timing==i)
rss += (train_Rap1a[index]-Rap1a)**2 + (train_GEFa[index]-GEFa)**2
return rss
In [45]:
def plotSimulate(ax1,ax2,param):
time,Rap1a_vals,Gefa_vals = Rap1model(param,both=True)
ax1.plot(time, Gefa_vals, color="red", alpha=0.5)
ax2.plot(time, Rap1a_vals, color="red", alpha=0.5)
return ax1,ax2
In [46]:
def plotRSS(ax3,ax4,generations,best_rss,mean_rss):
ax3.plot(generations, best_rss, color="blue", alpha=0.5)
ax4.plot(generations, mean_rss, color="blue", alpha=0.5)
return ax3,ax4
In [47]:
def plotInfo(train_t, train_GEFa, train_Rap1a, figsize=(12,8)):
fig = plt.figure(figsize=figsize)
ax1=fig.add_subplot(2,2,1)
ax1.set_xlabel("time")
ax1.set_title("GEFa")
ax1.scatter(train_t, train_GEFa, label="train data")
ax1.legend()
ax2 = fig.add_subplot(2,2,2)
ax2.set_xlabel("time")
ax2.set_title("Rap1")
ax2.scatter(train_t, train_Rap1a, label="train data")
ax2.legend()
ax3 = fig.add_subplot(2,2,3)
ax3.set_xlabel("generation")
ax3.set_ylabel("RSS")
ax3.set_title("The Best individual")
ax3.set_ylim(0,1.5)
ax4 = fig.add_subplot(2,2,4)
ax4.set_xlabel("generation")
ax4.set_ylabel("RSS")
ax4.set_title("The Mean of individuals")
ax4.set_ylim(0,6)
return ax1,ax2,ax3,ax4
In [48]:
import datetime
def now():
return datetime.datetime.now().strftime("%Y/%m/%d %H:%M:%S")
パラメタ定義¶
In [49]:
np.random.seed(2019) # シード固定(→再現性)
In [50]:
numGeneration = 250
ITERTION_COUNT = 10
In [51]:
K=4 # パラメタ数
N=100 # 一世代の個体数
indexes = np.array([i for i in range(N*2)]) # 各個体のインデックス
In [52]:
lb = 1e-3; ub = 1e2 # 初期パラメータの範囲
sigma_scale=0.1; # 正規乱数の分散パラメータ
In [53]:
rate=0.7 # スコア算出時に何割の個体と値を比べるか
q = int(N*2*rate)
EP実装¶
In [54]:
print("[START] {}".format(now()))
print("[PARAM] Iteration={}, numGeneration={}, numIndividual={}".format(ITERTION_COUNT, numGeneration,N))
ax1,ax2,ax3,ax4 = plotInfo(train_t_kadai3, train_GEFa, train_Rap1a)
for it in range(ITERTION_COUNT):
print("="*60)
params = np.random.uniform(lb,ub,(N,K)) # 初期パラメータ
best_rss = np.zeros(numGeneration+1) # 各世代の最も良い個体の値を格納
mean_rss = np.zeros(numGeneration+1) # 各世代の平均値を格納
generations = [i for i in range(numGeneration+1)]
for gen in range(numGeneration+1):
# NOTE:以下の値で新しいパラメタを設定する。
# 平均: 元のパラメタ
# 分散: sigma_scale × 元のパラメータ
new_params = np.random.normal(loc=params, scale=params*sigma_scale)
#=== パラメタを[lb,ub]の範囲に入れる ===
new_params = np.where(new_params>lb, new_params, lb)
new_params = np.where(new_params<ub, new_params, ub)
all_params = np.r_[params,new_params]
RSSs = np.array([CalRap1RSS_both(train_t_kadai3, train_GEFa, train_Rap1a, param, dt=0.1) for param in all_params])
best_rss[gen] = np.min(RSSs)
mean_rss[gen] = np.mean(RSSs)
scores = np.zeros(2*N)
for i in range(2*N):
rss = RSSs[i] # 今考えている個体のRSS
selected = np.random.choice(indexes[indexes!=i], q) # 考えている個体以外からq個体をランダムチョイス
selected_scores = RSSs[selected] # それらの個体のRSS
score = np.count_nonzero(selected_scores > rss) # 考えている個体よりもRSSが悪い（=大きい）個体の数
scores[i] = score # 記録
params = all_params[np.argsort(scores)[N:]] # スコアの良い半分を次の世代に残す。
best_params = all_params[np.argmax(scores)]
if gen%10==0:
k1,k2,k3,k4 = best_params
print("[{:02} | Gene{:03}] {}".format(it,gen,now()))
print("- best RSS: {:.5f} / k1={:.2f}, k2={:.2f}, k3={:.2f}, k4={:.2f}".format(best_rss[gen],k1,k2,k3,k4))
print("- mean RSS: {:.5f}".format(mean_rss[gen]))
print("-"*60)
ax1,ax2 = plotSimulate(ax1,ax2,best_params)
ax3,ax4 = plotRSS(ax3,ax4,generations,best_rss,mean_rss)
plt.tight_layout()
plt.show()
発展課題：Rap1のデータに対して、Rasモデルのパラメータ（$k_1,\ldots,k_6$）を推定せよ。¶
- RasモデルとRap1モデルの関係はどうなっているか？
- Rasモデルを使った場合とRap1モデルを使った場合で、パラメータや残差二乗和にどのような違いが現れたか？
- 残差二乗和が似ている最適化試行では、パラメータの推定値も似ているか？
予想：¶
$\frac{dGAP}{dt}=0$ になるようにパラメータがフィッティング（つまり、Rasモデルと本質的には変わらないということ）するのではないか？
もしくは変にオーバーフィッティング？？
In [ ]: