14-8. k-meansで細胞分裂シミュレーション

やること
k-meansについて
細胞分裂について
実行環境
細胞を用意する
点群をk-meansする関数
k-meansを可視化する関数
コードの改良
楕円形チャンバーで16細胞期まで
三角形チャンバーで16細胞期まで
L字型チャンバーで16細胞期まで
楕円切り欠きチャンバーで16細胞期まで
星型チャンバーで16細胞期まで
考察・展望（レポートかな？）

やること

k-means法を可視化していたら気がつきました。

細胞分裂後期っぽいなあ。

ということで今回はk-meansを応用して細胞分裂のシミュレーションを行います。

k-meansについて

k-means法のアルゴリズムと可視化についてはこちらの動画が分かりやすいです。

アルゴリズムがデータを分ける様子を可視化した【K-means】【クラスタリング】

こちらも同様に分かりやすい。

K-means 法を D3.js でビジュアライズしてみた

クラスタリングの定番アルゴリズム K-means 法（K平均法）の動作原理を理解するために、D3.js を使って可視化してみました。ステップ最初からN (ノード数):K (クラスター数):新規作成図をクリックするかボタンを押すと、1ステップずつ処理を行いますボタンを押すと、最初の状態に戻りますボタンを押すと、...

細胞分裂について

Nicolas Minc, et al., Influence of Cell Geometry on Division-Plane Positioning, Cell 144(3), 414-426, 2011をご参照ください。アブストをGoogle翻訳にぶっ込み、図をサーっと見ます。

媒精したウニ卵をさまざまな形状のマイクロチャンバーに閉じ込め、第一卵割面の決まり方を実験と数理モデルで検討した論文です。

実行環境

WinPython3.6をおすすめしています。

WinPython - Browse /WinPython_3.6/3.6.7.0 at SourceForge.net

Portable Scientific Python 2/3 32/64bit Distribution for Windows

Google Colaboratoryが利用可能です。

Google Colab

細胞を用意する

細胞に見立てた円状の点群を用意します。

import numpy as np
import numpy.random as nr
import matplotlib.pyplot as plt

#============================
# 設定（細胞の形状）
#============================
#円形
data = nr.rand(1000, 2)
mask = np.where(((data[:, 0] - 0.5)**2 + (data[:, 1] - 0.5)**2)**0.5 < 0.5)[0]
data = data[mask]

#============================
# 表示
#============================
plt.figure(figsize=(6, 6))
plt.scatter(data[:, 0], data[:, 1], c='k')
plt.xlim(0, 1); plt.ylim(0, 1)
plt.show(); plt.close(); print()

点群をk-meansする関数

sklearn.cluster.KMeans が使えればよいのですが、１ステップ毎に止めて可視化するような機能がなかったので、自分でk-meansを書きます。

ここでは「適当な重心」と「点群」を受け取って、k-meansを50ステップ進めてから「２つの重心の時系列データ」と「２つの点群の時系列データ」を返す関数を作りました。

#============================
# 重心と点群を受け取り、重心を２つに分裂させてk-meansする関数
#============================
def division(center, data):
    #重心の分裂
    center_1 = center + (nr.rand(2) * 0.01 - 0.005)
    center_2 = center + (nr.rand(2) * 0.01 - 0.005)
    
    #データ
    center_time_1 = []
    center_time_2 = []
    data_time_1 = []
    data_time_2 = []
    
    #k-meansを50ステップ進める
    for step in range(50):
        print('step:{}'.format(step))
        
        #２つの重心から各点までの距離計算
        dist_1 = ((data[:, 0] - center_1[0])**2 + (data[:, 1] - center_1[1])**2)**0.5
        dist_2 = ((data[:, 0] - center_2[0])**2 + (data[:, 1] - center_2[1])**2)**0.5
        
        #重心１に属するかどうかのbool
        bool_1 = dist_1 < dist_2
        
        #各重心に属する点群
        data_1 = data[bool_1, :]
        data_2 = data[~bool_1, :]
        
        #重心と所属点群を記録
        center_time_1.append(center_1)
        center_time_2.append(center_2)
        data_time_1.append(data_1)
        data_time_2.append(data_2)
        
        #重心の更新（配下の点群の重心に移動）
        center_1 = np.mean(data_1, axis=0)
        center_2 = np.mean(data_2, axis=0)
    
    #２つの重心と２つの点群の時系列データを返す（点群データのサイズはステップ毎に変わるのでndarrayにしないよう注意）
    return center_time_1, center_time_2, data_time_1, data_time_2

#適当な初期重心        
center = [nr.rand()*0.02 + 0.49, nr.rand()*0.02 + 0.49]

#分裂して２細胞期へ
center_time_1, center_time_2, data_time_1, data_time_2 = division(center, data)

print('len(center_time_1):\n{}'.format(len(center_time_1)))
print('len(data_time_1[0]):\n{}'.format(len(data_time_1[0])))
print('len(data_time_2[0]):\n{}'.format(len(data_time_2[0])))

len(center_time_1):
50
len(data_time_1[0]):
381
len(data_time_2[0]):
382

適当な初期重心と先ほどの点群を入れて、返り値のサイズを見てみました。重心１は(x, y)が50ステップ並んだ２次元配列です。１ステップ目に重心１に属した点群は381個、重心２に属した点群は382個でした。

k-meansを可視化する関数

何細胞期でも可視化できるように拡張性を意識して、可視化関数を作りました。

#============================
# 複数の重心、複数の点群の時系列データを受け取ってプロットする関数
#============================
def plot_frame(center_time_list, data_time_list):
    #50ステップ表示
    for step in range(50):
        #表示
        plt.figure(figsize=(6, 6))
        #各細胞の重心とデータをプロット
        for i in range(len(center_time_list)):
            #この細胞、このステップの情報
            center = center_time_list[i][step]
            data = data_time_list[i][step]
            #点群
            plt.scatter(data[:, 0], data[:, 1])
            #重心
            plt.scatter(center[0], center[1], s=100)
            #線（重い）
            #for j in range(len(data)):
            #    plt.plot([data[j, 0], center[0]], [data[j, 1], center[1]], '-k', alpha=0.2)
        
        plt.xlim(0, 1); plt.ylim(0, 1)
        plt.title('step:{}'.format(i))
        plt.show(); plt.close(); print()

#listに格納
center_time_list = [center_time_1, center_time_2]
data_time_list = [data_time_1, data_time_2]

#１→２細胞期を可視化
plot_frame(center_time_list, data_time_list)

いい感じです。

コードの改良

要領が分かりました。まず、重心の移動が急すぎるので、移動の速さを30%くらいに落とします。

#重心の更新（配下の点群の重心に移動）
#center_1 = np.mean(data_1, axis=0)
#center_2 = np.mean(data_2, axis=0)
center_1 = center_1*0.7 + np.mean(data_1, axis=0)*0.3
center_2 = center_2*0.7 + np.mean(data_2, axis=0)*0.3

16細胞期までシミュレーションするコードです。

#適当な初期重心        
center = [nr.rand()*0.02 + 0.49, nr.rand()*0.02 + 0.49]
#分裂して２細胞期へ
center_time_1, center_time_2, data_time_1, data_time_2 = division(center, data)
#listに格納
center_time_list = [center_time_1, center_time_2]
data_time_list = [data_time_1, data_time_2]
#１→２細胞期を可視化
plot_frame(center_time_list, data_time_list)


#分裂して４細胞期へ、listに格納
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[0][-1], data_time_list[0][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[1][-1], data_time_list[1][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
#２→４細胞期を可視化
plot_frame(center_time_list[-4:], data_time_list[-4:])


#分裂して８細胞期へ、listに格納
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[2][-1], data_time_list[2][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[3][-1], data_time_list[3][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[4][-1], data_time_list[4][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[5][-1], data_time_list[5][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
#４→８細胞期を可視化
plot_frame(center_time_list[-8:], data_time_list[-8:])


#分裂して16細胞期へ、listに格納
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[6][-1], data_time_list[6][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[7][-1], data_time_list[7][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[8][-1], data_time_list[8][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[9][-1], data_time_list[9][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[10][-1], data_time_list[10][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[11][-1], data_time_list[11][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[12][-1], data_time_list[12][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
center_time_1, center_time_2, data_time_1, data_time_2 = division(center_time_list[13][-1], data_time_list[13][-1])
center_time_list += [center_time_1, center_time_2]
data_time_list += [data_time_1, data_time_2]
#８→16細胞期を可視化
plot_frame(center_time_list[-16:], data_time_list[-16:])

楕円形チャンバーで16細胞期まで

#============================
# 設定（細胞の形状）
#============================
#楕円
data = nr.rand(1000, 2)
mask = np.where(((data[:, 0] - 0.5)**2 + (data[:, 1] - 0.5)**2)**0.5 < 0.5)[0]
data = data[mask]
data[:, 1] = data[:, 1]*0.7 + 0.15

疎な点群

密な点群

論文より比較用の画像を拝借します。

三角形チャンバーで16細胞期まで

#============================
# 設定（細胞の形状）
#============================
#三角
data = nr.rand(1000, 2)
mask = (data[:, 1] < 3.0*data[:, 0] - 0.5) * (data[:, 1] < -3.0*data[:, 0] + 2.5)
data = data[mask]

疎な点群

密な点群

論文より比較用の画像を拝借します。

L字型チャンバーで16細胞期まで

#============================
# 設定（細胞の形状）
#============================
#L字
data = nr.rand(1000, 2)
mask = (data[:, 0] > 0.5) * (data[:, 1] > 0.5)
data = data[~mask]

疎な点群

密な点群

楕円切り欠きチャンバーで16細胞期まで

#楕円切り欠き
data = nr.rand(1000, 2)
mask1 = np.where(((data[:, 0] - 0.5)**2 + (data[:, 1] - 0.5)**2)**0.5 < 0.5)[0]
data = data[mask1]
data[:, 1] = data[:, 1]*0.7 + 0.15
mask2 = (data[:, 1] < 4.0*data[:, 0] - 0.9) * (data[:, 1] < -4.0*data[:, 0] + 1.7)
data = data[~mask2]
mask3 = (data[:, 1] > 4.0*data[:, 0] - 0.7) * (data[:, 1] > -4.0*data[:, 0] + 1.9)
data = data[~mask3]

疎な点群

密な点群

論文によれば、切り欠きのところで割れるわけではない。

星型チャンバーで16細胞期まで

星型を作るのが面倒くさかったです。。

#星型
data = nr.rand(12000, 2)
mask1 = (data[:, 1] > 0.65) * (data[:, 1] > 3.08*data[:, 0] - 0.5)
data = data[~mask1]
mask2 = (data[:, 1] > 0.65) * (data[:, 1] > -3.08*data[:, 0] + 2.5)
data = data[~mask2]
mask3 = (data[:, 1] > 0.65) * (data[:, 1] > -3.08*data[:, 0] + 2.5)
data = data[~mask3]
mask4 = (data[:, 1] < -0.73*data[:, 0] + 0.65) * (data[:, 1] > 3.08*data[:, 0] - 0.5)
data = data[~mask4]
mask5 = (data[:, 1] < 0.73*data[:, 0] - 0.05) * (data[:, 1] > -3.08*data[:, 0] + 2.5)
data = data[~mask5]
mask6 = (data[:, 1] < 0.73*data[:, 0] - 0.05) * (data[:, 1] < -0.73*data[:, 0] + 0.65)
data = data[~mask6]

疎な点群