# Bringing ODEs to life

This simulation actually extends an example from the DifferentialEquations.jl documentation describing a growing cell population. Here we take those underlying model dynamics and combine them with CellularPotts.jl

We begin by loading in both CellularPotts and DifferentialEquations

`using CellularPotts, DifferentialEquations`

On the CellularPotts side, we need to create a new CellPotts model which requires a CellSpace, a CellState, and a list of penalties

The space will use is a 200×200 grid that defaults to periodic boundary conditions

`space = CellSpace(200,200)`

`200×200 Periodic 4-Neighbor CellSpace{Int64,2}`

In the CellState we specify one epithelial cell with a volume of 200 pixels The cell will be positioned at the halfway point within the space.

```
state = CellState(
names = :Epithelial,
volumes = 200,
counts = 1,
positions = size(space) .÷ 2
);
```

Keeping the model simple, we'll only include an Adhesion and Volume penalty

```
penalties = [
AdhesionPenalty([0 30;
30 30]),
VolumePenalty([5]),
]
```

```
2-element Vector{Penalty}:
Adhesion
Volume
```

Now that we have all the pieces, we can generate a new CPM model.

`cpm = CellPotts(space, state, penalties);`

By default, CellularPotts models to not record states as they change overtime to increase computional speed. To have the model record past states we can toggle the appropriate keyword.

`cpm.recordHistory = true;`

## DifferentialEquations.jl setup

From the DifferentialEquations example, a theoretical protein X was created for each cell that increases linearly in time with rate parameter α

`const α = 0.3;`

As ProteinX evolves over time for each cell, the CPM model also needs to step forward in time to try and minimize its energy. To facilitate this, we can use the callback feature from DifferentialEquations.jl. Here specifically we use the `PeriodicCallback`

function which will stop the ODE solve at regular time intervals and run some other function for us (Here it will be the `ModelStep!`

function).

```
function cpmUpdate!(integrator, cpm)
ModelStep!(cpm)
end
```

`cpmUpdate! (generic function with 1 method)`

This timeScale variable controls how often the callback is triggered. Larger timescales correspond to faster cell movement.

```
timeScale = 100
pcb = PeriodicCallback(integrator -> cpmUpdate!(integrator, cpm), 1/timeScale);
```

The ODE functions are taken directly from the DifferentialEquations example. Each cell is given the following differential equation

\[ \frac{\mathrm{d} X}{\mathrm{d} t} = \alpha X\]

```
function f(du,u,p,t)
for i in eachindex(u)
du[i] = α*u[i]
end
end
```

`f (generic function with 1 method)`

Also coming from the differential equations example, this callback is tiggered whenever a cell's ProteinX is greater than 1. Basically the cell will divide when when the ProteinX concentration is too large.

```
condition(u,t,integrator) = 1-maximum(u)
function affect!(integrator,cpm)
u = integrator.u
resize!(integrator,length(u)+1)
cellID = findmax(u)[2]
Θ = rand()
u[cellID] = Θ
u[end] = 1-Θ
#Adding a call to divide the cells in the CPM
CellDivision!(cpm, cellID)
return nothing
end
```

`affect! (generic function with 1 method)`

This will instantiate the ContinuousCallback triggering cell division

`ccb = ContinuousCallback(condition,integrator -> affect!(integrator, cpm));`

To pass multiple callbacks into DifferentialEquations, we need to collect them into a set.

`callbacks = CallbackSet(pcb, ccb);`

Define the ODE model and solve

```
u0 = [0.2]
tspan = (0.0,20.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), callback=callbacks);
```

## Visualization

We can replicate the plots from the original example

`using Plots, Printf`

Plot the total cell count over time

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

Plot ProteinX dynamics for a specific cell

```
ts = range(0, stop=20, length=1000)
plot(ts, first.(sol.(ts)), lw=3,
ylabel="Amount of X in Cell 1", xlabel="Time", legend=nothing)
```

Finally, we can create an animation of the CPM to see the cells dividing. I've dropped the first few frames because the first cell takes a while to divide.

```
proteinXConc = zeros(size(space)...)
anim = @animate for t in Iterators.drop(1:cpm.step.counter,5*timeScale)
currTime = @sprintf "Time: %.2f" t/timeScale
cpmt = cpm(t)
currSol = sol((t+1)/timeScale)
#Map protein concentrations to space, set :Medium to zero
for (i,cellID) in enumerate(cpmt.space.nodeIDs)
proteinXConc[i] = iszero(cellID) ? 0.0 : currSol[cellID]
end
plt = heatmap(
proteinXConc,
c = cgrad([:grey90, :grey, :gold], [0.1, 0.6, 0.9]),
clims = (0,1),
title=currTime,
titlelocation=:left,
titlefontsize = 32)
visualize!(plt,cpmt; colorby=:none)
end
gif(anim, "BringingODEsToLife.gif", fps = 30)
```

*This page was generated using Literate.jl.*