Help us understand the problem. What is going on with this article?

Pythonで粘菌ネットワーク

概説

粘菌数理モデルが最短経路を求める様子をgif動画にしてみた

ezgif.com-optimize.gif

ezgif.com-optimize (2).gif

粘菌ネットワークとは

粘菌は餌が近くにある時、自身の体で栄養を運ぶ管のネットワークを構築します。そのネットワークはしばしば餌場間の輸送効率の良いネットワークを構築することが知られています。そんな粘菌のふるまいにヒントを得てつくられたのが粘菌数理モデルです。中垣らはこのような粘菌に関する研究によって2008年と2010年に栄えあるイグノーベル賞を受賞しています。

粘菌数理モデル

はじめに、粘菌の体を$n$個の節(ノード)とそれどうしをつなぐ栄養管(エッジ)から構成されるグラフとみなし、ノードとエッジを以下のように表記します。

記号 意味
$N_i$ $i$番目のノード($i$=1,2,3,$\cdots$+$n$)
$M_{ij}$ $i$番目のノードと$j$番目のノードをつなぐエッジ

$M_{ij}$の中を$N_i$から$N_j$方向に流れる単位時間当たりの餌の流量$Q_{ij}$はHagen-Poiseuille流れを仮定すると以下のように書き表されます。

Q_{ij}=\frac{\pi r^4}{8\eta L_{ij}} (p_i-p_j)\tag{1}

物理量の意味は以下の表を参照してください。

物理量 意味
$\pi$ 円周率
$r$ 栄養管の半径
$\eta$ 栄養管を流れる流体の粘性率
$p_i$ $i$番目のノード点での圧力
$p_j$ $j$番目のノード点での圧力
$L_{ij}$ $M_{ij}$の経路長

ここで、

D_{ij}=\frac{\pi r^4}{8\eta}\tag{2}

となる変数$D_{ij}$を定義すると式1は以下のように書き換えられます。

Q_{ij}=\frac{D_{ij}}{L_{ij}} (p_i-p_j)\tag{3}

$D$は単位長さ当たりの流れやすさ(Conductivity)を意味します。($\eta$がエッジによらず一定であると仮定すれば)$D$の値はエッジがどれくらい太いかということを示します。$D$が大きければそのエッジは太く、小さければ細いです。$D$は時間がたつにつれて変化します。その変化量は以下の式に基づいて導出されます。

\frac{dD_{ij}}{dt}=f(|Q_{ij}|)-D_{ij}\tag{4}

右辺第1項は、あるエッジの流れやすさはそのエッジを流れる栄養の流量に応じて変化することを意味し、右辺第2項はエッジが細くなる速度を意味します。$f(|Q|)$の中身を

f(|Q|)=\frac{|Q|^\gamma}{1+|Q|^\gamma}\tag{5}

とおきます。式4,5から栄養の流量が多くなるほどそのエッジは太くなりさらに流れやすさがアップし、流量が少なく使われていないエッジは自ら細くなっていき栄養がより流れにくくなることがわかります。

エッジの長さや流れやすさは流れの向きによらないものとします。式3より$Q$は$i$,$j$を入れ替えると符号は反転します。

L_{ij}=L_{ji}\\
D_{ij}=D_{ji}\\
Q_{ij}=-Q_{ji}\tag{6}

ノードの中から餌場所となるノードを選択します。この説明では$N_1$を餌ノード(Source)、$N_2$をもうひとつの餌ノード(Sink)とします。$N_1$では単位時間あたり$I_0$の栄養が供給されます。$N_2$では単位時間当たり$I_0$の栄養が消費されます(=マイナスの栄養が供給される?)。これ以外から栄養の出入りはありません。キルヒホッフの法則(各ノードにおいて流入の総和と流出の総和は等しいという法則)を考えると、

\sum_j Q_{1j} = +I_0 (i=1)\\
\sum_j Q_{2j} = -I_0 (i=2)\\
\sum_j Q_{ij} = 0 (i\neq1,2)\\
\tag{7}

となります。この式は「単位時間当たりにある一つのノードについて他のノードから流入した分と他のノードへ流出した分の合計をとったとき、そのノードが$N_1$なら$I_0$、$N_2$なら$-I_0$、それ以外のノードなら0になっているよ」ということを言っています。

アルゴリズム

数理モデルのルールを一通り学んだところで、次は実際にどのようなステップを踏んで計算を進めていくかについて簡単なグラフを使いながら説明します。

step0.使うグラフ

図のような5個のノードからなる無向グラフについて考えます。

nenkin02.png

import networkx as nx
import matplotlib.pyplot as plt
import numpy as np

G=nx.Graph()
G.add_nodes_from([0,1,2,3,4])
G.add_edges_from([(0,2),(0,3),(2,3),(2,4),(1,4),(1,3),(3,4)])
pos = nx.spring_layout(G)
nx.draw_networkx(G,pos)
plt.axis('off')
plt.show()

このグラフは隣接行列を用いて以下のように書けます。隣接行列はノードどうしがつながっているかどうかを示します。つながっていれば1、つながっていなければ0です。隣接行列の一番上、左から3番目(0行2列目、Pythonではインデックスが0から始まることに注意してください)は「1」ですからノード0とノード2はつながっています。

\left(
\begin{matrix}
0&0&1&1&0 \\
0&0&0&1&1 \\
1&0&0&1&1 \\
1&1&1&0&1 \\
0&1&1&1&0
\end{matrix}
\right)\tag{8}
nbr_mtx=nx.to_numpy_matrix(G)

step.1 使う変数

各エッジの流れやすさ、長さ、各ノードの圧力、各エッジの流量をそれぞれ以下のように定義します。

\boldsymbol{D}=\left(
\begin{matrix}
0&0&D_{02}&D_{03}&0 \\
0&0&0&D_{13}&D_{14} \\
D_{20}&0&0&D_{23}&D_{24} \\
D_{30}&D_{31}&D_{32}&0&D_{34} \\
0&D_{41}&D_{42}&D_{43}&0
\end{matrix}
\right), \\
\quad\\
\boldsymbol{L}=\left(
\begin{matrix}
0&0&L_{02}&L_{03}&0 \\
0&0&0&L_{13}&L_{14} \\
L_{20}&0&0&L_{23}&L_{24} \\
L_{30}&L_{31}&L_{32}&0&L_{34} \\
0&L_{41}&L_{42}&L_{43}&0
\end{matrix}
\right), \\
\quad\\
\boldsymbol{p}=\left(
\begin{matrix}
p_0\\
p_1\\
p_2\\
p_3\\
p_4\\
\end{matrix}
\right), \\
\boldsymbol{Q}=\frac{\boldsymbol{D}}{\boldsymbol{L}}(\boldsymbol{p}^T-\boldsymbol{p})\tag{9}

step.2 一次連立方程式を立てて解き、圧力を求める

さて、ノード0が餌ノード($+I_0$)、ノード1をもうひとつの餌ノード($-I_0$)とおくと、式7より

\begin{align}
Q_{02}+Q_{03}&=I_0\\
Q_{13}+Q_{14}&=-I_0\\
Q_{20}+Q_{23}+Q_{24}&=0\\
Q_{30}+Q_{31}+Q_{32}+Q_{34}&=0\\
Q_{41}+Q_{42}+Q_{43}&=0
\end{align}\tag{10}

が成立します。式3からこの一次連立方程式は以下の行列形式に変換できます。

\left(
\begin{matrix}
\frac{D_{02}}{L_{02}}+\frac{D_{03}}{L_{03}}&
0&
-\frac{D_{02}}{L_{02}}&
-\frac{D_{03}}{L_{03}} &
0\\
0&
\frac{D_{13}}{L_{13}}+\frac{D_{14}}{L_{14}}&
0&
-\frac{D_{13}}{L_{13}}&
-\frac{D_{14}}{L_{14}} \\
-\frac{D_{20}}{L_{20}}&
0&
\frac{D_{20}}{L_{20}}+\frac{D_{23}}{L_{23}}+\frac{D_{24}}{L_{24}}&
-\frac{D_{23}}{L_{23}}&
-\frac{D_{24}}{L_{24}}\\
-\frac{D_{30}}{L_{30}}&
-\frac{D_{31}}{L_{31}}&
-\frac{D_{32}}{L_{32}}&
\frac{D_{30}}{L_{30}}+\frac{D_{31}}{L_{31}}+\frac{D_{32}}{L_{32}}+\frac{D_{34}}{L_{34}}&
-\frac{D_{34}}{L_{34}}\\
0&
-\frac{D_{41}}{L_{41}} &
-\frac{D_{42}}{L_{42}}&
-\frac{D_{43}}{L_{43}}&
\frac{D_{41}}{L_{41}}+\frac{D_{42}}{L_{42}}+\frac{D_{43}}{L_{43}}
\end{matrix}
\right)
\left(
\begin{matrix}
p_0\\p_1\\p_2\\p_3\\p_4
\end{matrix}
\right)
=
\left(
\begin{matrix}
I_0\\-I_0\\0\\0\\0
\end{matrix}
\right)\tag{11}

左辺の大きな行列(以後行列$\boldsymbol{A}$と呼びます)が逆行列を持つ場合連立方程式は解を持ち、方程式を解けば$\boldsymbol{p}$が求まります。行列$\boldsymbol{A}$の各成分をよく見ると以下のような法則が成り立っていることがわかります。この法則が分かれば、自分で行列$\boldsymbol{A}$を作ることができます。

(対角成分)\sum_{j=neighbors}\frac{D_{ij}}{L_{ij}}\\
(隣接あり)-\frac{D_{ij}}{L_{ij}}\\
(隣接なし)0

step.3 圧力から次世代の流れやすさを求める

$\boldsymbol{p}$が求まれば$\boldsymbol{Q}$を求めることができ、式4,5から次世代の$\boldsymbol{D}$が計算できます。

\begin{align}
D(t+dt)&=D(t)+dD\\
&=D(t)+(f(|Q(t)|)-D(t))dt\tag{12}
\end{align}

step.4 step.2と3を繰り返す

次世代の$\boldsymbol{D}$がわかれば次世代の$\boldsymbol{p}$が求まります。次世代の$\boldsymbol{p}$がわかれば次々世代の$\boldsymbol{D}$がわかります。あとはその繰り返しです。

計算

  • Anaconda+Jupyter notebookで実行
  • Python3.6.8
  • matplotlib 3.0.3
  • numpy 1.16.3
  • networkx 2.3

$D$の非ゼロ要素の初期値は全て1.0とします。今回各エッジの長さは常に1とします。

import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import PillowWriter, FuncAnimation

G=nx.Graph()
G.add_nodes_from([0,1,2,3,4])
G.add_edges_from([(0,2),(0,3),(2,3),(2,4),(1,4),(1,3),(3,4)])
pos = nx.spring_layout(G)

#変数宣言
nt=13
#今回はdt=1とする
nnodes=nx.number_of_nodes(G)
edge_list=list(G.edges)
nbr_mtx=nx.to_numpy_matrix(G)

conductivity=np.zeros((nt,nnodes,nnodes)) #D
length=np.zeros((nnodes,nnodes)) #L
pressure=np.zeros((nt,nnodes)) #p
flux=np.zeros(nnodes) #Q

#初期値と定数の設定
conductivity[0]=nbr_mtx #各エッジの流れやすさは最初1.0
length=nbr_mtx #今回は各エッジの長さは1.0とする
I0=2.0
flux[0]=I0
flux[1]=-I0
gamma=1.8

#ゼロ除算が発生しないように非ゼロ要素のみで割り算を行う
def divide_non_zero_element(D,L,num_nodes,list_edges):
    X=np.zeros((num_nodes,num_nodes))
    for i,j in list_edges:
        X[i,j]=D[i,j]/L[i,j]
        X[j,i]=D[j,i]/L[j,i]
    return X

#f(Q) for dD/dt 
def f(x):
    powered=x**gamma
    return powered/(powered+1)

#Dの時間変化量を求める
def dD(D,L,p,num_nodes,list_edges):
    X=divide_non_zero_element(D,L,num_nodes,list_edges)
    Q=np.multiply(X,np.expand_dims(p,axis=1)-p)
    ans=f(np.abs(Q))
    return ans

#一次連立方程式を解きpを求める
def deduce_p(D,L,B,num_nodes,list_edges):
    Y=divide_non_zero_element(D,L,num_nodes,list_edges)    
    A=np.diag(np.sum(Y,axis=1))-Y
    p=np.linalg.solve(A,B)
    return p

#pの初期値を求める
pressure[0]=deduce_p(conductivity[0],length,flux,nnodes,edge_list)

#繰り返し計算する
for t in range(0,nt-1):
    conductivity[t+1]=dD(conductivity[t],length,pressure[t],nnodes,edge_list)
    pressure[t+1]=deduce_p(conductivity[t+1],length,flux,nnodes,edge_list)

描画

%matplotlib nbagg
fig=plt.figure()

datatype=[('conductivity',float)]

def animate(i):
    A=np.matrix(conductivity[i],dtype=datatype)
    G=nx.from_numpy_matrix(A)
    weights=[G[u][v]['conductivity'] for u,v in G.edges()]

    plt.cla()
    nx.draw_networkx(G,pos,width=weights)
    plt.axis('off')
    plt.title('t=' + str(i))
    #plt.show()

anim = FuncAnimation(fig,animate,frames=nt,repeat=True)
#anim.save("nenkin01.gif", writer='pillow',fps=6)
plt.show()

ezgif.com-optimize.gif

もっと複雑なグラフで試す

はじめに適当なグラフを作ります。

G0=nx.triangular_lattice_graph(10,10)
pos0=nx.spring_layout(G0)
nx.draw_networkx(G0,pos0)
plt.show()

nenkin03.png

生成したグラフはノードのキー名が(0,0),(0,1),,,になっててこのままだと今回作ったプログラムでは使えないので、一度隣接行列に落とし、その隣接行列からグラフを再生成するというまわりくどいことをします。

nbr_mtx=nx.to_numpy_matrix(G0)
G=nx.from_numpy_matrix(nbr_mtx)
pos=nx.spring_layout(G)
nx.draw_networkx(G,pos)
plt.show()

nenkin04.png

次に、餌ノードと任意に選択します。ここではノード0,35を餌ノード(Source),ノード5,60,65をもうひとつの餌ノード(Sink)としています。餌ノードが複数になっても粘菌ネットワーク全体で餌の量が保存されるように餌ノードの生産量を調整しています。

計算

import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import PillowWriter, FuncAnimation

#変数宣言
nt=12
#dt=1とする
nnodes=nx.number_of_nodes(G)
edge_list=list(G.edges)
nbr_mtx=nx.to_numpy_matrix(G)

conductivity=np.zeros((nt,nnodes,nnodes)) #D
length=np.zeros((nnodes,nnodes)) #L
pressure=np.zeros((nt,nnodes)) #p
flux=np.zeros(nnodes) #Q

#初期値と定数の設定
conductivity[0]=nbr_mtx #各エッジの流れやすさは最初1.0
length=nbr_mtx #今回は各エッジの長さは1.0とする
I0=2.0
source_list=[0,35]
sink_list=[5,60,65]
source_len=len(source_list)
sink_len=len(sink_list)

for i in source_list:
    flux[i]=I0

#餌ノードの量が全体で保存されるように調整する
for i in sink_list:
    flux[i]=-I0*(source_len/sink_len)

gamma=1.5

#ゼロ除算が発生しないように非ゼロ要素のみで割り算を行う
def divide_non_zero_element(D,L,num_nodes,list_edges):
    X=np.zeros((num_nodes,num_nodes))
    for i,j in list_edges:
        X[i,j]=D[i,j]/L[i,j]
        X[j,i]=D[j,i]/L[j,i]
    return X

#f(Q) for dD/dt 
def f(x):
    powered=x**gamma
    return powered/(powered+1)

#Dの時間変化量を求める
def dD(D,L,p,num_nodes,list_edges):
    X=divide_non_zero_element(D,L,num_nodes,list_edges)
    Q=np.multiply(X,np.expand_dims(p,axis=1)-p)
    ans=f(np.abs(Q))
    return ans

#一次連立方程式を解きpを求める
def deduce_p(D,L,B,num_nodes,list_edges):
    Y=divide_non_zero_element(D,L,num_nodes,list_edges)    
    A=np.diag(np.sum(Y,axis=1))-Y
    p=np.linalg.solve(A,B)
    return p

#pの初期値を求める
pressure[0]=deduce_p(conductivity[0],length,flux,nnodes,edge_list)

#繰り返し計算する
for t in range(0,nt-1):
    conductivity[t+1]=dD(conductivity[t],length,pressure[t],nnodes,edge_list)
    pressure[t+1]=deduce_p(conductivity[t+1],length,flux,nnodes,edge_list)

描画

%matplotlib nbagg
fig=plt.figure()

datatype=[('conductivity',float)]

def animate(i):
    A=np.matrix(conductivity[i],dtype=datatype)
    G=nx.from_numpy_matrix(A)
    weights=[G[u][v]['conductivity'] for u,v in G.edges()]

    plt.cla()
    nx.draw_networkx(G,pos,width=weights,node_size=100)
    nx.draw_networkx_nodes(G,pos,nodelist=source_list,node_color='red')
    nx.draw_networkx_nodes(G,pos,nodelist=sink_list,node_color='purple')
    plt.axis('off')
    plt.title('t=' + str(i))
    #plt.show()

anim = FuncAnimation(fig,animate,interval=500,frames=nt,repeat=True)
#anim.save("nenkin02.gif", writer='pillow',fps=2)
plt.show()

ezgif.com-optimize (2).gif

粘菌数理モデルによる経路探索。赤色のノードは餌ノード(Source)、紫色のノードはもう一つの餌ノード(Sink)、青色のノードはそれ以外であることを示しています。$\gamma$の値が小さいと比較的輸送効率の悪い経路に対して寛容になります。

課題

計算回数を重ねるとLinAlgError: Singular matrixがでて止まります。おそらく経路にあまり関係のないエッジで成分がゼロに近づきすぎ、行列$\boldsymbol{A}$がランク落ちして逆行列を持たなくなってしまうことが原因と考えられます。もう使いそうにないエッジを切り捨てるとか改良が必要です。

今後やりたいこと

高速化(疎行列にしてみる、numbaとかscipyとか使う)
標高差の因子をいれたモデルを考える(低い位置には流れやすく高い位置には流れにくい)。

gif動画がループしないことについて

出力されたgif動画を他のプログラムで再生するとループせず最後のフレームで止まるのではないでしょうか。それはあなたのせいではありません。matplotlibのpillowWriterのバグです。matplotlib v3.1では修正される予定です。
https://github.com/t-makaro/animatplot/issues/24
ループするgif動画を保存したい場合はimagemagickをインストールしwriter='imagemagick`を使用してください。

参考文献

  • A. Tero, R. Kobayashi, and T. Nakagaki, "Physarum solver - a biologically inspired method of road-network navigation -", Physica A Statistical Mechanics and its Applications 363. 115(2006).

  • A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi, and T. Nakagaki, "Rules for biologically inspired adaptive network design", Science 327, 439(2010).

  • 加藤由人、大谷紀子、”粘菌ネットワークを用いた時間枠制約付き巡回セールスマン問題の解法”、情報処理学会第77回全国大会、4T-03(2015).

  • Python を使ってお手軽に数学 gif を作ってみよう

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした