source: project/release/4/nemo/trunk/examples/Golgi/Golgi.nemo @ 28024

Last change on this file since 28024 was 28024, checked in by Ivan Raikov, 8 years ago

nemo: overhaul of ionic concentrations interface

File size: 34.9 KB
Line 
1
2;; Reference: Theta-Frequency Bursting and Resonance in Cerebellar Granule Cells:Experimental
3;; Evidence and Modeling of a Slow K+-Dependent Mechanism
4;; Egidio D'Angelo,Thierry Nieus,Arianna Maffei,Simona Armano,Paola Rossi,Vanni Taglietti,
5;; Andrea Fontana and Giovanni Naldi
6
7(nemo-model Golgi
8
9  (
10   (input v celsius)
11
12   (const V_t = -35)
13
14   (component (type defaults)
15      (const celsius = 30)
16     
17      (output celsius)
18      )
19
20   (component (type membrane-capacitance)
21         (const C_m  = 1.0)
22         (output C_m  ))
23
24
25   (component (type decaying-pool) (name ca)
26
27      (input (ica from ion-currents))
28
29              (const  F       = 96485.3)
30              (const  d       = .2)
31              (const  cao     = 2.)
32              (const  cai0    = 50e-6)
33              (const  beta    = 1.3)
34
35              (d (ca) =  ((neg (ica) / (2 * F * d)) * 1e4 -
36                          (beta * (ca - cai0)))
37                         (initial cai0))
38             
39              (output ca cao)
40              )
41
42
43   (component (type decaying-pool) (name ca2)
44
45      (input (ica2 from ion-currents))
46
47      (const valence  = 2)
48      (const  F       = 96485.3)
49      (const  d       = .2)
50      (const  ca2o     = 2.)
51      (const  ca2i0    = 50e-6)
52      (const  beta    = 1.3)
53     
54      (d (ca2) =  ((neg (ica2) / (2 * F * d)) * 1e4 -
55                   (beta * (ca2 - ca2i0)))
56         (initial ca2i0))
57     
58      (output ca2 ca2o valence)
59      )
60
61   
62   (component (type ionic-current) (name CaHVA )
63             
64      (input
65       (cai from ion-pools)
66       (cao from ion-pools))
67             
68      (component (type gate)
69
70                 ;; rate constants
71                 
72                 (Q10 = (pow (3 ((celsius - 20) / 10))))
73                 
74                 (const Aalpha_s  = 0.04944)
75                 (const Kalpha_s  =  15.87301587302)
76                 (const V0alpha_s = -29.06)
77                 
78                 (const Abeta_s  = 0.08298)
79                 (const Kbeta_s  =  -25.641)
80                 (const V0beta_s = -18.66)
81                 
82                 (const Aalpha_u  = 0.0013)
83                 (const Kalpha_u  =  -18.183)
84                 (const V0alpha_u = -48)
85                 
86                 (const Abeta_u = 0.0013)
87                 (const Kbeta_u = 83.33)
88                 (const V0beta_u = -48)
89                 
90                 ;; rate functions
91                 
92                 (defun alpha_s (v Q10)
93                   (Q10 * Aalpha_s * exp((v - V0alpha_s) / Kalpha_s)))
94                 
95                 (defun beta_s (v Q10)
96                   (Q10 * Abeta_s * exp((v - V0beta_s) / Kbeta_s)))
97                 
98                 (defun alpha_u (v Q10)
99                   (Q10 * Aalpha_u * exp((v - V0alpha_u) / Kalpha_u)))
100                 
101                 (defun beta_u (v Q10)
102                   (Q10 * Abeta_u * exp((v - V0beta_u) / Kbeta_u)))
103
104                 
105                 (hh-ionic-gate
106                  (CaHVA  ;; ion name: exported variables will be of the form {ion}_{id}
107                   (initial-m  ((alpha_s (v Q10))/(alpha_s (v Q10) + beta_s (v Q10)) ))
108                   (initial-h  ((alpha_u (v Q10))/(alpha_u (v Q10) + beta_u (v Q10)) ))
109                   (m-power    2)
110                   (h-power    1)
111                   (m-inf      ((alpha_s (v Q10) / (alpha_s (v Q10) + beta_s (v Q10)) )))
112                   (m-tau      ((1 / (alpha_s (v Q10) + beta_s (v Q10)) )))
113                   (h-inf      ((alpha_u (v Q10) / (alpha_u (v Q10) + beta_u (v Q10)) )))
114                   (h-tau      ((1 / (alpha_u (v Q10) + beta_u (v Q10)) )))
115                   ))
116                 
117                 )
118             
119              (component (type pore)
120                         (const  gbar  = 460e-6)
121                         (output gbar ))
122             
123              (component (type permeating-ion) (name ca)
124                         (const F = 96485.3)
125                         (const R = 8.31342)
126                         (e = ((1e3) * (R * (celsius + 273.15)) / (2 * F) * ln (cao / cai)))
127                         (output e)
128                         )
129             
130              ) ;; end CaHVA current
131
132
133   (component (type voltage-clamp) (name CaHVA)
134             
135              (const vchold   = -71)
136              (const vcbase   = -69)
137              (const vcinc    = 10)
138              (const vcsteps  = 8)
139              (const vchdur   = 30)
140              (const vcbdur   = 100)
141             
142              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
143              )
144
145   
146   (component (type ionic-current) (name CaLVA )
147             
148      (input
149       (ca2i from ion-pools)
150       (ca2o from ion-pools))
151
152      (component (type gate)
153
154                 (const shift   = 2) ; screening charge for Ca_o = 2 mM
155                 
156                 (const v0_m_inf = -50)
157                 (const v0_h_inf = -78)
158                 (const k_m_inf  = -7.4)
159                 (const k_h_inf  = 5.0)
160       
161                 (const C_tau_m   = 3)
162                 (const A_tau_m   = 1.0)
163                 (const v0_tau_m1 = -25)
164                 (const v0_tau_m2 = -100)
165                 (const k_tau_m1  = 10)
166                 (const k_tau_m2 = -15)
167                 
168                 (const C_tau_h   = 85)
169                 (const A_tau_h   = 1.0)
170                 (const v0_tau_h1 = -46)
171                 (const v0_tau_h2 = -405)
172                 (const k_tau_h1  = 4)
173                 (const k_tau_h2  = -50)
174                 
175                         
176                 ;; rate functions
177                 
178                 (phi_m = (pow (5.0 ((celsius - 24) / 10))))
179                 (phi_h = (pow (3.0 ((celsius - 24) / 10))))
180                 
181                 (m_inf = (1.0 / ( 1 + exp ((v + shift - v0_m_inf) / k_m_inf)) ))
182                 (h_inf = (1.0 / ( 1 + exp ((v + shift - v0_h_inf) / k_h_inf)) ))
183                 
184                 (tau_m = ( (C_tau_m + A_tau_m / ( exp ((v + shift - v0_tau_m1) / k_tau_m1) +
185                                                       exp ((v + shift - v0_tau_m2) / k_tau_m2) ) ) / phi_m) )
186                 
187                 (tau_h = ( (C_tau_h + A_tau_h / ( exp ((v + shift - v0_tau_h1 ) / k_tau_h1) +
188                                                       exp ((v + shift - v0_tau_h2) / k_tau_h2) ) ) / phi_h) )
189                 
190
191                 (hh-ionic-gate
192                  (CaLVA  ;; ion name: exported variables will be of the form {ion}_{id}
193                   (initial-m  m_inf)
194                   (initial-h  h_inf)
195                   (m-power    2)
196                   (h-power    1)
197                   (m-inf      m_inf)
198                   (m-tau      tau_m)
199                   (h-inf      h_inf)
200                   (h-tau      tau_h)
201                   ))
202                 )
203     
204     
205      (component (type pore)
206                 (const  gbar  = 2.5e-4)
207                 (output gbar ))
208
209             
210      (component (type permeating-ion) (name ca2)
211                 (const valence = 2)
212                 (const F = 96485.3)
213                 (const R = 8.31342)
214                 (e = ((1e3) * (R * (celsius + 273.15)) / (2 * F) * ln (ca2o / ca2i)))
215                 (output e valence))
216
217             
218      ) ;; end CaLVA current
219
220
221
222   (component (type voltage-clamp) (name CaLVA)
223             
224              (const vchold   = -71)
225              (const vcbase   = -69)
226              (const vcinc    = 10)
227              (const vcsteps  = 8)
228              (const vchdur   = 200)
229              (const vcbdur   = 30)
230             
231              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
232              )
233
234
235     
236   (component (type ionic-current) (name HCN1)
237
238              (const Ehalf = -72.49)
239              (const c = 0.11305)
240             
241              (const rA = 0.002096)
242              (const rB = 0.97596)
243             
244              (defun r (potential)
245                ((rA * potential) + rB))
246             
247              (defun tau (potential t1 t2 t3)
248                (exp (((t1 * potential) - t2) * t3)))
249             
250              (defun o_inf (potential Ehalf c)
251                (1 / (1 + exp ((potential - Ehalf) * c))))
252             
253
254              (component (type gate)
255                         
256
257                         ;; rate constants
258       
259                         (const tCs = 0.01451)
260                         (const tDs = -4.056)
261                         (const tEs = 2.302585092)
262
263                         (o_slow_inf = ((1 - r (v)) * o_inf (v Ehalf c)))
264                       
265                         (tau_s =  (tau (v tCs tDs tEs)) )
266
267                         (d (o_slow) = ((o_slow_inf - o_slow) / tau_s)
268                            (initial o_slow_inf))
269                         
270                         (output o_slow)
271
272                         )
273             
274              (component (type gate)
275                         
276
277                         ;; rate constants
278       
279                         (const tCf = 0.01371)
280                         (const tDf = -3.368)
281                         (const tEf = 2.302585092)
282
283                         (o_fast_inf = (r (v) * o_inf (v Ehalf c)))
284                       
285                         (tau_f =  (tau (v tCf tDf tEf)) )
286
287                         (d (o_fast) = ((o_fast_inf - o_fast) / tau_f)
288                            (initial o_fast_inf))
289                         
290                         (output o_fast)
291
292                         )
293             
294              (component (type pore)
295                         (const  gbar  = 5e-5)
296                         (output gbar))
297             
298              (component (type permeating-ion) (name non-specific)
299                         (const e = -20)
300                         (output e ))
301             
302              ) ;; end HCN1 current
303
304
305             
306   (component (type voltage-clamp) (name HCN1)
307             
308              (const vchold   = -71)
309              (const vcbase   = -69)
310              (const vcinc    = 10)
311              (const vcsteps  = 8)
312              (const vchdur   = 30)
313              (const vcbdur   = 100)
314             
315              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
316              )
317
318
319     
320   (component (type ionic-current) (name HCN2)
321                         
322
323              (const Ehalf = -81.95)
324              (const c = 0.1661)
325             
326              ;; rate constants
327
328              (const rA = -0.0227)
329              (const rB = -1.4694)
330
331              (defun r (potential r1 r2)
332                (if (potential >= -64.70)
333                    then  0.0
334                    else (if (potential <= -108.70)
335                             then 1.0
336                             else (r1 * potential) + r2)))
337
338              (defun o_inf (potential Ehalf c)
339                (1 / (1 + exp((potential - Ehalf) * c))))
340             
341              (component (type gate)
342                         
343
344                         ;; rate constants
345
346                         (const tCs = 0.0152)
347                         (const tDs = -5.2944)
348                         (const tEs = 2.3026)
349                         
350                         (defun tau_slow (potential t1 t2 t3)
351                           (exp (t3 * ((t1 * potential) - t2))))
352
353                         (o_slow_inf = ((1 - r (v rA rB)) * o_inf (v Ehalf c)))
354
355                         (tau_s =  (tau_slow(v tCs tDs tEs)))
356
357                         (d (o_slow) = ((o_slow_inf - o_slow) / tau_s)
358                            (initial o_slow_inf))
359
360                         (output o_slow)
361
362                         )
363             
364              (component (type gate)
365                         
366
367                         ;; rate constants
368
369                         (const tCf = 0.0269)
370                         (const tDf = -5.6111)
371                         (const tEf = 2.3026)
372                         
373                         (defun tau_fast (potential t1 t2 t3)
374                           (exp (t3 * ((t1 * potential) - t2))))
375
376                         (o_fast_inf = (r (v rA rB) * o_inf (v Ehalf c)))
377
378                         (tau_f =  (tau_fast(v tCf tDf tEf)))
379
380                         (d (o_fast) = ((o_fast_inf - o_fast) / tau_f)
381                            (initial o_fast_inf))
382
383                         (output o_fast)
384
385                         )
386             
387              (component (type pore)
388                         (const  gbar  = 8e-5)
389                         (output gbar))
390             
391              (component (type permeating-ion) (name non-specific)
392                         (const e = -20)
393                         (output e ))
394             
395              ) ;; end HCN2 current
396
397
398
399   (component (type voltage-clamp) (name HCN2)
400             
401              (const vchold   = -71)
402              (const vcbase   = -69)
403              (const vcinc    = 10)
404              (const vcsteps  = 8)
405              (const vchdur   = 30)
406              (const vcbdur   = 100)
407             
408              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
409              )
410
411   
412
413   (component (type ionic-current) (name KA )
414
415              (defun sigm (x y) (1 / (exp (x / y) + 1)))
416             
417              (component (type gate)
418                         
419                         ;; rate constants
420                         
421                         (Q10 = (pow (3 ((celsius - 25.5) / 10))))
422
423                         (const Aalpha_a  = 0.8147)
424                         (const Kalpha_a  = -23.32708)
425                         (const V0alpha_a = -9.17203)
426                         (const Abeta_a   = 0.1655)
427                         (const Kbeta_a   = 19.47175)
428                         (const V0beta_a  = -18.27914)
429
430                         (const Aalpha_b  = 0.0368)
431                         (const Kalpha_b  = 12.8433)
432                         (const V0alpha_b = -111.33209)
433                         (const Abeta_b   = 0.0345)
434                         (const Kbeta_b   = -8.90123)
435                         (const V0beta_b  = -49.9537)
436
437                         (const V0_ainf = -38)
438                         (const  K_ainf = -17)
439
440                         (const V0_binf   = -78.8)
441                         (const K_binf    = 8.4)
442
443                         
444                         ;; rate functions
445
446                         (defun alpha_a (v Q10)
447                           (Q10 * Aalpha_a * sigm((v - V0alpha_a) Kalpha_a)))
448
449                         (defun beta_a (v Q10)
450                           (Q10 * (Abeta_a / exp((v - V0beta_a) / Kbeta_a))))
451
452                         (defun alpha_b (v Q10)
453                           (Q10 * Aalpha_b * sigm((v - V0alpha_b) Kalpha_b)))
454
455                         (defun beta_b (v Q10)
456                           (Q10 * Abeta_b * sigm((v - V0beta_b) Kbeta_b)))
457
458
459
460                         (a_inf = (1 / (1 + exp ((v - V0_ainf) / K_ainf))))
461                         (tau_a = (1 / (alpha_a (v Q10) + beta_a (v Q10)) ))
462                         (b_inf = (1 / (1 + exp ((v - V0_binf) / K_binf))))
463                         (tau_b = (1 / (alpha_b (v Q10) + beta_b (v Q10)) ))
464
465
466                         (hh-ionic-gate
467                          (KA  ;; ion name: exported variables will be of the form {ion}_{id}
468                           (initial-m  (a_inf))
469                           (initial-h  (b_inf))
470                           (m-power    3)
471                           (h-power    1)
472                           (m-inf      a_inf)
473                           (m-tau      tau_a)
474                           (h-inf      b_inf)
475                           (h-tau      tau_b)
476                           ))
477
478                         )
479
480                     
481              (component (type pore)
482                         (const  gbar  = 0.008)
483                         (output gbar ))
484
485             
486              (component (type permeating-ion) (name k)
487                         (const e = -84.69)
488                         (output e ))
489             
490              ) ;; end KA current
491             
492
493
494   (component (type voltage-clamp) (name KA)
495             
496              (const vchold   = -71)
497              (const vcbase   = -69)
498              (const vcinc    = 10)
499              (const vcsteps  = 8)
500              (const vchdur   = 30)
501              (const vcbdur   = 100)
502             
503              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
504              )
505
506   
507
508   (component (type ionic-current) (name KCa )
509
510      (input
511       (cai from ion-pools))
512             
513      (component (type gate)
514
515                 ;; rate constants
516                 (Q10 = (pow (3 ((celsius - 30) / 10))))
517                 
518                 (const Aalpha_c = 7)
519                 (const Balpha_c = 1.5e-3)
520                 
521                 (const Kalpha_c =  -11.765)
522                 
523                 (const Abeta_c = 1.)
524                 (const Bbeta_c = 0.15e-3)
525
526                 (const Kbeta_c = -11.765)
527                 
528                 ;; rate functions
529                 (defun alpha_c (v cai Q10)
530                   (Q10 * Aalpha_c / (1 + (Balpha_c * exp(v / Kalpha_c) / cai))))
531                 
532                 (defun beta_c (v cai Q10)
533                   (Q10 * Abeta_c / (1 + (cai / (Bbeta_c * exp (v / Kbeta_c))) )))
534
535
536                 (hh-ionic-gate
537                  (KCa  ;; ion name: exported variables will be of the form {ion}_{id}
538                   (initial-m  ((alpha_c (v cai Q10)) / (alpha_c (v cai Q10) + beta_c (v cai Q10)) ))
539                   (m-power    1)
540                   (h-power    0)
541                   (m-alpha    (alpha_c (v cai Q10)))
542                   (m-beta     (beta_c (v cai Q10)))
543                   ))
544                 
545                 )
546     
547      (component (type pore)
548                 (const  gbar  = 0.003)
549                 (output gbar ))
550     
551     
552      (component (type permeating-ion) (name k)
553                 (const e = -84.69)
554                 (output e ))
555     
556      (component (type modulating-ion) (name ca) )
557     
558      ) ;; end KCa current
559   
560   (component (type voltage-clamp) (name KCa)
561             
562              (const vchold   = -71)
563              (const vcbase   = -69)
564              (const vcinc    = 10)
565              (const vcsteps  = 8)
566              (const vchdur   = 30)
567              (const vcbdur   = 100)
568             
569              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
570              )
571
572
573   (component (type ionic-current) (name KM )
574             
575              (component (type gate)
576
577                         ;; rate constants
578                         (const Aalpha_n = 0.0033)
579
580                         (const Kalpha_n  = 40)
581                         (const V0alpha_n = -30)
582                         (const Abeta_n   = 0.0033)
583
584                         (const Kbeta_n  = -20)
585                         (const V0beta_n = -30)
586                         (const V0_ninf  = -35)
587                         (const   B_ninf = 6)
588                         
589                         (Q10 = (pow (3 ((celsius - 22) / 10))))
590                         
591                         ;; rate equations
592                         (defun alpha_n (v Q10)
593                           (Q10 * Aalpha_n * exp((v - V0alpha_n) / Kalpha_n) ))
594
595                         (defun beta_n (v Q10)
596                           (Q10 * Abeta_n * exp((v - V0beta_n) / Kbeta_n) ))
597
598                         (hh-ionic-gate
599                          (KM  ;; ion name: exported variables will be of the form {ion}_{id}
600                           (initial-m  ((1 / (1 + exp((neg (v - V0_ninf)) / B_ninf)))))
601                           (m-power    1)
602                           (h-power    0)
603                           (m-tau      ((1 / (alpha_n(v Q10) + beta_n (v Q10)) )))
604                           (m-inf      ((1 / (1 + exp((neg (v - V0_ninf)) / B_ninf)))))
605                           ))
606                         )
607             
608              (component (type pore)
609                         (const  gbar  = 0.001)
610                         (output gbar ))
611             
612              (component (type permeating-ion) (name k)
613                         (const e = -84.69)
614                         (output e ))
615             
616              ) ;; end KM current
617
618
619
620   (component (type voltage-clamp) (name KM)
621             
622              (const vchold   = -71)
623              (const vcbase   = -69)
624              (const vcinc    = 10)
625              (const vcsteps  = 8)
626              (const vchdur   = 30)
627              (const vcbdur   = 100)
628             
629              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
630              )
631
632   
633
634   (component (type ionic-current) (name KV )
635             
636              (defun linoid (x y)
637                (if ((abs (x / y)) < 1e-6)
638                    then (y * (1 - ((x / y) / 2)))
639                    else (x / (exp (x / y) - 1))
640                    ))
641
642
643              (component (type gate)
644
645                         ;; rate constants
646
647                         (const Aalpha_n = -0.01)
648                         (const Kalpha_n = -10)
649                         (const V0alpha_n = -26)
650                         (const Abeta_n = 0.125)
651       
652                         (const Kbeta_n = -80)
653                         (const V0beta_n = -36)
654
655                         ;; rate functions
656                         (Q10 = (pow (3 ((celsius - 6.3) / 10))))
657
658                         (defun alpha_n (v Q10)
659                           (Q10 * Aalpha_n * linoid ((v - V0alpha_n) Kalpha_n)))
660
661                         (defun beta_n (v Q10)
662                           (Q10 * Abeta_n * exp((v - V0beta_n) / Kbeta_n) ))
663
664                         (hh-ionic-gate
665                          (KV  ;; ion name: exported variables will be of the form {ion}_{id}
666                           (initial-m  ((alpha_n (v Q10)) / (alpha_n (v Q10) + beta_n (v Q10)) ))
667                           (m-power    4)
668                           (h-power    0)
669                           (m-alpha    (alpha_n(v Q10)))
670                           (m-beta     (beta_n(v Q10)))
671                           ))
672                         )
673
674              (component (type pore)
675                         (const  gbar  = 0.032)
676                         (output gbar ))
677             
678              (component (type permeating-ion) (name k)
679                         (const e = -84.69)
680                         (output e ))
681             
682              ) ;; end KV current
683             
684
685   (component (type voltage-clamp) (name KV)
686             
687              (const vchold   = -71)
688              (const vcbase   = -69)
689              (const vcinc    = 10)
690              (const vcsteps  = 8)
691              (const vchdur   = 30)
692              (const vcbdur   = 100)
693             
694              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
695              )
696
697
698
699   
700   (component (type ionic-current) (name SK2 )
701
702      (input
703       (cai from ion-pools))
704             
705      (component (type gate)
706
707                 (Q10 = (pow (3 ((celsius - 23) / 10))))
708                 
709                 (const diff = 3) ;; Diffusion factor
710                 
711                 ;; rates ca-independent
712                 (const invc1 = 80e-3)
713                 (const invc2 = 80e-3)
714                 (const invc3 = 200e-3)
715
716                 (const invo1 = 1)
717                 (const invo2 = 100e-3)
718                 (const diro1 = 160e-3)
719                 (const diro2 = 1.2)
720
721                 ;; rates ca-dependent
722                 (const dirc2 = 200)
723                 (const dirc3 = 160)
724                 (const dirc4 = 80)
725
726                 (invc1_t = (invc1 * Q10  ))
727                 (invc2_t = (invc2 * Q10 ))
728                 (invc3_t = (invc3 * Q10 ))
729                 (invo1_t = (invo1 * Q10 ))
730                 (invo2_t = (invo2 * Q10 ))
731                 (diro1_t = (diro1 * Q10 ))
732                 (diro2_t = (diro2 * Q10 ))
733                 (dirc2_t = (dirc2 * Q10 ))
734                 (dirc3_t = (dirc3 * Q10 ))
735                 (dirc4_t = (dirc4 * Q10 ))
736
737
738                 (dirc2_t_ca = (dirc2_t * (cai / diff) ))
739                 (dirc3_t_ca = (dirc3_t * (cai / diff) ))
740                 (dirc4_t_ca = (dirc4_t * (cai / diff) ))
741
742                 
743                 (reaction
744                  (SK2_z
745                   (transitions
746                    (<-> c1 c2 dirc2_t_ca invc1_t)
747                    (<-> c2 c3 dirc3_t_ca invc2_t)
748                    (<-> c3 c4 dirc4_t_ca invc3_t)
749                    (<-> c3 o1 diro1_t invo1_t)
750                    (<-> c4 o2 diro2_t invo2_t)
751                    )
752                   (conserve  ((1 = (c1 + c2 + c3 + c4 + o2 + o1))))
753                   (open o1 o2)
754                   (power 1)))
755                 
756                 (output SK2_z)
757                 
758                 )
759
760      (component (type pore)
761                 (const  gbar  = 0.038)
762                 (output gbar ))
763     
764     
765      (component (type permeating-ion) (name k)
766                 (const e = -84.69)
767                 (output e ))
768
769      (component (type modulating-ion) (name ca) )
770     
771      ) ;; end SK2 current
772
773
774
775   (component (type voltage-clamp) (name SK2)
776           
777              (const vchold   = -71)
778              (const vcbase   = -69)
779              (const vcinc    = 10)
780              (const vcsteps  = 8)
781              (const vchdur   = 30)
782              (const vcbdur   = 100)
783             
784              (output vchold vcbase vcinc vcsteps vchdur vcbdur)
785              )
786   
787   (component (type ionic-current) (name Na )
788             
789             
790              (defun linoid (x y)
791                (if ((abs (x / y)) < 1e-6)
792                    then (y * (1 - ((x / y) / 2)))
793                    else (x / (1 - exp (x / y) ))
794                    ))
795             
796             
797              (component (type gate)
798             
799              ;; rate constants
800                         
801                         (Q10 = (pow (3 ((celsius - 20) / 10))))
802
803                         (const Aalpha_u  = 0.3)
804                         (const Kalpha_u  = -10)
805                         (const V0alpha_u = -25)
806       
807                         (const Abeta_u  = 12)
808                         (const Kbeta_u  = -18.182)
809                         (const V0beta_u = -50)
810
811                         (const Aalpha_v  = 0.21)
812                         (const Kalpha_v  = -3.333)
813                         (const V0alpha_v = -50)
814 
815                         (const Abeta_v  = 3)
816                         (const Kbeta_v  = -5)
817                         (const V0beta_v = -17)
818
819                         ;; rate functions
820                         (defun alpha_u (v Q10)
821                           (Q10 * Aalpha_u * linoid((v - V0alpha_u) Kalpha_u) ))
822
823                         (defun beta_u (v Q10)
824                           (Q10 * Abeta_u * exp((v - V0beta_u) / Kbeta_u) ))
825
826                         (defun alpha_v (v Q10)
827                           (Q10 * Aalpha_v * exp((v - V0alpha_v) / Kalpha_v) ))
828                               
829                         (defun beta_v (v Q10)
830                           (Q10 * Abeta_v / (1 + exp((v - V0beta_v) / Kbeta_v) )))
831
832                         
833                         (hh-ionic-gate
834                          (Na  ;; ion name: exported variables will be of the form {ion}_{id}
835                           (initial-m  ((alpha_u (v Q10)) / (alpha_u (v Q10) + beta_u (v Q10)) ))
836                           (initial-h  ((alpha_v (v Q10)) / (alpha_v (v Q10) + beta_v (v Q10)) ))
837                           (m-power    3)
838                           (h-power    1)
839                           (m-alpha    (alpha_u (v Q10)))
840                           (m-beta     (beta_u (v Q10)))
841                           (h-alpha    (alpha_v (v Q10)))
842                           (h-beta     (beta_v (v Q10)))
843                           ))
844                         
845                         )
846             
847              (component (type pore)
848                         (const  gbar  = 0.048)
849                         (output gbar ))
850
851             
852              (component (type permeating-ion) (name na)
853                         (const e = 87.39)
854                         (output e ))
855             
856              ) ;; end Na current
857
858
859   (component (type voltage-clamp) (name Na)
860
861              (const vchold   = -71)
862              (const vcbase   = -60)
863              (const vcinc    = 10)
864              (const vcsteps  = 9)
865              (const vchdur   = 30)
866              (const vcbdur   = 100)
867           
868              (output vchold vcbase vcinc vcsteps vchdur vcbdur))
869
870   
871   (component (type ionic-current) (name NaP )
872             
873              (defun linoid (x y)
874                (if ((abs (x / y)) < 1e-6)
875                    then (y * (1 - ((x / y) / 2)))
876                    else (x / (exp (x / y) - 1))
877                    ))
878
879
880              (component (type gate)
881
882                         ;; rate constants
883                         ( Q10 = (pow (3 ((celsius - 30) / 10))))
884
885                         (const Aalpha_m  = -0.91)
886                         (const Kalpha_m  = -5)
887                         (const V0alpha_m = -40)
888                         (const Abeta_m   = 0.62)
889                         (const Kbeta_m   = 5)
890                         (const V0beta_m  = -40)
891                         (const V0_minf   = -43)
892                         (const B_minf    = 5)
893
894                         ;; rate functions
895
896                         (defun alpha_m (v Q10)
897                           (Q10 * Aalpha_m * linoid( ( v - V0alpha_m ) Kalpha_m) ))
898
899                         (defun beta_m (v Q10)
900                           (Q10 * Abeta_m * linoid( ( v - V0beta_m )  Kbeta_m ) ))
901
902                         (hh-ionic-gate
903                          (NaP  ;; ion name: exported variables will be of the form {ion}_{id}
904                           (initial-m  ((1 / (1 + exp ((neg (v - V0_minf)) / B_minf)))))
905                           (m-power    1)
906                           (h-power    0)
907                           (m-inf     ((1 / (1 + exp ((neg (v - V0_minf)) / B_minf)))))
908                           (m-tau     ((5 / (alpha_m (v Q10) + beta_m (v Q10)) )))
909                           ))
910                         )
911
912              (component (type pore)
913                         (const  gbar  = 0.00019)
914                         (output gbar ))
915
916              (component (type permeating-ion) (name na)
917                         (const e = 87.39)
918                         (output e ))
919             
920              ) ;; end NaP current
921
922
923   (component (type voltage-clamp) (name NaP)
924
925              (const vchold   = -71)
926              (const vcbase   = -60)
927              (const vcinc    = 10)
928              (const vcsteps  = 9)
929              (const vchdur   = 30)
930              (const vcbdur   = 100)
931           
932              (output vchold vcbase vcinc vcsteps vchdur vcbdur))
933
934   
935   (component (type ionic-current) (name NaR )
936             
937              (component (type gate)
938
939                         ;; rate constants
940                         (Q10 = (pow (3 ((celsius - 20) / 10))))
941
942                         (const Aalpha_s     = -0.00493)
943                         (const V0alpha_s    = -4.48754)
944                         (const Kalpha_s     = -6.81881)
945                         (const Shiftalpha_s = 0.00008)
946
947                         (const Abeta_s     = 0.01558)
948                         (const V0beta_s    = 43.97494)
949                         (const Kbeta_s     =  0.10818)
950                         (const Shiftbeta_s = 0.04752)
951
952                         (const Aalpha_f  = 0.31836)
953                         (const V0alpha_f = -80)
954                         (const Kalpha_f  = -62.52621)
955
956                         (const Abeta_f  = 0.01014)
957                         (const V0beta_f = -83.3332)
958                         (const Kbeta_f  = 16.05379)
959
960                         ;; rate functions
961                         (defun alpha_s (v Q10)
962                           (Q10 * (Shiftalpha_s + ((Aalpha_s * ((v + V0alpha_s) )) / (exp ((v + V0alpha_s) / Kalpha_s) - 1)))))
963
964                         (defun beta_s (v Q10)
965                           (let ((x1 ((v + V0beta_s) / Kbeta_s)))
966                             (Q10 * (Shiftbeta_s + Abeta_s * ((v + V0beta_s) / (exp((if (x1 > 200) then 200 else x1)) - 1))))))
967
968                         (defun alpha_f (v Q10)
969                           (Q10 * Aalpha_f * exp( ( v - V0alpha_f ) / Kalpha_f)))
970
971                         (defun beta_f (v Q10)
972                           (Q10 * Abeta_f * exp( ( v - V0beta_f ) / Kbeta_f )  ))
973
974                         (hh-ionic-gate
975                          (NaR  ;; ion name: exported variables will be of the form {ion}_{id}
976                           (initial-m  ((alpha_s (v Q10)) / (alpha_s (v Q10) + beta_s (v Q10)) ))
977                           (initial-h  ((alpha_f (v Q10)) / (alpha_f (v Q10) + beta_f (v Q10)) ))
978                           (m-power    1)
979                           (h-power    1)
980                           (m-alpha    (alpha_s (v Q10)))
981                           (m-beta     (beta_s (v Q10)))
982                           (h-alpha    (alpha_f (v Q10)))
983                           (h-beta     (beta_f (v Q10)))
984
985                           ))
986                         )
987
988              (component (type pore)
989                         (const  gbar  = 0.0017)
990                         (output gbar ))
991
992              (component (type permeating-ion) (name na)
993                         (const e = 87.39)
994                         (output e ))
995             
996              ) ;; end Nar current
997
998
999   (component (type voltage-clamp) (name NaR)
1000
1001              (const vchold   = -71)
1002              (const vcbase   = -60)
1003              (const vcinc    = 10)
1004              (const vcsteps  = 9)
1005              (const vchdur   = 30)
1006              (const vcbdur   = 100)
1007           
1008              (output vchold vcbase vcinc vcsteps vchdur vcbdur))
1009
1010     
1011   (component (type ionic-current) (name Lkg)
1012             
1013              (component (type pore)
1014                         (const  gbar  = (21e-6))
1015                         (output gbar))
1016             
1017              (component (type permeating-ion) (name non-specific)
1018                         (const e = -55)
1019                         (output e ))
1020             
1021              ) ;; end leak current
1022   )
1023
1024   ;; Following are templates for various driver scripts used to run this model
1025   (
1026    ("template.hoc" () 
1027#<<EOF
1028/*******Cerebellar Golgi Cell Model **********
1029
1030Developers:    Sergio Solinas & Egidio D'Angelo
1031Code contributors:  Thierry Neius, Shyam Diwakar, Lia Forti
1032Data Analysis: Sergio Solinas
1033
1034Work Progress: April 2004 - May 2007
1035
1036Developed At:  Università Degli Studi Di Pavia
1037               Dipartimento Di Scienze Fisiologiche
1038               Pavia - Italia
1039               
1040Model Published in:
1041             Sergio M. Solinas, Lia Forti, Elisabetta Cesana,
1042             Jonathan Mapelli, Erik De Schutter and Egidio D`Angelo (2008)
1043             Computational reconstruction of pacemaking and intrinsic
1044             electroresponsiveness in cerebellar golgi cells
1045             Frontiers in Cellular Neuroscience 2:2
1046
1047
1048********************************************/
1049
1050// In this configuration the ion channels
1051// were not corrected for the Liquid Junction potential.
1052// The ion reversal potential were corrected in agreement
1053// with the voltage shift.
1054// See Table 1 Solinas et al. 2008 Frontiers Neurosci 2:2
1055
1056begintemplate GoC
1057public soma,axon,elec,seal,dend
1058public SpikeTrain, RT_Vm, EL_Vm
1059
1060create soma
1061create axon
1062create elec,seal
1063create dend[3]
1064
1065objref SpikeTrain, netcon, nil
1066objref RT_Vm, EL_Vm
1067
1068proc init() {
1069
1070    RT_Vm = new Vector()
1071    EL_Vm = new Vector()
1072    create soma
1073    soma {
1074        nseg = 1
1075        diam = 27 // 22 pF Dieudonne98
1076        L = 27
1077        Ra = 100 // From Roth&Hausser2000
1078        celsius = 23
1079       
1080        insert Golgi
1081       
1082        cai0_ca_ion = 50e-6
1083        ca2i0_ca2_ion = cai0_ca_ion
1084       
1085        cai0_Golgi_Ca = cai0_ca_ion
1086        ca2i0_Golgi_Ca2 = cai0_ca_ion
1087       
1088        cai = cai0_ca_ion
1089        ca2i = cai
1090        ca2o = cao
1091       
1092        ena =  87.39
1093        ek  = -84.69
1094
1095        SpikeTrain = new Vector()
1096        netcon = new NetCon(&v(0.5),nil)
1097        netcon.threshold=-20
1098        netcon.record(SpikeTrain)
1099       
1100        RT_Vm.record(&v(0.5))
1101    }
1102   
1103    create dend[3]
1104    for i=0,2 {
1105        dend[i] {
1106            nseg = 10
1107            diam = 3
1108            L = 113
1109            Ra = 100
1110            celsius = 23
1111           
1112            insert Golgi_Lkg
1113        }
1114        connect dend[i](0), soma(1)     
1115    }
1116   
1117   
1118    create axon
1119    axon {
1120        nseg = 100
1121        diam = 2.4 // gives 90 pF to get to the 145 pF Forti06
1122        L = 1200
1123        Ra = 100
1124        celsius = 23
1125       
1126        insert Golgi_Lkg
1127    }
1128   
1129    connect axon(0), soma(0)
1130   
1131    create elec,seal
1132    elec {
1133        nseg = 1
1134        diam = 3
1135        L = 1000
1136        Ra = 36
1137        cm = 0.0015
1138        celsius = 23
1139        EL_Vm.record(&v(0.5))
1140
1141    }
1142       
1143    seal {
1144        nseg = 1
1145        diam = 3
1146        L = 1
1147        Ra = 1
1148        cm = 0.0001
1149        celsius = 23
1150       
1151    }
1152       
1153    connect seal(1), soma(1)
1154    connect elec(1), seal(0)
1155
1156}
1157
1158                   
1159endtemplate GoC
1160EOF
1161)
1162    ("protocol.hoc" () 
1163#<<EOF
1164load_file("nrngui.hoc")
1165load_file("Golgi_protocol.ses")
1166
1167// Load the Golgi cell template
1168xopen("Golgi_template.hoc")
1169objref goc
1170goc = new GoC()
1171
1172
1173prelength = 1000
1174mainlength = 2000
1175
1176
1177//**********************************************************************
1178proc simulate() { local preDT, mainDT, logsize  localobj logfile, tlog, Vlog, iklog, inalog, icalog, ica2log
1179   
1180    mainDT = 0.025
1181    preDT = 0.025
1182   
1183    dt = preDT
1184    tstop = prelength
1185    run()
1186   
1187   
1188    if ( stoprun == 1) return
1189   
1190    dt = mainDT
1191    continuerun(prelength + mainlength)
1192    if ( stoprun == 1) return
1193   
1194    logfile=$o1
1195    tlog=$o2
1196    Vlog=$o3
1197    iklog=$o4
1198    inalog=$o5
1199    icalog=$o6
1200    ica2log=$o7
1201   
1202    logsize = tlog.size()
1203   
1204    for i=0,tlog.size()-1 {
1205        logfile.printf("%g %g %g %g %g\n", tlog.x[i], Vlog.x[i], iklog.x[i], inalog.x[i], icalog.x[i], ica2log.x[i])
1206    }
1207
1208}
1209
1210
1211//*************User-Interface*******************************************
1212
1213nrnsecmenu(0.5, 1)
1214
1215xpanel("Spontaneous firing")
1216xvalue("Time for Initialization", "prelength")
1217xvalue("Main duration", "mainlength")
1218
1219xvalue("dt", "dt")
1220xvalue("t", "t")
1221xlabel("")
1222xbutton("Start", "simulate()")
1223xbutton("Stop", "stoprun = 1")
1224xpanel()
1225
1226vec_sizes = (prelength+mainlength)/dt + 1       // recorded traces are all this size
1227
1228
1229objref  iklog
1230iklog = new Vector(vec_sizes)
1231iklog.record (&goc.soma.ik(0.5))
1232
1233objref  inalog
1234inalog = new Vector(vec_sizes)
1235inalog.record (&goc.soma.ina(0.5))
1236
1237objref  icalog
1238icalog = new Vector(vec_sizes)
1239icalog.record (&goc.soma.ica(0.5))
1240
1241objref  ica2log
1242ica2log = new Vector(vec_sizes)
1243ica2log.record (&goc.soma.ica2(0.5))
1244
1245objref  Vlog
1246Vlog = new Vector(vec_sizes)
1247Vlog.record (&goc.soma.v(0.5))
1248
1249objref tlog
1250tlog = new Vector(vec_sizes,0)
1251tlog.record (&t)
1252
1253objref logfile
1254logfile = new File()
1255logfile.wopen ( "OFF_ON_OFF.dat" )
1256
1257simulate(logfile,tlog,Vlog,iklog,inalog,icalog,ica2log)
1258logfile.close()
1259
1260
1261quit()
1262EOF
1263)
1264    ("protocol.ses" () 
1265#<<EOF
1266{load_file("nrngui.hoc")}
1267objectvar save_window_, rvp_
1268objectvar scene_vector_[4]
1269objectvar ocbox_, ocbox_list_, scene_, scene_list_
1270{ocbox_list_ = new List()  scene_list_ = new List()}
1271{pwman_place(0,0,0)}
1272{
1273save_window_ = new Graph(0)
1274save_window_.size(0,5500,-80,50)
1275scene_vector_[2] = save_window_
1276{save_window_.view(0, -80, 5500, 130, 477, 9, 774.9, 232.3)}
1277graphList[0].append(save_window_)
1278save_window_.save_name("graphList[0].")
1279save_window_.addexpr("v(.5)", 1, 1, 0.8, 0.9, 2)
1280}
1281{doNotify()}
1282EOF
1283)
1284
1285    (".hoc" () 
1286#<<EOF
1287
1288create soma
1289access soma
1290
1291insert {{model}}
1292
1293soma {
1294        nseg = 1
1295
1296{% with diam = default(diam, 9.76), L = default(L, 9.76) %}
1297        diam = {{diam}}
1298        L = {{L}}
1299{% endwith %}
1300
1301        cm = 1
1302
1303        Ra = 100
1304    }
1305EOF
1306)
1307
1308    ("plot_neuron_original_iclamp.ses" () 
1309#<<EOF
1310
1311CaHVA0       = load ("original/CaHVA.dat");
1312CaHVA1       = load ("nemo_nmodl_generated/CaHVA.dat");
1313
1314CaLVA0       = load ("original/CaLVA.dat");
1315CaLVA1       = load ("nemo_nmodl_generated/CaLVA.dat");
1316
1317HCN10       = load ("original/HCN1.dat");
1318HCN11       = load ("nemo_nmodl_generated/HCN1.dat");
1319
1320HCN20       = load ("original/HCN2.dat");
1321HCN21       = load ("nemo_nmodl_generated/HCN2.dat");
1322
1323KA0       = load ("original/KA.dat");
1324KA1       = load ("nemo_nmodl_generated/KA.dat");
1325
1326KCa0       = load ("original/KCa.dat");
1327KCa1       = load ("nemo_nmodl_generated/KCa.dat");
1328
1329KM0      = load ("original/KM.dat");
1330KM1      = load ("nemo_nmodl_generated/KM.dat");
1331
1332KV0      = load ("original/KV.dat");
1333KV1      = load ("nemo_nmodl_generated/KV.dat");
1334
1335SK20      = load ("original/SK2.dat");
1336SK21      = load ("nemo_nmodl_generated/SK2.dat");
1337
1338Na0      = load ("original/Na.dat");
1339Na1      = load ("nemo_nmodl_generated/Na.dat");
1340
1341NaR0      = load ("original/NaR.dat");
1342NaR1      = load ("nemo_nmodl_generated/NaR.dat");
1343
1344NaP0      = load ("original/NaP.dat");
1345NaP1      = load ("nemo_nmodl_generated/NaP.dat");
1346
1347
1348subplot(3,4,1);
1349plot(CaHVA0(:,1),CaHVA0(:,2),CaHVA1(:,1),CaHVA1(:,2),'linewidth',2);
1350title ("CaHVA current");
1351
1352subplot(3,4,2);
1353plot(CaLVA0(:,1),CaLVA0(:,2),CaLVA1(:,1),CaLVA1(:,2),'linewidth',2);
1354title ("CaLVA current");
1355
1356subplot(3,4,3);
1357plot(HCN10(:,1),HCN10(:,2),HCN11(:,1),HCN11(:,2),'linewidth',2);
1358title ("HCN1 current");
1359
1360subplot(3,4,4);
1361plot(HCN20(:,1),HCN20(:,2),HCN21(:,1),HCN21(:,2),'linewidth',2);
1362title ("HCN2 current");
1363
1364subplot(3,4,5);
1365plot(KA0(:,1),KA0(:,2),KA1(:,1),KA1(:,2),'linewidth',2);
1366title ("KA current");
1367
1368subplot(3,4,6);
1369plot(KCa0(:,1),KCa0(:,2),KCa1(:,1),KCa1(:,2),'linewidth',2);
1370title ("KCa current");
1371
1372subplot(3,4,7);
1373plot(KM0(:,1),KM0(:,2),KM1(:,1),KM1(:,2),'linewidth',2);
1374title ("KM current");
1375
1376subplot(3,4,8);
1377plot(KV0(:,1),KV0(:,2),KV1(:,1),KV1(:,2),'linewidth',2);
1378title ("KV current");
1379
1380subplot(3,4,9);
1381plot(SK20(:,1),SK20(:,2),SK21(:,1),SK21(:,2),'linewidth',2);
1382title ("SK2 current");
1383
1384subplot(3,4,10);
1385plot(Na0(:,1),Na0(:,2),Na1(:,1),Na1(:,2),'linewidth',2);
1386title ("Na current");
1387
1388subplot(3,4,11);
1389plot(NaR0(:,1),NaR0(:,2),NaR1(:,1),NaR1(:,2),'linewidth',2);
1390title ("NaR current");
1391
1392subplot(3,4,12);
1393plot(NaP0(:,1),NaP0(:,2),NaP1(:,1),NaP1(:,2),'linewidth',2);
1394title ("NaP current");
1395
1396print  ("NEURON_Original_Iclamp.eps", "-depsc");
1397EOF
1398)
1399
1400
1401
1402    ("plot_neuron_original_vclamp.ses" () 
1403#<<EOF
1404
1405CaHVA0       = load ("original/CaHVA.dat");
1406CaHVA1       = load ("nemo_nmodl_generated/CaHVA.dat");
1407
1408CaLVA0       = load ("original/CaLVA.dat");
1409CaLVA1       = load ("nemo_nmodl_generated/CaLVA.dat");
1410
1411HCN10       = load ("original/HCN1.dat");
1412HCN11       = load ("nemo_nmodl_generated/HCN1.dat");
1413
1414HCN20       = load ("original/HCN2.dat");
1415HCN21       = load ("nemo_nmodl_generated/HCN2.dat");
1416
1417KA0       = load ("original/KA.dat");
1418KA1       = load ("nemo_nmodl_generated/KA.dat");
1419
1420KCa0       = load ("original/KCa.dat");
1421KCa1       = load ("nemo_nmodl_generated/KCa.dat");
1422
1423KM0      = load ("original/KM.dat");
1424KM1      = load ("nemo_nmodl_generated/KM.dat");
1425
1426KV0      = load ("original/KV.dat");
1427KV1      = load ("nemo_nmodl_generated/KV.dat");
1428
1429SK20      = load ("original/SK2.dat");
1430SK21      = load ("nemo_nmodl_generated/SK2.dat");
1431
1432Na0      = load ("original/Na.dat");
1433Na1      = load ("nemo_nmodl_generated/Na.dat");
1434
1435NaR0      = load ("original/NaR.dat");
1436NaR1      = load ("nemo_nmodl_generated/NaR.dat");
1437
1438NaP0      = load ("original/NaP.dat");
1439NaP1      = load ("nemo_nmodl_generated/NaP.dat");
1440
1441
1442subplot(3,4,1);
1443plot(CaHVA0(:,1),CaHVA0(:,2),CaHVA1(:,1),CaHVA1(:,2),'linewidth',2);
1444title ("CaHVA current");
1445
1446subplot(3,4,2);
1447plot(CaLVA0(:,1),CaLVA0(:,2),CaLVA1(:,1),CaLVA1(:,2),'linewidth',2);
1448title ("CaLVA current");
1449
1450subplot(3,4,3);
1451plot(HCN10(:,1),HCN10(:,2),HCN11(:,1),HCN11(:,2),'linewidth',2);
1452title ("HCN1 current");
1453
1454subplot(3,4,4);
1455plot(HCN20(:,1),HCN20(:,2),HCN21(:,1),HCN21(:,2),'linewidth',2);
1456title ("HCN2 current");
1457
1458subplot(3,4,5);
1459plot(KA0(:,1),KA0(:,2),KA1(:,1),KA1(:,2),'linewidth',2);
1460title ("KA current");
1461
1462subplot(3,4,6);
1463plot(KCa0(:,1),KCa0(:,2),KCa1(:,1),KCa1(:,2),'linewidth',2);
1464title ("KCa current");
1465
1466subplot(3,4,7);
1467plot(KM0(:,1),KM0(:,2),KM1(:,1),KM1(:,2),'linewidth',2);
1468title ("KM current");
1469
1470subplot(3,4,8);
1471plot(KV0(:,1),KV0(:,2),KV1(:,1),KV1(:,2),'linewidth',2);
1472title ("KV current");
1473
1474subplot(3,4,9);
1475plot(SK20(:,1),SK20(:,2),SK21(:,1),SK21(:,2),'linewidth',2);
1476title ("SK2 current");
1477
1478subplot(3,4,10);
1479plot(Na0(:,1),Na0(:,2),Na1(:,1),Na1(:,2),'linewidth',2);
1480title ("Na current");
1481
1482subplot(3,4,11);
1483plot(NaR0(:,1),NaR0(:,2),NaR1(:,1),NaR1(:,2),'linewidth',2);
1484title ("NaR current");
1485
1486subplot(3,4,12);
1487plot(NaP0(:,1),NaP0(:,2),NaP1(:,1),NaP1(:,2),'linewidth',2);
1488title ("NaP current");
1489
1490print  ("NEURON_Original_Vclamp.eps", "-depsc");
1491
1492EOF
1493)
1494
1495))
Note: See TracBrowser for help on using the repository browser.