Moran 過程の数値シ ミ ュ レ ーシ ョ ン
2011SE195

夏川将実

指導教員：小藤俊幸

1はじ めに有限集団での進化ダイナミクスの解析を始めよう．個体数は連続変数ではなく，整数で与えられる．その結果，進化ダイナミクスはもはや決定論的な微分方程式では記述さ(iv) 繁殖に B の 1 個体が，死亡に A の 1 個体が選択され，この事象が起こった後，A は 1 個体減り，変数 i はi − 1 となる．この事象の確率は i(N − i)/N 2 である．

3

遷移確率
xni を時刻 n で i 個である確率とし,れず，確率論的な定式化が必要である．
生物学的な問題を研究するための最良のアプローチは，i(N − i)
N2i(N − i) i(N − i)
Pi,i = 1 −−N2
N2i(N − i)Pi,i+1 =N2Pi,i−1 =まず決定論的な記述を試み，それが適切な側面を損なっているときのみ，確率論的な解析へ移行することである．この研究では有限サイズの集団における中立的浮動と一定の淘汰について調べ，ダーウィンの進化論に相対する中立説を指示していきたい．

2

中立的浮動
サイズ N で固定されたある集団を考えよう．個体の種

(1)
(2)
(3)

とおく. 時刻 n で i − 1 個のとき, 時刻 n + 1 で i 個となる
確率は Pi,i−1 xn
i−1 , 時刻 n で i 個のとき, 時刻 n + 1 で i 個となる確率は Pi,i xn
i , 時刻 n で i + 1 個のとき, 時刻 n + 1n
類は A と B の 2 種類ある．個体は同じ速度で繁殖する． で i 個となる確率は Pi+1,i xi+1 であるから,したがって，A と B は選択に関して中立な変異体である．

xn+1= Pi−1,i xni−1 + Pi,i xni + Pi+1,i xni+1i

すべての単位時間ステップにおいて，１個体がランダムに選択されて繁殖し，かつ１個体がランダムに選択されて除去される．ここでは復元抽出を用いよう．すなわち，同一の個体が繁殖と死亡に選ばれることも起こり得る．繁殖は突然変異することなく起こる．すなわち A は A を，B はB を生産する．

(4)となる.

4

プロ グラ ムMoran 過程をシミュレートするプログラムを作成した.
1 個体を繁殖に選択する

この確率過程はオーストラリアの集団遺伝子学者

P.A.P.Moran にちなんで名付けられており，1958 年に考案された．各時間ステップで常に 1 個体の出生と 1 個体の死亡が起こるため，集団サイズは必ず一定にあることが保証されてる．確率変数は個体 A の個体数のみであり，ここでは i と表そう．個体 B の個体数は N − i である．1変数の確率過程を調べることは，2 進数あるいはそれ以上の変数の過程よりも非常に簡単である．

Moran 過程は状態空間 i = 0, ...., N で定義される．個体 A を (出生または死亡で) 選択する確率は i/N で与えられる．個体 B を選択する確率は (N − i)/N で与えられる．
すべての単位時間ステップにおいて，次の 4 つの事象が起そして 1 個体を死亡に選択するこる可能性がある.
(i) 繁殖と死亡に A の 1 個体が選択され，この事象が起 2 番目に選択された個体は 1 番目に選択された個体の子供に置き換えられる
こった後の個体 A の数は，起こる前の数と同じであり，i2は変化しない．この事象の確率は (i/N ) である．
プログラムでは集団 A，B を数字の 1，0 に置き換え，乱(ii) 繁殖と脂肪に B の 1 個体が選択され，この事象が起こった後の個体 B の数は，起こる前の数と同じであり， 数を発生させるため rand 関数を用いた．
N − i は変化しない．この事象の確率は [(N − i)/N ]2 で 4.1 プロ グ ラ ムある．

(iii) 繁殖に A の 1 個体が，死亡に B の 1 個体が選択され，この事象が起こった後，A は 1 個体増え，変数 i はi + 1 となる．この事象の確率は i(N − i)/N 2 である．

以下にソースリストを示す.
main(){

n=8 ●●●●○●●●○○
n=9 ●●●●○○●○○○
この結果では n=70 ですべてが●になった．

int x[N];int i, n, p, q;
srand48(10000);5for(i = 0; i < N/2; i++){x[i] = 0;}
for(i = N/2; i < N; i++){x[i] = 1;}

比較
最初の 1 の数が半分（5 個）としたとき，100 回目に 1が i 個現れる確率（変数 x[i]) を遷移確率から計算するプ
ログラムと，実際に乱数を使いシミュレーションした結果を比較する．遷移確率での結果はこのようになった．

for(n = 0; n < 100; n++){p = floor(N*drand48());
q = floor((N - 1)*drand48());
if (q >= p) q = q + 1;
printf("n = %d, p = %d, q = %d\n", n, p, q);for(i = 0; i < N; i++){printf("%d ", x[i]);}printf("\n");x[q] = x[p];}図 1 シミュレーション結果}

このプログラムは初期値として，0，1 を 5 個体ずつ用意し，乱数により上記の (i)(ii)(iii)(iv) の事象を起こさせど

上記のプログラムを発展させ，n = 99 の結果のみを 100のように変化していくかを表したものである．このプログ個表示するプログラムを作った．更に結果の 1 の数を数えラムを用いて発展させていく．まずはこのプログラムを動るように組み，遷移確率と比較できるようにようにした．

かしどのようになるのかを検証する．

0 と 1 の数は 1:9 とし,1 の個数はデータの平均をとった.

4.2

乱数について1 の数

drand48 は区間 [0, 1) 上の一様乱数 (疑似乱数) を生成個数する関数である. m = 248 = 281474976710656 を法とす

確率る乗数 (multiplicand) a = 0xfdeece66d = 25214903917

1 の数加数 (addend) c = 0xb = 11 の線形合同法 (linear con-07.670.07

100.01

20.330.01

30.330.0140.66
0.01

500.01600.01の7800.01

90.66
0.01

10
90.33
0.87

個数0gruential method) (1) rn+1 = arn + c (mod m) により確率0.01生成される乱数を m で割って,[0, 1) 上の乱数としている.
注 0x は 16 進数を表す. 最初の 0x を除いた部分が「数字」
で a, b, c, d, e, f は, それぞれ,10, 11, 12, 13, 14, 15 に対 6 終わり に応する.
遷移確率と中立的浮動のシミュレーションの結果はほぼ4.3

結果

4.1 のプログラムの実行結果を示す. 0 と 1 では表しにくいので，0 を●，1 を○とする．
n=0 ●●●●●○○○○○
n=1 ●●●●●○●○○○
n=2 ○●●●●○●○○○
n=3 ○●●●●●●○○○
n=4 ●●●●●●●○○○
n=5 ●●●●○●●○○○
n=6 ●●●●○●●○○○
n=7 ●●●●○●●●○○

一致するという結果になった。これより、遺伝子の変化は自然淘汰に対して有利でも不利でもなく中立であるという中立説を指示する.

7

参考文献

参考文献
[1] MARTIN A. NOWAK[著] 竹内康博・佐藤一憲・巌佐庸・中岡慎治 [監訳]：『進化のダイナミクス』．共立出版，東京，2008.