医学のための数学~医学概論2014~その2

## hash-2.2.6 provided by Decision Patterns

plot of chunk unnamed-chunk-1

行列からグラフへ

4要素が次のような連立微分方程式を満たすような関係にあるときに

\[ \begin{equation} \frac{dx_1}{dt} = -x_2\\ \frac{dx_2}{dt} = -x_3\\ \frac{dx_3}{dt} = -x_4\\ \frac{dx_4}{dt} = x_1 \end{equation} \]

こんな増減変化が現れ

A <- matrix(c(0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 1, 0, 0, 0), byrow = TRUE, 
    4, 4)
eigen.out <- eigen(A)
V <- eigen.out[[2]]
U <- solve(V)
s <- eigen.out[[1]]
t <- seq(from = 0, to = 1, length = 1000) * 10
X <- matrix(0, length(t), 4)
x1 <- c(1, 3, 10, 100)
for (i in 1:length(t)) {
    X[i, ] <- V %*% diag(s^t[i]) %*% U %*% x1
}
matplot(Re(X), type = "l")

plot of chunk unnamed-chunk-2

そこには \[ \begin{equation} \begin{pmatrix}0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \\ 1 & 0 & 0 & 0 \end{pmatrix} \end{equation} \]

こんな行列が関わっているのでした。

この行列をよく見ると、\( x_i \)が\( x_j \)に影響を受けて減少するときには、i行j列に-1が立ち、逆に増加するときには、1が立つことを表しており、\( x_i,x_j \)の間に直接の関係がないときには0が立っていることを表しています。

生命現象を理解するときには、たくさんの要素の間の関係の有無を捉えて理解しようとすることがあります。

n <- 10
M <- diag(rep(1, n))
M <- M[c(2:n, 1), ]
library(igraph)
## Warning: package 'igraph' was built under R version 3.0.2
g <- graph.adjacency(M)
plot(g)

plot of chunk unnamed-chunk-3

これは代謝系に出てくる分子サイクルに見えます。

絵で描けばこのようになりますが、これを行列で表示すれば

M
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    0    1    0    0    0    0    0    0    0     0
##  [2,]    0    0    1    0    0    0    0    0    0     0
##  [3,]    0    0    0    1    0    0    0    0    0     0
##  [4,]    0    0    0    0    1    0    0    0    0     0
##  [5,]    0    0    0    0    0    1    0    0    0     0
##  [6,]    0    0    0    0    0    0    1    0    0     0
##  [7,]    0    0    0    0    0    0    0    1    0     0
##  [8,]    0    0    0    0    0    0    0    0    1     0
##  [9,]    0    0    0    0    0    0    0    0    0     1
## [10,]    1    0    0    0    0    0    0    0    0     0

となります。

2つのつながり具合の異なるネットワークを描いてみます。

g <- erdos.renyi.game(500, 1/200)
plot(g, vertex.size = 5, vertex.label = NA)

plot of chunk unnamed-chunk-5

g <- as.undirected(barabasi.game(500))
plot(g, vertex.size = 5, vertex.label = NA)

plot of chunk unnamed-chunk-5

どちらも、あっさりとした設定で作成したネットワークですが、見た印象はずいぶん違います。 生命現象を支えている分子ネットワークや機能ネットワークがどちらのパターンに近いのか、それともこれらとは異なるパターンなのかを明らかにすることは、構成要素の1つ1つを細かく調べることと同じくらい重要なことです。 1つ1つが分子である場合には、その構成要素を1つ1つ細かく調べるのは「分子生物学」です。 全体像を把握しに行くのは「システムバイオロジー」です。 全体像の把握にはそれに合ったアプローチがあり、ネットワーク理論やそれを支えるグラフ理論は重要です。

親子関係・伝承関係のグラフ

感染微生物はヒトの健康の大敵ですが、その微生物はたえず変化しています。その変遷の様子を捕まえるために、微生物の遺伝情報である核酸塩基配列を調べ、その継承・変遷関係を捕まえる必要があります。以下の例は、免疫低下を引き起こす疾患であるAIDSの原因ウイルス(HIVウイルス)の塩基配列変異に基づく系統分岐を描いたものです。

library(ape)
## Warning: package 'ape' was built under R version 3.0.2
## 
## Attaching package: 'ape'
## 
##  以下のオブジェクトはマスクされています (from 'package:igraph') : 
## 
##      edges
# example tree in NH format (a string)
data("hivtree.newick")
# generate file 'hivtree.phy' in working directory
cat(hivtree.newick, file = "hivtree.phy", sep = "\n")
tree.hiv <- read.tree("hivtree.phy")  # load tree
unlink("hivtree.phy")  # delete the file 'hivtree.phy'
plot(tree.hiv, show.tip.label = FALSE)

plot of chunk unnamed-chunk-6

ヒトの病気が家系内集積することがあり、それを引き起こしている遺伝因子や環境因子を突き止めることも重要で、家系図を解析することも有用です。

家系図もグラフの一種です。家系図では、夫婦に水平線を引き、その中点から子への線を引くので線が折れ線になりますが、父親と子、母親と子、の間に線を引くことにすれば、いわゆるグラフになります。

どちらのグラフも「絵」として眺めることも重要ですが、そこにどんな原理があるのかを解釈しようとするとき、「絵」の状態では難しいです。 しかしながら、グラフは「行列」表現できることをもう知っていますから、数学的に扱うことが可能です。

## Warning: package 'kinship2' was built under R version 3.0.2
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 3.0.2
## Loading required package: quadprog

plot of chunk unnamed-chunk-7

## Did not plot the following people: 2 5 6 7 8

少ないルールが作る木と言う構造

グラフの中でもぐるりと閉じた経路がないものに「木」というタイプのグラフがあります。

多細胞生物であるヒトの体の構造には分岐木状の構造があります。 身体に張り巡らされた血管系、口から肺へと通じる気道・気管系がその例です。 単一の細胞である受精卵からだんだんに特徴を変化させながら特別な形や機能を持った細胞に分化していくことも、「多細胞が作る機能的『分岐木』」です。 生物の進化を系統樹で表したものも『分岐木』です。

植物の木が成長する様子も、体内で木構造が成長する様子も基本は同じです。

plot of chunk unnamed-chunk-8

こんな感じでだんだんに複雑になっていきます。

しかしながら、この複雑な構造を作るルールは単純です。

5つの記号"F",“+”,“-”,“[”,“]"を用意し、 "F"をスタートとして、その後は \[ \begin{equation} "F" \to "FF-[-F+F+F]+[+F-F-F]" \end{equation} \] と入れ替える、という処理を繰り返すだけで、ある特徴的な"F,+,-"の文字列を作成することができます。

できた文字列を2次元に表示する決まりを与えると、ここに示したような木構造の成長が作成できます。

入れ替え文字列のパターンを変えると、木の様子も変わります。 以下のような入れ替えパターンで描いてみます。 \[ \begin{equation} "F" \to "F[+F]F[-F]F"\\ "F" \to "F[+F]F[-F][F]"\\ "F" \to "FF-[-F+F+F]+[+F-F-F]"\\ "F" \to "F[+X]F[-X]+X"\\ ("F","X") \to ("FF","F[+X]F[-X]+X")\\ ("F","X") \to ("FF","F[+X][-X]FX")\\ \end{equation} \] plot of chunk unnamed-chunk-9

L-systemと呼ばれる形式文法の1つで植物の成長プロセスを初めとした様々な自然物の構造を記述・表現できるアルゴリズムです。

限られた構成要素と少数の規則で複雑な構造を作り上げる仕組みの1つです。フラクタルなどもこの方法で構成することができます。

生命体・生命現象の多様性が限られた遺伝子の情報に詰め込むには、こんな仕組みも背景にありそうです。

library(fields)  # for tim.colors
## Warning: package 'fields' was built under R version 3.0.2
## Loading required package: spam
## Warning: package 'spam' was built under R version 3.0.2
## Loading required package: grid
## Spam version 0.40-0 (2013-09-11) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## 
## Attaching package: 'spam'
## 
##  以下のオブジェクトはマスクされています (from 'package:base') : 
## 
##      backsolve, forwardsolve 
## 
## Loading required package: maps
munit <- 10  # 解像度指定 (30は重い、20くらいが適当)
m = munit^2  # grid size
C = complex(real = rep(seq(-1.8, 0.6, length.out = m), each = m), imag = rep(seq(-1.2, 
    1.2, length.out = m), m))
C = matrix(C, m, m)
L <- 1
# 複素数zを0+0iから、1+0iに増やしたり
# |a+bi|一定で、偏角を0からpiに増やしたりする modul<-rep(0.2,L) #
# 偏角漸増(|a+bi|固定)
modul <- seq(from = 0, to = 1, length.out = L)  # 実軸に長くする
# args<-seq(from=0,to=pi,length.out=L) # 偏角漸増)
args <- rep(4 * pi/10, length.out = L)  # 偏角固定
Y <- array(0, c(m, m, L))
for (i in 1:L) {
    Z = complex(m = modul[i], a = args[i])  # z(n+1)=z(n)^2+Cのz
    K <- 20
    X = array(0, c(m, m, K))
    for (k in 1:K) {
        # X = array(0, c(m,m,20)) for (k in 1:20) {
        Z = Z^2 + C
        X[, , k] = exp(-abs(Z))
    }
    image(X[, , K], col = tim.colors(256))  # show final image in R
    Y[, , i] <- X[, , K]
}

plot of chunk unnamed-chunk-10

グラフを使って医療診断~ベイジアンネットワーク~

患者さんが病気の好発地域の出身者だったら、その病気かもしれない確率は高い。そんな患者さんが、その病気かもしれない症状で受診した。 検査をしてそれが異常か正常かで診断にたどり着きたい。 検査が異常になるか正常になるかは、病気の有無にも左右されるが、患者さんの生活習慣にも左右される。 また、同じような症状の別の病気でも同じような検査結果が出る。 さて、問診と検査で診断しよう。

これは普通の臨床診断プロセスです。

これにグラフを持ち込んでみる方法の1つがベイジアンネットワークです。

医療上の判断に計算機の支援を導入するのはまだ現実的とは考えられていませんが、膨大になる医学知識のすべてを活用することも1人の人間の能力を超えているのも確かです。

どのようにして計算機を活用していくかの入口を眺めておくことは決して悪いことではありません。

library(gRain)
## Warning: package 'gRain' was built under R version 3.0.2
## Loading required package: gRbase
## Warning: package 'gRbase' was built under R version 3.0.2
## Loading required package: graph
## 
## Attaching package: 'graph'
## 
##  以下のオブジェクトはマスクされています (from 'package:ape') : 
## 
##      edges 
## 
##  以下のオブジェクトはマスクされています (from 'package:igraph') : 
## 
##      degree, edges 
## 
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 3.0.2
yn <- c("yes", "no")
a <- cptable(~asia, values = c(1, 99), levels = yn)
t.a <- cptable(~tub + asia, values = c(5, 95, 1, 99), levels = yn)
s <- cptable(~smoke, values = c(5, 5), levels = yn)
l.s <- cptable(~lung + smoke, values = c(1, 9, 1, 99), levels = yn)
b.s <- cptable(~bronc + smoke, values = c(6, 4, 3, 7), levels = yn)
e.lt <- cptable(~either + lung + tub, values = c(1, 0, 1, 0, 1, 0, 0, 1), levels = yn)
x.e <- cptable(~xray + either, values = c(98, 2, 5, 95), levels = yn)
d.be <- cptable(~dysp + bronc + either, values = c(9, 1, 7, 3, 8, 2, 1, 9), 
    levels = yn)
plist <- compileCPT(list(a, t.a, s, l.s, b.s, e.lt, x.e, d.be))
pn <- grain(plist)
pn
## Independence network: Compiled: FALSE Propagated: FALSE 
##   Nodes: chr [1:8] "asia" "tub" "smoke" "lung" "bronc" "either" ...
summary(pn)
## Independence network: Compiled: FALSE Propagated: FALSE 
##  Nodes : chr [1:8] "asia" "tub" "smoke" "lung" "bronc" "either" ...
plot(pn)
## Loading required package: Rgraphviz

plot of chunk unnamed-chunk-11

pnc <- compile(pn, propagate = TRUE)

「アジア出身かどうかで結核(tub)のリスクが影響される」 「喫煙者かどうか(smoke)で肺がん(lung)かどうかのリスクが影響される」 「喫煙者かどうかで慢性気管支炎(bronc)かどうかのリスクが影響される」 「結核か肺がんのどちらか(either)が気になるが、そのどちらかは結核であるかどうか、肺がんであるかどうかが決める」 「結核か肺がんかのどちらか(either)は呼吸困難感(dysp)に影響を与え、慢性気管支炎も影響を与える」 「結核か肺がんかのどちらか(either)は胸部レントゲン写真(xray)に影響を与える」

という個々の情報をグラフにしたのがこの図です。

このネットワークでは、「アジア出身かどうか」「喫煙者かどうか」「レントゲン写真の所見」などの情報を投入すると、どの診断がそれらしいかを数字にして返してくれます。

このベイジアンネットワークは、このような診断の補助的な情報提供にも使えますし、実験データから意味を読み取るための手法としても使えます。