Sorry, I am no mathematician, so most of this is beyond my comprehension. But I have found this – the problem is to distinguish the bacteria (blue line) by adding phage (green line) in what is known as phage therapy. I calculated the minimal amount to add by the formula Vp = (mu/phi) * (1-(r_s0/kappa)) - omega/phi (which I obtained by equating dN/dt to zero end solving for P). I don’t get extinction using such amount, but If I add more and more phage, I achieve extinction:
So my question is: how do I calculate the right amount of phage and its time of delivery?
Perhaps the boundary problem is not the right tool, but the fact that if I add one billion phages I do reach extinction tells me there must be a solution. But Instead of randomly selecting values I would like to find the right amount straight away, even because the conditions will change in future models. Same thing for the time of administration. Here I subjectively chosen t=1000, but there should be a better time (specifically when bacteria are at a certain density, even if in this context the density is constant).
I tried the suggestions but:
julia> function bc!(residual, u, p, t)
residual[1] = (u[1][1] - N0)/N0 # keep the starting value constant
residual[2] = (u[1][2] - I0)/I0 # keep the starting value constant
residual[3] = u[end][1]/N0 # u[end][1] should be zero at the final timestep
end
bc! (generic function with 1 method)
julia> N0 = 1.0*10^7
1.0e7
julia> I0 = 1E5
100000.0
julia> P0 = 2*10^8
200000000
julia> pGuess = 1.0
1.0
julia> tspan = (0.0,1500.0)
(0.0, 1500.0)
julia> bvp0 = BVProblem(modelFun!, bc!, [N0,I0,P0], tspan)
BVProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 1500.0)
u0: [1.0e7, 100000.0, 2.0e8]
julia> sol0 = solve(bvp0, Shooting(Vern7()), dtmin=0.01, dt=0.5)
ERROR: MethodError: no method matching iterate(::DiffEqBase.NullParameters)
Closest candidates are:
iterate(::Base.EnvDict) at env.jl:119
iterate(::Base.EnvDict, ::Any) at env.jl:119
iterate(::JuliaInterpreter.ExprSplitter) at /home/gigiux/.julia/packages/JuliaInterpreter/LBEYA/src/construct.jl:535
...
Stacktrace:
[1] indexed_iterate(::DiffEqBase.NullParameters, ::Int64) at ./tuple.jl:84
[2] 𝔅!(::Array{Float64,1}, ::Array{Float64,1}, ::DiffEqBase.NullParameters, ::Float64) at /home/gigiux/Documents/model/util/funCollect.jl:7
[3] (::ODEFunction{true,typeof(modelFun!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing})(::Array{Float64,1}, ::Vararg{Any,N} where N) at /home/gigiux/.julia/packages/DiffEqBase/3iigH/src/diffeqfunction.jl:248
{etc}
Can I reach such a solution with boundary problems?
Thank you