A Mechanical Model of Rod-Shaped Bacteria

## A Mechanical Model of
                        ## Rod-Shaped Bacteria
                        #### Jonas Pleyer
                        ##### 03.02.2025

---

## Overview

| | |
                        |:---|---|
                        | 1. Rod-Shaped Bacteria | <div class="n-slides"> (3 slides)</div> |
                        | 2. Mathematical Model | <div class="n-slides"> (9 slides)</div>|
                        | 3. Computational Algorithms | <div class="n-slides"> (5 slides)</div>|
                        | 4. Parameter Estimation | <div class="n-slides"> (3 slides)</div>|
                        | 5. Outlook | <div class="n-slides"> (2 slide)</div>|

## Rod-Shaped Bacteria

- One of many bacterial
                          shapes
                            - We do not consider coccobacillus (elliptical shape)
                              for now
                            - Main interest: bacillus
                        - Length $\approx 1-10µm$, Diameter$\approx0.5-2µm$
                          
                        - Elongated shape (not simply point-particle) but flexible
                          
                          [(Amir et al., 2014)](https://doi.org/10.1073/pnas.1317497111)
                        - Asymmetric growth in direction of long axis
                          
                        - Why are we interested in this?<br>
                          <div style="color: orange;">$\Rightarrow$ We will get to that later</div>

</div><div class="fragment" data-fragment-index="1">
                        <img src="./Bacterial_morphology_diagram.svg" />
                        <div style="font-size: 0.6em">wikimedia.org/wiki/File:Bacterial_morphology_diagram.svg</div>
                        </div>
                        </div>

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` | <div style="width: 100%; text-align: center;">⎯⎯⎯ \" ⎯⎯⎯<div> |
                        | [(Martins et al., 2015)](https://doi.org/10.1109/ENBENG.2015.7088854) | Constructed model but no parameter estimations. |
                        | [(Winkle et al., 2017)](https://doi.org/10.1088/1478-3975/aa7bae) | Authors used game engine `Chipmunk3D` to simulate<br>rods but only studied global phenomena. |
                        | [(Doumic et al, 2020)](https://doi.org/10.3934/mbe.2020356) | Estimated distribution of some parameters (growth rate,<br>length) but no bending and not on individual level. |

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

</div>

<h2>Mechanics</h2>
                        <div style="display: grid; grid-template-columns: 45% 55%;">
                            <img src="../mechanics.png">
                            <div class="fragment" data-fragment-index="1">

- 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

1. 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}$$

2. 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}$$

</div>

## Growth

- Elongation is constant in length growth<br>
                          
                          [(Grover et al., 1977](https://doi.org/10.1016/0022-5193(77)90192-8)
                          [Rosenberger et al., 1978)](https://doi.org/10.1016/0022-5193(78)90132-7)
                        - 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}$$

</div><div>

## Cell Cycle

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

</div>

---

## 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}$$

```rust
                        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;
                        }
                        ```

</div>
                        </div>

---

## Interaction
                        #### Morse-Potential

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

</div>
                        <div class="fragment" data-fragment-index="2">

| | 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 |

</div>
                        </div>

---

## 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}

</div>
                        <div clas="fragment" data-fragment-index="3">

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

</div>
                        <div class="fragment" data-fragment-index="1">

</div>
                        <div class="fragment" data-fragment-index="2">

| | 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 |

</div>
                        </div>

</div>

---

## Comparison of Interaction Potentials

### Morse
                        ### Mie

<div>

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

</div>
                        <div>

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

</div>
                        </div>

---

<div style="display: grid; grid-template-columns: 45% 55%; width: 100%;">
                        <img
                            src="https://cellular-raza.com/logos/cellular_raza_dark_mode.svg"
                            style="margin: auto;"
                        />
                        <h2 style="margin: auto;">Variables/Parameters</h2>
                        <div>

```rust [1-2|3-5|6-8|9-13|1-13]
                        #[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,
                        }
                        ```
                        </div>
                        <div class="fragment">

```rust [1-2|3-6|7-16|1-17]
                        #[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,
                        }
                        ```

</div>
                        </div>
                        <div>
                        <div class="fragment">

| | 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 |

</div>
                        </div>

---

## Numerical Simulation Example

---

## Image Analysis

<div style="display: grid; grid-template-columns: 1fr 1fr;column-gap: 2em;">
                            <h4>Time Series of Microscopic Images</h4>
                            <h4></h4>
                            <div style="display: grid; grid-template-columns: 1fr 1fr 1fr">
                                <img width=100% src="growth-2/image001042.png"/>
                                <img width=100% src="growth-2/image001062.png"/>
                                <img width=100% src="growth-2/image001082.png"/>
                                <img width=100% src="growth-2/image001102.png"/>
                                <img width=100% src="growth-2/image001122.png"/>
                                <img width=100% src="growth-2/image001142.png"/>
                                <img
                                    width=100%
                                    src="growth-2-marked/image001042-marked.png"
                                    class="fragment"
                                    data-fragment-index="2"
                                />
                                <img
                                    width=100%
                                    src="growth-2-marked/image001062-marked.png"
                                    class="fragment"
                                    data-fragment-index="2"
                                />
                                <img
                                    width=100%
                                    src="growth-2-marked/image001082-marked.png"
                                    class="fragment"
                                    data-fragment-index="2"
                                />
                            </div>
                        <div>

<p style="color: orange;" class="fragment" data-fragment-index="1">
                            We need methods to map between images<br>
                            and our numerical simulation
                        </p>

1. Extract individual cells from image<br>
                           
                           ✅ (_Cell-Segmentation_)
                        2. Extract positions from individual cell<br>
                           
                           ❌ (_No existing algorithm_)
                        3. Advance simulation<br>
                           
                           ✅ (_see implementation of our model_)
                        4. Compare numerical and real-world results<br>
                           
                           ❓ (_what is the cost function?_)

</div>
                        </div>

---

## Image Analysis
                        #### Extracting Positions

<div style="display: grid; grid-template-columns: 1fr 1fr 1fr;">
                            <p class="fragment" data-fragment-index="1">Microscopic Image</p>
                            <p class="fragment" data-fragment-index="2">Cell Masks</p>
                            <p class="fragment" data-fragment-index="3">Skeletonization</p>
                            <img
                                width=100%
                                src="../fitting-methods/algorithm/image001042.png"
                                class="fragment"
                                data-fragment-index="1"
                            >
                            <img
                                width=100%
                                src="growth-2-marked/image001042-marked.png"
                                class="fragment"
                                data-fragment-index="2"
                            >
                            <img
                                width=100%
                                src="../fitting-methods/algorithm/mask-zoom.png"
                                class="fragment"
                                data-fragment-index="3"
                            >
                        </div>

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

---

## Image Analysis
                        #### Generating Images

- We can use the numerical simulation to generate<br>
                          
                          artificial data & images
                        - Create masks by projecting 3D image &<br>
                          
                          assigning unique color to each cell
                        - Use this data to test our position extraction algorithm<br>
                          
                          Advantage: known input positions

## Image Analysis
                        #### How to generate Near-Realistic Microscopic Images

- Noise
                          
                            - Camera Sensor Noise
                            - Statistic Noise from environmental influences
                        - Defects
                          
                            - Lens Aberration
                            - Dirt on Glass
                            - etc.
                        - Lighting (backlight most of the time)
                          
                        - Cellular properties: reflectivity, luminosity, tranlucency
                          
                        - Overlapping of bacteria
                          
                        - Dirt/Other material floating around

<p class="fragment" style="color: orange;">
                        $\Rightarrow$ Physically-Based Rendering!
                        </p>

</div><img width=100% src="../09395645494836445480/images/000000800.png">
                        </div>

---

## Image Analysis
                        #### Extracting Positions

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

</div>

---

## Image Analysis
                        #### Extracting Positions

![](../fitting-methods/displacement-calculations.png)

- 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

</div>
                        </div>

---

## Parameter Estimation
                        #### Constructing a Cost Function

- Assumption: We have no cell-division
                        - $\Rightarrow$ compare euclidean distance between vertices
                        - If we have cell-division, we need to do something else

## Parameter Estimation
                        #### Comparing Across Generations

</div></div>

---

## 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

</div>
                            <img
                                width=100%
                                src="../fitting-methods/estimate-parameters1/microscopic-images-2.png"
                            >
                            <img
                                width=100%
                                src="../fitting-methods/estimate-parameters1/microscopic-images-3.png"
                            >
                        </div>

---

## 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)

</div>

---

## Parameter Estimation
                        #### Optimization Strategies

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

---

### Python Package

```bash
                        pip install cr_mech_coli
                        ```

```python [1|3-5|7-11|13-19|1-19]
                        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
                        )
                        ```

</div>
                        <div>
                        <h3>Documentation</h3>
                        <iframe
                            src="https://jonaspleyer.github.io/cr_mech_coli/"
                            style="width: 1100px; height: 980px; max-width: 100%; max-height: 100%;"
                        ></iframe></div>
                        </div>

---

## Outlook

1. 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
                        2. Distribution of Parameters
                          
                            - Can we fit enough cells individually simultanously?
                            - How are parameters of cells distributed?
                        3. Inheritance of Parameters
                          
                            - Will our algorithms work over multiple bacteria generations?
                            - How are parameters inherited to daughter cells?

<p style="color: orange;" class="fragment">
                        We are meeting up with Alexander Rohrbach<br>
                        this week to discuss how to take things further.
                        </p>