3D Cell Sorting

Cell Sorting is a naturally occuring phenomenon which drives many biological processes [1,2]. While the underlying biological reality can be quite complex, it is rather simple to describe such a system in its most basic form. The underlying principle is that interactions between cells are specific.

Mathematical Description

We assume that cells are spherical objects which interact via force potentials. The two positions of cells are $x_i,x_j$ and their distance is $r=|x_i-x_j|$.

$$\begin{align} \sigma_{i,j} &= \frac{r}{R_i + R_j}\\ U(\sigma_{i,j}) &= V_0 \left(\frac{1}{3\sigma_{i,j}^3} - \frac{1}{\sigma_{i,j}}\right) \end{align}$$

The values $R_i,R_j$ are the radii of the cells ($i\neq j$) interacting with each other. For simplification, we can assume that they are identical $R_i=R_j=R$.

Furthermore, we assume that the equation of motion is given by

$$\begin{equation} \partial^2_t x = F - \lambda \partial_t x \end{equation}$$

where the first term is the usual force term $F = - \nabla V$ obtained by differentiating the given potential and the second term is a damping term which arises due to the cells being immersed inside a viscuous fluid.

ℹ️

Note that we set the mass to $m=1$ thus rescaling the involved parameters. This means, that units of $V_0$ and $\lambda$ are changing and they incorporate this property.

We can assume that interactions between cells are restricted to close ranges and thus enforce a cutoff $\xi\geq R_i+R_j$ for the interaction where the resulting force is identical to zero. We further assume that cells of different species do not attract each other but do repel. To describe this behaviour, we set the potential to zero when $r>R_i+R_j$ (i.e., $\sigma_{i,j}>1$) and both cells have distinct species type $s_i$. In total we are left with

$$\begin{equation} V(\sigma_{i,j}) = \begin{cases} 0 &\text{ if } \sigma_{i,j}\geq\xi/(R_i+R_j)\\ 0 &\text{ if } s_i\neq s_j \text{ and } \sigma_{i,j}\geq 1\\ U(\sigma_{i,j}) &\text{ else } \end{cases}. \end{equation}$$

Parameters

In total, we are left with only 4 parameters to describe our system.

Parameter	Symbol	Value
Cell Radius	$R_i$	$6.0 \text{ µm}$
Potential Strength	$V_0$	$2\text{ µm}^2\text{ }/\text{ min}^2$
Damping Constant	$\lambda$	$2\text{ min}^{-1}$
Interaction Range	$\xi$	$1.5 (R_i+R_j)=3R_i$

The following table shows additional values which are used to initialize the system. In total, 1600 cells with random initial positions and zero velocity were placed inside the domain.

Property	Symbol	Value
Time Stepsize	$\Delta t$	$0.2\text{ min}$
Time Steps	$N_t$	$10'000$
Domain Size	$L$	$110\text{ µm}$
Cells Species 1	$N_{C,1}$	$800$
Cells Species 2	$N_{C,2}$	$800$

The chosen total simulated time is thus $2000\text{ min}=33.33\text{ h}$.

Results

Initial State

Cells are initially placed randomly inside the cuboid simulation domain.

Movie

Final State

After the simulation has finished, the cells have assembled into connected regions of the same species.

Code

The code for this simulation and the visualization can be found in the examples folder of cellular_raza.

Full Code

cellular_raza-examples/cell_sorting/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


use cellular_raza::prelude::*;

use nalgebra::Vector3;
use num::Zero;
use rand::{Rng, SeedableRng};
use rand_chacha::ChaCha8Rng;
use serde::{Deserialize, Serialize};

pub const N_CELLS_1: usize = 800;
pub const N_CELLS_2: usize = 800;

pub const CELL_DAMPENING: f64 = 2.0;
pub const CELL_RADIUS: f64 = 6.0;

pub const CELL_MECHANICS_RELATIVE_INTERACTION_RANGE: f64 = 1.5;
pub const CELL_MECHANICS_POTENTIAL_STRENGTH: f64 = 2.0;

pub const DT: f64 = 0.2;
pub const N_TIMES: u64 = 10_000;
pub const SAVE_INTERVAL: u64 = 20;

pub const N_THREADS: usize = 4;

pub const DOMAIN_SIZE: f64 = 110.0;

#[derive(Clone, Debug, Deserialize, PartialEq, Serialize)]
enum Species {
    RedCell,
    BlueCell,
}

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

impl Interaction<Vector3<f64>, Vector3<f64>, Vector3<f64>, (f64, Species)>
    for CellSpecificInteraction
{
    fn calculate_force_between(
        &self,
        own_pos: &Vector3<f64>,
        _own_vel: &Vector3<f64>,
        ext_pos: &Vector3<f64>,
        _ext_vel: &Vector3<f64>,
        ext_info: &(f64, Species),
    ) -> Result<(Vector3<f64>, Vector3<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((Vector3::zeros(), Vector3::zeros()));
                        }
                        false => (own_pos - ext_pos).normalize(),
                    };
                    let r = self.cell_radius * min_relative_distance_to_center;
                    (r, dir)
                }
            };
        let (ext_radius, species) = ext_info;
        // 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;

        if *species == self.species {
            Ok((
                -repelling_force - attracting_force,
                repelling_force + attracting_force,
            ))
        } else {
            Ok((-repelling_force, repelling_force))
        }
    }

    fn get_interaction_information(&self) -> (f64, Species) {
        (self.cell_radius, self.species.clone())
    }
}

#[derive(CellAgent, Clone, Deserialize, Serialize)]
struct Cell {
    #[Interaction]
    interaction: CellSpecificInteraction,
    #[Mechanics]
    mechanics: NewtonDamped3D,
}

fn main() -> Result<(), SimulationError> {
    // Define the seed
    let mut rng = ChaCha8Rng::seed_from_u64(1);

    let cells = (0..N_CELLS_1 + N_CELLS_2)
        .map(|n| {
            let pos = Vector3::from([
                rng.gen_range(0.0..DOMAIN_SIZE),
                rng.gen_range(0.0..DOMAIN_SIZE),
                rng.gen_range(0.0..DOMAIN_SIZE),
            ]);
            Cell {
                mechanics: NewtonDamped3D {
                    pos,
                    vel: Vector3::zero(),
                    damping_constant: CELL_DAMPENING,
                    mass: 1.0,
                },
                interaction: CellSpecificInteraction {
                    species: match n <= N_CELLS_1 {
                        true => Species::BlueCell,
                        false => Species::RedCell,
                    },
                    potential_strength: CELL_MECHANICS_POTENTIAL_STRENGTH,
                    relative_interaction_range: CELL_MECHANICS_RELATIVE_INTERACTION_RANGE,
                    cell_radius: CELL_RADIUS,
                },
            }
        })
        .collect::<Vec<_>>();

    let domain = CartesianCuboid3New::from_boundaries_and_interaction_ranges(
        [0.0; 3],
        [DOMAIN_SIZE; 3],
        [CELL_MECHANICS_RELATIVE_INTERACTION_RANGE * CELL_RADIUS * 2.0; 3],
    )?;

    let time = cellular_raza::core::time::FixedStepsize::from_partial_save_steps(
        0.0,
        DT,
        N_TIMES,
        SAVE_INTERVAL,
    )?;
    let storage_builder = StorageBuilder::new().location("out/cell_sorting");

    let settings = cellular_raza::core::backend::chili::Settings {
        n_threads: N_THREADS.try_into().unwrap(),
        time,
        storage: storage_builder,
        show_progressbar: true,
    };

    run_simulation!(
        domain,
        agents: cells,
        settings,
        aspects: [Mechanics, Interaction]
    )?;
    Ok(())
}

References

[1] M. S. Steinberg, “Reconstruction of Tissues by Dissociated Cells,” Science, vol. 141, no. 3579. American Association for the Advancement of Science (AAAS), pp. 401–408, Aug. 02, 1963. doi: 10.1126/science.141.3579.401.

[2] F. Graner and J. A. Glazier, “Simulation of biological cell sorting using a two-dimensional extended Potts model,” Physical Review Letters, vol. 69, no. 13. American Physical Society (APS), pp. 2013–2016, Sep. 28, 1992. doi: 10.1103/physrevlett.69.2013.