cellular_raza - Think like a cell!

Jonas Pleyer

- Think like a cell!

05.02.2024

Think like a cell!

05.02.2024

Jonas Pleyer

Outline

Agent-Based Modelling

What? Why? How?
Example Workflow

From theory to practice
Implementation details

## Cellular Properties
            | Property | Examples |
            |---|---|
            | Mechanics | Stochastic Motion, Chemotaxis |
            | Forces | Collision, Friction, Adhesion |
            | Intracellular Dynamics | Reaction Networks, Metabolics |
            | Cycle | Division, Death, Differentiation |
            <div style="margin-left: auto; margin-right: auto; font-size: 0.75em;">
              <img src="media/E._coli_Bacteria_(16578744517).jpg" width="600px">
              <p>Credit: 
                <a href="https://upload.wikimedia.org/wikipedia/commons/d/d5/E._coli_Bacteria_%2816578744517%29.jpg">en.m.wikipedia.org/</a>
                (under CC 2.0 License)
              </p>
            </div>

---

## Domain Properties
            | Property | Examples |
            |---|---|
            | Physical Shape | Petri Dish, Mother-Machine |
            | Fluid-Dynamics | Diffusion, Flow |
            | Custom | Removal/Insertion of cells |
            <div style="margin-left: auto; margin-right: auto; font-size: 0.75em;">
              <img src="media/202004_petri_dish_with_cultured_bacteria_blue.svg" width="600px">
              <p>Credit: 
                <a href="https://upload.wikimedia.org/wikipedia/commons/7/73/202004_petri_dish_with_cultured_bacteria_blue.svg">en.m.wikipedia.org/</a>
                (under CC 4.0 License)
              </p>
            </div>

---

## Example - Bacterial Growth Model
            Include delay due to metabolism in typical growth model

\\begin{alignat}{7}
              \dot{L} &= -&&\lambda R L\\\\
              \dot{G} &= &&\lambda R L &+ \alpha R &G\\\\
              \dot{R} &= &&&-\frac{\alpha}{G_\text{max}} R &G
            \\end{alignat}

| Variable | Description |
            | --- | --- |
            | $L$ | Cells in lag-phase |
            | $G$ | Cells in growth-phase
            | $R$ | Resource (nutrients) |
            | $\lambda$ | Transition rate from $L$ to $G$ |
            | $\alpha$ | Growth rate |

---

## Two-Species Model
            ### with inhibition

\\begin{alignat}{7}
              \dot{L_A} &= -&&\lambda_A R L_A\\\\
              \dot{G_A} &= &&\lambda_A R L_A &+ \alpha_A R &G_A\\\\
              \dot{L_B} &= -&&\lambda_B R L_B\\\\
              \dot{G_B} &= &&\lambda_B R L_B &&&+ \frac{\alpha_B}{1 + \mu_I I} &R G_B\\\\
              \dot{R} &= &&&-\frac{\alpha_A}{G_\text{max}} R &G_A &-\frac{\alpha_B}{(1 + \mu_I I) G_\text{max}} &R G_B \\\\
              \dot{I}   &= &&\kappa G_A
            \\end{alignat}

- Inhibitor $I$ is produced with $\kappa$ by species $A$

---

## Observations
            - Variables $L,G,R,I$ are global
            - Parameters $\lambda,\alpha,\kappa,\mu$ correspond to global processes
            - $G_{max}$ is heuristic

---

## Think like a cell
            - Identify processes
            - Link to cellular properties
            - Find mathematical description
            - Choose physical domain

---

## Processes

| Process | Cellular Property | Description | 
            | --- | --- | --- |
            | $L \longrightarrow G$ | Cycle | Stochastic event |
            | Proliferation | Cycle | Division rule |
            | Shape | Forces + Mechanics | Spherical Potential |
            | Nutrient consumption | Intracellular Dynamics | coupled ODE |

<br>

### Physical Domain
            Cartesian Cuboid with extracellular diffusion

---

### Reactions
            \\begin{align}
              \dot{V_c} &= \frac{u}{1+\nu_I [I]} \sigma [R]V_c\\\\
              \dot{[R]} &= \sum\limits_c -\frac{u}{1+\nu_I [I]} \frac{V_c}{V_D} [R] \delta(x-x_c) + D_R\Delta [R]\\\\
              \dot{[I]} &= \sum\limits_c b\frac{V_c}{V_D}\delta(x-x_c) + D_I\Delta [I]
            \\end{align}

### Forces + Mechanics
            \\begin{equation}
              m \ddot{\vec{x}} = \vec{F} - \lambda \dot{\vec{x}}
            \\end{equation}

### Cycle
            - Custom update rules for agents

---

<h2>Snapshots</h2>
            <div style="display: grid; grid-template-columns: 40% 40%; margin-left: 20%;">
              <img src="media/pool_model-snapshot_00006000.png" style="width: 100%; margin: 0;">
              <img src="media/pool_model-snapshot_00012000.png" style="width: 100%; margin: 0;">
              <img src="media/pool_model-snapshot_00018000.png" style="width: 100%; margin: 0;">
              <img src="media/pool_model-snapshot_00024000.png" style="width: 100%; margin: 0;">
            </div>

---

## Comments
            - only mechanistic parameters
            - coarse-grained case is readily obtained
              - same amount of parameters
            - even just the mindset is helpful
              - no numerical simulation required

---

- Transfer approach to computational model
            - Require Flexibility in model design
            - User should be in control
            - Aim for good performance

## Etymology

<div>

`cellular_raza` = Cell + "tabula rasa"

</div>

- Start from "clean slate"

---

## Implementation

- Use generic programming
              - Allows for flexibility
              - Retains performance

```rust
            // Rust example
            fn update<T>(element: &mut T) {...}
            //                         ↑
            //         This is a generic argument
            ```

```cpp
            // C++ example
            template<class T>
            void update(T &element) {...}
            //          ↑
            //         This is a template argument
            ```

</div>

- Build all update functions with generics

---

## Inherent Assumptions
            - All interactions are local (finite range)
              - Cells only interact with next neighbours
              - and environment
              - No global rules
            - No explicit time-dependence in system

---

## Mechanics
            ```rust
            pub trait Mechanics<Pos, Vel, For, Float = f64> {
              fn pos(&self) -> Pos;
              fn velocity(&self) -> Vel;
              fn set_pos(&mut self, pos: &Pos);
              fn set_velocity(&mut self, velocity: &Vel);
              fn calculate_increment(&self, force: For) -> Result<(Pos, Vel), CalcError>;
            }
            ```

| Name | Dim | | | Equations |
            | --- |:---:| ---:|:--- | --- |
            | `Brownian` | 1, 2, 3 | $\dot{\vec{x}}$ | $=$ | $-\frac{D}{k_B T}\nabla V(x) + \sqrt{2D}R(t)$ |
            | `Langevin` | 1, 2, 3 | $M \ddot{\mathbf{X}}$ | $=$ | $- \mathbf{\nabla} U(\mathbf{X})$ <br> $- \gamma M\dot{\mathbf{X}} + \sqrt{2 M \gamma k_{\rm B} T}\mathbf{R}(t)$ |
            | `NewtonDamped` | 1, 2, 3 | $m \ddot{\vec{x}}$ | $=$ | $\vec{F} - \lambda \dot{\vec{x}}$ |
            | `VertexMechanics` | 2 | | | <p style="color: rgba(250, 250, 250, 0.5);" >(Too much space needed)</p> |

---

## Forces
            ```rust
            pub trait Interaction<Pos, Vel, Force, Inf = ()> {

fn get_interaction_information(&self) -> Inf;

fn calculate_force_between(
                  &self,
                  own_pos: &Pos,
                  own_vel: &Vel,
                  ext_pos: &Pos,
                  ext_vel: &Vel,
                  ext_info: &Inf
              ) -> Option<Result<Force, CalcError>>;
            }
            ```

| Name | Dim | Equations |
            | --- |:---:|:---:|
            | `BoundLennardJones` | 1, 2, 3 | $U(r) = 4\epsilon\left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right]$ <br> $V(r)=\min(U(r), \beta)\theta(r-\xi)$ |
            | `BacteriaReactions` | 2 | $\vec{z} = \vec{x_1} - \vec{x_2} \hspace{1cm} r = \|\|\vec{z}\|\|$ <br> $\sigma = \frac{r}{r_1+r_2} \hspace{1.375cm} q=0.2$ <br> $\vec{F} = \epsilon \frac{1+q}{q+\sigma} \theta(\sigma - 1) \frac{1}{r}\vec{z}$ |

---

## Cellular Reactions

```rust
            pub trait CellularReactions<CI, CE = CI> {
              fn get_intracellular(&self) -> CI;

fn set_intracellular(&mut self, concentration_vector: CI);

fn calculate_intra_and_extracellular_reaction_increment(
                  &self,
                  internal_concentration_vector: &CI,
                  external_concentration_vector: &CE
              ) -> Result<(CI, CE), CalcError>;
            }
            ```

| Name | Equations |
            | --- | --- |
            | `BacteriaReactions` | $\dot{V_c} = \frac{u}{1+\nu_I [I]} \sigma [R]V_c$ <br> $\dot{[R]} = \sum\limits_c -\frac{u}{1+\nu_I [I]} \frac{V_c}{V_D} [R] \delta(x-x_c) + D_R\Delta [R]$ <br> $\dot{[I]} = \sum\limits_c b\frac{V_c}{V_D}\delta(x-x_c) + D_I\Delta [I]$ |
            | `NoReactions` | |

---

## Cycle

```rust
              pub trait Cycle<Cell, Float = f64> {

fn update_cycle(
                    rng: &mut ChaCha8Rng,
                    dt: &Float,
                    cell: &mut Cell
                ) -> Option<CycleEvent>;

fn divide(
                    rng: &mut ChaCha8Rng,
                    cell: &mut Cell
                ) -> Result\<Cell, DivisionError>;
              ```

```rust
                // Default method
                // Can be overwritten
                fn update_conditional_phased_death(
                  rng: &mut ChaCha8Rng,
                  dt: &Float,
                  cell: &mut Cell
                ) -> Result<bool, DeathError> {
                  Ok(true)
                }
              }
              ```

</div>

Usually implemented individually for every new Cell-Agent.
            <div class="fragment">

### Example
              ```rust
              impl Cycle<Cell> for Cell {
                  fn update_cycle(rng: &mut ChaCha8Rng, dt: &f64, cell: &mut Cell) -> Option<CycleEvent> {
                      // Increase the current age of the cell
                      cell.current_age += dt;

// If the cell is older than the current age, return a division event
                      if cell.current_age > cell.maximum_age {
                          return Some(CycleEvent::Division)
                      }
                      None
                  }
                  // ...
              }
              ```

</div>

---

## Physical Domain

#### Domain

#### Subdomain

```rust
              pub trait Domain<C, S> {
                type SubDomainIndex;
                type VoxelIndex;

fn get_all_voxel_indices(
                  &self
                ) -> Vec<Self::VoxelIndex>;
              ```

```rust
              pub trait SubDomain<C> {
                type VoxelIndex;

fn get_voxel_index_of(
                    &self,
                    cell: &C
                ) -> Result<Self::VoxelIndex, BoundaryError>;
              ```

```rust
                fn decompose(
                  self,
                  n_subdomains: NonZeroUsize,
                  cells: Vec<C>
                ) -> Result<
                  DecomposedDomain<Self::SubDomainIndex, S, C>,
                  DecomposeError
                >
                  where S: SubDomain<C>;
              }
              ```
              ```rust
                fn get_neighbor_voxel_indices(
                  &self,
                  voxel_index: &Self::VoxelIndex
                ) -> Vec<Self::VoxelIndex>;

fn apply_boundary(
                  &self,
                  cell: &mut C
                ) -> Result<(), BoundaryError>;

fn get_all_indices(
                  &self
                ) -> Vec<Self::VoxelIndex>;
              }
              ```

</div>

| Name | Dim | Description |
            | --- |:---:| --- |
            | `CartesianCuboid` | 1, 2, 3 | Cuboid $[a_0, b_0] \times \dots \times [a_d, b_d] \sub \mathbb{R}^d$ <br> with reflective boundary conditions. |

---

## Techniques
            - Domain decomposition
              - Split simulation domain into multiple subdomains
              - Execute upate steps in parallel
              - Possible due to locality of interactions
            - Storage Manager
              - Export parts or store complete simulation
              - xml, json or sled (embedded) database

---

## Features
            - Fine-grained control over parameters
            - Combine building-blocks to new Agents
            - Multiple backends possible
            - Clear separatin between concepts and Implementation
              - Defines "Language of ABMs"

Documentation

Cell Sorting

Start

End

Species-specific interaction
Only forces and mechanics
Just 4 Parameters

Radius $r$, Damping $\lambda$,
Interaction Range $\sigma$,
Potential Strength $V_0$

Think like a cell!

05.02.2024

Jonas Pleyer

Outline

Agent-Based Modelling

Documentation

Cell Sorting

Start

End

Semi-Vertex Model

(Still in Testing phase)

Bacterial Branching

Performance