A Mechanical Model of Rod-Shaped Bacteria

Overview


1. Rod-Shaped Bacteria	(3 slides)
2. Mathematical Model	(9 slides)
3. Computational Algorithms	(5 slides)
4. Parameter Estimation	(3 slides)
5. Outlook	(2 slide)

Rod-Shaped Bacteria

Growth of E.Coli

Institut für den Wissenschaftlichen Film - Göttingen 1981

Rod-Shaped Bacteria

State of current models

ABM/Paper	State
Biocellion	Only cylindrical potentials
BSim	Supports Rod-Shaped Bacteria; No parameters estimated.
`gro`	⎯⎯⎯ " ⎯⎯⎯
(Martins et al., 2015)	Constructed model but no parameter estimations.
(Winkle et al., 2017)	Authors used game engine `Chipmunk3D` to simulate rods but only studied global phenomena.
(Doumic et al, 2020)	Estimated distribution of some parameters (growth rate, length) but no bending and not on individual level.

$\Rightarrow$ No literature on measuring properties of individual cells!

Mechanics

Vertices $\vec{x}_i$ connected by springs
Exert force between two adjacent vertices

$\vec{c}_i = \vec{x}_{i}-\vec{x}_{i-1}$ $\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}$

Mechanics

Bending force acts on vertices $\vec{x}_{i-1},\vec{x}_i,\vec{x}_{i+1}$
Proportional to curvature $\kappa_i$ at vertex $\vec{x}_i$

$\begin{align} \kappa_i &= 2\tan\left(\frac{\alpha_i}{2}\right)\\ \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}$

Mechanics

Sum all 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}$

Integrate equations of motion

$\begin{align} \partial_t^2 \vec{x} &= \partial_t\vec{x} + \sqrt{2D}\vec{\xi}\\ \partial_t\vec{x} &= \vec{F}_\text{total} - \lambda \partial_t\vec{x} \end{align}$

Growth

Elongation is constant in length growth
(Grover et al., 1977 Rosenberger et al., 1978)
Model as linear growth

$\begin{equation} \partial_t l = \mu \end{equation}$

For larger collections: introduce phenomenological saturation term depending on number of neighbors $n$

$\begin{equation} \partial_t l = \mu\left(1 - \frac{\text{min}(n, N)}{N}\right) \end{equation}$

Cell Cycle

Cell divides in the middle when reaching threshold length $l_\text{max}$
Calculate new spring lengths
Calculate new positions of vertices for both cells (interpolate new vertices between vertices of mother-cell)
(Optional) Pick new growth rate (or other parameters)

Interaction

Calculate interaction for each vertex $\vec{x}_i$
Pointwise to nearest point on ther cells polygon

$\vec{p} = h\vec{w}_j + (1-h)\vec{w}_{j+1}$

for (i, vi) in vertices_cell_1.iter().enumerate() {
    // Calculate nearest point, index and segment as above
    let (p, j, h) = get_nearest_point_index_segment(&vi, &cell2)?;

    // Calculate force between these two points
    let calculated_force = cell1.calculate_force_between(&p)?;

    // Store calculated force for cell1 and cell2
    forces_cell_1[i] +=         calculated_force;
    forces_cell_2[j] -=    h  * calculated_force;
    forces_cell_2[j] -= (1-h) * calculated_force;
}

Interaction

Morse-Potential

$\begin{equation} V(r) = V_0\left(1 - e^{-\lambda(r-R)}\right)^2 \end{equation}$

	Name	Type
$R$	Radius	`f32` $\in\mathbb{R}^+$	Fit
$\lambda$	Stiffness	`f32` $\in\mathbb{R}^+$	Fit
$\zeta$	Cutoff	`f32` $\in\mathbb{R}^+$	Assumption
$V_0$	Strength	`f32` $\in\mathbb{R}^+$	Fit

Interaction

Mie-Potential

$\begin{equation} U(r) = C\epsilon\left[ \left(\frac{\sigma}{r}\right)^n - \left(\frac{\sigma}{r}\right)^m\right] \end{equation}$

$C = \frac{n}{n-m}\left(\frac{n}{m}\right)^{\frac{n}{n-m}}$

	Name	Type
$R$	Radius	`f32` $\in\mathbb{R}^+$	Fit
$\epsilon$	Strength	`f32` $\in\mathbb{R}^+$	Fit
$\beta$	Bound	`f32` $\in\mathbb{R}^+$	Assumption
$\zeta$	Cutoff	`f32` $\in\mathbb{R}^+$	Assumption
$n,m$	Exponents	`f32` $\in\mathbb{R}^+$	Fit

Comparison of Interaction Potentials

Morse

Mie

Less steep for $r\rightarrow0$
Width of the well affects range and steepness
Less parameters

Can represent more types of interactions
Needs an upper bound for $r\rightarrow0$

Variables/Parameters

 #[derive(CellAgent, Clone, Debug, Deserialize, Serialize)]
pub struct RodAgent {
    /// Determines mechanical properties of the agent
    #[Mechanics]
    pub mechanics: RodMechanics<f32, 3>,
    /// Determines interaction between agents
    #[Interaction]
    pub interaction: RodInteraction<PhysicalInteraction>,
    /// Rate with which the cell grows
    pub growth_rate: f32,
    /// Threshold at which the cell will divide
    pub division_length: f32,
} #[derive(CellAgent, Clone, Debug, Deserialize, Serialize)]
pub struct RodAgent {
    /// Determines mechanical properties of the agent
    #[Mechanics]
    pub mechanics: RodMechanics<f32, 3>,
    /// Determines interaction between agents
    #[Interaction]
    pub interaction: RodInteraction<PhysicalInteraction>,
    /// Rate with which the cell grows
    pub growth_rate: f32,
    /// Threshold at which the cell will divide
    pub division_length: f32,
} #[derive(CellAgent, Clone, Debug, Deserialize, Serialize)]
pub struct RodAgent {
    /// Determines mechanical properties of the agent
    #[Mechanics]
    pub mechanics: RodMechanics<f32, 3>,
    /// Determines interaction between agents
    #[Interaction]
    pub interaction: RodInteraction<PhysicalInteraction>,
    /// Rate with which the cell grows
    pub growth_rate: f32,
    /// Threshold at which the cell will divide
    pub division_length: f32,
} #[derive(CellAgent, Clone, Debug, Deserialize, Serialize)]
pub struct RodAgent {
    /// Determines mechanical properties of the agent
    #[Mechanics]
    pub mechanics: RodMechanics<f32, 3>,
    /// Determines interaction between agents
    #[Interaction]
    pub interaction: RodInteraction<PhysicalInteraction>,
    /// Rate with which the cell grows
    pub growth_rate: f32,
    /// Threshold at which the cell will divide
    pub division_length: f32,
} #[derive(CellAgent, Clone, Debug, Deserialize, Serialize)]
pub struct RodAgent {
    /// Determines mechanical properties of the agent
    #[Mechanics]
    pub mechanics: RodMechanics<f32, 3>,
    /// Determines interaction between agents
    #[Interaction]
    pub interaction: RodInteraction<PhysicalInteraction>,
    /// Rate with which the cell grows
    pub growth_rate: f32,
    /// Threshold at which the cell will divide
    pub division_length: f32,
}

 #[derive(Clone, Debug, PartialEq)]
pub struct RodMechanics<f32, 3> {
    /// The current position
    pub pos: Matrix<f32, Dyn, Const<3>, _>,
    /// The current velocity
    pub vel: Matrix<f32, Dyn, Const<3>, _>,
    /// Controls stochastic motion
    pub diffusion_constant: f32,
    /// Spring tension between vertices
    pub spring_tension: f32,
    /// Stiffness between two edges
    pub rigidity: f32,
    /// Target spring length
    pub spring_length: f32,
    /// Daming constant
    pub damping: f32,
} #[derive(Clone, Debug, PartialEq)]
pub struct RodMechanics<f32, 3> {
    /// The current position
    pub pos: Matrix<f32, Dyn, Const<3>, _>,
    /// The current velocity
    pub vel: Matrix<f32, Dyn, Const<3>, _>,
    /// Controls stochastic motion
    pub diffusion_constant: f32,
    /// Spring tension between vertices
    pub spring_tension: f32,
    /// Stiffness between two edges
    pub rigidity: f32,
    /// Target spring length
    pub spring_length: f32,
    /// Daming constant
    pub damping: f32,
} #[derive(Clone, Debug, PartialEq)]
pub struct RodMechanics<f32, 3> {
    /// The current position
    pub pos: Matrix<f32, Dyn, Const<3>, _>,
    /// The current velocity
    pub vel: Matrix<f32, Dyn, Const<3>, _>,
    /// Controls stochastic motion
    pub diffusion_constant: f32,
    /// Spring tension between vertices
    pub spring_tension: f32,
    /// Stiffness between two edges
    pub rigidity: f32,
    /// Target spring length
    pub spring_length: f32,
    /// Daming constant
    pub damping: f32,
} #[derive(Clone, Debug, PartialEq)]
pub struct RodMechanics<f32, 3> {
    /// The current position
    pub pos: Matrix<f32, Dyn, Const<3>, _>,
    /// The current velocity
    pub vel: Matrix<f32, Dyn, Const<3>, _>,
    /// Controls stochastic motion
    pub diffusion_constant: f32,
    /// Spring tension between vertices
    pub spring_tension: f32,
    /// Stiffness between two edges
    pub rigidity: f32,
    /// Target spring length
    pub spring_length: f32,
    /// Daming constant
    pub damping: f32,
}

	Property	Type
$\mu$	Growth Rate	`f32` $\in\mathbb{R}^+$	Fit
$l_\text{max}$	Division Length	`f32` $\in\mathbb{R}^+$	Fit
$\vec{x}_i$	Position	`Matrix` $\in\mathbb{R}^{d\times 3}$	I.V.
$\vec{v}_i$	Velocity	`Matrix` $\in\mathbb{R}^{d\times 3}$	I.V.
-	Diffusion	`f32` $\in\mathbb{R}^+$	$0$
$l$	Spr. Length	`f32` $\in\mathbb{R}^+$	I.V.
$\gamma$	Spr. Tension	`f32` $\in\mathbb{R}^+$	Fit
$\eta$	Rigidity	`f32` $\in\mathbb{R}^+$	Fit
$\lambda$	Damping	`f32` $\in\mathbb{R}^+$	Fit

Numerical Simulation Example

Image Analysis

Time Series of Microscopic Images

We need methods to map between images
and our numerical simulation

Extract individual cells from image
✅ (Cell-Segmentation)
Extract positions from individual cell
❌ (No existing algorithm)
Advance simulation
✅ (see implementation of our model)
Compare numerical and real-world results
❓ (what is the cost function?)

Image Analysis

Extracting Positions

Microscopic Image

Cell Masks

Skeletonization

Take skeleton of cell-mask and interpolate to polygon with $N$ vertices

Image Analysis

Extracting Positions

Crosses show extracted positions
Lines show actual polygon
Seems to match quite well

Image Analysis

Extracting Positions

Measure difference between rod lengths
Calculate average vertex difference from extracted to known position
Good agreement with given positions
In the beginning small underestimation
In the end small overestimation

Parameter Estimation

Results

2 Pairs are matching well
Left pair is not aligned very well

$\Rightarrow$ this could be due to the MorsePotentialF32 interaction potential

Parameter Estimation

Results

Only one parameter for all agents (not individual)
Only one time-step is fitted
Duration: $\approx 15\text{min}$ (depending on your hardware)

Parameter Estimation

Optimization Strategies

Estimate parameters for each time-step
Run simulation over multiple time-steps and compare results
- option to rebase every $n$ time-steps and start with new initial conditions
Global vs Local optimization
- unfortunately in general not possible to calculate Jacobian

Python Package

pip install cr_mech_coli

 import cr_mech_coli as crm
 
# Configure inputs
config = crm.Configuration()
agent_settings = crm.AgentSettings()
 
# Run simulation
cell_container = crm.run_simulation(
    config,
    agent_settings,
)
 
# Render Images
render_settings = crm.RenderSettings()
crm.store_all_images(
    cell_container,
    config.domain_size,
    render_settings
) import cr_mech_coli as crm
 
# Configure inputs
config = crm.Configuration()
agent_settings = crm.AgentSettings()
 
# Run simulation
cell_container = crm.run_simulation(
    config,
    agent_settings,
)
 
# Render Images
render_settings = crm.RenderSettings()
crm.store_all_images(
    cell_container,
    config.domain_size,
    render_settings
) import cr_mech_coli as crm
 
# Configure inputs
config = crm.Configuration()
agent_settings = crm.AgentSettings()
 
# Run simulation
cell_container = crm.run_simulation(
    config,
    agent_settings,
)
 
# Render Images
render_settings = crm.RenderSettings()
crm.store_all_images(
    cell_container,
    config.domain_size,
    render_settings
) import cr_mech_coli as crm
 
# Configure inputs
config = crm.Configuration()
agent_settings = crm.AgentSettings()
 
# Run simulation
cell_container = crm.run_simulation(
    config,
    agent_settings,
)
 
# Render Images
render_settings = crm.RenderSettings()
crm.store_all_images(
    cell_container,
    config.domain_size,
    render_settings
) import cr_mech_coli as crm
 
# Configure inputs
config = crm.Configuration()
agent_settings = crm.AgentSettings()
 
# Run simulation
cell_container = crm.run_simulation(
    config,
    agent_settings,
)
 
# Render Images
render_settings = crm.RenderSettings()
crm.store_all_images(
    cell_container,
    config.domain_size,
    render_settings
)

Documentation

Outlook

Cell-Segmentation & Cell-Tracking
- Generate Dataset of near-realistic images with masks
- Asses if this can improve existing Cell-Segmentation tools
- This allows us to easily study Cell-Tracking solutions
Distribution of Parameters
- Can we fit enough cells individually simultanously?
- How are parameters of cells distributed?
Inheritance of Parameters
- Will our algorithms work over multiple bacteria generations?
- How are parameters inherited to daughter cells?

We are meeting up with Alexander Rohrbach
this week to discuss how to take things further.

	#[derive(CellAgent, Clone, Debug, Deserialize, Serialize)]
	pub struct RodAgent {
	/// Determines mechanical properties of the agent
	#[Mechanics]
	pub mechanics: RodMechanics<f32, 3>,
	/// Determines interaction between agents
	#[Interaction]
	pub interaction: RodInteraction<PhysicalInteraction>,
	/// Rate with which the cell grows
	pub growth_rate: f32,
	/// Threshold at which the cell will divide
	pub division_length: f32,
	}

	#[derive(Clone, Debug, PartialEq)]
	pub struct RodMechanics<f32, 3> {
	/// The current position
	pub pos: Matrix<f32, Dyn, Const<3>, _>,
	/// The current velocity
	pub vel: Matrix<f32, Dyn, Const<3>, _>,
	/// Controls stochastic motion
	pub diffusion_constant: f32,
	/// Spring tension between vertices
	pub spring_tension: f32,
	/// Stiffness between two edges
	pub rigidity: f32,
	/// Target spring length
	pub spring_length: f32,
	/// Daming constant
	pub damping: f32,
	}

	import cr_mech_coli as crm

	# Configure inputs
	config = crm.Configuration()
	agent_settings = crm.AgentSettings()

	# Run simulation
	cell_container = crm.run_simulation(
	config,
	agent_settings,
	)

	# Render Images
	render_settings = crm.RenderSettings()
	crm.store_all_images(
	cell_container,
	config.domain_size,
	render_settings
	)

A Mechanical Model of

Rod-Shaped Bacteria

Jonas Pleyer

03.02.2025