# Growing cell population using callback

**URL:** <https://discourse.julialang.org/t/growing-cell-population-using-callback/84970>\
**Category:** Numerics\
**Tags:** diffeq\
**Created:** [July 29, 2022, 9:59am UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970 "2022-07-29T09:59:34Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [July 29, 2022, 9:59am UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/1 "2022-07-29T09:59:35Z")

</div>

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](https://diffeq.sciml.ai/stable/features/callback_functions/#Example-3:-Growing-Cell-Population)

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.

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/f/0/f0e939b3aa9897c6ada0279abf606a0fdd352127.png)

and

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/4/5/45aaa2d7abebb8707042a88dfed701d62f4bc378.png)

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

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [July 29, 2022, 11:09am UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/2 "2022-07-29T11:09:46Z")

</div>

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

> [@hbooth](#):
>
> `idx = findall(y-> y == 1,u)`

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.

---

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [July 29, 2022, 11:15am UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/3 "2022-07-29T11:15:42Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [July 29, 2022, 11:32am UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/4 "2022-07-29T11:32:30Z")

</div>

For me it seems to work with

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

```

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

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [July 29, 2022, 11:42am UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/5 "2022-07-29T11:42:03Z")

</div>

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:

```julia
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);

```

---

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [July 29, 2022, 11:57am UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/6 "2022-07-29T11:57:21Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/b/3/b3004d658b68eb85fbfd09e7985880be845c1162.png)

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](https://global.discourse-cdn.com/julialang/original/3X/5/f/5fad6eedab5cbc36a94942fe198b275c5e1405eb.png)

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

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [July 29, 2022, 12:10pm UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/7 "2022-07-29T12:10:29Z")

</div>

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…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 29, 2022, 12:30pm UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/8 "2022-07-29T12:30:03Z")

</div>

> [@hbooth](#):
>
> You can see that the concentrations drop (division) before it actually hits one on the second division. Any ideas?

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](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)

---

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [July 29, 2022, 1:05pm UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/9 "2022-07-29T13:05:15Z")

</div>

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:

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/f/1/f12e601eb834763b0b1d4a74b108131f1bc52783.png)  
and

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/b/e/beb4d7004278689e170a4338a77838951a01a710.png)

Sorry if I’m missing something simple here.

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [July 29, 2022, 2:05pm UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/10 "2022-07-29T14:05:32Z")

</div>

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:

```julia
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](https://diffeq.sciml.ai/stable/basics/plot/#Density)

---

<div class="post-metadata">

**Author:** ![hbooth](https://avatars.discourse-cdn.com/v4/letter/h/e19b73/32.png) [@hbooth](https://discourse.julialang.org/u/hbooth)\
**Post date:** [July 29, 2022, 2:11pm UTC](https://discourse.julialang.org/t/growing-cell-population-using-callback/84970/11 "2022-07-29T14:11:33Z")

</div>

![image](https://global.discourse-cdn.com/julialang/original/3X/9/3/932dcb0cab549b3925d6a5d0290b4aa201bc2846.png)

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