Bacterial Branching Patterns

Spatio-temporal patterns of bacterial growth such as in Bacillus Subtilis have been studied for numerous years [1,2]. They are typically modeled by a system of PDEs (Partial Differential Equations) which describe intracellular reactions, growth and spatial processes (usually via Diffusion), such as diffusion of nutrients and movement of the cells. Here, we only consider two variables: the spatial-temporal distribution $n(x,t)$ of the available nutrients and the bacterial population density $b(x,t)$. The rescaled set of coupled partial differential equations read:

$$\begin{alignat}{5} \dot{n} &= &&\nabla^2n &- &f(n, b)\\ \dot{b} &= d&&\nabla^2b &+ \theta &f(n, b) \end{alignat}$$

The function $f(n,b)$ describes the nutrient consumption by the bacterial metabolism and $d$ is the ratio of diffusion constants. The parameter $\theta$ is the “gain” in bacterial mass per nutrient volume resulting from growth and division.

One critique of these models is that the pattern will diffuse over the course of time, thus not creating a persistent pattern. This comes from modeling cellular motility via a diffusion equation. However, stable patterns could be achieved if an equilibrium state exists where cells remain at their locations.

Mathematical Description

To formulate the above equations in an agent-based approach, we need to define cellular behaviour on an individual-based level.

Mechanics & Interaction

We represent cells as soft spheres with dynamics determined by NewtonDamped2DF32 and their interaction given by the MorsePotentialF32 type.

Intra- & extracellular Reactions

The nutrient resource is freely diffusible throughout the simulation domain. Individual cell-agents take up the extracellular nutrient resource and grow proportionally. Only a fraction of the nutrients is converted to actual growth of cell volume.

$$\begin{align} \dot{V}_c &= \alpha u n_e\\ \dot{n}_e(x) &= D\Delta n_e - u \sum\limits_{c=1}^N n_e(x)\delta(x-x_c) \end{align}$$

The components of these PDEs describe the extracellular nutrients as well as the change in volume of the individual cells. $n_e$ is the spatially distributed extra-cellular nutrient concentration which undergoes diffusion with the diffusion constant $D$ while $V_c$ is the volume of cell $c$ positioned at $x_c$. The parameter $u$ is the uptake rate of the nutrient while $\alpha$ describes the conversion of the nutrient by the cellular metabolism resulting in an increase of the volume $V_c$.

Cycle

Once cells have reached a threshold $\tau$ in size (measured in multiple of given initial radius), they will divide. The newly created agents inherit all parameter values and thus the individual behaviour of their mother cell. They will continue to take up nutrients, process them and divide. Usually, this cycle of producing new generations ceases due to depletion of nutrients after a few division events. For simplicity we ignore cell death in our simulation.

Parameters

The parameter values have been chosen such that our simulation yields realistic results.

Parameter	Symbol	Value
Cell Radius	$R$	$6.0 \text{ µm}$
Potential Strength	$V_0$	$2\text{ µm}^2\text{ }/\text{ min}^2$
Potential Stiffness	$\lambda$	$0.15\text{ µm}^{-1}$
Damping Constant	$\lambda$	$1\text{ min}^{-1}$
Interaction Range	$\xi$	$1.0 R$
Uptake Rate	$u$	$1.0 \text{ min}^{-1}$
Growth Rate	$\alpha$	$13.0 \text{ µm}^3\text{ l }/ \text{ µg}$
Division threshold	$\tau$	$2.0R$
Diffusion Constant	$D$	$80 \text{ µm}^2\text{ }/\text{ min}$

Initial State

Property	Symbol	Value
Time Stepsize	$\Delta t$	$0.12\text{ min}$
Time Steps	$N_t$	$20'000$
Domain Size	$L$	$3000\text{ µm}$
Centered Starting Domain Size	$L_0$	$300 \text{ µm}$
Number of cells	$N_0$	$400$
Initial Nutrients	$n_e$	$10 \text{ µg }/\text{ l}$

Results

Figure 1: Final snapshot of the fully grown bacterial colony.

ℹ️

The picture shown above was generated with modified parameters on a larger domain size of $L=10000\mu m$ and with an increased number of simulation steps $t_\text{max}=6000$. It is recommended to use multiple threads --threads XY to calculate this result.

Movie

Code

The simulation can be executed from a CLI tool which allows users to specify their own parameters.

CLI Arguments

# cargo run --bin cr_bacteria_branching -- --help
Usage: cr_bacteria_branching [OPTIONS]

Options:
  -h, --help     Print help
  -V, --version  Print version

Bacteria:
  -n, --n-bacteria-initial <N_BACTERIA_INITIAL>
          [default: 5]
  -r, --radius <RADIUS>
          [default: 6]
      --division-threshold <DIVISION_THRESHOLD>
          Multiple of the radius at which the cell will divide [default: 2]
      --potential-stiffness <POTENTIAL_STIFFNESS>
          [default: 0.15]
      --potential-strength <POTENTIAL_STRENGTH>
          [default: 2]
      --damping-constant <DAMPING_CONSTANT>
          [default: 1]
  -u, --uptake-rate <UPTAKE_RATE>
          [default: 1]
  -g, --growth-rate <GROWTH_RATE>
          [default: 13]

Domain:
  -d, --domain-size <DOMAIN_SIZE>
          Overall size of the domain [default: 3000]
      --voxel-size <VOXEL_SIZE>
          Size of one voxel containing individual cells.
          This value should be chosen `>=3*RADIUS`. [default: 30]
      --domain-starting-size <DOMAIN_STARTING_SIZE>
          Size of the square for initlal placement of bacteria [default: 100]
      --reactions-dx <REACTIONS_DX>
          Discretization of the diffusion process [default: 20]
      --diffusion-constant <DIFFUSION_CONSTANT>
          [default: 80]
      --initial-concentration <INITIAL_CONCENTRATION>
          [default: 10]

Time:
      --dt <DT>                        [default: 0.1]
      --tmax <TMAX>                    [default: 2000]
      --save-interval <SAVE_INTERVAL>  [default: 200]

Other:
      --threads <THREADS>  Meta Parameters to control solving [default: 2]

Simulation

Cargo

cellular_raza-examples/bacterial_branching/Cargo.toml

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


[package]
name = "cr_bacteria_branching"
version = "0.1.0"
edition = "2021"
authors = ["Jonas Pleyer <jonas.pleyer@fdm.uni-freiburg.de>"]

[dependencies]
serde = { workspace = true, features=["rc"] }
rand = { workspace = true, features=["small_rng"] }
rand_chacha = { workspace = true }
nalgebra = { workspace = true }
num = { workspace = true }
cellular_raza = { path="../../cellular_raza", features=["default"] }
plotters = { workspace = true }
rayon = "1.10.0"
ndarray = { workspace = true, features = ["blas", "serde", "serde-1"] }
clap = { version = "4.5.31", features = ["derive"] }

ℹ️

The dependencies which are derived from the workspace either via workspace = true or path="../../" should be replaced with the versions used in the workspace Cargo.toml.

Bacterial Properties

cellular_raza-examples/bacterial_branching/src/bacterial_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


use core::f32::consts::{PI, SQRT_2};

use cellular_raza::prelude::*;
use serde::{Deserialize, Serialize};

use crate::ReactionVector;

#[derive(Clone, Serialize, Deserialize, CellAgent)]
pub struct MyAgent {
    #[Mechanics]
    pub mechanics: NewtonDamped2DF32,
    #[Interaction]
    pub interaction: MorsePotentialF32,
    pub uptake_rate: f32,
    pub division_radius: f32,
    pub growth_rate: f32,
}

impl Cycle<MyAgent, f32> for MyAgent {
    fn update_cycle(
        _rng: &mut rand_chacha::ChaCha8Rng,
        _dt: &f32,
        cell: &mut MyAgent,
    ) -> Option<CycleEvent> {
        // If the cell is not at the maximum size let it grow
        if cell.interaction.radius > cell.division_radius {
            return Some(CycleEvent::Division);
        }
        None
    }

    fn divide(
        rng: &mut rand_chacha::ChaCha8Rng,
        c1: &mut MyAgent,
    ) -> Result<MyAgent, DivisionError> {
        // Clone existing cell
        let mut c2 = c1.clone();

        let r = c1.interaction.radius;

        // Make both cells smaller
        // Also keep old cell larger
        c1.interaction.radius /= SQRT_2;
        c2.interaction.radius /= SQRT_2;

        // Generate cellular splitting direction randomly
        use rand::Rng;
        let alpha = rng.random_range(0.0..2.0 * PI);
        let dir_vec = nalgebra::Vector2::from([alpha.cos(), alpha.sin()]);

        // Define new positions for cells
        // It is randomly chosen if the old cell is left or right
        let offset = dir_vec * r / SQRT_2;
        let old_pos = c1.pos();

        c1.set_pos(&(old_pos + offset));
        c2.set_pos(&(old_pos - offset));

        Ok(c2)
    }
}

// COMPONENT DESCRIPTION
// 0         CELL AREA
impl Intracellular<ReactionVector> for MyAgent {
    fn set_intracellular(&mut self, intracellular: ReactionVector) {
        self.interaction.radius = (intracellular[0] / PI).powf(0.5);
    }

    fn get_intracellular(&self) -> ReactionVector {
        vec![PI * self.interaction.radius.powf(2.0)].into()
    }
}

impl ReactionsExtra<ReactionVector, ReactionVector> for MyAgent {
    fn calculate_combined_increment(
        &self,
        _intracellular: &ReactionVector,
        extracellular: &ReactionVector,
    ) -> Result<(ReactionVector, ReactionVector), CalcError> {
        let extra = extracellular;
        let u = self.uptake_rate;

        let uptake = u * extra;

        let incr_intra: ReactionVector = vec![self.growth_rate * uptake[0]].into();
        let incr_extra = -uptake;

        Ok((incr_intra, incr_extra))
    }
}

Domain

cellular_raza-examples/bacterial_branching/src/subdomain.rs

Main Simulation

cellular_raza-examples/bacterial_branching/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
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
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218


use cellular_raza::prelude::*;

use clap::{Args, Parser};
use nalgebra::Vector2;
use num::Zero;

use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha8Rng;

use serde::{Deserialize, Serialize};

mod bacteria_properties;

use bacteria_properties::*;

pub(crate) type ReactionVector = nalgebra::DVector<f32>;

#[derive(Clone, Args, Debug)]
#[group()]
#[clap(next_help_heading = Some("Bacteria"))]
struct BacterialParameters {
    #[arg(short, long, default_value_t = 5)]
    n_bacteria_initial: u32,
    #[arg(short, long, default_value_t = 6.0)]
    radius: f32,
    /// Multiple of the radius at which the cell will divide
    #[arg(long, default_value_t = 2.0)]
    division_threshold: f32,
    #[arg(long, default_value_t = 0.15)]
    potential_stiffness: f32,
    #[arg(long, default_value_t = 2.0)]
    potential_strength: f32,
    #[arg(long, default_value_t = 1.0)]
    damping_constant: f32,
    #[arg(short, long, default_value_t = 1.0)]
    uptake_rate: f32,
    #[arg(short, long, default_value_t = 13.0)]
    growth_rate: f32,
}

#[derive(Clone, Args, Debug)]
#[group()]
#[clap(next_help_heading = Some("Domain"))]
struct DomainParameters {
    /// Overall size of the domain
    #[arg(short, long, default_value_t = 3000.0)]
    domain_size: f32,
    #[arg(
        long,
        default_value_t = 30.0,
        help = "\
        Size of one voxel containing individual cells.\n\
        This value should be chosen `>=3*RADIUS`.\
    "
    )]
    voxel_size: f32,
    /// Size of the square for initlal placement of bacteria
    #[arg(long, default_value_t = 100.0)]
    domain_starting_size: f32,
    /// Discretization of the diffusion process
    #[arg(long, default_value_t = 20.0)]
    reactions_dx: f32,
    #[arg(long, default_value_t = 80.0)]
    diffusion_constant: f32,
    #[arg(long, default_value_t = 10.0)]
    initial_concentration: f32,
}

#[derive(Clone, Args, Debug)]
#[group()]
#[clap(next_help_heading = Some("Time"))]
struct TimeParameters {
    #[arg(long, default_value_t = 0.1)]
    dt: f32,
    #[arg(long, default_value_t = 2000.0)]
    tmax: f32,
    #[arg(long, default_value_t = 200)]
    save_interval: usize,
}

#[derive(Clone, Parser, Debug)]
#[command(version, about, long_about = None)]
struct Parameters {
    #[command(flatten)]
    bacteria: BacterialParameters,

    #[command(flatten)]
    domain: DomainParameters,

    #[command(flatten)]
    time: TimeParameters,

    #[clap(help_heading = Some("Other"))]
    /// Meta Parameters to control solving
    #[arg(long, default_value_t = 2)]
    threads: usize,
}

fn main() -> Result<(), SimulationError> {
    let parameters = Parameters::parse();
    run_sim(parameters)
}

fn run_sim(parameters: Parameters) -> Result<(), SimulationError> {
    let Parameters {
        bacteria:
            BacterialParameters {
                n_bacteria_initial,
                radius: cell_radius,
                division_threshold,
                potential_stiffness,
                potential_strength,
                damping_constant,
                uptake_rate,
                growth_rate,
            },
        domain:
            DomainParameters {
                domain_size,
                voxel_size: domain_voxel_size,
                domain_starting_size,
                reactions_dx,
                diffusion_constant,
                initial_concentration,
            },
        time:
            TimeParameters {
                dt,
                tmax: t_max,
                save_interval,
            },
        threads: n_threads,
    } = parameters;

    let ds = domain_size / 2.0;
    let dx = domain_starting_size / 2.0;

    // Fix random seed
    let mut rng = ChaCha8Rng::seed_from_u64(2);

    let cells = (0..n_bacteria_initial)
        .map(|_| {
            let x = rng.random_range(ds - dx..ds + dx);
            let y = rng.random_range(ds - dx..ds + dx);

            let pos = Vector2::from([x, y]);
            MyAgent {
                mechanics: NewtonDamped2DF32 {
                    pos,
                    vel: Vector2::zero(),
                    damping_constant,
                    mass: 1.0,
                },
                interaction: MorsePotentialF32 {
                    radius: cell_radius,
                    potential_stiffness,
                    cutoff: 2.0 * division_threshold * cell_radius,
                    strength: potential_strength,
                },
                uptake_rate,
                division_radius: division_threshold * cell_radius,
                growth_rate,
            }
        })
        .collect::<Vec<_>>();

    let cond = dt - 0.5 * reactions_dx / diffusion_constant;
    if cond >= 0.0 {
        println!(
            "❗❗❗WARNING❗❗❗\n\
            The stability condition \
            dt <= 0.5 dx^2/D for the integration \
            method is not satisfied. This can \
            lead to solving errors and inaccurate \
            results."
        );
    }

    if domain_voxel_size < division_threshold * cell_radius {
        println!(
            "❗❗❗WARNING❗❗❗\n\
            The domain_voxel_size {domain_voxel_size} has been chosen \
            smaller than the length of the interaction {}. This \
            will probably yield incorrect results.",
            division_threshold * cell_radius,
        );
    }

    let domain = CartesianDiffusion2D {
        domain: CartesianCuboid::from_boundaries_and_interaction_range(
            [0.0; 2],
            [domain_size, domain_size],
            domain_voxel_size,
        )?,
        reactions_dx: [reactions_dx; 2].into(),
        diffusion_constant,
        initial_value: ReactionVector::from(vec![initial_concentration]),
    };

    let storage = StorageBuilder::new().priority([StorageOption::SerdeJson]);
    let time = FixedStepsize::from_partial_save_freq(0.0, dt, t_max, save_interval)?;
    let settings = Settings {
        n_threads: n_threads.try_into().unwrap(),
        time,
        storage,
        show_progressbar: true,
    };

    let _storager = run_simulation!(
        agents: cells,
        domain: domain,
        settings: settings,
        aspects: [Mechanics, Interaction, ReactionsExtra, Cycle],
        parallelizer: Rayon,
        zero_reactions_default: |_| nalgebra::DVector::zeros(1),
    )?;
    Ok(())
}

Plotting

Plotting

cellular_raza-examples/bacterial_branching/src/plotting.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
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
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236


import json
from pathlib import Path
from glob import glob
import os
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import tqdm
import multiprocessing as mp
import itertools


def get_last_output_path(search_dir: Path | None = Path("out")) -> Path:
    """
    Parameters
    ----------
    search_dir: Path
        Directory to search for results. Defaults to "out".

    Returns
    -------
    output_path: Path
        Most recent folder in the given search_dir.
    """
    return Path(sorted(list(glob(str(search_dir) + "/*")))[-1])


def _get_all_iteration_files(output_path: Path = get_last_output_path()) -> list[Path]:
    return [Path(p) for p in sorted(glob(str(output_path) + "/cells/json/*"))]


def _iteration_to_file(iteration: int, output_path: Path, cs: str = "cells") -> Path:
    return output_path / "{}/json/{:020}".format(cs, iteration)


def get_all_iterations(output_path: Path = get_last_output_path()) -> list[int]:
    iterations_files = _get_all_iteration_files(output_path)
    return [int(os.path.basename(it)) for it in iterations_files]


def load_cells_at_iteration(
    iteration: int,
    output_path: Path,
):
    iteration_file = _iteration_to_file(iteration, output_path, "cells")
    data = []
    for filename in glob(str(iteration_file) + "/*"):
        file = open(filename)
        di = json.load(file)["data"]
        data.extend([b["element"][0] for b in di])
    df = pd.json_normalize(data)
    for key in [
        "cell.mechanics.pos",
        "cell.mechanics.vel",
    ]:
        df[key] = df[key].apply(lambda x: np.array(x, dtype=float))
    return df


def load_subdomains_at_iteration(
    iteration: int, output_path: Path = get_last_output_path()
) -> pd.DataFrame:
    iteration_file = _iteration_to_file(iteration, output_path, "subdomains")
    data = []
    for filename in glob(str(iteration_file) + "/*"):
        file = open(filename)
        di = json.load(file)["element"]
        data.append(di)
    df = pd.json_normalize(data)
    for key in [
        "subdomain.domain_min",
        "subdomain.domain_max",
        "subdomain.min",
        "subdomain.max",
        "subdomain.dx",
        "subdomain.voxels",
        "reactions_min",
        # "reactions_max",
        "reactions_dx",
        "extracellular.data",
        "ownership_array.data",
    ]:
        df[key] = df[key].apply(lambda x: np.array(x, dtype=float))
    return df


def plot_iteration(
    iteration: int,
    intra_bounds: tuple[float, float],
    extra_bounds: tuple[float, float],
    output_path: Path = get_last_output_path(),
    save_figure: bool = True,
    figsize: int = 32,
) -> matplotlib.figure.Figure | None:
    dfc = load_cells_at_iteration(iteration, output_path)
    dfs = load_subdomains_at_iteration(iteration, output_path)

    # Set size of the image
    domain_min = dfs["subdomain.domain_min"][0]
    domain_max = dfs["subdomain.domain_max"][0]
    fig, ax = plt.subplots(figsize=(figsize, figsize))
    ax.set_xlim([domain_min[0], domain_max[0]])
    ax.set_ylim([domain_min[1], domain_max[1]])

    # Plot background
    max_size = np.max([dfsi["index_max"] for _, dfsi in dfs.iterrows()], axis=0)
    all_values = np.zeros(max_size)
    for n_sub, dfsi in dfs.iterrows():
        values = dfsi["extracellular.data"].reshape(dfsi["extracellular.dim"])[:, :, 0]
        filt = dfsi["ownership_array.data"].reshape(dfsi["ownership_array.dim"])
        filt = filt[1:-1, 1:-1]

        index_min = np.array(dfsi["index_min"])
        slow = index_min
        shigh = index_min + np.array(values.shape)
        all_values[slow[0] : shigh[0], slow[1] : shigh[1]] += values * filt
    ax.imshow(
        all_values.T,
        vmin=extra_bounds[0],
        vmax=extra_bounds[1],
        extent=(domain_min[0], domain_max[0], domain_min[1], domain_max[1]),
        origin="lower",
    )

    # Plot cells
    points = np.array([p for p in dfc["cell.mechanics.pos"]])
    radii = np.array([r for r in dfc["cell.interaction.radius"]])
    radii_div = np.array([r for r in dfc["cell.division_radius"]])
    s = np.clip(
        (radii / radii_div - intra_bounds[0]) / (intra_bounds[1] - intra_bounds[0]),
        0,
        1,
    )

    color_high = np.array([233, 170, 242]) / 255
    color_low = np.array([129, 12, 145]) / 255
    color = np.tensordot((1 - s), color_low, 0) + np.tensordot(s, color_high, 0)

    # Plot cells as circles
    from matplotlib.patches import Circle
    from matplotlib.collections import PatchCollection

    collection = PatchCollection(
        [
            Circle(
                points[i, :],
                radius=radii[i],
            )
            for i in range(points.shape[0])
        ],
        facecolors=color,
    )
    ax.add_collection(collection)
    ax.text(
        0.05,
        0.05,
        "Agents: {:9}".format(points.shape[0]),
        transform=ax.transAxes,
        fontsize=14,
        verticalalignment="center",
        bbox=dict(boxstyle="square", facecolor="#FFFFFF"),
    )

    ax.set_axis_off()
    if save_figure:
        os.makedirs(output_path / "images", exist_ok=True)
        fig.savefig(
            output_path / "images/cells_at_iter_{:010}".format(iteration),
            bbox_inches="tight",
            pad_inches=0,
        )
        plt.close(fig)
        return None
    else:
        return fig


def __plot_all_iterations_helper(args_kwargs):
    iteration, kwargs = args_kwargs
    plot_iteration(iteration, **kwargs)


def plot_all_iterations(
    intra_bounds: tuple[float, float],
    extra_bounds: tuple[float, float],
    output_path: Path = get_last_output_path(),
    n_threads: int | None = None,
    **kwargs,
):
    pool = mp.Pool(n_threads)
    kwargs["intra_bounds"] = intra_bounds
    kwargs["extra_bounds"] = extra_bounds
    kwargs["output_path"] = output_path
    iterations = get_all_iterations(output_path)
    args = zip(
        iterations,
        itertools.repeat(kwargs),
    )
    print("Plotting Results")
    _ = list(
        tqdm.tqdm(pool.imap(__plot_all_iterations_helper, args), total=len(iterations))
    )


def generate_movie(opath: Path | None = None, play_movie: bool = True):
    if opath is None:
        opath = get_last_output_path(opath)
    bashcmd = f"ffmpeg\
        -v quiet\
        -stats\
        -y\
        -r 30\
        -f image2\
        -pattern_type glob\
        -i '{opath}/images/*.png'\
        -c:v h264\
        -pix_fmt yuv420p\
        -strict -2 {opath}/movie.mp4"
    os.system(bashcmd)

    if play_movie:
        print("Playing Movie")
        bashcmd2 = f"firefox ./{opath}/movie.mp4"
        os.system(bashcmd2)


if __name__ == "__main__":
    output_path = get_last_output_path()
    plot_all_iterations(
        (0.5, 1),
        (0, 10.0),
        output_path,
    )

    generate_movie(output_path)

ℹ️

We are currently working on a rewrite of this example using the chili backend under cellular_raza-examples/bacterial_branching.

References

[1] K. Kawasaki, A. Mochizuki, M. Matsushita, T. Umeda, and N. Shigesada, “Modeling Spatio-Temporal Patterns Generated byBacillus subtilis,” Journal of Theoretical Biology, vol. 188, no. 2. Elsevier BV, pp. 177–185, Sep. 1997. doi: 10.1006/jtbi.1997.0462.

[2] M. Matsushita et al., “Interface growth and pattern formation in bacterial colonies,” Physica A: Statistical Mechanics and its Applications, vol. 249, no. 1–4. Elsevier BV, pp. 517–524, Jan. 1998. doi: 10.1016/s0378-4371(97)00511-6.