2015-02-21

多スイッチオンオフモデル〜遺伝子による確率論的決定機構

スイッチのオンオフはスイッチ数pのとき、 $2^p$ 個の頂点を持つp次元単位格子上の状態推移として考えられる
この酔歩にあって、ノードに存在確率の濃淡を作るには、あるノードには、流入しやすいけれど流出しにくい、ということを実現することが必要だ
何がそのようなオンオフの制御をするかは、後付けで考えるとして、このシミュレーションをするには
格子グラフを作る
流出路を絶つ
ということをすればよい
また、オン・オフのどちらがより可能性が高いか、ということに差をつければ、同じグラフにおいて、定常状態のノード確率ベクトルに違いが生じる
推移確率行列に固有ベクトルが得られるわけだけれど、そのどれに落ちるかは、初期値次第(または乱雑項の影響でドリフトするかも)
ロバスト性は、初期値が同じ・ほぼ同じときに、ある程度の乱雑項があっても同じ固有値・固有ベクトルに落ちてくれるようなことでもあるし
「低確率すぎるもの」には期待しない〜「０の状態」と「あるべき状態」との間で生起確率に峻別できる差を持たせる、というようなことでもある
実際には、「存在させる」という作業と、その状態の個数の確保とは別プロセスの方が無難だろう。後者は「増殖力」という仕組みで対応可能な感じがする
ほとんど起こりえないオンオフ状態の細胞がどうなるかも確率的には決まっている。どの状態にどれくらいの確率で収束するかによって、「状態間の距離」を定義できるだろう
以下のソースは、多次元単位格子上の酔歩

p <- 4
n <- 2^p
X <- as.matrix(expand.grid(rep(list(0:1),p)))
x <- apply(t(X) * 2^(0:(p-1)),2,sum)

D <- as.matrix(dist(X,"manhattan"))
d <- 0  + (D == 1)
n.edge <- sum(d)

edge.list <- which(d==1,arr.ind=TRUE)

edge.dir <- rep(0,length(n.edge))
edge.sign <- rep(0,length(n.edge))
edge.val <- rep(0,length(n.edge))
for(i in 1:n.edge){
	tmp <- X[edge.list[i,1],] - X[edge.list[i,2],]
	tmp2 <- which(abs(tmp)==1)
	tmp3 <- sign(tmp[tmp2])
	edge.val[i] <- tmp2
	edge.sign[i] <- tmp3
	edge.dir[i] <- tmp3 * tmp2
}

affected.item <- 1:n
edge.eff <- sample(affected.item,p)


subs <- matrix(0,n,n)
for(i in 1:n){
	subs[i,] <- apply(-t(X)+X[i,],2,min) >= 0
}
d2 <- d
for(i in 1:p){
	this.item <- affected.item[i]
	tmp <- edge.eff[i]
	tmp2 <- which(subs[,tmp]==1)
	affected.edge.id <- which(edge.val==this.item)
	affected.edge.list <- edge.list[affected.edge.id,]
	for(j in 1:length(affected.edge.list[,1])){
		if(affected.edge.list[j,1] %in% tmp2){
			d2[affected.edge.list[j,1],affected.edge.list[j,2]] <- 0
		}
	}
}

my.lattice <- function(p){
	n <- 2^p
	X <- as.matrix(expand.grid(rep(list(0:1),p)))
	x <- apply(t(X) * 2^(0:(p-1)),2,sum)

	D <- as.matrix(dist(X,"manhattan"))
	d <- 0  + (D == 1)
	return(list(X=X,d=d))
}
my.self.ctrl.lattice <- function(p,affected.item=1:2^p,edge.eff=sample(affected.item,p)){
	n <- 2^p
	tmp.lattice <- my.lattice(p)
	X <- tmp.lattice$X
	d <- tmp.lattice$d
	#X <- as.matrix(expand.grid(rep(list(0:1),p)))
	x <- apply(t(X) * 2^(0:(p-1)),2,sum)

	#D <- as.matrix(dist(X,"manhattan"))
	#d <- 0  + (D == 1)

	n.edge <- sum(d)

	edge.list <- which(d==1,arr.ind=TRUE)

	#edge.dir <- rep(0,length(n.edge))
	#edge.sign <- rep(0,length(n.edge))
	edge.val <- rep(0,length(n.edge))
	for(i in 1:n.edge){
		tmp <- X[edge.list[i,1],] - X[edge.list[i,2],]
		tmp2 <- which(abs(tmp)==1)
		tmp3 <- sign(tmp[tmp2])
		edge.val[i] <- tmp2
		#edge.sign[i] <- tmp3
		#edge.dir[i] <- tmp3 * tmp2
	}

	subs <- matrix(0,n,n)
	for(i in 1:n){
		subs[i,] <- apply(-t(X)+X[i,],2,min) <= 0
	}
	d2 <- d
	if(length(edge.eff)>0){
		for(i in 1:length(edge.eff)){
			this.item <- affected.item[i]
			tmp <- edge.eff[i]
			tmp2 <- which(subs[,tmp]==1)
			affected.edge.id <- which(edge.val==this.item)
			affected.edge.list <- edge.list[affected.edge.id,]
			for(j in 1:length(affected.edge.list[,1])){
				if(affected.edge.list[j,1] %in% tmp2){
					d2[affected.edge.list[j,1],affected.edge.list[j,2]] <- 0
				}
			}
		}
	}
	return(list(X=X,d.lat=d,d=d2,el=edge.list,e.val=edge.val,subs=subs,affected.item=affected.item,edge.eff=edge.eff))
}


my.rw.graph <- function(d,n.step,pr=1,x.init=c(1,rep(0,ncol(d)-1)),trans=FALSE){
	if(trans){
		d <- t(d)
	}
	n <- ncol(d)
	ret <- matrix(0,n.step+1,n)
	ret[1,] <- x.init
	tmp <- apply(d,1,sum)
	tmp2 <- tmp
	tmp2[which(tmp2==0)] <- 1
	d.pr <- d/tmp2*pr
	
	tmp.diag <- rep(1-pr,n)
	tmp.diag[which(tmp==0)] <- 1
	diag(d.pr) <- tmp.diag
	d.pr <- t(t(d.pr)/apply(d.pr,1,sum))
	for(i in 1:n.step){
		ret[i+1,] <- t(d.pr) %*% ret[i,]
	}
	ret
}

p <- 8
affected.item=1:2^p
k <- 0
edge.eff=sample(affected.item,p-k)
out <- my.self.ctrl.lattice(p,affected.item,edge.eff)
n.step <- 1000
pr <- 0.8
out2 <- my.rw.graph(out$d,n.step,pr)
par(mfcol=c(1,2))
matplot(out2,type="l")

plot(sort(out2[n.step+1,]),type="h")
par(mfcol=c(1,1))

p <- 8
d.lat <- my.lattice(p)$d
tmp <- d.lat[nrow(d.lat),]
tmp1 <- which(tmp==1)[1:4]
d.lat. <- d.lat
d.lat.[nrow(d.lat),] <- 0
d.lat.[tmp1,] <- 1

x <- my.rw.graph(d=d.lat.,n.step=100,pr=1)
matplot(x,type="l")

d.lat. <- d.lat

tmp2 <- sample((2^p/2):(2^p),50)
tmp2 <- c(tmp2,2^p)
for(i in 1:length(tmp2)){
	d.lat.[tmp2[i],1:tmp2[i]] <- 0
}
tmp3 <- sample(tmp2,length(tmp2)/2)
for(i in 1:length(tmp3)){
	d.lat.[tmp3[i],tmp3[i]:(2^p)] <- 0
}

x <- my.rw.graph(d=d.lat.,n.step=100,pr=1)
matplot(x,type="l")

k <- 1+runif(1)*0.5
d.lat.off <- d.lat.on <- d.lat.
d.lat.off[upper.tri(d.lat.off)] <- d.lat.off[upper.tri(d.lat.off)] * k
d.lat.on[!upper.tri(d.lat.on)] <- d.lat.on[!upper.tri(d.lat.on)] * k

x.off <- my.rw.graph(d=d.lat.off,n.step=100,pr=0.5)
x.on <- my.rw.graph(d=d.lat.on,n.step=100,pr=0.5)

par(mfcol=c(1,2))
matplot(x.off,type="l")
matplot(x.on,type="l")
par(mfcol=c(1,1))

f <- 0.95
mjr <- which(x.off[nrow(x.off),]+x.on[nrow(x.on),] > quantile(x.off[nrow(x.off),]+x.on[nrow(x.on),],1-50/(2^p)))
par(mfcol=c(1,2))
matplot(cbind(x.off[nrow(x.off),mjr],x.on[nrow(x.on),mjr]),type="l")
plot(x.off[nrow(x.off),mjr],x.on[nrow(x.on),mjr])
abline(0,1,col=2)
par(mfcol=c(1,1))

library(MCMCpack)
d.lat. <- d.lat

tmp2 <- sample((2^p/2):(2^p),15)
tmp2 <- c(tmp2,2^p)
for(i in 1:length(tmp2)){
	d.lat.[tmp2[i],1:tmp2[i]] <- 0
}
tmp3 <- sample(tmp2,length(tmp2)/2)
for(i in 1:length(tmp3)){
	d.lat.[tmp3[i],tmp3[i]:(2^p)] <- 0
}


x.inits <- rdirichlet(9,rep(1,2^p))
par(mfcol=c(3,3))
xs <- list()
for(i in 1:length(x.inits[,1])){
	xs[[i]] <- my.rw.graph(d=d.lat.,n.step=1000,pr=0.5,x.init=x.inits[i,])
	matplot(xs[[i]],type="l")
}
par(mfcol=c(1,1))
e.out <- eigen(d.lat.)
plot(Arg(e.out[[1]]))
plot(e.out[[1]])

何回のトライアルをすれば、すべてのタイプが少なくとも１回は起きるか、そして起きてほしくないタイプは１回も起きないか…