- モルフォゲンの産生、拡散、分解によって濃度勾配がつくられることをシミュレートする。
- 離散的
- 陽解法
- Julia
考える系
一次元のモルフォゲン濃度分布を考えます。位置$x$、時刻$t$におけるモルフォゲン濃度を$u(x,t)$とおき、空間を$dx$、時刻を$dt$で区切ることで離散化します。濃度空間分布の長さを1とすれば、区切られた領域の数は$1/dx$です。
初期濃度
u_0=[0,0,0,…0,0,0]
拡散項
拡散速度は、ある領域と隣の領域との濃度勾配および接触面積に比例します。
\frac{D_u \Biggl(\frac{u(x-dx)-u(x)}{dx} \Biggl)dx+D_u \Biggl(\frac{u(x+dx)-u(x)}{dx} \Biggl)dx}{dx^2}\\
=D_u \Biggl(\frac{u(x-dx)-u(x)+u(x+dx)-u(x)}{dx^2} \Biggl)
ここで$D_u$は拡散係数。
1番目の式の分子の第1項が1つ左の領域とのやり取り、第2項が1つ右の領域とのやり取りをあらわします。
Sourceのみのモデル
産生項
左から1番目の領域を産生部(source)とします。そこでは、モルフォゲンが単位時間ごとに+1産生されます。
u_p=[1,0,0…0,0,0]
境界条件
このモデルでは、右端と左端では外に物質のやり取りがないと考えます(zero-flux condition)。その場合、$u(x-dx),u(x+dx)$はイテレータを使って下図のように計算されます。
計算
using Plots
gr()
dx=0.05
dt=0.1
u0=zeros(1/dx) #初期濃度
up=zeros(1/dx);up[1]=1.0 #産生濃度
Du=0.01 #拡散係数
function left(arr)
prepend!(collect(Iterators.take(arr,length(arr)-1)),arr[1])
end
function right(arr)
append!(collect(Iterators.drop(arr,1)),arr[end])
end
function diffusion(arr)
Du*(left(arr)-arr+right(arr)-arr)/dx/dx
end
function dudt(arr)
arr+dt*(up+diffusion(arr))
end
u=u0
anim = @animate for g=0:150
plot(u,
xticks=0:1:20,
ylim=(0.0,2.5),
yticks=0.0:0.1:2.5,
seriestype=:bar,
title="gen = $g")
u=dudt(u)
end
@time gif(anim, "tmp/diff01.gif", fps = 5)
79.689804 seconds (1.74 k allocations: 60.438 KiB)
モルフォゲンは増え続け定常状態に落ち着いてくれません。
Source-Sinkモデル
左端の領域では濃度が常に1になるように産生が行われ(source)、右端の領域では分解によって濃度は常に0である(sink)という固定境界条件を置きます。
計算
using Plots
gr()
dx=0.05
dt=0.1
u0=zeros(1/dx)
Du=0.01
#固定境界条件
function fixConc(arr)
v=arr
v[1]=1.0
v[end]=0.0
return v
end
function left(arr)
prepend!(collect(Iterators.take(arr,length(arr)-1)),arr[1])
end
function right(arr)
append!(collect(Iterators.drop(arr,1)),arr[end])
end
function diffusion(arr)
Du*(left(arr)-arr+right(arr)-arr)/dx/dx
end
function dudt(arr)
fixConc(arr+dt*diffusion(arr))
end
u=u0
anim = @animate for g=0:250
plot(u,
xticks=0:1:20,
ylim=(0.0,1.5),
yticks=0.0:0.1:1.5,
seriestype=:bar,
title="gen = $g")
u=dudt(u)
end
gif(anim, "tmp/diff02.gif", fps = 5)
200回ほど計算したところで一次関数の形をした定常状態になるのが分かります。
SDDモデル
はじめのモデルにsinkを付け加える代わりに全ての領域においてモルフォゲンは自然分解する(degradation)と仮定します。分解速度degはその領域における濃度に比例します。
計算
#Source-Diffusion-Degradation model
using Plots
gr()
dx=0.05
dt=0.1
Du=0.01
deg=0.5 #分解速度
u0=zeros(1/dx) #初期濃度
up=zeros(1/dx);up[1]=1.0 #産生濃度
function left(arr)
prepend!(collect(Iterators.take(arr,length(arr)-1)),arr[1])
end
function right(arr)
append!(collect(Iterators.drop(arr,1)),arr[end])
end
function diffusion(arr)
Du*(left(arr)-arr+right(arr)-arr)/dx/dx
end
function dudt(arr)
arr+dt*(up-deg*arr+diffusion(arr)) #分解項を追加
end
u=u0
anim = @animate for g=0:250
plot(u,
xticks=0:1:20,
ylim=(0.0,1.5),
yticks=0.0:0.1:1.5,
seriestype=:bar,
title="gen = $g")
u=dudt(u)
end
gif(anim, "tmp/diff03.gif", fps = 5)
指数関数的な定常状態に到達しました。この定常状態の式を求めるには微分方程式を解く必要があります。