Unformatted text preview: Theoretical Foundations of Neural Modeling: BME 580.681 Raimond L. Winslow The Departments of Biomedical Engineering and Computer Science The Johns Hopkins University School of Medicine and Whiting School of Engineering 5505090 (O ce Phone) [email protected] (Internet) September 15, 1992 Contents
1 Linear Cable Theory
1.1 The Cable Equation : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1.1.1 De nition of Speci c Membrane Parameters : : : : : : : : : : : : : 1.1.2 Derivation of the Cable Equation : : : : : : : : : : : : : : : : : : : 1.1.3 Physical Meaning of Rm Cm : : : : : : : : : : : : : : : : : : : : : : r 1.1.4 Physical Meaning of Rm 2 : : : : : : : : : : : : : : : : : : : : : : : : Ri 1.1.5 Alternative Forms of the Cable Equation : : : : : : : : : : : : : : : 1.1.6 Input Conductance of an In nite Passive Cable : : : : : : : : : : : 1.1.7 Further Properties of In nite Cables : : : : : : : : : : : : : : : : : 1.2 Finite Length Dendritic Cable Models : : : : : : : : : : : : : : : : : : : : 1.2.1 Input Conductance of Finite Cylinders : : : : : : : : : : : : : : : : 1.2.2 SteadyState Response to Current Injection : : : : : : : : : : : : : 1.2.3 TimeDomain Solutions of Finite Length Cable Problems: SealedEnd Boundary Conditions : : : : : : : : : : : : : : : : : : : : : : : 1.2.4 TimeDomain Solutions of Finite Length Cable Problems: RC Soma Model : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1.3 Idealized Branched Cable Models : : : : : : : : : : : : : : : : : : : : : : : 1.3.1 Equivalent Cylinder Models: Peripheral Spread of Current : : : : : 1.3.2 Equivalent Cylinder Models: Central Spread of Current : : : : : : : 1.3.3 Examples of Peripheral and Central Current Spread : : : : : : : : : 1.3.4 Estimating the Number of Synaptic Inputs to Equivalent Cylinder Models : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 1.3.5 Estimating the Site of Synaptic Inputs to Equivalent Cylinder Models 2.1 Subdivision of Dendritic Trees into Discrete Compartments : : : : : : : 2.2 Membrane Models : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 2.3 Compartmental Model State Equations in Response to Current Inputs 2.3.1 State Stability of the Homogeneous Equation : : : : : : : : : : 2.3.2 SteadyState Response to Maintained Current Inputs : : : : : : 1 3 3 5 7 7 8 8 9 9 10 12 12 23 30 30 33 35 38 38 3 2 Compartmental Modeling : : : : : : : : : : 44 44 46 47 49 51 2.3.3 Transient Response to Current Inputs : : : : : : : : : : : : : : : 2.3.4 Example: Compartmental Model Impulse Response : : : : : : : 2.4 Compartmental Model State Equations in Response to Synaptic Inputs 2.4.1 A Simpli ed Case : : : : : : : : : : : : : : : : : : : : : : : : : : 2.5 Implementation of Compartmental Models : : : : : : : : : : : : : : : : 2.5.1 Linear Models : : : : : : : : : : : : : : : : : : : : : : : : : : : : 2.5.2 Nonlinear Models : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 52 54 54 57 58 58 59 3 Design Principles 4 References 3.1 Synaptic Saturation : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 61 3.1.1 Single Synaptic Inputs : : : : : : : : : : : : : : : : : : : : : : : : : 61 61 63 2 1 Linear Cable Theory
1.1 The Cable Equation
In this section, we derive a partial di erential equation which can be solved for the distribution of voltage along a uniform dendritic cable with passive membrane properties as a function of both space and time. The equation is known as the cable equation. Reviews which should de nitely be consulted by the reader are Rall (1977) and Jack et al. (1975). 1.1.1 De nition of Speci c Membrane Parameters Consider a small dendritic cylinder with length x and radius r, as shown in Fig. 1.1. Assume that the cylinder has passive, uniform membrane properties. We would like to 2r x
Figure 1.1: Uniform Dendritic Segment derive expressions for the axial resistance ri (resistance to current ow along the cylinder axis), membrane resistance rm (resistance to current ow through the cell membrane), and membrane capacitance cm of the cylinder. To do this, it is useful to de ne three speci c membrane parameters: 3 Ri The resistance to current ow along the axis of a uniform dendritic cylinder 1 cm long, with endface surface area of 1 cm2. Ri is known as speci c axial resistance, and is expressed in units of cm. The inverse quantity Gi is axial conductance, and is expressed in units of Siemens per cm. Rm The conductance of a patch of membrane with area 1 cm2, Gm, is known as speci c membrane conductance, and is expressed in units of Siemens per cm2. The inverse quantity Rm, speci c membrane resistance, is expressed in units of ; cm2. Cm The capacitance of a 1 cm2 patch of membrane is known as speci c membrane capacitance, and is expressed in Farads per cm2. This is usually treated as a biological constant with a value of 1 uF/cm2.
Using these de nitions, the axial resistance ri and conductance gi for a dendritic cylinder of length x can be de ned as r2 G gi = x i ri = rx Ri 2 (1.1) (1.2) An explanation of the numerator of Eq. 1.1 is straightforward. Assuming a uniform distribution of charge carriers throughout the cylinder core, then increasing the crosssectional area of the cylinder increases the number of charge carriers. The ux of charge along the cylinder core at any given voltage (the axial current) therefore increases with increasing crosssectional area. The resistance of a medium is determined by the probability of collision between charge carriers and other molecules. Increasing the length of the cylinder (the denominator term in Eq. 1.1) therefore increases the probability of such collisions, leading to a lower conductance (higher resistance). The membrane conductance, gm, is proportional to the number of ion channels in the membrane. Assuming that channels are distributed uniformly within the membrane, then gm is proportional to the cylinder surface area: gm = 2 r xGm (1.3) rm = 2 Rm x (1.4) r Finally, total membrane capacitance cm is proportional to the cylinder surface area: cm = 2 r xCm
4 (1.5) We begin derivation of the cable equation by decomposing an in nite passive dendritic cable into a number of subcylinders of in nitesimal length xi, as shown in Fig. 1.2. Figure 1.3 shows an equivalent electrical circuit. The resistor and capacitor values in this
Vi1 Ii1 xi1 Vi Ii Vi+1 Ii+1 1.1.2 Derivation of the Cable Equation xi xi+1 Figure 1.2: Partitioning of a uniform cable into subcylinders circuit are those described for the cylinder segment in Eqs. 1.11.5. For simplicity, we have assumed zero extracellular resistance (making the extracellular resistance nonzero does not alter the form of the cable equation see Rall, 1977 and Jack et al, 1975 for derivations). We will apply two physical principles to derive the cable equation: 1) Ohm's
Ii1 Ii1
m Vi1 I
c Inside Ii Vi Ii
m c m Ii+1 I Ii+1 Vi+1 I
c Outside Figure 1.3: Equivalent electrical circuit of a partitioned uniform cable law and 2) Kircho 's current law. Referring to Fig. 1.3 and applying Ohm's law at node i, we obtain Taking the limit as xi goes to zero yields an expression for the axial current Ii: 5 Vi = Vi ; Vi;1 = ;Ii Ri r2xi Vi = Vi;1 ; Ii( Ri r2xi ) (1.6) (1.7) V lim 0 xi = @Vi = ;Ii Ri2 (1.8) xi ! @xi r i Simpli cation yields an expression for the axial current Ii: r2 @Vi (1.9) Ii = ; R @x i i This is a fundamentally important relationship which will be used frequently in our study of cable theory. Equation 1.9 is intuitively reasonable the axial current is large when the voltage gradient along the cable is large. We next apply Kircho 's current law at node i, from which we obtain Ii = Iim + I c + Ii+1 Ii = Ii+1 ; Ii = ;ViGm 2 r xi ; Cm2 r xi @Vi @t Taking the limit as xi goes to zero gives I @I lim 0 xi = @xi = ;Vi Gm 2 r ; Cm2 r @Vi xi ! @t i i
(1.10) (1.11) (1.12) @I An alternative expression for @xii in terms of Vi can be obtained by di erentiation of Eq. 1.9 with respect to xi. Substitution of the result into Eq. 1.12 yields: @ 2Vi = ; @Ii Ri @x2 @xi r2 i = Vi(Gm 2 r) Ri2 + Cm2 r @Vi Ri2 r @t r Ri 2 + C R 2 @Vi = Vi R r m i r @t m
r Multiplication of both sides of the above equation by 2 Rm gives the cable equation: Ri (1.13) (1.14) (1.15) 2 r ( 2 Rm ) @ V2i = Vi + RmCm @Vi R @x @t
i i (1.16) 6 1.1.3 Physical Meaning of Rm Cm Suppose the dendritic cylinder is spaceclamped at a uniform voltage V0 8x 2 (;1 1) at time t = t0, so that @V = 0. Under this assumption, the distribution of voltage along @x the cable remains uniform for all time following release of the voltageclamp, and the cable equation reduces to V + RmCm @V = 0 @t V (t0) = Vo
This is a linear rst order di erential equation in V , with solution
1 (1.17) (1.18) V (t) = V0 e; RmCm (t;t0) 8t t0 (1.19) V (t0) = V0 (1.20) Thus, when a spatially uniform initial voltage V0 is imposed on the cylinder at time t = t0, potential decays exponentially with time constant m at every point along the cylinder following release of the voltageclamp, where the membrane time constant is de ned as
m = RmCm (1.21) In this section, we derive the steadystate response of a passive in nite cable to current injection at x = 0. In the steadystate, @V = 0. Substitution into Eq. 1.16 yields @t
2 ( Rm r ) @ V2 = V Ri 2 @x This is a linear constant coe cient second order ODE, the solution of which is 1.1.4 Physical Meaning of Rm r2 Ri (1.22) (1.23) V = c1e 1 x + c2e; 1 x 8x r = Rm 2 (1.24) Ri Intuitively, it is clear that the solution is symmetric about x = 0, so we will only consider solutions for x 0. In order that the solution remain bounded at x = 1, c1 = 0. This yields 7 s V = c2e; 1 x (1.25) V (0) = c2 (1.26) Equation 1.25 states that in the steadystate, voltage decays exponentially as a function r of distance from the site of current injection, with space constant = Rm 2 . This result Ri holds only for an in nite cable. As we will see later, boundary conditions play a critical role in determining decay of steadystate membrane potential in nite length cylinders. q The membrane time and space constants m and provide natural units in which to express time and distance when solving passive cable problems. Using the de nitions of the time and space constants given by Eqs. 1.24 and 1.21, Eq. 1.16 may be rewritten as 2 2 @ V (x t) = V (x t) + m @V (x t) (1.27) @x2 @t De ne the variables X and T as (1.28) X = x T = t (1.29) Substitution into Eq. 1.27 yields @ 2V (X T ) = V (X T ) + @V (X T ) (1.30) @X 2 @T This is the form of the cable equation in terms of the dimensionless variables X and T . We have seen that the quantities and m are important normalizing parameters when working with properties of the cable equation. There is one more normalizing parameter that will prove to be useful. This parameter is the input conductance of a semiin nite cable , known as G1 (a semiin nite cable is one which begins at location x = 0 and extends to in nity). We have previously derived an expression for the steadystate distribution of voltage V (x) in response to current injection at x = 0 in an in nite cable (Eq. 1.25). We have also derived an expression for the axial current I (x) (Eq. 1.9). Substitution of Eq. 1.25 into Eq. 1.9 yields 8
m 1.1.5 Alternative Forms of the Cable Equation 1.1.6 Input Conductance of an In nite Passive Cable 2 2 I (x) = cR r e; 1 x i The input conductance of a semiin nite cable is therefore given by 2 I 2 2 G1 = V (0) = cR r c1 = Rr (0) i 2 i (1.31) (1.32) We previously derived an expression for axial current at location x (Eq. 1.9). If we let X = x= , then we can rewrite this expression as 2 (X (X (1.33) I (X T ) = ; Rr @V @X T ) = ;G1 @V @X T ) i An alternative expression for G1 can be derived by incorporating the de nition of the space constant (Eq. 1.24) into the above, and de ning the dendritic diameter d = 2r: d3=2 (1.34) G1 = 2(R R )1=2 i m We will not discuss any additional properties of in nite cables. There are three reasons for this. First, the solution of in nite cable problems is quite complicated in the timedomain. Thus, Laplace transform techniques are typically used to solve these problems in the frequency domain, an approach which provides little physical intuition into the nature of cable problems. Second, timedomain solutions for nite cable problems are far more tractable, and an understanding of the role played by boundary and initial conditions in shaping the nature of solutions provides a great deal of insight into the nature of cable problems. Finally, with the development of more and more accurate techniques for recording from single cells, estimates of the electrotonic length of cell dendrites is dropping. Twenty years ago, dendrites were believed to be electrotonically extensive (say 3 4 space constants, approximately equivalent to an in nite cylinder). More recent estimates in a variety of neuron types provide electrotonic length estimates typically less than one space constant. For those interested, extensive descriptions of in nite cable problems are given in the book by Jack et al (1975). 1.1.7 Further Properties of In nite Cables 1.2 Finite Length Dendritic Cable Models
Figure 1.4 shows the simplest possible model of a neuron a membrane cylinder. The cylinder has uniform speci c membrane parameters Cm, Rm, and Ri and nite length l. 9 I0 V0 x=0 Il Vl x=l Figure 1.4: Finite length uniform cable model The variables I0 and V0 are axial current and voltage at location x=0. Variables Il and Vl are corresponding values at location x=l. Despite its simplicity, there are a number of ways in which this structure is relevant to neural modeling: A. The model accurately describes properties of short dendritic segments with passive membrane properties. Understanding the input output properties of such cables will aid in understanding the responses of passive cell models, which are typically represented as the concatenation of dendritic segments B. Under certain conditions, the response properties of branched dendritic trees are equivalent to those of this single cylinder model C. The e ects of certain changes to the model, such as adding a soma at location x=0, can be understood by modi cation of the boundary conditions at x=0. 1.2.1 Input Conductance of Finite Cylinders The following derivation parallels that in Rall (1959). We have previously derived the steadystate solution of the cable equation: 8x2 0 l] (1.35) V (x) = c1e x + c2e; x This equation describes the steadystate distribution of voltage along either nite or in nitelength cables. Di erences in the solutions of these two problems are determined by di erences in the boundary conditions, as expressed by I0 and Il. It is convenient to rewrite Eq. 1.35 in the following form (note that the variables c1 and c2 uniquely specify the variables A and B used below): V (x) = A + B e l e; x + A ; B e; l e x (1.36) 2 2 = A (e l;x + e; l;x ) + B (e l;x ; e; l;x ) (1.37) 2 2 = Acosh( l ; x ) + Bsinh( l ; x ) (1.38) 10 Writing the equations in this form makes it easier to apply the boundary condition at x = l (exponents are all zero at x = l). We previously derived an expression for the axial current at location x in terms of V (x) (Eq. 1.9). Di erentiation of Eq. 1.38 with respect to x, and substitution into Eq. 1.9 yields an expression for I0: 2 r2 (x I0 = ; R @[email protected] t) jx=0 = Rr Asinh( l ) + Bcosh( l )] (1.39) i i Letting L = l , an expression for the input conductance Gin can be obtained from Eqs. 1.39 and 1.38: 2 I (L) + Bcosh Gin = V0 = Rr Asinh(L) + Bsinh(L) (1.40) (L) 0 i Acosh Dividing top and bottom by the term Acosh(L) gives B + tanh(L) (1.41) Gin = G1 ( 1A + B tanh(L) ) A The conclusion is that the input conductance of a nite cable is speci ed given the speci c membrane parameters Ri and Rm , the cylinder length and radius l and r (these quantities specify G1 ), and an as yet unspeci ed parameter B . What is B ? A A Previously we noted that Eq. 1.35 describes the steadystate distribution of voltage along nite length cables as long as the appropriate boundary conditions are applied. We have already applied the boundary condition at x = 0. We will learn about the parameter B by applying the boundary condition at x = l to obtain an expression for the output A conductance. Evaluation of Eqs. 1.9 and 1.38 at x = l yields:
2 r2 (x Il = ; R @[email protected] t) jx=l = Rr B i i Vl = A De ning a new constant Gout , we have 2 I Gout = Vl = Rr B l i A = G1 B A (1.42) (1.43) (1.44) (1.45) (1.46) 11 This yields B = Gout A G1 is therefore the terminating conductance Gout normalized by the input conductance of an in nite version of the same cable, G1. Input conductance of the nite length cable is therefore determined by intrinsic cable properties and the end boundary condition: Gout + tanh(L) Gin = G1 ( G1 Gout ) (1.47) 1 + G1 tanh(L) Using Eq. 1.38, we can express the steadystate voltage distribution along the cable in units of electrotonic distance X as B A 1.2.2 SteadyState Response to Current Injection
V (X ) = Acosh(L ; X ) + Bsinh(L ; X ) V (0) Acosh(L) + Bsinh(L) cosh(L ; X ) + Gout sinh(L ; X ) G1 = Gout sinh(L) cosh(L) + G1 (1.48) (1.49) Thus, the distribution of voltage along a nite length cable is also determined by the terminating boundary condition as speci ed by the value of Gout . 1.2.3 TimeDomain Solutions of Finite Length Cable Problems: SealedEnd Boundary Conditions
Our next task is to solve the cable equation for the distribution of voltage V (X T ) along a passive, uniform cable 8 X 2 0 L], and 8 T T0 in response to an arbitrary current input applied at location X = Xi at time T0. More formally, we seek to solve @ 2V (X T ) = V (X T ) + @V (X T ) (1.50) @X 2 @T subject to boundary conditions (X I0 = ;G1 @V dX T ) jX =0 (1.51) (X IL = ;G1 @V dX T ) jX =L (1.52) and initial condition V (X0 T0) = E (X ) 8X 2 0 L] The solution to this problem involves three steps: 12 (1.53) A. Solution of the cable equation subject to arbitrary boundary conditions at X = 0 and L, and an arbitrary initial voltage distribution (the "initial condition") at time T = T0 B. Specialization of the initial condition so that V (X T ) is the impulse response K (X T Xi ) of the cable to a Diracdelta current pulse applied at location X = Xi C. Use of the convolution integral to obtain the response of the cable to an arbitrary current input given the impulse response function K (X T Xi). Assume a solution of Eq. 1.50 exists with form V (X T ) = F (X )G(T ) (1.54) This solution technique is known as the method of separation of variables, since it assumes that the solution V (X T ) can be decomposed into a component varying only in X and a component varying only in T. As we will see, this assumption reduces the problem of solving a partial di erential equation to one of solving two coupled ordinary di erential equations. Whether or not the method of separation of variables can be applied to any particular partial di erential equation depends on the form of the equation, and on the form and shape of the boundary/initial conditions (Weinberger, 1965). These conditions are met for the cable problems we will study. Substitution of Eq. 1.54 into the cable equation (Eq. 1.50) yields F 00(X )G(T ) = F (X )G(T ) + F (X )G0(T ) (1.55) Rearranging terms gives F 00(X ) = G(T ) + G0(T ) 8X 2 0 L] 8T T (1.56) 0 F (X ) G(T ) What does the above equation tell us? A. Equation 1.56 holds for all T. Set T = T1. Then the righthand side of Eq. 1.56 is a constant. The lefthand side must equal this constant regardless of the value of X. B. Equation 1.56 also holds for all X. Choose X = X1. Then the lefthand side of Eq. 1.56 is a constant. The right hand side must be equal to this constant for all values of T. In particular, equality must hold for T = T1. Items A and B therefore imply that the left and righthand sides of Eq. 1.56 equal the same constant. This yields 13 General Solution of the Finite Length Cable Equation F 00(X ) = G(T ) + G0(T ) = ;u2 (1.57) F (X ) G(T ) As will be seen momentarily, setting the unknown constant equal to ;u2 is purely a matter of convenience. We will assume for generality that u is a complex number given by u = +j . Equation 1.57 yields two equations which may be solved separately for F(X) and G(T): F 00(X ) + u2F (X ) = 0 (1.58) 0 (T ) + (1 + u2)G(T ) = 0 G (1.59) Equation 1.58 is a second order linear constant coe cient ODE with solution F (X ) = AejuX + Be;juX 8X 2 0 L] (1.60) Equation 1.59 is a rst order linear ODE with solution G(T ) = Ce;(1+u2)T 8T T0 (1.61) Therefore, the general form of the solution is V (X T ) = c1ejuX + c2e;juX ]e;(1+u2)T 8X 2 0 L] T T0 (1.62) The three unknowns c1, c2, and u must be determined from the three boundary and initial conditions given by Eqs. 1.51 1.53. In order to proceed, we will assume that the sealedend boundary conditions apply at X = 0 and X = L. We will then derive the particular form of the solution to Eq. 1.62 subject to these boundary conditions (also see (Rall, 1969)). These boundary conditions state that (X (1.63) I0 = ;G1 @V dX T ) jX =0 = 0 (X IL = ;G1 @V dX T ) jX =L = 0 (1.64) Solution of the Finite Length Cable Equation: SealedEnd Boundary Conditions In words, these boundary conditions state that no charge can leak across the membrane at the end of the cable, hence the name sealedend boundary conditions. The boundary condition at X = 0 (Eq. 1.63) implies that @V (X T ) j (1.65) @X X =0 = 0 14 Application of this boundary condition to Eq. 1.62 yields @V (X T ) j juX ;juX (1.66) @X X =0 = ju c1e ; c2e ]jX =0 = ju c1 ; c2] = 0 (1.67) One possible solution to Eq. 1.67 is u = 0. However, substitution into Eq. 1.62 shows that this eliminates the possibility of spatiotemporal variation of the solution, and violates our assumption of a general initial condition E (X ). The only remaining solution of Eq. 1.67 is that c1 = c2. The resulting solution therefore simpli es to V (X T ) = c2 ejuX + e;juX ]e;(1+u2)T 8X 2 0 L] T T0 (1.68) Next, we apply the boundary condition at X = L in conjunction with Eq. 1.68 to obtain @V (X T ) j juL ;juL (1.69) @X X =L = juc2 e ; e ] = 0 Once again, u = 0 is a possible solution, but eliminates the possibility of spatiotemporal variation of the solution and leads to a possible violation of the initial condition. Setting c2 = 0 yields the trivial solution V = 0. This leaves ejuL ; e;juL = 0 (1.70) Substitution of u = + j into the above, and simpli cation yields cos( L) e; L ; e L] + jsin( L) e; L + e L] = 0 (1.71) In order for this equality to hold, both real and imaginary terms in Eq. 1.71 must be zero. Note that the term e; L + e L is always greater than or equal to two. Thus, the imaginary part of Eq. 1.71 is only zero when sin( L) = 0. Under this constraint, cos( L) is always 1. Therefore, equality in Eq. 1.71 holds if and only if sin( L) = 0 ) = n (1.72) L e; L ; e L = 0 ) = 0 (1.73) The values for given above are known as separation constants. Since there are an in nite number, we will denote them as n. Since cos( nX ) = cos(; n X ), we have the following solution for V (X T ): = n L Note that in Eq. 1.74, at any time T = T1, the term
n V (X T ) = n=0 1 X Ancos( nX )e;(1+ 2 )T n 8X 2 0 L] 8T T0 (1.74) (1.75) 15 F (X ) = is the Fourier Series representation of the voltage distribution along the cable, with the n = 0 term giving the DC component of the voltage, and terms for N 1 giving terms of successively higher spatial frequency. Inspection of Eq. 1.74 also shows that the solution of the cable equation subject to the sealedend boundary conditions in response to any current input is a weighted sum of exponentials. It is important to note that the time constants associated with each of these exponentials are speci ed by application of the boundary conditions. These time constants are de ned as follows: m (1.77) n = 1 + 2 n = 0 1 2 3 ::: For the particular case of the sealedend boundary condition,
n n=0 1 X Ancos( nX ) (1.76) n = 0 1 2 ::: (1.78) 1 + ( nL )2 Note that m = 0, the membrane time constant. Time constants n for n 1 are known as equalizing time constants. Equalizing time constants are less than the membrane time constant. Thus, terms expressing the contribution of high spatial frequencies to the solution decay faster than do terms with lower spatial frequency. The weighting factors An for each term in Eq. 1.74 are determined by application of the initial condition in Eq. 1.53. To do this, we will use the fact that the set of functions fcos( n L) n = 0 1 2 :::g are orthogonal over the interval 0 L]. In this sense, the following derivation parallels that for the Fourier series coe cients. For simplicity, set T0 = 0. Application of the initial condition in Eq. 1.53 yields
n = m V (X 0) = Multiplication of both sides of Eq. 1.79 by cos( k X ), and integration from 0 to L yields n=0 1 X Ancos( nX ) (1.79)
1 X An Z L cos( nX )cos( k X )dX Z L V (X 0)cos(
0 k X )dX = Consider the term n=0 0 (1.80) (1.81) Z L cos(
0 n X )cos( k X )dX 16 When n 6= k, the integral evaluates to 2 2 ;1 L n sin( n X )cos( k X ) ; k cos( n X )sin( k X )]( n ; k ) j0 = 2 2 ;1 n sin( n L)cos( k L) ; k cos( n L)sin( k L)]( n ; k ) We have shown previously that the boundary conditions I0 = I1 = 0 imply that sin( 0. Therefore, the integral in Eq. 1.83 is zero when n 6= k. When n = k, L cos2( nX )dX = X + sin(2 n X ) jL 2 4 n 0 0 n = L + sin(2 n L) = L 1 + sin(2 LL) ] 2 4 2 2 Z (1.82) (1.83) n L) = (1.84) (1.85) Incorporation of this result into Eq. 1.80, and application of the initial condition V (X 0) = E (X ) yields an expression for the coe cients An: ;1 L 2 n E (X )cos( n X )dX (1.86) An = L 1 + sin(2 LL) ] 2 n 0 Since sin(2 n L) = sin(2n ) = 0, the expression for An simpli es to 1 L A0 = L E (X )dX (1.87) 0 2 L (1.88) An = L E (X )cos( n X )dX 8n 1 L n n Z Z Z We have therefore obtained the complete solution to Eqs. 1.50 1.53: 1 (1.89) V (X T ) = Ancos( n X )e;(1+( nL )2)T L n=0 1 L A0 = L E (X )dX (1.90) 0 2 L (1.91) An = L 0 E (X )cos( n X )dX 8n 1 L A simple intuitive example highlights the role of the initial condition in shaping the distribution of voltage along the cable for T T0. Assume the initial condition E (X ) = E (uniform initial voltage along the cable). Then Eq. 1.91 reduces to An = 0 when n 6= 0, and An = E otherwise. Therefore, V (X T ) = Ee;(T ;T0) 8X 2 0 L] 8T T0 (1.92) The conclusion is that under the assumption of sealedend boundary conditions and uniform initial conditions, voltage decays exponentially at the same rate everywhere along the 0 X Z Z 17 cable. Put another way, all equalizing time constants are zero except for 0 which is equal to the membrane time constant. The equalizing time constants are nonzero (for the case of sealedend boundary conditions) only when the initial voltage distribution is nonuniform. In this case, the initial condition has a Fourier series representation with terms corresponding to n 1 having sucessively higher spatial frequencies. Subsequent to application of the initial condition, charge is dissipated not only through the membrane, but also by lateral currents owing along the core of the cable. These axial currents are of course predicted when the gradient term in Eq. 1.9 is nonzero (as it is only when V(X,T) is nonuniform). As will be shown in a subsequent section, voltage responses have more rapid transient components when the initial condition is nonuniform simply because there are more pathways by which charge may dissipate at any point along the cable. Impulse Response Function Under SealedEnd Boundary Conditions At this point, we have accomplished our rst goal of solving the cable equation subject to arbitrary boundary conditions (we in fact specialized the boundary conditions, but the procedure is similar for other boundary conditions). Our next goal is to compute voltage along a linear passive cable in response to an arbitrary current input applied at any point. To do this, we must derive the impulse response of the cable (Rinzel and Rall, 1974). The approach to this problem will be illustrated by the following simple example. The
V(t) I(t) R C Figure 1.5: Parallel RC Circuit Model input to the parallel RC circuit shown in Fig. 1.5 is a current I(t), and output is the voltage V(t). V(t) can be computed for any I(t) if we know the impulse response of the system. Suppose a Diracdelta current function is applied at time t=0. The e ect of this input is to place a unit charge on the capacitor so that V (0) = C ;1. If we apply Kircho 's current law to the circuit of Fig. 1.5, we obtain a rst order linear constant coe cient ODE which can be solved subject to the initial condition V (0) = c;1:
1 (1.93) V (t) = C ;1e; RC t If V(t) is in fact the impulse response K(t) of the parallel RC circuit, then we should be able to use K(t) de ned by Eq. 1.93 to compute voltage response to an arbitrary current 18 input. Assume that the input I(t) is a step current input u(t) applied at time t=0. Then
1 V (t) = u( )C ;1e; RC (t; )d (1.94) 0 1 = R(1 ; e; RC t) (1.95) In other words, the unitstep response of this system is a voltage V(t) which rises with time constant RC to a steadystate value of R intuitively the correct result. The purpose of this example is to show that the voltage response of the simple parallel RC circuit to a Diracdelta current input can be determined by placing an instantaneous charge on the membrane capacitance at time t=0, thereby setting the initial voltage across the capacitor. The resulting voltage response of the system for time t 0 seconds is equal to the impulse response. Why is this important? Eqs. 1.891.91 can be used to solve for V(X,T) given the sealedend boundary conditions and an initial voltage distribution E(X). We will construct an E(X) which corresponds to the application of a Diracdelta current function at X = X1. The resulting voltage response will be the impulse response of the cable, K (X T X1). By analogy to our previous example, let the applied current be Q0 (T ) (X ; X ) (1.96) I (X T ) = Zt m {z} Charge per Unit Time 1 Setting Q0 = m corresponds to application of a unit amplitude Diracdelta current input at location X = X1 at time T = 0. The initial voltage distribution corresponding to this input can be determined as follows. Since charge is only placed at X = X1 , we have Q(X 0) = Q0 (X ; X1) (1.97) V (X 0) = V0 (X ; X1) (1.98) The total charge on a cable of length L is given by Z L Q(X 0)dX = Z L CV (X 0)dX
0 0 (1.99) where C is the capacitance of a length of dendritic cylinder. Using the fact that C = Cm d and substitution of the de nitions in Eqs. 1.971.98 into Eq. 1.99 yields (1.100) Q0 = Cm d V0 ) V0 = Rdm = R1 The initial voltage distribution corresponding to a unit delta dirac current input applied at location X = X1 is therefore V (X 0) = R1 (X ; X1 ) (1.101) 19 The voltage response V(X,T) to the initial voltage distribution given by Eq. 1.101 is the impulse response of the cable to a Diracdelta current input applied at location X = X1 denoted by K (X T X1). Substitution of this initial condition into Eqs. 1.891.91 yields the form of the impulse response: R1 e;T + 1 2R1 cos( n X )cos( n X )e; 1+( nL )2]T K (X T X1) = L (1.102) L 1 L n=1 L X Step Response Under SealedEnd Boundary Conditions The result derived in the previous section has theoretical signi cance since it can be used to predict the voltage response of a cable model to an arbitrary current input. Impulse responses can't be measured physiologically since we can't apply Diracdelta current pulses to real neurons. However, we can provide reasonable approximations of step current inputs to real cells. Analyses of these responses will enable us to deduce properties of the neuron being studied. In this section, we will derive the unitstep response of nite length cables subject to sealedend boundary conditions. Since there is a linear relationship between injected current and voltage along the cable, we can use the convolution integral to relate the current injected at X = X1, denoted by I (X1 T ), to voltage V (X T ): V (X T ) = Z T I (X
0 1 S )K (X T ; S X1)dS (1.103) Letting I (X1 T ) = u(T ) and substitution of the form of K (X T X1) into Eq. 1.103 yields 1 2 V (X T ) = R1 (1;e;T )+ 2R1 cos( n X1)cos( n X ) L2 +L(n )2 ] 1;e; 1+( nL )2]T ](1.104) L L L n=1 L Consider the voltage response at point X = 0 in response to current injection at X = 0: 1 2 V (0 T ) = f R1 + 2R1 L2 +L(n )2 ]g (1.105) L n=1 L 1 2 (1.106) ; f R1 e;T + 2R1 L2 +L(n )2 ]e; 1+( nL )2]T g L n=1 L The rst term in the above equation is the steadystate response, whereas the second is the transient response. The input resistance of a nite length cable with sealedend boundary conditions is given by Eq. 1.41: R1 Rin = tanh(L) (1.107) X X X 20 Since we are applying a unit current input, the above expression also equals the steadystate voltage at location X = 0. The step response normalized to a peak value of 1 is therefore L2 ]e; 1+( nL )2]T V (0 T ) = 1 ; tanh(L) e;T ; 1 2tanh(L) (1.108) V (0 1) L L L2 + (n )2 n=1 X Inspection of Eqs. 1.104 and 1.108 shows that the transient component of the unit step response is a weighted sum of exponentials. The rst time constant, m, is the membrane time constant and is equal to RmCm. The remaining time constants are all equalizing time constants, which are typically much faster than m. The ratio of the membrane to equalizing time constants is 2 m = 1 + (n ) (1.109) L n As an example, when L = :5, the ratio of the rst two time constants m is 40. This means 1 that if the membrane time constant is 10 msec, the rst equalizing time constant is 250 usec. Under these assumptions, the step response of the cable is well described by a single exponential for t :250 usec. The di erence in membrane and equalizing time constants can be highlighted by considering only the transient portion of Eq. 1.108: 1 2tanh(L) 2 ;1 V (0 T 1 + ( n ) ] e; 1+( nL )2 ]T (1.110) 1 ; V (0 1)) = 1 tanh(L) e;T + L L L n=1 Plots of Eq. 1.110, frequently referred to as charging curves, are shown in Fig. 1.6. The solid and dotted lines plot the natural logarithm of Eq. 1.110 for L = 2 and :5, respectively. The dashed line plots the natural logarithm of a pure exponential charging curve of the form e;T . A logarithmic plot of a pure exponential yields a straight line with slope of 1 and yintercept 0. The natural log of the charging curve for L = :5 (dotted line) appears to be well described by a straight line of slope 1 for T :05, but is not for smaller values of T. The very brief initial transient seen for values of T < :05 represents a contribution to the shape of the charging curve by the second exponential term in Eq. 1.110, a term 0 with time constant 1. When L = 2, the ratio 1 is 3.45. Continuing with our example, if 0 = m is 10 msec, then 1 = 2:9 msec and 2 = 920 usec. The data shown by the solid line in Fig. 1.6 show clear evidence of multiexponential decay of the charging curve. For large values of T, the curve asymptotically approaches a straight line of slope 1. The rapid initial transient is determined by the more rapidly decaying exponential terms governed by the equalizing time constants 1 and 2. These data suggest that for this particular cable model, both membrane and equalizing time constants can be estimated by tting exponential terms to appropriate regions of the Estimating Membrane and Equalizing Time Constants X 21 ln[1V/Vinf] 1.5
0.0 T 1.0 0.5 0.0 0.2 0.4 0.6 0.8 1.0 Figure 1.6: Natural Log of the Normalized Cable Step Response. charging curve. This procedure for estimating time constants, known as the Rall Peeling Procedure (see Rall, 1960), is illustrated in Fig. 1.7. Panel A shows a least squares t to the late exponential tail of the charging curve for L = 2. The estimated time constant is ^0 = 10:3 msec (assuming an actual value of 10 msec). Panel B shows the residual data, de ned as the original charging curve data minus the least squares estimate, plotted on logarithmic coordinates. A least squares t of a straight line to this data yields an estimate of 1, ^1 = 2:9 msec (the actual value). Once 0 0 = m is known, speci c membrane resistance can be estimated as the ratio Cm = Rm . If both 0 and 1 are known, then the electrotonic length of the cylinder L can be estimated using Eq. 1.109, thus specifying all parameters of the model (including all membrane time constants for n 2). 22 A. B. 0.0 0.5 ln[1V/Vinf] 1.0 ln[1V/Vinf]
0.0 T 0.5 1.0 1.5 2.0 2.0 1.5 2.5 10
0.0 8 6 4 2 0.5 1.0 T 1.5 2.0 Figure 1.7: Illustration of the Rall Peeling Procedure. 1.2.4 TimeDomain Solutions of Finite Length Cable Problems: RC Soma Model
In the previous section, we derived an expression for voltage as a function of both space and time along a passive, uniform, nite length cable subject to the sealedend boundary conditions and an arbitrary initial voltage distribution. In this section, we will construct a more realistic model of a neuron by altering the boundary conditions at one end of the cable. Rather than imposing the sealedend boundary condition at both ends, we will terminate one end using a simple parallel RC circuit. The remaining end will be terminated with a sealedend. The resulting neuronal structure is shown in Fig. 1.8. This model was developed originally by Rall (1969) , and is referred to as either an RCsoma model or the Rall motoneuron model. The model has played a fundamental role in the development of the modern concept of the neuron. The key assumptions of this model are: 23 I0 IL Rs Cs
X=0 X=L Figure 1.8: RCSoma Neuron Model A. the soma is su ciently small that it can be modeled as an isopotential structure with only passive membrane B. the dendritic tree of the neuron being modeled can be represented using a single passive cable with uniform properties. The isopotential assumption A is reasonable given typical major axis length of a typical elliptical neuronal soma (1020 um see Pickard, 1971). The precise conditions under which a dendritic tree can be represented using a single equivalent cylinder will be described later. As always, the assumption of linear membrane properties inherent in these passive membrane models may only be valid about some neighborhood of the neuron's voltage operating point. As before, we seek a solution of the cable equation subject to the relevant boundary and initial conditions: @ 2V (X T ) = V (X T ) + @V (X T ) (1.111) @X 2 @T I0 L = ;G1 @V (X T ) jX =0 L (1.112) dX V (X T0) = E (X ) (1.113) The boundary condition at X = L is straightforward and yields @V (X T ) j (1.114) @X X =L = 0 The boundary condition at X = 0 is somewhat more complicated. In order to express this boundary condition as a condition on V (X T ), it is convenient to switch momentarily to the use of the variables l and t rather than the dimensionless variables X and T . Letting Rs and Cs denote total somatic resistance and capacitance, we obtain r2 @V (x t) j = V (x t) + C @V (x t) j (1.115) s Ri @x x=0 Rs @t x=0 24 Substitution of X = x and T = tm yields (X G1 @V @X T ) jX =0 = V (X T ) + Cs @V (X T ) jX =0 (1.116) Rs @T m Eq. 1.116 can be simpli ed further by making an additional assumption (Rall, 1969). The assumption is that Rm and Cm are the same everywhere in the soma and dendrite (this assumption is known to be invalid for certain neurons see Shelton, 1985). De ne the quantity as follows: (1.117) = Gd = G1 tanh(L) Gs Gs where Gd is the input conductance of the cable (dendrite) at X = 0 and Gs is the soma input conductance. The parameter is known as the dendritic dominance ratio. Substitution of Eq. 1.117 into Eq. 1.116 yields @V (X T ) j V (X T ) + Cs @V (X T ) ]j (1.118) X =0 = tanh(L) @X Rs Gs @X X =0 m Gs However, Rs Gs = 1, and Cs = Cm As = 1 (1.119) Rm Cm R1m As m Gs where As is the soma area. Therefore, the boundary condition at X = 0 can nally be written as @V (X T ) j @V (X T ) ]tanh(L) (1.120) X =0 = V (X T ) + @X @T In order to simplify application of the boundary condition at X = L, it is convenient to write the solution V(X,T)=F(X)G(T) so that F(X) has the form of Eq. 1.37: V (X T ) = Aeju(L;X ) + Be;ju(L;X )]e;(1+u2 )T (1.121) where u = + j , as before. Application of the boundary condition at X = L yields @V (X T ) j ;ju(L;X ) ; Aeju(L;X )]jX =L = 0 (1.122) X =L = ju Be @X As in our previous analysis, u = 0 is a solution, however it eliminates the possibility of arbitrary initial conditions, so we reject it. This leaves A = B as the solution. The resulting form of V(X,T) is V (X T ) = A eju(L;X ) + e;ju(L;X )]e;(1+u2)T (1.123)
25 Application of the boundary condition at X = 0 by substitution of Eq. 1.123 into Eq. 1.120 yields the following equation: ju e;juL ; ejuL] = ;u2 ejuL + e;juL ]tanh(L) (1.124) Solving Eq. 1.124 for the separation constant u is straightforward, but algebraically tedious. The procedure is to let u = + j , apply Euler's indentity, and equate real and imaginary terms. Doing so demonstrates that = 0 and satis es the following transcendental equation: (1.125) sin( L) = ; cos( L)tanh(L) ) tan( L) = ; tanh(L) ( L) L Thus, we obtain the result that V (X T ) = Acos (L ; X )]e;(1+ 2)T (1.126) (1.127) tan( L) = ; tanh(L) ( L) L V (X 0) = E (X ) (1.128) Figure 1.9: Numerical Evaluation of the Equalizing Time Constants for the RCSoma Model Figure 1.9 shows a plot of the function y = tan(x) (solid lines) superimposed on a plot of y = ;mx (dashed line), where 26 m = tanh(L) (1.129) L The abscissa values corresponding to the points of intersection of these two curves (the points xn = n L) specify the equalizing time constants of the RCsoma model. These values may be determined numerically using Newton's method or the method of bisection. Using the symmetry of the roots of Eq. 1.127, we may rewrite the solution as V (X T ) = tan( nL) = ; tanh(L) ( n L) (1.131) L V (X 0) = E (X ) (1.132) Therefore, as in the previous case, the charging curve of the RC soma model is a weighted sum of exponential functions. The rst exponential term of the series has a time constant o = m = Rm Cm . The remaining time constants are determined by the slope m, and are equalizing time constants. Their values are determined by the dendritic dominance ratio and the electrotonic length of the cable L. Two cases are of special interest. First, ! 0 implies that either Gd ! 0 or Gs ! 1. The case Gd ! 0 can occur in several ways: a) the cable length x ! 0 b) speci c membrane resistance Rm in the dendrite ! 1 and c) cable diameter r ! 0. The case Gs ! 1 can occur if the cable is terminated by a short circuit to ground at X = 0. In either case, the slope m of the dashed line in Fig. 1.9 approaches 1 and solutions of Eq. 1.127 are those values of nL such that tan( nL) = 1. Therefore
n n=0 1 X Ancos n (L ; X )]e ;(1+ n 2 )T (1.130) There is one subtlety here, however. When Gd ! 0, the distribution of membrane potential becomes increasingly uniform. Therefore, as described in Sect. 1.2.3, the coe cients An for n 1 vanish, and membrane potential is described by a single exponential term (as it would be for a pure soma, or single compartment model). When, Gs ! 1, all higher order exponential terms contribute to the form of the charging curve (assuming L greater than zero). The second case of interest is when = 1. In this case, Gd ! 1 (cable diameter d ! 1 speci c axial resistance Ri ! 0 speci c membrane resistance Rm ! 0) or Gs ! 0 (the sealedend boundary condition). In either case, the slope m is zero, and the dashed line of Fig. 1.9 is a horizontal line through the origin. The points of intersection of this line with the tangent function in Fig. 1.9 are those for the sealedend boundary condition given previously. The case Gd ! 1 is uninteresting since it assures that voltage is zero everywhere along the cable. 27 = (2n ; 1) 2L (1.133) As described previously in Sect. 1.2.3, the initial condition V (X 0) = E (X ) can be used to solve for the remaining coe cients An in Eq. 1.130. The approach is identical to that shown in Eqs. 1.791.84 that is, both sides of Eq. 1.130 are multiplied by cos k (L ; X )] and integrated with respect to X over the range of 0 to L. The solution is: V (X T ) = 2 1 + sin(2 nL) ];1 L E (X )cos (L ; X )]dX An = L n 2 nL 0 tan( nL) = ; tanh(L) ( n L) L V (X 0) = E (X ) n=0 1 X Ancos n (L ; X )]e Z ;(1+ n 2 )T (1.134) (1.135) (1.136) (1.137) Impulse and Step Response of the RCSoma Model Derivation of the unit impulse and step response functions for the RCsoma model is based on Eqs. 1.1341.137 following the methods described previously in Sect. 1.2.3. Estimating Membrane and Equalizing Time Constants, Electrotonic Length, and Dendritic Dominance In Sect. 1.2.3, we saw that the voltage response of nite length cables with sealedend boundary conditions is completely speci ed by two parameters, m and L. The membrane time constant could be estimated by application of the Rall Peeling Procedure to obtain the decay rate of the rst exponential term governing the shape of model charging curves. Estimation of the rst equalizing time constant 1 and calculation of the ratio 0 yielded 1 an estimate of the electrotonic length L (see Eq. 1.109). There are three parameters specifying the Rall motoneuron model , L, and m. One approach to solving for these parameters uses the assumption that L is su ciently small (less than about .4, see Rall, 1977). In this case, tanh(L) L so that tanh Ad = G1A G (L) (1.138) As s m Thus, is constrained if the surface areas of both the neuron dendrite (Ad) and soma (As) are known. Such information can be obtained from reconstruction of HRP lled cells. Next, since 0 = 0 is a root of Eq. 1.136, the rst exponential term of the model charging curve has time constant m. Therefore, the Rall Peeling Procedure can be used to estimate this term to obtain an estimate of Rm under the assumption that Cm is 1uF=cm2. This leaves one remaining parameter, L. Figure 1.10 shows a plot of the function 28 m /1 Electrotonic Length L Figure 1.10: Plots of the ratio of
m m 1 versus electrotonic length L for the RCsoma model (1.139) versus L for various values of ( = 1 2 5 20 bottom to top). L can be estimated from these curves given (which may be determined from the geometry of the neuron when L is su ciently small) and the rst two time constants m and 1 (which may be estimated using the Rall Peeling Procedure). A more general approach to estimating membrane parameters can be used when both charging curve measurements (yielding estimates of 0 = Rm Cm, 1, and Gin ) and soma membrane surface area are available (As). In this case, is computed as ; Rm = Cm1 0 (1.140) Gin ; 1 = AG (1.141) and L solves
s m 1 = 1 + 2 n where tan( 1L) = ; tanh(L) ( 1L) L
1 (1.142) (1.143) = s0
1 ;1 Note that estimated values of L will be more accurate when L < 1. This follows from the fact that the slope of the curves in Fig. 1.10 are large in this region. Therefore estimation errors in the rst two equalizing time constants map to smaller errors in estimates of L. When L is large, even small errors in estimating the time constant map to large errors in estimated electrotonic length. 29 1.3 Idealized Branched Cable Models
In this section, we will extend the nite cable model of neuronal dendrites developed previously to include simple models of dendritic branching. This will lead to development of what has become known as the "equivalent cylinder conditions" the conditions under which a branched dendritic tree can be represented as a single unbranched cylinder (Rall, 1962 a far more readable treatment is given in the book by Jack et al, 1975). Along the way, we will also develop an iterative algorithm for computing wholecell input resistance, one of the rst techniques used to generate estimates of speci c membrane resistance Rm (Rall, 1959). 1.3.1 Equivalent Cylinder Models: Peripheral Spread of Current
G21 L21 Figure 1.11 shows the simplest model of a branched dendritic tree a primary dendrite Go L11 G11 G22 L22 Figure 1.11: Example of a binary dendritic tree which bifurcates to form two daughter branches. De ne the following quantities: G0 input conductance of the dendritic tree Gjk conductance looking out of branch k of order j G1 input conductance of an in nite extension of the kth branch of order j jk Bjk the quantity Gjk G1 jk
30 Ljk the electrotonic length of branch k of order j Assume that speci c membrane resistance Rm is uniform throughout the dendritic tree (or can be written as a function of a constant and the branch order), and that the terminating conductances G21 and G22 and speci c core resistance Ri are known. Then Eq. 1.47 can be used to iteratively compute the input conductance of the dendritic tree as follows: (1.144) G11 = G1 1B21 + tanh((L21)) ] + G1 1B22 + tanh((L22)) ] 22 + B tanh L 21 + B21tanh L21 22 22 G0 = G1 1B11 + tanh((L11)) ] (1.145) 11 + B tanh L
11 11 It can be seen that G0 is a function of the single unknown quantity Rm (this follows from the de nition of the term Bjk ). Thus, if the input resistance of the dendrite is known through experimental measurement, then an estimate of Rm can be obtained by choosing that value yielding the measured input conductance. The procedure for computing input resistance of a simple binary dendritic tree can be applied in an iterative manner to compute input resistance of an arbitrarily complex binary dendritic tree. This procedure was in fact the rst technique used to estimate Rm in real neurons (Rall, 1959).The accuracy of this estimate of Rm will of course depend on whether or not electrode penetration results in the creation of any signi cant shunt resistance. An equivalent cylinder representation of the branched dendritic tree shown in Fig. 1.11 can be derived in the following fashion. Assume that Rm and Ri are uniform throughout the dendritic tree. Then G11 is given by (1.146) G11 = 2 (Rm Ri);1=2 1B21 + tanh((L21)) ]d3=2 + B21tanh L21 21 + 2 (Rm Ri);1=2 1B22 + tanh((L22)) ]d3=2 (1.147) + B tanh L 22
22 22 where d21 and d22 are the diameters of the second order branches. Assume that B21 = B22 = B2, and L21 = L22 = L2. The rst condition states that the boundary conditions on the second order branches are the same. The second condition asserts that the electrotonic length of the branches are identical. Under these two assumptions, Eq. 1.147 yields G11 = 2 (Rm Ri);1=2 1B2 + tanh((L2)) ](d3=2 + d3=2) (1.148) 22 + B2tanh L2 21 G11 is therefore the input conductance of a cable with electrotonic length L2, diameter 3= d2 = (d212 + d3=2)2=3, and terminating boundary condition B2. It follows that a current 22 injected at the branch point will produce the same voltage along such a cable as it will along either of the two second order branches in the original model. We can therefore 31 Go L11 , d11 G11 L 2 , d2 B2 Figure 1.12: A "collapsed" version of the binary dendritic tree shown previously collapse the model shown in Fig. 1.11 to the simpli ed model shown in Fig. 1.12. Current injected at X = 0 in this model will also produce the same voltage distribution along the two cables as would be observed in the branched model. Next, if d3=2 = d3=2 = (d3=2 + d3=2) (1.149) 11 2 21 22 then the branched dendritic tree of Fig. 1.11 can be collapsed to a single cylinder with electrotonic length L11 + L2, diameter d11, and boundary condition B2, as shown in Fig. 1.13. Go L11 + L2 , d11 B2 Figure 1.13: An equivalent cylinder representation of the binary dendritic tree shown previously This cylinder model of the branched dendritic tree is known as the equivalent cylinder model, and the condition speci ed by Eq. 1.149 is known as the "3/2 power branching rule". In summary, the conditions which must be met in order for a branched binary dendritic tree to have an equivalent cylinder representation are that: A. Rm and Ri are uniform throughout the dendritic tree (this constraint can be relaxed see Jack, Noble, and Tsien, 1975) B. Daughter branches of the same primary dendrite must have identical electrotonic length and boundary conditions C. the 3/2 power rule must hold at each branch point in the dendritic tree 32 1.3.2 Equivalent Cylinder Models: Central Spread of Current From what we have derived so far, the equivalent cylinder has the property that current injection at X = 0 produces the same voltage as would be observed at the same electrotonic distance from the point of current injection in the branched model. Distribution of voltage in response to somatic current injection is the same in branched dendrites and equivalent cylinder models because the models have identical output conductances at every branch point. We will now consider how the distribution of voltage in response to injection of current into the distal dendrites of branched cylinder models compares to that in equivalent cylinder models. To do this, we will consider a specialized case illustrated in Fig. 1.14. We will make the following simplifyingassumptions:
Branch A L21, d21, G21 Branch C L11, d11, G11 I/2 C X=L 2 Branch B L22 , d22 , G22 X=0 I/2 X=L 2 Branch C L11, d11, G11 X=0 Branch A I C L2 , d2 , G2 X=L 2 X=0 Figure 1.14: Central Spread of Current in Branched and Equivalent Cylinder Models A. The branched model satis es the equivalent cylinder conditions 33 B. Sealedend boundary conditions apply on all cylinders C. Branching is symmetric (d21 = d22). Given these assumptions, we may compute the voltage at a point X along branch A in response to symmetric injection of current (magnitude I/2) at the distal tips of branches A and B. To do this, we use the linear superposition principle to calculate the voltage along branch A due to current injection at the end of branch A, and add this to the voltage along branch A due to current injection at the end of branch B. Let VA (X ) be the voltage at location X on branch A due to current injection at the distal end of branch A alone. Then tanh L L21 + L VA (X ) = I G21 1B+ B21tanh((L 21)) cosh(cosh; X ) + B21sinh(L21); X ) (1.150) 2 11 21 + (L21) B21sinh( 21 21 Letting L11 = L1 and L22 = L21 = L2, and de ning the boundary condition B21 as G1 tanh(L2) + G1 tanh(L1) = tanh(L ) + 2tanh(L ) 11 (1.151) B21 = 22 2 1 1 G21 Then L2) + 2 (L1)] ; VA (X ) = I G11 1 + tanh(tanh(L tanhtanh(tanh(L2 ; X ) cosh(L2L )X ) (1.152) 2 11 L2) cosh( 2 1) + In particular, the voltage at branch point C is 1 VC = VA (L2) = I G11 tanh(L ) + tanh(L ) cosh1(L ) (1.153) 2 11 1 2 2 The voltage at location X on branch A due to current injection at the distal end of branch B alone is cosh(X (1.154) VB (X ) = VA (L2) cosh(L )) 2 Using the linear superposition principle, the voltage at location X on either branch A or B due to both current sources is therefore 1 V (X ) = I G11 tanh(L ) + tanh(L ) cosh1(L ) 2 11 2 1 2 cosh(X cosh(L2 ; X ) + tanh(L2) + 2tanh(L1)]sinh(L2 ; X ) + cosh(L )) ](1.155) After application of the appropriate hyperbolic function identities, this reduces to tanh( tanh(L ; V (X ) = I G11 1 +tanh(LL1)+ tanh(2L;)X ) cosh(L2L )X ) (1.156) cosh( 2 1) 2 11 34
2 Next, we will compute voltage at position X along the equivalent cylinder model in response to injection of current with magnitude I at the end of branch A. B tanh L L2 ; X VA (X ) = I G11 1B+ + 1tanh((L 2)) cosh(cosh(L ) + B1sinh(L2); X ) (1.157) 2 2 ) + B1sinh(L2 11 1 where 1 Gin 11 (1.158) B1 = G1 = G11tanh(L11) = tanh(L11) G1 11 11 Incorporation of Eq. 1.158 into Eq. 1.157, and simpli cation, yields tanh( tanh(L ; V (X ) = I G11 1 +tanh(LL1)+ tanh(2L;)X ) cosh(L2L )X ) (1.159) cosh( 2 1) 2 11 The agreement of Eqs. 1.156 and 1.159 shows that current injection, as illustrated in Fig. 1.14, produces the same voltage at equal electrotonic locations in the two models. Other conclusions concerning central spread of current in idealized branched models and their equivalent cylinder representations are described in (Jack, Noble and Tsien, 1975). 1.3.3 Examples of Peripheral and Central Current Spread A very simple example can be used to highlight the di erence in peripheral versus central current spread. Assume an idealized neuron model consisting of a passive cable terminated at X=0 by an RC soma and at X=L by the sealedend boundary condition. Let the physical length and diameter of the cable be 200 and .2 um, respectively. Let the soma diameter be 20 um. Further, assume a speci c core resistance value (Ri) of 100 ; cm, a dendritic speci c membrane resistance Rm d of 100,000 ;cm2 , a somatic speci c membrane resistance of 10,000 ; cm2, and speci c membrane capacitance of 1uF=cm2. While these values of speci c membrane parameters may seem extreme, they are comparable to somatic Rm values estimated for the cerebellar Purkinje neuron by Shelton , and to dendritic Rm values estimated in this same cell (Shelton, 1985) and in the retinal ganglion cell (Coleman and Miller, 1989). Under these assumptions, the dendritic space constant is 707 um. The dendrite therefore has an electrotonic length of .3 (comparable to the electrotonic length of dendrites in mudpuppy retinal ganglion cells), and the wholecell input resistance is 200M . Calculations show that the dendritic dominance ratio = 6x10;4 . Therefore, the charging curve in response to current injection at the soma is indistinguishable from a pure exponential, with a time constant of 10 msec. One of the most misunderstood aspects of linear cable theory is the notion that measurement of an exponential charging curve in response to current injection at the soma 35 4.0 V(x)/IR
Bi I
Bs Bi P 3.0
Gp Soma Pb Bc1 2.0
Bc2 P Bs 1.0
Gp Pb Soma .25 .5 .75 Bc1 Bc2 1.0 X Figure 1.15: Steadystate voltage values in a branched dendritic tree obeying the 3/2 power rule (with symmetric branching) in response to current injection at point Bi. Abscissa is electrotonic distance from the soma. Ordinate is voltage response scaled by that measured in a semi in nite cylinder with membrane properties identical to the fourth order branch. 36 implies that the neuron in question can be modeled using a single isopotential compartment. The error of this idea can be seen by computing the distribution of voltage along the dendritic cable of the above example in response to current injection at the denritic tip. Calculations show that the ratio V ((X =:3) , voltage at the site of current injection divided V X =0) by the soma voltage, is .007. In other words, there is an extremely large voltage gradient along the cable. When viewed from the dendritic end, the cable appears to be terminated by a "shortcircuit" to ground (the soma). The large voltage gradient is therefore due to the voltage drop produced by axial current ow through the cable core conductance. Thus, despite the fact that the model neuron exhibits an exponential charging curve, it most de nitely may not be modeled using a single isopotential compartment. An additional example rst shown by Rall (1973, Fig. 4) demonstrates a functional consequence of the 3/2 power branching rule. Rall's data are reproduced in Fig. 1.15. The inset to this gure shows a branched neuron model with four orders of branching. Branching is assumed to obey the 3/2 power rule, and is symmetric (that is, dx1 = dx2). Dendrites are terminated with a sealedend boundary condition, and all branches have an electrotonic length of .25. The ordinate shows scaled voltage in response to current injection at the end of a single fourth order branch at the point marked Bi as a function of distance from the soma (abscissa). Voltage is normalized by the term IR1, the voltage that would be observed in response to current injection into a semiin nite version of the fourth order branch. The ordinate values are therefore dimensionless, and independent of the particular value of I. These data show that there is a very large voltage drop between the site of current injection (point Bi) and the third order branch point (point P). This voltage drop is due primarily to current ow through the core resistance of the fourth order branch. The line with endpoints labeled P and Bs shows voltage decay along the sibling fourth order branch which emerges from the branch point P. There is very little voltage drop along this branch. Put another way, very little current ows into the sibling fourth order branch. The reason for this is that the input impedance of this branch is high relative to that of the third order parent branch at point P because the fourth order branch has a smaller diameter than the third order branch and is terminated with a sealedend. Similarly, these data show that at every lower order branch point, charge is preferentially distributed towards the soma. This is an inevitable consequence of the 3/2 power branching rule the diameter of the parent branch is always greater than that of the daughter branches. Therefore, current injected into the dendrites tends to be diverted towards the soma. Neurons with dendritic trees obeying the 3/2 power branching rule are well suited for funneling charge generated by synaptic inputs to the soma. Surprisingly, there have been relatively few studies which have critically examined the branching patterns of particular classes of neurons. The study by Lux et al (1970) has 3=2 shown that the ratio d3=2 dx 1d3=2 has a value of about one in motoneurons. The same +
x+1 1 x+1 2 37 is true for retinal ganglion cells of the mudpuppy (R. F. Miller, unpublished results). Thus, the dendritic trees of these cells appear to be well suited for funneling charge to the soma. A very notable exception to the 3/2 power branching rule has been studied by Miller and Bloom eld (1983). They demonstrated that retinal "starburst" amacrine cells have dendrites which are very thin near the soma and which become larger at higher branch orders. Cable analyses showed that rather than funneling charge to the soma, the dendritic architecture of these cells appeared to be specialized for distributing charge within the dendritic tree. This was consistent with older hypotheses that the soma of these cells plays only a nutrative, maintenance role, with local signal processing occuring within the dendrites. Therefore, these two cells (amacrine and ganglion), both elements of the same tissue, are perfect examples illustrating the role play of di erences in dendritic geometry in determining dendritic sugnal processing. 1.3.4 Estimating the Number of Synaptic Inputs to Equivalent Cylinder Models One surprising consequence of analyzing current spread in idealized branched dendritic trees is that it is impossible to decide the number of sites at which current is injected (or synaptic input occurs) within a dendritic tree unless recordings can be made within the di erent branches of the tree. This idea is illustrated schematically in Fig. 1.16. This gure shows the case where a current I is injected into the distal extremity of a single branch, versus the case where the current is divided equally and injected into the distal dendrites of both second order branches. The linear superposition principle can be used to compute voltage at the branch point for each of the models. In both cases, the same voltage is produced at the branch point. Therefore, the voltage distribution along the rst order branch is identical in the two cases (although voltage distribution is di erent within the second order branches). We therefore conclude that the number of synaptic inputs within the second order branch cannot be determined by measuring electrical events in the rst order branch. 1.3.5 Estimating the Site of Synaptic Inputs to Equivalent Cylinder Models
Under certain, rather ideal conditions, it is possible to estimate the elctrotonic distance from the soma to sites of synaptic input by measuring excitatory postsynaptic potentials (EPSP's) at the soma. In order to understand the basis of this technique, we will rst consider the form of the voltage response of a semiin nite cable when distance between the recording and current injection site is increased. The data in Fig. 1.17 are reproduced from Figure 3.17 of Jack et al (1975). 38 I
V1 V2 I/2
V1 V2 I/2 = +
V3 2V1 V3 > 2V2 2V3 I/2
V2 V1 I/2
V3 V2 = I/2 +
V3 I/2
V1 V1+V3 > 2V2 V1+V3 Figure 1.16: Numbers of synaptic inputs to second order branches cannot be estimated by recording voltage responses in lower order branches. 39 Figure 1.17: Potential change at varying electrotonic locations in response to square wave m current of duration 10 applied at X=0 40 These data show plots of V(X) at X=0, .2, .4, ... in response to a rectangular current m pulse with duration 10 applied at X=0 of a semiin nite cable. There are four important features of these data. First, the amplitude of the voltage response decreases with increasing distance from the site of current injection.This is due to leakage of charge through the dendritic membrane. Second, there is an onset delay which increases with increasing distance from the injection site. Since the voltage across the membrane capacitance cannot change instantaneously, this delay represents the nite charging time of the membrane capacitance. Third, since axial current is decreasing with increasing distance from the site of injection, and since the magnitude of this current determines the rate of charging of the membrane capacitance and therefore membrane potential, risetime slows with increasing distance from the site of current injection. Fourth, the width of the voltage transient increases with increasing distance from the soma. This occurs because at su ciently long times after o set of the stimulus, membrane potential has become more equalized through electrotonic spread of current. Potential decay is therefore less in uenced by the smaller equalizing time constants, and is more heavily determined by the larger membrane time constant. Figure 1.18: Risetime and halfwidth plots for alpha function current injection in a semiin nite cable for di erent values. In principle, all four response characteristics (decay of peak response, delay time, risetime, response width) could be used to estimate the electrotonic distance between the input and output sites. However, use of the relative amplitude of peak response requires that we know voltages at both the site of current injection and the recording site information that is not usually available. Use of the delay time (or travel time along the cable) also requires 41 knowledge of the exact time at which the synaptic input is activated again, information not readily available. This leaves the remaining two response characteristics as criteria for determining the location of single synaptic inputs, since they are independent of knowledge of precisely when the synaptic input occurs or the amplitude of the EPSP at the site of input. Figure 1.18 shows plots of both the rise time and response half width as a function of distance between the recording and current injection sites for the semiin nite cable model described above. In this case, the injected current had the form I (T ) = 2Te; T (1.160) This function reaches its peak value at T = 1= , and approximates the shape of real synaptic currents fairly accurately. Risetime of response voltage was de ned as the time (normalized by the membrane time constant) required to go from 10% to 90% of the voltage transient. Halfwidth was de ned as the time between succesive positive and negative crossings of a threshold voltage set half way between the baseline response and the peak voltage. The data in Fig. 1.18 show that both risetime and half width increase with increasing distance from the site of current injection. Both are reasonably independent of the duration of the current pulse (as re ected by di erent values of ). Similar computations may be performed using the Rall motoneuron model described in Sect. 1.2.4. The results are shown in Figures 7.19 (halfwidth) and 7.20 (rise time) of Jack et al (1975). However, in this case, both response measures are functions of m, L, and . Jack et al (1975) have shown that the dependence on is weak, but that uncertainty in m and L can signi cantly e ect estimates of input location. These results indicate that if: a) a neuron conforms to the equivalent cylinder model b) if L and m are known c) it is known a priori that the neuron receives temporally coincident synaptic input within a highly restricted region of the dendritic tree then the electrotonic location of synaptic inputs can be estimated by measuring EPSP risetime and half width. That this is the case has been demonstrated very elegantly by Redman and Walmsley (1983). These investigators electrically stimulated individual Ia a erents to motoneurons, and recorded EPSP's at the soma. EPSP's were averaged over many trials, and the risetime and halfwidths estimated. The motoneuron membrane time constant was also measured by "peeling" charging curves in response to current injection at the soma and by estimating the nal EPSP time constant of decay. Speci c axial resistance was set to 70 ; cm from previous work, and membrane capacitance to the standard value of 1 uF/cm2. Both the Ia a erents and the motoneuron were HRP lled and reconstructed. Once the motoneuron branch diameters were known, it was possible to compute the electrotonic length of each branch, and to compute the electrotonic distance of the Ia inputs from the soma. In addition, an equivalent cylinder model of the neuron was constructed, 42 single function current inputs were applied at various locations in the denritic tree, and the location yielding the best t to measured data determined. The results of this study showed that in those cases where the anatomy indicated that excitatory projections were highly localized, there was excellent agreeement between estimates of synaptic location based on the morphological data and the electrophysiological estimates. Single  function current inputs were not capable of reproducing EPSP data for those neurons in which synaptic inputs were shown to be spatially distributed, or when localized inputs were not activated at the same time. In light of all the highly restrictive assumptions which must be assumed in the absence of morphological data, halfwidth analyses have very limited application to the problem of estimating synaptic location. 43 2 Compartmental Modeling
So far, we have considered highly simpli ed models of neurons which take the form of partial di erential equations. We will now replace these distributed parameter models with spatially discrete versions, or lumped parameter models, known as compartmental models. The theory of compartmental modeling was developed early in 1964 by Rall. However, it saw very limited application for the next fteen to twenty years. The reasons for this were twofold. The rst and primary reason was that compartmental modeling by its very nature requires a complete description of the branching pattern of the neuron under study, including measurements of dendritic lengths and diameters. This kind of information has become widely available only recently with the development of HRP and ourescence staining techniques. The second reason for the limited initial interest was that numerical integration of the resulting systems of di erential equations was, and still is, a computationally demanding task. The computing resources needed to solve such problems were not as commonly available to researchers in 1964 as they are today. Our analyses of compartmental models will be restricted once again to consideration of passive models, that is, models in which there are no voltagegated ion channels in the cell membrane. The analyses fall naturally into two domains: a) models in which currents are applied as inputs and voltage is the output (linear models) and b) models in which synaptic conductance changes are inputs and voltages are outputs (such models are inherently nonlinear). 2.1 Subdivision of Dendritic Trees into Discrete Compartments
The rst step in developing a compartmental model of a cell or neuron is to decompose the dendritic tree into small compartments within which spatial nonuniformities of all state 44 variables (including not just membrane potential, but possibly ionic concentrations) can be neglected. Thus, nonuniformity of state variables is represented only in terms of di erences between, and not within, compartments. For simplicity, we will illustrate the principles of compartmental modeling by restricting our attention to an unbranched passive cable in which there is a single state variable associated with each compartment (membrane potential). Compartmental modeling methods are easily generalized to more complicated binary dendritic trees and systems with a greater diversity of state variables (including active channels). A cable, and its decomposition into a number of isopotential compartments was shown previously in Fig.1.2. The rst question that must be answered is "how short must the segments be in order for the isopotential assumption to hold?" This question can be answered by applying what we have learned so far about cable theory. Assume that each dendritic segment is .05 space constants in length (L=.05). Note that in order to know the physical length that corresponds to this electrotonic length, we must know Rm, Ri, and the dendritic diameter d. Compartmental models can only be constructed for those cells for which we have complete morphological data, so that we can assume the dendritic diameter is known for each segment. The range of Ri values is quite restricted, and is usually assumed to be about 110 ; cm (see Shelton, 1985, for an excellent discussion of experimental and theoretical estimates of Ri). The remaining parameter Rm usually is not known apriori. It is therefore necessary to assume a range of values for Rm . Given an electrotonic length per segment we can compute the voltage decay along that segment in response to current injection at one end subject to particular boundary conditions using Eq. 1.49: V (X ) = cosh(L ; X ) + Gout sinh(L ; X ) G1 (2.1) V (0) cosh(L) + Gout sinh(L) G1 Assuming the sealed end boundary condition, voltage at the end of the cylinder in response to current injection at the opposite end decays by only .1% the segment is essentially isopotential. Decay depends on the particular boundary conditions assumed. For example, if we assume that each segment is terminated by a conductance G1, then potential decays by 5% over the length of the cable. The actual boundary conditions are of course determined by the structure of the neuron being modeled. It has been determined empirically that as long as the electrotonic length of each segment is less than about .1 (using the sealedend condition), then numerical solutions of cable problems obtained using compartmental models agree well with analytic solutions. Note that for systems in which there are synaptic inputs, care must be taken in to choose L :1 when these inputs are maximally activated. 45 2.2 Membrane Models
Membrane segments are modeled as described by Rall (1964), and as shown in Fig. 2.1.
gj gm j cj
e i Vj (t) g j(t)
e gj(t) i Er Ej Ej Figure 2.1: Electrical model of an isopotential dendritic compartment. The values of circuit elements are determined by the radius rj and length lj of the j th compartment, and are de ned as follows: gj is axial conductance of the j th compartment, gj = 2 rj lj Gi cj is total membrane capacitance of the j th compartment, cj = 2 rj lj Cm gjm is total passive membrane conductance, gjm = 2 rj lj Gm Er is resting membrane potential in the absence of any synaptic input Eje is the reversal potential of the ion(s) whose membrane permeability is increased during excitatory synaptic input to compartment j. Its value is given by the Nernst potential for single ions, or by the GoldmanHodgkinKatz equation for multiple species. Eji is the reversal potential of the ion(s) whose membrane permeability is increased during inhibitory synaptic input to compartment j. Its value is determined by the Nernst or GHK equations. gje (t) is the membrane conductance change produced by excitatory synaptic input (Ee > Er )
46 gji (t) is the membrane conductance change produced by inhibitory synaptic input (Ei Er ) The functional form of the synaptic conductance changes must be assumed. Once again, the socalled function is used, so that: t (2.2) gje i(t) = ( tt )e(1; tp )Ge i j Ge i j
p = Ge i De i2 rj lj (2.3) In the above equations, tp is the time at which the synaptic conductance reaches a peak value, Ge i is the peak conductance change due to excitatory or inhibitory synaptic input, j Ge i is the per channel conductance, De i is channel density in units of channels per um2, rj and lj are the segment radius and length in um, respectively. 2.3 Compartmental Model State Equations in Response to Current Inputs
For simplicity, we will begin by deriving state equations for a compartmental model of a uniform, passive unbranched cylinder in which the input is assumed to be injected current. Thus, we will ignore complications posed by the timevarying conductances shown in Fig. 2.1 which model synaptic input. An electrical model of an unbranched cable is shown in Fig. 2.2. The elements gj and Vj are axial conductance and voltage of the j th comg j1 Vj1 z j1 gj Vj zj g j+1 Vj+1 z j+1 Figure 2.2: Electrical model of a dendritic cable. partment. The term zj represents the parallel combination of the membrane capacitance and resting membrane conductance for the j th compartment (zero synaptic conductance). State equations for the j th compartment are derived by applying Kircho 's current law at each node. For a linear cable partitioned into k compartments, this yields: j( vj;1(t);vj (t)]gj + vj+1(t);vj (t)]gj+1+Ij (t) = cj dvdt t) + vj (t);Er]gjm 1 < j < k(2.4) 47 1( m v2(t) ; v1(t)]g2 + I1(t) = c1 dvdt t) + v1(t) ; Er ]g1 j = 1 (2.5) k m vk;1(t) ; vk (t)]gk + Ik (t) = ck dvdt(t) + vk (t) ; Er ]gk j = k (2.6) where Ij (t) is current injected into the j th compartment. De ne a normalized voltage vj (t): vj (t) = vj (t) ; Er (2.7) so that all voltages are expressed relative to the resting potential. This yields dvj = vj;1 gj ; vj gjm + gj + gj+1] c1 + vj+1 gjc+1 + Ijc(t) 1 < j < k (2.8) dt cj j j j dv1 = ;v gm + g ] 1 + v g2 + I1(t) j = 1 (2.9) 2c 1 1 2 c1 dt c1 1 dvk = v gk ; v gm + g ] 1 + Ik (t) j = k (2.10) k k ;1 c k k dt ck ck k These equations can be written in matrix form. De ne v (t) = v1(t) v2 (t) : : : vk (t)]T (2.11) (2.12) u(t) = I1(t) I2(t) : : : Ik (t)]T B = diag 1=c1 1=c2 : : : 1=ck ] (2.13) where diag ] denotes a diagonal matrix, in this case, with dimension k x k. Note that B is a fullrank, invertible matrix. Finally, de ne the k x k matrix A as ; c11 g12 0 0 0 : : : 0 c g3 g2 2 ; c2 c2 0 0 : : : 0 c2 g3 ; 3 g4 0 : 0 : : 0 c3 c3 c3 0 : : : : : : 0 A= 0 : : : : : : : : : gk;1 ; k;1 gk 0 0 0 0 0 0 ck;1 ck;1 ck;1 gk ; ckk 0 0 0 0 0 0 0 ck where m 1<j<k (2.14) j = gj + gj +1 + gj m (2.15) 1 = g2 + g1 m (2.16) k = gk + gk Then 48 v (t) = Av (t) + Bu(t) _ {z  velocity field } (2.17) Equation 2.17 has the form of a timeinvariant linear system with k x 1 state vector v (t), k x k state transition matrix A, and k x 1 input vector u(t). Elements of the matrix A are determined by the connectivity of the compartments modeling the neuron. Put another way, A contains information on the geometry of the neuron. For the case of an unbranched cylinder, A is a tri diagonal matrix. This follows from the fact that in this case, every compartment connects at most to a right and a lefthand neighbor. For a general binary dendritic tree, there will be no more than one additional entry in each row of A. The location of this entry will depend on the ordering of the compartments within the dendritic tree. The nonzero elements of the input vector u(t) correspond to compartments receiving current injection. We will now investigate stability of Eq. 2.17 when the input vector u(t) = 0 (state stability). We will do this to obtain a better understanding of the general dynamic properties of compartmental models. We are speci cally interested in knowing things like how many critical (equilibrium) points such models have, and whether or not behavior such as oscillation can be observed. The zero input assumption yields v (t) = Av (t) _ (2.18) If all the eigenvalues of the matrix A have negative real part, then the system de ned by Eq. 2.18 has a single critical, or equilibrium point which is a stable global attractor. Thus, our task is to prove that A has eigenvalues with negative real part. In order to do this, it is useful to represent A in terms of a new basis set in which the representation of A is a real symmetric matrix, since there are a number of useful theorems describing properties of eigenvalues/eigenvectors of such matrices. De ne the k x k fullrank diagonal matrix M as
; ; M = diag c1 2 c; 2 : : : ck 2 ] (2.19) 2 Next, express the state vector v (t) as a linear combination of the columns of the matrix M. In other words, the set of k, k x 1 vectors which are the columns of M = m1 m2 : : : mk ] form a new basis in which the representation of v (t) is y(t): v (t) = My(t) ) y(t) = M ;1v (t) (2.20) It follows that the representation of Eq. 2.18 in this new basis is
1 1 1 2.3.1 State Stability of the Homogeneous Equation 49 y(t) = M ;1 AMy(t) _ (2.21) Since M is a fullrank matrix, the systems de ned by Eqs.2.18 and 2.21 are equivalent in that there is a 1:1 mapping between trajectories of the two systems. In addition, the matrices A and M ;1 AM are related by what is known as a similarity transformation. Why is this important? Theorem 2.1 Let M be a fullrank (invertible) k x k matrix. Let A be a k x k matrix. Then the matrices A and M ;1AM have identical eigenvalues Let y i = M ;1 xi (the representation of the eigenvalues of A in the basis M). Then (2.23) (A ; iI )Myi = 0 ;1 (A ; I )My = 0 (2.24) M i i ;1 AMy = M (2.25) iyi i Therefore, yi and i are eigenvectors and eigenvalues of M ;1 AM . Proof 2.1 Let i and xi be the eigenvalues and corresponding eigenvectors of A so that (A ; i I )xi = 0 (2.22) The advantage of this similarity transformation is that it yields a basis in which the state transition matrix M ;1 AM of Eq. 2.21 is a real symmetric matrix: ; c11 pcg12c2 0 0 0 : : : 0 g2 2 pg3 pc1 c2 ; c2 0 0 : : : 0 c2 c3 g3 3 pg4 pc2 c3 ; c3 0 : : 0 c3 c4 0 : ;1 AM = 0 M 0 : : : : : : 0 : : : : : : : : : gk;1 k;1 p gk 0 0 0 0 0 0 pck;2 ck;1 ; ck;1 ck;1ck k pckg;1 ck 0 0 0 0 0 0 0 ; ckk Real symmetric matrices have real eigenvalues. Since A and M ;1AM have identical eigenvalues, it follows that A has real eigenvalues. If we can now show that all the real eigenvalues of A are negative, then we will have demonstrated that Eq. 2.18 is stable. To prove this, we will use Gershgorin's Circle Theorem: Theorem 2.2 Every eigenvalue of a matrix A with elements ajk lies in at least one of the discs in the complex plane with center ajj and radius rj , where rj = j ajk j (2.26)
k6=j X 50 Proof 2.2 See pgs. 263264 of Quinney, 1985 .
Applying Gershgorin's Circle Theorem to the matrix A, we obtain (2.27) ajj = ;(gj + gj+1 + gjm ) c1 1 < j < k j rj = gj + gj+1 1 < j < k (2.28) cj m a11 = ;(g2c+ g1 ) (2.29) 1 (2.30) r1 = g2 c1 m ;(gk + gk ) (2.31) akk = ck rk = gk (2.32) ck We already know that the eigenvalues of A are purely real. Using the above result, we can now show that every eigenvalue of A must lie in at least one of the intervals dj , where gm ;gjm (2.33) dj = ;2(gj c+ gj+1 ) ; cj cj ] 1 < j < k j j m m d1 = ;c2g2 ; gc1 ;cg1 ] (2.34) 1 1 1 m m (2.35) dk = ;c2gk ; gk ;cgk ] ck k k Since each compartment has nonzero surface area, each one of these intervals is bounded away from zero. Therefore, since every interval is strictly less than zero, every eigenvalue must be less than zero. We conclude that since the eigenvalues of A are real and negative, the homogeneous system given by Eq. 2.18 has a single globally attracting stable critical point at v (t) = 0. Since the system is linear, this critical point is also exponentially stable. In other words, when started at any point in the state space, the system state vector decays exponentially to the zero state. The system does not exhibit more complicated motion such as periodic oscillation (unless the input is itself oscillatory). The steadystate solution v to maintained current input u is trivially obtained as v = ;A;1Bu (2.36) 51 2.3.2 SteadyState Response to Maintained Current Inputs Additional quantities such as pointtopoint DC transfer functions (the ratio of current applied in compartment j to voltage in compartment i) and attenuation factors (the ratio of voltage in compartment j to voltage in compartment i) can be derived in a straightforward fashion from the above equation (Perkel and Mulloney, 1978a Perkel and Mulloney, 1978b). A standard result from linear system theory is that the solution to Eq. 2.17 is 2.3.3 Transient Response to Current Inputs Z t eA(t; )Bu( )d v (t) =  {z 0 + t } {z } Unforced Response  Forced Response
eA(t;t0)v
0 (2.37) In this section, we will show that responses of compartmental models to arbitrary current inputs can be expressed as a linear combination of input functions convolved with a nite number of decaying exponential terms (as opposed to an in nite number of terms for the cable equation). For simplicity, we can set v0 = 0 and t0 = 0 without loss of generality. Then v(t) is given by the forced response in Eq. 2.37. Thus, an analytic solution form, the well known convolution integral (but in a multidimensional system), exists. The question is "What is eAt?" Formally, eAt has a power series expansion: 2 3 eAt = I + At + A t2 + A t3 + : : : (2.38) 2! 3! However, in practical terms, it becomes di cult to evaluate such an in nite series unless the matrix A is diagonal. Our problem will be simpli ed if we choose a new basis set in which the representation of the state transition matrix A is a diagonal matrix ^. We previously de ned an equivalent, or conjugate system to that given in Eq. 2.17: y(t) = M ;1AMy(t) + M ;1 Bu(t) _ (2.39) v (t) = My(t) We know that M ;1AM is a real symmetric matrix. It therefore has a set of k orthonormal eigenvectors zi, corresponding to eigenvalues i. De ne a socalled modal matrix D of M ;1 AM : (2.40) D = d1 d2 : : : dk ] Express the vector y(t) using D as the new basis set: y(t) = Dr(t) (2.41) 52 Then substitution of Eq. 2.41 into Eq. 2.39 yields r(t) = DT (M ;1 AM )Dr(t) + (DT M ;1B )u(t) _ (2.42) = ^r(t) + Fu(t) (2.43) where ^ = diag 1 2 : : : k ] (2.44) T M ;1 B F = D (2.45) T = D;1 . In this new basis set, the system de ned by Eq.2.43 is uncoupled. and D Theorem 2.3 Let B be a k x k matrix with eigenvectors di and eigenvalues i . De ne the modal matrix D to be the matrix with columns equal to the eigenvectors of B. Then DT BD is a diagonal matrix with diagonal elements equal to the eigenvalues i. Proof 2.3 Let D be the modal matrix for B. Then (2.46) BD = B d1 d2 : : : : dk ] = 1d1 2d2 : : : k dk ] (2.47) dT 1 . 1d : : : k d ] T BD = DT d . D (2.48) 1 1 2 d2 : : : k dk ] = . 1 k T dk = diag 1 : : : k ] (2.49) We can now write the solution to Eq. 2.43 (assuming r0 = 0 and t0 = 0): t r(t) = 0 e^(t; )Fu( )d (2.50) where e^t = diag e 1t e 2t : : : e kt] (2.51) y(t) = Dr(t) v (t) = My(t) ) v (t) = MDr(t) (2.52) Inspection of Eq. 2.50 shows that the product of e^(t; )F is the matrix F with the ith row multiplied by e i(t; ). The ith entry of the solution vector therefore has the general form t (2.53) ri(t) = Fij e i(t; )uj ( )d Z X Z
j Equation 2.53 demonstrates an important result. It shows that the response of a passive compartmental model to arbitrary current inputs is a weighted combination of exponentials convolved with the current inputs. The time constants of the exponentials are given by the eigenvalues of the original state transition matrix A. If the inputs are reasonably simple in form, then closed form solutions for the convolution integral are readily obtained, and an exact response can be computed (as opposed to direct numerical integration of the underlying state equations). 53 0 Set uj (t) = 0 8j 6= i, and ui(t) = (t). The current input is therefore a Diracdelta current pulse applied to compartment i. For simplicity, set the initial condition r0 = 0 and t0 = 0. Then substitution of the above expression for u(t) into Eq. 2.43 yields the impulse response k(t): e 1tF1i k(t) = MD ... (2.54) t Fki ek The response of the compartmental model to a Diracdelta current input applied to the ith compartment is therefore a weighted sum of k exponentials, with the exponentials being the eigenvalues of the original state transition matrix A. A similar approach can be used to obtain the unit step response. 2.3.4 Example: Compartmental Model Impulse Response 2.4 Compartmental Model State Equations in Response to Synaptic Inputs
In the previous section, we saw how the state equations for compartmental models subject to current input are those of a timeinvariant multi dimensional linear system. In this section, we will examine these same models when inputs are synaptic conductance changes rather than injected currents. We will show that in this case, the resulting state equations are no longer linear, but rather, are members of a class of nonlinear systems known as bilinear systems. There are no useful analytic solutions for such state equations, so they can only be solved using numerical techniques. Once again, we will derive the state equations for a compartmental model of an unbranched dendritic cylinder as shown in Fig. 2.2. However, in this case the term zj represents the parallel combination of the membrane capacitance, resting membrane conductance, and synaptic conductances as shown in the membrane model in Fig. 2.1. State equations for the j th compartment are derived by applying Kircho 's current law at each node. This yields: j( vj;1(t) ; vj (t)]gj + vj+1(t) ; vj (t)]gj+1 = cj dvdt t) + vj (t) ; Er ]gjm (2.55) + vj (t) ; Eje ]gje(t) + vj (t) ; Eji ]gji (t) (2.56) 1( m e e i i v2(t) ; v1(t)]g2 = c1 dvdt t) + v1(t) ; Er ]g1 + v1(t) ; E1 ]g1 (t)+ v1(t) ; E1]g1(t)(2.57) k m e e i i vk;1(t);vk(t)]gk = ck dvdt(t) + vk(t);Er]gk + vk(t);Ek]gk (t)+ vk(t);Ek]gk (t)(2.58) 54 De ne a normalized voltage vj (t): j ;E vj (t) = vE(et); E i r (2.59) j j The numerator term is deviation of membrane potential in the j th compartment, vj (t), from the resting potential Er . The denominator term is the maximum voltage change that can be observed in any compartment. Values of vj (t) therefore lie in the interval (;1 1). Then the state equations reduce to dvj = vj;1 gj ; vj (gjm + gje(t) + gji (t)) + gj + gj+1 ] c1 + vj+1 gjc+1 (2.60) dt cj j j Eje ; Er gje(t) Eji ; Er gji (t) 1<j <k (2.61) + Ee ; Ei c + Ee ; Ei c j j j j j j e e i i dv1 = ;v (gm +ge (t)+gi (t))+g ] 1 + v g2 + E1 ; Er g1 (t) + E1 ; Er g1(t) (2.62) 2c 1 1 1 1 2 c1 E e ; E i c1 e i dt E1 ; E1 c1 1 1 1
e e i i dvk = v gk ; v (gm+ge (t)+gi (t))+g ] 1 + Ek ; Er gk (t) + Ek ; Er gk (t) (2.63) k k;1 c k k k k e i e i dt ck Ek ; Ek ck Ek ; Ek ck k These equations may once again be written in matrix form. Let v (t) be a k x 1 state vector with elements equal to the normalized voltage in each of the k compartments of the model. Let u(t) be a 2k x 1 input vector with elements the excitatory and inhibitory synaptic inputs to each compartment: e i e i e i u(t) = g1 (t) g1(t) g2 (t) g2(t) : : :gk (t) gk (t)]T (2.64) Let A be the k x k matrix de ned previously in Eq. 2.3, and de ne a k x k diagonal matrix G(t) as e + i e + i (2.65) G(t) = diag ; g1(t) c g1(t) : : : ; gk (t) c gk (t) ] 1 k Finally, de ne a k x 2k matrix B as 0 B= 0 0 0
c1 E1 ;E1
e 1 E1 ;Er e i c1 E1 ;E1 i 1 E1 ;Eri e 0 0 0 0 c2 E2 ;E2 e 1 E2 ;Er e i 0 0 0 0 c2 E2 ;E2 i 1 E2 ;Eri e 0 0 0 0 0 0 : 0 0 : : : 0 0 : : 0 : : 0 : : 0 : : 0 e i Ek ;Er 1 Ek ;Er 0 c1k Eke;Eki ck Ek ;Eki e 55 Then the state equations describing the time rate of change of voltage along the cable can be written as v (t) = A + G(t)]v (t) + Bu(t) _ (2.66) In the above equation, A is the same matrix which was de ned previously as containing information about the topology of the neuron being modeled. The matrix G(t) represents perturbations applied to the diagonal terms of A by synaptic conductance changes. The state transition matrix of Eq. 2.66 is therefore timevarying. This in and of itself is not a problem there is a very well developed theory of the dynamics of linear, timevarying multidimensional systems. The di culty stems from the fact that these time varying perturbations to the state transition matrix are caused by the input vector u(t). This can be made clear by rewriting the above equation: v (t) = Av (t) + Bu(t) + G(t)v (t) _ (2.67) =  {z } constant parameter linear system
Av (t) + Bu(t) + In this equation, element njj , the j th diagonal element of the k x k diagonal matrix Ni, is i c;1 if the j th compartment receives synaptic input, and is zero otherwise. Therefore each j matrix Ni corresponds to one of the 2k synaptic inputs and has a single non zero element marking the compartment which receives this input. The term ui(t) is the ith element of the input vector u(t). The equation makes it clear that if any of the Ni are nonzero, then there is a multiplicative interaction between the input and the state vector. The consequence of this interaction is that the state equation is nonlinear. We will prove this assertion by contradiction. Suppose that v0 and v1 are solutions of Eq. 2.67 corresponding to inputs u0 and u1. Then the system is linear i v = 0v0 + 1v1 is a solution in response to the input u = 0u0 + 1u1. Let u0i and u1i denote the ith elements of vectors u0 and u1, respectively. Then it follows that 2k X Niv (t)ui(t) i=1 {z terms } nonlinear (2.68) vo + v1 = Av0 + Bu0 + _ _
i=1 2k 2k X Niv0u0i + Av1 + Bu1 + X Niv1u1i i=1 i=1 2k 2k X X + Niv0u1i + Niv1u0i
i=1 (2.69) (2.70) Since v0 and v1 are solutions in response to inputs u0 and u1, it follows that ; Niv0u1i =
i=1 2k X 2k X Niv1u0i
i=1 Using the de nition of the matrices Ni, we can expand these equations to obtain 56 v11(u01 + u02) v12(u03 + u04) = (2.71) . . . . . . ;v0k(u1 2k;1 + u1 2k) v1k(u0 2k;1 + u0 2k ) Equation 2.71 is true if and only if the system is linear. If we can nd even a single contradiction to Eq. 2.71, then we will show that the system is nonlinear. Let u0 = u1 ) v0 = v1. Then from the above equation we have that ;v0 = v1, which is a contradiction the system of equations 2.67 is nonlinear. Equation 2.67 de nes a class of nonlinear systems known as bilinear systems. These systems have the property that the unforced response is described by a constant parameter linear state equation, however, the forced response adds a nonlinear term to the velocity eld. While the above proof demonstrates nonlinearity, it adds little physical intuition as to how this nonlinearity arises. The nonlinearity is an immediate consequence of a multiplicative interaction between the state and the input. This multiplicative nonlinearity arises from terms in the state equation which compute the synaptic current Ijsyn into each compartment: Ijsyn = vj (t) ; Esyn ]gjsyn (2.72) ;v01(u11 + u12) ;v02(u13 + u14) 2.4.1 A Simpli ed Case e i e i dv1 = ;v ( + g1 (t) + g1(t) )] ;1+v 1 ;1+ E e ; Er g1(t) ;1+ E i ; Er g1(t) ;1(2.75) 1 1 gm 2 L2 dt gm E e ; E i gm E e ; E i gm e i e i dvk = v 1 ;1;v ( + gk (t) + gk (t) )] ;1+ E e ; Er gk (t) ;1+ E i ; Er gk (t) ;1 (2.76) k;1 k 1 gm dt L2 gm E e ; E i gm E e ; E i gm where gm is the total nonsynaptic membrane conductance in each compartment, and 1 (2.77) 1 = 1 + L2 2 (2.78) 2 = 1 + L2 We have assumed that the excitatory and inhibitory synaptic current reversal potentials, E e and E i are the same in each compartment. De ning The state equations derived previously can be simpli ed considerably when modeling a passive, uniform cable partitioned into k segments of electrotonic length L: gje (t) gji (t) dvj = vj;1 1 2 ;1 ; vj ( 2 + gm + gm )] ;1 + vj+1 1 2 ;1 (2.73) dt L L E e ; Er gje(t) ;1 + E i ; Er gji (t) ;1 + E e ; E i gm (2.74) E e ; E i gm 57 e i e i e i 1 1 2 2 k k u(t) = gg(t) gg(t) gg(t) gg(t) : : : gg(t) gg(t) ]T m m m m m m e i e i G(t) = diag ; g1(t)g+ g1(t) : : : ; gk (t)g+ gk (t) ] m m E e;Er E i ;Er 0 0 0 : : i i E e;E1 E e ;E1 E e;Er E i ;Er 0 0 0 Ee;Ei Ee;Ei : : 0 0 0 0 : : : B= 0 0 0 0 0 0 : : : : : : e: : E ;Er E i ;Er 0 0 0 0 0 Ee;Ei Ee;Ei 0 0 0 0 0 0 0 0 ; 1 12 0 0 0 : : : 1 2 ;L 2 1 2 0 0 : : : 0 L 1 2 ;L 2 1 2 0 : : 0 : 0 L L A= 0 0 : : : : : : 0 : : : : : : : : : 12 ; 2 12 0 0 0 0 0 0 L L 0 0 0 0 0 0 0 12 ; 1 L we obtain the following state equations: v (t) = ;1 A + G(t)]v (t) + ;1Bu(t) _ (2.79) (2.80) E e;Er E e;E i : : : : : 0 E i ;Er E e ;E i 0 0 0 0 : 0 (2.81) 2.5 Implementation of Compartmental Models
2.5.1 Linear Models
When input to compartmental models is injected current, the resulting state equations are those of a linear timeinvariant system with velocity eld given by (Eq.2.17). Direct numerical integration of the velocity eld can be used to solve for the state vector as a function of time. However, equation 2.53 shows that it is possible to nd closed form solutions for the convolution integral in Eq.2.53 when current inputs are su ciently simple. In such cases, evaluation of Eq.2.53 will only require numerical computation of eigenvalues and eigenvectors of the state transition matrix A. This approach can lead to a considerable savings in computation time. Alternatives to the compartmental approach have also been used. Butz and Cowan (1974) developed a geometrical calculus for modeling passive dendritic trees with arbitrary branching patterns. Dendritic trees were assumed to be modeled as a concatenation 58 of uniform, passive cables, and impedance functions Kij (w) relating current injected at point j to voltage at location i were derived in the frequency domain. This approach was generalized by Koch and Poggio (1985) to include arbitrary linear membrane models. When the inputs to compartmental models are synaptic conductance changes, numerical integration of state equations must be employed. The righthand sides of equations 2.17 and 2.67 are vectors describing the time rate of change of membrane potential for the cases of current and synaptic inputs, respectively. Given an initial condition vector and these velocity elds, a variety of numerical techniques can be used to integrate the state equations. Mathematical software subroutine packages with implementations of RungeKutta, Adams, Gear, and BurlischStoer integration algorithms include the International Mathematics and Statistics Library (IMSL) and the Naval Surface Weapons Center (NSWC) libraries. A variety of alternative software tools are available for integrating the state equations. This includes general purpose mathematical software such as Mathematica, interactive software tools specialized for the study of nonlinear systems such as INSITE, KAOS, and Phaser, and special purpose tools designed speci cally for studying neuronal compartmental models (Genesis and Geppeto). Geppeto (Tiller et al., 1992) has the advantage that it is a software tool that automatically builds compartmental models of cells which have been reconstructed using a Eutectic Neuronal Reconstruction System. One of the rst approaches for simulating properties of compartmental models involved the use of SPICE, an electrical circuit simulation package (Segev et al., 1985). In order to use SPICE, it was necessary to calculate the values of all the resistive and capacitive circuit elements in the model, and then describe the interconnection of all these elements using the SPICE modeling language. Synaptic function inputs were built by designing an electrical circuit whose node voltage satis ed a di erential equation with solution equal to the  function. The voltage at this node was then used to control a voltage dependent conductance element available in the SPICE language to model a synaptic input. While the use of SPICE was an important early contribution to the eld of compartmental modeling, it has proven to be extremely slow. In addition, it is very di cult to model voltagedependent synaptic conductances in the SPICE language. One important numerical issue that must be considered when working with compartmental models has to do with the eigenvalues of the state transition matrix A. In general, the time rate of change of voltage in each compartment is determined by the eigenvalue i for that compartment. We previously derived upper and lower bounds on these eigenvalues in Eqs.2.332.35 for the case of a compartmental model subject to current inputs. A similar approach can be taken for the model in which inputs are membrane conductance changes. Assume that all speci c membrane parameters are uniform throughout the dendritic tree. Let = Rm Cm be the membrane time constant. Then it is possible to show, using the 59 2.5.2 Nonlinear Models approach developed in Sect.2.4.1, that the system eigenvalues are constrained to lie in the intervals ge (t) + gi (t) ge (t) + gi (t) g (2.82) dj = ;1 ;(1 + j gm j ) ; 2( gj +m j+1 ) ;(1 + j gm j )] gj j j e i e i g2 (2.83) d1 = ;1 ;(1 + g1 (t)g+ g1(t) ) ; 2 gm ) ;(1 + g1(t)g+ g1(t) )] m m 1 1 1 e i e i gk dk = ;1 ;(1 + gk (t)g+ gk (t) ) ; 2( gm ) ;(1 + gk (t)g+ gk (t) )] (2.84) m m This simpli es to
k k k (t) 2 (1 + Lj 2 Aj+1 ) ;(1 + gje(t) + gji (2.85) )] 2 2 Aj m gj Lj Lj+1 e i 2 L12 A2 ;(1 + g1 (t) + g1(t) )] d1 = (2.86) m g1 L12 L22 A1 e i 2 ;(1 + gk (t) + gk (t) )] dk = (2.87) m gk Lk 2 where Lj and Aj are the electrotonic length and surface area of the j th compartment, respectively. Therefore, in general, model time constants are determined by both the relative electrotonic lengths and surface areas of compartments, as well as the magnitude of the conductance changes produced by synaptic inputs. Di culties arise when the eigenvalues begin to take on radically di erent values. This can happen in several ways. First, because of the way in which the dendritic tree is partitioned, electrotonic lengths and surface areas of compartments could di er. A second factor plays a role when large conductance changes occur in some compartments, but not others. This will make the lower bound on certain eigenvalues increasingly negative. The characteristic time constants associated with these eigenvalues (the inverse of the eigenvalue) will become smaller (decay faster), thereby forcing the use of very short integration time steps. Systems with this property are referred to as "'sti "' systems. Numerical integration of such systems is time consuming because the largest integration time step that may be employed is determined by the most rapidly decaying characteristic time constant. In such cases, it is necessary to employ a sti integration routine such as Gear's algorithm. This problem is not present when simulating a cable of uniform dimension partitioned into compartments of equal electrotonic length. In this case, it is possible to derive explicit solutions for the eigenvalues of A. These eigenvalues are all clustered near 1. dj = e i ;1 ;(1 + gj (t) + gj (t) ) ; gjm e i ;1 ;(1 + g1 (t) + g1 (t) ) ; m g1 e i ;1 ;(1 + gk (t) + gk (t) ) ; m gk 60 3 Design Principles
3.1 Synaptic Saturation
Referring back to Fig. 2.1, we can derive the relationship between synaptic conductance and membrane potential for the case of an isopotential compartment. For simplicity, lets ignore the inhibitory synaptic input. Application of Kircho 's current law yields a di erential equation for v(t). De ning a normalized voltage v (t) = v(t) ; Er , letting _ rm = g1 , and setting v (t) to zero yields an expression for the DC potential in response to m a maintained synaptic conductance input ge : 3.1.1 Single Synaptic Inputs germ = (E ; E ) ggm v = (Ee ; Er ) 1 + ge rm (3.1) e e r 1 + ggm A plot of the steadystate normalized membrane potential v as a function of the base 10 logarithm of synaptic conductance ge relative to the passive membrane conductance gm is shown in Fig. 3.1. We can determine a region where there is an approximate linear e mapping from normalized conductance ggm to membrane potential v by expanding Eq. 3.1 in a Taylor series about 0: ge ge ge ge (E e ; E r ) ( gm ) ; ( gm )2 + ( gm )3 + 0( gm )4] (3.2) The rst term of this function is shown by the dashed curve in Fig.3.1. The relationship between synaptic conductance and voltage is approximately linear as long as the ratio of synaptic conductance ge to membrane conductance gm is small (less than about .1). Within e this regime, cinductance inputs can be modeled as current inputs with Iin = (Ee ; Er ) ggm .
61
ge gm (Ee ; Er ) ge 1 + gm e v* 1 0.8 0.6 0.4 0.2 3 2 1 1 2 3 Log10(Ge/Gm) Figuree 3.1: Plot of steadystate normalized membrane potential v (ordinate) versus log10 ggm (abscissa). Solid line is the function given by Eq. Dashed line is the linear approximation. 62 4 References
Butz, E. G. and Cowan, J. D. (1974). Transient potentials in dendritic systems of arbitrary geometry. Biophys. J., 14:661{689. Coleman, P. A. and Miller, R. F. (1989). Meausrement of passive membrane parameters with wholecell recordings from neurons in the intact amphibian retina. J. Neurophysiol., 61:218{230. Jack, J. J. B., Noble, D., and Tsien, R. W. (1975). Electric current ow in excitable cells. Oxford Science Publications, Clarendon Press, Oxford, UK. Koch, C. and Poggio, T. (1985). A simple algorithm for solving the cable equation in dendritic trees of arbitrary geometry. J. Neurosci. Meth., 12:303{315. Lux, H. D., Schubert, P., and Kreutzberg, G. W. (1970). Direct matching of morphological and electrophysiological data in cat spinal motoneurons. In andJ. K. S. Jansen, P. A., editor, Excitatory Synaptic Mechanisms, pages 189{198, Oslo:Universitetsforlaget. Miller, R. F. and Bloom eld, S. A. (1983). Electroanatomy of a unique amacrine cell in the rabbit retina. Proc. Natl. Acad. Sci. USA, 80:3069{3073. Perkel, D. and Mulloney, B. (1978a). Calibrating compartmental models of neurons. Am. J. Physiol., 235:R93{R98. Perkel, D. and Mulloney, B. (1978b). Electrotonic properties of neurons: Steadystate compartmental models. J. Neurophysiol., 41:627{639. Pickard, W. F. (1971). The spatial variation of plasmalemma potential in a spherical cell polarized by a small current source. Math. Biosci., 10:307{328. Quinney, D. (1985). An Introduction to the Numerical Solution of Di erential Equations. Research Studies Press, Letchworth, UK. 63 Rall, W. (1959). Branching dendritic trees and motoneuron membrane resistivity. Exp. Neurol., 1:491{527. Rall, W. (1960). Membrane potential transients and membrane time constant of motoneurons. Exp. Neurol., 2:503{532. Rall, W. (1962). Theory of physiological properties of dendrites. Ann. N.Y. Acad. Sci., 96:1071{1092. Rall, W. (1964). Theoretical signi cance of dendritic trees for neuronal inputoutput relations. In Reiss, R. F., editor, Neural Theory and Modeling, pages 73{97, Stanford, CA. Stanford University Press. Rall, W. (1969). Time constants and electrotonic length of membrane cylinders and neurons. Biophys. J, 9(12):1483{1508. Rall, W. (1973). Branch input resistance and steady attenuation for input to one branch of a dendritic neuron model. Biophys. J, 13:648{688. Rall, W. (1977). Core conductor theory and cable properties of neurons. In Handbook of Physiology: The Nervous System. Cellular Biology of Neurons, volume 1, chapter 3, pages 39{97. Am. Physiol. Soc., Bethesda, MD. Redman, S. and Walmsley, B. (1983). The time course of synaptic potentials evoked in cat spinal motoneurones at identi ed group ia synapses. J. Physiol., 343:117{133. Rinzel, J. and Rall, W. (1974). Transient response in a dendritic neuron model for current injected at one branch. Biophys. J., 14:759{790. Segev, I., Fleshman, J. W., Miller, J. P., and Bunow, B. (1985). Modeling the electrical behaviour of anatomically complex neurons with a network analysis program: Passive membrane. Biol. Cybern., 53:27{40. Shelton, D. P. (1985). Membrane resistivity estimated for the purkinje neuron by means of a passive computer model. Neurosci., 14:111{131. Tiller, M., Jayaram, S., Winslow, R. L., and Miller, R. F. (1992). Geppeto: A graphics application for generating and manipulating compartmental models. J. Neurosci. Meth., page submitted. Weinberger, H. F. (1965). A First Course in Partial Di erential Equations with Complex Variables and Transform Methods. Blaisdell Publishing Company, New York. 64 ...
View
Full
Document
This note was uploaded on 04/15/2008 for the course BME 580.406 taught by Professor Wang during the Spring '08 term at Johns Hopkins.
 Spring '08
 Wang
 Biomedical Engineering

Click to edit the document details