1. No category

第2回数理デザイン道場@IBM天城ホームステッド
ベイズ統計と生細胞画像解析
Terumasa Tokunaga [1][2], Osamu Hirose [2][3], Yu Toyoshima [2][4], !
Takayuki Teramoto [2][5], Sayuri Kuge [2][5], Takeshi Ishihara [2][5], !
Yuichi Iino [2][4], Yuishi Iwasaki [2][7] and Ryo Yoshida [1][2][6][8]
[1] The Institute of Statistics Mathematics, Research Organization and Systems, [2] JST, CREST,!
[3] Kanazawa University, [4] The University of Tokyo, [5] Kyushu University, !
[6] The Graduate University for Advanced Studies, [7] Ibaraki University, [8] JST, ERATO
第2回数理デザイン道場@IBM天城ホームステッド
1
課題の設定: 動物の神経系の理解
線虫の神経系
神経細胞数: 302（雌雄同体）.!
感覚/介在/運動ニューロン.!
ASEL
化学シナプス結合とギャップ結合.
ASEL
AIAL
神経回路モデル
特徴
全ての神経間結合が既知.!
ゲノミクスが進んでおり, 突然変異体の作成
が容易.
多数の細胞特異的プロモーターが知られており,
遺伝子導入により特定の細胞での遺伝子発現
が可能.
＊
城大•岩崎先生の講義資料より
第2回数理デザイン道場@IBM天城ホームステッド
2
課題の設定
Caイオンイメージングによる神経活動度の可視化
20-24
＊九大•石原研究室より提供
神経活動の指標となるカルシウムイオンを鋭敏に検出できる, 蛍光性カルシウム
イオンセンサー（蛍光タンパク質）が開発.
特定の神経細胞に発現. 外部刺激に対する神経系の応答を直接計測.
第2回数理デザイン道場@IBM天城ホームステッド
課題の設定: 線虫神経系のまるごと計測
mCherryによる多数の神経細胞核の形状の可視化
＊九大•石原研究室より提供
およそ150個の神経細胞を同時にモニタリング可能な計測技術が確立.
Caイオンイメージングにより, 神経細胞活動度も同時に計測.
3
第2回数理デザイン道場@IBM天城ホームステッド
課題の設定: ロードマップ
神経活動「まるごと」計測データと神経回路モデルを統合: データ同化
➡
神経回路の動作原理の包括的な理解
4
第2回数理デザイン道場@IBM天城ホームステッド
4Dイメージングデータからの!
神経活動度の自動定量化
第2回数理デザイン道場@IBM天城ホームステッド
神経活動度の自動定量化
Data characteristics
Cells form to be globular or more general.
Sizes and shapes of cells are very similar.
Number of cells = ~ 100-160 (closely-spaced)
Move rapidly and coordinately.
5
第2回数理デザイン道場@IBM天城ホームステッド
6
神経活動度の自動定量化
Tasks for automated quantification
➡
➡
Cell detection
Multiple cell tracking
Outline of the proposed method
(a)
(b)
Detection
Conversion
（kernel density estimation）
(d)
(c)
Tracking
Segmentation
第2回数理デザイン道場@IBM天城ホームステッド
7
Step 1. Conversion of 3D image
Original data
Coordinate of voxel: xi (i = 1, !, 512 × 256 × 20)
Normalized intensity: wi ≥ 0,
n
∑w
i
=1
i=1
Estimated density (continuous function)
Original data (digital)
p=0
Conversion
Kernel density estimation
n
p(x) = ∑ wi k(x − xi )
i=1
Intensity
{
}
k(x − xi ) ∝ exp −(x − xi )T Σ −1 (x − xi )
• Center: position of each voxel!
• Mixing rate: intensity of each voxel
x ∈! 3+
第2回数理デザイン道場@IBM天城ホームステッド
8
Step 1. Conversion of 3D image
Assumption:
Each cell lies around a local mode of the density function.
Cells (
p=0 )
Cell tracking → tracking of local modes changing over time.
detect
t =1
track
t=2
t=3
Local modes can be found by using optimization techniques designed for
a continuous function, for instance, gradient descent algorithm.
Time
第2回数理デザイン道場@IBM天城ホームステッド
Step 2. Cell counting and detection
Hill-climbing algorithm for detecting local modes
Cell detection = Peak detection
Optimization methods for continuous objective functions are applicable.
We show the hill-climbing method derived upon an idea of EM-algoritnm.
Recurrent formula (EM-algorithm)
ψ new = ∑ ui (ψ old )xi
n
i=1
(
)
(
with ui ψ old ∝ wi k ψ old − xi
)
(i) Start from an arbitrary initial position.!
(ii) Each process monotonically climbs a
hill.
(iii) The search process converges at the
nearest local mode.
initial
points
9
第2回数理デザイン道場@IBM天城ホームステッド
10
Step 2. Cell counting and detection
Extension of hill-climbing method
Independent search
Repulsive parallel hill-climbing method
initial points
initial points
improve
Many processes starting from
different initial points
➡ tends to converge to a same local
mode.
repulsion
Different processes explore
different local modes.
We can detect much more
diverse objects with one parallel run.
第2回数理デザイン道場@IBM天城ホームステッド
Step 2. Cell counting and detection
Detection test using real imaging data
search process
movement of search process
45 objects are included in the image.
80 initial points were allocated in a conner of the 3D space.
Repulsive force severity is controlled (lowered gradually until zero).
45 local peaks were detected by one parallel optimization.
11
第2回数理デザイン道場@IBM天城ホームステッド
12
Step 3. Cell tracking
Tracking of multiple objects
Many cells are distributed densely.!
This causes two kind of tracking errors;
correct tracking
t =1
t=2
(i) merging
t =1
(ii) switching
t=2
Errors:
(i) more than one trackers merge to the same target.!
(ii) trackers switch their targeting local modes.
t =1
t=2
第2回数理デザイン道場@IBM天城ホームステッド
13
Step 3. Cell tracking
Transition model based on an Markov Random Field (MRF)
Spatiotemporal covariation of the movement of cells are modeled.!
Each cell tends to move in conjunction with neighboring cells.
ψ t , j − ψ t−1, j =
∑ (ψ
t,k
− ψ t−1, k + ηk, j
k∈C j
)
with ηk, j ~ N(0, λ k, j )
C j : neighbor set of the j-th tracker
clique
ψ1
ψ2
ψ3
t −1
ψ1
ψ2
ψ3
t
第2回数理デザイン道場@IBM天城ホームステッド
14
Step 3. Cell tracking
Maximization of posterior distribution
Likelihood
maxmize p ( Ψ t | Ψ t−1 , I t )
image force
n
⎛ 1
⎞
∝ ∏⎜
log ∑ k(ψ t , j − xi )wt , i ⎟
⎠
j=1 ⎝ 2γ
i=1
g
ψ t, j
ψ t−1, j
×
likelihood (image force)
(
) (
⎛ 1
× ∏ exp ⎜ −
ψ t , j − ψ t−1, j − ψ t , k − ψ t−1, k
⎝ 2λ
k∈C j
prior distribution (constraint force)
λ k, j ∝ ψ t−1, j − ψ t−1, k
)
⎞
⎟
Σ⎠
2
t −1
t
Prior distribution
−2
Σ
constraint force
dispersion of prior distribution
t −1
t
第2回数理デザイン道場@IBM天城ホームステッド
15
Step 3. Cell tracking
Tracking result
tracker
minimum spanning tree
第2回数理デザイン道場@IBM天城ホームステッド
16
自動定量化に向けた課題
Cell detection
Multiple cell tracking
Segmentation
Image alignment
Skelton estimation
線虫が激しく動くデータでの正確な追跡の実現.
一部の細胞が一時的にフレームからはみ出す場合の対処.
追跡中に発生したswitchingの検出法の構築.
単一細胞（cell body＋dendrite＋axon）の正確な追跡法の開発.
アノテーション付きデータとの自動マッチング法の開発.
計算速度.
開発中
Kernel density estimation
運用中
カーネル密度推定に基づく, 神経活動度の自動定量化法を開発.
第2回数理デザイン道場@IBM天城ホームステッド
神経回路モデルと計測データの統合!
: データ同化
第2回数理デザイン道場@IBM天城ホームステッド
実験データと神経回路モデルの統合
データ同化研究の事例
統計数理研究所データ同化研究開発センター
http://daweb.ism.ac.jp/contents/index.php
宇宙空間の電磁現象の解明
細胞内流動の可視化
台風の進路予測
48 hrs. prediction
48 hrs. prediction
24 hrs. prediction
24 hrs. prediction
24 hrs. observation data
are accumulated…
データ同化：地球科学分野で発展した, データと数理モデルの統合解析手法.!
モデルに含まれる物理パラメータの推定.!
適切な初期条件や境界条件の推定.!
部分的な観測データからシステム全体の内部状態の推定.!
現在は工学や経済学, 生物学など広範な応用.
17
第2回数理デザイン道場@IBM天城ホームステッド
19
実験データと神経回路モデルの統合
ベイズの定理
p(x | θ ) × p(θ )
p(θ | x) =
∝ p(x | θ ) × p(θ )
p(x)
事後確率
θ
x
尤度
事前分布
: 潜在変数, モデルのパラメータ
: 観測データ
!
!
!
第2回数理デザイン道場@IBM天城ホームステッド
実験データと神経回路モデルの統合
津波データ同化による海底地形の推定
15
Nakamura et al.,
シミュレーションモデル
⎧ : 浅水波方程式を用いた津波
シミュレータ
⎨
⎩ 観測データ
: 北海道西南沖地震時の
海面変位（潮位データ）
海底地形がわかれば浅水波方程式により津波の到達時間は高精度で予測可能.
海底地形は正確にわかっていない（主要なデータベース間で平均5%の差異）
データ同化により, 津波の伝搬プロセスから海底地形（境界条件）を補正.
日本海の海底地形は, 既存データベースの地形よりも浅いと推定.
➡
地球科学の問題では, validationが本質的に困難
第2回数理デザイン道場@IBM天城ホームステッド
20
実験データと神経回路モデルの統合
神経回路モデルの例: コンダクタンスベースモデル
神経細胞を点素子として取り扱う.
i 番目の神経細胞の膜電位 Vi (t) に関する方程式;
dVi (t)
1
chem
gap gap
stim
C
= − m (Vi (t) − Vi leak ) + ∑ nichem
I
(t)
+
n
I
(t)
+
I
(t).
∑
j
ij
ij ij
i
dt
Ri
j
j
m
i
I
chem
ij
(t) =
gichem
j
1+ exp(− β i j (Vj (t) − V ))
eq
j
chem
ij
(E
− Vi (t)), E
gap
I igap
(t)
=
g
j
i j (V j (t) − Vi (t)).
⎧
⎨
⎩
Vi (t) :
I
chem
ij
I
gap
ij
膜電位
(t) : 化学シナプス結合を通して流入する電流
(t) : ギャップ結合を通して流入する電流
I istim (t) : それ以外の要因による膜電流
nichem
:
j
化学シナプス結合数 →既知
nigapj :
ギャップ結合数
→既知
chem
i j
⎧⎪ EEPSP for excitatory,
=⎨
⎪⎩ E IPSP for excitatory
V j (t )
I ijchem (t )
Vi (t )
I ijgap (t )
第2回数理デザイン道場@IBM天城ホームステッド
21
実験データと神経回路モデルの統合
シミュレーションモデル: 点素子モデル
⎧
⎨
⎩ 観測データ:
内部状態: V(t) = (V1 (t), ... , V302 (t)) ←302個
輝度（カルシウムイオン濃度）:
y(t) = (y1 (t), ..., y m (t)) ←m < 302
状態空間モデル表現
⎧
⎨
⎩
システムモデル: V(t) = f(V(t− 1), θ ) + η (t)
システムノイズ
部分的な!
観測データ
観測モデル: y(t) = g(V(t)) + ε (t)
観測ノイズ
内部状態
計測される量
Vt
It
膜電位
カルシウム応答
nCa + F ! Can F
2+
y(t)
計算
事後分布
p(V | y(t))
推定
全内部状態
V(t)
第2回数理デザイン道場@IBM天城ホームステッド
22
課題
モデル内の未知パラメータ（活性化係数など）の推定法.
統計的な手法（データ同化）により推定
摂動実験によって予め推定
➡
➡
観測されていない神経細胞の膜電位の推定
➡
➡
統計的な手法（データ同化）により推定
Validationはどうするか？
V j (t )
I ijchem (t )
I ijgap (t )
Vi (t )
モデリング
計測
バイオロジーの分野では, 適切な実験を
行うことでvalidationができる！
データサイエンスと計測の専門家が連携し,
効率良く知識の獲得に結びつく計測の形を
デザインして行くことが重要.
統合
検証
計測デザイン
データ同化
検証
モデルの改良