-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathOriginOfTheHighFrequencyExtraCellularSignal.html
1712 lines (1448 loc) · 73.3 KB
/
OriginOfTheHighFrequencyExtraCellularSignal.html
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
<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<!-- 2016-02-18 jeu. 14:27 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Origin Of The (High Frequency) Extra-cellular Signal</title>
<meta name="generator" content="Org-mode" />
<meta name="author" content="Christophe Pouzat" />
<style type="text/css">
<!--/*--><![CDATA[/*><!--*/
.title { text-align: center;
margin-bottom: .2em; }
.subtitle { text-align: center;
font-size: medium;
font-weight: bold;
margin-top:0; }
.todo { font-family: monospace; color: red; }
.done { font-family: monospace; color: green; }
.priority { font-family: monospace; color: orange; }
.tag { background-color: #eee; font-family: monospace;
padding: 2px; font-size: 80%; font-weight: normal; }
.timestamp { color: #bebebe; }
.timestamp-kwd { color: #5f9ea0; }
.org-right { margin-left: auto; margin-right: 0px; text-align: right; }
.org-left { margin-left: 0px; margin-right: auto; text-align: left; }
.org-center { margin-left: auto; margin-right: auto; text-align: center; }
.underline { text-decoration: underline; }
#postamble p, #preamble p { font-size: 90%; margin: .2em; }
p.verse { margin-left: 3%; }
pre {
border: 1px solid #ccc;
box-shadow: 3px 3px 3px #eee;
padding: 8pt;
font-family: monospace;
overflow: auto;
margin: 1.2em;
}
pre.src {
position: relative;
overflow: visible;
padding-top: 1.2em;
}
pre.src:before {
display: none;
position: absolute;
background-color: white;
top: -10px;
right: 10px;
padding: 3px;
border: 1px solid black;
}
pre.src:hover:before { display: inline;}
pre.src-sh:before { content: 'sh'; }
pre.src-bash:before { content: 'sh'; }
pre.src-emacs-lisp:before { content: 'Emacs Lisp'; }
pre.src-R:before { content: 'R'; }
pre.src-perl:before { content: 'Perl'; }
pre.src-java:before { content: 'Java'; }
pre.src-sql:before { content: 'SQL'; }
table { border-collapse:collapse; }
caption.t-above { caption-side: top; }
caption.t-bottom { caption-side: bottom; }
td, th { vertical-align:top; }
th.org-right { text-align: center; }
th.org-left { text-align: center; }
th.org-center { text-align: center; }
td.org-right { text-align: right; }
td.org-left { text-align: left; }
td.org-center { text-align: center; }
dt { font-weight: bold; }
.footpara { display: inline; }
.footdef { margin-bottom: 1em; }
.figure { padding: 1em; }
.figure p { text-align: center; }
.inlinetask {
padding: 10px;
border: 2px solid gray;
margin: 10px;
background: #ffffcc;
}
#org-div-home-and-up
{ text-align: right; font-size: 70%; white-space: nowrap; }
textarea { overflow-x: auto; }
.linenr { font-size: smaller }
.code-highlighted { background-color: #ffff00; }
.org-info-js_info-navigation { border-style: none; }
#org-info-js_console-label
{ font-size: 10px; font-weight: bold; white-space: nowrap; }
.org-info-js_search-highlight
{ background-color: #ffff00; color: #000000; font-weight: bold; }
/*]]>*/-->
</style>
<script type="text/javascript">
/*
@licstart The following is the entire license notice for the
JavaScript code in this tag.
Copyright (C) 2012-2013 Free Software Foundation, Inc.
The JavaScript code in this tag is free software: you can
redistribute it and/or modify it under the terms of the GNU
General Public License (GNU GPL) as published by the Free Software
Foundation, either version 3 of the License, or (at your option)
any later version. The code is distributed WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU GPL for more details.
As additional permission under GNU GPL version 3 section 7, you
may distribute non-source (e.g., minimized or compacted) forms of
that code without the copy of the GNU GPL normally required by
section 4, provided you include this license notice and a URL
through which recipients can access the Corresponding Source.
@licend The above is the entire license notice
for the JavaScript code in this tag.
*/
<!--/*--><![CDATA[/*><!--*/
function CodeHighlightOn(elem, id)
{
var target = document.getElementById(id);
if(null != target) {
elem.cacheClassElem = elem.className;
elem.cacheClassTarget = target.className;
target.className = "code-highlighted";
elem.className = "code-highlighted";
}
}
function CodeHighlightOff(elem, id)
{
var target = document.getElementById(id);
if(elem.cacheClassElem)
elem.className = elem.cacheClassElem;
if(elem.cacheClassTarget)
target.className = elem.cacheClassTarget;
}
/*]]>*///-->
</script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
displayAlign: "center",
displayIndent: "0em",
"HTML-CSS": { scale: 100,
linebreaks: { automatic: "false" },
webFont: "TeX"
},
SVG: {scale: 100,
linebreaks: { automatic: "false" },
font: "TeX"},
NativeMML: {scale: 100},
TeX: { equationNumbers: {autoNumber: "AMS"},
MultLineWidth: "85%",
TagSide: "right",
TagIndent: ".8em"
}
});
</script>
<script type="text/javascript"
src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML"></script>
</head>
<body>
<div id="content">
<h1 class="title">Origin Of The (High Frequency) Extra-cellular Signal</h1>
<div id="table-of-contents">
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#orgheadline1">1. Introduction</a></li>
<li><a href="#orgheadline4">2. Relation between membrane potential and extracellular potential</a>
<ul>
<li><a href="#orgheadline2">2.1. Basic equations</a></li>
<li><a href="#orgheadline3">2.2. Membrane current density</a></li>
</ul>
</li>
<li><a href="#orgheadline13">3. Numerical integration of the H & H equation</a>
<ul>
<li><a href="#orgheadline8">3.1. A standardized form for the non-linear reaction-diffusion equations</a>
<ul>
<li><a href="#orgheadline5">3.1.1. The heat equation</a></li>
<li><a href="#orgheadline6">3.1.2. Adding the reaction term: Lee's method</a></li>
<li><a href="#orgheadline7">3.1.3. Boundary conditions</a></li>
</ul>
</li>
<li><a href="#orgheadline9">3.2. Python code doing the job</a></li>
<li><a href="#orgheadline10">3.3. Some functions definitions</a></li>
<li><a href="#orgheadline11">3.4. Checking against Cooley and Dodge (1966) results</a></li>
<li><a href="#orgheadline12">3.5. Checking against the traveling wave solution</a></li>
</ul>
</li>
<li><a href="#orgheadline15">4. Numerical investigation of the radius effect on the extracellular potential</a>
<ul>
<li><a href="#orgheadline14">4.1. A fast way to do the job</a></li>
</ul>
</li>
</ul>
</div>
</div>
<div id="outline-container-orgheadline1" class="outline-2">
<h2 id="orgheadline1"><span class="section-number-2">1</span> Introduction</h2>
<div class="outline-text-2" id="text-1">
<p>
We want to explore here the basic properties of the <i>extra-cellular</i> potential generated by a uniform active—that is, able to propagate an action potential (or <i>spike</i>)—cable or axon. We are going to use the conductance model of <a href="http://onlinelibrary.wiley.com/doi/10.1113/jphysiol.1952.sp004764/abstract">Hodgkin and Huxley (1952)</a> together with the cable model making the "full" H & H model.
</p>
<p>
Remember that H & H <i>did not</i> solve their full model in their <i>opus magnum</i>, remember also that the <i>mechanical calculator</i> they had at this time was far less powerful that any of the smartphones everyone has nowadays in his/her pocket. They used a "trick" looking at the propagation of a waveform without deformation at a constant speed \(\theta\), that is, a spike or action potential. In this way the spatial derivatives of the membrane potential can be expressed as time derivatives (as we will see bellow) and the partial differential equation (PDE) of the full model can be replaced by ordinary differential equations (ODE).
</p>
<p>
So we are going to start by deriving the expression of the extra-cellular potential generated by a "cable like neurite"—a neurite with a large length to radius ratio—that can approximated by a <i>line source</i>. We will follow a classical development that is very clearly explained in the book of Plonsey and Barr (2007) <i>Bioelectricity. A Quantitative Approach</i>, published by Springer. This development will lead us to an equation relating the extra-cellular potential to the integral of the weighted second partial derivative of the membrane potential with respect to space, \(\partial^2 V_m(x,t) / \partial x^2\). The H & H model will give us actual values for this derivative but that will require a numerical solution. We will then explore the effect of the axonal diameter on the extra-cellular potential.
</p>
</div>
</div>
<div id="outline-container-orgheadline4" class="outline-2">
<h2 id="orgheadline4"><span class="section-number-2">2</span> Relation between membrane potential and extracellular potential</h2>
<div class="outline-text-2" id="text-2">
</div><div id="outline-container-orgheadline2" class="outline-3">
<h3 id="orgheadline2"><span class="section-number-3">2.1</span> Basic equations</h3>
<div class="outline-text-3" id="text-2-1">
<p>
The electrostatic potential \(\Phi_e\) [mV] generated by a constant <b>point source</b> of intensity \(I_0\) [mA] is given by:
</p>
\begin{align}\label{eq:stat}\tag{1} \Phi_e = \frac{1}{4 \pi \sigma_e} \frac{I_0}{r} \, ,\end{align}
<p>
where \(\sigma_e\) [S/cm] is the conductivity of the extracellular medium assumed homogeneous and \(r\) [cm] is the distance between the source and the electrode (Plonsey and Barr, 2007, <i>Bioelectricity: A Quantitative Approach</i>, p. 29).
</p>
<p>
For an extended source with a large length to diameter ratio (a cable) that can be approximated by a <b>line source</b>; the generalization of the previous equation for a continuous line source along the x axis (between \(x_{min}\) and \(x_{max}\)) of the 3D Euclidean space equipped with Cartesian coordinates when the electrode is located at \((X,Y,Z)\) is:
</p>
\begin{align}\label{eq:stat1}\tag{2} \Phi_e(X,Y,Z) = \frac{1}{4 \pi \sigma_e} \int_{x_{min}}^{x_{max}} \frac{i_m(x)}{r(x)} dx \, ,\end{align}
<p>
where:
</p>
\begin{align}\tag{3} r(x) \doteq \sqrt{(x-X)^2+Y^2+Z^2}\;,\end{align}
<p>
and \(i_m(x)\) [mA/cm] is the <i>current density</i> at position \(x\) [cm] along the cable.
</p>
</div>
</div>
<div id="outline-container-orgheadline3" class="outline-3">
<h3 id="orgheadline3"><span class="section-number-3">2.2</span> Membrane current density</h3>
<div class="outline-text-3" id="text-2-2">
<p>
We get an expression for \(i_m(x)\) by considering a small piece of cable of radius \(a\) [cm] and of length \(\Delta x\) [cm] (Plonsey and Barr, 2007).
</p>
<p>
If the intracellular potential at position \(x\) is written \(\Phi_i(x)\), then Ohm's law—the current equals the potential drop multiplied by the conductance—implies that the <i>axial current</i> \(I_i(x)\) [mA] is given by (\(\sigma_i\) [S/cm] is the intracellular conductivity):
</p>
\begin{align}
I_i(x) &= -\pi a^2 \sigma_i \frac{\Phi_i(x+\Delta x) -
\Phi_i(x)}{\Delta x} \nonumber \\
&\xrightarrow[\Delta x \to 0]{ } -\pi a^2 \sigma_i \frac{d \Phi_i(x)}{dx} \, . \label{eq:stat2}\tag{4}
\end{align}
<p>
Then the charge conservation implies that the membrane current density \(i_m(x)\) (positive for
an outgoing current) is given by:
</p>
\begin{align}
I_i(x+\Delta x) - I_i(x) &= -i_m(x)\, \Delta{}x \nonumber \\
\frac{d I_i(x)}{dx} &= -i_m(x). \label{eq:stat3}\tag{5}
\end{align}
<p>
Combining equation 4 and equation 5 we get:
</p>
\begin{align}
\label{eq:stat4}\tag{6}
i_m(x) &= \pi a^2 \sigma_i \frac{d^2 \Phi_i(x)}{d x^2}\, .
\end{align}
<p>
Now, writing the membrane potential \(V_m = \Phi_i - \Phi_e\) we have:
</p>
\begin{align}
\label{eq:stat5}\tag{7}
i_m(x) &= \pi a^2 \sigma_i \frac{d^2 V_m(x)}{dx^2} \,.
\end{align}
<p>
This allows us to rewrite equation 2 as:
</p>
\begin{align}
\label{eq:stat6}\tag{8}
\Phi_e(X,Y,Z) = \frac{a^2 \sigma_i}{4 \sigma_e} \int_{x_{min}}^{x_{max}} \frac{1}{\sqrt{(x-X)^2+Y^2+Z^2}}
\frac{d^2 V_m(x)}{dx^2} dx \,.
\end{align}
<p>
The quasi-static approximation (Plonsey, 1967, <i>The bulletin of mathematical biophysics</i> <b>29</b>:657-664; Nicholson and Freeman, 1975, <i>Journal of Neurophysiology</i> <b>38</b>: 356-368)—that
amounts to considering the extracellular medium as purely resistive— leads to
a more general, <b>time dependent</b>, version of equation 8:
</p>
\begin{align}
\label{eq:stat7}\tag{9}
\Phi_e(X,Y,Z,t) = \frac{a^2 \sigma_i}{4 \sigma_e} \int_{x_{min}}^{x_{max}} \frac{1}{\sqrt{(x-X)^2+Y^2+Z^2}}
\frac{\partial^2 V_m(x,t)}{\partial x^2} dx \,.
\end{align}
<p>
<b>Notice that the derivation of equations 8 and 9 does not assume anything about the origin of the membrane potential
non-uniformity.</b>
</p>
<p>
If the membrane potential deviation with respect to rest, \(\Delta{}V_m\), and its derivatives are null at the boundaries of the integration domain, then two rounds of integration by part give (with \(X=0\) and \(h = \sqrt{Y^2+Z^2}\)):
</p>
\begin{align}
\label{eq:statPart}\tag{10}
\Phi_e(h) = \frac{a^2 \sigma_i}{4 \sigma_e} \int_{x_{min}}^{x_{max}} \left(\frac{3 u^2}{(u^2+h^2)^{5/2}} - \frac{1}{(u^2+h^2)^{3/2}}\right) \Delta{}V_m(u) du \, .
\end{align}
<p>
At that stage, in order to go further, we need an explicit expression or value for \(\Delta{}V_m\). We are going to solve numerically the H & H equations for that.
</p>
</div>
</div>
</div>
<div id="outline-container-orgheadline13" class="outline-2">
<h2 id="orgheadline13"><span class="section-number-2">3</span> Numerical integration of the H & H equation</h2>
<div class="outline-text-2" id="text-3">
<p>
We follow here the exposition of Tuckwell (1988) <i>Introduction to theoretical neurobiology. Volume 2.</i> CUP, pp 54-70<sup><a id="fnr.1" class="footref" href="#fn.1">1</a></sup>. We want to solve the following set of equations:
</p>
\begin{align}
C_m \, \frac{\partial V_m}{\partial t} &= \frac{a \sigma_i}{2} \frac{\partial^2 V_m}{\partial x^2} + \overline{g}_K n^4 (V_K-V_m) + \overline{g}_{Na} m^3 h (V_{Na}-V_m) + g_l (V_l - V_m) + I_A \, , \label{eq:HH-PDE}\tag{11}\\
\frac{\partial n}{\partial t} &= \alpha_n(V_m) (1-n) - \beta_n(V_m) n \, , \label{eq:HH-n}\tag{12}\\
\frac{\partial m}{\partial t} &= \alpha_m(V_m) (1-m) - \beta_m(V_m) m \, , \label{eq:HH-m}\tag{13}\\
\frac{\partial h}{\partial t} &= \alpha_h(V_m) (1-h) - \beta_h(V_m) h \, , \label{eq:HH-h}\tag{14}
\end{align}
<p>
where \(C_m\) is the membrane capacitance per unit area [\(\mu{}F/cm^2\)]; \(\overline{g}_K\), \(\overline{g}_{Na}\) and \(g_l\) are the potassium, sodium and leak conductances per unit area [mS/cm\(^2\)]; \(V_K\), \(V_{Na}\) and \(V_l\) are the potassium, sodium and leak currents reversal potentials [mV] and \(I_a\) is "externally" applied current per unit area [\(\mu{}A/cm^2\)] and all the \(\alpha\) and \(\beta\) functions are measured in [1/ms]. <b>You have to be careful here with the units since</b> \(\sigma_i\) <b>is usually given in [S/cm] leading to a first term on the right hand side of equation 11 in</b> [mA/cm\(^2\)] <b>while the left hand side as well as the remaining terms on the right hand side are in</b> [\(\mu{}A/cm^2\)].
</p>
</div>
<div id="outline-container-orgheadline8" class="outline-3">
<h3 id="orgheadline8"><span class="section-number-3">3.1</span> A standardized form for the non-linear reaction-diffusion equations</h3>
<div class="outline-text-3" id="text-3-1">
<p>
We will consider a <i>reaction-diffusion</i> system with the form:
</p>
\begin{align}
\mathbf{u}_t = \mathbf{D} \, \mathbf{u}_{xx} + \mathbf{F}(\mathbf{u}) \, , \label{eq:reaction-diffusion}\tag{15}
\end{align}
<p>
where the \(t\) subscript stands for the partial derivative with respect to time, the \(xx\) subscripts stands for the second partial derivative with respect to position, \(\mathbf{u} = \left(u_1(x,t),\ldots,u_n(x,t)\right)^T \in \mathbb{R}^n\), \(\mathbf{D}\) is a diagonal \(n \times n\) matrix of diffusion coefficients \(\left(D_1,\ldots,D_n\right)\) and \(\mathbf{F}(\cdot) = \left(F_1(\cdot),\ldots,F_n(\cdot)\right)^T\) is a vector-valued function. The corresponds with the above H & H equations is obtained by setting: \(\mathbf{u} = \left(V_m,n,m,h\right)^T\); \(\left(D_1,D_2,D_3,D_4\right) = \left(\frac{a \sigma_i}{2 C_m},0,0,0\right)\), \(F_1(\mathbf{u}) = \left(\overline{g}_K n^4 (V_K-V_m) + \overline{g}_{Na} m^3 h (V_{Na}-V_m) + g_l (V_l - V_m) + I_A\right)/C_m\), \(F_2(\mathbf{u}) \equiv F_2(V_m,n)\), \(F_3(\mathbf{u}) \equiv F_3(V_m,m)\) and \(F_4(\mathbf{u}) \equiv F_4(V_m,h)\) are given by equations 12, 13 and 14.
</p>
</div>
<div id="outline-container-orgheadline5" class="outline-4">
<h4 id="orgheadline5"><span class="section-number-4">3.1.1</span> The heat equation</h4>
<div class="outline-text-4" id="text-3-1-1">
<p>
Let us consider a simpler problem, the <i>heat equation</i>:
</p>
\begin{align}
u_t = D \, u_{xx} \, , \label{eq:heat-equation}\tag{16}
\end{align}
<p>
where \(u(x,t)\) is a scalar. A numerical integration procedure is possible by <i>finite differencing</i>. Here, the heat equation (16) is replaced by a finite difference equation whose solution <i>approximates</i> the one of the heat equation. We discretize the \(x\) axis using \(m+1\) equally spaced points (with a step \(\Delta{}x\)) and the \(t\) axis using \(n+1\) equally spaced times (with a step \(\Delta{}t\)). We write the approximate solution as:
</p>
\begin{align}
U_{i,j} = u(i \Delta{}x,j \Delta{}t) \quad i = 0,\ldots,m \; i = 0,\ldots,n \, . \label{eq:discrete-u}\tag{17}
\end{align}
<p>
The finite difference approximations of the required derivatives are:
</p>
\begin{align}
u_t(x,t) &\approx \frac{U_{i,j+1}-U_{i,j}}{\Delta{}t} \, , \label{eq:u_t}\tag{18} \\
u_x(x,t) &\approx \frac{U_{i+1,j}-U_{i,j}}{\Delta{}x} \, , \label{eq:u_x}\tag{19} \\
u_{xx}(x,t) &\approx \frac{u_x(x,t)-u_x(x-\Delta{}x,t)}{\Delta{}x} \, , \nonumber \\
&\approx \frac{U_{i+1,j}-2 \, U_{i,j} + U_{i-1,j}}{\Delta{}x^2} \, . \label{eq:u_xx}\tag{20} \\
\end{align}
<p>
The numerical integration of the heat equation with the finite difference equation is obtained by establishing a relation between the \(U_{i,j+1}\) and the \(U_{i,j}\). One methods approximates the second spatial derivative at \(t\) by the one at \(t+\Delta{}t\) giving the scheme:
</p>
\begin{align}
\frac{U_{i,j+1}-U_{i,j}}{\Delta{}t} = \frac{D}{\Delta{}x^2} \left(U_{i+1,j+1}-2 \, U_{i,j+1} + U_{i-1,j+1}\right)\, . \label{eq:Ames-scheme}\tag{21}
\end{align}
<p>
<a href="https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method">Crank and Nicolson</a> used the average of the approximations to the second space derivatives at the \(jth\) and \((j+1)th\) time points to get:
</p>
\begin{align}
\frac{U_{i,j+1}-U_{i,j}}{\Delta{}t} = \frac{D}{2 \Delta{}x^2} \left(U_{i+1,j+1}-2 \, U_{i,j+1} + U_{i-1,j+1} + U_{i+1,j}-2 \, U_{i,j} + U_{i-1,j}\right)\, . \label{eq:Crank-Nicolson}\tag{22}
\end{align}
<p>
More generally a weight factor \(\lambda\) can be used with weight \(\lambda\) for the \((j+1)th\) time points and weight \((1-\lambda)\) for the \(jth\) with \(0 \le \lambda \le 1\). Then with:
</p>
\begin{align}
r \doteq \frac{D \Delta{}t}{\Delta{}x^2} \, , \label{eq:step-ratio}\tag{23}
\end{align}
<p>
we have:
</p>
\begin{align}
-r \lambda U_{i-1,j+1} + (1+2 r \lambda) U_{i,j+1} -r \lambda U_{i+1,j+1} = r (1-\lambda) U_{i-1,j} + \left(1-2 r (1-\lambda)\right) U_{i,j} + r (1-\lambda) U_{i+1,j}\, , \label{eq:general-Crank-Nicolson}\tag{23}
\end{align}
<p>
where all the unknown terms in \(j+1\) are on the left side. Since \(i = 0,1,\ldots,m\) there are \(m+1\) equations with \(m+1\) unknown. This integration scheme is called <i>implicit</i> because a linear system must be solved to obtain the values of \(u(x,t)\) at the next time step. The system defined by equation 23 is <i>tridiagonal</i> and can be solved without matrix inversion. In <code>Python</code>, the <a href="http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html">scipy.linalg</a> sub-module provides the <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html#scipy.linalg.solve_banded">solve_banded</a> function to work efficiently with linear systems exhibiting a banded structure.
</p>
</div>
</div>
<div id="outline-container-orgheadline6" class="outline-4">
<h4 id="orgheadline6"><span class="section-number-4">3.1.2</span> Adding the reaction term: Lee's method</h4>
<div class="outline-text-4" id="text-3-1-2">
<p>
We now add a <i>reaction term</i> \(F(u)\) to the scalar heat equation:
</p>
\begin{align}
u_t = D \, u_{xx} + F(u) \, . \label{eq:heat-equation-plus-reaction}\tag{24}
\end{align}
<p>
In the Crank-Nicolson method the second space derivative is approximated by the average of its finite-difference approximations at time points \(j\) and \(j+1\). A similar estimate is needed for \(F(u)\); in other words we need \(F(U_{i,j+^1/_2})\) and we approximate \(U_{i,j+^1/_2}\) by:
</p>
\begin{align}
U_{i,j+^1/_2} &\approx U_{i,j} + (U_{i,j} - U_{i,j-1})/2 \nonumber \\
&\approx \frac{3}{2} U_{i,j} - \frac{1}{2} U_{i,j-1} \, . \label{eq:mid-point}\tag{25}
\end{align}
<p>
And Lees' modification of the Crank-Nicolson method gives the tridiagonal system (remember that \(\lambda\) in equation 23 equals \(^1/_2\) for the Crank-Nicolson method):
</p>
\begin{align}
-\frac{r}{2} U_{i-1,j+1} + (1+r) U_{i,j+1} -\frac{r}{2} U_{i+1,j+1} = \frac{r}{2} U_{i-1,j} + (1-r) U_{i,j} + \frac{r}{2}U_{i+1,j} + \Delta{}t F\left(\frac{3}{2} U_{i,j} - \frac{1}{2} U_{i,j-1}\right)\, . \label{eq:Lees-method}\tag{26}
\end{align}
<p>
Clearly this last equation can only be used if \(j>0\) so for \(j=0\) we use an <i>explicit</i> version (we use \(j\) instead of \(j+1\) in the right hand side of equation 21):
</p>
\begin{align}
U_{i,1} = r \left(U_{i-1,0} -2 U_{i,0} + U_{i+1,0}\right) + \Delta{}t F\left(U_{i,0}\right) + U_{i,0}\, . \label{eq:Lees-method-explicit}\tag{27}
\end{align}
</div>
</div>
<div id="outline-container-orgheadline7" class="outline-4">
<h4 id="orgheadline7"><span class="section-number-4">3.1.3</span> Boundary conditions</h4>
<div class="outline-text-4" id="text-3-1-3">
<p>
There is still one problem to consider before starting writing our code: the <i>boundary conditions</i>, that is what happens at the ends of the cable. There two "extreme" possibilities (and a third one in between the two). The first possibility consists in imposing the voltage at both ends, this leads to the <i>Dirichlet conditions</i>:
</p>
\begin{align}
u(0,t) &= \alpha \, , \label{eq:Dirichlet-0}\tag{28} \\
u(L,t) &= \beta \, . \label{eq:Dirichlet-L}\tag{29}
\end{align}
<p>
The finite difference version is:
</p>
\begin{align}
U_{0,j} &= \alpha \, , \quad j=0,1,\ldots \, , \label{eq:Dirichlet-0-discrete}\tag{30} \\
U_{m,j} &= \beta \, , \quad j=0,1,\ldots \, . \label{eq:Dirichlet-L-discrete}\tag{31}
\end{align}
<p>
These conditions reduce the number of unknown in our linear system by 2, from \(m+1\) to \(m-1\) and correspond to voltage-clamping the ends of the cable.
</p>
<p>
The more common conditions in simulation studies are the <i>Neumann conditions</i> where the values of the space derivatives of the potential are imposed at the ends:
</p>
\begin{align}
u_x(0,t) &= \alpha \, , \label{eq:Neumann-0}\tag{32} \\
u_x(L,t) &= \beta \, . \label{eq:Neumann-L}\tag{33}
\end{align}
<p>
The common values chosen are \(\alpha = \beta = 0\) often referred to as the "sealed ends" conditions—the ones we are going to choose in our numerical implementation. To get the finite difference version, a quick solution would be using \(u_x(x,t) \approx \left(U_{i+1,j}-U_{i,j}\right) / \Delta{}x\), but we can do better—in term of the approximation of the space derivative by its finite difference version at fixed \(\Delta{}x\) using:
</p>
\begin{align}
u_x(i \Delta{}x,j \Delta{}t) &\approx \frac{U_{i+1,j} - U_{i-1,j}}{2 \Delta{}x} \, . \label{eq:central-difference}\tag{34}
\end{align}
<p>
Can you see why? Then the Neumann conditions become:
</p>
\begin{align}
U_{-1,j} &= -2 \alpha \Delta{}x + U_{1,j}\, , \label{eq:Neumann-0-discrete}\tag{35} \\
U_{m+1,j} &= 2 \beta \Delta{}x + U_{m-1,j} \, . \label{eq:Neumann-L-discrete}\tag{36}
\end{align}
<p>
This amounts to introducing "false boundaries" and substituting 35 in 26, the first equation becomes (for \(j>0\)):
</p>
\begin{align}
(1+r) U_{0,j+1} -r U_{1,j+1} = - 2 r \alpha \Delta{}x + (1-r) U_{0,j} + r U_{1,j} + \Delta{}t F\left(\frac{3}{2} U_{0,j} - \frac{1}{2} U_{0,j-1}\right)\, . \label{eq:Lees-left}\tag{37}
\end{align}
<p>
At \(j=0\) the substitution in equation 27 leads to:
</p>
\begin{align}
U_{0,1} = 2 r \left(U_{1,0} - U_{0,0} - \alpha \Delta{}x\right) + \Delta{}t F\left(U_{0,0}\right) + U_{0,0}\, . \label{eq:Lees-left-at-0}\tag{38}
\end{align}
<p>
At the other end we get for \(j>0\):
</p>
\begin{align}
-r U_{m-1,j+1} + (1+r) U_{m,j+1} = 2 r \beta \Delta{}x + r U_{m-1,j} + (1-r) U_{m,j} + \Delta{}t F\left(\frac{3}{2} U_{m,j} - \frac{1}{2} U_{m,j-1}\right)\, , \label{eq:Lees-right}\tag{39}
\end{align}
<p>
while for \(j=0\) we have:
</p>
\begin{align}
U_{m,1} = 2 r \left(U_{m-1,0} - U_{m,0} + \beta \Delta{}x\right) + \Delta{}t F\left(U_{m,0}\right) + U_{m,0}\, . \label{eq:Lees-right-at-0}\tag{40}
\end{align}
</div>
</div>
</div>
<div id="outline-container-orgheadline9" class="outline-3">
<h3 id="orgheadline9"><span class="section-number-3">3.2</span> Python code doing the job</h3>
<div class="outline-text-3" id="text-3-2">
<p>
We are going to solve the standard H & H model using the Neumann boundary conditions with \(\alpha = \beta = 0\) ("sealed ends"). We start by an <code>IPython</code> session—but it wokrs as well with a classical <code>Python</code> session—loading the two main modules we are going on a regular basis, <code>numpy</code> and <code>pylab</code> a sub-module of <code>matplotlib</code>:
</p>
<div class="org-src-container">
<pre class="src src-python">import numpy as np
import matplotlib.pylab as plt
plt.ion()
#%matplotlib inline
plt.style.use('ggplot')
</pre>
</div>
<p>
The three last commands give us <i>interactive</i> graphics (<code>plt.ion</code>) or <i>inline</i> graphics when using the <code>jupyter notebook</code> (in that case, comment the previous line with "#" and uncomment the following one) and a nicer default style for the graphs (<code>plt.style.use('ggplot')</code>). We then assign a few variables considering an axon with a radius \(a\) of 1 \(\mu{}m\) that is \(10^{-4}\) cm (for quantitative data on CNS axons diameters, see <a href="http://www.jneurosci.org/content/32/2/626.abstract">Perge et al (2013)</a>):
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock1">a = 1e-4
</pre>
</div>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock2">Cm = 1.0 # H & H 1952 [μF / cm^2]
rho = 35.4 # H & H 1952, rho is the inverse of σi [Ω cm]
D = a / (2.0 * rho * Cm) # the "Diffusion" constant
D
</pre>
</div>
<pre class="example">
1.4124293785310736e-06
</pre>
<p>
Notice that with this choice of units <code>D</code> is measured in cm\(^2\) / \(\mu{}s\). We define next, for each activation variable, \(n, m, h\) the \(\alpha(v)\) and \(\beta(v)\) functions—where the formal parameter \(v\) stands for the <b>deviation of the membrane voltage with respect to rest</b>—as well as a function returning the steady-state value of the variable at a given voltage. We start with the \(n\) activation variable—the <code>numpy</code> module must have been previously imported with the alias <code>np</code> (<code>import numpy as np</code>)—:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock3">def alpha_n(v):
if np.abs(v-10.0) < 1e-10:
return 0.1
else:
return 0.01*(10.0 - v)/(np.exp((10.0-v)/10.0)-1.0)
def beta_n(v):
return 0.125*np.exp(-0.0125*v)
n_inf = np.vectorize(lambda v: alpha_n(v)/(alpha_n(v) + beta_n(v)))
</pre>
</div>
<p>
Notice that we took care of the special case \(v=10\) using the limit to avoid the undefined expression \(0/0\). The <code>n_inf</code> function has been defined in a <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.vectorize.html#numpy.vectorize">vectorized form</a> since our definition of <code>alpha_n</code> works only with scalar arguments. Having defined these functions it is always a good idea to make a couple of graphs to make sure that we did things properly (we should get figures 4 and 5, p 511 of H & H 1952; <i>don't forget that the membrane voltage convention at that time was the opposite of the one now used</i>):
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock4">vv = np.linspace(-50,110,201)
plt.plot(vv,np.vectorize(alpha_n)(vv),lw=2)
plt.plot(vv,np.vectorize(beta_n)(vv),lw=2)
plt.plot(vv,n_inf(vv),lw=2)
</pre>
</div>
<div class="figure">
<p><img src="figsL1/n_activation.png" alt="n_activation.png" />
</p>
<p><span class="figure-number">Figure 1:</span> \(\alpha_n\) (red), \(\beta_n\) (blue) and \(n_{\infty}\) (violet) as a function of the membrane voltage deviation with respect to rest.</p>
</div>
<p>
We can now define function <code>F_n</code> corresponding to the \(F_2\) of equation 15 whose expression is given by equation 12; this function takes two formal parameters: the membrane potential (deviation) <code>v</code> and the activation variable <code>n</code>:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock5">def F_n(v,n):
if np.abs(v-10.0) < 1e-10:
alpha = 0.1
else:
alpha = 0.01*(10.0 - v)/(np.exp((10.0-v)/10.0)-1.0)
beta = 0.125*np.exp(-0.0125*v)
return alpha*(1-n)-beta*n
vF_n = np.vectorize(F_n)
</pre>
</div>
<p>
It is again a good idea to use these newly defined functions to make sure that nothing "too pathological" happens:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock6">F_n(20,0.6)
</pre>
</div>
<pre class="example">
0.0048690095444177128
</pre>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock7">vF_n([-10,0,10,20,30],[0.1,0.2,0.3,0.4,0.5])
</pre>
</div>
<pre class="example">
array([ 0.01400882, 0.02155814, 0.03690637, 0.05597856, 0.07269618])
</pre>
<p>
Notice that we "redefine" <code>alpha_n</code> and <code>beta_n</code> inside <code>F_n</code>, this is to gain execution time by avoiding function calls. We also define a vectorized version <code>vF_n</code> that will take two formal parameters, <code>v</code> and <code>n</code>, that can be vectors. We proceed in the same way with the \(m\) activation variable:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock8">def alpha_m(v):
if np.abs(v-25.0) < 1e-10:
return 1.0
else:
return 0.1*(25.0 - v)/(np.exp((25.0 - v)/10.0)-1.0)
def beta_m(v):
return 4*np.exp(-.0555*v)
m_inf = np.vectorize(lambda v: alpha_m(v)/(alpha_m(v) + beta_m(v)))
</pre>
</div>
<p>
The graphs (not shown) giving figures 7 and 8 pp 515-516 are obtained with:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock9">vv = np.linspace(-50,110,201)
plt.plot(vv,np.vectorize(alpha_m)(vv),lw=2)
plt.plot(vv,np.vectorize(beta_m)(vv),lw=2)
plt.plot(vv,m_inf(vv)*10,lw=2)
plt.xlim(-10,110)
plt.ylim(0,10)
</pre>
</div>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock10">def F_m(v,m):
if np.abs(v-25.0) < 1e-10:
alpha = 1.0
else:
alpha = 0.1*(25.0 - v)/(np.exp((25.0 - v)/10.0)-1.0)
beta = 4*np.exp(-.0555*v)
return alpha*(1-m)-beta*m
vF_m = np.vectorize(F_m)
</pre>
</div>
<p>
A quick check gives:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock11">vF_m([-10,0,10,20,30],[0.1,0.2,0.3,0.4,0.5])
</pre>
</div>
<pre class="example">
array([-0.59869277, -0.62114902, -0.38730895, -0.06484611, 0.2569922 ])
</pre>
<p>
And for the \(h\) activation variable:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock12">def alpha_h(v):
return 0.07*np.exp(-0.05*v)
def beta_h(v):
return 1.0/(np.exp((30.0 - v)/10.0) + 1.0)
def h_inf(v):
return alpha_h(v)/(alpha_h(v) + beta_h(v))
</pre>
</div>
<p>
Notice that since <code>alpha_h</code> is already (implicitly) vectorized, there is no need to use <code>np.vectorize</code> when defining function <code>h_inf</code>. The graphs (not shown) giving figures 9 and 10 pp 517-518 are obtained with:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock13">vv = np.linspace(-50,110,201)
plt.plot(vv,np.vectorize(alpha_h)(vv),lw=2)
plt.plot(vv,np.vectorize(beta_h)(vv),lw=2)
plt.plot(vv,h_inf(vv),lw=2)
</pre>
</div>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock14">def F_h(v,h):
return 0.07*np.exp(-0.05*v)*(1-h)-1.0/(np.exp((30.0 - v)/10.0) + 1.0)*h
vF_h = np.vectorize(F_h)
</pre>
</div>
<p>
A quick check gives:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock15">vF_h([-10,0,10,20,30],[0.1,0.2,0.3,0.4,0.5])
</pre>
</div>
<pre class="example">
array([ 0.10207082, 0.04651483, -0.00604087, -0.09212563, -0.24219044])
</pre>
<p>
We define next <code>F_V</code> corresponding to the \(F_1\) of equation 15. This function takes 5 formal parameters: <code>v</code>, <code>n</code>, <code>m</code>, <code>h</code> and <code>Ia</code> the injected current. The maximal conductances [mS / cm\(^2\)] and reversal potentials [mV] from H & H (1952) are assigned to local variables in the function. A vectorized version is also defined:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock16">def F_V(v,n,m,h,Ia):
GNa, GK, GL = 120.0, 36.0, 0.3 # H & H 1952
ENa, EK, EL = 115.0, -12.0, 10.5987 # H & H 1952
return (GK*n**4*(EK-v)+GNa*m**3*h*(ENa-v)+GL*(EL-v)+Ia)/Cm
vF_V = np.vectorize(F_V)
</pre>
</div>
<p>
We can now make a first (explicit) step. We are going to consider a thin cable with a 1 \(\mu{}m\) radius and we start by getting its length constant: \(\lambda = \sqrt{a/2 \rho_i \sigma_m}\). We already set \(\rho_i = 35.4\) [\(\Omega{}\) cm], we get the resting value of \(\sigma_m\) [S / cm\(^2\)] by getting the activation variables values at resting level (don't forget that the conductance densities given by H & H are in [mS]):
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock17">sigma_m_rest = (36*n_inf(0)**4+120*m_inf(0)**3*h_inf(0)+0.3)/1000
sigma_m_rest
</pre>
</div>
<pre class="example">
0.00067725364844574128
</pre>
<p>
This gives us a length constant at rest in cm:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock18">lambda_rest = np.sqrt(1e-4/2/rho/sigma_m_rest)
lambda_rest
</pre>
</div>
<pre class="example">
0.045667548060889344
</pre>
<p>
So our length constant is roughly 500 \(\mu{}m\). We will pick a space discretization step of 50 \(\mu{}m\) (5 \(\times 10^{-3}\) cm) equal to a tenth of the length constant and choose a cable length of 20000 \(\mu{}m\) (2 cm), forty times the length constant. We then choose our time discretization step such that the value \(r\) defined by equation 23 is not too large, say 2 (the reason for using an implicit method like the <a href="https://en.wikipedia.org/wiki/Crank%E2%80%93Nicolson_method">Crank-Nicolson</a> method instead of an explicit one in that the latter is stable only if \(r \le 0.5\)). That gives us for \(\Delta{}t\) (remember that our <code>D</code> above is in cm\(^2\) / \(\mu{}s\) and we want a result in \(ms\)):
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock19">Delta_x = 5e-3
r = 2
Delta_t = r*Delta_x**2/D/1000
Delta_t
</pre>
</div>
<pre class="example">
0.0354
</pre>
<p>
To be on the safe side, we will pick a \(\Delta{}t\) of 0.025 ms:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock20">Delta_t = 0.025
</pre>
</div>
<p>
We now need 4 vectors containing the membrane voltage (deviation) and the value of each activation variable at each discrete location along our cable:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock21">L = 2
M = L/Delta_x
v_0 = np.zeros(M+1)
n_0 = np.ones(M+1)*n_inf(0)
m_0 = np.ones(M+1)*m_inf(0)
h_0 = np.ones(M+1)*h_inf(0)
</pre>
</div>
<p>
We also need a vector of the same length with the injected current density at each point along the axon:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock22">Ia_0 = np.zeros(M+1)
Ia_0[0] = 1000.0
</pre>
</div>
<p>
We can now define a function performing a single time step with the explicit method using equations 27, 38 and 40:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock23">def explicit_step(v,n,m,h,Ia):
v_new = np.copy(v)
n_new = np.copy(n)
m_new = np.copy(m)
h_new = np.copy(h)
reaction_term = Delta_t * vF_V(v,n,m,h,Ia)
diffusion_term = np.zeros(len(v))
diffusion_term[1:-1] = (v[0:-2]-2*v[1:-1]+v[2:])*r
diffusion_term[0] = 2*r*(v[1]-v[0])
diffusion_term[-1] = 2*r*(v[-2]-v[-1])
v_new += diffusion_term + reaction_term
n_new += Delta_t*vF_n(v,n)
m_new += Delta_t*vF_m(v,m)
h_new += Delta_t*vF_h(v,h)
return v_new,n_new,m_new,h_new
</pre>
</div>
<p>
We perform one explicit step with:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock24">v_1, n_1, m_1, h_1 = explicit_step(v_0,n_0,m_0,h_0,Ia_0)
</pre>
</div>
<p>
The general time step using Lees' method is an implicit one and requires a banded matrix (containing the voltage factor on the right hand side of equations 28, 37 and 39) to be define that's what do now:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock25">A = np.zeros((3,M+1))
A[0,2:] = -r/2.0 # upper diagonal
A[0,1] = -r # upper diagonal
A[1,:] = 1.0 + r # diagonal
A[2,:-3] = -r/2.0 # lower diagonal
A[2,-2] = -r # lower diagonal
</pre>
</div>
<p>
We now define a function doing one Lees' step. The function needs the present and previous (or old) values of v, n, m and h as well as Ia. The function assumes that the banded matrix <code>A</code> above is already available in the environment and loads function <code>solve_banded</code> from <code>scipy.linalg</code> sub-module:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock26">def lees_step(v_old,n_old,m_old,h_old,Ia_old,
v_present,n_present,m_present,h_present,Ia_present):
from scipy.linalg import solve_banded
v_extra = 1.5*v_present-0.5*v_old # extrapolated mid-point value
n_extra = 1.5*n_present-0.5*n_old # extrapolated mid-point value
m_extra = 1.5*m_present-0.5*m_old # extrapolated mid-point value
h_extra = 1.5*h_present-0.5*h_old # extrapolated mid-point value
Ia_extra = 1.5*Ia_present-0.5*Ia_old # extrapolated mid-point value
n_new = np.copy(n_present)+Delta_t*vF_n(v_extra,n_extra)
m_new = np.copy(m_present)+Delta_t*vF_m(v_extra,m_extra)
h_new = np.copy(h_present)+Delta_t*vF_h(v_extra,h_extra)
reaction_term = Delta_t*vF_V(v_extra,n_extra,m_extra,h_extra,Ia_extra)
diffusion_term = (1-r)*np.copy(v_present)
diffusion_term[1:-1] += (v_present[0:-2] + v_present[2:])*r/2.0
diffusion_term[0] += r*v_present[1]
diffusion_term[-1] += r*v_present[-2]
v_new = solve_banded((1,1),A,reaction_term+diffusion_term)
return v_new, n_new, m_new, h_new
</pre>
</div>
<p>
We make one step with:
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock27">v_2,n_2,m_2,h_2 = lees_step(v_0,n_0,m_0,h_0,Ia_0,v_1,n_1,m_1,h_1,Ia_0)
</pre>
</div>
<p>
Now 2000 more steps stopping the stimulation after 2 ms or 80 steps (this take a few seconds on my slow laptop):
</p>
<div class="org-src-container">
<pre class="src src-python" id="orgsrcblock28">v_M = np.zeros((2002,int(M+1)))
v_M[0] = v_0
v_M[1] = v_1
n_M = np.zeros((2002,int(M+1)))
n_M[0] = n_0
n_M[1] = n_1
m_M = np.zeros((2002,int(M+1)))
m_M[0] = m_0
m_M[1] = m_1
h_M = np.zeros((2002,int(M+1)))
h_M[0] = h_0
h_M[1] = h_1
for i in range(2,2002):
if i < 80:
v_M[i,:],n_M[i,:],m_M[i,:],h_M[i,:] = lees_step(v_M[i-2,:],n_M[i-2,:],m_M[i-2,:],h_M[i-2,:],Ia_0,
v_M[i-1,:],n_M[i-1,:],m_M[i-1,:],h_M[i-1,:],Ia_0)
if i == 80:
v_M[i,:],n_M[i,:],m_M[i,:],h_M[i,:] = lees_step(v_M[i-2,:],n_M[i-2,:],m_M[i-2,:],h_M[i-2,:],Ia_0,
v_M[i-1,:],n_M[i-1,:],m_M[i-1,:],h_M[i-1,:],0)
if i > 80:
v_M[i,:],n_M[i,:],m_M[i,:],h_M[i,:] = lees_step(v_M[i-2,:],n_M[i-2,:],m_M[i-2,:],h_M[i-2,:],0,
v_M[i-1,:],n_M[i-1,:],m_M[i-1,:],h_M[i-1,:],0)
</pre>
</div>
<p>
We can graph the spatial profile of the membrane potential deviation at different times like every 40 time steps or every ms for the first 10 ms: