講義全体の流れ¶
ノートに記録するのは初めてなので、今回の実習の講義日程も記録しておく。
目的
次世代シークエンス解析関連の実験操作と解析を学ぶ
| # | 日付 | 内容 |
|---|---|---|
| 1 | 6/25(火) | はじめに/Total RNAの調製とBioanalyzerによる濃度測定 |
| 2 | 6/26(水) | ナノポアシークエンサーMinION RNA-seq 1日目 |
| 3 | 6/27(木) | ナノポアシークエンサーMinION RNA-seq 2日目 |
| 4 | 6/28(金) | イルミナシークエンサー RNA-seq 1日目 |
| 5 | 7/2(火) | イルミナシークエンサー RNA-seq 2日目 |
| 6 | 7/3(水) | 数理モデル講義・演習（社会連携講座） |
| 7 | 7/4(木) | RNA-seq解析 |
| 8 | 7/5(金) | RNA-seq解析とまとめ |
解析解¶
$$ \begin{aligned} \frac{dx(t)}{dt} &= r\left(1-\frac{x(t)}{K}\right)x(t)\\ \frac{dx(t)}{dt} &= \frac{r}{K}\left(K-x(t)\right)x(t)\\ \int\left(\frac{1}{K-x(t)}+\frac{1}{x(t)}\right)dx(t) &= \int rdt\\ \ln\left(\frac{x(t)}{K-x(t)}\right) &= Crt\qquad(C:\mathrm{const.})\\ \frac{x(t)}{K-x(t)} &= Ae^{rt}\qquad(A=e^C)\\ \end{aligned} $$ここで、$t=0$ で $x(t)=x_0$ だから、$A = \frac{x_0}{K-x_0}$ となるので、これを代入して、
$$x(t) = \frac{K}{1 + \left(\frac{K}{x_0}-1\right)e^{-rt}}$$シミュレーション（数値計算）¶
In [1]:
import numpy as np
import matplotlib.pyplot as plt
In [2]:
xmin=0; xmax=10
K = 1000
r = 1.2
x0 = 10
In [3]:
# 厳密解
def Logistic_equation(t,r,K,x0):
return K/(1+ (K/x0-1)*(np.e**(-r*t)))
t_exact = np.linspace(xmin, xmax, 1000)
exact = np.vectorize(Logistic_equation)(t_exact,r,K,x0)
In [4]:
# 数値計算
def d_Logisitc(x,r,K):
return r*(1-x/K)*x
t_numerical = [xmin+i for i in range(xmax-xmin+1)]
numerical = []
x = x0
for i in range(xmax+1):
x = x + d_Logisitc(x,r,K)
numerical.append(x)
In [5]:
plt.scatter(t_numerical, numerical, label="numerical")
plt.plot(t_exact, exact, label="exact")
plt.title("Logistic equation: $r={},\ K={}, x_0={}$".format(r,K,x0))
plt.xlabel("$t$")
plt.ylabel("$x(t)$")
plt.legend()
plt.show()
組織の中で1細胞生まれたがん細胞が拡がる確率を計算する：¶
単位時間経過後にがん細胞の数が $i$ から $j$ になる確率を $P_{ij}$ とすると、
$$ \begin{cases} P_{i\rightarrow i+1} &= \frac{N-i}{N}\frac{ri}{(N-i) + ri} \\ P_{i\rightarrow i-1} &= \frac{i}{N}\frac{N-i}{(N-i) + ri} \\ P_{i\rightarrow i} &= 1-P_{i\rightarrow i+1}-P_{i\rightarrow i+1} = \frac{1}{N}\frac{(N-i)^2+ri^2}{(N-i) + ri}\\ \end{cases} $$したがって、$i$ 細胞からスタートしたがん細胞が最終的に固定する確率を $V_i$ とすると、
$$ \begin{aligned} V_i &= V_{i+1}P_{i\rightarrow i+1} + V_{i-1}P_{i\rightarrow i-1} + V_iP_{i\rightarrow i}\\ &= V_{i+1}\frac{N-i}{N}\frac{ri}{(N-i) + ri} + V_{i-1}\frac{i}{N}\frac{N-i}{(N-i) + ri} + V_i\frac{1}{N}\frac{(N-i)^2+ri^2}{(N-i) + ri}\\ V_{i+1}-V_i &= \frac{1}{r}(V_i-V_{i-1}), \quad (V_0=0, V_N=1)\\ V_i &= V_1\frac{1-(1/r)^i}{1-1/r} \end{aligned} $$$i=N$ を代入して、
$$V_1 = \frac{1-1/r}{1-(1/r)^N}$$シミュレーション（数値計算）¶
In [6]:
rs = [i/10 for i in range(5,21)]
Ns = [10,100]
In [7]:
# 厳密解
def Fixation_probability(N,r):
if r!=1: return (1-1/r) / (1-(1/r)**N)
else: return 1/N
exact10,exact100 = [[Fixation_probability(n,r) for r in rs] for n in Ns]
In [8]:
ITER_MAX = 10000
COUNT = 1000
# シミュレーション
def Numeric_analysis(N,r,seed=0):
np.random.seed(seed)
fixed = 0
for c in range(COUNT):
mu=1 # mutation
no=N-1 # normal
for it in range(ITER_MAX):
rand = np.random.rand()
if (N-mu)/N * r*mu/(N-mu + r*mu) > rand: mu+=1; no-=1
elif mu/N * (N-mu)/(N-mu+r*mu) > (1-rand): mu-=1; no+=1
else: continue
if no == N: break
if mu == N: fixed += 1; break
return fixed/COUNT
numerical10, numerical100 = [[Numeric_analysis(n,r) for r in rs] for n in Ns]
In [9]:
plt.scatter(rs, exact10, c="#80160e", marker="s", label="exact, $N=10$"); plt.plot(rs, exact10,c="#80160e")
plt.scatter(rs, exact100, c="#80160e", marker="o", label="exact, $N=100$"); plt.plot(rs, exact100,c="#80160e")
plt.scatter(rs, numerical10, c="#6da3fc", marker="s", label="numerical, $N=10$"); plt.plot(rs, numerical10,c="#6da3fc")
plt.scatter(rs, numerical100, c="#6da3fc", marker="o", label="numerical, $N=100$"); plt.plot(rs, numerical100,c="#6da3fc")
plt.title("The probability of the fixation probability $V_1$")
plt.xlabel("fitness $r$")
plt.ylabel("$V_1$")
plt.grid()
plt.legend()
plt.show()
$N$ が小さければ、適応度 $r$ が $1$ より小さくても固定する可能性がある！
ということが言いたかったらしいのですが、当たり前やなあ。という感じですね。
課題3: 競争系のダイナミクス¶
$$ \begin{cases} \frac{dx}{dt} = (a-bx-cy)x\\ \frac{dy}{dt} = (d-ex-fy)y \end{cases}\ \cdots (1) $$【仮定】
- 種 $x,y$ はそれぞれ $x,y$ がいなければロジスティック成長（$b,e$ が種内競争係数を意味。）
- 種 $x$ と種 $y$ はお互いに増殖を抑制し合う（$c,f$ が種間競争係数を意味。）
1. 式$(1)$ の平衡点を解く¶
解の平衡点は、$\frac{dx}{dt}=0$ かつ $\frac{dy}{dt}=0$ を満たす点なので、以下の $4$ 点である。
$$ \begin{cases} x = 0 & y = 0\\ x = 0 & y = \frac{d}{f}\\ x = \frac{a}{b} & y = 0\\ x = \frac{af-cd}{bf-ce} & y = \frac{ae-bd}{ce-bf} \end{cases} $$2. 種 $x$ と種 $y$ が共に正の値をとる平衡点が出現するパラメータ $(a,b,c,d,e,f)$ の関係式を示す¶
以下の平衡点 $(x,y)$ が共に正の値を取れば良い。
$$x = \frac{af-cd}{bf-ce}, y = \frac{bd-ae}{bf-ce}$$3. アイソクライン法で問2の平衡点の安定性を調べる¶
In [10]:
ITER_MAX = 1000
def plot_nurcline(a,b,c,d,e,f,inix,iniy):
plt.clf()
plt.figure(figsize=(6,6))
xmax=int(a/b)+1
ymax=int(d/f)+1
x=np.linspace(0,xmax,10)
y=np.linspace(0,ymax,10)
# ヌルクライン
X,Y = np.meshgrid(x, y)
dX = (a - b*X - c*Y)*X
dY = (d - e*X - f*Y)*Y
plt.quiver(X, Y, dX, dY)
# 直線
func1 = a/c - b/c*x
func2 = d/f - e/f*x
plt.plot(x,func1, color="#721940", label="line1")
plt.plot(x,func2, color="#F1E3A7", label="line2")
# 平衡点
plt.scatter(0,0, s=500, color="#1A1406", label="bp1")
plt.scatter(0,d/f, s=500, color="#721940", label="bp2")
plt.scatter(a/b,0, s=500, color="#F1E3A7", label="bp3")
plt.scatter((a*f-c*d)/(b*f-c*e), (a*e-b*d)/(c*e-b*f), s=500, color="#FFF2DE", label="bp4")
# シミュレーション
x=inix; y=iniy
lst_x=[]; lst_y=[]
for it in range(ITER_MAX):
lst_x.append(x); lst_y.append(y);
x+=1/10*(a-b*x-c*y)*x; y+=1/10*(d-e*x-f*y)*y
plt.plot(lst_x,lst_y, color="black", label="simulation: {}".format(i))
plt.scatter(inix,iniy, s=50, color="black", label="initial state: {}".format(i), marker="s")
plt.title("The simulation results.")
plt.xlim(0,xmax)
plt.ylim(0,ymax)
plt.legend()
plt.grid()
plt.show()
In [11]:
plot_nurcline(a=1,b=1,c=2,d=1,e=2,f=1,inix=0.8,iniy=0.3) # af<cd, bd<ae
In [12]:
plot_nurcline(a=1,b=3,c=2,d=1,e=2,f=1,inix=0.8,iniy=0.3) # af<cd, bd>ae
In [13]:
plot_nurcline(a=1,b=1,c=2,d=1,e=2,f=3,inix=0.8,iniy=0.3) # af>cd, bd<ae
In [14]:
plot_nurcline(a=1,b=3,c=2,d=1,e=2,f=3,inix=0.8,iniy=0.3) # af>cd, bd>ae
