2010-12-10

中庸

こちらで、セルオートマトンを扱った
決定論的な推移が観察できて、ある形が、そのままの形で２次元平面を移動するパターンを作った
これはライフゲームと呼ばれるものである(Wiki記事)
Wiki記事にもあるように、周辺の濃度過多と濃度過少とが死をもたらし、そうでないときに、生存する。また、死の条件があるので、個体数が減りすぎないように誕生の条件がある、と読める。
- 誕生:死んでいるセルに隣接する生きたセルがちょうど3つあれば、次の世代が誕生する。
- 生存:生きているセルに隣接する生きたセルが2つか3つならば、次の世代でも生存する。
- 過疎:生きているセルに隣接する生きたセルが1つ以下ならば、過疎により死滅する。
- 過密:生きているセルに隣接する生きたセルが4つ以上ならば、過密により死滅する。
多型の生死を考えよう
多型箇所が残る「理由」として、ヘテロが有利な場合、というのがあるのは、良く知られたこと
- SNPなどの２アレル型多型で言えば、あるアレルのホモと別のアレルのホモよりも、ヘテロの方が有利なとき、多型性は消失しにくい
これを多座位の複合形質に適用すると、多座位のもたらす効果がある方向に偏っているときと、別の方向に偏っているときに、不利で、そうでない場合が有利、と
このアナロジーで行けば、ライフゲームを繁殖・生殖に使おうとすれば、「過密・過疎」を、「１染色体上の機能性アレルの１方向への蓄積ともう１方向への蓄積」に対応させて、それを「死」の条件とするなどすることが可能だろうか
また、ここ数日のシリーズ記事にあるように、空間に波があって、それが伝播するようにしたいので、なにがしかの条件を入れないといけないのだが、まだ、どうすると「定形波」の進行になるのか、その生物的意味は何かが、ぼんやりしているので、そこは、置いておくとして、データをうまく図示できるようにしておくことは、常に、ポジティブな役割を果たすので、そんなことをやってみる
なお、多型が、存続しがち、ということを取り扱うときに、ドリフトアウトが多発すると面倒なので、空間を、閉じる(１次元空間なら、円にする)ことで、その点の手間を省くことも、有効な気がするので、そんな風にしてみたい
ハプロイドで扱っているので、(当然のことながら)HWEを仮定している。ディプロイドに持ち上げるには、まだ、色々、ある。こちらでコメントしたように、ハプロイドに空間上の位置を入れると、HWEを仮定してもHWDを生じさせることはできるが…。
色々ありすぎて、ライフゲームやソリトンのモデルが使えなくなると、面白くないから工夫が必要か…

Nx<-200 # 空間　１次元格子
Nm<-100 # 座位数
Nt<-500 # 世代数

# 染色体のアレル格納
A<-array(0,c(Nx,Nm,Nt))

# 初期値発生
Xrange<-round(Nx*0.45,0):round(Nx*0.55,0)
Xrange<-1:Nx
A[Xrange,,1]<-sample(c(0,1),length(Xrange)*Nm,replace=TRUE,prob=c(0.95,0.05))

Xrange1<-1:round(Nx*0.1,0)
Xrange2<-round(Nx*0.55,0):round(Nx*0.6,0)
A[Xrange1,,1]<-sample(c(0,1),length(Xrange1)*Nm,replace=TRUE,prob=c(0.1,0.9))
A[Xrange2,,1]<-sample(c(0,1),length(Xrange2)*Nm,replace=TRUE,prob=c(0.3,0.7))

# メイティング条件(近すぎるところにある染色体同士は組合わせない(禁　近親婚)
d<-3

# 極端なアレル集積を認めない条件
fUp<-1
fLo<-0

UpL<-Nm*fUp
LoL<-Nm*fLo

# 極端なアレル集積ハプロイドができたら、それは生まれずに致死として扱い、次の子をなすことにするときの、最大試行回数
MaxLoop<-10

# ハプロイドごとの、「機能性アレル数」の集計用
SumM<-matrix(0,Nx,Nt)
SumM[,1]<-apply(A[,,1],1,sum)
for(i in 2:Nt){
	for(j in 1:Nx){
		loop<-TRUE
		countLoop<-0
		while(loop){
			# 親ハプロイドは周辺から
			tmpP<-(j-d)%%Nx
			tmpM<-(j+d)%%Nx
			# ハプロイド空間を円周にする条件
			if(tmpP==0)tmpP<-Nx
			if(tmpM==0)tmpM<-Nx
			#if(tmpP<1)tmpP<-Nx+tmpP
			#if(tmpM>Nx)tmpM<-tmpP-Nx
			# 全座位は組換え自由(異なる染色体上)の仮定
			tmpDip<-cbind(A[tmpP,,i-1],A[tmpM,,i-1])
			selection<-sample(c(1,2),Nm,replace=TRUE)
			tmp<-tmpDip[cbind(1:Nm,selection)]
			# 極端ハプロイドの排除
			sumtmp<-sum(tmp)
			if(countLoop==MaxLoop){
				A[j,,i]<-tmp
				loop<-FALSE
			}else if(sumtmp>=LoL & sumtmp<=UpL){
				A[j,,i]<-tmp
				loop<-FALSE
			}
			countLoop<-countLoop+1
		}

	}
	SumM[,i]<-apply(A[,,i],1,sum)
}

image(SumM)

persp(SumM,phi=60,theta=50,xlab="chromosomes",ylab="generations",zlab="No.alleles")

tmpxyz<-cbind(c(row(SumM)),c(col(SumM)),c(SumM))

library(rgl)

geta<-1

xx<-(geta+sqrt(tmpxyz[,3]))*cos(tmpxyz[,1]/Nx*2*pi)
yy<-(geta+sqrt(tmpxyz[,3]))*sin(tmpxyz[,1]/Nx*2*pi)
zz<-tmpxyz[,2]

maxxy<-max(abs(xx),abs(yy))

xylim<-c(-maxxy,maxxy)

plot3d(xx,yy,zz,xlim=xylim,ylim=xylim,type="l",xlab="SpaceX",ylab="SpaceY",zlab="Generation",main="Circular Space x time",col="blue")


plot3d(row(SumM),col(SumM),SumM,type="l")
plot3d(row(SumM),col(SumM),SumM,xlab="Space",ylab="Generation",zlab="Value")