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
//! Model description: cell types, the contact matrix and the constraints.
/// A cell type's own parameters. Type 0 is the medium and takes no volume
/// constraint.
///
/// Every term is off at zero, so `CellType::default()` is the medium and a
/// type states only what it uses. New terms therefore leave existing
/// descriptions and existing callers alone.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct CellType {
/// Target volume in lattice sites.
pub target_volume: f64,
/// Volume constraint strength.
pub lambda_volume: f64,
/// Target surface in boundary bonds.
pub target_surface: f64,
/// Surface constraint strength.
pub lambda_surface: f64,
/// Volume at which a cell divides. Zero never divides.
pub division_volume: f64,
/// Probability per step that a cell of this type dies. Zero never dies.
pub death_rate: f64,
/// Target major axis in sites. Zero leaves the length unconstrained.
pub target_length: f64,
/// Length constraint strength.
pub lambda_length: f64,
/// Whether a copy that would locally pinch this cell is refused.
pub connected: bool,
/// Steps of memory a site keeps after this cell takes it. Zero leaves the
/// cell with no persistence.
pub max_activity: f64,
/// Strength of the persistence.
pub lambda_activity: f64,
/// A constant drift, as a force per axis. The work is the displacement
/// along it, so a cell with one travels that way whatever its neighbours
/// do.
pub external: [f64; 3],
}
/// The whole model: what a Blueprint would state, in one struct.
#[derive(Clone, Debug)]
pub struct Model {
/// Diffusing species, in the order a description lists them.
pub species: Vec<crate::field::Species>,
/// Secretion and uptake per cell type, index 0 the medium.
pub exchange: Vec<crate::field::Exchange>,
/// Chemotactic sensitivity per cell type and species, index 0 the medium.
/// Positive climbs a gradient.
pub chemotaxis: Vec<Vec<f64>>,
/// Lattice width in sites.
pub width: usize,
/// Lattice height in sites.
pub height: usize,
/// Lattice depth in sites. One is a plane, and every neighbourhood then
/// drops its out-of-plane offsets.
pub depth: usize,
/// Contact energies, row-major over `(type_a, type_b)`, symmetric.
pub contact: Vec<f64>,
/// One entry per type, index 0 the medium.
pub types: Vec<CellType>,
/// Metropolis fluctuation amplitude, the Potts temperature.
pub temperature: f64,
/// 1 for the four von Neumann neighbours, 2 for the eight Moore ones.
pub neighbour_order: u8,
/// Random seed.
pub seed: u64,
}
impl Model {
/// Number of cell types, counting the medium.
#[must_use]
pub fn n_types(&self) -> usize {
self.types.len()
}
/// How many dimensions the lattice spans.
#[must_use]
pub fn dimensions(&self) -> usize {
if self.depth > 1 { 3 } else { 2 }
}
/// Whether any type keeps a memory of where it has been.
#[must_use]
pub fn has_motility(&self) -> bool {
self.types
.iter()
.any(|t| t.lambda_activity != 0.0 && t.max_activity > 0.0)
}
/// Whether any type drifts.
#[must_use]
pub fn has_external_potential(&self) -> bool {
self.types
.iter()
.any(|t| t.external.iter().any(|&v| v != 0.0))
}
/// Whether any type refuses a copy that would pinch it.
#[must_use]
pub fn has_connectivity(&self) -> bool {
self.types.iter().any(|t| t.connected)
}
/// Whether any type constrains its length, so a run without one skips the
/// moment bookkeeping.
#[must_use]
pub fn has_length_constraint(&self) -> bool {
self.types.iter().any(|t| t.lambda_length != 0.0)
}
/// Whether any type follows a gradient.
#[must_use]
pub fn has_chemotaxis(&self) -> bool {
self.chemotaxis
.iter()
.any(|row| row.iter().any(|&l| l != 0.0))
}
/// Whether any type divides or dies, so a run without either skips the
/// per-step scan over labels.
#[must_use]
pub fn has_population_events(&self) -> bool {
self.types
.iter()
.any(|t| t.division_volume > 0.0 || t.death_rate > 0.0)
}
/// Contact energy between two types.
#[must_use]
pub fn contact_energy(&self, a: u8, b: u8) -> f64 {
self.contact[a as usize * self.n_types() + b as usize]
}
/// Check the shapes agree before a run starts.
///
/// # Errors
/// A message naming the first inconsistency found.
pub fn validate(&self) -> Result<(), String> {
let n = self.n_types();
if n == 0 {
return Err("a model needs at least the medium type".into());
}
if self.contact.len() != n * n {
return Err(format!(
"contact matrix is {} entries for {n} types, want {}",
self.contact.len(),
n * n
));
}
for a in 0..n {
for b in 0..n {
let (ab, ba) = (self.contact[a * n + b], self.contact[b * n + a]);
if (ab - ba).abs() > 1e-12 {
return Err(format!("contact matrix is asymmetric at ({a}, {b})"));
}
}
}
if self.width == 0 || self.height == 0 || self.depth == 0 {
return Err("the lattice has a zero dimension".into());
}
if !matches!(self.neighbour_order, 1..=3) {
return Err(format!(
"neighbour order {} is outside the range one to three",
self.neighbour_order
));
}
if self.temperature <= 0.0 {
return Err("the temperature must be positive".into());
}
if !self.exchange.is_empty() && self.exchange.len() != n {
return Err(format!(
"exchange has {} entries for {n} types",
self.exchange.len()
));
}
if !self.chemotaxis.is_empty() && self.chemotaxis.len() != n {
return Err(format!(
"chemotaxis has {} entries for {n} types",
self.chemotaxis.len()
));
}
for (index, row) in self.chemotaxis.iter().enumerate() {
if !row.is_empty() && row.len() != self.species.len() {
return Err(format!(
"type {index} states {} chemotactic sensitivities for {} species",
row.len(),
self.species.len()
));
}
}
for (index, spec) in self.exchange.iter().enumerate() {
for (what, list) in [("secretion", &spec.secretion), ("uptake", &spec.uptake)] {
if !list.is_empty() && list.len() != self.species.len() {
return Err(format!(
"type {index} states {} {what} rates for {} species",
list.len(),
self.species.len()
));
}
}
}
Ok(())
}
}
impl Default for Model {
/// A one-site medium-only lattice at a temperature that runs. It exists so
/// a caller states the terms it cares about and leaves the rest, which is
/// what keeps a new energy term from touching every call site.
fn default() -> Self {
Self {
species: Vec::new(),
exchange: Vec::new(),
chemotaxis: Vec::new(),
width: 1,
height: 1,
depth: 1,
contact: vec![0.0],
types: vec![CellType::default()],
temperature: 10.0,
neighbour_order: 2,
seed: 0,
}
}
}