forked from jlubo/arbor_network_consolidation
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy patharborNetworkConsolidation.py
executable file
·1733 lines (1415 loc) · 87.3 KB
/
arborNetworkConsolidation.py
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
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/bin/python3
# Arbor simulation of memory consolidation in recurrent spiking neural networks, consisting of leaky integrate-and-fire neurons that are
# connected via current-based, plastic synapses (early phase based on calcium dynamics and late phase based on synaptic tagging and capture).
# Intended to reproduce the results of Luboeinski and Tetzlaff, Commun. Biol., 2021 (https://doi.org/10.1038/s42003-021-01778-y). For further
# details, see https://doi.org/10.48550/arXiv.2411.16445.
# Copyright 2021-2024 Jannik Luboeinski
# License: Apache-2.0 (http://www.apache.org/licenses/LICENSE-2.0)
# Contact: mail[at]jlubo.net
import arbor
from arbor import units as U
import numpy as np
import gc
import json
import argparse
import os
import subprocess
import traceback
import copy
from outputUtilities import getTimestamp, getFormattedTime, getTimeDiff, \
setDataPathPrefix, getDataPath, \
openLog, writeLog, writeAddLog, flushLog, closeLog, \
getPresynapticId
setDataPathPrefix("data_")
###############################################################################
# completeProt
# Takes a, possibly incomplete, stimulation protocol dictionary and returns a complete
# protocol after adding keys that are missing. This allows to leave out unnecessary keys
# when defining protocols in the JSON config file.
# - prot: stimulation protocol, possibly incomplete
# - return: complete stimulation protocol
def completeProt(prot):
if prot.get("scheme") is None:
prot["scheme"] = ""
if prot.get("time_start") is None:
prot["time_start"] = 0
if prot.get("duration") is None:
prot["duration"] = 0
if prot.get("freq") is None:
prot["freq"] = 0
if prot.get("N_stim") is None:
prot["N_stim"] = 0
if prot.get("I_0") is None:
prot["I_0"] = 0
if prot.get("sigma_WN") is None:
prot["sigma_WN"] = 0
if prot.get("explicit_input") is None \
or prot["explicit_input"].get("receivers") is None \
or prot["explicit_input"].get("stim_times") is None:
prot["explicit_input"] = { "receivers": [ ], "stim_times": [ ] }
return prot
###############################################################################
# protFormatted
# Takes a complete stimulation protocol dictionary and creates a formatted string for
# human-readable output.
# - prot: stimulation protocol
# - return: formatted string with information about the protocol
def protFormatted(prot):
if not prot['scheme']:
fmstr = "none"
elif prot['scheme'] == "EXPLICIT":
fmstr = f"{prot['scheme']}, pulses at {prot['explicit_input']['stim_times']} ms, to neurons {prot['explicit_input']['receivers']}"
elif prot['scheme'] == "STET":
fmstr = f"{prot['scheme']} (Poisson) to neurons {prot['explicit_input']['receivers']}, starting at {prot['time_start']} s"
elif prot['I_0'] != 0:
fmstr = f"{prot['scheme']} (OU), I_0 = {prot['I_0']} nA, sigma_WN = {prot['sigma_WN']} nA s^1/2"
elif prot['duration'] != 0:
fmstr = f"{prot['scheme']} (OU), {prot['N_stim']} input neurons at {prot['freq']} Hz, starting at {prot['time_start']} s, lasting for {prot['duration']} s"
else:
fmstr = f"{prot['scheme']} (OU), {prot['N_stim']} input neurons at {prot['freq']} Hz, starting at {prot['time_start']} s"
return fmstr
###############################################################################
# getInputProtocol
# Defines specific input/stimulation protocols. Returns start and end time. If 'label' is provided,
# also returns event generator instances (may either implement an explicit schedule, or a protocol for
# an 'ou_input' mechanism with the following behavior: if the value of the 'weight' parameter is 1,
# is switched on, else if it is -1, stimulation is switched off).
# - protocol: protocol that is being used for stimulation
# - runtime: full runtime of the simulation (with unit)
# - dt : duration of one timestep (with unit)
# - explicit_input_times: list of times (in ms) at which explicit input pulses shall be given
# - label [optional]: label of the mechanism that receives the the stimulation
# - return: a list of 'arbor.event_generator' instances, the start time of the stimulus, and the end time of the stimulus in ms
def getInputProtocol(protocol, runtime, dt, explicit_input_times, label = ""):
prot_name = protocol['scheme'] # name of the protocol (defining its structure)
start_time = protocol['time_start']*U.s # time at which the stimulus starts
if prot_name == "RECT":
duration = protocol['duration']*U.s # duration of the stimulus
end_time = start_time + duration + dt # time at which the stimulus ends
if label:
# create regular schedules to implement a stimulation pulse that lasts for 'duration'
stim_on = arbor.event_generator(
label,
1,
arbor.regular_schedule(start_time, dt, start_time + dt))
stim_off = arbor.event_generator(
label,
-1,
arbor.regular_schedule(start_time + duration, dt, start_time + duration + dt))
return [stim_on, stim_off], start_time, end_time
return [], start_time, end_time
elif prot_name == "ONEPULSE":
end_time = start_time + 100*U.ms + dt # time at which the stimulus ends
if label:
# create regular schedules to implement a stimulation pulse that lasts for 0.1 s
stim_on = arbor.event_generator(
label,
1,
arbor.regular_schedule(start_time, dt, start_time + dt))
stim_off = arbor.event_generator(
label,
-1,
arbor.regular_schedule(start_time + 100*U.ms, dt, end_time))
return [stim_on, stim_off], start_time, end_time
return [], start_time, end_time
elif prot_name == "TRIPLET":
end_time = start_time + 1100*U.ms + dt # time at which the stimulus ends in ms
if label:
# create regular schedules to implement pulses that last for 0.1 s each
stim1_on = arbor.event_generator(
label,
1,
arbor.regular_schedule(start_time, dt, start_time + dt))
stim1_off = arbor.event_generator(
label,
-1,
arbor.regular_schedule(start_time + 100*U.ms, dt, start_time + 100*U.ms + dt))
stim2_on = arbor.event_generator(
label,
1,
arbor.regular_schedule(start_time + 500*U.ms, dt, start_time + 500*U.ms + dt))
stim2_off = arbor.event_generator(
label,
-1,
arbor.regular_schedule(start_time + 600*U.ms, dt, start_time + 600*U.ms + dt))
stim3_on = arbor.event_generator(
label,
1,
arbor.regular_schedule(start_time + 1000*U.ms, dt, start_time + 1000*U.ms + dt))
stim3_off = arbor.event_generator(
label,
-1,
arbor.regular_schedule(start_time + 1100*U.ms, dt, end_time))
return [stim1_on, stim1_off, stim2_on, stim2_off, stim3_on, stim3_off], start_time, end_time
return [], start_time, end_time
elif prot_name == "STET":
duration = protocol['duration']*U.s # duration of the stimulus (for "standard" STET: 1200 ms)
last_pulse_start_time = start_time + duration
t_start = np.linspace(start_time, last_pulse_start_time, num=3, endpoint=True) # start times of the pulses (in ms)
t_end = np.linspace(start_time + 1000*U.ms, last_pulse_start_time + 1000*U.ms, num=3, endpoint=True) # end times of the pulses (in ms)
end_time = t_end[-1]
freq = protocol['freq']*U.Hz # average spike frequency in Hz
seed = int(datetime.now().timestamp() * 1e6)
# average number of spikes (random number drawn for every timestep, then filtered with probability):
stimulus_times_exc = np.array([])
rng = np.random.default_rng(seed)
num_timesteps = np.int_(np.round_((t_end[0]-t_start[0])/dt))
for i in range(len(t_start)):
spike_mask = rng.random(size=num_timesteps) < freq*dt/1000.
timestep_values = np.linspace(t_start[i], t_end[i], num=num_timesteps, endpoint=False)
spikes = timestep_values[spike_mask]
stimulus_times_exc = np.concatenate([stimulus_times_exc, spikes])
if label:
stim_explicit = arbor.event_generator(
label,
0.,
arbor.explicit_schedule(stimulus_times_exc))
return [stim_explicit], start_time, end_time
return [], start_time, end_time
elif prot_name == "FULL":
end_time = runtime # time at which the stimulus ends
if label:
# create a regular schedule that lasts for the full runtime
stim_on = arbor.event_generator(
label,
1,
arbor.regular_schedule(start_time, dt, start_time + dt))
return [stim_on], start_time, end_time
return [], start_time, end_time
elif prot_name == "EXPLICIT" and explicit_input_times.any():
# depends on input pulses explicitly defined by 'explicit_input_times'
start_time = np.min(explicit_input_times)*U.ms # time at which the stimulation starts
end_time = np.max(explicit_input_times)*U.ms # time at which the stimulation ends
if label:
# create a schedule with explicitly defined pulses
stim_explicit = arbor.event_generator(
label,
0.,
arbor.explicit_schedule(explicit_input_times*U.ms))
return [stim_explicit], start_time, end_time
return [], start_time, end_time
return [], np.nan*U.ms, np.nan*U.ms
#####################################
# NetworkRecipe
# Implementation of Arbor simulation recipe
class NetworkRecipe(arbor.recipe):
# constructor
# Sets up the recipe object; sets parameters from config dictionary (deep-copies the non-immutable structures to enable their modification)
# - config: dictionary containing configuration data
# - adap_dt: duration of one adapted timestep
# - plasticity_mechanism: the plasticity mechanism to be used
# - store_latest_state: specifies if latest synaptic state shall be probed and stored
def __init__(self, config, adap_dt, plasticity_mechanism, store_latest_state):
# The base C++ class constructor must be called first, to ensure that
# all memory in the C++ class is initialized correctly. (see https://github.com/tetzlab/FIPPA/blob/main/STDP/arbor_lif_stdp.py)
arbor.recipe.__init__(self)
self.s_desc = config['simulation']['short_description'] # short description of the simulation
self.N_exc = int(config["populations"]["N_exc"]) # number of neurons in the excitatory population
self.N_inh = int(config["populations"]["N_inh"]) # number of neurons in the inhibitory population
self.N_CA = int(config["populations"]["N_CA"]) # number of neurons in the cell assembly
self.N_tot = self.N_exc + self.N_inh # total number of neurons (excitatory and inhibitory)
self.p_c = config["populations"]["p_c"] # probability of connection
self.p_r = config["populations"]["p_r"] # probability of recall stimulation
self.props = arbor.neuron_cable_properties() # initialize the cell properties to match Neuron's defaults
# (cf. https://docs.arbor-sim.org/en/v0.5.2/tutorial/single_cell_recipe.html)
cat = arbor.load_catalogue("./custom-catalogue.so") # load the catalogue of custom mechanisms
cat.extend(arbor.default_catalogue(), "") # add the default catalogue
cat.extend(arbor.stochastic_catalogue(), "") # add the stochastic catalogue
self.props.catalogue = cat
self.plasticity_mechanism = plasticity_mechanism #+ "/sps_=sps" # set the plasticity mechanism and the particle for diffusion
self.runtime = config["simulation"]["runtime"]*U.s # runtime of the simulation in biological time
self.std_dt = config["simulation"]["dt"]*U.ms # duration of one standard timestep (for rich computation)
self.adap_dt = adap_dt # duration of one adapted timestep
self.neuron_config = config["neuron"]
self.syn_config = config["synapses"]
self.syn_plasticity_config = config["synapses"]["syn_exc_calcium_plasticity"]
self.h_0 = self.syn_plasticity_config["h_0"]*U.mV
self.D = 1*U.m2/U.s # diffusivity
self.w_ei = config["populations"]["w_ei"]*U.mV
self.w_ie = config["populations"]["w_ie"]*U.mV
self.w_ii = config["populations"]["w_ii"]*U.mV
self.learn_prot = completeProt(config["simulation"]["learn_protocol"]) # protocol for learning stimulation as a dictionary with the keys "time_start" (starting time), "scheme" (scheme of pulses), and "freq" (stimulation frequency)
self.recall_prot = completeProt(config["simulation"]["recall_protocol"]) # protocol for recall stimulation as a dictionary with the keys "time_start" (starting time), "scheme" (scheme of pulses), and "freq" (stimulation frequency)
self.bg_prot = completeProt(config["simulation"]["bg_protocol"]) # protocol for background input as a dictionary with the keys "time_start" (starting time), "scheme" (scheme of pulses), "I_0" (mean), and "sigma_WN" (standard deviation)
self.sample_gid_list = config['simulation']['sample_gid_list'] # list of the neurons that are to be probed (given by number/gid)
self.sample_pre_list = config['simulation']['sample_pre_list'] # list of the presynaptic neurons for probing synapses (one value for each value in `sample_gid`; -1: no synapse probing)
if config['populations']['conn_file']: # if a connections file is specified, load the connectivity matrix from that file
self.conn_matrix = np.loadtxt(config['populations']['conn_file']).T
else: # there is no pre-defined connectivity matrix -> generate one
rng = np.random.default_rng() # random number generator
self.conn_matrix = rng.random((self.N_tot, self.N_tot)) <= self.p_c # two-dim. array of booleans indicating the existence of any incoming connection
self.conn_matrix[np.identity(self.N_tot, dtype=bool)] = 0 # remove self-couplings
np.savetxt(getDataPath(self.s_desc, "connections.txt"), self.conn_matrix.T, fmt="%.0f") # save the connectivity matrix
## set diffusing particles
self.props.set_ion("sps_", valence=1, int_con=0*U.mM, ext_con=0*U.mM, rev_pot=0*U.mV) #diff=self.D) # signal triggering protein synthesis
self.props.set_ion("pc_", valence=1, int_con=0*U.mM, ext_con=0*U.mM, rev_pot=0*U.mV) # common pool of plasticity-related proteins to diffuse
# name, charge, internal_concentration, external_concentration, reversal_potential, reversal_potential_method, diffusivity
# for testing
self.exc_neurons_counter = 0
self.inh_neurons_counter = 0
self.connections_counter = 0
# states to store/load
self.h = None
self.z = None
self.p = None
self.max_num_trace_probes = 0
self.max_num_latest_state_probes = None
self.store_latest_state = store_latest_state
# loadState
# Loads the relevant variables of the latest state of all neurons and synapses (from the end of the previous simulation phase)
# - prev_phase: number specifying the previous phase from which to load the data
def loadState(self, prev_phase):
state_path = getDataPath(self.s_desc, f"state_after_phase_{prev_phase}")
# loading the adjacency matrices for h, z, and p
self.h = np.loadtxt(os.path.join(state_path, "h_synapses.txt")).T
self.z = np.loadtxt(os.path.join(state_path, "z_synapses.txt")).T
self.p = np.loadtxt(os.path.join(state_path, "p_compartments.txt"))
# cell_kind
# Defines the kind of the neuron given by gid
# - gid: global identifier of the cell
# - return: type of the cell
def cell_kind(self, gid):
return arbor.cell_kind.cable # note: implementation of arbor.cell_kind.lif is not ready to use yet
# cell_description
# Defines the morphology, cell mechanism, etc. of the neuron given by gid
# - gid: global identifier of the cell
# - return: description of the cell
def cell_description(self, gid):
# cylinder morphology
tree = arbor.segment_tree()
radius_µm = self.neuron_config["radius"] # radius of cylinder in µm
height_µm = 2*radius_µm # height of cylinder in µm
tree.append(arbor.mnpos,
arbor.mpoint(-1/2 * height_µm, 0, 0, radius_µm),
arbor.mpoint(1/2 * height_µm, 0, 0, radius_µm),
tag=1)
labels = arbor.label_dict({"center": "(location 0 0.5)"})
area = 2 * np.pi * radius_µm*U.um * height_µm*U.um # surface area of the cylinder (excluding the circle-shaped ends, since Arbor does not consider current flux there)
area_cm2 = 2 * np.pi * radius_µm * height_µm * (1e-4)**2 # surface area of the cylinder in cm^2
volume = np.pi * radius_µm**2*U.um2 * height_µm*U.um # volume of the cylinder
self.volume = volume
#i_factor = (1/area) * U.nA/(U.mA/U.cm2) # current to current density conversion factor (nA to mA/cm^2; necessary for point neurons)
i_factor = (1e-9/1e-3) / area_cm2 # current to current density conversion factor (nA to mA/cm^2; necessary for point neurons)
C_mem = self.neuron_config["C_mem"]*U.F # membrane capacitance
c_mem = C_mem / area # specific membrane capacitance in F/m^2, computed from absolute capacitance of a point neuron
# initialize decor
decor = arbor.decor()
# set CV policy
cv_policy = arbor.cv_policy("(max-extent 100)")
#decor.discretization(cv_policy) # since v0.10.1-dev-76120d1, this has to be done in `cable_cell()` (see below)
# neuronal dynamics
decor.set_property(Vm=self.neuron_config["V_init"]*U.mV, cm=c_mem)
mech_neuron = arbor.mechanism(self.neuron_config["mechanism"])
R_leak = self.neuron_config["R_leak"]*U.MOhm
R_reset = self.neuron_config["R_reset"]*U.MOhm
tau_mem = R_leak * C_mem # membrane time constant
V_rev = self.neuron_config["V_rev"]*U.mV
V_reset = self.neuron_config["V_reset"]*U.mV
V_th = self.neuron_config["V_th"]*U.mV
t_ref = self.neuron_config["t_ref"]*U.ms
U_substance_amount = U.nano*U.nano*U.mol #1e-18*U.mol
theta_pro = self.syn_plasticity_config["theta_pro"]*U_substance_amount
p_max = self.syn_plasticity_config["p_max"]*U.mM
mech_neuron.set("R_leak", R_leak.value_as(U.MOhm))
mech_neuron.set("R_reset", R_reset.value_as(U.MOhm))
mech_neuron.set("I_0", 0) # set to zero (background input is applied via OU process ou_bg)
mech_neuron.set("i_factor", i_factor)
mech_neuron.set("V_rev", V_rev.value_as(U.mV))
mech_neuron.set("V_reset", V_reset.value_as(U.mV))
mech_neuron.set("V_th", V_th.value_as(U.mV))
mech_neuron.set("t_ref", t_ref.value_as(U.ms))
mech_neuron.set("theta_pro", theta_pro.value_as(U_substance_amount))
mech_neuron.set("p_max", p_max.value_as(U.mM))
# diffusion of particles/signals
decor.set_ion("sps_", int_con=0*U.mM, diff=self.D) # signal to trigger protein synthesis
decor.set_ion("pc_", int_con=0*U.mM, diff=self.D) # proteins
# excitatory (postsynaptic) neurons
if gid < self.N_exc:
# set neuronal state variables to values loaded from previous state of the network
if self.p is not None:
mech_neuron.set('pc_init', self.p[gid])
# parameter output
if gid == 0:
writeAddLog("area =", area, "µm^2")
writeAddLog("volume =", volume, "µm^3")
writeAddLog("i_factor =", i_factor, "(mA/cm^2)/nA")
writeAddLog("c_mem =", c_mem, "F/m^2")
writeAddLog("tau_mem =", tau_mem, "ms")
# initialize plastic exponential mechanism for excitatory
mech_expsyn_exc = arbor.mechanism(self.plasticity_mechanism)
# set standard parameters
mech_expsyn_exc.set('h_0', self.h_0.value_as(U.mV))
mech_expsyn_exc.set('theta_tag', self.syn_plasticity_config["theta_tag"])
#mech_expsyn_exc.set('area', area.value_as(U.um**2))
if not self.plasticity_mechanism[-3:] == "_ff":
mech_expsyn_exc.set('R_mem', R_leak.value_as(U.MOhm))
mech_expsyn_exc.set('tau_syn', self.syn_config["tau_syn"])
mech_expsyn_exc.set('Ca_pre', self.syn_plasticity_config["Ca_pre"])
mech_expsyn_exc.set('Ca_post', self.syn_plasticity_config["Ca_post"])
mech_expsyn_exc.set('theta_p', self.syn_plasticity_config["theta_p"])
mech_expsyn_exc.set('theta_d', self.syn_plasticity_config["theta_d"])
mech_expsyn_exc.set('sigma_pl', self.syn_plasticity_config["sigma_pl"])
exc_exc_connections = np.array(self.conn_matrix[gid][0:self.N_exc], dtype=bool) # array of booleans defining all incoming excitatory connections
exc_exc_pre_neurons = np.arange(self.N_exc)[exc_exc_connections] # array of excitatory presynaptic neurons given by their gid
# place exc.->exc. synapses; set synaptic state variables to values loaded from previous state of the network
for presyn_gid in exc_exc_pre_neurons:
if (self.h is not None and self.z is not None):
mech_expsyn_exc.set('h_init', self.h[gid][presyn_gid])
mech_expsyn_exc.set('z_init', self.z[gid][presyn_gid])
if self.N_exc <= 4: # only for small networks
writeAddLog(f"Setting loaded values for synapse {presyn_gid}->{gid}...\n" +
f" h = {self.h[gid][presyn_gid]}\n" +
f" z = {self.z[gid][presyn_gid]}")
# add synapse
decor.place('"center"', arbor.synapse(mech_expsyn_exc), f"syn_ee_calcium_plasticity_{presyn_gid}_{gid}") # place synapse at the center of the soma
# (because of point neuron approximation)
# initialize non-plastic exponential mechanism for inhibitory synapse
mech_expsyn_inh = arbor.mechanism('expsyn_curr')
mech_expsyn_inh.set('w', -self.w_ie.value_as(U.mV) * self.h_0.value_as(U.mV))
mech_expsyn_inh.set('R_mem', R_leak.value_as(U.MOhm))
mech_expsyn_inh.set('tau', self.syn_config["tau_syn"])
# place inh.->exc. synapses
inh_exc_connections = np.array(self.conn_matrix[gid][self.N_exc:self.N_tot], dtype=bool) # array of booleans defining all incoming inhibitory connections
inh_exc_pre_neurons = np.arange(self.N_exc, self.N_tot)[inh_exc_connections] # array of inhibitory presynaptic neurons given by their gid
for presyn_gid in inh_exc_pre_neurons:
decor.place('"center"', arbor.synapse(mech_expsyn_inh), f"syn_ie_{presyn_gid}_{gid}") # place synapse at the center of the soma (because: point neuron)
self.exc_neurons_counter += 1 # for testing
# inhibitory (postsynaptic) neurons
else:
# initialize non-plastic exponential mechanism for excitatory synapse
mech_expsyn_exc = arbor.mechanism('expsyn_curr')
mech_expsyn_exc.set('w', self.w_ei.value_as(U.mV) * self.h_0.value_as(U.mV))
mech_expsyn_exc.set('R_mem', R_leak.value_as(U.MOhm))
mech_expsyn_exc.set('tau', self.syn_config["tau_syn"])
# place exc.->inh. synapses
exc_inh_connections = np.array(self.conn_matrix[gid][0:self.N_exc], dtype=bool) # array of booleans defining all incoming excitatory connections
exc_inh_pre_neurons = np.arange(self.N_exc)[exc_inh_connections] # array of inhibitory presynaptic neurons given by their gid
for presyn_gid in exc_inh_pre_neurons:
decor.place('"center"', arbor.synapse(mech_expsyn_exc), f"syn_ei_{presyn_gid}_{gid}") # place synapse at the center of the soma (because: point neuron)
# initialize non-plastic exponential mechanism for inhibitory synapse
mech_expsyn_inh = arbor.mechanism('expsyn_curr')
mech_expsyn_inh.set('w', -self.w_ii.value_as(U.mV) * self.h_0.value_as(U.mV))
mech_expsyn_inh.set('R_mem', R_leak.value_as(U.MOhm))
mech_expsyn_inh.set('tau', self.syn_config["tau_syn"])
# place inh.->inh. synapses
inh_inh_connections = np.array(self.conn_matrix[gid][self.N_exc:self.N_tot], dtype=bool) # array of booleans defining all incoming inhibitory connections
inh_inh_pre_neurons = np.arange(self.N_exc, self.N_tot)[inh_inh_connections] # array of inhibitory presynaptic neurons given by their gid
for presyn_gid in inh_inh_pre_neurons:
decor.place('"center"', arbor.synapse(mech_expsyn_inh), f"syn_ii_{presyn_gid}_{gid}") # place synapse at the center of the soma (because: point neuron)
self.inh_neurons_counter += 1 # for testing
# learning stimulation to all neurons in the assembly core ('as' and 'ans' subpopulations); input current described by Ornstein-Uhlenbeck process
# accounting for a population of neurons
if gid < self.N_CA:
mech_ou_learn_stim = arbor.mechanism('ou_input')
ou_mu = self.learn_prot['N_stim'] * self.learn_prot['freq']*U.Hz * 1*U.s * self.h_0 / R_leak
#ou_sigma = np.sqrt(self.learn_prot['N_stim'] * self.learn_prot['freq']*U.Hz / (2 * self.syn_config['tau_syn']*U.ms)) * 1*U.s * self.h_0/R_leak
ou_sigma = np.sqrt(1000.0 * self.learn_prot['N_stim'] * self.learn_prot['freq'] / (2 * self.syn_config['tau_syn'])) * self.h_0.value_as(U.mV)/R_leak.value_as(U.MOhm) * U.nA
mech_ou_learn_stim.set('mu', ou_mu.value_as(U.nA))
#mech_ou_learn_stim.set('sigma', ou_sigma.value_as(U.nA / np.sqrt(U.s))) # fractional units not supported!
mech_ou_learn_stim.set('sigma', ou_sigma.value_as(U.nA)) # technically incorrect unit (see line above)
mech_ou_learn_stim.set('tau', self.syn_config["tau_syn"])
decor.place('"center"', arbor.synapse(mech_ou_learn_stim), "ou_learn_stim")
# recall stimulation to some neurons in the assembly core ('as' subpopulation); input current described by Ornstein-Uhlenbeck process
# accounting for a population of neurons
if gid < self.p_r*self.N_CA:
mech_ou_recall_stim = arbor.mechanism('ou_input')
ou_mu = self.recall_prot['N_stim'] * self.recall_prot['freq']*U.Hz * 1*U.s * self.h_0 / R_leak
ou_sigma = np.sqrt(1000.0 * self.recall_prot['N_stim'] * self.recall_prot['freq'] / (2 * self.syn_config['tau_syn'])) * self.h_0.value_as(U.mV)/R_leak.value_as(U.MOhm) * U.nA
mech_ou_recall_stim.set('mu', ou_mu.value_as(U.nA))
#mech_ou_recall_stim.set('sigma', ou_sigma.value_as(U.nA / np.sqrt(U.s))) # fractional units not supported!
mech_ou_recall_stim.set('sigma', ou_sigma.value_as(U.nA)) # technically incorrect unit (see line above)
mech_ou_recall_stim.set('tau', self.syn_config["tau_syn"])
decor.place('"center"', arbor.synapse(mech_ou_recall_stim), "ou_recall_stim")
# background input current to all neurons; input current described by Ornstein-Uhlenbeck process
# with given mean and standard deviation
mech_ou_bg = arbor.mechanism('ou_input')
ou_mu = self.bg_prot["I_0"]*U.nA
ou_sigma = self.bg_prot['sigma_WN'] * np.sqrt(1000. / (2 * self.syn_config['tau_syn'])) * U.nA
mech_ou_bg.set('mu', ou_mu.value_as(U.nA))
#mech_ou_bg.set('sigma', ou_sigma.value_as(U.nA / np.sqrt(U.s))) # fractional units not supported!
mech_ou_bg.set('sigma', ou_sigma.value_as(U.nA)) # technically incorrect unit (see line above)
mech_ou_bg.set('tau', self.syn_config["tau_syn"])
decor.place('"center"', arbor.synapse(mech_ou_bg), "ou_bg")
# additional excitatory delta synapse for explicit input
mech_deltasyn_exc = arbor.mechanism('deltasyn')
mech_deltasyn_exc.set('g_spike', 1000*(V_th.value_as(U.mV)-V_reset.value_as(U.mV))*np.exp(self.std_dt.value_as(U.ms)/tau_mem.value_as(U.ms))) # choose sufficently large increase in conductance
decor.place('"center"', arbor.synapse(mech_deltasyn_exc), "explicit_input")
# place spike detector
decor.place('"center"', arbor.threshold_detector(V_th), "spike_detector")
# paint neuron mechanism
decor.paint('(all)', arbor.density(mech_neuron))
return arbor.cable_cell(tree, decor, labels, discretization=cv_policy)
# connections_on
# Defines the list of incoming synaptic connections to the neuron given by gid
# - gid: global identifier of the cell
# - return: connections to the given neuron
def connections_on(self, gid):
connections_list = []
rr = arbor.selection_policy.round_robin
rr_halt = arbor.selection_policy.round_robin_halt
connections = self.conn_matrix[gid]
assert connections[gid] == 0 # check that there are no self-couplings
self.connections_counter += np.sum(connections)
exc_connections = np.array(connections*np.concatenate((np.ones(self.N_exc, dtype=np.int8), np.zeros(self.N_inh, dtype=np.int8)), axis=None), dtype=bool) # array of booleans indicating all incoming excitatory connections
inh_connections = np.array(connections*np.concatenate((np.zeros(self.N_exc, dtype=np.int8), np.ones(self.N_inh, dtype=np.int8)), axis=None), dtype=bool) # array of booleans indicating all incoming inhibitory connections
assert not np.any(np.logical_xor(np.logical_or(exc_connections, inh_connections), connections)) # test if 'exc_connections' and 'inh_connections' together yield 'connections' again
exc_pre_neurons = np.arange(self.N_tot)[exc_connections] # array of excitatory presynaptic neurons indicated by their gid
inh_pre_neurons = np.arange(self.N_tot)[inh_connections] # array of inhibitory presynaptic neurons indicated by their gid
assert np.logical_and(np.all(exc_pre_neurons >= 0), np.all(exc_pre_neurons < self.N_exc)) # test if the excitatory neuron numbers are in the correct range
assert np.logical_and(np.all(inh_pre_neurons >= self.N_exc), np.all(inh_pre_neurons < self.N_tot)) # test if the inhibitory neuron numbers are in the correct range
# delay constants
#d0 = self.syn_config["t_ax_delay"] # delay time of the postsynaptic potential in ms
d0 = max(self.adap_dt.value_as(U.ms), self.syn_config["t_ax_delay"])*U.ms # delay time of the postsynaptic potential in ms
#d1 = self.syn_plasticity_config["t_Ca_delay"] # delay time of the calcium increase in ms (only for plastic synapses)
d1 = max(self.adap_dt.value_as(U.ms), self.syn_plasticity_config["t_Ca_delay"])*U.ms # delay time of the calcium increase in ms (only for plastic synapses)
# excitatory (postsynaptic) neurons
if gid < self.N_exc:
# incoming excitatory synapses
for src in exc_pre_neurons:
syn_target = f"syn_ee_calcium_plasticity_{src}_{gid}"
connections_list.append(arbor.connection((src, "spike_detector"), (syn_target, rr_halt), 1, d0)) # for postsynaptic potentials
connections_list.append(arbor.connection((src, "spike_detector"), (syn_target, rr), -1, d1)) # for plasticity-related calcium dynamics
# incoming inhibitory synapses
for src in inh_pre_neurons:
connections_list.append(arbor.connection((src,"spike_detector"), (f"syn_ie_{src}_{gid}", rr), 1, d0))
# inhibitory (postsynaptic) neurons
else:
# incoming excitatory synapses
for src in exc_pre_neurons:
connections_list.append(arbor.connection((src,"spike_detector"), (f"syn_ei_{src}_{gid}", rr), 1, d0))
# incoming inhibitory synapses
for src in inh_pre_neurons:
connections_list.append(arbor.connection((src,"spike_detector"), (f"syn_ii_{src}_{gid}", rr), 1, d0))
if self.N_exc <= 4: # only for small networks
writeAddLog(f"Setting connections for gid = {gid}...")
for conn in connections_list:
writeAddLog(f" {conn}")
return connections_list
# event_generators
# Event generators for input to synapses
# - gid: global identifier of the cell
# - return: events generated from Arbor schedule
def event_generators(self, gid):
inputs = []
# background input current to all neurons
stim, _, _ = getInputProtocol(self.bg_prot, self.runtime, self.std_dt,
np.array(self.bg_prot["explicit_input"]["stim_times"]), "ou_bg")
inputs.extend(stim)
# explicitly specified pulses for learning stimulation to defined receiving neurons
if (self.learn_prot["scheme"] == "EXPLICIT" and self.learn_prot["explicit_input"]["receivers"]) \
or self.learn_prot["scheme"] == "STET":
if gid in self.learn_prot["explicit_input"]["receivers"]:
stim, _, _ = getInputProtocol(self.learn_prot, self.runtime, self.std_dt,
np.array(self.learn_prot["explicit_input"]["stim_times"]), "explicit_input")
inputs.extend(stim)
# Ornstein-Uhlenbeck learning stimulation to cell assembly neurons ('as' and 'ans' subpopulations)
elif gid < self.N_CA:
stim, _, _ = getInputProtocol(self.learn_prot, self.runtime, self.std_dt,
np.array([]), "ou_learn_stim")
inputs.extend(stim)
# explicitly specified pulses for recall stimulation to defined receiving neurons
if self.recall_prot["scheme"] == "EXPLICIT" and self.recall_prot["explicit_input"]["receivers"] \
or self.recall_prot["scheme"] == "STET":
if gid in self.recall_prot["explicit_input"]["receivers"]:
stim, _, _ = getInputProtocol(self.recall_prot, self.runtime, self.std_dt,
np.array(self.recall_prot["explicit_input"]["stim_times"]), "explicit_input")
inputs.extend(stim)
# Ornstein-Uhlenbeck recall stimulation to cell assembly neurons ('as' subpopulation)
elif gid < self.p_r*self.N_CA:
stim, _, _ = getInputProtocol(self.recall_prot, self.runtime, self.std_dt,
np.array([]), "ou_recall_stim")
inputs.extend(stim)
return inputs
# global_properties
# Sets properties that will be applied to all neurons of the specified kind
# - gid: global identifier of the cell
# - return: the cell properties
def global_properties(self, kind):
assert kind == arbor.cell_kind.cable # assert that all neurons are technically cable cells
return self.props
# num_cells
# - return: the total number of cells in the network
def num_cells(self):
return self.N_tot
# set_max_num_trace_probes
# Sets the maximal number of probes to measure the traces of a particular neuron or synapse
def set_max_num_trace_probes(self, num):
if num > self.max_num_trace_probes:
self.max_num_trace_probes = num
# get_max_num_trace_probes
# - return: the number of probes to measure the traces of a particular neuron or synapse
def get_max_num_trace_probes(self):
return self.max_num_trace_probes
# set_max_num_latest_state_probes
# Sets the maximal number of probes to measure the latest synaptic state of a neuron
def set_max_num_latest_state_probes(self, num):
if self.max_num_latest_state_probes is None \
or num > self.max_num_latest_state_probes:
self.max_num_latest_state_probes = num
# get_max_num_latest_state_probes
# - return: the number of probes to measure the latest synaptic state of a neuron
def get_max_num_latest_state_probes(self):
if self.max_num_latest_state_probes is None:
return 0
return self.max_num_latest_state_probes
# probes
# Sets the probes to measure neuronal and synaptic state
# - gid: global identifier of the cell
# - return: the probes on the given cell
def probes(self, gid):
if self.N_exc <= 4: # only for small networks
writeAddLog(f"Setting probes for gid = {gid}...")
# loop over all potential excitatory presynaptic neurons to set probes for latest state
latest_state_probes = []
if self.store_latest_state and gid < self.N_exc:
latest_state_probes.extend([arbor.cable_probe_density_state('"center"', self.neuron_config["mechanism"], "pc", "tag_pc_latest")])
for presyn_gid in range(self.N_exc):
if self.conn_matrix[gid][presyn_gid]:
syn_target = f"syn_ee_calcium_plasticity_{presyn_gid}_{gid}"
latest_state_probes.extend([arbor.cable_probe_point_state(syn_target, self.plasticity_mechanism, "h", f"tag_h_latest_{presyn_gid}"),
arbor.cable_probe_point_state(syn_target, self.plasticity_mechanism, "z", f"tag_z_latest_{presyn_gid}")])
#arbor.cable_probe_point_state(syn_target, self.plasticity_mechanism, "pc", f"tag_pc_{presyn_gid}")])
self.set_max_num_latest_state_probes(len(latest_state_probes))
# loop over all synapses whose traces are to be probed (for each 'gid' in 'sample_gid_list')
all_trace_probes = []
for i in range(len(self.sample_gid_list)):
trace_probes = []
if self.sample_gid_list[i] == gid:
# for every neuron to be probed (upon first specification): get membrane potential, total current, and external input currents
if self.sample_gid_list[:i+1].count(gid) == 1:
trace_probes.extend([arbor.cable_probe_membrane_voltage('"center"', "tag_v"),
arbor.cable_probe_total_ion_current_cell("tag_I_tot"),
arbor.cable_probe_point_state_cell("ou_input", "I_ou", "tag_I_ou")])
# gets the synapse identifier for the corresponding neuron identifier in 'sample_gid_list'
presyn_gid = getPresynapticId(self.sample_pre_list, i)
# for synapses to be probed: additionally get synaptic calcium concentration, early- and late-phase weight, and concentration of PRPs
if presyn_gid >= 0:
syn_target = f"syn_ee_calcium_plasticity_{presyn_gid}_{gid}"
trace_probes.extend([arbor.cable_probe_point_state(syn_target, self.plasticity_mechanism, "Ca", f"tag_Ca_{presyn_gid}"),
arbor.cable_probe_point_state(syn_target, self.plasticity_mechanism, "h", f"tag_h_{presyn_gid}"),
arbor.cable_probe_point_state(syn_target, self.plasticity_mechanism, "z", f"tag_z_{presyn_gid}"),
#arbor.cable_probe_ion_diff_concentration('"center"', "sps_", f"tag_sps__{presyn_gid}"),
#arbor.cable_probe_point_state(syn_target, self.plasticity_mechanism, "sps", f"tag_sps_{presyn_gid}"),
arbor.cable_probe_density_state('"center"', self.neuron_config["mechanism"], "spsV", f"tag_spsV_{presyn_gid}"),
#arbor.cable_probe_point_state(syn_target, self.plasticity_mechanism, "pc", f"tag_pc_{presyn_gid}"),
#arbor.cable_probe_ion_diff_concentration('"center"', "pc_", f"tag_pc__{presyn_gid}"),
arbor.cable_probe_density_state('"center"', self.neuron_config["mechanism"], "pc", f"tag_pc__{presyn_gid}")])
if self.N_exc <= 4: # only for small networks
writeAddLog(f" set probes for synapse {presyn_gid}->{gid}")
self.set_max_num_trace_probes(len(trace_probes))
all_trace_probes.extend(trace_probes)
if self.N_exc <= 4: # only for small networks
writeAddLog(" max_num_trace_probes =", self.get_max_num_trace_probes())
writeAddLog(" max_num_latest_state_probes =", self.get_max_num_latest_state_probes())
#return [*all_trace_probes, *latest_state_probes] # see the warning above
return [*latest_state_probes, *all_trace_probes]
#####################################
# runSimPhase
# Run a simulation phase given a recipe, runtime, and timestep.
# - phase: number of the current simulation phase (should equal to 1 upon first call)
# - recipe: the Arbor recipe to be used
# - config: configuration of model and simulation parameters (as a dictionary from JSON format)
# - rseed: random seed
# - t_final_red: final simulated time (reduced for the current phase)
# - adap_dt: the adapted timestep
# - platform: the hardware platform to be used (options: "CPU", "GPU")
# - plotting: specifies whether plots should be created or not
# - return: spike data and trace data
def runSimPhase(phase, recipe, config, rseed, t_final_red, adap_dt, platform, plotting):
###############################################
# initialize
s_desc = config['simulation']['short_description'] # short description of the simulation, including hashtags for benchmarking
sample_gid_list = config['simulation']['sample_gid_list'] # list of the neurons that are to be probed (given by number/gid)
sample_pre_list = config['simulation']['sample_pre_list'] # list of the presynaptic neurons for probing synapses (one value for each value in `sample_gid`; -1: no synapse probing)
sample_curr = config['simulation']['sample_curr'] # pointer to current data (0: total membrane current, 1: OU stimulation current, 2: OU background noise current)
output_period = config['simulation']['output_period'] # sampling size in timesteps (every "output_period-th" timestep, data will be recorded for plotting)
loc = 0 # for single-compartment neurons, there is only one location
###############################################
# prepare domain decomposition and simulation
writeLog("CPU:", subprocess.check_output("lscpu | grep \"^Model name:\"", shell=True).decode()
.replace('\n', ' ')
.replace('Model name:', '')
.strip())
if platform == "CPU":
alloc = arbor.proc_allocation(threads=1, gpu_id=None) # select one thread and no GPU (default; cf. https://docs.arbor-sim.org/en/v0.7/python/hardware.html#arbor.proc_allocation)
hint = arbor.partition_hint()
hint.cpu_group_size = 1 #recipe.num_cells()
context = arbor.context(alloc, mpi=None) # constructs a local context without MPI connection
#context = arbor.context(threads="avail_threads", gpu_id=0)
writeLog("GPU: none")
elif platform == "GPU":
#alloc = arbor.proc_allocation(threads=1, gpu_id=None) # select one thread and no GPU (default; cf. https://docs.arbor-sim.org/en/v0.7/python/hardware.html#arbor.proc_allocation)
#alloc = arbor.proc_allocation(threads="avail_threads", gpu_id=0) # select one thread and no GPU (default; cf. https://docs.arbor-sim.org/en/v0.7/python/hardware.html#arbor.proc_allocation)
hint = arbor.partition_hint()
hint.prefer_gpu = True
hint.gpu_group_size = 500 #recipe.num_cells()
#context = arbor.context(alloc, mpi=None) # constructs a local context without MPI connection
context = arbor.context(threads=1, gpu_id=0)
writeLog("GPU:", subprocess.check_output("nvidia-smi -L", shell=True).decode()
.replace('\n', '\n ')
.strip())
else:
raise ValueError(f"Unknown platform: '{platform}'")
meter_manager = arbor.meter_manager()
meter_manager.start(context)
hints = {arbor.cell_kind.cable: hint}
domains = arbor.partition_load_balance(recipe, context, hints) # constructs a domain_decomposition that distributes the cells in the model described by an arbor.recipe over the distributed and local hardware resources described by an arbor.context (cf. https://docs.arbor-sim.org/en/v0.5.2/python/domdec.html#arbor.partition_load_balance)
meter_manager.checkpoint('load-balance', context)
sim = arbor.simulation(recipe, context, domains, seed = rseed)
# output, to be printed only once
if phase == 1:
num_conn_str = f"Number of connections: {int(recipe.connections_counter)}" + \
f" (expected value: {round(recipe.p_c * (recipe.N_tot**2 - recipe.N_tot), 1) if recipe.p_c > 0 else 'n/a'}"
if config['populations']['conn_file']:
num_conn_str += f", loaded from '{config['populations']['conn_file']}')"
else:
num_conn_str += ", generated)"
writeLog(num_conn_str)
writeLog(context)
writeLog(hint)
writeLog(domains)
# set metering checkpoint
meter_manager.checkpoint('simulation-init', context)
# create sampling schedules
writeLog(f"Creating schedule with adap_dt = {adap_dt.value_as(U.ms)} ms, output_period = {output_period}.")
reg_sched = arbor.regular_schedule(0*U.ms, output_period*adap_dt, t_final_red) # create schedule for recording traces (of either rich or fast-forward dynamics)
final_sched = arbor.explicit_schedule([t_final_red-adap_dt]) # create schedule for recording the final timestep
#sampl_policy = arbor.sampling_policy.lax # use exact policy, just to be safe!?
###############################################
# set handles
handle_mem, handle_tot_curr, handle_stim_curr, handle_Ca_specific, handle_h_specific, handle_z_specific, handle_sps_specific, handle_p_specific = [], [], [], [], [], [], [], []
handle_h_syn_latest, handle_z_syn_latest, handle_p_comp_latest = [], [], []
synapse_mapping = [[] for i in range(recipe.N_exc)] # list that contains a list of matrix indices indicating the presynaptic neurons for each (excitatory) neuron -> sparse representation
max_num_latest_state_probes = recipe.get_max_num_latest_state_probes() # the number of latest-state probes before specific sampling probes begin (for excitatory neurons)
writeAddLog(f"Setting handles with max_num_latest_state_probes = {max_num_latest_state_probes}.")
# loop over all compartments and all potential synapses to set handles for latest state
if max_num_latest_state_probes > 0: # only if latest state probes exist
for i in range(recipe.N_exc):
handle_p_comp_latest.append(sim.sample(i, "tag_pc_latest", final_sched)) # neuronal protein amount
if recipe.N_exc <= 4: # only for small networks
writeAddLog(f"Setting handles for latest state of neuron {i}...\n" +
f" handle(p) = {handle_p_comp_latest[-1]} " + str(sim.probe_metadata(i, "tag_pc_latest")))
for j in range(recipe.N_exc):
# check if synapse is supposed to exist
if recipe.conn_matrix[i][j]:
#num_sample_syn = len(synapse_mapping[i]) # number of synapses that have so far been found for postsynaptic neuron 'i'
handle_h_syn_latest.append(sim.sample(i, f"tag_h_latest_{j}", final_sched)) # early-phase weight
handle_z_syn_latest.append(sim.sample(i, f"tag_z_latest_{j}", final_sched)) # late-phase weight
if recipe.N_exc <= 4: # only for small networks
writeAddLog(f"Setting handles for latest state of synapse {j}->{i}...\n" +
f" handle(h) = {handle_h_syn_latest[-1]} " + str(sim.probe_metadata(i, f"tag_h_latest_{j}")) + "\n" +
f" handle(z) = {handle_z_syn_latest[-1]} " + str(sim.probe_metadata(i, f"tag_z_latest_{j}")))
synapse_mapping[i].append(j)
# loop over elements in 'sample_gid_list' (and thereby, 'sample_pre_list') to set handles for specific sampling
for i in range(len(sample_gid_list)):
# retrieve the current neuron index
sample_gid = sample_gid_list[i]
# retrieve the synapse index by the corresponding neuron identifier
sample_presyn_gid = getPresynapticId(sample_pre_list, i)
writeAddLog(f"Setting handles for specific sample #{i} (neuron {sample_gid}, " +
f"synapse {sample_presyn_gid}->{sample_gid})...")
# for all neurons: get membrane potential and current(s)
handle_mem.append(sim.sample(sample_gid, "tag_v", reg_sched)) # membrane potential
handle_tot_curr.append(sim.sample(sample_gid, "tag_I_tot", reg_sched)) # total current
handle_stim_curr.append(sim.sample(sample_gid, "tag_I_ou", reg_sched)) # stimulus current
writeAddLog(f" handle(V) = {handle_mem[-1]} " + str(sim.probe_metadata(sample_gid, "tag_v")) +
f"\n handle(I_tot) = {handle_tot_curr[-1]} " + str(sim.probe_metadata(sample_gid, "tag_I_tot")) +
f"\n handle(I_stim) = {handle_stim_curr[-1]} " + str(sim.probe_metadata(sample_gid, "tag_I_ou")))
# for excitatory neurons with synapses to be probed: get synapse data
if sample_presyn_gid >= 0: # synapse probes exist only if this condition is true
handle_Ca_specific.append(sim.sample(sample_gid, f'tag_Ca_{sample_presyn_gid}', reg_sched)) # calcium amount
handle_h_specific.append(sim.sample(sample_gid, f'tag_h_{sample_presyn_gid}', reg_sched)) # early-phase weight
handle_z_specific.append(sim.sample(sample_gid, f'tag_z_{sample_presyn_gid}', reg_sched)) # late-phase weight
handle_sps_specific.append(sim.sample(sample_gid, f'tag_spsV_{sample_presyn_gid}', reg_sched)) # signal triggering protein synthesis
handle_p_specific.append(sim.sample(sample_gid, f'tag_pc__{sample_presyn_gid}', reg_sched)) # neuronal protein amount
writeAddLog(f" handle(Ca) = {handle_Ca_specific[-1]} " + str(sim.probe_metadata(sample_gid, f"tag_Ca_{sample_presyn_gid}")) +
f"\n handle(h) = {handle_h_specific[-1]} " + str(sim.probe_metadata(sample_gid, f"tag_h_{sample_presyn_gid}")) +
f"\n handle(z) = {handle_z_specific[-1]} " + str(sim.probe_metadata(sample_gid, f"tag_z_{sample_presyn_gid}")) +
f"\n handle(sps) = {handle_sps_specific[-1]} " + str(sim.probe_metadata(sample_gid, f"tag_spsV_{sample_presyn_gid}")) +
f"\n handle(p) = {handle_p_specific[-1]} " + str(sim.probe_metadata(sample_gid, f"tag_pc__{sample_presyn_gid}")))
writeAddLog("Number of set handles:" +
f"\n len(handle_h_syn_latest) = {len(handle_h_syn_latest)}" +
f"\n len(handle_z_syn_latest) = {len(handle_z_syn_latest)}" +
f"\n len(handle_p_comp_latest) = {len(handle_p_comp_latest)}" +
f"\n len(handle_mem) = {len(handle_mem)}" +
f"\n len(handle_tot_curr) = {len(handle_tot_curr)}" +
f"\n len(handle_stim_curr) = {len(handle_stim_curr)}" +
f"\n len(handle_Ca_specific) = {len(handle_Ca_specific)}" +
f"\n len(handle_h_specific) = {len(handle_h_specific)}" +
f"\n len(handle_z_specific) = {len(handle_z_specific)}" +
f"\n len(handle_sps_specific) = {len(handle_sps_specific)}" +
f"\n len(handle_p_specific) = {len(handle_p_specific)}")
meter_manager.checkpoint('handles-set', context)
#arbor.profiler_initialize(context) # can be enabled if Arbor has been compiled with -DARB_WITH_PROFILING=ON
###############################################
# run the simulation
if plotting: # show the progress banner only if plotting is enabled (not on cluster)
sim.progress_banner()
sim.record(arbor.spike_recording.all)
#sim.set_binning_policy(arbor.binning.regular, adap_dt)
sim.run(tfinal=t_final_red, dt=adap_dt) # recall
meter_manager.checkpoint('simulation-run', context)
# output of metering and profiler summaries
writeAddLog("Metering summary:\n", arbor.meter_report(meter_manager, context))
#writeAddLog("Profiler summary:\n", arbor.profiler_summary()) # can be enabled if Arbor has been compiled with -DARB_WITH_PROFILING=ON
###############################################
# get and store the latest state of all compartments and synapses (for the next simulation phase)
if max_num_latest_state_probes > 0: # only if latest state probes exist
h_synapses = np.zeros((recipe.N_exc, recipe.N_exc))
z_synapses = np.zeros((recipe.N_exc, recipe.N_exc))
p_compartments = np.zeros(recipe.N_exc)
state_path = getDataPath(s_desc, f"state_after_phase_{phase}")
if not os.path.isdir(state_path): # if this directory does not exist yet
os.mkdir(state_path)
# loop over all compartments and synapses and set handles to get latest state
n = 0 # 1-dim of synapse (for sparse representation)
for i in range(recipe.N_exc):
p_samples = sim.samples(handle_p_comp_latest[i])
if len(p_samples) > 0:
tmp, _ = p_samples[loc]
p_compartments[i] = tmp[:, 1].squeeze()
else:
raise Exception(f"No data for handle_p_comp_latest[{i}] = {handle_p_comp_latest[i]}")
# loop over synapses
for num_sample_syn in range(len(synapse_mapping[i])):
j = synapse_mapping[i][num_sample_syn]
if recipe.N_exc <= 4: # only for small networks
writeAddLog(f"Getting latest state of synapse {j}->{i}...")
if recipe.conn_matrix[i][j]:
h_samples = sim.samples(handle_h_syn_latest[n])
z_samples = sim.samples(handle_z_syn_latest[n])
if len(h_samples) > 0:
tmp, _ = h_samples[loc]
h_synapses[i][j] = tmp[:, 1].squeeze()
else:
raise Exception(f"No data for handle_h_syn_latest[{n}] = {handle_h_syn_latest[n]}")
if len(z_samples) > 0:
tmp, _ = z_samples[loc]
z_synapses[i][j] = tmp[:, 1].squeeze()
else:
raise Exception(f"No data for handle_z_syn_latest[{n}] = {handle_z_syn_latest[n]}")
n += 1
else:
raise Exception(f"Entry with i={i}, j={j} not found in connectivity matrix.")
np.savetxt(os.path.join(state_path, "h_synapses.txt"), h_synapses.T, fmt="%.4f")
np.savetxt(os.path.join(state_path, "z_synapses.txt"), z_synapses.T, fmt="%.4f")
np.savetxt(os.path.join(state_path, "p_compartments.txt"), p_compartments, fmt="%.4f")
###############################################
# get data traces
sample_gid_list = config['simulation']['sample_gid_list'] # list of the neurons that are to be probed (given by number/gid)
sample_curr = config['simulation']['sample_curr'] # pointer to current data (0: total membrane current, 1: OU stimulation current, 2: OU background noise current)
data, times = [], []
data_mem, data_curr, data_Ca, data_h = [], [], [], []
data_z, data_sps, data_p = [], [], []
dimin = 0
# loop over elements in 'sample_gid_list' (and thereby, 'sample_pre_list') to get traces
for i in range(len(sample_gid_list)):
# retrieve the current neuron index
sample_gid = sample_gid_list[i]
# retrieve the presynaptic neuron identifier by the corresponding index
sample_presyn_gid = getPresynapticId(sample_pre_list, i)
writeLog(f"Getting traces of sample #{i} (neuron {sample_gid}, " +
f"synapse {sample_presyn_gid}->{sample_gid})...")
# get neuronal data
data, _ = sim.samples(handle_mem[i])[loc]
if len(sim.samples(handle_mem[i])) <= 0:
raise ValueError(f"No data for handle_mem[{i}] (= {handle_mem[i]}).")
times = data[:, 0]
data_mem.append(data[:, 1])
writeAddLog(f" Read from handle {handle_mem[i]} (V, {bool(sim.samples(handle_mem[i]))})")
if sample_curr > 0: # get OU stimulation current specified by sample_curr