Edit on GitHub

EpithelialDynamics1D

Documentation for EpithelialDynamics1D.

This is a package for simulating epithelial dynamics in one dimension, supporting cell migration and cell proliferation, implementing the model in Baker et al. (2019). In this model, cells are represented as intervals between points, with the $i$th cell given by $(x_i, x_{i+1})$, $i=1,\ldots,n-1$. The node positions $x_i$ are governed by the differential equation

\[\eta\dfrac{\mathrm dx_i}{\mathrm dt} = F\left(x_i - x_{i-1}\right) - F\left(x_{i+1} - x_i\right), \quad i=2,\ldots,n-1,\]

assuming $x_1 < x_2 < \cdots < x_n$. If the left boundary is fixed, then $\mathrm dx_1/\mathrm dt = 0$, otherwise $\eta\mathrm dx_1/\mathrm dt = -F(x_2 - x_1)$. Similarly, if the right boundary is fixed then $\mathrm dx_n/\mathrm dt = 0$, otherwise $\eta\mathrm dx_n/\mathrm dt = F(x_n - x_{n-1})$. The parameter $\eta$ is called the damping constant or the viscosity coefficient and $F$ is the force law. In this model, cells are modelled as springs, whose forces are governed by this force law $F$.

The proliferation mechanism that we use assumes that at most one cell can divide at a time. Let $C_i(t)$ denote the event that the $i$th cell divides in the time interval $[t, t+\mathrm dt)$. We assume that $\Pr[C_i(t)] = G_i\,\mathrm dt$, where $G_i = G(|x_{i+1} - x_i|)$ for some proliferation law $G$. When this division occurs we place a new position at $(x_i + x_{i+1})/2$, adjusting the indices of the $x_i$ accordingly so that they remain sorted. We implement this mechanism by attempting a proliferation event every $\mathrm dt = \Delta t$ units of time, so that the probability of a proliferation event occuring at the time $t$ is $\sum_{i=1}^{n-1} G_i\Delta t$, and the probability that the $i$th cell proliferates, given that a proliferation event does occur, is $G_i/\sum_{j=1}^{n-1} G_j$. DiffEqCallbacks.jl is used to implement the callback that attempts this event periodically while simulating the problem.

Examples of how to simulate these systems are given in the sidebar.

EpithelialDynamics1D.CellProblem — Type
CellProblem{F,DF,DP,G,GP}

Struct for representing a cell simulation problem.

Fields

  • force_law::F: The force law, e.g. (δ, p) -> p.k * (p.s - δ).
  • force_law_parameters::Fp: The parameters for the force law.
  • proliferation_law::G = (δ, p) -> zero(δ): The proliferation law, e.g. (δ, p) -> p.β.
  • proliferation_law_parameters::Gp = nothing: The parameters for the proliferation law.
  • proliferation_period::Float64 = 1e-2: How often to attempt a proliferation event.
  • initial_time::Float64 = 0.0: The initial time.
  • final_time::Float64: The final time.
  • fix_left::Bool = true: Whether to fix the left endpoint.
  • fix_right::Bool = true: Whether to fix the right endpoint.
  • damping_constant::Float64: The damping constant.
  • initial_condition::Vector{Float64}: The initial cell positions. Must be sorted.

Solving

You can solve a CellProblem like you would solve a problem from DifferentialEquations.jl, e.g. with 

solve(prob, Tsit5(), saveat = 0.1)

The other kwargs for solve that we provide are:

  • proliferation::Bool = false: Whether to include proliferation.
  • rng::AbstractRNG = Random.default_rng(): The random number generator to use for proliferation.
  • sort::Bool = false: Whether to sort the cells after each step. This is useful for preventing the cells from becoming unsorted, which can happen due to numerical errors or force laws like F(q) = k/q^(n+1), n ≥ 0.

For proliferation problems, EnsembleProblem can be useful - see the docs.

source
EpithelialDynamics1D.SteadyCellProblem — Type
SteadyCellProblem{M}

Defines a steady-state problem for a CellProblem. Only has a single field, prob::CellProblem.

You can solve this problem as you would a NonlinearProblem, e.g. with 

solve(prob, alg)

where alg is e.g. DynamicSS(TRBDF2())) from SteadyStateDiffEq.jl.

source
FiniteVolumeMethod1D.FVMProblem — Type
FVMProblem(prob::CellProblem,
    mesh_points=copy(prob.initial_condition);
    diffusion_function=:continuum,
    diffusion_parameters=nothing,
    reaction_function=:continuum,
    reaction_parameters=nothing,
    proliferation=false)

Constructs an FVMProblem from a given CellProblem.

source
MovingBoundaryProblems1D.MBProblem — Type
MBProblem(prob::CellProblem,
    mesh_points=copy(prob.initial_condition);
    diffusion_function=:continuum,
    diffusion_parameters=nothing,
    reaction_function=:continuum,
    reaction_parameters=nothing,
    moving_boundary_function=:continuum,
    moving_boundary_parameters=nothing,
    rhs_function = :continuum,
    rhs_parameters=nothing,
    proliferation=false)

Constructs an MBProblem from a given CellProblem.

source
EpithelialDynamics1D.cell_densities — Method
cell_densities(cell_positions::AbstractVector{T}) where {T<:Number}

Compute the cell densities from the cell positions. The returned vector q has length(q) = length(cell_positions) - 1, with q[i] the density of the ith cell (cell_positions[i], cell_positions[i+1]) given by 1/(cell_positions[i+1] - cell_positions[i]).

If you want estimates of the densities at each node rather than for a cell, see node_densities.

source
EpithelialDynamics1D.cell_midpoints — Method
cell_midpoints(cell_positions::AbstractVector{T}) where {T<:Number}

Compute the midpoints of the cells from the cell positions. The returned vector x has length(x) = length(cell_positions) - 1, with x[i] the midpoint of the ith cell (cell_positions[i], cell_positions[i+1]) given by 0.5(cell_positions[i] + cell_positions[i+1]).

source
EpithelialDynamics1D.cell_numbers — Method
cell_numbers(sol::EnsembleSolution; indices = eachindex(sol), alpha=0.05)

Computes summary statistics for the cell numbers from an EnsembleSolution to a CellProblem.

Arguments

  • sol::EnsembleSolution: The ensemble solution to a CellProblem.

Keyword Arguments

  • indices = eachindex(sol): The indices of the cell simulations to consider.
  • alpha::Float64 = 0.05: The significance level for the confidence intervals.

Outputs

  • N::Vector{Vector{Int}}: The cell numbers for each cell simulation.
  • means::Vector{Float64}: The mean cell numbers for each cell simulation.
  • lowers::Vector{Float64}: The lower bounds of the confidence intervals for the cell numbers for each cell simulation.
  • uppers::Vector{Float64}: The upper bounds of the confidence intervals for the cell numbers for each cell simulation.
source
EpithelialDynamics1D.continuum_limit — Function
continuum_limit(prob::CellProblem,
    mesh_points=copy(prob.initial_condition);
    proliferation=false)

Returns the FVMProblem (if the right boundary is fixed) or the MBProblem (if the right boundary is free) that defines the continuum limit to the CellProblem prob. The argument mesh_points defines the mesh points to use for the PDE problem, which can be a grid of points or a number of points to use. If the problem should be considered to have proliferation, set proliferation = true.

source
EpithelialDynamics1D.get_knots — Function
get_knots(sol, num_knots = 500; indices = eachindex(sol), stat=maximum, parallel=false)

Computes knots for each time, covering the extremum of the cell positions across all cell simulations. You can restrict the simultaions to consider using the indices. The knots are obtained by applying stat, a tuple of functions for the left and right sides, to the vector of extrema at each time. For example, if stat=maximum then, at each time, the knots range between the smallest position observed and the maximum position observed across each simulation at that time.

source
EpithelialDynamics1D.leading_edges — Method
leading_edges(sol::EnsembleSolution; indices = eachindex(sol), alpha=0.05)

Computes summary statistics for the leading edges from an EnsembleSolution to a CellProblem.

Arguments

  • sol::EnsembleSolution: The ensemble solution to a CellProblem.

Keyword Arguments

  • indices = eachindex(sol): The indices of the cell simulations to consider.
  • alpha::Float64 = 0.05: The significance level for the confidence intervals.

Outputs

  • L::Vector{Vector{Float64}}: The leading edges for each cell simulation.
  • means::Vector{Float64}: The mean leading edges for each cell simulation.
  • lowers::Vector{Float64}: The lower bounds of the confidence intervals for the leading edges for each cell simulation.
  • uppers::Vector{Float64}: The upper bounds of the confidence intervals for the leading edges for each cell simulation.
source
EpithelialDynamics1D.node_densities — Method
node_densities(sol::EnsembleSolution;
    indices=eachindex(sol),
    num_knots=500,
    stat=(minimum, maximum),
    knots=get_knots(sol, num_knots; indices, stat),
    alpha=0.05,
    interp_fnc=(u, t) -> LinearInterpolation{true}(u, t),
    smooth_boundary=true)

Computes summary statistics for the node densities from an EnsembleSolution to a CellProblem. Any negative values are set to zero.

Arguments

  • sol::EnsembleSolution: The ensemble solution to a CellProblem.

Keyword Arguments

  • indices = eachindex(sol): The indices of the cell simulations to consider.
  • num_knots::Int = 500: The number of knots to use for the spline interpolation.
  • stat = (minimum, maximum): How to summarise the knots for get_knots.
  • parallel = false: Whether to use multithreading for the loops.
  • knots::Vector{Vector{Float64}} = get_knots(sol, num_knots; indices, stat, parallel): The knots to use for the spline interpolation.
  • alpha::Float64 = 0.05: The significance level for the confidence intervals.
  • interp_fnc = (u, t) -> LinearInterpolation{true}(u, t): The function to use for constructing the interpolant.
  • smooth_boundary::Bool = true: Whether to use the smooth boundary node densities.
  • extrapolate = false: Whether to allow for extrapolation when considering evaluating outside the range of cells for a given time and a given simulation.

Outputs

  • q::Vector{Vector{Vector{Float64}}}: The node densities for each cell simulation.
  • r::Vector{Vector{Vector{Float64}}}: The cell positions for each cell simulation.
  • means::Vector{Vector{Float64}}: The mean node densities for each cell simulation.
  • lowers::Vector{Vector{Float64}}: The lower bounds of the confidence intervals for the node densities for each cell simulation.
  • uppers::Vector{Vector{Float64}}: The upper bounds of the confidence intervals for the node densities for each cell simulation.
  • knots::Vector{Vector{Float64}}: The knots used for the spline interpolation.
source
EpithelialDynamics1D.node_densities — Method
node_densities(cell_positions::AbstractVector{T}; smooth_boundary=true) where {T<:Number}

Compute the cell densities from the cell positions, assigning a density to each node. If smooth_boundary is true, then

\[q_i(t) = \begin{cases} \dfrac{2}{x_2(t)-x_1(t)} - \dfrac{2}{x_3(t)-x_1(t)} & i=1, \\[6pt] \dfrac{2}{x_{i+1}(t)-x_{i-1}(t)} & i=2,\ldots,n-1, \\[6pt] \dfrac{2}{x_n(t) - x_{n-1}(t)} - \dfrac{2}{x_n(t)-x_{n-2}(t)} & i=n. \end{cases}\]

Otherwise,

\[q_i(t) = \begin{cases} \dfrac{1}{x_2(t)-x_1(t)} & i=1, \\[6pt] \dfrac{2}{x_{i+1}(t)-x_{i-1}(t)} & i=2,\ldots,n-1, \\[6pt] \dfrac{1}{x_n(t) - x_{n-1}(t)} & i=n. \end{cases}\]

If you want estimates of the densities for each cell rather than at each node, see cell_densities.

source