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.
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
.
Cargo
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
[package]
name = "cr_cell_sorting"
version = "0.1.0"
edition = "2021"
authors = ["Jonas Pleyer <jonas.pleyer@fdm.uni-freiburg.de>"]
[dependencies]
serde = { workspace = true, features=["rc"] }
rand = { workspace = true, features=["small_rng"] }
rand_chacha = { workspace = true }
nalgebra = { workspace = true }
cellular_raza = { path="../../cellular_raza" }
num = { workspace = true }
[[bin]]
name = "cr_cell_sorting_default"
path = "src/main.rs"
[[bin]]
name = "cr_cell_sorting_langevin"
path = "src/main-langevin.rs"
[[bin]]
name = "cr_cell_sorting_brownian"
path = "src/main-brownian.rs"
workspace = true
or
path="../../"
should be replaced with the versions used in the workspace
Cargo.toml
.
Full Code
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.