- 確率的に何かが起きる Consider a stochastic process
- 次に何が起きるかは「今」の影響を受けるが、「前」の影響を受けない "Next event" depends on "Now" but not on "Past".
ある確率で、機能性細胞が死滅するとする。その確率はpである。初め、細胞は全部でn個あった。生存細胞数をシミュレーションする。
Assume that cells die stochasitically. The probability of death is p. Simulate for n cells at the begininng.
```{r}
n <- 1000
p <- 0.1
s.or.d <- sample(0:1,n,replace=TRUE,prob=c(1-p,p))
s.or.d[1:100]
sum(s.or.d) # 死細胞数 No. dead cells
sum(s.or.d)/n # 死細胞割合 Fraction of dead cells
```
細胞は確率的に死ぬが、生き残った細胞は増殖する。死亡確率をp、増殖(１細胞が２細胞になる)確率をqとする。初め、細胞は全部でn個あった。
Cells may die but living cells will divide into two. Both are stochastic, with probabilities p and q, respectively.
At the begininng there are n cells.
いくつの細胞が死に、いくつの細胞が分裂するかは、その時々の生細胞数に依存するので、「今」の影響を受けている。
しかしながら、一段階前に何個の細胞が生きていたかには影響を受けていない。
No. death/division depend on the No. cells "NOW", but not on "One period previous".
```{r}
p <- 0.05
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
n.sv.1 <- rep(0,length(weeks))
n.sv.1[1] <- 100
for(i in 2:length(weeks)){
# death
tmp1 <- sample(0:1,n.sv.1[i-1],replace=TRUE,prob=c(1-p,p))
n.death <- sum(tmp1)
# No. living cells for now
tmp2 <- n.sv.1[i-1] - n.death
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living cells for now
n.sv.1[i] <- tmp2 + tmp4
}
plot(n.sv.1)
```
生きている細胞がm 個のとき、1/mの仕事をしなくてはならず、仕事量が多いと、それに応じて死亡率が上がるとする。分裂確率は仕事割合によって変わらないとする。
When there are m cells, each cell has to do 1/m job and the death rate is positively correlated with quantity of job. Division rate is independent from job.
```{r}
p <- 0.05 # p depends on number of living cells
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
n.sv.2 <- rep(0,length(weeks))
n.sv.2[1] <- 100
for(i in 2:length(weeks)){
# death rate
p. <- p + 1/n.sv.2[i-1]
# p. <= 1
p. <- min(p.,1)
# death
tmp1 <- sample(0:1,n.sv.2[i-1],replace=TRUE,prob=c(1-p.,p.))
n.death <- sum(tmp1)
# No. living cells for now
tmp2 <- n.sv.2[i-1] - n.death
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living cells for now
n.sv.2[i] <- tmp2 + tmp4
}
plot(n.sv.2)
```
加速なしの場合と並べてプロットしてみる。
Plot two simulation results together.
```{r}
matplot(cbind(n.sv.1,n.sv.2),type="b")
```
分裂して生まれたばかりの「赤ん坊細胞」は、すぐには分裂できないとするが、次のタイミングでは分裂可能とする。
Baby cells just born do not divide but they become dividable one period later.
この場合は、生細胞を「赤ん坊」か「そうでないか」の２通りで記録する必要がある。
You have to distinct baby cells and non-baby cells.
```{r}
p <- 0.05
q <- 0.03
weeks <- 1:300
# No. living cell from 1st through 100th week
# baby and non-baby
n.sv.3.nb <- n.sv.3.bb <- rep(0,length(weeks))
n.sv.3.nb[1] <- 100
for(i in 2:length(weeks)){
# death of non-baby
tmp1.nb <- sample(0:1,n.sv.3.nb[i-1],replace=TRUE,prob=c(1-p,p))
n.death.nb <- sum(tmp1.nb)
# death of baby
tmp1.bb <- sample(0:1,n.sv.3.bb[i-1],replace=TRUE,prob=c(1-p,p))
n.death.bb <- sum(tmp1.bb)
# No. living non-baby cells for now
tmp2 <- n.sv.3.nb[i-1] - n.death.nb
# Cell division
tmp3 <- sample(0:1,tmp2,replace=TRUE,prob=c(1-q,q))
# No. new born cells
tmp4 <- sum(tmp3)
# No. living non-baby cells. Living baby cells are now non-baby
n.sv.3.nb[i] <- tmp2 + n.sv.3.bb[i-1]-n.death.bb
# No. newborns
n.sv.3.bb[i] <- tmp4
}
plot(n.sv.3.nb+n.sv.3.bb)
matplot(cbind(n.sv.3.nb,n.sv.3.bb,n.sv.3.nb+n.sv.3.bb),type="l")
```
ある疾患には、４つの状態があるという。急性軽症・急性重症・慢性・癌化の４つである。この４状態に、健康・死亡の２状態を加えて、全部で６状態を考える。
Assume a disease has 4 statuses, acureLight,acuteSevere,chronic,cancer. Along these 4, two more statuses, healthy and dead, should be considered.
状態に番号を振る。健康：１、急性軽症：２、急性重症：３、慢性：４、癌化：５、死亡：６とする。
Call 4+6 statuses with 1,2,...,6 as 1:healthy, 2:acuteLight, 3:acuteSevere, 4:chronic, 5:cancer, 6:dead.
健康からは、急性軽症、急性重症、死亡(この疾患とは無関係の死亡)に変わりうるとして、その確率をp12,p13,p16とする。
Healthy individual can become acuteLight/acuteSevere/Dead(dead from some reason away from this disease). Let p12,p13,p16 denote their probability.
変わりえないのは、慢性・癌化の２状態であり、その確率p14=p15=0である。
今、健康なまま、という確率はp11とする。
Healthy individual can not get "chronic/cancer" directry, meaning p14=p15=0. (S)he may stay healthy with probability p11.
このとき、p11 + p12 + p13 + p14 + p15 + p16=1であって、p1iはすべて0以上である。これを、状態１からの推移確率ベクトルと呼ぶ。
p = (p11,...,p16) is called the transition vector from status 1 and p11 + ... + p16 = 1 and p1i >= 0.
同様にして、1-6のそれぞれに推移確率ベクトルを作ることができる。
Each status has its own transition vector.
6つの状態の状態推移ベクトルを合わせると、6x6行列が作れる。これを状態推移行列という。
With these 6 transition vectors, you can make a 6x6 matrix, called "transition matrix".
状態推移行列を使って、6状態にある確率の推移を見てみる。
Using this transition matrix, history of probabilities in 6 statuses can be simulated.
推移確率は１年間での状態変化確率とする。
Transition probability is per-year one.
まず行列を作る。Make a matrix.
```{r}
# From 1:healthy
p1 <- c(0.9,0.08,0.0001,0,0,0.0199)
# From 2:acuteLight
p2 <- c(0.95,0,0,0.04,0,0.01) # 治るか慢性化するか、死亡するか
# From 3:acuteSevere
p3 <- c(0.8,0,0,0.01,0,0.19) # 死亡率0.19
# From 4:chronic
p4 <- c(0.01,0,0,0.95,0.03,0.01) # 癌化する
# From 5:cancer
p5 <- c(0,0,0,0,0.9,0.1)
# From 6: dead
p6 <- c(0,0,0,0,0,1)
P <- cbind(p1,p2,p3,p4,p5,p6)
P
apply(P,2,sum)
```
初期状態は、1:Healthyに居たとする。状態確率ベクトルが(1,0,0,0,0,0)である、という。
Initially you are healthy(1), then your status vector is (1,0,0,0,0,0).
```{r}
v1 <- c(1,0,0,0,0,0)
# １年後 One-year later
v2 <- P %*% v1
v2
# ２年後、two-year later
v3 <- P %*% v2
v3
x <- 1:50
v.all <- matrix(0,length(x),length(v1))
v.all[1,] <- v1
for(i in 2:length(x)){
v.all[i,] <- P %*% v.all[i-1,]
}
matplot(v.all,type="l")
```
慢性状態を治癒に持ち込む確率が、0.01だが、これを0.1に上げる治療A。
Assume a treatment A that increases p41 from 0.01 to 0.1.
```{r}
PA <- P
PA[1,4] <- 0.1
PA[4,4] <- PA[4,4]-(0.1-0.01)
PA
```
```{r}
v.all.A <- matrix(0,length(x),length(v1))
v.all.A[1,] <- v1
for(i in 2:length(x)){
v.all.A[i,] <- PA %*% v.all.A[i-1,]
}
# 治療法Aの効果を死亡状態の割合で比較する
matplot(x,cbind(v.all[,6],v.all.A[,6]),type="l")
```
癌を治す治療Bとして、治癒成功率を0.1にしてみよう。
Assume a treatment B that cures cancer, p51 = 0.1.
```{r}
PB <- P
PB[1,5] <- 0.1
PB[5,5] <- PB[5,5]-0.1
PB
```
```{r}
v.all.B <- matrix(0,length(x),length(v1))
v.all.B[1,] <- v1
for(i in 2:length(x)){
v.all.B[i,] <- PB %*% v.all.B[i-1,]
}
# 治療法Aの効果を死亡状態の割合で比較する
matplot(x,cbind(v.all[,6],v.all.A[,6],v.all.B[,6]),type="l")
```
make.kame.hyoushi <- function(ttl="",n.kame=30,kame.col="blue",sub=""){
t <- seq(from=0,to=1,length=1000)*2*pi
kora <- cbind(cos(t),sin(t))
s <- 1/4
teashi <- kora*s
thetas <- c(0,1,2,4,5)/6*2*pi
teashi.ctr <- cbind(cos(thetas),sin(thetas))
teashi.coords <- matrix(0,nrow=0,ncol=2)
for(i in 1:length(thetas)){
teashi.coords <- rbind(teashi.coords,cbind(teashi.ctr[i,1]+teashi[,1],teashi.ctr[i,2]+teashi[,2]))
}
teashi.out <- teashi.coords[which(apply(teashi.coords^2,1,sum)>1),]
shippo.st <- c(cos(pi),sin(pi))
shippo.s <- 0.1
shippo.l <- 0.4
shippo.x <- seq(from=0,to=-shippo.l,length=100)
shippo.y <- shippo.s * sin(shippo.x*2*pi/shippo.l)
shippo.coords <- cbind(shippo.st[1] + shippo.x,shippo.st[2]+shippo.y)
kame <- rbind(kora,teashi.out,shippo.coords)
b <- 0.9
kame2 <- kame
kame2[,2] <- kame2[,2]*b
X <- matrix(0,nrow=0,ncol=2)
for(i in 1:n.kame){
theta <- runif(1)*2*pi
R <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),2,2)
new.kame <- t(R %*% t(kame2))
mv <- runif(2)*n.kame/2
new.kame <- cbind(new.kame[,1]+mv[1],new.kame[,2]+mv[2])
X <- rbind(X,new.kame)
}
jpeg(filename = paste(ttl,".jpeg"), width = 5000, height = 5000,
quality = 100, pointsize = 180,bg = "white", res = NA,
restoreConsole = TRUE)
plot(X,col=grey(1),axes=FALSE,xlab="",ylab="",main=paste(ttl,"\n",sub))
points(X,cex=0.5,pch=20,xlim=range(X),ylim=range(X),col=kame.col)
dev.off()
}
n.kame <- 15
kame.col <- "orange"
ttl <- "状態推移シミュレーション"
subttl <- "〜臨床的例を用いて〜"
make.kame.hyoushi(ttl,n.kame,kame.col,subttl)
