Constraints in a stochastic differential equation problem

Specific Domains Modelling & Simulations
,
1

I am trying to implement constraints into a spatial chemical langevin equation. I have an array where every two indices corresponds to a spatial location. Every odd/even index corresponds to the number of agents of type w/m at the given location. I would like to restrict the number of agents of types w and m to be between 0 and K at every location, where K is some number called the carrying capacity. I have figured out how to simulate my model (The code is below), but have not yet figured out how to implement constraints at every location in my array. How does one do this?

# CLE Parameters/Set-Up:
K = 100.0 # Carrying capacity
M = 1000 # Number of locations - boundaries
T = 10.0 # Final time
dt = 1e-9 # Time-step 
r_w = 0.1 # Wild-type growth rate
r_m = 0.2 # Mutant growth rate
N = 10 # Number of trajectories

S_init = repeat([K, 0.0],M+2) # Initial wave (esentially a down-step)
S_init[M+2:end] .= 0.0 # ... initial wave setup continued

# Deterministic Portion
function ϕ(dS,S,p,t)
    dS[1] = p[1] * S[1] * (p[3] - S[1] - S[2]) + p[3] * (S[3] - S[1])
    dS[2] = p[2] * S[2] * (p[3] - S[1] - S[2]) + p[3] .* (S[4] - S[2])
    dS[3:2:end-3] .= p[1].*S[3:2:end-3].*(p[3] .- S[3:2:end-3] .- S[4:2:end-2]) .+ p[3] .* (S[1:2:end-5] .- 2.0 .* S[3:2:end-3] .+ S[5:2:end-1])
    dS[4:2:end-2] .= p[2].*S[4:2:end-2].*(p[3] .- S[3:2:end-3] .- S[4:2:end-2]) .+ p[3] .* (S[2:2:end-4] .- 2.0 .* S[4:2:end-2] .+ S[6:2:end])
    dS[end-1] = p[1] * S[end-1] * (p[3] - S[end-1] - S[end]) + p[3] * (S[end-3] - S[end-1])
    dS[end] = p[2] * S[end] * (p[3] - S[end-1] - S[end]) + p[3] * (S[end-2] - S[end])
end

# Noise Term
function ξ(dS,S,p,t)
    dS[1] = sqrt(clamp((2.0 - p[1]) * S[1] * (p[3] - S[1] - S[2]) + 2.0 * S[1] * (S[2] - S[3]) + p[3] * (S[1] + S[3]), 0.0, p[3]^2))
    dS[2] = sqrt(clamp((2.0 - p[2]) * S[2] * (p[3] - S[1] - S[2]) + 2.0 * S[2] * (S[1] - S[4]) + p[3] * (S[2] + S[4]), 0.0, p[3]^2))
    dS[3:2:end-3] .= sqrt.(clamp.(6.0 .* S[3:2:end-3] .* (p[3] .- S[3:2:end-3]) .+ (p[3] .- 2.0 .* S[3:2:end-3]) .* (S[1:2:end-5] .- 2.0 .* S[3:2:end-3] .+ S[5:2:end-1]) .- p[1].*S[3:2:end-3].*(p[3] .- S[3:2:end-3] .- S[4:2:end-2]), 0.0, p[3]^2))
    dS[4:2:end-2] .= sqrt.(clamp.(6.0 .* S[4:2:end-2] .* (p[3] .- S[4:2:end-2]) .+ (p[3] .- 2.0 .* S[4:2:end-2]) .* (S[2:2:end-4] .- 2.0 .* S[4:2:end-2] .+ S[6:2:end]) .- p[2].*S[4:2:end-2].*(p[3] .- S[3:2:end-3] .- S[4:2:end-2]), 0.0, p[3]^2))
    dS[end-1] = sqrt(clamp((2.0 - p[1]) * S[end-1] * (p[3] - S[end-1] - S[end]) + 2.0 * S[end-1] * (S[end] - S[end-3]) + p[3] * (S[end-1] + S[end-3]), 0.0, p[3]^2))
    dS[end] = sqrt(clamp((2.0 - p[2]) * S[end] * (p[3] - S[end-1] - S[end]) + 2.0 * S[end] * (S[end-1] - S[end-2]) + p[3] * (S[end] + S[end-2]), 0.0, p[3]^2))
end

prob_sde = SDEProblem(ϕ,ξ,S_init,(0.0,T),[r_w,r_m,K]) # params p = [r_w,r_m,K])
ensemble_prob_sde = EnsembleProblem(prob_sde)
sol = solve(ensemble_prob_sde,RKMil(),EnsembleThreads(),trajectories=N,dt=dt)
2

You’d have to design the equation so the constraint is by design. Otherwise it’s not an SDE.

1 Like