# Growing cell population using callback

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")
`````` and

``````ts = range(0, stop=10, length=100)
plot(ts,map((x)->x,sol.(ts)),lw=3,
ylabel="Amount of X in Cell 1",xlabel="Time")
`````` 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

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.

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…

For me it seems to work with

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

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

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

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 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? You can see that the concentrations drop (division) before it actually hits one on the second division. Any ideas?

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…

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)

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,sol.(ts)),lw=3,
ylabel="Amount of X in Cell 1",xlabel="Time",plotdensity=1000)
``````

gives and

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

gives Sorry if I’m missing something simple here.

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,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 perfect - that works! That’s what happens when you copy the plotting code and don’t think about what its doing . Thanks again.

1 Like