Growing cell population using callback

Specific Domains Numerics
1

Hi everyone,

I am building a model of a growing cell population. I want to make use of callback to evaluate the ‘division’ event, similar to the example in this tutorial - Event Handling and Callback Functions · DifferentialEquations.jl

However unlike in the tutorial, I want to maintain a spatial relationship between parent and daughter cells - if u represents a 1-D array of cells and their protein concentration, then when the protein concentration of a specific cell idx hits 1, it divides and its daughter cell is placed next to it at u[idx+1], shifting all the other cells one space along (as if we called insert!(u,idx,val)), rather than at the end of the array as in the tutorial.

In addition, I want the possibility for multiple cells to be able to divide - suppose we start with two cells with equal starting concentrations, then when one reaches a protein concentration of 1 the other will have also (assuming identical protein dynamics), meaning that 2 new cells should be created rather than just 1 (as is the behaviour in the tutorial’s model).

I have tried the below code, and it doesn’t produce the expected behaviour.

const α = 0.3

function f(du,u,p,t)
    for i in 1:length(u)
        du[i] = α*u[i]
    end
end

function condition(u,t,integrator)
    1 - maximum(u) # has any cell hit 1? 
end

function affect!(integrator)
    
    u = integrator.u
    
    idx = findall(y-> y == 1,u) # find all cells which have hit concentration == 1
    
    resize!(integrator,length(u)+length(idx)) # resize solution to accomodate all new divisions

    for id in reverse(idx) # for each cell hitting 1, divide:
        for i in reverse(id:length(u)-1) # this just replicates the behaviour of insert!(u,idx+1,u[idx]/0.5) and u[idx] = 0.5*u[idx]
            if i == id + 1
                u[i+1] = u[i]
                u[i] = 0.5*u[id]
            elseif i == id
                u[i] = 0.5*u[id]
            else
                u[i+1] = u[i]
            end
        end
    end
    
    nothing
end

callback = ContinuousCallback(condition,affect!)
u0 = [0.2]
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5(),callback=callback);

plot(sol.t,map((x)->length(x),sol[:]),lw=3,
     ylabel="Number of Cells",xlabel="Time")

image

and

ts = range(0, stop=10, length=100)
plot(ts,map((x)->x[1],sol.(ts)),lw=3,
     ylabel="Amount of X in Cell 1",xlabel="Time")

image

I might have the wrong approach here - I did think about using VectorContinuousCallback where each event in out was a particular cell dividing. However, this would require the length of out to be dynamic (the same as u) which I don’t think is supported.

I am new to this part of diffeq and so any help would be appreciated! Thanks

2

I didn’t run the code, but could it be that this line

should be
idx = findall(y-> y >= 1.0, u)

The point is that the integrator might overshoot just a little bit so that it doesn’t have exact equality.
Maybe setting rootfind=SciMLBase.RightRootFind could also help.

3

Thanks for the suggestion - I just tried changing it to idx = findall(y-> y >= 1.0, u), but unfortunately that doesn’t seem to fix it. I’ll have a look at rootfind…

4

For me it seems to work with

callback = ContinuousCallback(condition,affect!, rootfind=SciMLBase.RightRootFind)

and idx = findall(y-> y >= 1.0, u)…

5

I think the insert function had a small bug too, e.g. the second cell would be sometimes set to zero instead of 50% of the parent.

This seems to work:

using OrdinaryDiffEq

function f(du,u,p,t)
    for i in 1:length(u)
        du[i] = p.α*u[i]
    end
end

function condition(u,t,integrator)
    return 1.0 - maximum(u) # has any cell hit 1? 
end

function affect!(integrator)
    
    u = integrator.u
    
    idx = findall(y-> y >= 1.0, u) # find all cells which have hit concentration == 1
    
    resize!(integrator,length(u)+length(idx)) # resize solution to accomodate all new divisions

    for id in reverse(idx) # for each cell hitting 1, divide:
        u[id+2:end] .= @view u[id+1:end-1]
        u[id] *= 0.5
        u[id+1] = u[id]
    end
    
    nothing
end

callback = ContinuousCallback(condition,affect!, rootfind=SciMLBase.RightRootFind)
u0 = [0.2]
tspan = (0.0,10.0)
p = (; α = 0.3)
prob = ODEProblem(f,u0,tspan, p)
sol = solve(prob,Tsit5(),callback=callback);
3 Likes
6

Yes, I was just about to report back - the missing link seems to be SciMLBase.RightRootFind. Also, thanks for the bug spot and the far more elegant replacement code!

This now has the following cell division behaviour

image

If we look at the concentration dynamics of the cell in position 1 (which shouldn’t move as its the starting cell - it doesn’t get shifted along) you can see that there maybe seems to be an error in the plotting?

image

You can see that the concentrations drop (division) before it actually hits one on the second division. Any ideas?

7

Can you check the actual numbers, e.g. sol[:,1]? Maybe it is just something with the plot, it could be that it is confused since the time values are the same for the solution before and after the callback. On my computer the solution does actually hit 1.0.
The default option for save_positions=(true,true) is such that it should save the result before and after the callback…

8

That’s likely just the accuracy of the points plotted with the interpolation. You’d have to increase the number of plotted points (plotdensity) to make “spikes” more exact. A case where this is even more pronounced is the Van Der Pol stiff case:

https://benchmarks.sciml.ai/html/StiffODE/VanDerPol.html

The spikes are all the same there, but there are not “enough” plotted points to always hit the same part of the spike.

As @SteffenPL says, the values should be exactly at the dot.

(Maybe we should in the interpolation always plot the solution values too, but then if the solution is a million points, it’ll plot a million + 10_000 by default instead of the 10_000 and would cause crashing in Plots. We should think about a better heuristic though)

9

The solution does hit 1.0 and so the problem as @ChrisRackauckas says seems to just be the plotting accuracy. Although, playing around with plotdensity values doesn’t seem to have an effect for me:

ts = range(0, stop=10, length=100)
plot(ts,map((x)->x[1],sol.(ts)),lw=3,
     ylabel="Amount of X in Cell 1",xlabel="Time",plotdensity=1000)

gives

image
and

ts = range(0, stop=10, length=100)
plot(ts,map((x)->x[1],sol.(ts)),lw=3,
     ylabel="Amount of X in Cell 1",xlabel="Time",plotdensity=1000000)

gives

image

Sorry if I’m missing something simple here.

10

Ah, you give explicit time points to evaluate the interpolation of the solution. So you need to change ts and the plot density doesn’t have any effect in your case.

In you case this should work:

ts = range(0, stop=10, length=1000)
plot(ts,map((x)->x[1],sol.(ts)),lw=3, ylabel="Amount of X in Cell 1",xlabel="Time")

Otherwise you would need to call the plot function directly applied to sol, like here Plot Functions · DifferentialEquations.jl

11

image

perfect - that works! That’s what happens when you copy the plotting code and don’t think about what its doing :laughing:. Thanks again.

1 Like