-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathnewnetclamp.hoc
1079 lines (899 loc) · 38.6 KB
/
newnetclamp.hoc
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
/************************************************************
For more information, consult ModelDoc.pdf
************************************************************/
loadstart = startsw() // record the start time of the set up
/***********************************************************************************************
I. LOAD LIBRARIES
***********************************************************************************************/
{load_file("nrngui.hoc")} // Standard definitions - NEURON library file
{load_file("netparmpi.hoc")} // Contains the template that defines the properties of
// the ParallelNetManager class, which is used to set
// up a network that runs on parallel processors
{load_file("./setupfiles/ranstream.hoc")} // Contains the template that defines a RandomStream
// class used to produce random numbers
// for the cell noise (what type of noise?)
{load_file("./setupfiles/CellCategoryInfo.hoc")} // Contains the template that defines a
// CellCategoryInfo class used to store
// celltype-specific parameters
{load_file("./setupfiles/SynStore.hoc")} // Contains the template that defines a
{load_file("./setupfiles/defaultvar.hoc")} // Contains the proc definition for default_var proc
{load_file("./setupfiles/parameters.hoc")} // Loads in operational and model parameters that can
// be changed at command line
{load_file("./setupfiles/set_other_parameters.hoc")}// Loads in operational and model parameters
// that can't be changed at command line
//{load_file("./setupfiles/check_for_files.hoc")} // Checks that the necessary chosen files exist
// (datasets, connectivity, and stimulation)
default_var("gidOI", 0) // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("cellindOI", 0) // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("spkflag", 0) // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("stimflag", 0) // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("resultsfolder", "00001") // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("runname", "None") // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("origRunName", "CutSma_040_11") // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("celltype2use", "axoaxoniccell") // Run Comments // -c 'strdef RunName ' 'RunName = "Woo"'
default_var("relpath"," /cygdrive/c/Users/maria_000/Documents/repos/ca1") //"/home/casem/repos/ca1")
default_var("resultspath","networkclamp_results")
default_var("sl","/")
celltypeOI=cellindOI
myi_flag=1
//gidOI=0
i=0
ij=0
gid=0
/*strdef origRunName
origRunName = "CutSma_040_11"
PrintConnDetails=2*/
/***********************************************************************************************
II. SET MODEL SIZE, CELL DEFINITIONS
***********************************************************************************************/
celsius=34
{load_file("./setupfiles/load_cell_category_info.hoc")} // Reads the 'cells2include.hoc' file and
// loads the information into one
// 'CellCategoryInfo' object for each cell
// type (bio or art cells?). Info includes
// number of cells, gid ranges, type name
strdef path2ConnData
sprint(path2ConnData, "../networkclamp_results/%s/%s/conndata.dat", origRunName, resultsfolder)
strdef path2SynData
sprint(path2SynData, "../networkclamp_results/%s/%s/syndata.dat", origRunName, resultsfolder)
{load_file("./setupfiles/load_cell_conns.hoc")} // Load in the cell connectivity info
{load_file("./setupfiles/load_cell_syns.hoc")} // Load in the cell connectivity info
strdef tempFileStr // Define a string reference to store the name of the
// current cell template file
proc loadCellTemplates(){local i // Proc to load the template that defines (each) cell class
//for i=0, numCellTypes-1 { // Iterate over each cell type in cells2include (and art cells?)
i = celltypeOI
sprint(tempFileStr,"./cells/class_%s.hoc",celltype2use) // Concatenate the
//sprint(tempFileStr,"./cells/class_%s.hoc",cellType[i].technicalType) // Concatenate the
// path and file
load_file(tempFileStr) // Load the file with the template that defines the class
// for each cell type
//}
if (stimflag==1) {
load_file("./cells/class_sintrain.hoc") // Load the file with the template that defines the class
//tstop = 100
//SimDuration = tstop
print "Using PhasicData, SimDuration = ", SimDuration
} else {
load_file("./cells/class_ppvec.hoc") // Load the file with the template that defines the class
}
}
loadCellTemplates() // Run the newly defined proc
proc calcNetSize(){local i // Calculate the final network size (after any cell death)
//cellType[0].numCells = cellType[numCellTypes-1].cellEndGid - cellType[0].cellEndGid + 4*2 // just added now for spontburst
//cellType[0].updateGidRange(0)
totalCells = 0 // Initialize totalCells (which counts the number of 'real'
// cells) so we can add to it iteratively in the 'for' loop
ncell = cellType[0].numCells // Initialize ncell (which counts all 'real' and 'artificial'
// cells) so we can add to it iteratively in the 'for' loop
for i=1,numCellTypes-1 { // Run the following code for 'real' cell types only - need a different way of singling out real cells?
if (cellType[i].layerflag==0) { // For cell types in layer 0, which are susceptible to death
cellType[i].numCells = int(cellType[i].numCells * ((100-PercentCellDeath)/100))
// Calculate the number of cells surviving after cell loss
if (cellType[i].numCells == 0) {cellType[i].numCells = 1} // If all cells of one type
// are killed, let 1 live
}
cellType[i].updateGidRange(cellType[i-1].cellEndGid+1) // Update the gid range for each
// cell type
totalCells = totalCells + cellType[i].numCells // Update the total number of cells
// after sclerosis, not including
// artificial cells
ncell = ncell + cellType[i].numCells // Update the total number of cells
// after sclerosis, including
// artificial cells
}
}
calcNetSize()
proc calcBinSize(){local NumGCells
for i=0, numCellTypes-1 { // Using the specified dimensions of the network (in um) and
// the total number of cells of each type, set the number
// of bins in X, Y, Z dimensions such that the cells will be
// evenly distributed throughout the allotted volume
// just changed this so even the stim cells will be allotted, as now we have some
// stimulation protocols that incorporate stim cell position
cellType[i].setBins(scLongitudinalLength,scTransverseLength,LayerVector.x[cellType[i].layerflag])
// For the z length, use the height of the layer in which the
// cell somata are found for this cell type
}
}
calcBinSize()
/***********************************************************************************************
III.SET UP PARALLEL CAPABILITY AND WRITE OUT RUN RECEIPT
***********************************************************************************************/
objref pnm, pc, nc, nil
proc parallelizer() {
pnm = new ParallelNetManager(ncell) // Set up a parallel net manager for all the cells
pc = pnm.pc
pnm.round_robin() // Incorporate all processors - cells 0 through ncell-1
// are distributed throughout the hosts
// (cell 0 goes to host 0, cell 1 to host 1, etc)
}
parallelizer()
iterator pcitr() {local i2, startgid // Create iterator for use as a standard 'for' loop
// throughout given # cells usage:
// for pcitr(&i1, &i2, &gid, it_start, it_end) {do stuff}
// it_start and it_end let you define range over
// which to iterate
// i1 is the index of the cell on the cell list for that host
// i2 is the index of that cell for that cell type on that host
numcycles = int($4/pc.nhost)
extra = $4%pc.nhost
addcycle=0
if (extra>pc.id) {addcycle=1}
startgid=(numcycles+addcycle)*pc.nhost+pc.id
i1 = numcycles+addcycle // the index into the cell # on this host.
i2 = 0 // the index of the cell in that cell type's list on that host
if (startgid<=$5) {
for (i3=startgid; i3 <= $5; i3 += pc.nhost) { // Just iterate through the cells on
// this host(this simple statement
// iterates through all the cells on
// this host and only these cells because
// the roundrobin call made earlier dealt
// the cells among the processors in an
// orderly manner (like a deck of cards)
$&1 = i1
$&2 = i2
$&3 = i3
// if (i3 == $5) {iterator_statement}
iterator_statement
i1 += 1
i2 += 1
}
}
}
iterator fastitr() {local i2, startgid // Create iterator for use as a standard 'for' loop
// throughout given # cells usage:
// for fastitr(&i1, &i2, &gid, it_start, it_end) {do stuff}
// it_start and it_end let you define range over
// which to iterate
// i1 is the index of the cell on the cell list for that host
// i2 is the index of that cell for that cell type on that host
numcycles = int($4/pc.nhost)
extra = $4%pc.nhost
addcycle=0
if (extra>pc.id) {addcycle=1}
startgid=(numcycles+addcycle)*pc.nhost+pc.id
i1 = numcycles+addcycle // the index into the cell # on this host.
i2 = 0 // the index of the cell in that cell type's list on that host
if (startgid<=$5) {
for (i3=startgid; i3 <= $5; i3 += pc.nhost) { // Just iterate through the cells on
// this host(this simple statement
// iterates through all the cells on
// this host and only these cells because
// the roundrobin call made earlier dealt
// the cells among the processors in an
// orderly manner (like a deck of cards)
$&1 = i1
$&2 = i2
$&3 = i3
if (i3 == $5) {iterator_statement}
// iterator_statement
i1 += 1
i2 += 1
}
}
}
objref strobj
strobj = new StringFunctions()
strdef direx
if (strcmp(UID,"0")==0 && pc.id==0) {
type = unix_mac_pc() // 1 if unix, 2 if mac, 3 if mswin, or 4 if mac osx darwin
if (type==3) {
{system("cscript //NoLogo setupfiles/uuid.vbs", direx)} // pc
} else {
{system("uuidgen", direx)} // unix or mac
strobj.left(direx, strobj.len(direx)-1)
}
UID = direx
}
strdef cmdstr, cmd
//{load_file("./setupfiles/save_run_info.hoc")}
loadtime = startsw() - loadstart // Calculate the set up time (now - recorded start time) in seconds
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (set up)\n************\n", loadtime)}
createstart = startsw() // Record the start time of the cell creation
/***********************************************************************************************
IV. CREATE, UNIQUELY ID, AND POSITION CELLS
***********************************************************************************************/
objref memfile, topfile //, pccc
//{pccc = new ParallelContext()}
{zzz=0}
if (pc.id==0) {
{sprint(cmd,"./networkclamp_results/%s/%s/memory.dat", runname, resultsfolder)}
{memfile = new File(cmd)}
{memfile.wopen()}
{sprint(cmd,"./networkclamp_results/%s/%s/topoutput.dat", runname, resultsfolder)}
{topfile = new File(cmd)}
{topfile.wopen()}
}
strdef direx1, direx2, TopCommand
func mallinfo() {local m // this function is a wrapper for the system function mallinfo
m = nrn_mallinfo(0)
if (pc.id == 0) {
if (PrintTerminal>1) {printf("Memory (host 0) - %s: %ld. Since last call: %ld\n", $s2, m, m - $1)}
memfile.printf("%s\t%f\t%ld\t%ld\n", $s2, startsw()-loadstart, m, m - $1) // Print the spike time and spiking cell gid
sprint(TopCommand,"top -p `pgrep %s | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1", TopProc)
//print "TopCommand = ", TopCommand
if (strcmp(TopProc,"")!=0) { system(TopCommand, direx1)}
//system("top -p `pgrep special | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1",direx1)
//system("top -p `pgrep nrniv | head -n20 | tr \"\\n\" \",\" | sed 's/,$//'` -n1 -b | tail -n2 | head -n1",direx2)
//if (strcmp(direx1,direx2)<0) {
// topfile.printf("%s\t%s\n", $s2, direx2) // Print the spike time and spiking cell gid
//} else {
topfile.printf("%s\t%s\n", $s2, direx1) // Print the spike time and spiking cell gid
//}
}
return m
}
{zzz = mallinfo(zzz, "Prior to creating cells")}
objref cells, ransynlist, ranstimlist, raninitlist
cells = new List()
ransynlist = new List()
ranstimlist = new List()
raninitlist = new List()
{load_file("./setupfiles/clamp/netclamp_create_cells_pos.hoc")} // Creates each cell on its assigned host
// and sets its position using the algorithm
// defined above
strdef cmd
//new
strdef mname, tmpstr
objref strobj
objref mechstring[9], mechlength
proc print_mechs(){local i, k localobj mt, cell // this code courtesy of Jose Ambros-Ingerson via the NEURON forum
mechlength = new Vector(1)
strobj = new StringFunctions()
cell = pc.gid2cell($1)
access cell.soma {
mt = new MechanismType(0)
objref mechstring[mt.count()]
k = 0
for i=0, mt.count()-1 {
mt.select( i )
mt.selected(mname)
if( ismembrane(mname)) {
if (strcmp(mname,"capacitance")!=0 && strcmp(mname,"morphology")!=0 && strcmp(mname,"iconc_CaZ")!=0 && strcmp(mname,"iconc_Ca")!=0 && strcmp(mname,"iconcCa")!=0 && strcmp(mname,"ccanl")!=0 && strcmp(mname,"cad")!=0 && strcmp(mname,"pas")!=0 && strcmp(mname,"vmax")!=0 && strcmp(mname,"extracellular")!=0 ) { //
if (strobj.substr(mname,"_ion")==-1) {
//printf("myi_%s \n", mname)
sprint(tmpstr, "myi_%s", mname) // "cell.soma.%s(0.5)", tmpstr
mechstring[k] = new String()
mechstring[k].s = tmpstr
{k = k+1}
}
}
}
}
}
//{printf("mt count = %g, mechlength0 = %g mechstring = %s\n", mt.count(), k-1, mechstring[k-1].s)}
{mechlength.x[0] = k}
//{printf("'bout to leave. mechstring count = %g \n", mechstring.count)}
//return mechstring
}
//newend
createtime = startsw() - createstart // Calculate time taken to create the cells
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (created cells)\n************\n", createtime)}
connectstart = startsw() // Grab start time of cell connection
oldtimeout = pc.timeout(0)
zzz = mallinfo(zzz, "After creating cells")
/***********************************************************************************************
V. CONNECT THE MODEL CELLS AND CONNECT THE STIMULATION CELLS TO SOME MODEL CELLS
***********************************************************************************************/
celsius=34
{load_file("./setupfiles/nc_append_functions.hoc")} // Defines nc_append and nc_appendo, which
// are procs (?) needed to create the netcon
// object associated with each connection
// (synapse)
{pc.barrier()}
print "got here?"
{load_file("./setupfiles/clamp/list_randfaststim_connections.hoc")} // Loads in a particular definition of the connectCells function,
// depending on the value of the Connectivity string
print "or here?"
zzz = mallinfo(zzz, "After connecting cells")
{pc.barrier()}
if (stimflag==0) {
{load_file("./setupfiles/clamp/netclamp_stimulation.hoc")} // Loads in a particular stimulation protocol, including the spike
// vectors for the artificial cells and their connections to the
// 'real' cells (where do we define the # of art cells and create them?)
} else {
if (stimflag==1) {
{load_file("./setupfiles/clamp/netclamp_oscillation_stimulation.hoc")}
}
}
zzz = mallinfo(zzz, "After defining stimulation")
print "got there!"
{load_file("./setupfiles/write_conns_functions.hoc")} // Defines two funcs (procs?): printNumConFile()
// writes the number of synapses between every
// combination of pre- and post-synaptic cell type
// and tracenet() lists the presynaptic cell type,
// postsynaptic cell type, and synapse type for
// every connection
if (PrintConnSummary==1) {printNumConFile()} // Write "numcons.dat": a QUICK, SMALL summary file of the # of
// connections by pre and post cell types
/*
if (PrintConnDetails==2) {allCellConnsParallel()} // Write "connections.dat": a VERY SLOW, LARGE file of all cell
// connections (pre and post gids, synapse type, host on which
// (target?) cell exists)
if (PrintConnDetails==1) {allCellConnsParallel()} // Write the file of all cell connections (pre and post gids, synapse type,
// host on which cell exists), "connections.dat"
if ((PrintConnDetails==0) && (PrintVoltage==2)) {recordedCellConnsParallel()} // If detailed connections are not being written out, at least write them out for cells that are being traced
*/
connecttime = startsw() - connectstart // calculate time taken to connect the cells
if (pc.id == 0) {printf("\nTIME HOST 0: %g seconds (connected cells)\n************\n", connecttime)}
/***********************************************************************************************
VI. INITIALIZE AND RUN NETWORK, OUTPUT RESULT FILES
***********************************************************************************************/
//new
strdef mname, headstr, outfile, cmdstr
objref mt, cell, strobj, f
strobj = new StringFunctions()
cellmechs=1
if (cellmechs==1) {
cellind = cellindOI
if (pc.gid_exists(gidOI)) {
// add recording locations to cell:
cell = pc.gid2cell(gidOI)
if (cell.is_art==0) { resultsfolder
sprint(outfile, "./networkclamp_results/%s/%s/celldata_%s.dat", origRunName, resultsfolder, cellType[cellind].cellType_string)
//sprint(outfile, "%s%s%s%scelldata_%s.dat", relpath, sl, resultspath, sl, cellType[cellind].cellType_string)
f = new File(outfile)
f.wopen()
access cell.soma // access cell.soma[0]
distance() // set the origin for measuring the distance
sprint(headstr,"secname\tpos\tx\ty\tz\tdist\tdiam\tcm\tRa")
mt = new MechanismType(0)
for i=0, mt.count()-1 {
mt.select( i )
mt.selected(mname)
if (strobj.substr(mname,"ch_")>-1) {
sprint(headstr,"%s\tgmax_%s", headstr, mname)
}
}
sprint(headstr,"%s\n", headstr)
f.printf(headstr)
forsec cell.all { // forsec cell.all // forall
for (x) {
xdist = distance(x)
sprint(headstr,"%s\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g", secname(), x, x3d(0), y3d(0), z3d(0), xdist, diam, cm, Ra)
for i=0, mt.count()-1 {
mt.select( i )
mt.selected(mname)
if (strobj.substr(mname,"ch_")>-1) {
g=0
e=0
if ( ismembrane(mname)) {
sprint(cmdstr, "g = gmax_%s", mname)
execute1(cmdstr)
}
sprint(headstr,"%s\t%g", headstr, g)
}
}
sprint(headstr,"%s\n", headstr)
f.printf(headstr)
}
}
f.close()
}
}
}
//newend
objref mytrace, myvrec, cell
objref lfpkmatrix
lfpkmatrix = new Vector()
cell = pc.gid2cell(gidOI)
/*********************
inserted stuff right below here:
************/
/*
Computes xyz coords of nodes in a model cell whose topology & geometry are defined by pt3d data.
Code by Ted Carnevale.
*/
proc interpxyz() { local ii, nn, kk, xr, nsegs localobj xx, yy, zz, xint, yint, zint, ll, range
// get the data for the section
nn = $1
nsegs = $2
xx = $o3
yy = $o4
zz = $o5
ll = $o6
xint = $o7
yint = $o8
zint = $o9
// to use Vector class's .interpolate()
// must first scale the independent variable
// i.e. normalize length along centroid
ll.div(ll.x[nn-1])
// initialize the destination "independent" vector
range = new Vector(nsegs+2)
range.indgen(1/nsegs)
range.sub(1/(2*nsegs))
range.x[0]=0
range.x[nsegs+1]=1
// length contains the normalized distances of the pt3d points
// along the centroid of the section. These are spaced at
// irregular intervals.
// range contains the normalized distances of the nodes along the
// centroid of the section. These are spaced at regular intervals.
// Ready to interpolate.
xint.interpolate(range, ll, xx)
yint.interpolate(range, ll, yy)
zint.interpolate(range, ll, zz)
}
proc setup_npole_lfp() { local i, j, k, m, n, r, h, s, avglfp, typei, celltype, ex, ey, ez, sx, sy, sz, lfptype, ii, nn localobj myr, dists, lst, xx, yy, zz, ll, xint, yint, zint strdef ks
rho = 333.0 // extracellular resistivity, [ohm cm]
fdst = 0.1 // percent of distant cells to include in the computation
ex = mypoint.x[0] //$1
ey = mypoint.x[1] //$2
ez = mypoint.x[2] //$3
//printf ("host %d: entering setup_npole_lfp\n", pc.id)
// Vector for storing longitudinal and perpendicular distances
// between recording electrode and compartments
dists = new Vector(2)
m = 0
n = 0
// Relative to the recording electrode position
forsec cell.all {
insert extracellular
n = n + nseg
}
lfpkmatrix = new Vector(n+1)
//printf ("host %d: npole_lfp: lfpkmatrix.nrow = %d lfpkmatrix.ncol = %d\n", pc.id, lfpkmatrix.nrow, lfpkmatrix.ncol)
// Iterates over each compartment of the cell
j = 0
forsec cell.all {
if (ismembrane("extracellular")) {
nn = n3d()
xx = new Vector(nn)
yy = new Vector(nn)
zz = new Vector(nn)
ll = new Vector(nn)
for ii = 0,nn-1 {
xx.x[ii] = x3d(ii)
yy.x[ii] = y3d(ii)
zz.x[ii] = z3d(ii)
ll.x[ii] = arc3d(ii)
}
xint = new Vector(nseg+2)
yint = new Vector(nseg+2)
zint = new Vector(nseg+2)
interpxyz(nn,nseg,xx,yy,zz,ll,xint,yint,zint)
j = 0
for (x,0) {
sx = xint.x[j]
sy = yint.x[j]
sz = zint.x[j]
// l = L/nseg is compartment length
// r is the perpendicular distance from the electrode to a line through the compartment
// h is longitudinal distance along this line from the electrode to one end of the compartment
// s = l + h is longitudinal distance to the other end of the compartment
l = L/nseg
r = sqrt((ex-sx)*(ex-sx) + (ey-sy)*(ey-sy) + (ez-sz)*(ez-sz))
h = l/2
s = l + h
k = 0.0001 * area(x) * (rho / (4 * PI * l)) * abs(log((sqrt(h*h + r*r) - h) / (sqrt(s*s + r*r) - s)))
sprint(ks,"%g",k)
if ((strcmp(ks,"nan") == 0) || (strcmp(ks,"-nan") == 0)) {
k = 0
}
lfpkmatrix.x[j] = k
j = j + 1
}
}
}
}
setup_npole_lfp()
/************** done inserting *********/
func npole_lfp() { local vlfp
vlfp = 0 // Initialize the LFP variable
forsec cell.all {
if (ismembrane("extracellular")) {
j = 0
for (x,0) {
vlfp = vlfp + (i_membrane(x) * lfpkmatrix.x[j])
j = j + 1
}
}
} // vlfp = cell.apical[cell.NumApical-1].v(0.5) - cell.basal[0].v(0.5)
return vlfp
}
objref lfplist, lfptracelist, vec
lfplist = new List()
lfptracelist = new List()
lfp_dt = 0.5 // sampling interval for LFP
objref mytraceidxlist, lfptracelist
proc sample_lfp() { local i, avglfp, startgid, endgid, gid, celltype localobj vec, lst
avglfp = npole_lfp(MaxEDist) // x, y, z position of electrode in microns
if (pc.id == 0) {
vec = new Vector()
vec.append(t, avglfp)
lfplist.append(vec.c)
}
cvode.event(t + lfp_dt, "sample_lfp()") // execute sample_lfp lfp_dt ms from now
}
strdef filestr
sprint(filestr, "./networkclamp_results/%s/%s/allsyns.dat", origRunName, resultsfolder)
objref fAll, fSyn
fAll = new File(filestr)
fAll.wopen()
sprint(filestr, "./networkclamp_results/%s/%s/synlist.dat", origRunName, resultsfolder)
fSyn = new File(filestr)
fSyn.wopen()
strdef cmdstr, tempFileStr, sname
objref newSecRef
if (pc.gid_exists(gidOI)) {
cell = pc.gid2cell(gidOI)
print "cell is located at: ", cell.x, ", ", cell.y, ", ", cell.z
mytrace = new Vector((tstop-tstart)/dt)
mytrace.record(&cell.soma.v(0.5))
//new
if (myi_flag==1) {
cell.soma{
sname = secname()
print_mechs(gidOI)
cellind = cellindOI
for rt = 0, mechlength.x[0]-1 {
sprint(cmdstr, "{objref %s_%s}", cellType[cellind].cellType_string, mechstring[rt].s)
execute1(cmdstr)
sprint(cmdstr, "{%s_%s = new Vector(%g)}", cellType[cellind].cellType_string, mechstring[rt].s, (tstop-tstart)/dt)
execute1(cmdstr)
sprint(cmdstr, "{%s_%s.record(&%s.%s(0.5))}", cellType[cellind].cellType_string, mechstring[rt].s, sname, mechstring[rt].s)
execute1(cmdstr)
}
}
}
//newend
ist = 0
ien = 0
myi = 0
access cell.soma[0]
{distance()}
for precellType = 0, numCellTypes-1 {
ist = myi
for r=0, cellType[celltypeOI].SynList[precellType].count()-1 {
sprint(cmdstr,"newSecRef = cell.%s",cellType[celltypeOI].SynList[precellType].object(r).SecRefStr)
execute(cmdstr) // sets newSecRef
forsec newSecRef {
for (x,0) {
//y=0
execute(cellType[celltypeOI].SynList[precellType].object(r).CondStr)
if (y==1) {
fSyn.printf("%s %s %g %s\n", cellType[celltypeOI].cellType_string, cellType[precellType].cellType_string, myi, secname())
myi = myi + 1
}
}
}
}
ien = myi - 1
if (ien>=ist) {fAll.printf("%s %s %g %g\n", cellType[celltypeOI].cellType_string, cellType[precellType].cellType_string, ist, ien)}
}
fAll.close()
fSyn.close()
/*forsec cell.all {
insert extracellular
insert xtra
}
myvrec = new Vector((tstop-tstart)/dt)
myvrec.record(&lfp_xtra)*/
} else {
print "couldnt make cell ", gidOI
}
/*
load_file("clamp/interpxyz.hoc") // only interpolates sections that have extracellular
load_file("clamp/setpointers.hoc") // automatically calls grindaway() in interpxyz.hoc
load_file("clamp/calcrxc.hoc") // computes transfer r between segments and recording electrodes
*/
proc init() { local dtsav, temp, secsav, secondordersav // Redefine the proc that initializes the
// simulation (why?)
dtsav = dt // Save the desired dt value to reset it after temporarily
// changing dt to run a quick presimulation to allow the
// network to reach steady state before we start 'recording'
secondordersav = secondorder // Save the desired secondorder value to reset it after
// temporarily changing secondorder to run a quick presimulation
// to allow the network to reach steady state before we start
// 'recording' its results
finitialize(v_init) // Call finitialize from within the custom init proc, just as the default
// init proc does. Note that finitialize will call the INITIAL block for
// all mechanisms and point processes inserted in the sections and set the
// initial voltage to v_init for all sections
t = -200 // Set the start time for (pre) simulation; -500 ms to allow the network to
// reach steady state before t = 0, when the real simulation begins
dt= 10 // Set dt to a large value so that the presimulation runs quickly
secondorder = 0 // Set secondorder to 0 to set the solver to the default fully implicit backward
// euler for numerical integration (see NEURON ref)
temp= cvode.active() // Check whether cvode, a type of solver (?) is on
if (temp!=0) {cvode.active(0)} // If cvode is on, turn it off while running the presimulation
while(t<-100) { fadvance() if (PrintTerminal>2) {print t}} // Run the presimulation from t = -500
// to t = -100 (why not 0?) to let the
// network and all its components reach
// steady state. Integrate all section
// equations over the interval dt,
// increment t by dt, and repeat until
// t at -100
if (temp!=0) {cvode.active(1)} // If cvode was on and then turned off, turn it back on now
t = tstart // Set t to the start time of the simulation
dt = dtsav // Reset dt to the specified value for the simulation
secondorder = secondordersav // Reset secondorder to the specified value for the simulation
if (cvode.active()){
cvode.re_init() // If cvode is active, initialize the integrator
} else {
fcurrent() // If cvode is not active, make all assigned variables
// (currents, conductances, etc) consistent with the
// values of the states
}
frecord_init() // see email from ted - fadvance() increments the recorder, so we need to fix the index it ends up at
}
{load_file("./setupfiles/write_results_functions.hoc")} // Defines the procs spikeout() and timeout()
// which write out the spikeraster (spikeout)
// and the time taken to run each section of
// the code (timeout)
use_cache_efficient=1
get_spike_hist=0
use_bin_queue=0
use_spike_compression=0
if (use_spike_compression==1) {
maxstepval = 2.5
} else {
maxstepval = 10
}
proc rrun(){ // Run the network simulation and write out the results
local_minimum_delay = pc.set_maxstep(maxstepval) // Set every machine's max step size to minimum delay of
// all netcons created on pc using pc.gid_connect, but
// not larger than 10
if (pc.id==0 && use_spike_compression==1) {print "Host ", pc.id, " has local_minimum_delay=", local_minimum_delay}
stdinit() // Call the init fcn (which is redefined in this code) and
// then make other standard calls (to stop watch, etc)
runstart = startsw() // grab start time of the simulation
pc.psolve(tstop) // Equivalent to calling cvode.solve(tstop) but for parallel NEURON;
// solve will be broken into steps determined by the result of
// set_maxstep
runtime = startsw() - runstart // Calculate runtime of simulation
// Print a time summary message to screen
writestart = startsw()
comptime = pc.step_time
avgcomp = pc.allreduce(comptime, 1)/pc.nhost
maxcomp = pc.allreduce(comptime, 2)
if (maxcomp>0) {
if (pc.id == 0) { printf("load_balance = %.3f\n", avgcomp/maxcomp)}
if (pc.id == 0) { printf("exchange_time = %.2f ms\n", runtime - maxcomp) }
} else {
if (pc.id == 0) { printf("no load balance info available\nno spike exchange info available\n")}
}
/*
spikeoutfast()
TimeOutSerial() // Write out a file of run times for each code section
if (PrintVoltage>0) {voltageout()} // Write out any voltages that were recorded during simulation
if (PrintConnDetails==1) {allCellConnsParallel()} // Write the file of all cell connections (pre and post gids, synapse type,
// host on which cell exists), "connections.dat"
if ((PrintConnDetails==0) && (PrintVoltage==1)) {recordedCellConnsParallel()} // If detailed connections are not being written out, at least write them out for cells that are being traced
SummaryOut(runtime) // Write out a file with the total number of cells, spikes, and connections (including contributions from artificial cells)
*/
writetime = startsw() - writestart // Calculate write time of program
totaltime = startsw() - loadstart
allottedtime = JobHours*3600
if (pc.id == 0) {printf("****\nTIME SUMMARY for host 0\nset up in %g seconds\ncreated cells in %g seconds\nconnected cells in %g seconds\nran simulation in %g seconds\nwrote results in %g seconds\nTOTAL TIME = %g seconds\nALLOTTED TIME = %g seconds\n************\nThis run is called: %s %s\n************\n", loadtime, createtime, connecttime, runtime, writetime, totaltime, allottedtime, origRunName, resultsfolder)}
}
objref fih
// ComputeNpoleLFP=1 ComputeDipoleLFP=0
if (ComputeNpoleLFP > 0) {
// execute sample_lfp() at t = 0,
// right after the mechanism INITIAL blocks have been executed
fih = new FInitializeHandler("sample_lfp()")
}
objref fihw
fihw = new FInitializeHandler(2, "midbal()")
walltime = startsw()
strdef cmd, cmdo
newtstop = tstop
warningflag=0
proc midbal() {local wt, thisstep
wt = startsw()
if (t>0) {
thisstep = wt - walltime
simleft = tstop - t
compleft = JobHours*3600 - (startsw() - loadstart)
//print "t=",t,", tstop=",tstop,", simleft=",simleft,", jobtime=",JobHours*3600,", timeused=",startsw() - loadstart,", compleft=", compleft
if (warningflag==0 && (simleft/StepBy*thisstep+EstWriteTime)>compleft && pc.id == 0) {
newtstop = int((compleft-EstWriteTime)/thisstep)*StepBy + t
print "Not enough time to complete ", tstop, " ms simulation, simulation will likely stop around ", newtstop, " ms"
warningflag=1
}
if (pc.id == 0) { printf("%g ms interval at t=%g ms was %g s\n", StepBy, t, thisstep) }
if ((thisstep+EstWriteTime)>compleft && t<tstop) { // not enough time for another step, end now
{pc.barrier()}
sprint(cmd,"cat results/%s/runreceipt.txt | sed -e 's/^SimDuration = [^ ]*/SimDuration = %g;/' > x ; mv x results/%s/runreceipt.txt", RunName, t, RunName)
if (pc.id == 0) { {system(cmd,cmdo)} }
if (pc.id == 0) { print "simulation stopped early at t=", t, " ms"}
tstop = t
stoprun=1
}
}
walltime = wt
cvode.event(t+StepBy, "midbal(0)")
}
{cvode.cache_efficient(use_cache_efficient)} // always double check that this addition does not affect the spikeraster (via pointers in mod files, etc)
if (use_bin_queue==1) {
use_fixed_step_bin_queue = 1 // boolean
use_self_queue = 0 // boolean - this one may not be helpful for me, i think it's best for large numbers of artificial cells that receive large numbers of inputs
{mode = cvode.queue_mode(use_fixed_step_bin_queue, use_self_queue)}
}
if (use_spike_compression==1) {
nspike = 3 // compress spiketimes or not
gid_compress = 0 //only works if fewer than 256 cells on each proc
{nspike = pc.spike_compress(nspike, gid_compress)}
}
objref spkhist
if (get_spike_hist==1) {
spkhist = new Vector(pc.nhost)
if (pc.id==0) {
pc.max_histogram(spkhist)
}
}
zzz = mallinfo(zzz, "Before running simulation")
{rrun()} // Run the network simulation (in proc rrun)
zzz = mallinfo(zzz, "After running simulation")
if (pc.id==0) {
{memfile.close()}
{topfile.close()}
}
objref f4
strdef cmd
if (pc.id==0 && get_spike_hist==1) {
sprint(cmd,"./results/%s/spkhist.dat", RunName)
f4 = new File(cmd)
f4.wopen()
// Open for appending to file
for i=0, pc.nhost-1 {
f4.printf("%d\t%d\n", i, spkhist.x[i]) // Print the spike time and spiking cell gid
}
f4.close()
}
strdef outfile, umfile
objref f, umf
sprint(outfile, "./networkclamp_results/%s/%s/mytrace_%d_soma.dat", runname, resultsfolder, gidOI)
f = new File(outfile)
f.wopen()
f.printf("t\tv\n")
for i=0, (tstop-tstart)/dt-1 {
f.printf("%g\t%g\n", i*dt, mytrace.x[i])
}
f.close()
strdef outfile, umfile
objref f, umf
sprint(outfile, "./networkclamp_results/%s/%s/mytrace_%d_syns.dat", runname, resultsfolder, gidOI)
f = new File(outfile)
f.wopen()
f.printf("t")
for i = 0, numCellTypes-1 {
f.printf("\t%s", cellType[i].cellType_string)
}
f.printf("\n")
for ti=0, (tstop-tstart)/dt-1 {
f.printf("%g", ti*dt)
for i = 0, numCellTypes-1 {
totcur=0
for m=0,pretyperec[i].count()-1 {
totcur=totcur+pretyperec[i].o(m).x[ti]
}
f.printf("\t%g", totcur)
}
f.printf("\n")
}
f.close()
sprint(outfile, "./networkclamp_results/%s/%s/othertrace_%d_soma.dat", runname, resultsfolder, gidOI)
f = new File(outfile)
f.wopen()
f.printf("t\tv")
print_mechs(gidOI) // mechstring =
if (myi_flag==1) {
for rt = 0, mechlength.x[0]-1 {
f.printf("\t%s", mechstring[rt].s)
}
}
f.printf("\n")
for i=0, (tstop-tstart)/dt-1 {
f.printf("%g\t%g", i*dt, mytrace.x[i])
if (myi_flag==1) {
for rt = 0, mechlength.x[0]-1 {
sprint(cmdstr, "{f.printf(\"\\t%%g\", %s_%s.x[i])}", cellType[cellindOI].cellType_string, mechstring[rt].s)
execute1(cmdstr)
}
}
{f.printf("\n")}
}
f.close()
proc SpikeOutSerial() {local i, rank localobj f // Write out a spike raster (cell, spike time)
pc.barrier() // Wait for all ranks to get to this point
sprint(cmd,"./networkclamp_results/%s/%s/spikeraster.dat", runname, resultsfolder)
f = new File(cmd)
if (pc.id == 0) { // Write header to file 1 time only
f.wopen()
f.close()
}