# Patterning

In this model we implement the paper from Corson et al. (2017).

using CellBasedModels
using GLMakie #Can be changes to CairoMakie
using Distributions
Makie.inline!(true)
true

## Create model

First we define some functions that we will use for then model

fσ(x) = (1+tanh(2*x))/2
fs0(x,t,l,S0,τg,L) = S0*fσ(1-t/τg)*(exp(-x^2/(2*L^2))+exp(-(1-x)^2/(2*L^2))) + fσ(t/τg-1)*(exp(-x^2/(2*l^2))+exp(-(1-x)^2/(2*l^2)))
fs0 (generic function with 1 method)

We create the ABM model.

model = ABM(2,

agent = Dict(
:s0 => Float64,
:u => Float64,
:s => Float64
),

model = Dict(
:a0=>Float64,
:a1=>Float64,
:τ=>Float64,
:l=>Float64,
:D=>Float64,
:S0=>Float64,
:τg=>Float64,
:L=>Float64
),

agentODE = quote

s0 = fs0(x,t,l,S0,τg,L)

s = 0
@loopOverNeighbors i2 begin
d = minimum(
[
(x-x[i2])^2+(y-y[i2])^2,
(x-x[i2]+1)^2+(y-y[i2])^2,
(x-x[i2]-1)^2+(y-y[i2])^2,
(x-x[i2])^2+(y-y[i2]+1)^2,
(x-x[i2])^2+(y-y[i2]-1)^2,
(x-x[i2]+1)^2+(y-y[i2]+1)^2,
(x-x[i2]+1)^2+(y-y[i2]-1)^2,
(x-x[i2]-1)^2+(y-y[i2]+1)^2,
(x-x[i2]-1)^2+(y-y[i2]-1)^2
]
)
s += exp(-d/(2*l^2)) * u[i2]*(a0 + 3*u[i2]^3*a1/(1+u[i2]^2))
end

dt( u ) = fσ(2*(u-s-s0))/τ -u/τ

end,

agentSDE = quote
dt(u) = D
end,

agentAlg=DifferentialEquations.EM()
);

## Initialize the community

Lx = 1
Ly = 1
Nx = 18
Ny = 9

com = Community(model,
N=2*Nx*Ny,
dt=0.001,
)

#Global parameters
λ = 5*10^-6; a0 = .05; a1 = 1 - a0; τ = 1/2; l = 0.085#1.75*λ;
D = 5*10^-5; S0 = 2; τg = 1; L = .2; N = 324
com[:a0] = a0
com[:a1] = a1
com[:τ] = τ
com[:τg] = τg
com[:l] = l
com[:D] = D
com[:s0] = S0
com[:L] = L

#Positions
dist = Uniform(-1,1)
posx = zeros(2*Nx*Ny); posy = zeros(2*Nx*Ny)
for i in 1:Nx
for j in 1:Ny
posx[Nx*(j-1)+i] = Lx*(i-0.5)/Nx +0.01*rand(dist)
posy[Nx*(j-1)+i] = Ly*(j-0.5)/Ny +0.01*rand(dist)
end
end
for i in 1:Nx
for j in 1:Ny
posx[Nx*Ny+Nx*(j-1)+i] = Lx*(i)/Nx +0.01*rand(dist)
posy[Nx*Ny+Nx*(j-1)+i] = Ly*(j-0.5+cos(pi/3))/Ny +0.01*rand(dist)
end
end

com[:x] = posx; com[:y] = posy;
#Concentration
u0 = [fs0(i,0.,l,S0,τg,L) for i in posx]
com[:u] = u0;
com[:s0] = u0;

## Evolution

evolve!(com,steps=4000,saveEach=10,saveCurrentState=true);

## Plotting results

comOut = getParameter(com,[:u])

fig = Figure(resolution=(1500,300))

for (i,time) in enumerate(1:round(Int64,length(com)/4):length(com))
ax = Axis(fig[1,i])
meshscatter!(ax,com[:x],com[:y],markersize=3*10^-2,color=comOut[:u][time])
end

display(fig)