Computer Programs in Biomedicine 7 (1977) 171-178 © Elsevier/North-Holland Biomedical Press

171

A METHOD FOR EVALUATION OF THE ELECTRIC POTENTIAL FIELD OF A NEURON WITH TIME-DEPENDENT SURFACE POTENTIAL Ivan DVORAK Mathematical Center of the Institutes of Biological Research, Czechoslovak Academy of Sciences, Budd]ovickh 1083, Prague 4, 142 20, Czechoslovakia and Lurrt(r SRCH Hybrid Computation Laboratory, Institute of Hematology and Blood Transfusion, U Nemocnice 5, Prague 2, 128 20, Czechoslovakia Results of the mathematical simulation of the extraceUular electric field of the neuron are presented, yielding the formula describing the time and space dependence of the extracellular potential. The general method of the evaluation of the formulae is presented, which is applicable for any computation of the electric field surrounding a body having time-dependent surface potential. The implementation of the resulting formulae in FORTRAN IV for the special case of neuronal model is presented, including the list of input and output data. The possible applications of the suggested method are studied in concluding discussion. Neurophysiology

Potential field

Hybrid computer

Model of neuron

Laplace equation

1. Introduction

2. The model of the neuron

The contemporary neurophysiology, thanks to the steady improvement in the recording technique, pays serious attention to the time and space dependence of the extracellular potential of an active neuron. The understanding of the recordings of spike waveforms requires a framework dealing with the model of the extracellular field. From the general point of view this model leads to the evaluation of electric potential in the surroundings of a body having time-dependent surface potential. This evalution is a common problem appearing often in different disciplines of contemporary biology. Therefore a common model enabling the computations of a wider class of the models is desirable. The time and space dependence of the electric field are usually complicated thus some simplifications of the model are necessary, which can influence its accuracy. Methods presented in this paper suggest a possible way how the electric field can be described and materialization of this method within the model o f a neuron together with digital computer program shows the way how acceptable agreement between computed and experimental data can be reached.

Rosenfalck [1] and Golovtchinskij [2,3] have proved that the extracellular electric field of a neuron can be satisfactorily approximated by a quasistationary field. The time and space dependence of the extracellular potential ~ representing the spike activity can therefore be obtained by the solution of Laplace equation for time variable boundary conditions: A~=O

(1)

This equation is quite common for various biological models. If the boundary conditions fulfil some reasonable constraints the solution of eq. (1) can be expressed in an integral form. For the purposes of modeling neuronal electric field the formula was derived by Plonsey [4,5] as follows 1 f[~b

-,

3

1

in(7, t) "

1

t

17-

' (2)

where ~(~*, t) represents the deviation of the extracellu-

172

L Dvo~dk, L. Srch, Evaluation o f the potentialft'eld o f a neuron

lar potential (at point ~ of external space in the time instant t) from its resting value, ep(r~ t) the deviation of the surface potential (at point 7 of the neuronal surface membrane in time instant t) from its resting value, tn(r, t) the normal component of the transmembrane current density, o the conductivity of the extracellular space, S surface of the neuronal membrane and O/On the normal differentiation to this surface directed outward in the extracellular space. The formula (2) expresses the dependence of the extracellular potential on time and on the coordinates of the measuring point as well as on the boundary conditions of the formula (1). The main problem imbedded in modeling of the extracellular field of the active neuron is the expression of the boundary conditions precisely enough. The detailed model of the surface properties of a neuron has been suggested by one of the authors [6]. It approximates the neuronal shape of the neuron by the rotationally symmetrical body (dendrites being neglected) and assumes the spreading of the depolarization originated in one point upon axon and soma a s well (see fig. 1), the time courses of the surface potential and the normal component of the transmembrane current density being approximated by a simple function in accordance with the results published by Hodgkin-Huxley [7]. The more detailed description will be published elsewhere. The main goal of the model is description of the surface properties of the neuron by several parameters,

namely by: (i) dimensions of the neuron (ii) maximal values of the surface potential and the normal component of the transmembrane current density (iii) the conductivity of the external medium• To achieve the conveniently warting model of the extracellular field describing the waveform of the potential as the function of parameters mentioned above, the modification of the formula (2) with respect to the used type of computer is desired. Having an opportunity to materialize the model of the neuron on EAI-690 hybrid computer, the authors were not restricted by the computational technique. They had the possibility of choice among analog, hybrid and digital ways of solution. The pure analog solution of the eq. (1) (not even of its modification (2)) with our boundary conditions is not possible, because of the great number of variables and parameters. The speed of computations and the possibility of the continuous mapping of variables resuiting in great accuracy of computations are the advantages of hybrid solution when compared with the pure digital one. However, in our case these advantages are not so expressive as - with regard to the great number of variables and parameters - some discretization is necessary in the hybrid solution as well as in the digital one.

On the other hand the great scale of the application of digital computers in life sciences speaks in favour of

z /

/ SSIBLE HEASURING POINT

.....

\

Fig. 1. The shape of the body approximating the neuron.

........

..........

....

L Dvofdk, L. Srch, Evaluation of the potentialfieM of a neuron application of digital solution, as the digital program can be easily used in the common use. In the preliminary analysis of the problem a way of modification of eq. (2) was found, simplifying considerably the computations. This, together with greater applicability of digital computers, was the main reason why the digital solution was chosen finally.

integration in this formula over surface S can be at every time instant t substituted by integration over ~2(t) representing a part of S where ~p(t') and in(t') differ from zero• If we introduce the notion I

___~---~ - G(r, ~ )

(6)

iT0 1 -+~. On 1 7 - ~[ = F( r, ~ )

3. Description of the computational method

173

(7)

the eq. (2) can be rewritten According to the previous paragraph the formula (2) must be modified into the form allowing effective digital computations. The following method - primarily developed for the model of the extracellular electric field of the neuron - is applicable for any simulation of the electric field around a body having timedependent surface potential if the following conditions are fulfilled: (i) The extracellular potential can be expressed in form (2). (ii) The time and space dependence of surface variables Op(7, t) and in(7, t) can be expressed in the forms (3) and (4)

%(L t) = •

(3)

--..).

in(r, t) = in(t'),

(4)

where t' = t - r ( 7 ) .

1

, ~ in(t') -+ Ckp(t)F(r,-~)+ o G ( r , g ) d a

O(ff, t)=~-~ f

a(t)

(8)

Functions (~p(t') and in(t' ) depend on (t') in the same way as on t at the initial point, where relation t' = t is valid. This dependence can be derived from the Hodgkin-Huxley membrane model [7]. The waveforms of mentioned functions differ from zero only in finite time interval (O, T), Tbeing about 2 - 4 ms; thus we consider in our model only this interval to be impor. tant. For the sake of the discretization of the formula (8) we approximate continuous functions (gp(t') and f in(t' ) in interval (O, T) by step functions composed from N constant segments of equal length. Dividing int t t terval (O, T) into N intervals (tk, tk+ 1), where t k is given by relations (9) and (10) t

tk = kX

k = O, 1,2,3 .... N - 1

(5)

Function r(-~) in the eq. (5) is a growing function of 171 equal to the delay of time courses of functlons • 4~p(r, -+ t) and in( ' ~r, t) at the surface point ~" in comparison with those at initial point r 0. • --)(iii) The surface variables ~p(-~, t) and in( r, t) are identically zero for negative values of t'. In the case of the mentioned model of a neuron these conditions are fulfilled, see [6]. Suggested modifications of eq. (2) consist of two steps: (i) Discretizations of continuous variables or their approximation by constant-in-parts functions. (ii) Approximation of the surface integral by the double sum. As follows from the formula (2) a part of surface S where q~p(t') = O, in(t') = 0 has no influence on the waveform of the extracellular potential. That is why

tN = T, N

(9)

I~

(10)

integer + 1 ,

and introducing the formulas

pk=~lu--tik+l (~p(t')dt'

(11 )

~k+l i.(t')dt',

(12)

k _

1

q -2/6, _/0 3'k(t') - ~ 1

t'~(t'k,

tlc+ 1)

(13)

t' E (t~, t~+ ?)

where U and I are maximal values of 4~p(t') and in(t') respectively in the interval (0, T), formulas (3) and

174

L Dvofdk, L. Srch, Evaluation o f the potential fieM o f a neuron

(4) are simplified into

(21) into (8) and denoting ¢/(~) = ¢(~, ti) results in

N

Ct,(t') = U ~

7k(t')p k

i

(14)

E f

k=O N

in(t')=I

~

477" m = l

7k(t')q k .

(15)

k=O

ff_mF(r

The time course of function q~(ff,t) at chosen point is the desired result of computations. Experimental data prove that it is a continuous function of time, thus without loss of information, $(~, t) can be computed in several discrete time instants only. This simplifies computations immensely. We chose sampling instants equidistantly distributed with interval k between consecutive samples. Denoting the sampling instants as t/, we can write tj =/k

/ = 0,1,2 ... Z .

(16)

We shall denote by symbol ~,n the part of S covered by the spreading depolarization during time interval ((m - 1)A, mk). In contrast to the previous considerations we are speaking about real time now. The area of depolarized membrane in time t i = jA is then f2(ti). According to condition (ii) (paragraph 4), the following relations are true

G ( r-~,

d~ M .

(22)

The indexes k,/, m are not independent. If we indicate as M the largest value of m for which ~m is not an empty region, the relation (23)

Z =N +M

is fulfilled. Taking (23) into consideration the summation in (22) should be performed only over the range (max(1 , j - N ) , min(1, M)) since other terms are equal to zero. Denoting un(

T) = f

F(r,

~)da

rn ,

(24)

$2m

Gm(~) = J G(r, ~ ) d ~ rn ,

(25)

12m

eq. (22) can be rewritten in the form (26) min(j,M)

J f2(t/) = U ~2rn ,

(17)

rn=l

~2m=i n~2m=°=O

, ~ ) + _1q l - .m o

i 4~v,

i,v=l,2...

(18)

Because of these relations the integral over g2(ti) can be divided into the sum of the integrals, each of them being taken only over one special I2m. Let us study the function 3'k(t') at the time instants t i. Taking into consideration only points/(-7) belonging to one of above mentioned ~m, the relation

r(~) E ((m - 1)k, m•) ¢,/(-7) E ~rn

(19)

is valid. Eqs. (19) and (7) result in the relation t'(P, ti) E ( q - re)X, (j - m + l )k) ¢~ / ( 7 ) E ~2rn . (20)

For the function 7k(t') from (13) and (20) follows 7k = i - m (ti - r(P)) = 8 m s f o r / ( 7 ) e ~s.

(21 )

Substituting the formulas (14), (15), (16), (17) and

1 E ¢i(5 %. m--max(1,;-a [Upi_mFmiT)

I . m m ~''1 +oq'G (~)J.

(26)

The formula (26) can be used for numerical computations; the only problem arises with the evaluation of integrals (24) and (25). For given value of ~ functions F(-F, ~) and G(F, ~--')are continuous on S and their first differentiations as well, thus an integral over small segment of S can be substituted by the value of integrated function in one point situated in the centre of this segment multiplied by the area of the segment. The integral over whole S can thus be substituted by the sum of integrals figured out by this approximative method. It can be shown that when S is divided into a sufficient number of segments with comparable sizes arising inaccuracy can be neglected. Applying this method we cover all surface S by the net of knot points. One segment is related to each point according to the rule that only one knot point is situated in every segment, all segments together

L DvoTdk, L. Srch, Evaluation of the potentialfield of a neuron cover the whole surface S and there is no common point belonging to two different segments. Indicating the knot points by the index e and denoting the area of the segment described by the knot point r e by Se the formulas (24) and (25) are transformed into Fm(~) : ~

F(re, ~)Se,

(27)

Gm(~) = ~ a ( r-e' >, - +~ ) S e , e~ Tm

(28)

eG 7~n

where 7*n represents the set of knot points belonging to the region ~m. Substituting (27) and (28) into (26) we get the final fornmla

min(j,M) ~b/.(~) - 4n m=max(l,j--N) e E T

Se

e

[u#-mF(re, 7) +-qj-mG(re,-~) I 1, 0

(29)

enabling the evaluation of the potential ~.(~) in sampling instants tj at any chosen point ~ of the extracellular space.

4. Program NEURON The formula (29) was implemented by the program NEURON, its partial terms being evaluated with the aid of the formulas given in [6]. The full description of the program is beyond the scope of this paper, therefore only short characteristics are given.

5. Hardware and software specifications

The program NEURON coded in FORTRAN IV is designed for EAI-690 hybrid computer. The main program is using only the digital part (EAI-640) of the computer, the analog part (EAI-680) is used only when output drawing of the results is desired. Memory required 12 KB. Peripherals used: line printer, card reader, teletype, paper tape punch, X - Y plotter. Number of cards in source desk: ~350. Card punching code: EAL Paper tape punching code ASCII.

175

6. Description of the program The program NEURON is composed from the main program with eight subroutines. The main program makes all 1/0 operations, drives the strategy of the computations, writes the teletype messages to the operator and performs the identification and punching of the results to the paper tape. If desired it starts the subroutine HYPLOT performing the interpolation of results and drawing of the resulting curves on the X - Y plotter. In the other subroutines the computations of implemented functions are performed. 6.1. Input data The input data of the program NEURON are read on the card reader. They are divided into three independent groups: (i) Dimensions and properties of the neuron (see fig. 1). 1st card: a , b , f , d , u , R 2nd card: r, L,.Ia, K Here/2 means the velocity of the spreading depolarization over the surface of the neuron and K is the constant characteristic for the shape of the neuron (for detailed description se6 [6]). (ii) Properties of external medium 3rd card: U, I, o (see formula (29)) (iii) Cylindrical coordinates of measuring point 4th card: z, p, (see fig. 1). All three groups of data are read only at the beginning of the task. When the computations are repeated, it is possible to read only the group containing the changed data, the other two remaining unchanged. On the teletype the messages are written for the operator according to which he can operate the changes of the input data during repeated computations. 6.2. Output data The main output of the program NEURON is achieved by the line printer. All output data are collected in output record (see sample output record), which contains all input data and the table of the values of the extracellular potential of the neuron at the measuring instants. This record is printed even if a

176

1. Dvo~dk, L. Sreh, Evaluation o f the potential field o f a neuron

\ START

Reading I U,~

_~Reading dimensions

h

/

of a neuron

_ _ ~

the c o o r ~ a t e s of a measuring Point

j

/ Eh~D

Punching

Teletype

results into paper tape

Computation of p a r t i a l

Readi~.

command

results on the X-Y

Reading stored results

of output

of the

record

p o t e n t i a l ~)

plotter

from paper tapd

Flowchart of the program NEURON part o f the data is changed. When desired, by the special order from the teletype the resulting tables of the values of the extracellular potential are punched to the paper tape. This output enables another compilation of the results by the program HYPLOT.

690. It is an independent unit which enables interpolation of the discrete results of the NEURON program into continuous curve and its drawing on the X - Y

~b [ m Y ]

6.3. Program H Y P L O T

0.50

There is a special subroutine named HYPLOT among the subroutines of the NEURON program. This subroutine does not connect directly with the model of a neuron; it is also the only one which employs the analog part (EAI-680) of the computer EAI-

0.2!

• 1

°.

°(

Oi}

/,x~'{

~,U

J o..];~T ~ii>7 / , ( : L!. L&;s SPACE C~¢.~ :,£D /~,,,F(!fP i S I.TS 5.,.i l