Bacterial Rods
Mathematical Description
To model the spatial mechanics of elongated bacteria, we represent them as a collection of auxiliary vertices $\{\vec{v}_i\}$ which are connected by springs in ascending order. Furthermore, we assume that the cells are flexible described by their stiffness property. A force $\vec{F}$ interacting between cellular agents determines the radius (thickness) of the rods and an attractive component can model adhesion between cells.
Mechanics
In principle we can assign individual lengths $\{l_i\}$ and strengths $\{\gamma\}_i$ to each spring. The internal force acting on vertex $\vec{v}_i$ can be divided into 2 contributions coming from the 2 springs pulling on it. In the case when $i=0,N_\text{vertices}$, this is reduced to only one internal component. We denote with $\vec{c}_{i}$ the connection between two vertices $$\begin{align} \vec{c}_i = \vec{v}_{i}-\vec{v}_{i-1} \end{align}$$ and can write down the resulting force $$\begin{align} \vec{F}_{i,\text{springs}} = &-\gamma_i\left(1 - \frac{l_i}{\left|\vec{c}_i\right|}\right) \vec{c}_i\\ &+ \gamma_{i+1}\left(1 - \frac{l_{i+1}}{\left|\vec{c}_{i+1}\right|}\right) \vec{c}_{i+1} \end{align}$$
In addition to springs between individual vertices $\vec{v}_i$, we assume that each angle at a vertex between two other is subject to a stiffening force. Assuming that $\alpha_i$ is the angle between the connections and $\vec{d}_i=\vec{c}_i/|\vec{c}_i|$ is the normalized connection, we can write down the forces acting on vertices $\vec{v}_i,\vec{v}_{i-1},\vec{v}_{i+1}$ $$\begin{align} \vec{F}_{i,\text{stiffness}} &= \eta_i\left(\pi-\alpha_i\right) \frac{\vec{d}_i - \vec{d}_{i+1}}{|\vec{d}_i-\vec{d}_{i+1}|}\\ \vec{F}_{i-1,\text{stiffness}} &= -\frac{1}{2}\vec{F}_{i,\text{stiffness}}\\ \vec{F}_{i+1,\text{stiffness}} &= -\frac{1}{2}\vec{F}_{i,\text{stiffness}} \end{align}$$ where $\eta_i$ is the angle stiffness at vertex $\vec{v}_i$. We can see that the stiffening force does not move the overall center of the cell in space. The total force is the sum of external and interal forces. $$\begin{equation} \vec{F}_{i,\text{total}} = \vec{F}_{i,\text{springs}}+ \vec{F}_{i,\text{stiffness}} + \vec{F}_{i,\text{external}} \end{equation}$$ and are integrated via $$\begin{align} \partial_t^2 \vec{x} &= \partial\vec{x} + D\vec{\xi}\\ \partial_t\vec{x} &= \vec{F}_\text{total} \end{align}$$ where $D$ is the diffusion constant and $\vec{\xi}$ is the wiener process (compare with brownian motion).
Interaction
When calculating forces acting between the cells, we can use a simplified model to circumvent the numerically expensive integration over the complete length of the rod. Given a vertex $\vec{v}_i$ on one cell, we calculate the closest point $\vec{p}$ on the polygonal line given by the vertices $\{\vec{w}_j\}$ of the interacting cell. Furthermore we determine the value $q\in[0,1]$ such that $$\begin{equation} \vec{p} = (1-q)\vec{w}_j + q\vec{w}_{j+1} \end{equation}$$ for some specific $j$. The force is then calculated between the points $\vec{v}_i$ and $\vec{p}_i$ and acts on the vertex $\vec{w}_i,\vec{w}_{i+1}$ with relative strength $(1-q)$ and $q$. $$\begin{align} \vec{F}_{i,\text{External}} = \vec{F}(\vec{v}_i,\vec{p}) \end{align}$$ For this example, we reused the interaction shape of the cell-sorting example ignoring the species aspect.
Cycle
To simulate proliferation, we introduce a growth term for the spring lengths $l_i$ $$\begin{equation} \partial_t l_i = \mu \end{equation}$$ which will increase the length of the cell indefenitely unless we impose a condition for the division event. We define a threshold (in our case double of the original length) for the total length of the cell at which it divides. To construct a new cell, we cannot simply copy the existing one twice, but we also need to adjust internal parameters in the process. The following actions need to be taken for the old and new agent.
- Assign a new growth rate (pick randomly from uniform distribution in $[0.8\mu_0,1.2\mu_0]$ where $\mu_0$ is some fixed value)
- Assign new positions
- Calculate new spring lengths $$\tilde{l}_i = l_i\left(\frac{1}{2} - \frac{r}{\sum\limits_i l_i}\right)$$
- Calculate middle of old cell $$\vec{m} = \frac{1}{N_\text{vertices}}\sum\limits_i\vec{v}_i$$
- Calculate positions of new vertices $\vec{w}_i$ $$\begin{align} q_i &= \frac{i}{N_\text{vertices}}\\ \vec{w}_{i,\text{new},\pm} &= (1-q_i)\vec{m} + q_i(\vec{v}_{\pm\text{start}} - \vec{m}) \end{align}$$
Domain
In this simulation example, the domain plays an important role. The domain consists of an elongated box with reflective boundary conditions. Cells are initially placed in the left part. Due to their repulsive potential at short distances, they begin to push each other into the remaining space.
Results
Initial Snapshot
Movie
Final Snapshot
Code
The code is part of the examples and can be found in the official github repository under bacterial-rods.
Custom Domain
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
use cellular_raza::building_blocks::{CartesianCuboid, CartesianSubDomain};
use cellular_raza::concepts::*;
use serde::Serialize;
use crate::Agent;
#[derive(Clone, Domain)]
pub struct MyDomain<const D2: usize> {
#[DomainRngSeed]
pub cuboid: CartesianCuboid<f64, D2>,
pub gravity: f64,
pub damping: f64,
}
impl<const D2: usize> DomainCreateSubDomains<MySubDomain<D2>> for MyDomain<D2> {
type SubDomainIndex = usize;
type VoxelIndex = [usize; D2];
fn create_subdomains(
&self,
n_subdomains: std::num::NonZeroUsize,
) -> Result<
impl IntoIterator<Item = (Self::SubDomainIndex, MySubDomain<D2>, Vec<Self::VoxelIndex>)>,
DecomposeError,
> {
let gravity = self.gravity;
let damping = self.damping;
let subdomains = self.cuboid.create_subdomains(n_subdomains)?;
Ok(subdomains
.into_iter()
.map(move |(subdomain_index, subdomain, voxels)| {
(
subdomain_index,
MySubDomain {
subdomain,
force: MySubDomainForce { gravity, damping },
},
voxels,
)
}))
}
}
impl<const D1: usize, const D2: usize> SortCells<Agent<D1, D2>> for MyDomain<D2> {
type VoxelIndex = [usize; D2];
fn get_voxel_index_of(&self, cell: &Agent<D1, D2>) -> Result<Self::VoxelIndex, BoundaryError> {
let pos = cell.pos();
let index = (pos.row_mean().transpose() - self.cuboid.get_min())
.component_div(&self.cuboid.get_dx());
let res: [usize; D2] = index.try_cast::<usize>().unwrap().into();
Ok(res)
}
}
#[derive(Clone, SubDomain, Serialize)]
pub struct MySubDomain<const D2: usize> {
#[Base]
pub subdomain: CartesianSubDomain<f64, D2>,
#[Force]
pub force: MySubDomainForce,
}
impl<const D1: usize, const D2: usize>
SubDomainMechanics<nalgebra::SMatrix<f64, D1, D2>, nalgebra::SMatrix<f64, D1, D2>>
for MySubDomain<D2>
{
fn apply_boundary(
&self,
pos: &mut nalgebra::SMatrix<f64, D1, D2>,
vel: &mut nalgebra::SMatrix<f64, D1, D2>,
) -> Result<(), BoundaryError> {
// TODO refactor this with matrix multiplication!!!
// This will probably be much more efficient and less error-prone!
// For each position in the springs Agent<D1, D2>
pos.row_iter_mut()
.zip(vel.row_iter_mut())
.for_each(|(mut p, mut v)| {
// For each dimension in the space
for i in 0..p.ncols() {
// Check if the particle is below lower edge
if p[i] < self.subdomain.get_domain_min()[i] {
p[i] = 2.0 * self.subdomain.get_domain_min()[i] - p[i];
v[i] = v[i].abs();
}
// Check if the particle is over the edge
if p[i] > self.subdomain.get_domain_max()[i] {
p[i] = 2.0 * self.subdomain.get_domain_max()[i] - p[i];
v[i] = -v[i].abs();
}
}
});
// If new pos is still out of boundary return error
for j in 0..pos.nrows() {
let p = pos.row(j);
for i in 0..pos.ncols() {
if p[i] < self.subdomain.get_domain_min()[i]
|| p[i] > self.subdomain.get_domain_max()[i]
{
return Err(BoundaryError(format!(
"Particle is out of domain at pos {:?}",
pos
)));
}
}
}
Ok(())
}
}
#[derive(Clone, Serialize)]
pub struct MySubDomainForce {
pub damping: f64, // 0.75 / SECOND,
pub gravity: f64,
}
impl<const D1: usize, const D2: usize>
SubDomainForce<
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
> for MySubDomainForce
{
fn calculate_custom_force(
&self,
_pos: &nalgebra::SMatrix<f64, D1, D2>,
vel: &nalgebra::SMatrix<f64, D1, D2>,
) -> Result<nalgebra::SMatrix<f64, D1, D2>, cellular_raza::concepts::CalcError> {
let mut force = nalgebra::SMatrix::<f64, D1, D2>::zeros();
// Gravity
force.column_mut(D2 - 1).add_scalar_mut(-self.gravity);
// Damping force
force -= self.damping * vel;
Ok(force)
}
}
impl<const D1: usize, const D2: usize> SortCells<Agent<D1, D2>> for MySubDomain<D2> {
type VoxelIndex = [usize; D2];
fn get_voxel_index_of(&self, cell: &Agent<D1, D2>) -> Result<Self::VoxelIndex, BoundaryError> {
let pos: [f64; D2] = cell.pos().row_mean().transpose().into();
let index = self.subdomain.get_index_of(pos)?;
Ok(index)
}
}
Main Simulation
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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
use std::ops::Div;
use cellular_raza::building_blocks::{nearest_point_from_point_to_line, CartesianCuboid};
use cellular_raza::concepts::{CalcError, CellAgent, Cycle, CycleEvent, Position, Velocity};
use cellular_raza::core::backend::chili;
use rand::SeedableRng;
use serde::{Deserialize, Serialize};
mod custom_domain;
use custom_domain::*;
pub const METRE: f64 = 1.0;
pub const MILI_METRE: f64 = 1e-3;
pub const MICRO_METRE: f64 = 1e-6;
pub const SECOND: f64 = 1.0;
pub const MINUTE: f64 = 60.0 * SECOND;
pub const HOUR: f64 = 24.0 * MINUTE;
/// # Mechanics
/// The cell consists of two 2D points which are coupled via a spring.
/// We represent this by a matrix
/// ```
/// # use nalgebra::Matrix5x2;
/// let cell_pos = Matrix5x2::new(
/// 1.0, -2.0,
/// -2.0, 1.0
/// );
/// let cell_pos_point_0 = cell_pos.row(0);
/// let cell_pos_point_1 = cell_pos.row(1);
/// ```
#[derive(CellAgent, Clone, Deserialize, Serialize)]
pub struct Agent<const D1: usize, const D2: usize> {
// Mechanics
pos: nalgebra::SMatrix<f64, D1, D2>,
vel: nalgebra::SMatrix<f64, D1, D2>,
random_velocity: nalgebra::SMatrix<f64, D1, D2>,
diffusion_constant: f64,
spring_tension: f64,
angle_stiffness: f64,
spring_length: f64,
// Interaction
interaction_potential: f64,
radius: f64,
interaction_range: f64,
// Cycle
growth_rate: f64,
max_spring_length: f64,
}
impl<const D1: usize, const D2: usize>
cellular_raza::concepts::Mechanics<
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
> for Agent<D1, D2>
{
fn calculate_increment(
&self,
force: nalgebra::SMatrix<f64, D1, D2>,
) -> Result<
(
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
),
CalcError,
> {
use core::ops::AddAssign;
// Calculate internal force between the two points of the Agent
let mut total_force = force;
// Calculate force exerted by spring action between individual vertices
let dist_internal =
self.pos.rows(0, self.pos.nrows() - 1) - self.pos.rows(1, self.pos.nrows() - 1);
dist_internal.row_iter().enumerate().for_each(|(i, dist)| {
let dir = dist.normalize();
let force_internal = -self.spring_tension * (dist.norm() - self.spring_length) * dir;
total_force.row_mut(i).add_assign(0.5 * force_internal);
total_force.row_mut(i + 1).add_assign(-0.5 * force_internal);
});
// Calculate force exerted by angle-contributions
use itertools::Itertools;
dist_internal
.row_iter()
.tuple_windows::<(_, _)>()
.enumerate()
.for_each(|(i, (conn1, conn2))| {
let angle = conn1.angle(&-conn2);
let force_direction = (conn1.normalize() - conn2.normalize()).normalize();
let force = self.angle_stiffness * (std::f64::consts::PI - angle) * force_direction;
total_force.row_mut(i).add_assign(-0.5 * force);
total_force.row_mut(i + 1).add_assign(force);
total_force.row_mut(i + 2).add_assign(-0.5 * force);
});
Ok((self.vel.clone() + self.random_velocity, total_force))
}
fn get_random_contribution(
&self,
rng: &mut rand_chacha::ChaCha8Rng,
dt: f64,
) -> Result<
(
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
),
cellular_raza::prelude::RngError,
> {
let distr = match rand_distr::Normal::new(0.0, dt.sqrt()) {
Ok(e) => Ok(e),
Err(e) => Err(cellular_raza::concepts::RngError(format!("{e}"))),
}?;
let dpos = std::f64::consts::SQRT_2
* self.diffusion_constant
* nalgebra::SMatrix::<f64, D1, D2>::from_distribution(&distr, rng)
/ dt;
let dvel = nalgebra::SMatrix::<f64, D1, D2>::zeros();
Ok((dpos, dvel))
}
}
impl<const D1: usize, const D2: usize> Position<nalgebra::SMatrix<f64, D1, D2>> for Agent<D1, D2> {
fn pos(&self) -> nalgebra::SMatrix<f64, D1, D2> {
self.pos.clone()
}
fn set_pos(&mut self, pos: &nalgebra::SMatrix<f64, D1, D2>) {
self.pos = pos.clone();
}
}
impl<const D1: usize, const D2: usize> Velocity<nalgebra::SMatrix<f64, D1, D2>> for Agent<D1, D2> {
fn velocity(&self) -> nalgebra::SMatrix<f64, D1, D2> {
self.vel.clone()
}
fn set_velocity(&mut self, velocity: &nalgebra::SMatrix<f64, D1, D2>) {
self.vel = velocity.clone();
}
}
impl<const D1: usize, const D2: usize>
cellular_raza::concepts::Interaction<
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
f64,
> for Agent<D1, D2>
{
fn calculate_force_between(
&self,
own_pos: &nalgebra::SMatrix<f64, D1, D2>,
_own_vel: &nalgebra::SMatrix<f64, D1, D2>,
ext_pos: &nalgebra::SMatrix<f64, D1, D2>,
_ext_vel: &nalgebra::SMatrix<f64, D1, D2>,
ext_radius: &f64,
) -> Result<
(
nalgebra::SMatrix<f64, D1, D2>,
nalgebra::SMatrix<f64, D1, D2>,
),
CalcError,
> {
use core::ops::AddAssign;
let mut force_own = nalgebra::SMatrix::<f64, D1, D2>::zeros();
let mut force_ext = nalgebra::SMatrix::<f64, D1, D2>::zeros();
use itertools::Itertools;
for (i, p1) in own_pos.row_iter().enumerate() {
for (j, (p2_n0, p2_n1)) in ext_pos.row_iter().tuple_windows::<(_, _)>().enumerate() {
// Calculate the closest point of the external position
let (_, nearest_point, rel_length) = nearest_point_from_point_to_line(
&p1.transpose(),
&(p2_n0.transpose(), p2_n1.transpose()),
);
let dist = p1 - nearest_point.transpose();
let r = dist.norm();
if r < ext_radius + self.radius + self.interaction_range {
let sigma = r / (self.radius + ext_radius);
let strength = (1.0 / sigma.powf(4.0) - 1.0 / sigma.powf(2.0)).min(0.2);
let force_strength = self.interaction_potential * strength * dist.normalize();
force_own.row_mut(i).add_assign(-force_strength);
force_ext
.row_mut(j)
.add_assign(-(1.0 - rel_length) / 2.0 * force_strength);
force_ext
.row_mut((j + 1) % D1)
.add_assign(-rel_length / 2.0 * force_strength);
}
}
}
Ok((-force_own, force_ext))
}
fn get_interaction_information(&self) -> f64 {
self.radius
}
fn is_neighbor(
&self,
own_pos: &nalgebra::SMatrix<f64, D1, D2>,
ext_pos: &nalgebra::SMatrix<f64, D1, D2>,
_ext_inf: &f64,
) -> Result<bool, CalcError> {
for own_point in own_pos.row_iter() {
for ext_point in ext_pos.row_iter() {
if (own_point - ext_point).norm() < 2.0 * self.radius {
return Ok(true);
}
}
}
Ok(false)
}
fn react_to_neighbors(&mut self, neighbors: usize) -> Result<(), CalcError> {
let base_rate = 3.0 * MICRO_METRE / MINUTE;
if neighbors > 0 {
self.growth_rate = (base_rate * (10.0 - neighbors as f64) / 10.0).max(0.0);
} else {
self.growth_rate = base_rate;
}
Ok(())
}
}
impl<const D1: usize, const D2: usize> Cycle<Agent<D1, D2>> for Agent<D1, D2> {
fn update_cycle(
_rng: &mut rand_chacha::ChaCha8Rng,
dt: &f64,
cell: &mut Agent<D1, D2>,
) -> Option<cellular_raza::prelude::CycleEvent> {
if cell.spring_length < cell.max_spring_length {
cell.spring_length += dt * cell.growth_rate;
None
} else {
Some(CycleEvent::Division)
}
}
fn divide(
rng: &mut rand_chacha::ChaCha8Rng,
cell: &mut Agent<D1, D2>,
) -> Result<Agent<D1, D2>, cellular_raza::prelude::DivisionError> {
use rand::Rng;
let c1 = cell;
let mut c2 = c1.clone();
let base_rate = 1.5 * MICRO_METRE / MINUTE;
c1.growth_rate = rng.gen_range(0.8 * base_rate..1.2 * base_rate);
c2.growth_rate = rng.gen_range(0.8 * base_rate..1.2 * base_rate);
let n_rows = c1.pos.nrows();
// Calculate the fraction of how much we need to scale down the individual spring lengths
// in order for the distances to still work.
let div_factor = 0.5 - c1.radius / (n_rows as f64 * c1.spring_length);
// Shrink spring length
c1.spring_length *= div_factor;
c2.spring_length *= div_factor;
// Define new positions
let middle = if D1 % 2 == 0 {
1.0 * c1.pos.row(D1.div(2))
} else {
0.5 * (c1.pos.row(D1.div(2)) + c1.pos.row(D1.div(2) + 1))
};
let p1 = c1.pos.row(0).to_owned();
let p2 = c1.pos.row(n_rows - 1).to_owned();
let d1 = (p1 - middle) * (1.0 - div_factor);
let d2 = (p2 - middle) * (1.0 - div_factor);
c1.pos.row_iter_mut().for_each(|mut r| {
r -= middle;
r *= div_factor;
r += middle;
r += d1
});
c2.pos.row_iter_mut().for_each(|mut r| {
r -= middle;
r *= div_factor;
r += middle;
r += d2;
});
Ok(c2)
}
}
fn main() -> Result<(), chili::SimulationError> {
// Define the dimensionality of the problem
const D1: usize = 5;
const D2: usize = 3;
// Define initial random seed
use rand::Rng;
let mut rng = rand_chacha::ChaCha8Rng::seed_from_u64(5);
// Give agent default values
let agent = Agent {
pos: nalgebra::SMatrix::<f64, D1, D2>::zeros(),
vel: nalgebra::SMatrix::<f64, D1, D2>::zeros(),
random_velocity: nalgebra::SMatrix::<f64, D1, D2>::zeros(),
diffusion_constant: 0.1 * MICRO_METRE.powf(2.0) / SECOND,
spring_tension: 20.0 / SECOND.powf(2.0),
angle_stiffness: 2.0 * MICRO_METRE / SECOND.powf(2.0),
interaction_potential: 3e6 * MICRO_METRE.powf(2.0) / SECOND.powf(2.0),
spring_length: 3.0 * MICRO_METRE,
max_spring_length: 6.0 * MICRO_METRE,
radius: 3.0 * MICRO_METRE,
interaction_range: 0.25 * MICRO_METRE,
growth_rate: 3.0 * MICRO_METRE / MINUTE,
};
// Place agents in simulation domain
let domain_size = 50.0 * MICRO_METRE;
let delta_x = agent.spring_length * D1 as f64;
let agents = (0..5).map(|_| {
let mut new_agent = agent.clone();
new_agent.spring_length = rng.gen_range(1.5..2.5) * MICRO_METRE;
let mut pos = nalgebra::SMatrix::<f64, D1, D2>::zeros();
pos[(0, 0)] = rng.gen_range(delta_x..2.0 * delta_x);
pos[(0, 1)] = rng.gen_range(delta_x / 3.0..delta_x * 2.0 / 3.0);
pos[(0, 2)] = rng.gen_range(delta_x..2.0 * delta_x);
let theta = rng.gen_range(0.0..2.0 * std::f64::consts::PI);
for i in 1..pos.nrows() {
let phi =
theta + rng.gen_range(-std::f64::consts::FRAC_PI_8..std::f64::consts::FRAC_PI_8);
let mut direction = nalgebra::SVector::<f64, D2>::zeros();
if D2 > 0 {
direction[0] = phi.cos();
}
if D2 > 1 {
direction[1] = phi.sin();
}
let new_pos = pos.row(i - 1) + agent.spring_length * (direction).transpose();
use core::ops::AddAssign;
pos.row_mut(i).add_assign(new_pos);
}
new_agent.set_pos(&pos);
new_agent
});
// Domain Setup
let domain_sizes = [4.0 * domain_size, delta_x, 3.0 * delta_x];
let domain_segments = [8, 1, 2];
let domain = MyDomain {
cuboid: CartesianCuboid::from_boundaries_and_n_voxels(
[0.0; D2],
domain_sizes,
domain_segments,
)?,
// TODO enable this
gravity: 0.0 * 2e-7 * 9.81 * METRE / SECOND.powf(2.0),
damping: 1.5 / SECOND,
};
// Storage Setup
let storage_builder = cellular_raza::prelude::StorageBuilder::new().location("./out");
// Time Setup
let t0 = 0.0 * MINUTE;
let dt = 0.05 * SECOND;
let save_interval = 2.0 * SECOND;
let t_max = 30.0 * MINUTE;
let time_stepper = cellular_raza::prelude::time::FixedStepsize::from_partial_save_interval(
t0,
dt,
t_max,
save_interval,
)?;
let settings = chili::Settings {
n_threads: 8.try_into().unwrap(),
time: time_stepper,
storage: storage_builder,
show_progressbar: true,
};
println!("Running Simulation");
let _storage = chili::run_simulation!(
domain: domain,
agents: agents,
settings: settings,
aspects: [Mechanics, Interaction, DomainForce, Cycle],
)?;
Ok(())
}
Plotting
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
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
import pyvista as pv
import numpy as np
import argparse
import itertools
import tqdm
import concurrent.futures
from pathlib import Path
import glob
import json
import os
import pandas as pd
def get_last_output_path():
return Path(sorted(glob.glob("out/*"))[-1])
def get_all_iterations(output_path: Path | None = None):
if output_path is None:
output_path = get_last_output_path()
folders = glob.glob(str(output_path / "cells/json") + "/*")
iterations = [int(Path(fi).name) for fi in folders]
return output_path, iterations
def load_cells_from_iteration(output_path: Path, iteration: int):
load_path = Path(output_path) / "cells/json/{:020}".format(iteration)
results = []
for file in glob.glob(str(load_path) + "/*"):
f = open(file)
elements = json.load(f)["data"]
elements = [element["element"][0] for element in elements]
results.extend(elements)
df = pd.json_normalize(results)
if len(df) > 0:
# Format individual entries for easier use later on
df["identifier"] = df["identifier"].apply(lambda x: tuple(x))
df["cell.pos"] = df["cell.pos"].apply(lambda x: np.array(x, dtype=float).reshape(3, -1))
df["cell.vel"] = df["cell.vel"].apply(lambda x: np.array(x, dtype=float))
return df
def get_cell_meshes(iteration: int, path: Path):
cells = load_cells_from_iteration(path, iteration)
positions = np.array([x for x in cells["cell.pos"]], dtype=float)
radii = np.array([x for x in cells["cell.radius"]], dtype=float)
growth_rate = np.array([x for x in cells["cell.growth_rate"]], dtype=float)
cell_surfaces = []
for i, p in enumerate(positions):
meshes = []
# Add the sphere at the first point
meshes.append(pv.Sphere(center=p[:,0], radius=radii[i]))
for j in range(max(p.shape[1]-1,0)):
# Add a sphere for each following point
meshes.append(pv.Sphere(center=p[:,j+1], radius=radii[i]))
# Otherwise add cylinders
pos1 = p[:,j]
pos2 = p[:,j+1]
center = 0.5 * (pos1 + pos2)
direction = pos2 - center
radius = radii[i]
height = float(np.linalg.norm(pos1 - pos2))
# sphere = pv.Sphere(center=pos, radius=radii[i])
# meshes.append(sphere)
cylinder = pv.Cylinder(center, direction, radius, height)
meshes.append(cylinder)
combined = pv.MultiBlock(meshes).combine()
if combined is None:
return None
merged = combined.extract_surface()
if merged is None:
return None
merged = merged.clean()
cell_surfaces.append((merged, growth_rate[i]))
return cell_surfaces
def plot_spheres(
iteration: int,
path: Path = Path("./"),
opath: Path | None = None,
overwrite: bool = False,
transparent_background: bool = False,
):
if opath is None:
opath = path / "images/{:010}.png".format(iteration)
opath.parent.mkdir(parents=True, exist_ok=True)
if os.path.isfile(opath) and overwrite is False:
return None
cell_meshes = get_cell_meshes(iteration, path)
cell_meshes = cell_meshes if cell_meshes is not None else []
# General Settings
plotter = pv.Plotter(off_screen=True)
plotter.set_background([100, 100, 100])
# Draw box around everything
dx = 3e-6
box = pv.Box(bounds=(-dx, 200e-6+dx, 0-dx, 15e-6+dx, 0-dx, 45e-6+dx))
plotter.add_mesh(box, style="wireframe", line_width=15)
color_min = np.array([69, 124, 214])
color_max = np.array([82, 191, 106])
for (cell, growth_rate) in cell_meshes:
# TODO MAGIC NUMBERS
q = growth_rate / (3e-6 / 60)
color = (1-q) * color_min + q * color_max
plotter.add_mesh(
cell,
show_edges=False,
color=(int(color[0]), int(color[1]), int(color[2])),
)
# Define camera
# TODO MAGIC NUMBERS
plotter.camera.position = (100e-6, -250e-6, -250e-6)
plotter.camera.focal_point = (100e-6, 25e-6, 22.5e-6)
plotter.enable_ssao(radius=12)
plotter.enable_anti_aliasing()
img = plotter.screenshot(opath, transparent_background=transparent_background)
plotter.close()
return img
def __plot_spheres_helper(args):
(args, kwargs) = args
plot_spheres(*args, **kwargs)
def plot_all_spheres(
path: Path,
n_threads: int | None = None,
overwrite:bool=False,
transparent_background: bool = False
):
iterations = [it for it in get_all_iterations(path)[1]]
if n_threads is None:
n_cpu = os.cpu_count()
if n_cpu is not None and n_cpu > 2:
n_threads = n_cpu -2
else:
n_threads = 1
args = [(it,) for it in iterations]
kwargs = {
"path": path,
"overwrite": overwrite,
"transparent_background": transparent_background,
}
with concurrent.futures.ProcessPoolExecutor(max_workers=n_threads) as executor:
_ = list(tqdm.tqdm(executor.map(
__plot_spheres_helper,
zip(args, itertools.repeat(kwargs))
), total=len(iterations)))
# TODO actually use this
if __name__ == "_main__":
parser = argparse.ArgumentParser(
prog="Plotter",
description="A plotting CLI for the springs example",
epilog="For suggestions or bug-reports please refer to\
https://github.com/jonaspleyer/cellular_raza/"
)
parser.add_argument(
"-n",
"--iteration",
help="Specify which iteration to plot.\
Multiple arguments are accepted as a list.\
If not specified, plot all iterations."
)
parser.add_argument(
"-i",
"--input",
help="Input path of files. If not given,\
it will be determined automatically by searching in './out/'"
)
parser.add_argument(
"-o",
"--output-path",
help="Path where to store generated images."
)
args = parser.parse_args()
if __name__ == "__main__":
print("Plotting Individual Snapshots")
output_path = get_last_output_path()
plot_all_spheres(output_path, transparent_background=True)
print("Generating Movie")
bashcmd = f"ffmpeg\
-v quiet\
-stats\
-y\
-r 30\
-f image2\
-pattern_type glob\
-i '{output_path}/images/*.png'\
-c:v h264\
-pix_fmt yuv420p\
-strict -2 {output_path}/movie.mp4"
os.system(bashcmd)
print("Playing Movie")
bashcmd2 = f"firefox {output_path}/movie.mp4"
os.system(bashcmd2)