How to set a systems of ODEs to zero with Boundary Value Problem?

Specific Domains Optimization (Mathematical)
, , ,
1

I have a system of ordinary differential equations (ODEs) that model the interaction between a bacterium and phage:

dN/dt = μN(1-(N+I)/κ) -φNP -ωN
dI/dt = φNP -ηI -ωI
dP/dt = βηI -φNP -ωP

I am looking for the values of P that bring N to zero in a stable term. As is, the addition of P to the system reduces N, but then N and P go into oscillation and then they reach an equilibrium where N >> 0.

treatment_1
treatment_11920×1440 111 KB

I understand this might be a problem for Boundary Value Problem but how do I set the system?
I understand BVP requires a function to work upon, which probably is the dN/dt I indicated, hence:

function fun!(du, u, p, t) 
    μ, κ, φ, ω = p
    #=
    du[1] = susceptible
    =#
    du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * u[1] * u[3]) - (ω * u[1])
end

(do I need to include dI/dt?) This is derived from the system of ODEs I set for the model:

function modelFun!(du, u, p, t) # base infection model with logistic term
    μ, κ, φ, ω, η, β = p
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = phages
    =#
    du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * u[1] * u[3]) - (ω * u[1])
    du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
    du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])
end

I then need an implicit function bc, but what would that be?
How is the setting of this kind of problem?
The parameters I used are:

κ = 2.2 ⋅ 10^7
μ = 0.47
φ = 10^-9
β = 50
η = 1.0
ω = 0.05
N0 = 1.0*10^7
P0 = 2*10^8

Thank you

2

I tried with:

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
  ...