Linear concentration ramp with obstacle

We replicate the chemotaxis example with the linear ramp, with the addition of an obstacle: a vertical wall that hinders the drift along the gradient direction, with only a little hole for bacteria to pass through. The wall is made of two boxes, and bacteria slide along its surfaces.

using MicrobeAgents
using Plots

@inline function concentration_field(microbe, model)
    C₀ = model.C₀
    C₁ = model.C₁
    Lx = first(spacesize(model))
    pos = position(microbe)
    concentration_field(pos,Lx,C₀,C₁)
end
@inline concentration_field(pos,Lx,C₀,C₁) = C₀ + (C₁-C₀)*pos[1]/Lx

@inline function concentration_gradient(microbe, model)
    C₀ = model.C₀
    C₁ = model.C₁
    Lx = first(spacesize(model))
    pos = position(microbe)
    concentration_gradient(pos,Lx,C₀,C₁)
end
@inline concentration_gradient(pos,Lx,C₀,C₁) = SVector{length(pos)}(i==1 ? (C₁-C₀)/Lx : 0.0 for i in eachindex(pos))

Lx, Ly = 3000, 1500 # domain size (μm)
periodic = false
space = ContinuousSpace((Lx,Ly); periodic)
Δt = 0.05 # timestep (s)

# model setup
C₀, C₁ = 0.0, 20.0 # μM
wall = [
    ObstacleBox((1400, 0), (1600, 650)),
    ObstacleBox((1400, 850), (1600, Ly)),
]
obstacle_sets = (
    Obstacles(space, wall; contact=Slide()),
    Walls(space; contact=Reflect()),
)
properties = Dict(
    :C₀ => C₀,
    :C₁ => C₁,
    :chemoattractant => ChemicalField{2}(;
        concentration_field, concentration_gradient
    ),
    :obstacles => obstacle_sets,
)
model = StandardABM(Microbe{2}, space, Δt; properties, container=Vector)
n = 20 # number of microbes
for i in 1:n
    # start all bacteria in left half of the domain
    x = rand(abmrng(model)) * 700
    y = rand(abmrng(model)) * Ly
    pos = SVector{2,Float64}(x,y)
    motility = RunTumble([30.0], 0.67, Isotropic(2))
    rotational_diffusivity = 0.1
    add_agent!(pos, model; motility, rotational_diffusivity,
        behaviors = (chemotaxis = BrownBerg(),))
end
model

T = 300 # simulation time (s)
nsteps = round(Int, T/Δt)
adata = [position]
adf, mdf = run!(model, nsteps; adata, when=5)

traj = Analysis.adf_to_matrix(adf, :position)
x = first.(traj)
y = last.(traj)

ts = unique(adf.time) .* Δt
lw = eachindex(ts) ./ length(ts) .* 3
xmesh = range(0,Lx,length=100)
ymesh = range(0,Ly,length=100)
c = concentration_field.(Iterators.product(xmesh,ymesh),Lx,C₀,C₁)
heatmap(xmesh, ymesh, c', cbar=false, c=:Greens,
    ratio=1, axis=false, grid=false, xlims=(0,Lx), ylims=(0,Ly)
)
foreach(b -> plot!(Shape(
        [b.center[1] .+ b.halfwidth[1] .* (-1, 1, 1, -1)...],
        [b.center[2] .+ b.halfwidth[2] .* (-1, -1, 1, 1)...]
    ); c=:black, lab=false), wall)
plot!(x, y, lab=false, lw=lw, lc=(1:n)', palette=palette(:acton,size(x,2)))
scatter!(x[end,:], y[end,:], lab=false, m=:c, mc=1:n, msw=0.5, ms=8)