Adjust solution obtained from DifferentialEquations.jl

[Luigi_Marongiu]

Hello,
I am using DifferentialEquations.jl to solve these ordinary differential equations:

where dN/dt is bacteria susceptible, dI/dt is bacteria infected and dV/dt is phages. I wrote the following code:

function dynamo!(du, u, p, t)
    μ, κ, φ, ω, η, β = 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

# set parameters
mu    = 0.16        # maximum growth rate
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate
beta  = 50.0        # burst size
omega = 0.05        # outflow
eta   = 1.0         # lyse rate
tmax  = 4000.0      # time span 0-tmax
s0    = 50000.0     # initial susceptible population
i0    = 0.0         # initial infected population
v0    = 80.0        # initial phage population

u0 = [s0, i0, v0]
p = [mu, kappa, phi, omega, eta, beta]
tspan = (0.0, tmax)

# set solver
prob = ODEProblem(dynamo!, u0, tspan, p)
soln = solve(prob)

and it works, but the plot I get is this:

instead of this:

I guess the difference is in the algorithm chosen by the solver (unless there are flaws in the conversion of the equations).
Can I tune the solver in order to get the expected results? and if I change the algorithm, won’t I introduce a bias in the results, for they will be what they should be but what they are supposed to be? In this case, I have a certain solution available, but what for real life data?
Thank you

[BLI]

It is a little confusing that you use different symbols in your code compared to the model picture you include.

Anyway, observe that your plot (I don’t know which package you used) uses linear scaling of the ordinate axis (“y axis”), while the plot you compare with uses logarithmic scaling of the ordinate axis.

If you use the Plots package, you can change to logarithmic scaling by adding keyword yscale = :log10 in the plot function (if my memory serves me well).

[Luigi_Marongiu]

I though adjustment of the scale was automatic, sorry. But you were right: it was a problem of scale. However, with yscale = :log10 I got

With ax = gca(); ax.set_yscale("log") I got:

Which is consistent with the reference figure. I don’t know why the colours of the lines are different from those of the legend though. Thank you

[BLI]

I’ve modified your code slightly… First, I use the following packages:

using Plots, LaTeXStrings; pyplot()
using DifferentialEquations

Next, I leave your function as it is, but change 3 of the following values:

# set parameters
mu    = 0.16        # maximum growth rate
kappa = 2.2e7    # maximum population density
phi   = 1e-9     # adsorption rate
beta  = 50.0        # burst size
omega = 0.05        # outflow
eta   = 1.0         # lyse rate
tmax  = 4000.0      # time span 0-tmax
s0    = 50000.0     # initial susceptible population
i0    = 1e-6         # initial infected population
v0    = 80.0        # initial phage population

I’ve changed the following:

1. kappa = 2.2*10^7 → kappa = 2.2e7. The latter is the proper way to write 2.2\cdot 10^7. Your method is potentially dangerous for the following reason: 10^7 is treated as an integer in Julia. When you multiply floating point value 2.2 by integer 10^7 you do get a floating point number. So, in this case: no problem. But consider that 10^18 becomes 1000000000000000000 =10^{18}, while 10^19 becomes -8446744073709551616 \neq 10^{19}… what happens here is that 10^{19} is larger than the largest integer that Julia normally handles, and what you see is integer overflow.
2. In general, xen where x is an integer or * floating point* number, n is a positive or negative integer, and e is … e – creates a floating point number x\cdot 10^n. xen is the proper Julia way to write x\cdot 10^n, while x*10^n is potentially dangerous.
3. I have thus also changed phi = 10.0^-7 → phi = 1e-7.
4. I have also changed i0 = 0 → i0 = 1e-6. The reason is that with i0 = 0, you get an error message with logarithmic scaling because the logarithm of 0 obviously is -\infty.

With the Plots package and LaTeXStrings, plotting can be done as follows:

plot(soln,label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9))
plot!(xlabel="time [h]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

leading to:

Personally, I’d have scaled the time axis to make the plot easier to “read”. I don’t know how to do that with the data structure of soln (Maybe Chris R has a trick up his sleeve). One way to do it is to take out the solutions into a matrix:

Soln = reduce(hcat,soln.u)';

and plot the solution as follows:

plot(soln.t/24/7,Soln,lw=2,label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9))
plot!(xlabel="time [week]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

leading to:

The key advantage of the first way to plot the solution is that the solution soln includes an interpolation function which is used to generate a smooth plot, while in the second approach (Soln) plot uses linear interpolation between the data points. (Perhaps not so easy to see here…)

[BLI]

The trick up the sleeve… how to scale the time axis (or other things) while keeping the DifferentialEquations interpolator for the solution…

So – keep the soln, i.e., solution as described above, and forget about the Soln matrix unpacking.

First, I pick out which variables in the solution I want to plot:

plot(soln,vars=[(0,1),(0,2),(0,3)],label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9))
plot!(xlabel="time [h]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

Here, this gives the same plot as above (plot(soln, label=...)). But I have specified that the first plot involves the mapping s_0 \rightarrow s_1 where s is my shorthand notation for soln and s_0 = time (independent variable) while s_1 is dependent variable 1; indicated by tuple (0,1). Next, (0,2) is s_0 \rightarrow s_2, while (0,3) is s_0 \rightarrow s_3. Thus, vars=[(0,1), (0,2), (0,3)] implies to plot each variable (1,2,3) as a function of time (0).

But I can also include a mapping of the plots. Suppose I create a function for scaling hour to week. Thus:

t_h2w(t,u) = (t/24/7,u)

takes a pair of time and solution ((t,u)) and returns the same pair with the first argument (time) scaled from hour to week.

I then inform the plot command about this transformation as follows:

plot(soln,vars=[(t_h2w,0,1),(t_h2w,0,2),(t_h2w,0,3)],label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9))
plot!(xlabel="time [w]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

With this strategy, I get a weird result…

So… what happened is that plots changed the time axis to weeks (or rather: scaled t \rightarrow t/24/7, but did not change the time axis limit. To fix this, I simply need to do as follows:

 plot(soln,vars=[(t_h2w,0,1),(t_h2w,0,2),(t_h2w,0,3)],label=["Susceptible bacteria" "Infected bacteria" "Bacteriophage"])
plot!(yscale=:log10,ylim=(1e2,1e9),xlim=(0,tmax/24/7))
plot!(xlabel="time [w]", ylabel="concentration [1/mL]",title=L"\eta = 1.0")

and presto!

The advantage of this approach over my previous suggestion of extracting data into matrix Soln is that we have used the interpolating functions of soln, thus getting a smoother solution.