# How to test multiple parameters with SDE?

Hello,
I would like to find the minimal amount of phages required to distinguish a bacterial population. I have run this model:

``````using DifferentialEquations
function growth!(du, u, p, t)
μ, κ, φ, ω, η, β = p
#=
du = susceptible
du = infected
du = phages
=#
du = ((μ * u) * (1 - ((u+u)/κ))) - (φ * u * u) - (ω * u)
du = (φ * u * u) - (η * u) - (ω * u)
du = (β * η * u) - (φ * u * u) - (ω * u)
end
# parameters
mu    = 0.47        # maximum growth rate susceptible (B. longum)
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate resident commensal phage
eta   = 1.0         # lyse rate resident commensal phage
beta  = 50.0        # burst size resident commensal phage
omega = 0.05        # outflow
tmax  = 4000.0      # time span 0-tmax
r_s0  = 1.0*10^7    # initial resident susceptible population u
r_i0  = 0.0         # initial resident infected population u
r_v0  = 1000.0      # initial resident phage population u
T_inf    = 1000     # time infection
amount_inf = 500         # amount infection
Vϕ = (mu/phi) * (1-(r_s0/kappa)) - omega/phi
Tϕ = T_inf
# apply
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, 0]
parms = [mu, kappa, phi, omega, eta, beta]
# inoculum
condition1(u, t, integrator) = t==T_inf
affect1!(integrator) = integrator.u += amount_inf
cb1 = DiscreteCallback(condition1, affect1!)
# extintion
condition2(u, t, integrator) = u-1
function affect2!(integrator)
integrator.u = 0
integrator.u = 0
end
cb2 = ContinuousCallback(condition2, affect2!)
condition3(u, t, integrator) = u-1
function affect3!(integrator)
integrator.u = 0
end
cb3 = ContinuousCallback(condition3, affect3!)
# run
modification = CallbackSet(cb1, cb2, cb3)
prob = ODEProblem(growth!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[T_inf])
``````

where I have a condition that simulates extinction when the number of elements falls below 1. I tried to calculate the best amount of phage with the formula `Vϕ = (mu/phi) * (1-(r_s0/kappa)) - omega/phi` which gives 2·10⁸ phages, but the value is not good because the bacterium does not go extinct:

Essentially, I would like to try several values for both `Vϕ` (amount of phages) and `Tϕ` (time of inoculum), so Parallel Ensemble seems a good tool for the job. But how do I implement it?
In the example reported on the SDE page, the Lotka-Volterra equation is modified by a function simulating noise:

``````function g(du, u, p, t)
du = p * u
du = p * u
end
``````

but in my case, I need to modify Vϕ and Tϕ in growth’s u. I also need to define the range of value to test with the additional function:

``````function prob_func(prob, I, repeat)
x = 0.3*rand(2)
remake(prob, p=[p[1:2]; x]
end
``````

but what is the role of this function (that in the example simulates random numbers for function `f` implementing the Lotka-Volterra model.
What is the syntax in my case?
Let’s say I would like to test a range of Vϕ from 1 to 10⁹ with steps of 10 millions and Tϕ=tspan with steps of 1 hour. How would I do that?
Thank you

I tried to implement the modulation manually with a for loop:

``````k = 1
I = [1.0, 1000.0, 10000.0, 100000.0, 1000000.0, 9000000.0, 9400000.0,
9450000.0, 9460000.0, 9470000.0, 9480000.0, 9490000.0,
9500000.0,
10000000.0, 80000000.0, 90000000.0, 100000000.0, 1000000000.0,
10000000000.0, 100000000000.0]
for  i in I
println(k, " > ", i)
# inoculum
condition1(u, t, integrator) = t==Tp
affect1!(integrator) = integrator.u += i
cb1 = DiscreteCallback(condition1, affect1!)
modification = CallbackSet(cb1, cb2, cb3)
prob = ODEProblem(growth!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp])
# plot
clf()
myPlot = plot(soln.t, getindex.(soln.u, 1)+getindex.(soln.u, 2),
linestyle = "-", color = "blue")
myPlot = plot(soln.t, getindex.(soln.u, 3),
linestyle = "-", color = "green")
ax = gca()
ax.set_yscale("log")
ax.set_xlim([1, tmax])
ax.set_ylim([1e0, 1e10])
legend(["Bacteria (total)", "Phages"],
loc="lower right")
xlabel("Time (h)")
ylabel("Concentration (1/mL)")
Title = string("Inoculation of ", round(i; digits = 1), " units")
title(Title)
fig = string("/home/gigiux/Documents/model/plot/treat/treat_", k, ".png")
savefig(fig,
dpi = 300, format = "png", transparent = false)
k = k+1
end
``````

and I got this breakpoint:

So by giving more than 9 490 000 phages, the bacteria go extinct. But I could proceed to check the between 9 490 000 and 9 485 000 and so forth. This approach is a very time-consuming way of proceeding.
Can I do the same (actually better) with SDE? Also because I would like to take into account also the time of inoculation…