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
use crate::{
err::Result,
fit::Attribute,
sim::{CVPolicy, GlobalProperties, IClamp, ModelType, Node, Probe, Simulation},
Map,
};
use anyhow::{anyhow, bail};
use serde::Serialize;
/// Resources to store in the output.
/// (src_gid, src_tag, tgt_tag, weight, delay)
type ConnectionData = (usize, usize, usize, f64, f64);
/// (location, variable, tag)
type ProbeData = (Option<String>, String, usize);
/// (location, mech, params, tag)
type SynapseData = (String, String, Map<String, f64>, usize);
/// (location, delay, duration, current, tag)
type IClampData = (String, f64, f64, f64, usize);
/// Metadata about cell,
/// mostly info discarded during generation
#[derive(Debug, Serialize)]
pub struct CellMetaData {
/// cell kind
pub kind: String,
/// population name
pub population: String,
/// cell type label
pub type_id: u64,
}
impl CellMetaData {
pub fn from(node: &Node) -> Self {
let kind = match node.node_type.model_type {
ModelType::Biophysical { .. } => String::from("biophys"),
ModelType::Single { .. } => String::from("single"),
ModelType::Point { .. } => String::from("point"),
ModelType::Virtual { .. } => String::from("virtual"),
};
Self {
population: node.pop.to_string(),
kind,
type_id: node.node_type.type_id,
}
}
}
#[derive(Debug, Serialize)]
pub struct CableGlobalProperties {
pub celsius: f64,
pub v_init: f64,
}
#[derive(Debug, Serialize)]
pub struct Bundle {
pub time: f64,
pub time_step: f64,
pub size: usize,
pub max_cv_length: Option<f64>,
/// gid to morphology and acc ids
/// this works as an index into the next two fields.
pub cell_bio_ids: Map<usize, (usize, usize)>,
/// acc index to name
pub morphology: Vec<String>,
/// morphology index to name
pub decoration: Vec<String>,
/// cell kinds, 0 = cable, 1 = lif, 2 = spike source, ...
pub cell_kind: Vec<u64>,
/// synapse data, cross-linked with incoming connections.
/// Location, Synapse, Parameters, Tag
/// May only be set iff kind==cable
pub synapses: Map<usize, Vec<SynapseData>>,
/// stimuli; May only be set iff kind==cable
/// location, delay, duration, amplitude, tag
pub current_clamps: Map<usize, Vec<IClampData>>,
/// List of data exporters
/// location, variable, tag. NOTE _could_ make variable an u64?!
pub probes: Map<usize, Vec<ProbeData>>,
/// Incoming connections as (src_gid, src_tag, tgt_tag, weight, delay)
pub incoming_connections: Map<usize, Vec<ConnectionData>>,
/// Spiking threshold
pub spike_threshold: f64,
/// sparse map of gids to LIF cell descrption. Valid iff kind(gid) == LIF
pub gid_to_lif: Map<usize, Map<String, f64>>,
/// sparse map of gids to virtual cell spike trains. Valid iff kind(gid) == Virtual
/// Will generate SpikeSource cells in Arbor
pub gid_to_vrt: Map<usize, Vec<f64>>,
/// dense map of gids to metadata
pub metadata: Vec<CellMetaData>,
/// cable cell global settings
pub cable_cell_globals: Option<CableGlobalProperties>,
/// cell counts by kind
pub count_by_kind: [usize; 3],
}
const KIND_CABLE: u64 = 0;
const KIND_LIF: u64 = 1;
const KIND_SOURCE: u64 = 2;
const SYNAPSES: &[(&str, &str)] = &[("exp2syn", "Exp2Syn")];
fn fudge_synapse_dynamics(old: &str) -> String {
if let Some(rep) = SYNAPSES.iter().find(|p| p.0 == old) {
rep.1
} else {
old
}
.to_string()
}
impl Bundle {
pub fn new(sim: &Simulation) -> Result<Self> {
// Reverse lookup tables, used internally for uniqueness and index generation.
let mut acc_to_cid = Map::new();
let mut mrf_to_mid = Map::new();
// Look up tables to write out
let mut gid_to_meta = Vec::new();
let mut cell_bio_ids = Map::new();
let mut morphology = Vec::new();
let mut decoration = Vec::new();
let mut cell_kind = Vec::new();
let mut incoming_connections = Map::new();
let mut synapses = Map::new();
let mut current_clamps = Map::new();
let mut probes = Map::new();
let mut gid_to_lif = Map::new();
let mut gid_to_vrt = Map::new();
let mut count_by_kind = [0; 3];
for gid in 0..sim.size {
let node = sim.reify_node(gid)?;
gid_to_meta.push(CellMetaData::from(&node));
if !node.incoming_edges.is_empty() {
if matches!(node.node_type.model_type, ModelType::Biophysical { .. }) {
let mut inc = Vec::new();
let mut syn = Vec::new();
for (ix, edge) in node.incoming_edges.iter().enumerate() {
inc.push((
edge.src_gid as usize,
0, // in our SONATA model, there is _one_ source on each cell.
ix,
edge.weight,
edge.delay,
));
let mech = edge.mech.as_ref().ok_or(anyhow!("Edge has no mechanism"))?;
let mech = fudge_synapse_dynamics(mech);
let loc = format!(
"(on-components {} (segment {}))",
edge.target.1, edge.target.0
);
syn.push((loc, mech, edge.dynamics.clone(), ix));
}
incoming_connections.insert(gid, inc);
synapses.insert(gid, syn);
} else {
let inc = node
.incoming_edges
.iter()
.map(|e| {
(
e.src_gid as usize,
0, // in our SONATA model, there is _one_ source on each cell.
0,
e.weight,
e.delay,
)
})
.collect::<Vec<_>>();
incoming_connections.insert(gid, inc);
};
}
match &node.node_type.model_type {
ModelType::Biophysical {
model_template,
attributes,
} => {
cell_kind.push(KIND_CABLE);
count_by_kind[KIND_CABLE as usize] += 1;
match model_template.as_ref() {
"ctdb:Biophys1.hoc" => {
let mid = if let Some(Attribute::String(mrf)) =
attributes.get("morphology")
{
if !mrf_to_mid.contains_key(mrf) {
let mid = morphology.len();
morphology.push(mrf.to_string());
mrf_to_mid.insert(mrf.to_string(), mid);
}
mrf_to_mid[mrf]
} else {
bail!("GID {gid} is a biophysical cell, but has no morphology.");
};
let cid = if let Some(Attribute::String(fit)) =
attributes.get("dynamics_params")
{
let acc = fit;
if !acc_to_cid.contains_key(acc) {
let cid = decoration.len();
decoration.push(acc.to_string());
acc_to_cid.insert(acc.to_string(), cid);
}
acc_to_cid[acc]
} else {
bail!(
"GID {gid} is a biophysical cell, but has no dynamics_params."
);
};
cell_bio_ids.insert(gid, (mid, cid));
}
acc if acc.starts_with("nml:") => {
let mid = if let Some(Attribute::String(mrf)) =
attributes.get("morphology")
{
if !mrf_to_mid.contains_key(mrf) {
let mid = morphology.len();
morphology.push(mrf.to_string());
mrf_to_mid.insert(mrf.to_string(), mid);
}
mrf_to_mid[mrf]
} else {
bail!("GID {gid} is a biophysical cell, but has no morphology.");
};
let acc = acc.strip_prefix("nml:").unwrap();
if !acc_to_cid.contains_key(acc) {
let cid = decoration.len();
decoration.push(acc.to_string());
acc_to_cid.insert(acc.to_string(), cid);
};
let cid = acc_to_cid[acc];
cell_bio_ids.insert(gid, (mid, cid));
}
t => bail!("Unknown model template <{t}> for gid {gid}"),
}
}
ModelType::Virtual { .. } => {
// The fields are largely irrelevant here
let data: &mut Vec<f64> = gid_to_vrt.entry(gid).or_default();
if let Some(group) = sim.virtual_spikes.get(&node.pop) {
if let Some(ts) = group.get(&node.node_id) {
data.append(&mut ts.clone());
}
}
cell_kind.push(KIND_SOURCE);
count_by_kind[KIND_SOURCE as usize] += 1;
}
ModelType::Point { model_template, .. } => {
cell_kind.push(KIND_LIF);
count_by_kind[KIND_LIF as usize] += 1;
match model_template.as_ref() {
"nrn:IntFire1" => {
// Taken from nrn/IntFire1.mod and adapted to Arbor.
let mut params = Map::from([
("cm".to_string(), 1.0),
("U_neutral".to_string(), 0.0),
("U_reset".to_string(), 0.0),
("U_th".to_string(), 1.0), // NOTE IntFire1 do be weird.
("U_0".to_string(), 0.0),
("t_ref".to_string(), 5.0),
("tau".to_string(), 10.0),
]);
for (k, v) in node.dynamics.iter() {
match k.as_ref() {
"tau" =>
params.insert("tau".to_string(), *v),
"refrac" =>
params.insert("t_ref".to_string(), *v),
_ => bail!("Unknown parameter <{k}> for template IntFire1 at gid {gid}")
};
}
gid_to_lif.insert(gid, params);
}
t => bail!("Unknown model template <{t}> for gid {gid}"),
}
}
mt => bail!("Cannot write ModelType {mt:?}"),
}
}
for (gid, ics) in &sim.iclamps {
let mut stim = Vec::new();
for IClamp {
amplitude_nA,
delay_ms,
duration_ms,
tag,
location,
} in ics
{
stim.push((
location.clone(),
*delay_ms,
*duration_ms,
*amplitude_nA,
*tag,
));
}
current_clamps.insert(*gid as usize, stim);
}
for (gid, sim_probes) in &sim.reports {
let mut prbs = Vec::new();
for (ix, probe) in sim_probes.iter().enumerate() {
match probe {
Probe::CableVoltage(ls) => {
for l in ls {
prbs.push((Some(l.clone()), "voltage".into(), ix));
}
}
Probe::Lif => {
prbs.push((None, "voltage".into(), ix));
}
Probe::CableIntConc(ion, ls) => {
for l in ls {
prbs.push((Some(l.clone()), ion.clone(), ix));
}
}
Probe::CableState(var, ls) => {
for l in ls {
prbs.push((Some(l.clone()), var.clone(), ix));
}
}
Probe::CableExtConc(ion, ls) => {
for l in ls {
prbs.push((Some(l.clone()), ion.clone(), ix));
}
}
}
}
probes.insert(*gid as usize, prbs);
}
let max_cv_length = match &sim.cv_policy {
CVPolicy::Default => None,
CVPolicy::MaxExtent(l) => Some(*l),
};
let cable_cell_globals =
if let Some(GlobalProperties { celsius, v_init }) = sim.global_properties {
Some(CableGlobalProperties { celsius, v_init })
} else {
None
};
// Sort all spike sources. Just to be sure...
gid_to_vrt.values_mut().for_each(|ts| {
ts.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))
});
Ok(Bundle {
time: sim.tfinal,
time_step: sim.dt,
max_cv_length,
size: sim.size,
cell_bio_ids,
morphology,
decoration,
synapses,
cable_cell_globals,
probes,
incoming_connections,
cell_kind,
count_by_kind,
current_clamps,
spike_threshold: sim.spike_threshold,
gid_to_lif,
gid_to_vrt,
metadata: gid_to_meta,
})
}
}