Bacterial Rods
Bacteria exist in many physical shapes such as spheroidal, rod-shaped or spiral [1,2]. We set out to define a simple yet effective coarse-grained mechanical cellular model to describe the motion of rod-shaped bacteria.
Mathematical Description
To model the spatial mechanics of elongated bacteria [3], we represent them as a collection of auxiliary vertices $\{\vec{v}_i\}$ which are connected by springs in ascending order. Furthermore, we assume that the cells are flexible which is described by their curvature property. A force $\vec{F}$ interacting between cellular agents determines the radius (thickness) of the rods and an attractive component can model adhesion between cells.
Mechanics
The internal force acting on vertex $\vec{v}_i$ can be divided into 2 contributions coming from the 2 springs pulling on it. We assume that the segment lengths $l$ and tension $\gamma$ are spatially independent. In the case when $i=0,N_\text{vertices}$, this is reduced to only one internal component. We denote with $\vec{c}_{i}$ the edge between two vertices
$$\begin{align} \vec{c}_i = \vec{v}_{i}-\vec{v}_{i-1} \end{align}$$
and can write down the resulting force
$$\begin{align} \vec{F}_{i,\text{springs}} = &-\gamma\left(1 - \frac{l}{\left|\vec{c}_i\right|}\right) \vec{c}_i\\ &+ \gamma\left(1 - \frac{l}{\left|\vec{c}_{i+1}\right|}\right) \vec{c}_{i+1} \end{align}$$
In addition to springs between individual vertices $\vec{v}_i$, we assume that each angle at a vertex between edges is subject to a force indiced by curvature. We assume that we can model the mechanical properties of the bacterium as an elastic rod. We define $\alpha_i = \sphericalangle(\vec{c}_{i-1},\vec{c}_i)$ as the angle between adjacent edges. The bending force can be assumed to be proportional to the curvature $\kappa_i$ at each vertex $\vec{v}_i$
$$\begin{equation} \kappa_i = 2\tan\left(\frac{\alpha_i}{2}\right). \end{equation}$$
The resulting force acts along the angle bisector which can be calculated from the edge vectors. The forces acting on vertices $\vec{v}_i,\vec{v}_{i-1},\vec{v}_{i+1}$ are given by
$$\begin{align} \vec{F}_{i,\text{curvature}} &= \eta\kappa_i \frac{\vec{c}_i - \vec{c}_{i+1}}{|\vec{c}_i-\vec{c}_{i+1}|}\\ \vec{F}_{i-1,\text{curvature}} &= -\frac{1}{2}\vec{F}_{i,\text{curvature}}\\ \vec{F}_{i+1,\text{curvature}} &= -\frac{1}{2}\vec{F}_{i,\text{curvature}} \end{align}$$
where $\eta_i$ is the rigidity at vertex $\vec{v}_i$ (see Figure 1). We can see that the curvature force does not move the overall center of the rod in space. The total force is the sum of external and interal forces.
$$\begin{equation} \vec{F}_{i,\text{total}} = \vec{F}_{i,\text{springs}}+ \vec{F}_{i,\text{curvature}} + \vec{F}_{i,\text{external}} \end{equation}$$
and are integrated via
$$\begin{align} \partial_t^2 \vec{x} &= \partial_t\vec{x} + \sqrt{2D}\vec{\xi}\\ \partial_t\vec{x} &= \vec{F}_\text{total} \end{align}$$
where $D$ is the diffusion constant and $\vec{\xi}$ is the wiener process (compare with
brownian motion).
These mechanics are available under the
RodMechanics building block.
Figure 1: One cell-agent consists of multiple vertices which are connected by springs.
Any angle which deviates from $\pi$ between the adjacent edges introduces a curvature force.
The area indicated around the polygon is the collections of points with distance less than $R$
to some point $\vec{p}$ on the segments.
It visualizes the interaction range.
Interaction
When calculating forces acting between the cells, we can use a simplified model to circumvent the numerically expensive integration over the complete length of the rod. Given a vertex $\vec{v}_i$ on one cell, we calculate the closest point $\vec{p}$ on the polygonal line given by the vertices $\{\vec{w}_j\}$ of the interacting cell. Furthermore we determine the value $q\in[0,1]$ such that
$$\begin{equation} \vec{p} = (1-q)\vec{w}_j + q\vec{w}_{j+1} \end{equation}$$
for some specific $j$. The force is then calculated between the points $\vec{v}_i$ and $\vec{p}_i$ and acts on the vertex $\vec{w}_i,\vec{w}_{i+1}$ with relative strength $(1-q)$ and $q$.
$$\begin{align} \vec{F}_{i,\text{External}} = \vec{F}(\vec{v}_i,\vec{p}) \end{align}$$
We provide the RodInteraction
struct to convert a point-wise interaction to a rod-rod interaction potential.
In our case, we use the
MorsePotential building block.
We detect neighbors who are in range of $2R$ where $R$ is the radius of interaction.
Cycle
To simulate proliferation, we introduce a growth term for the spring lengths $l$. We reduce the growth rate with increasing number of neighbors $n$.
$$\begin{equation} \partial_t l = \mu \left(1 - \frac{\min(n,N)}{N}\right) \end{equation}$$
which will increase the length of the cell indefenitely unless we impose a condition for the division event. We define a threshold (in our case double of the original length) for the total length of the cell at which it divides. To construct a new cell, we cannot simply copy the existing one twice, but we also need to adjust internal parameters in the process. The following actions need to be taken for the old and new agent.
- Assign a new growth rate (pick randomly from uniform distribution in $[0.8\text{ }\mu_0,1.2\text{ }\mu_0]$ where $\mu_0$ is some fixed value)
- Assign new positions
- Calculate new spring lengths $$\tilde{l} = l\left(\frac{1}{2} - \frac{r}{\sum\limits_i l}\right)$$
- Calculate middle of old cell $$\vec{m} = \frac{1}{N_\text{vertices}}\sum\limits_i\vec{v}_i$$
- Calculate positions of new vertices $\vec{w}_i$ $$\begin{align} q_i &= \frac{i}{N_\text{vertices}}\\ \vec{w}_{i,\text{new},\pm} &= (1-q_i)\vec{m} + q_i(\vec{v}_{\pm\text{start}} - \vec{m}) \end{align}$$
Domain
In this simulation example, the domain plays an important role. The domain consists of an elongated box with reflective boundary conditions. Cells are initially placed in the left part. Due to their repulsive potential at short distances, they begin to push each other into the remaining space.
Parameters
| Parameter | Symbol | Value |
|---|---|---|
RodMechanics |
||
| Spring Tension | $\gamma$ | $10\text{ min}^{-2}$ |
| Diffusion Constant | $D$ | $0$ $\text{µm}^2\text{min}^{-2}$ |
| Damping | $\lambda$ | $1.5\text{min}^{-1}$ |
| Spring Length | $l$ | $3\text{ µm}$ |
| Rigidity | $\eta$ | $2$ $\text{µm min}^{-1}$ |
MorsePotential |
||
| Interaction Radius | $R$ | $3\text{ µm}$ |
| Potential Stiffness | $\lambda$ | $0.5\text{ µm}$ |
| Cutoff | $5\text{ µm}$ | |
| Strength | $V_0$ | $0.1\text{ µm}^2\text{min}^{-2}$ |
| Cycle | ||
| Division Threshold | $6$ µm | |
| Growth Rate | $\mu$ | $0.1\text{ µm min}^{-1}$ |
CartesianCuboidRods |
||
| Domain Length | $L$ | $50\text{ µm}$ |
| Time | ||
| Stepsize | $\Delta t$ | $0.1\text{ min}$ |
| Save Interval | $\Delta t_\text{save}$ | $2.5\text{ min}$ |
| Total Time | $T$ | $25\text{ h}$ |
Results
Initial Snapshot
Movie
Final Snapshot
Code
The code is part of the examples and can be found in the official github repository under bacterial-rods.
Main Simulation
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
use cellular_raza::prelude::*;
use rand::SeedableRng;
use serde::{Deserialize, Serialize};
pub const METRE: f64 = 1.0;
pub const MILI_METRE: f64 = 1e-3;
pub const MICRO_METRE: f64 = 1e-6;
pub const MINUTE: f64 = 1.0;
pub const HOUR: f64 = 60.0 * MINUTE;
pub const GROWTH_BASE_RATE: f64 = 0.1 * MICRO_METRE / MINUTE;
#[derive(CellAgent, Clone, Deserialize, Serialize)]
pub struct Agent {
#[Mechanics]
mechanics: RodMechanics<f64, 3>,
// Interaction
interaction: RodInteraction<MorsePotential>,
// Cycle
growth_rate: f64,
spring_length_threshold: f64,
}
impl
cellular_raza::concepts::Interaction<
nalgebra::MatrixXx3<f64>,
nalgebra::MatrixXx3<f64>,
nalgebra::MatrixXx3<f64>,
f64,
> for Agent
{
fn calculate_force_between(
&self,
own_pos: &nalgebra::MatrixXx3<f64>,
own_vel: &nalgebra::MatrixXx3<f64>,
ext_pos: &nalgebra::MatrixXx3<f64>,
ext_vel: &nalgebra::MatrixXx3<f64>,
ext_info: &f64,
) -> Result<(nalgebra::MatrixXx3<f64>, nalgebra::MatrixXx3<f64>), CalcError> {
self.interaction
.calculate_force_between(own_pos, own_vel, ext_pos, ext_vel, ext_info)
}
fn get_interaction_information(&self) -> f64 {
self.interaction.0.radius
}
fn is_neighbor(
&self,
own_pos: &nalgebra::MatrixXx3<f64>,
ext_pos: &nalgebra::MatrixXx3<f64>,
_ext_inf: &f64,
) -> Result<bool, CalcError> {
for own_point in own_pos.row_iter() {
for ext_point in ext_pos.row_iter() {
if (own_point - ext_point).norm() < 2.0 * self.interaction.0.radius {
return Ok(true);
}
}
}
Ok(false)
}
fn react_to_neighbors(&mut self, neighbors: usize) -> Result<(), CalcError> {
if neighbors > 0 {
self.growth_rate = (GROWTH_BASE_RATE * (8.0 - neighbors as f64) / 8.0).max(0.0);
} else {
self.growth_rate = GROWTH_BASE_RATE;
}
Ok(())
}
}
impl Cycle<Agent> for Agent {
fn update_cycle(
_rng: &mut rand_chacha::ChaCha8Rng,
dt: &f64,
cell: &mut Agent,
) -> Option<cellular_raza::prelude::CycleEvent> {
cell.mechanics.spring_length += cell.growth_rate * dt;
if cell.mechanics.spring_length > cell.spring_length_threshold {
Some(CycleEvent::Division)
} else {
None
}
}
fn divide(
_rng: &mut rand_chacha::ChaCha8Rng,
cell: &mut Agent,
) -> Result<Agent, cellular_raza::prelude::DivisionError> {
let c2_mechanics = cell.mechanics.divide(cell.interaction.0.radius)?;
let mut c2 = cell.clone();
c2.mechanics = c2_mechanics;
Ok(c2)
}
}
fn main() -> Result<(), SimulationError> {
// Define the dimensionality of the problem
let nrows: usize = 5;
// Define initial random seed
use rand::Rng;
let mut rng = rand_chacha::ChaCha8Rng::seed_from_u64(5);
// Give agent default values
let agent = Agent {
mechanics: RodMechanics {
pos: nalgebra::MatrixXx3::zeros(nrows),
vel: nalgebra::MatrixXx3::zeros(nrows),
diffusion_constant: 0.0 * MICRO_METRE.powf(2.0) / MINUTE,
spring_tension: 10.0 / MINUTE.powf(2.0),
rigidity: 1.0 * MICRO_METRE / MINUTE.powf(2.0),
damping: 1.5 / MINUTE,
spring_length: 3.0 * MICRO_METRE,
},
interaction: RodInteraction(MorsePotential {
radius: 3.0 * MICRO_METRE,
potential_stiffness: 0.5 / MICRO_METRE,
strength: 0.1 * MICRO_METRE.powf(2.0) / MINUTE.powf(2.0),
cutoff: 5.0 * MICRO_METRE,
}),
spring_length_threshold: 6.0 * MICRO_METRE,
growth_rate: GROWTH_BASE_RATE,
};
// Place agents in simulation domain
let domain_size = 50.0 * MICRO_METRE;
let delta_x = agent.mechanics.spring_length * nrows as f64;
let agents = (0..5).map(|_| {
let mut new_agent = agent.clone();
new_agent.mechanics.spring_length = rng.gen_range(1.5..2.5) * MICRO_METRE;
let mut pos = nalgebra::MatrixXx3::zeros(nrows);
pos[(0, 0)] = rng.gen_range(delta_x..2.0 * delta_x);
pos[(0, 1)] = rng.gen_range(delta_x / 3.0..delta_x * 2.0 / 3.0);
pos[(0, 2)] = rng.gen_range(delta_x..2.0 * delta_x);
let theta = rng.gen_range(0.0..2.0 * std::f64::consts::PI);
for i in 1..pos.nrows() {
let phi =
theta + rng.gen_range(-std::f64::consts::FRAC_PI_8..std::f64::consts::FRAC_PI_8);
let mut direction = nalgebra::Vector3::zeros();
direction[0] = phi.cos();
direction[1] = phi.sin();
let new_pos = pos.row(i - 1) + agent.mechanics.spring_length * (direction).transpose();
use core::ops::AddAssign;
pos.row_mut(i).add_assign(new_pos);
}
new_agent.mechanics.set_pos(&pos);
new_agent
});
// Domain Setup
let domain_sizes = [4.0 * domain_size, delta_x, 3.0 * delta_x];
let domain_segments = [8, 1, 2];
let domain = CartesianCuboidRods {
domain: CartesianCuboid::from_boundaries_and_n_voxels(
[0.0; 3],
domain_sizes,
domain_segments,
)?,
};
// Storage Setup
let storage_builder = cellular_raza::prelude::StorageBuilder::new().location("./out");
// Time Setup
let t0 = 0.0 * MINUTE;
let dt = 0.1 * MINUTE;
let save_interval = 2.5 * MINUTE;
let t_max = 25.0 * HOUR;
let time_stepper = cellular_raza::prelude::time::FixedStepsize::from_partial_save_interval(
t0,
dt,
t_max,
save_interval,
)?;
let settings = Settings {
n_threads: 8.try_into().unwrap(),
time: time_stepper,
storage: storage_builder,
show_progressbar: true,
};
println!("Running Simulation");
run_simulation!(
domain: domain,
agents: agents,
settings: settings,
aspects: [Mechanics, Interaction, Cycle],
zero_force_default: |c: &Agent| {
nalgebra::MatrixXx3::zeros(c.mechanics.pos.nrows())
},
)?;
Ok(())
}
Plotting
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
import pyvista as pv
import numpy as np
import argparse
import itertools
import tqdm
import concurrent.futures
from pathlib import Path
import glob
import json
import os
import pandas as pd
def get_last_output_path():
return Path(sorted(glob.glob("out/*"))[-1])
def get_all_iterations(output_path: Path | None = None):
if output_path is None:
output_path = get_last_output_path()
folders = glob.glob(str(output_path / "cells/json") + "/*")
iterations = [int(Path(fi).name) for fi in folders]
return output_path, iterations
def load_cells_from_iteration(output_path: Path, iteration: int):
load_path = Path(output_path) / "cells/json/{:020}".format(iteration)
results = []
for file in glob.glob(str(load_path) + "/*"):
f = open(file)
elements = json.load(f)["data"]
elements = [element["element"][0] for element in elements]
results.extend(elements)
df = pd.json_normalize(results)
if len(df) > 0:
# Format individual entries for easier use later on
df["identifier"] = df["identifier"].apply(lambda x: tuple(x))
df["cell.mechanics.pos"] = df["cell.mechanics.pos"].apply(
lambda x: np.array(x[0], dtype=float).reshape(3, -1)
)
df["cell.mechanics.vel"] = df["cell.mechanics.vel"].apply(
lambda x: np.array(x[0], dtype=float)
)
return df
def get_cell_meshes(iteration: int, path: Path):
cells = load_cells_from_iteration(path, iteration)
positions = np.array([x for x in cells["cell.mechanics.pos"]], dtype=float)
radii = np.array([x for x in cells["cell.interaction.radius"]], dtype=float)
growth_rate = np.array([x for x in cells["cell.growth_rate"]], dtype=float)
cell_surfaces = []
for i, p in enumerate(positions):
meshes = []
# Add the sphere at the first point
meshes.append(pv.Sphere(center=p[:,0], radius=radii[i]))
for j in range(max(p.shape[1]-1,0)):
# Add a sphere for each following point
meshes.append(pv.Sphere(center=p[:,j+1], radius=radii[i]))
# Otherwise add cylinders
pos1 = p[:,j]
pos2 = p[:,j+1]
center = 0.5 * (pos1 + pos2)
direction = pos2 - center
radius = radii[i]
height = float(np.linalg.norm(pos1 - pos2))
cylinder = pv.Cylinder(center, direction, radius, height)
meshes.append(cylinder)
combined = pv.MultiBlock(meshes).combine()
if combined is None:
return None
merged = combined.extract_surface()
if merged is None:
return None
merged = merged.clean()
cell_surfaces.append((merged, growth_rate[i]))
return cell_surfaces
def plot_spheres(
iteration: int,
path: Path = Path("./"),
opath: Path | None = None,
overwrite: bool = False,
transparent_background: bool = False,
):
if opath is None:
opath = path / "images/{:010}.png".format(iteration)
opath.parent.mkdir(parents=True, exist_ok=True)
if os.path.isfile(opath) and overwrite is False:
return None
cell_meshes = get_cell_meshes(iteration, path)
cell_meshes = cell_meshes if cell_meshes is not None else []
# General Settings
plotter = pv.Plotter(off_screen=True, window_size=[1280, 1280])
plotter.set_background([100, 100, 100])
# Draw box around everything
dx = 3e-6
box = pv.Box(bounds=(-dx, 200e-6+dx, 0-dx, 15e-6+dx, 0-dx, 45e-6+dx))
plotter.add_mesh(box, style="wireframe", line_width=15)
color_min = np.array([69, 124, 214])
color_max = np.array([82, 191, 106])
for (cell, growth_rate) in cell_meshes:
# TODO MAGIC NUMBERS
q = growth_rate / 1e-7
color = (1-q) * color_min + q * color_max
plotter.add_mesh(
cell,
show_edges=False,
color=(int(color[0]), int(color[1]), int(color[2])),
diffuse=0.5,
ambient=0.5,
)
# Define camera
# TODO MAGIC NUMBERS
plotter.camera.position = (100e-6, -300e-6, -300e-6)
plotter.camera.focal_point = (100e-6, 25e-6, 22.5e-6)
plotter.enable_ssao(radius=12)
plotter.enable_anti_aliasing()
img = plotter.screenshot(opath, transparent_background=transparent_background)
plotter.close()
return img
def __plot_spheres_helper(args):
(args, kwargs) = args
plot_spheres(*args, **kwargs)
def plot_all_spheres(
path: Path,
n_threads: int | None = None,
overwrite:bool=False,
transparent_background: bool = False
):
iterations = [it for it in get_all_iterations(path)[1]]
if n_threads is None:
n_cpu = os.cpu_count()
if n_cpu is not None and n_cpu > 2:
n_threads = n_cpu -2
else:
n_threads = 1
args = [(it,) for it in iterations]
kwargs = {
"path": path,
"overwrite": overwrite,
"transparent_background": transparent_background,
}
with concurrent.futures.ProcessPoolExecutor(max_workers=n_threads) as executor:
_ = list(tqdm.tqdm(executor.map(
__plot_spheres_helper,
zip(args, itertools.repeat(kwargs))
), total=len(iterations)))
# TODO actually use this
if __name__ == "_main__":
parser = argparse.ArgumentParser(
prog="Plotter",
description="A plotting CLI for the springs example",
epilog="For suggestions or bug-reports please refer to\
https://github.com/jonaspleyer/cellular_raza/"
)
parser.add_argument(
"-n",
"--iteration",
help="Specify which iteration to plot.\
Multiple arguments are accepted as a list.\
If not specified, plot all iterations."
)
parser.add_argument(
"-i",
"--input",
help="Input path of files. If not given,\
it will be determined automatically by searching in './out/'"
)
parser.add_argument(
"-o",
"--output-path",
help="Path where to store generated images."
)
args = parser.parse_args()
if __name__ == "__main__":
print("Plotting Individual Snapshots")
output_path = get_last_output_path()
plot_all_spheres(output_path, transparent_background=True, overwrite=True)
print("Generating Movie")
bashcmd = f"ffmpeg\
-v quiet\
-stats\
-y\
-r 30\
-f image2\
-pattern_type glob\
-i '{output_path}/images/*.png'\
-c:v h264\
-pix_fmt yuv420p\
-strict -2 {output_path}/movie.mp4"
os.system(bashcmd)
print("Playing Movie")
bashcmd2 = f"firefox {output_path}/movie.mp4"
os.system(bashcmd2)
References
[1] A. Zapun, T. Vernet, and M. G. Pinho, “The different shapes of cocci,” FEMS Microbiology Reviews, vol. 32, no. 2. Oxford University Press (OUP), pp. 345–360, Mar. 2008. doi: 10.1111/j.1574-6976.2007.00098.x
[2] K. D. Young, “The Selective Value of Bacterial Shape,” Microbiology and Molecular Biology Reviews, vol. 70, no. 3. American Society for Microbiology, pp. 660–703, Sep. 2006. doi: 10.1128/mmbr.00001-06.
[3] C. Billaudeau et al., “Contrasting mechanisms of growth in two model rod-shaped bacteria,” Nature Communications, vol. 8, no. 1. Springer Science and Business Media LLC, Jun. 07, 2017. doi: 10.1038/ncomms15370.