Semi-Vertex Model

Vertex models are a very popular choice in describing multicellular systems. They are actively being used in great variety such as to describe mechanical properties of plant cells [1] or organoid structures of epithelial cells [2,3].

Mathematical Description

In this model, we are only concerned with cellular forces and their representation in space. One single cell-agent can be described by a collection of (ordered) vertices which in turn also allows for a dual description in terms of edges.

$$\begin{align} \{\vec{v}_i\}_{i=0\dots n}\\ \vec{v}_i = \begin{bmatrix}v_{i,0}\\v_{i,1}\end{bmatrix} \end{align}$$

In the following text, we assume that vertices are always ordered (clockwise or anti-clockwise) and this ordering is identical for every cell in our simulation.

Mechanics

Every vertex is connected to its next neighbours in order via springs with an associated length $d$ and spring constant $\gamma$. The potential used to calculate the force $F_i$ acting along the edges of the cell between vertex $i$ and $i+1$ is given by

$$\begin{align} \vec{F}_{\text{edges},i} &= - \gamma \left(|\vec{v}_i - \vec{v}_{i+1}| - d\right) \frac{\vec{v}_i - \vec{v}_{i+1}}{|\vec{v}_i - \vec{v}_{i+1}|}\\ %V_\text{edges} &= \sum\limits_{i=0}^n \frac{\gamma}{2}\left(d_i - d\right)^2 \end{align}$$

where $d_i = |\vec{v}_i - \vec{v}_{i+1}|$ is the distance between individual vertices.

From the length of the individual edges, we can determine the total 2D volume $V$ of the cell when the equilibrium configuration of a perfect hexagon is reached.

$$\begin{equation} V = d^2\sum\limits_{i=0}^{n-1}\frac{1}{2\sin(\pi/n)} \end{equation}$$

However, since the individual vertices are mobile, we require an additional mechanism which simulates a central pressure $P$ depending on the currently measured volume $\tilde{V}$. This area can be calculated by summing over the individual areas of the triangles given by two adjacent vertices and the center point $\vec{c}=\sum_i\vec{v}_i/(n+1)$. They can be calculated by using the parallelogramm formula

$$\begin{align} \tilde{V}_i &= \det\begin{vmatrix} \vec{v}_{i+1} - \vec{c} & \vec{v}_i - \vec{c} \end{vmatrix}\\ &= \det\begin{pmatrix} (\vec{v}_{i+1} - \vec{c})_0 & (\vec{v}_{i} - \vec{c})_0\\ (\vec{v}_{i+1} - \vec{c})_1 & (\vec{v}_{i} - \vec{c})_1 \end{pmatrix}\\ \tilde{V} &= \sum\limits_{i=0}^{n-1}\tilde{V}_i \end{align}$$

The resulting force then points from the center of the cell $\vec{c}$ towards the individual vertices $\vec{v}_i$.

$$\begin{align} \vec{F}_{\text{pressure},i} = P\left(V-\tilde{V}\right)\frac{\vec{v}_i - \vec{c}}{|\vec{v}_i - \vec{c}|} \end{align}$$

These mechanical considerations alone are enough to yield perfect hexagonal configurations for individual cells without any interactions. If we also take into account an external force acting on the cell, the total force acting on the individual vertices $\vec{v}_i$ can be calculated via

$$\begin{equation} \vec{F}_{\text{total},i} = \vec{F}_{\text{external},i} + \vec{F}_{\text{edges},i} + \vec{F}_{\text{pressure},i} \end{equation}$$

Interaction

Cell-agents are interacting via forces $\vec{F}(\vec{p},\vec{q})$ which are dependent on two points $\vec{p}$ and $\vec{q}$ in either cell. The mechanical model we are currently using does not fully capture the essence of these cellular interactions. In principle, we would have to calculate the total force $\vec{F}’$ by integrating over all points either inside the cell or on its boundary but for the sake of simplicity we consider a different approach.

Let us denote the vertices of the two cells in question with $\{\vec{v}_i\}$ and $\{\vec{w}_j\}$.

Given a vertex $v_i$ of the first cell and the other cells vertices $\{w_j\}$ we asses if the vertex $v_i$ is contained in the polygon $\{w_j\}$ and calculate either the inside or a outside interaction force. In this way, by iterating over all vertices $v_i$, we can calculate the total external force $\vec{F}_\text{external}$ as the sum of all individual contributions. The same procedure when switching $v_i$ and $w_j$ results in a symmetric interaction.

Case 1: Outside Interaction

In this case, we assume that the vertex $\vec{v}_i$ in question is not inside the other cell. We make the simplified assumption that each vertex $\vec{v}_i$ is interacting with the closest point on the outer edge of the other cell. Given these sets of vertices, we calculate for each vertex $\vec{v}_i$ the closest point $$\begin{equation}\vec{p} = (1-q)\vec{w}_j + q\vec{w}_{j+1}\end{equation}$$ (assuming that we set $\vec{w}_{j+1}=\vec{w}_1$ when $j=N_\text{vertices}$) on the edge and then the force acting on this vertex can be calculated $$\begin{equation}\vec{F}_{\text{outside},i} = \vec{V}(\vec{v}_i, \vec{p})\end{equation}$$ by applying $\vec{F}$ on them. The force acting on the other cell acts on the vertices $j$ and $j+1$ with relative strength $1-q$ and $q$ respectively. $$\begin{alignat}{5} &\vec{F}_{\text{outside},j} &=& - &(1-q)&\vec{V}(\vec{v}_i,\vec{p})\\ &\vec{F}_{\text{outside},j+1} &=& &-q&\vec{V}(\vec{v}_i,\vec{p}) \end{alignat}$$

ℹ️

In the actual implementation of this approach, we additionally use a bounding box around each cell to quickly identify if a given vertex is outside the other cell.

Case 2: Inside Interaction

In the second case, a vertex of the other cell $\vec{w}_j$ has managed to move inside. Here, a different force $\vec{W}$ acts which is responsible for pushing the vertex outwards. The force is calculated between the center of our cell $$\begin{equation}\vec{v}_c = \frac{1}{N_\text{vertices}}\sum\limits_i \vec{v}_i\end{equation}$$ and the external vertex in question. The force which is calculated this way acts in equal parts on all vertices. $$\begin{alignat}{5} &\vec{F}_{\text{inside},j} &=& &&\frac{1}{N_\text{vertices}}\vec{W}(\vec{v}_c,\vec{w}_j)\\ &\vec{F}_{\text{inside},i} &=&-&&\frac{1}{N_\text{vertices}}\vec{W}(\vec{v}_c,\vec{w}_j) \end{alignat}$$

Cycle

We introduce an additional mechanism by which every cells area grows linearly over time. The growth parameter is chosen from a uniform distribution initially.

$$\begin{equation} \dot{V} = \alpha \end{equation}$$

By this process, cells will start to push on each other and thus expand the whole tissue structure.

Parameters

Parameter	Symbol	Value
Area	$V$	$500$
Spring Tension	$\gamma$	$2.0$
Central Pressure	$P$	$0.5$
Interaction Range	$\beta$	$0.5$
Potential Strength	$V_0$	$10.0$
Damping	$\lambda$	$0.2$
Growth Rate	$\alpha$	$5.0$
Domain Size	$L$	$800$
Simulation Steps	$N_\text{step}$	$20000$
Time Increment	$\Delta t$	$0.005$

Results

Initial State

Cells are placed in a perfect hexagonal grid such that edges and vertices align. The only distinction between the cells is their varying growth rate as described earlier.

Movie

The cells have grown and pushed on each other thus creating small spaces in between them.

Final State

Code

Cell Properties

cellular_raza-examples/semi_vertex/src/cell_properties.rs

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101


use cellular_raza::prelude::*;

use nalgebra::{Unit, Vector2};
use serde::{Deserialize, Serialize};

#[derive(Serialize, Deserialize, Clone, core::fmt::Debug)]
pub struct DirectedSphericalMechanics {
    pub pos: Vector2<f64>,
    pub vel: Vector2<f64>,
    pub orientation: Unit<Vector2<f64>>,
}

#[derive(Serialize, Deserialize, Clone, core::fmt::Debug)]
pub struct OutsideInteraction {
    pub potential_strength: f64,
    pub interaction_range: f64,
}

#[derive(Serialize, Deserialize, Clone, core::fmt::Debug)]
pub struct InsideInteraction {
    pub potential_strength: f64,
    pub average_radius: f64,
}

#[derive(Serialize, Deserialize, CellAgent, Clone, core::fmt::Debug)]
pub struct MyCell<const D: usize> {
    #[Mechanics]
    pub mechanics: VertexMechanics2D<D>,
    #[Interaction]
    pub interaction: VertexDerivedInteraction<OutsideInteraction, InsideInteraction>,
    pub growth_rate: f64,
}

impl Interaction<Vector2<f64>, Vector2<f64>, Vector2<f64>> for OutsideInteraction {
    fn calculate_force_between(
        &self,
        own_pos: &Vector2<f64>,
        _own_vel: &Vector2<f64>,
        ext_pos: &Vector2<f64>,
        _ext_vel: &Vector2<f64>,
        _ext_info: &(),
    ) -> Result<(Vector2<f64>, Vector2<f64>), CalcError> {
        // Calculate distance and direction between own and other point
        let z = ext_pos - own_pos;
        let r = z.norm();
        if r == 0.0 {
            return Ok(([0f64; 2].into(), [0f64; 2].into()));
        }
        let dir = z.normalize();

        // Introduce Non-dimensional length variable
        let sigma = r / (self.interaction_range);
        let spatial_cutoff = if r > self.interaction_range { 0.0 } else { 1.0 };

        // Calculate the strength of the interaction with correct bounds
        let strength = self.potential_strength * (1.0 - sigma);

        // Calculate only attracting and repelling forces
        let force = -dir * strength * spatial_cutoff;
        Ok((-force, force))
    }

    fn get_interaction_information(&self) -> () {}
}

impl Interaction<Vector2<f64>, Vector2<f64>, Vector2<f64>> for InsideInteraction {
    fn calculate_force_between(
        &self,
        own_pos: &Vector2<f64>,
        _own_vel: &Vector2<f64>,
        ext_pos: &Vector2<f64>,
        _ext_vel: &Vector2<f64>,
        _ext_info: &(),
    ) -> Result<(Vector2<f64>, Vector2<f64>), CalcError> {
        // Calculate direction between own and other point
        let z = own_pos - ext_pos;
        let r = z.norm();
        let dir = z.normalize();

        let force = self.potential_strength * dir / (0.5 + 0.5 * r / self.average_radius);
        Ok((-force, force))
    }

    fn get_interaction_information(&self) -> () {}
}

impl<const D: usize> Cycle for MyCell<D> {
    fn update_cycle(
        _rng: &mut rand_chacha::ChaCha8Rng,
        dt: &f64,
        cell: &mut Self,
    ) -> Option<CycleEvent> {
        let a = cell.mechanics.get_cell_area();
        cell.mechanics.set_cell_area(a + dt * cell.growth_rate);
        None
    }

    fn divide(_rng: &mut rand_chacha::ChaCha8Rng, _cell: &mut Self) -> Result<Self, DivisionError> {
        todo!()
    }
}

Custom Domain

cellular_raza-examples/semi_vertex/src/custom_domain.rs

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108


use cellular_raza::building_blocks::{CartesianCuboid, CartesianSubDomain};
use cellular_raza::concepts::*;

use serde::Serialize;

use crate::MyCell;

#[derive(Clone, Domain)]
pub struct MyDomain {
    #[DomainRngSeed]
    pub cuboid: CartesianCuboid<f64, 2>,
}

impl cellular_raza::concepts::DomainCreateSubDomains<MySubDomain> for MyDomain {
    type SubDomainIndex = usize;
    type VoxelIndex = [usize; 2];

    fn create_subdomains(
        &self,
        n_subdomains: std::num::NonZeroUsize,
    ) -> Result<
        impl IntoIterator<Item = (Self::SubDomainIndex, MySubDomain, Vec<Self::VoxelIndex>)>,
        DecomposeError,
    > {
        let subdomains = self.cuboid.create_subdomains(n_subdomains)?;
        Ok(subdomains
            .into_iter()
            .map(move |(subdomain_index, subdomain, voxels)| {
                (subdomain_index, MySubDomain { subdomain }, voxels)
            }))
    }
}

impl<const D: usize> cellular_raza::concepts::SortCells<MyCell<D>> for MyDomain {
    type VoxelIndex = [usize; 2];

    fn get_voxel_index_of(&self, cell: &MyCell<D>) -> Result<Self::VoxelIndex, BoundaryError> {
        let pos = cell.pos().row_mean().transpose();
        self.cuboid.get_voxel_index_of_raw(&pos)
    }
}

#[derive(Clone, SubDomain, Serialize)]
pub struct MySubDomain {
    #[Base]
    pub subdomain: CartesianSubDomain<f64, 2>,
}

impl<const D: usize> cellular_raza::concepts::SortCells<MyCell<D>> for MySubDomain {
    type VoxelIndex = [usize; 2];

    fn get_voxel_index_of(&self, cell: &MyCell<D>) -> Result<Self::VoxelIndex, BoundaryError> {
        let pos = cell.pos().row_mean().transpose();
        self.subdomain.get_index_of(pos)
    }
}

impl<const D: usize>
    cellular_raza::concepts::SubDomainMechanics<
        nalgebra::SMatrix<f64, D, 2>,
        nalgebra::SMatrix<f64, D, 2>,
    > for MySubDomain
{
    fn apply_boundary(
        &self,
        pos: &mut nalgebra::SMatrix<f64, D, 2>,
        vel: &mut nalgebra::SMatrix<f64, D, 2>,
    ) -> Result<(), BoundaryError> {
        // TODO refactor this with matrix multiplication!!!
        // This will probably be much more efficient and less error-prone!

        // For each position in the springs MyCell<D>
        pos.row_iter_mut()
            .zip(vel.row_iter_mut())
            .for_each(|(mut p, mut v)| {
                // For each dimension in the space
                for i in 0..p.ncols() {
                    // Check if the particle is below lower edge
                    if p[i] < self.subdomain.get_domain_min()[i] {
                        p[i] = 2.0 * self.subdomain.get_domain_min()[i] - p[i];
                        v[i] = v[i].abs();
                    }

                    // Check if the particle is over the edge
                    if p[i] > self.subdomain.get_domain_max()[i] {
                        p[i] = 2.0 * self.subdomain.get_domain_max()[i] - p[i];
                        v[i] = -v[i].abs();
                    }
                }
            });

        // If new pos is still out of boundary return error
        for j in 0..pos.nrows() {
            let p = pos.row(j);
            for i in 0..pos.ncols() {
                if p[i] < self.subdomain.get_domain_min()[i]
                    || p[i] > self.subdomain.get_domain_max()[i]
                {
                    return Err(BoundaryError(format!(
                        "Particle is out of domain at pos {:?}",
                        pos
                    )));
                }
            }
        }
        Ok(())
    }
}

Main Simulation

cellular_raza-examples/semi_vertex/src/main.rs

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101


use backend::chili;
use cellular_raza::prelude::*;

use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha8Rng;
use serde::{Deserialize, Serialize};

// Mechanical parameters
pub const CELL_MECHANICS_AREA: f64 = 500.0;
pub const CELL_MECHANICS_SPRING_TENSION: f64 = 2.0;
pub const CELL_MECHANICS_CENTRAL_PRESSURE: f64 = 0.5;
pub const CELL_MECHANICS_INTERACTION_RANGE: f64 = 5.0;
pub const CELL_MECHANICS_POTENTIAL_STRENGTH: f64 = 10.0;
pub const CELL_MECHANICS_DAMPING_CONSTANT: f64 = 0.2;
pub const CELL_MECHANICS_DIFFUSION_CONSTANT: f64 = 0.0;

// Parameters for domain
pub const DOMAIN_SIZE_X: f64 = 800.0;
pub const DOMAIN_SIZE_Y: f64 = 800.0;

// Time parameters
pub const N_TIMES: u64 = 20_001;
pub const DT: f64 = 0.005;
pub const T_START: f64 = 0.0;
pub const SAVE_INTERVAL: u64 = 50;

// Meta Parameters to control solving
pub const N_THREADS: usize = 10;

mod cell_properties;
mod custom_domain;

use cell_properties::*;
use custom_domain::*;
use time::FixedStepsize;

fn main() -> Result<(), chili::SimulationError> {
    // Fix random seed
    let mut rng = ChaCha8Rng::seed_from_u64(2);

    // Define the simulation domain
    let domain = MyDomain {
        cuboid: CartesianCuboid::from_boundaries_and_interaction_range(
            [0.0; 2],
            [DOMAIN_SIZE_X, DOMAIN_SIZE_Y],
            2.0 * VertexMechanics2D::<6>::inner_radius_from_cell_area(CELL_MECHANICS_AREA),
        )?,
    };

    // Define cell agents
    let models = VertexMechanics2D::fill_rectangle_flat_top(
        CELL_MECHANICS_AREA,
        CELL_MECHANICS_SPRING_TENSION,
        CELL_MECHANICS_CENTRAL_PRESSURE,
        CELL_MECHANICS_DAMPING_CONSTANT,
        CELL_MECHANICS_DIFFUSION_CONSTANT,
        [
            [0.1 * DOMAIN_SIZE_X, 0.1 * DOMAIN_SIZE_Y].into(),
            [0.9 * DOMAIN_SIZE_X, 0.9 * DOMAIN_SIZE_Y].into(),
        ],
    );
    println!("Generated {} cells", models.len());

    let growth_rate = 5.0;
    let cells = models
        .into_iter()
        .map(|model| MyCell {
            mechanics: model,
            interaction: VertexDerivedInteraction::from_two_forces(
                OutsideInteraction {
                    potential_strength: CELL_MECHANICS_POTENTIAL_STRENGTH,
                    interaction_range: CELL_MECHANICS_INTERACTION_RANGE,
                },
                InsideInteraction {
                    potential_strength: 1.5 * CELL_MECHANICS_POTENTIAL_STRENGTH,
                    average_radius: CELL_MECHANICS_AREA.sqrt(),
                },
            ),
            growth_rate: rng.gen_range(0.0..1.0) * growth_rate,
        })
        .collect::<Vec<_>>();

    // Define settings for storage and time solving
    let settings = chili::Settings {
        time: FixedStepsize::from_partial_save_steps(0.0, DT, N_TIMES, SAVE_INTERVAL)?,
        n_threads: N_THREADS.try_into().unwrap(),
        show_progressbar: true,
        storage: StorageBuilder::new()
            .location("out/semi_vertex")
            .priority([StorageOption::SerdeJson]),
    };

    // Run the simulation
    let _storager = chili::run_simulation!(
        agents: cells,
        domain: domain,
        settings: settings,
        aspects: [Mechanics, Interaction, Cycle],
    )?;
    Ok(())
}

Plotting

cellular_raza-examples/semi_vertex/plot.py

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185


import matplotlib
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import numpy as np
import json
import pandas as pd
from pathlib import Path
from glob import glob
import tqdm
from threading import Thread
import os

def get_last_output_path() -> Path:
    return Path(sorted(list(glob("out/semi_vertex/*")))[-1])

class Loader:
    def __init__(self, opath: Path = get_last_output_path()):
        self.opath = opath
        folders = glob(str(self.opath / "cells/json/*"))
        self.__iterations = np.sort(np.array([int(Path(f).name) for f in folders]))
        self.__all_files = {
            iteration: list(glob(str(self.opath / "cells/json/{:020}/*.json".format(iteration))))
            for iteration in self.__iterations
        }

    def load_cells(self, iteration: int) -> pd.DataFrame:
        # Load all json files
        results = []
        for file in self.__all_files[iteration]:
            f = open(file)
            batch = json.load(f)
            results.extend([b["element"][0] for b in batch["data"]])
        df = pd.json_normalize(results)
        df["cell.mechanics.points"] = df["cell.mechanics.points"].apply(
            lambda x: np.array(x, dtype=float).reshape((2, -1)).T
        )
        self.cells = df
        return df

    @property
    def iterations(self) -> np.ndarray:
        return np.copy(self.__iterations)

def run_plotter_iterations(instance: int, opath: Path, iterations, **kwargs):
    plotter = Plotter(opath)
    iterations = tqdm.tqdm(iterations, total=len(iterations), position=instance)
    plotter.plot_cells_at_iterations(iterations, **kwargs, progress_bar=False)

class Plotter:
    loader: Loader

    def __init__(self, opath: Path = get_last_output_path()):
        self.fig, self.ax = plt.subplots(figsize=(12, 12))
        self.loader = Loader(opath)
        self.cm = matplotlib.colormaps.get_cmap("YlGn")

    def plot_cells(self, df_cells):
        polys = [Polygon(pos) for pos in df_cells["cell.mechanics.points"]]
        growth_rates = df_cells["cell.growth_rate"]
        rate_max = np.max(growth_rates)
        t = [(0.5 + 0.5 * rate / rate_max) for rate in growth_rates]
        facecolors = [self.cm(ti) for ti in t]
        pc = PatchCollection(
            polys,
            facecolors=facecolors,
            edgecolors='white'
        )
        self.pc = self.ax.add_collection(pc)

    def plot_cells_at_iter(
            self,
            iteration: int,
            save_path: Path | None = None,
            overwrite: bool = False,
            transparent: bool = False,
        ):
        if save_path is None:
            save_path = self.loader.opath / "images"
        save_path.mkdir(parents=True, exist_ok=True)
        save_path = save_path / "snapshot-{:020}.png".format(iteration)
        if save_path.exists() and overwrite is False:
            return save_path
        cells = self.loader.load_cells(iteration)

        self.ax.set_xlim(0, 800)
        self.ax.set_ylim(0, 800)
        self.ax.set_axis_off()
        self.plot_cells(cells)
        self.fig.tight_layout()
        self.fig.savefig(save_path, transparent=transparent, dpi=100)
        self.pc.remove()
        return save_path

    def plot_cells_at_iterations(
            self,
            iterations: list[int] | np.ndarray,
            save_path: Path | None = None,
            overwrite: bool = False,
            transparent: bool = True,
            progress_bar: bool = True,
        ):
        iterations = tqdm.tqdm(iterations) if progress_bar else iterations
        for iteration in iterations:
            self.plot_cells_at_iter(iteration, save_path, overwrite, transparent)

def plot_cells_at_all_iterations(
    opath: Path = get_last_output_path(),
    save_path: Path | None = None,
    overwrite: bool = False,
    transparent: bool = False,
    n_threads: int | None = None,
):
    plt.close('all')
    matplotlib.use('svg')
    if n_threads is None:
        n_threads = os.cpu_count()
    results = []
    for i in range(0, n_threads):
        loader = Loader(opath)
        chunksize = int(len(loader.iterations) / n_threads)
        left_over = len(loader.iterations) % n_threads
        lower = i * chunksize
        upper = lower + chunksize
        if left_over > 0:
            upper += 1
            left_over -= 1
        iterations = loader.iterations[lower:upper]
        t = Thread(
            target=run_plotter_iterations,
            args=[i, loader.opath, iterations],
            kwargs={
                "save_path":save_path,
                "overwrite":overwrite,
                "transparent":transparent
            },
        )
        results.append(t)

    for t in results:
        t.start()

def plot_growth(opath: Path = get_last_output_path(), iteration: int | None = None):
    loader = Loader(opath)
    iteration = int(loader.iterations[0]) if iteration is None else iteration
    fig, ax = plt.subplots(1, 2, figsize=(12, 6))
    cells = loader.load_cells(iteration)
    values = np.array(cells["cell.growth_rate"])
    _, bins, patches = ax[0].hist(
        values,
        bins=int(len(values) / 10),
        linewidth=1,
        facecolor=(0.75, 0.75, 0.75),
        edgecolor='k',
        fill=True,
    )
    bin_centers = 0.5 * (bins[:-1] + bins[1:])

    col = bin_centers - min(bin_centers)
    col /= max(col)

    cm = matplotlib.colormaps.get_cmap("YlGn")
    for c, p in zip(col, patches):
        plt.setp(p, "facecolor", cm(0.5 + 0.5*c))
    ax[0].set_xlabel("Growth rate [area/simulation step]")
    ax[0].set_title("Distribution of growth-rates")
    sub_iterations = loader.iterations[::50]
    areas = np.array([
        loader.load_cells(iteration)["cell.mechanics.cell_area"]
        for iteration in sub_iterations
    ])
    ax[1].errorbar(
        sub_iterations,
        np.mean(areas, axis=1),
        np.std(areas, axis=1),
        color="k",
    )
    ax[1].set_xlabel("Simulation Step")
    ax[1].set_title("Average cell area")
    fig.tight_layout()
    fig.savefig("{}/growth.png".format(loader.opath))

if __name__ == "__main__":
    plot_cells_at_all_iterations(transparent=True, n_threads=1)
    plot_growth()

The code for this simulation and the visualization can be found in the examples folder of cellular_raza.

References

[1] R. M. H. Merks, M. Guravage, D. Inzé, and G. T. S. Beemster, “VirtualLeaf: An Open-Source Framework for Cell-Based Modeling of Plant Tissue Growth and Development” Plant Physiology, vol. 155, no. 2. Oxford University Press (OUP), pp. 656–666, Feb. 01, 2011. doi: 10.1104/pp.110.167619.

[2] A. G. Fletcher, M. Osterfield, R. E. Baker, and S. Y. Shvartsman, “Vertex Models of Epithelial Morphogenesis,” Biophysical Journal, vol. 106, no. 11. Elsevier BV, pp. 2291–2304, Jun. 2014. doi: 10.1016/j.bpj.2013.11.4498.

[3] D. L. Barton, S. Henkes, C. J. Weijer, and R. Sknepnek, “Active Vertex Model for cell-resolution description of epithelial tissue mechanics,” PLOS Computational Biology, vol. 13, no. 6. Public Library of Science (PLoS), p. e1005569, Jun. 30, 2017. doi: 10.1371/journal.pcbi.1005569.