Initialize Celani at steady state

src/chemotaxis/celani.jl

@@ -13,7 +13,7 @@ tuned through a `chemotactic_precision` factor inspired by
 Default parameters:
 - `motility = RunTumble(speed = [30.0])`
 - `turn_rate = 1.49` Hz
-- `state = 1`
+- `state = 0`
 - `rotational_diffusivity = 0.26` rad²/s
 - `gain = 50.0`
 - `memory = 1` s
@@ -42,7 +42,7 @@ mutable struct Celani{D} <: AbstractMicrobe{D}
         speed::Real = rand_speed(motility),
         turn_rate::Real = 1/0.67, # 1/s
         markovian_variables::Vector{<:Real} = zeros(3),
-        state::Real = 1.0,
+        state::Real = 0.0,
         rotational_diffusivity::Real = 0.26, # rad²/s
         gain::Real = 50.0, # 1
         memory::Real = 1.0, # s
@@ -67,11 +67,10 @@ function _affect!(microbe::Celani, model)
     M = rand(Normal(c,σ)) # measurement
     λ = 1/microbe.memory
     W = microbe.markovian_variables
-    S = microbe.state
     W[1] += (-λ*W[1] + M)*Δt
     W[2] += (-λ*W[2] + W[1])*Δt
     W[3] += (-λ*W[3] + 2*W[2])*Δt
-    S = λ^2*(W[2] - λ*W[3]/2)
+    microbe.state = λ^2*(W[2] - λ*W[3]/2)
     return nothing
 end # function
 
@@ -80,4 +79,39 @@ function _turnrate(microbe::Celani, model)
     β = microbe.gain
     S = microbe.state
     return ν₀*(1-β*S) # modulated turn rate
-end # function
+end # function
+
+# Celani requires a custom add_agent! method
+# to initialize the markovian variables at steady state
+# depending the concentration field at the agent position
+function Agents.add_agent!(
+    pos::Agents.ValidPos,
+    A::Type{Celani{D}},
+    model::AgentBasedModel,
+    properties...;
+    vel = nothing,
+    speed = nothing,
+    kwargs...
+) where D
+    id = nextid(model)
+    microbe = A(id, pos, properties...; vel=ntuple(zero,D), speed=0, kwargs...)
+    microbe.vel = isnothing(vel) ? rand_vel(model.rng, D) : vel
+    microbe.speed = isnothing(speed) ? rand_speed(model.rng, microbe.motility) : speed
+    initialize_markovian_variables!(microbe, model)
+    add_agent_pos!(microbe, model)
+end
+
+"""
+    initialize_markovian_variables!(microbe::Celani, model)
+Initialize the internal markovian variables of `microbe` to be at steady state
+with respect to the `concentration_field` defined by `model`.
+"""
+function initialize_markovian_variables!(microbe::Celani, model)
+    W = microbe.markovian_variables
+    λ = 1/microbe.memory
+    M = model.concentration_field(microbe.pos, model)
+    W[1] = M/λ
+    W[2] = W[1]/λ
+    W[3] = 2W[2]/λ
+    nothing
+end

test/microbe-creation.jl

@@ -111,7 +111,21 @@ using LinearAlgebra: norm
             affect!(m, model)
             # concentration_field is null, so state should be unchanged
             @test m.markovian_variables == zeros(3)
-            @test m.state == 1.0
+            @test m.state == 0.0
+
+            # if the model has a non-zero concentration field
+            # markovian variables and internal state should be
+            # initialized from the appropriate steady state values
+            C = 2.0
+            concentration_field(pos, model) = C
+            properties = Dict(:concentration_field => concentration_field)
+            model = StandardABM(Celani{D}, ntuple(_->100,D), 0.1; properties)
+            add_agent!(model)
+            m = model[1]
+            p = m.pos
+            λ = 1/m.memory
+            @test m.state == 0.0
+            @test m.markovian_variables == [C/λ, C/λ^2, 2C/λ^3]
         end
     end
     @testset "Xie" begin