モルフォゲンの濃度勾配によるパターン形成

* モルフォゲンの産生、拡散、分解によって濃度勾配がつくられることをシミュレートする。
* 離散的
* 陽解法
* Julia

## 考える系
一次元のモルフォゲン濃度分布を考えます。位置$x$、時刻$t$におけるモルフォゲン濃度を$u(x,t)$とおき、空間を$dx$、時刻を$dt$で区切ることで離散化します。濃度空間分布の長さを1とすれば、区切られた領域の数は$1/dx$です。

## 初期濃度

```math
u_0=[0,0,0,…0,0,0]
```

## 拡散項
拡散速度は、ある領域と隣の領域との濃度勾配および接触面積に比例します。

```math
\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産生されます。

```math
u_p=[1,0,0…0,0,0]
```

### 境界条件
このモデルでは、右端と左端では外に物質のやり取りがないと考えます(zero-flux condition)。その場合、$u(x-dx),u(x+dx)$はイテレータを使って下図のように計算されます。

![SODmodel01.jpg](https://qiita-image-store.s3.amazonaws.com/0/253018/a7dee7cb-8340-e970-5971-a60ed37bdab0.jpeg)


### 計算
```julia
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)
![diff01.gif](https://qiita-image-store.s3.amazonaws.com/0/253018/b7464248-af45-a140-cd37-8c2456ff29c1.gif)

モルフォゲンは増え続け定常状態に落ち着いてくれません。


## Source-Sinkモデル
左端の領域では濃度が常に1になるように産生が行われ(source)、右端の領域では分解によって濃度は常に0である(sink)という固定境界条件を置きます。

### 計算

```julia
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)
```

![diff02.gif](https://qiita-image-store.s3.amazonaws.com/0/253018/0ca8697b-9a56-7f13-7adf-08c9b1cf8827.gif)

200回ほど計算したところで一次関数の形をした定常状態になるのが分かります。

## SDDモデル
はじめのモデルにsinkを付け加える代わりに全ての領域においてモルフォゲンは自然分解する(degradation)と仮定します。分解速度degはその領域における濃度に比例します。

### 計算

```julia
#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)
```
![diff03.gif](https://qiita-image-store.s3.amazonaws.com/0/253018/1ff96e5b-3706-87e9-40eb-86cf8a1f244f.gif)
指数関数的な定常状態に到達しました。この定常状態の式を求めるには微分方程式を解く必要があります。

## 参考文献
 * [三浦岳、発生の数理、京都大学学術出版会、pp35-53、2015](
https://www.amazon.co.jp/%E7%99%BA%E7%94%9F%E3%81%AE%E6%95%B0%E7%90%86-%E4%B8%89%E6%B5%A6-%E5%B2%B3/dp/4876988870)