-
Notifications
You must be signed in to change notification settings - Fork 0
/
brumley.jl
69 lines (64 loc) · 2.35 KB
/
brumley.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
export MicrobeBrumley, brumley_affect!, brumley_turnrate
"""
MicrobeBrumley{D} <: AbstractMicrobe{D}
Model of chemotactic bacterium from 'Brumley et al. (2019) PNAS'.
The model is optimized for simulation of marine bacteria and accounts
for the presence of (gaussian) sensing noise in the chemotactic pathway.
Default parameters:
- motility = RunReverseFlick(speed = Degenerate(46.5))
- turn_rate = 2.22 Hz → '1/τ₀'
- state = 0.0 → 'S'
- rotational_diffusivity = 0.035 rad²/S
- adaptation_time = 1.3 s → 'τₘ'
- receptor_gain = 50.0 μM⁻¹ → 'κ'
- motor_gain = 50.0 → 'Γ'
- chemotactic_precision = 6.0 → 'Π'
- radius = 0.5 μm → 'a'
"""
Base.@kwdef mutable struct MicrobeBrumley{D} <: AbstractMicrobe{D}
id::Int
pos::NTuple{D,Float64} = ntuple(zero, D)
motility = RunReverseFlick(speed = Degenerate(46.5))
vel::NTuple{D,Float64} = rand_vel(D) .* rand(motility.speed)
turn_rate::Float64 = 1/0.45 # 1/s
state::Float64 = 0.0
rotational_diffusivity::Float64 = 0.035 # rad²/s
adaptation_time::Float64 = 1.3 # s
receptor_gain::Float64 = 50.0 # 1/μM
motor_gain::Float64 = 50.0 # 1
chemotactic_precision::Float64 = 6.0 # 1
radius::Float64 = 0.5 # μm
end # struct
function brumley_affect!(microbe, model)
Δt = model.timestep
Dc = model.compound_diffusivity
τₘ = microbe.adaptation_time
α = exp(-Δt/τₘ) # memory persistence factor
a = microbe.radius
Π = microbe.chemotactic_precision
κ = microbe.receptor_gain
u = model.concentration_field(microbe.pos, model)
∇u = model.concentration_gradient(microbe.pos, model)
∂ₜu = model.concentration_time_derivative(microbe.pos, model)
# gradient measurement
μ = dot(microbe.vel, ∇u) + ∂ₜu # mean
σ = Π * sqrt(3*u / (π*a*Dc*Δt^3)) # noise
M = rand(Normal(μ,σ)) # measurement
# update internal state
S = microbe.state
microbe.state = α*S + (1-α)*κ*τₘ*M
return nothing
end # function
function brumley_turnrate(microbe, model)
ν₀ = microbe.turn_rate # unbiased
Γ = microbe.motor_gain
S = microbe.state
return (1 + exp(-Γ*S)) * ν₀/2 # modulated turn rate
end # function
function microbe_step!(microbe::MicrobeBrumley, model)
microbe_step!(
microbe, model;
affect! = brumley_affect!,
turnrate = brumley_turnrate
)
end # function