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 = { version = "0.16.0", 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 {
        [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 = [self.growth_rate * uptake[0]].into();
        let incr_extra = -uptake;

        Ok((incr_intra, incr_extra))
    }
}
Domain
cellular_raza-examples/bacterial_branching/src/subdomain.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
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393

pub const N_REACTIONS: usize = 1;
pub type ReactionVector = nalgebra::SVector<f32, N_REACTIONS>;

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

#[derive(Domain, Clone)]
pub struct CartesianDiffusion2D {
    #[DomainRngSeed]
    #[DomainPartialDerive]
    #[SortCells]
    pub domain: CartesianCuboid<f32, 2>,
    /// The discretization must be a multiple of the voxel size.
    /// This quantity will be used as an initial estimate and rounded to the nearest candidate.
    pub reactions_dx: nalgebra::Vector2<f32>,
    pub diffusion_constant: f32,
    pub initial_value: ReactionVector,
}

#[derive(Deserialize, SubDomain, Clone, Serialize, Debug)]
pub struct CartesianDiffusion2DSubDomain {
    #[Base]
    #[SortCells]
    #[Mechanics]
    pub subdomain: CartesianSubDomain<f32, 2>,
    pub reactions_min: nalgebra::Vector2<f32>,
    pub reactions_max: nalgebra::Vector2<f32>,
    pub reactions_dx: nalgebra::Vector2<f32>,
    pub index_min: nalgebra::Vector2<usize>,
    pub index_max: nalgebra::Vector2<usize>,
    pub extracellular: ndarray::Array3<f32>,
    pub ownership_array: ndarray::Array2<bool>,
    pub diffusion_constant: f32,
    increments: [ndarray::Array3<f32>; 3],
    increments_start: usize,
    helper: ndarray::Array3<f32>,
}

pub struct BorderInfo {
    min_sent: nalgebra::Vector2<usize>,
    max_sent: nalgebra::Vector2<usize>,
}

impl CartesianDiffusion2DSubDomain {
    fn assign_neighbor(&mut self, neighbor: NeighborValue) {
        use ndarray::*;
        let NeighborValue { min, max, values } = neighbor;

        // Cast everything to isize to avoid overflows
        let min = min.cast::<isize>();
        let max = max.cast::<isize>();
        let index_min = self.index_min.cast::<isize>();
        let index_max = self.index_max.cast::<isize>();

        // Legend
        // o = own array
        // x = neighbor value
        //
        // x x x x x x . . . .
        // x x x x x x . . . .
        // x x x o o o o o o o -- shared_min[1]
        // x x x o o o o o o o
        // x x x o o o o o o o -- shared_max[1]
        // . . . o o o o o o o
        // . . . o o o o o o o
        //       |   |
        //       |   shared_max[0]
        //       |
        //       shared_min[0]
        let shared_min = min.sup(&index_min);
        let shared_max = max.inf(&index_max);

        // x x x x x x . . . .                   o o o o o o . . . .
        // x x h h h h . . . . -- helper_min[1]  o o o o o o . . . .
        // x x h o o o o o o o                   o o o o o o h x x x -- helper_min[1]
        // x x h o o o o o o o                   o o o o o o h x x x
        // x x h o o o o o o o -- helper_max[1]  o o o o o o h x x x
        // . . . o o o o o o o                   . . . h h h h x x x -- helper_max[1]
        // . . . o o o o o o o                   . . . x x x x x x x
        //     |     |                                 |     |
        //     |     helper_max[0]                     |     helper_max[0]
        //     |                                       |
        //     heper_min[0]                            helper_min[0]
        let helper_min = shared_min.add_scalar(-1).sup(&min);
        let helper_max = shared_max.add_scalar(1).inf(&max);

        let nmin = helper_min - min;
        let nmax = helper_max - min;
        let hmin = (helper_min - index_min).add_scalar(1);
        // let hmax = (helper_max - index_min).inf(&shared_max).add_scalar(1);
        let hmax = nmax - nmin + hmin;

        Zip::from(
            self.helper
                .slice_mut(s![hmin[0]..hmax[0], hmin[1]..hmax[1], ..])
                .lanes_mut(Axis(2)),
        )
        .and(values.slice(s![nmin[0]..nmax[0], nmin[1]..nmax[1]]))
        .and(
            self.ownership_array
                .slice(s![hmin[0]..hmax[0], hmin[1]..hmax[1]]),
        )
        .for_each(|mut w, v, t| {
            if let (false, Some(vi)) = (*t, v) {
                w.assign(&vi);
            }
        });
    }
}

pub struct NeighborValue {
    min: nalgebra::Vector2<usize>,
    max: nalgebra::Vector2<usize>,
    values: ndarray::Array2<Option<ndarray::Array1<f32>>>,
}

impl SubDomainReactions<nalgebra::SVector<f32, 2>, ReactionVector, f32>
    for CartesianDiffusion2DSubDomain
{
    type BorderInfo = BorderInfo;
    type NeighborValue = NeighborValue;

    fn treat_increments<I, J>(
        &mut self,
        neighbors: I,
        sources: J,
    ) -> Result<(), cellular_raza::concepts::CalcError>
    where
        I: IntoIterator<Item = Self::NeighborValue>,
        J: IntoIterator<Item = (nalgebra::SVector<f32, 2>, ReactionVector)>,
    {
        use core::ops::AddAssign;
        use ndarray::*;
        let dx2 = self.reactions_dx[0].powf(-2.0);
        let dy2 = self.reactions_dx[1].powf(-2.0);
        let dd2 = -2.0 * (dx2 + dy2);

        // Helper variable to store current concentrations
        let co = &self.extracellular;

        // Use helper array which is +2 in every spatial dimension larger than the original array
        // We do this to seamlessly incorporate boundary conditions
        // Fill inner part of the array
        self.helper.fill(0.0);
        // _ _ _ _ _ _ _
        // _ x x x x x _
        // _ x x x x x _
        // _ _ _ _ _ _ _
        self.helper.slice_mut(s![1..-1, 1..-1, ..]).assign(&co);

        // First assume that we obey neumann-boundary conditions and fill outer parts accordingly
        // _ x x x x x _
        // x _ _ _ _ _ x
        // x _ _ _ _ _ x
        // _ x x x x x _
        self.helper
            .slice_mut(s![0, 1..-1, ..])
            .assign(&co.slice(s![0, .., ..]));
        self.helper
            .slice_mut(s![-1, 1..-1, ..])
            .assign(&co.slice(s![-1, .., ..]));
        self.helper
            .slice_mut(s![1..-1, 0, ..])
            .assign(&co.slice(s![.., 0, ..]));
        self.helper
            .slice_mut(s![1..-1, -1, ..])
            .assign(&co.slice(s![.., -1, ..]));

        for neighbor in neighbors {
            self.assign_neighbor(neighbor);
        }

        // Set increment to next time-step to 0.0 everywhere
        let start = self.increments_start;
        self.increments[start].fill(0.0);

        let dc = self.diffusion_constant;
        // - 2u[i,j] /dx^2 - 2u[i,j]/dy^2
        self.increments[start].add_assign(&(dd2 * dc * &self.helper.slice(s![1..-1, 1..-1, ..])));
        // + u[i-1,j]/dx^2
        self.increments[start].add_assign(&(dx2 * dc * &self.helper.slice(s![..-2, 1..-1, ..])));
        // + u[i+1,j]/dx^2
        self.increments[start].add_assign(&(dx2 * dc * &self.helper.slice(s![2.., 1..-1, ..])));
        // + u[i,j-1]/dy^2
        self.increments[start].add_assign(&(dy2 * dc * &self.helper.slice(s![1..-1, ..-2, ..])));
        // + u[i,j+1]/dy^2
        self.increments[start].add_assign(&(dy2 * dc * &self.helper.slice(s![1..-1, 2.., ..])));

        for (pos, dextra) in sources {
            let index = self.get_extracellular_index(&pos)?;
            let dextra: [f32; N_REACTIONS] = dextra.into();
            self.increments[start]
                .slice_mut(ndarray::s![index[0], index[1], ..])
                .scaled_add(1.0, &ndarray::Array1::<f32>::from_iter(dextra));
        }

        Ok(())
    }

    fn update_fluid_dynamics(&mut self, dt: f32) -> Result<(), cellular_raza::concepts::CalcError> {
        use core::ops::AddAssign;
        let start = self.increments_start;
        let n_incr = self.increments.len();

        // Adams-Bashforth 3rd order
        let k1 = 5.0 / 12.0;
        let k2 = 8.0 / 12.0;
        let k3 = -1.0 / 12.0;
        self.extracellular.add_assign(
            &(k1 * dt * &self.increments[start]
                + k2 * dt * &self.increments[(start + 1) % n_incr]
                + k3 * dt * &self.increments[(start + 2) % n_incr]),
        );
        self.extracellular.map_inplace(|x| *x = x.max(0.0));

        // TODO DEBUGGING
        // if start == 0 {
        //     todo!();
        // }

        self.increments_start = (self.increments_start + 1) % n_incr;
        Ok(())
    }

    fn get_extracellular_at_pos(
        &self,
        pos: &nalgebra::SVector<f32, 2>,
    ) -> Result<ReactionVector, cellular_raza::concepts::CalcError> {
        let index = self.get_extracellular_index(pos)?;
        let res = ReactionVector::from_iterator(
            self.extracellular
                .slice(ndarray::s![index[0], index[1], ..])
                .to_owned()
                .into_iter(),
        );
        Ok(res)
    }

    fn get_neighbor_value(&self, border_info: Self::BorderInfo) -> Self::NeighborValue {
        use ndarray::*;
        // Calculate shared indices plus padding of one
        let BorderInfo { min_sent, max_sent } = border_info;
        let min = min_sent.map(|x| x.saturating_sub(1)).sup(&self.index_min);
        let max = max_sent.map(|x| x.saturating_add(1)).inf(&self.index_max);

        let omin = min - self.index_min;
        let omax = max - self.index_min;
        let values = Zip::from(
            self.extracellular
                .slice(s![omin[0]..omax[0], omin[1]..omax[1], ..])
                .lanes(Axis(2)),
        )
        .and(
            self.ownership_array
                .slice(s![omin[0] + 1..omax[0] + 1, omin[1] + 1..omax[1] + 1]),
        )
        .map_collect(|v, &o| if o { Some(v.to_owned()) } else { None })
        .to_owned();

        NeighborValue { min, max, values }
    }

    fn get_border_info(&self) -> Self::BorderInfo {
        Self::BorderInfo {
            min_sent: self.index_min,
            max_sent: self.index_max,
        }
    }
}

impl CartesianDiffusion2DSubDomain {
    fn get_extracellular_index(
        &self,
        pos: &nalgebra::Vector2<f32>,
    ) -> Result<nalgebra::Vector2<usize>, CalcError> {
        let index = (pos - self.reactions_min).component_div(&self.reactions_dx);
        if index
            .iter()
            .enumerate()
            .any(|(n, &x)| x < 0.0 || x > self.index_max[n] as f32)
        {
            return Err(CalcError(format!(
                "Could not find index for position {:?}",
                pos
            )));
        }
        Ok(index.map(|x| x.floor() as usize))
    }
}

impl DomainCreateSubDomains<CartesianDiffusion2DSubDomain> for CartesianDiffusion2D {
    type SubDomainIndex = usize;
    type VoxelIndex = [usize; 2];

    fn create_subdomains(
        &self,
        n_subdomains: core::num::NonZeroUsize,
    ) -> Result<
        impl IntoIterator<
            Item = (
                Self::SubDomainIndex,
                CartesianDiffusion2DSubDomain,
                Vec<Self::VoxelIndex>,
            ),
        >,
        DecomposeError,
    > {
        let dx = self.reactions_dx;
        let dx_domain = self.domain.get_dx();
        let n_diffusion = dx_domain
            .component_div(&dx)
            .map(|x| (x.round() as usize).max(1));
        let dx = dx_domain.component_div(&n_diffusion.cast::<f32>());

        let diffusion_constant = self.diffusion_constant;

        // Calculate lattice points for subdomains
        // Introduce padding on the outside of the simulated domain.
        let [nrows, ncols] = n_diffusion
            .component_mul(&self.domain.get_n_voxels())
            .into();

        let extracellular_total =
            ndarray::Array3::from_shape_fn((nrows, ncols, N_REACTIONS), |(_, _, n)| {
                self.initial_value[n]
            });

        Ok(self
            .domain
            .create_subdomains(n_subdomains)?
            .into_iter()
            .map(move |(index, subdomain, voxels)| {
                let max_domain = [nrows, ncols].into();
                let mut min: nalgebra::Vector2<usize> = max_domain;
                let mut max: nalgebra::Vector2<usize> = [0; 2].into();
                for vox in subdomain.get_voxels() {
                    min = min.inf(&vox.into());
                    max = max.sup(&vox.into());
                }

                // Multiply with number of voxels in each dimension
                let min = min.component_mul(&n_diffusion);
                // Here we need to add one more step since this is supposed to be an upper limit
                let max = max.component_mul(&n_diffusion) + n_diffusion;
                let max_domain = max_domain.component_mul(&n_diffusion);

                let extracellular = extracellular_total
                    .slice(ndarray::s![min[0]..max[0], min[1]..max[1], ..])
                    .into_owned();

                let reactions_min = min.cast::<f32>().component_mul(&dx);
                let reactions_max = max.cast::<f32>().component_mul(&dx);

                // Has entry `true` if the given point is owned by this subdomain.
                let d = extracellular.dim();
                let mut ownership_array =
                    ndarray::Array2::<bool>::from_elem((d.0 + 2, d.1 + 2), false);
                for v in subdomain.get_voxels() {
                    let one = nalgebra::Vector2::from([1; 2]);
                    let v: nalgebra::Vector2<usize> = v.into();
                    let vox = v.component_mul(&n_diffusion);
                    let voxp1 = (v + one).component_mul(&n_diffusion);
                    let lower = (vox - min).add_scalar(1);
                    let upper = (voxp1 - min).inf(&max_domain).add_scalar(1);
                    ownership_array
                        .slice_mut(ndarray::s![lower[0]..upper[0], lower[1]..upper[1]])
                        .fill(true);
                }

                let sh = extracellular.shape();
                let increment = ndarray::Array3::zeros((sh[0], sh[1], sh[2]));
                let helper = ndarray::Array3::zeros((sh[0] + 2, sh[1] + 2, sh[2]));
                (
                    index,
                    CartesianDiffusion2DSubDomain {
                        subdomain,
                        reactions_min,
                        reactions_max,
                        reactions_dx: dx,
                        index_min: min,
                        index_max: max,
                        extracellular,
                        ownership_array,
                        diffusion_constant,
                        increments_start: 0,
                        increments: [increment.clone(), increment.clone(), increment.clone()],
                        helper,
                    },
                    voxels,
                )
            }))
    }
}
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

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;
mod subdomain;

use bacteria_properties::*;
use subdomain::*;

#[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([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,
    )?;
    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.