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 described by a system of PDEs (Partial Differential Equations) which describe intracellular reactions, growth and spatial contributions (typically via Diffusion processes) which describe diffusion of nutrients and movement of the cells. We consider two variables: the spatial-temporal distribution $n$ of the available nutrient and the bacterial population density $b$. Together with the ratio of diffusion of both components $d$ the 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. 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

Mechanics & Interaction

We represent cells as soft spheres with the same interaction potential as in the cell-sorting example but omitting the species-specificity.

Intra- & extracellular Reactions

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

$$\begin{align} \dot{n}_i^c &= u(n_e(x_c) - n_i^c) - (\alpha + \sigma)n_i^c\\ \dot{V}_c &= \alpha n_i^c V_c \left(1 - \frac{V_c}{V_\text{max}}\right)\\ \dot{n}_e &= D\Delta n_e - u \sum\limits_{c=1}^N \delta(x-x_c)(n_e-n_i) \end{align}$$

The components of these PDEs describe the extracellular and intracellular nutrients as well as the volume of the individual cell. $n_e$ is the spatially distributed extra-cellular nutrient concentration which undergoes diffusion with the diffusion constant $D$ while $n^c_i$ is the intra-cellular nutrient concentration of cell $c\in\{1,\dots,N\}$ positioned at $x_c$. The volumeo f cell $c$ is given by $V_c$. The parameter $u$ is the uptake rate of the nutrient while $\alpha$ describes the consumption of the nutrient by the cellular metabolism resulting in an increase of the volume $V_c$. In contrast, $\sigma$ degrades the intracellular nutrients «««< HEAD After reaching $90\%$ of the maximum volume $V_\text{max}$ the cell divides into two equally sized cells with volume $V_\text{max}/2$. The gain ΞΈ is approximately given by $\theta\approx \log(2) \alpha u/ (u + \alpha)$. If the uptake rate is very large the gain is limited by the nutrient processing rate $\alpha$ (the nutrient metabolism), but in the case where the nutrient metabolism is very fast (efficient) the total gain is limited by the uptake rate $u$.

Cycle

Once cells have reached a minimum age $\tau$ and nutrient threshold $n_t$, they will divide. The newly created agents behave exactly the same as their parent and they will continue to take up nutrients, process them and divide. This also means, we do not alter their internal parameters during the division process. Usually, this cycle of producing new generations ceases due to low concentration of the nutrient resource after a few division events. Our simulation duration is chosen short enough such that the colony remains in its equilibrium state and we thus do not need to model cell death.

Parameters

Parameter Symbol Value
Cell Radius $R$ $6.0\hspace{0.25em}\mu \text{m}$
Potential Strength $V_0$ $2\hspace{0.25em}\mu\text{m}^2\hspace{0.25em}\text{min}^{-2}$
Damping Constant $\lambda$ $2\hspace{0.25em}\text{min}^{-1}$
Interaction Range $\xi$ $1.5 \hspace{0.25em} R$
Turnover Rate $\sigma$ $0.025 \hspace{0.25em}\text{min}^{-1}$
Uptake Rate $u$ $0.05 \hspace{0.25em}\text{min}^{-1}$
Growth Rate $\alpha$ $0.1\hspace{0.25em}\mu m^2\hspace{0.25em}\mu g^{-1}\hspace{0.25em}l$
Nutrient Division Threshold $n_t$ $0.8\hspace{0.25em}\mu g\hspace{0.25em} l^{-1}$
Age threshold $\tau$ $65\hspace{0.25em}\text{min}$
Diffusion Constant $D$ $12\hspace{0.25em}\mu m^2 \hspace{0.25em}\text{min}^{-1}$

Initial State

Property Symbol Value
Time Stepsize $\Delta t$ $0.25\hspace{0.25em}\text{min}$
Time Steps $N_t$ $20'000$
Domain Size $L$ $3000\hspace{0.25em}\mu\text{m}$
Centered Starting Domain Size $L_0$ $300\hspace{0.25em}\mu m$
Number of cells $N_0$ $400$
Initial Intracellular Nutrients $n_i^c$ $1.0 \hspace{0.25em}\mu g\hspace{0.25em} \hspace{0.25em}l^{-1}$
Initial Extracellular Nutrients $n_e$ $25\hspace{0.25em}\mu g\hspace{0.25em} l^{-1}$

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\approx 30000\hspace{0.25em}\mu m$ and with an increased number of simulation steps. Unfortunately the exact numbers have been lost.

Code

⚠️
The code for this example is implemented via the cpu_os_threads backend. This backend is not receiving new updates or further development. We are currently working at porting the existing code to the chili backend.

Bacterial Properties
cellular_raza-examples/bacteria_population/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
 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

use cellular_raza::concepts::reactions_old::*;
use cellular_raza::prelude::*;

use nalgebra::Vector2;
use num::Zero;
use rand::Rng;
use serde::{Deserialize, Serialize};

pub const NUMBER_OF_REACTION_COMPONENTS: usize = 1;
pub type ReactionVector = nalgebra::SVector<f64, NUMBER_OF_REACTION_COMPONENTS>;
pub type MyCellType =
    ModularCell<NewtonDamped2D, CellSpecificInteraction, OwnCycle, OwnReactions, GradientSensing>;

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

impl Interaction<Vector2<f64>, Vector2<f64>, Vector2<f64>, f64> for CellSpecificInteraction {
    fn get_interaction_information(&self) -> f64 {
        self.cell_radius
    }

    fn calculate_force_between(
        &self,
        own_pos: &Vector2<f64>,
        _own_vel: &Vector2<f64>,
        ext_pos: &Vector2<f64>,
        _ext_vel: &Vector2<f64>,
        ext_radius: &f64,
    ) -> Result<(Vector2<f64>, Vector2<f64>), CalcError> {
        let min_relative_distance_to_center = 0.3162277660168379;
        let (r, dir) =
            match (own_pos - ext_pos).norm() < self.cell_radius * min_relative_distance_to_center {
                false => {
                    let z = own_pos - ext_pos;
                    let r = z.norm();
                    (r, z.normalize())
                }
                true => {
                    let dir = match own_pos == ext_pos {
                        true => {
                            return Ok((Vector2::zeros(), Vector2::zeros()));
                        }
                        false => (own_pos - ext_pos).normalize(),
                    };
                    let r = self.cell_radius * min_relative_distance_to_center;
                    (r, dir)
                }
            };
        // Introduce Non-dimensional length variable
        let sigma = r / (self.cell_radius + ext_radius);
        let bound = 4.0 + 1.0 / sigma;
        let spatial_cutoff = (1.0
            + (self.relative_interaction_range * (self.cell_radius + ext_radius) - r).signum())
            * 0.5;

        // Calculate the strength of the interaction with correct bounds
        let strength = self.potential_strength
            * ((1.0 / sigma).powf(2.0) - (1.0 / sigma).powf(4.0))
                .min(bound)
                .max(-bound);

        // Calculate only attracting and repelling forces
        let attracting_force = dir * strength.max(0.0) * spatial_cutoff;
        let repelling_force = dir * strength.min(0.0) * spatial_cutoff;

        Ok((
            -repelling_force - attracting_force,
            repelling_force + attracting_force,
        ))
    }
}

#[derive(Serialize, Deserialize, Debug, Clone)]
pub struct OwnCycle {
    age: f64,
    pub division_age: f64,
    pub maximum_cell_radius: f64,
    pub growth_rate: f64,
    pub food_threshold: f64,
    food_growth_rate_multiplier: f64,
    food_division_threshold: f64,
}

impl OwnCycle {
    pub fn new(
        division_age: f64,
        maximum_cell_radius: f64,
        growth_rate: f64,
        food_threshold: f64,
        food_growth_rate_multiplier: f64,
        food_division_threshold: f64,
    ) -> Self {
        OwnCycle {
            age: 0.0,
            division_age,
            maximum_cell_radius,
            growth_rate,
            food_threshold,
            food_growth_rate_multiplier,
            food_division_threshold,
        }
    }
}

impl Cycle<MyCellType> for OwnCycle {
    fn update_cycle(
        rng: &mut rand_chacha::ChaCha8Rng,
        dt: &f64,
        cell: &mut MyCellType,
    ) -> Option<CycleEvent> {
        // If the cell is not at the maximum size let it grow
        if cell.interaction.cell_radius < cell.cycle.maximum_cell_radius {
            let growth_difference = (cell.cycle.maximum_cell_radius * cell.cycle.growth_rate * dt)
                .min(cell.cycle.maximum_cell_radius - cell.interaction.cell_radius);
            cell.cellular_reactions.intracellular_concentrations[0] -=
                cell.cycle.food_growth_rate_multiplier * growth_difference
                    / cell.cycle.maximum_cell_radius;
            cell.interaction.cell_radius += growth_difference;
        }

        // Increase the age of the cell and divide if possible
        cell.cycle.age += dt;

        // Calculate the modifier (between 0.0 and 1.0) based on food threshold
        let relative_division_food_level = ((cell.get_intracellular()[0]
            - cell.cycle.food_division_threshold)
            / (cell.cycle.food_threshold - cell.cycle.food_division_threshold))
            .clamp(0.0, 1.0);

        if
        // Check if the cell has aged enough
        cell.cycle.age > cell.cycle.division_age &&
            // Check if the cell has grown enough
            cell.interaction.cell_radius >= cell.cycle.maximum_cell_radius &&
            // Random selection but chance increased when significantly above the food threshold
            rng.gen_range(0.0..1.0) < relative_division_food_level
        {
            return Some(CycleEvent::Division);
        }
        None
    }

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

        let r = c1.interaction.cell_radius;

        // Make both cells smaller
        // ALso keep old cell larger
        let relative_size_difference = 0.2;
        c1.interaction.cell_radius *= (1.0 + relative_size_difference) / std::f64::consts::SQRT_2;
        c2.interaction.cell_radius *= (1.0 - relative_size_difference) / std::f64::consts::SQRT_2;

        // Generate cellular splitting direction randomly
        let angle_1 = rng.gen_range(0.0..2.0 * std::f64::consts::PI);
        let dir_vec = nalgebra::Rotation2::new(angle_1) * Vector2::from([1.0, 0.0]);

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

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

        // Decrease the amount of food in the cells
        c1.cellular_reactions.intracellular_concentrations *=
            (1.0 + relative_size_difference) * 0.5;
        c2.cellular_reactions.intracellular_concentrations *=
            (1.0 - relative_size_difference) * 0.5;

        // New cell is completely new so set age to 0
        c2.cycle.age = 0.0;

        Ok(c2)
    }
}

#[derive(Serialize, Deserialize, Clone, Debug)]
pub struct OwnReactions {
    pub intracellular_concentrations: ReactionVector,
    pub turnover_rate: ReactionVector,
    pub production_term: ReactionVector,
    pub degradation_rate: ReactionVector,
    pub secretion_rate: ReactionVector,
    pub uptake_rate: ReactionVector,
}

impl CellularReactions<ReactionVector> for OwnReactions {
    fn calculate_intra_and_extracellular_reaction_increment(
        &self,
        internal_concentration_vector: &ReactionVector,
        external_concentration_vector: &ReactionVector,
    ) -> Result<(ReactionVector, ReactionVector), CalcError> {
        let mut increment_extracellular = ReactionVector::zero();
        let mut increment_intracellular = ReactionVector::zero();

        for i in 0..NUMBER_OF_REACTION_COMPONENTS {
            let uptake = self.uptake_rate[i] * external_concentration_vector[i];
            let secretion = self.secretion_rate[i] * internal_concentration_vector[i];
            increment_extracellular[i] = secretion - uptake;
            increment_intracellular[i] = self.production_term[i]
                - increment_extracellular[i]
                - self.turnover_rate[i] * internal_concentration_vector[i];
        }

        Ok((increment_intracellular, increment_extracellular))
    }

    fn get_intracellular(&self) -> ReactionVector {
        self.intracellular_concentrations
    }

    fn set_intracellular(&mut self, concentration_vector: ReactionVector) {
        self.intracellular_concentrations = concentration_vector;
    }
}

#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct GradientSensing {}

impl
    InteractionExtracellularGradient<
        MyCellType,
        nalgebra::SVector<Vector2<f64>, NUMBER_OF_REACTION_COMPONENTS>,
    > for GradientSensing
{
    fn sense_gradient(
        _cell: &mut MyCellType,
        _gradient: &nalgebra::SVector<Vector2<f64>, NUMBER_OF_REACTION_COMPONENTS>,
    ) -> Result<(), CalcError> {
        Ok(())
    }
}
Main Simulation
cellular_raza-examples/bacteria_population/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
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

use cellular_raza::building_blocks::{*, cartesian_cuboid_n_old::*};
use cellular_raza::concepts::domain_old::*;
use cellular_raza::concepts::*;
use cellular_raza::core::backend::cpu_os_threads::*;
use cellular_raza::core::storage::*;

use nalgebra::Vector2;

use num::Zero;

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

use serde::{Deserialize, Serialize};

// Number of cells to put into simulation in the Beginning
pub const N_BACTERIA_INITIAL: u32 = 400;

// Mechanical parameters
pub const BACTERIA_MECHANICS_RADIUS: f64 = 6.0;
pub const BACTERIA_MECHANICS_RELATIVE_INTERACTION_RANGE: f64 = 1.6;
pub const BACTERIA_MECHANICS_POTENTIAL_STRENGTH: f64 = 2.0;
pub const BACTERIA_MECHANICS_VELOCITY_REDUCTION: f64 = 2.0;

// Reaction parameters of the cell
pub const BACTERIA_FOOD_INITIAL_CONCENTRATION: f64 = 1.0;
pub const BACTERIA_FOOD_TURNOVER_RATE: f64 = 0.025;
pub const BACTERIA_FOOD_UPTAKE_RATE: f64 = 0.05;

// Parameters for cell cycle
pub const BACTERIA_CYCLE_DIVISION_AGE_MIN: f64 = 60.0;
pub const BACTERIA_CYCLE_DIVISION_AGE_MAX: f64 = 70.0;
pub const BACTERIA_CYCLE_GROWTH_RATE: f64 = 0.1;
pub const BACTERIA_CYCLE_FOOD_THRESHOLD: f64 = 2.0;
pub const BACTERIA_CYCLE_FOOD_GROWTH_RATE_MULTIPLIER: f64 = 10.0;
pub const BACTERIA_CYCLE_FOOD_DIVISION_THRESHOLD: f64 = BACTERIA_FOOD_INITIAL_CONCENTRATION * 0.8;

// Parameters for domain
pub const DOMAIN_SIZE: f64 = 3_000.0;
pub const DOMAIN_MIDDLE: Vector2<f64> = nalgebra::vector![DOMAIN_SIZE / 2.0, DOMAIN_SIZE / 2.0];

// Where will the cells be placed initially
pub const STARTING_DOMAIN_X_LOW: f64 = DOMAIN_SIZE / 2.0 - 150.0;
pub const STARTING_DOMAIN_X_HIGH: f64 = DOMAIN_SIZE / 2.0 + 150.0;
pub const STARTING_DOMAIN_Y_LOW: f64 = DOMAIN_SIZE / 2.0 - 150.0;
pub const STARTING_DOMAIN_Y_HIGH: f64 = DOMAIN_SIZE / 2.0 + 150.0;

// Parameters for Voxel Reaction+Diffusion
pub const VOXEL_FOOD_DIFFUSION_CONSTANT: f64 = 25.0;
pub const VOXEL_FOOD_INITIAL_CONCENTRATION: f64 = 12.0;

// Time parameters
pub const N_TIMES: usize = 20_001;
pub const DT: f64 = 0.25;
pub const T_START: f64 = 0.0;
pub const SAVE_INTERVAL: usize = 250;

// Meta Parameters to control solving
pub const N_THREADS: usize = 1;
pub const N_PLOTTING_THREADS: usize = 14;

mod bacteria_properties;
mod plotting;

use bacteria_properties::*;
use plotting::*;

fn voxel_definition_strategy(voxel: &mut CartesianCuboidVoxel2<NUMBER_OF_REACTION_COMPONENTS>) {
    voxel.diffusion_constant = ReactionVector::from([VOXEL_FOOD_DIFFUSION_CONSTANT]);
    voxel.extracellular_concentrations = ReactionVector::from([VOXEL_FOOD_INITIAL_CONCENTRATION]);
    voxel.degradation_rate = ReactionVector::zero();
    voxel.production_rate = ReactionVector::zero();
}

fn create_domain() -> Result<CartesianCuboid2, CalcError> {
    CartesianCuboid2::from_boundaries_and_interaction_ranges(
        [0.0; 2],
        [DOMAIN_SIZE, DOMAIN_SIZE],
        [BACTERIA_MECHANICS_RADIUS * BACTERIA_MECHANICS_RELATIVE_INTERACTION_RANGE; 2],
    )
}

#[derive(Clone, Serialize, Deserialize)]
pub struct CellNumberController {
    target_cell_number: i64,
    stored_ids: std::collections::HashSet<(u64, u64)>,
    full: bool,
}

type Observable = Option<(i64, Vec<(u64, u64)>)>;

impl Controller<MyCellType, Observable> for CellNumberController {
    fn measure<'a, I>(&self, cells: I) -> Result<Observable, CalcError>
    where
        I: IntoIterator<Item = &'a CellAgentBox<MyCellType>> + Clone,
    {
        if !self.full {
            let mut n_cells = 0_i64;
            let positions = cells
                .into_iter()
                .map(|c| {
                    n_cells += 1;
                    c.get_id()
                })
                .collect();
            Ok(Some((n_cells, positions)))
        } else {
            Ok(None)
        }
    }

    fn adjust<'a, 'b, I, J>(&mut self, measurements: I, cells: J) -> Result<(), ControllerError>
    where
        Observable: 'a,
        MyCellType: 'b,
        I: Iterator<Item = &'a Observable>,
        J: Iterator<Item = (&'b mut CellAgentBox<MyCellType>, &'b mut Vec<CycleEvent>)>,
    {
        // If we are not full, we
        if !self.full {
            let mut total_cell_number: i64 = 0;
            let all_ids: std::collections::HashSet<_> = measurements
                .into_iter()
                .filter_map(|meas| {
                    if let Some((n_cells, ids)) = meas {
                        total_cell_number += n_cells;
                        Some(ids.into_iter())
                    } else {
                        None
                    }
                })
                .flatten()
                .map(|&id| id)
                .collect();

            if total_cell_number > self.target_cell_number {
                self.stored_ids = all_ids;
                self.full = true;
            }
        }
        if self.full {
            // Kill all cells which do not match ids
            for (cell, cell_events) in cells.into_iter() {
                if !self.stored_ids.contains(&cell.get_id()) {
                    cell_events.push(CycleEvent::Remove);
                }
            }
        }
        Ok(())
    }
}

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

    // ###################################### DEFINE SIMULATION DOMAIN ######################################
    let domain = create_domain().unwrap();

    // ###################################### DEFINE CELLS IN SIMULATION ######################################
    let cells = (0..N_BACTERIA_INITIAL)
        .map(|_| {
            let x = rng.gen_range(STARTING_DOMAIN_X_LOW..STARTING_DOMAIN_X_HIGH);
            let y = rng.gen_range(STARTING_DOMAIN_Y_LOW..STARTING_DOMAIN_Y_HIGH);

            let pos = Vector2::from([x, y]);
            ModularCell {
                mechanics: NewtonDamped2D {
                    pos,
                    vel: Vector2::zero(),
                    damping_constant: BACTERIA_MECHANICS_VELOCITY_REDUCTION,
                    mass: 1.0,
                },
                interaction: CellSpecificInteraction {
                    potential_strength: BACTERIA_MECHANICS_POTENTIAL_STRENGTH,
                    relative_interaction_range: BACTERIA_MECHANICS_RELATIVE_INTERACTION_RANGE,
                    cell_radius: BACTERIA_MECHANICS_RADIUS,
                },
                interaction_extracellular: GradientSensing {},
                cycle: OwnCycle::new(
                    rng.gen_range(0.0..BACTERIA_CYCLE_DIVISION_AGE_MAX),
                    BACTERIA_MECHANICS_RADIUS,
                    BACTERIA_CYCLE_GROWTH_RATE,
                    BACTERIA_CYCLE_FOOD_THRESHOLD,
                    BACTERIA_CYCLE_FOOD_GROWTH_RATE_MULTIPLIER,
                    BACTERIA_CYCLE_FOOD_DIVISION_THRESHOLD,
                ),
                cellular_reactions: OwnReactions {
                    intracellular_concentrations: ReactionVector::from([
                        BACTERIA_FOOD_INITIAL_CONCENTRATION,
                    ]),
                    turnover_rate: ReactionVector::from([BACTERIA_FOOD_TURNOVER_RATE]),
                    production_term: ReactionVector::zero(),
                    degradation_rate: ReactionVector::zero(),
                    secretion_rate: ReactionVector::zero(),
                    uptake_rate: ReactionVector::from([BACTERIA_FOOD_UPTAKE_RATE]),
                },
                volume: (DOMAIN_SIZE / 151.0).powf(2.0), //2.0*std::f64::consts::PI*BACTERIA_MECHANICS_RADIUS.powf(2.0),
            }
        })
        .collect::<Vec<_>>();

    // ###################################### CREATE SUPERVISOR AND RUN SIMULATION ######################################
    let setup = SimulationSetup::new(
        domain,
        cells,
        TimeSetup {
            t_start: 0.0,
            t_eval: (0..N_TIMES)
                .map(|i| (T_START + DT * i as f64, i % SAVE_INTERVAL == 0))
                .collect::<Vec<(f64, bool)>>(),
        },
        SimulationMetaParams {
            n_threads: N_THREADS,
            ..Default::default()
        },
        StorageBuilder::new()
            .location("out/bacteria_population")
            .init(),
        CellNumberController {
            target_cell_number: 15_000,
            stored_ids: std::collections::HashSet::new(),
            full: false,
        },
    );

    let strategies = Strategies {
        voxel_definition_strategies: &voxel_definition_strategy,
    };

    let mut supervisor = SimulationSupervisor::initialize_with_strategies(setup, strategies);

    let mut simulation_result = supervisor.run_full_sim().unwrap();

    // ###################################### PLOT THE RESULTS ######################################
    simulation_result.plotting_config = PlottingConfig {
        n_threads: Some(N_PLOTTING_THREADS),
        image_size: 1500,
        image_type: ImageType::BitMap,
        ..Default::default()
    };

    simulation_result
        .plot_spatial_all_iterations_custom_cell_voxel_functions(plot_modular_cell, plot_voxel)
        .unwrap();
}
Plotting
cellular_raza-examples/bacteria_population/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

use cellular_raza::building_blocks::cartesian_cuboid_n_old::*;
use cellular_raza::concepts::domain_old::ExtracellularMechanics;
use cellular_raza::concepts::reactions_old::*;
use cellular_raza::prelude::*;

use plotters::{
    backend::BitMapBackend,
    coord::types::RangedCoordf64,
    prelude::{Cartesian2d, Circle, DrawingArea, ShapeStyle},
    style::colors::colormaps::{ColorMap, DerivedColorMap, ViridisRGB},
    style::RGBColor,
};

use crate::bacteria_properties::*;

pub fn plot_voxel(
    voxel: &CartesianCuboidVoxel2<NUMBER_OF_REACTION_COMPONENTS>,
    root: &mut DrawingArea<BitMapBackend, Cartesian2d<RangedCoordf64, RangedCoordf64>>,
) -> Result<(), DrawingError> {
    // Define lower and upper bounds for our values
    let lower_bound = 0.0;
    let upper_bound = 12.0;
    let concentration = voxel.get_total_extracellular()[0];

    // This should give a nice colormap
    let voxel_color = ViridisRGB::get_color_normalized(concentration, lower_bound, upper_bound);
    let rectangle = plotters::prelude::Rectangle::new(
        [
            (voxel.get_min()[0], voxel.get_min()[1]),
            (voxel.get_max()[0], voxel.get_max()[1]),
        ],
        Into::<ShapeStyle>::into(&voxel_color).filled(),
    );
    root.draw(&rectangle)?;
    Ok(())
}

pub fn plot_modular_cell(
    modular_cell: &MyCellType,
    root: &mut DrawingArea<BitMapBackend, Cartesian2d<RangedCoordf64, RangedCoordf64>>,
) -> Result<(), DrawingError> {
    let cell_border_color = plotters::prelude::BLACK;

    let relative_border_thickness = 0.25;

    // Plot the cell border
    let dx = root.get_x_range().end - root.get_x_range().start;
    let dx_pix = root.get_x_axis_pixel_range().end - root.get_x_axis_pixel_range().start;

    let s = modular_cell.interaction.cell_radius / dx * dx_pix as f64;
    let cell_border = Circle::new(
        (modular_cell.pos().x, modular_cell.pos().y),
        s,
        Into::<ShapeStyle>::into(&cell_border_color).filled(),
    );
    root.draw(&cell_border)?;

    let lower_bound = 0.0;
    let upper_bound = modular_cell.cycle.food_threshold / 10.0;

    // Define colormap
    let derived_colormap = DerivedColorMap::new(&[RGBColor(102, 52, 83), RGBColor(247, 126, 201)]);

    // Plot the inside of the cell
    let cell_inside_color = derived_colormap.get_color_normalized(
        modular_cell.get_intracellular()[0],
        lower_bound,
        upper_bound,
    );

    // let cell_inside_color = Life::get_color_normalized(modular_cell.get_intracellular()[1], 0.0, modular_cell.cellular_reactions.intracellular_concentrations_saturation_level[1]);
    let cell_inside = Circle::new(
        (modular_cell.pos().x, modular_cell.pos().y),
        s * (1.0 - relative_border_thickness),
        Into::<ShapeStyle>::into(&cell_inside_color).filled(),
    );
    root.draw(&cell_inside)?;
    Ok(())
}

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.