辅导voltage、讲解CS/python, C/C++编程语言、辅导PDE留学生

2018.11.25 - 首页 >> Python编程

Propagation of voltage in a neuron: the cable equation

Connection with course material: During the semester, we have spent a substantial amount of time studying

the heat equation:

u(x, t)

t = k

2u(x, t)

. (1)

Indeed, Eq. (1) is a foundational model in mathematical physics. Not only is Eq. (1) relevant for the

transfer of heat via thermal conduction, but it is also mathematically describes the process of diffusion

whereby molecules move from a region of high to low concentration. Furthermore, Eq. (1) is at the heart of

classical cable theory, which describes the propagation of electric current along the membrane of neurons

(the cellular building blocks of the brain). However, one must incorporate additional terms, described by the

function f(u), into the partial differential equation (1) to account for ion channels, which open and close to

let ions in and out of the cell in response to changes in a neuron’s voltage:

v(x, t)

t =2

v(x, t)

+ f(v(x, t)), (2)

where we have set k = 1. In this project, you will incorporate two different types of channel models into

Eq. (2), analyze the resulting models, and study the predictions of the associated solutions.

Background: A.L. Hodgkin and A.F. Huxley received the 1963 Nobel Prize in Physiology/Medicine for

their groundbreaking work describing the initiation and propagation of action potentials in the squid giant

axon [1]. Central to their work was the observation that each part of the axonal membrane could be treated

as a component in an electric circuit (See Fig. 1). An ordinary differential equation model describing the

dynamics of the given electric circuit produces action potentials given the right parameter values. Action

potentials, also called “spikes,” are abrupt events in which the voltage of a neuron rapidly rises and falls.

Such activations are the primary units of information in the brain, since they allow communication in between

cells via synapses. While the genesis of spikes is interesting in its own right, this project will be

concerned with how these spikes can propagate from a neuron’s main body (soma) down an axon.

What follows is a derivation of a partial differential equation (PDE) model describing the propagation

of potential along the length of a neuronal membrane. To start, you will consider a passive membrane so

that ions simple “leak” out of the cell if there is a voltage gradient. To gain an intuitive understanding of the

resulting equation’s behavior, you will find the stationary solutions of the model to start. After this, you will

solve for the Green’s function of the PDE and use it to determine the impulse response of the membrane

to inputs. Then, a nonlinear model will be introduced, which accounts for the “two state” nature of some

neural ion channels. Rather than analyzing all possible solutions, you will solve the resulting equation for

traveling wave solutions, given appropriate boundary conditions.

Deriving the model: The current per unit length Jm passing through the membrane of a neuron is given in

terms of voltage v by:

Jm(x) = x 1

R(x)

vx,

also known as the heat or diffusion equation. Spatial transfer of potential v throughout the neural “cable”

is determined by the inverse resistance: R(x) = Ri(x) + Re(x), the sum of the intra- and extracelluar

resistences (See Fig. 2). On an infinitesimal segment [x, x + dx] with intra- and extracelular resistences

Ridx and Redx, denote the intra- and extracellular voltages at either end as vi(x), ve(x), vi(x + dx), and

Nernst-Planck, and Nernst equations. In particular, Section 2.2.4 has a good explanation of spacecharge

neutrality. The GHK equations are derived, and there are examples of how to use the

concepts, including a worked-out resting potential exercise on p. 30. Chapter 3 has equivalent

circuit descriptions.

3.3 The Hodgkin-Huxley (H-H) equations

The Hodgkin-Huxley (H-H) model of the squid giant axon [126] is a triumph of mathematical

modeling. Their formulation of an ODE model that reproduces the action potential, and extension

to a partial differential equation (PDE) to explain its propagation along the axon, marks the

beginning of mathematical neuroscience. The model was built on more than 15 years of beautiful

and painstaking experimental work, interrupted by World War II, that culminated in a remarkable

series of papers [127, 128, 124, 123, 125, 126]. This gained Hodgkin and Huxley a Nobel prize in

1963, along with J.C. Eccles (for his work on synapses and discovery of excitatory and inhibitory

post-synaptic potentials: see §4.1). Here we give only the essence of the final mathematical model

from [126]; for an introduction with some historical notes, see [150, §5.1]. See [122] for a personal

story that also has general relevance for scientific research.

Exercise 15. A possible journal club team project: Read and present the experimental H-H(-K)

papers [127, 128, 124, 123, 125].

Figure 3.5: The equivalent circuit for the giant axon of squid, reproduced from [126, Fig. 1].

The leak conductance gL is assumed constant, but the sodium and potassium conductances vary,

indicated by the arrow crossing the resistor symbol. Batteries in series with each resistor represent

reversal potentials. The external current I flows inward in the convention of [126], in contrast to

our convention of outward flow (Fig. 3.4).

The H-H model comprises four coupled nonlinear first-order differential equations, listed below,

Figure 1: The equivalent circuit for the giant axon of squid [1]. The leak conductance gL is assumed

constant, due to the fact that chloride ions constantly flow through the channel, depending only on the

gradient across the cell membrane. The sodium and potassium conductances gNa and gK vary with changes

in the voltage v, indicated by the arrow crossing the resistor symbol. Batteries in series with each resistor

represent reversal potentials. The external current I flows inward in the convention of [1].

ve(x + dx) (See Fig. 2a). Ohm’s law states that the internal and external currents Ii(x) and Ie(x) are related

to the voltages as

vi(x + dx vi(x) = Ii(x)Ridx, ve(x + dxve(x) = e(x)Redx.

Dividing by dx and taking the limit as dx → 0 yields

Ii(x) =

Ri(x)

x and Ie(x) =1

Re(x)

vex . (3)

As demonstrated in the current diagram (Fig. 2b), Kirchoff’s laws relate the current across the membrane

dx, having intra- and extracellular resistences Ridx, Redx, and denote the intra- and extracellular

voltages at either end as vi(x), ve(x) and vi(x + dx), ve(x + dx): Fig. 3.15(a). With the convention

of positive current flow from left to right (x to x + dx), Ohm’s law implies that

vi(x + dx) vi(x) = Ii(x)Ridx and ve(x + dxve(x) = Ie(x)Redx, (3.44)

where Ii, Ie denote the internal and external currents in the axial direction. Dividing by dx and

taking te limit dx → 0, as in calculus, gives

Ii(x) = 1

and Ie(x) = 1

x . (3.45)

Figure 3.15: a) Current flow along an axon. b) Schematic for the Kirchhoff law current balance.

We next balance the axial currents Ii, Ie at x and x + dx with the transmembrane current,

appealing to Kirchhoff’s laws. If Im is the transmembrane current per unit cable length, this gives

Ie(x + dx)Ie(x) = Imdx = Ii(x)Ii(x + dx) : (3.46)

in words: what is lost in external current is gained in internal current. See Fig. 3.15(b) and recall

that outward transmembrane currents are positive: §3.2. Dividing by dx and taking the limit again

gives

Im(x) = Ie

x = Ii

x . (3.47)

We now use the facts that v = vi ? ve and the total axial current Ia = Ii + Ie is constant to

eliminate explicit reference to the internal and external currents and voltages. From (3.45) and

ve = vi v we get

Ia = Ii + Ie =1

x 1

(vi v)

x =

!Ri + Re

RiRe

(3.48)

Since Ia/x =0(Ia constant), this gives

# 1

(3.49)

Figure 2: (a) Current flowing along a neuronal axon. (b) Schematic for Kirchoff law current balance.

per unit length Jm(x) to the Ie(x) and Ii(x) as:

Ie(x + dx) Ie(x) = Jm(x)dx = Ii(x) Ii(x + dx).

Divide these by dx and taking the limit dx → 0 yields:

Jm(x) = Ie . (4)

Using the fact that v = vi?ve is the voltage across the cell membrane and the total axial current Ia = Ii+Ie

is constant, we can eliminate explicit reference to the internal and external currents and voltages:

Ia = Ii + Ie =

ReIa

Ri + Re

x +

Ri + Re

and since Ia is constant, this gives

Ri(x) + Re(x)

Re(x)

Ri(x) + Re(x)

Finally, assuming constant resistances (Re(x) ≡ Rˉ

e, Ri(x) ≡ Rˉ

i), and applying Eq. (3) and (4), we derive

the expression

Jm(x) =

The neural membrane as a whole behaves like a resistor-capacitor (RC) circuit, and the nonlinear properties

of the resistor-like ion channels (See Fig. 1) yield a nonlinear differential equation for the evolution of

the transmembrane voltage v in time t. Considering the time variable as well now, the current across the

membrane per unit length Jm(x, t) is given by the sum of the capacitive current JC(x, t) = C

, the outward

ionic JION (v(x, t), t), and inward applied current JA(x, t) per unit length [1]:

Jm(x, t) = JC(x, t) + JION (v(x, t), t) JA(x, t) (5)

so upon plugging in Eq. (5), and applying the expression for capacitive current yields:

v(x, t)

t = JION (v(x, t), t) + 1

v(x, t)

+ JA(x, t). (6)

Changing variables x 7→ y/√

R and t 7→ Cs in Eq. (6) cancels the capacitance C and resistance R parameters.

Furthermore, we assume the ionic currents are given by a nonlinear function f(v(y, s)), which is only

a function of the rescaled voltage and accounts for the multiple possible stable states of the voltage [2]. The

justification of the form for f(v(y, s)) will be discussed in more detail with each example. This yields the

cable equation:

v(y, s)

s =2

v(y, s)

+ f(v(y, s)) + Jext(x, t).

For familiarity’s sake, we write this in terms of the original variables (x, t):

v(x, t)

t =2v(x, t)

+ f(v(x, t)) + Jext(x, t). (7)

1. Passive membrane. The ionic currents described by the function f(v(x, t)) in Eq. (7) can in principle

comprise many different types of ion channel. For simplicity, we begin by considering a passive

membrane current. This means that the current flows according to Ohm’s law VION = IION RION ,

which we can write IION = VION /RION . For Eq. (7), assuming resistances can be rescaled to unity,

this means f(v) = v. Therefore, we are working with the inhomogeneous PDE:

v(x, t)

t =2

v(x, t)

v(x, t) + Jext(x, t). (8)

Voltage changes at location x are determined by: (i) diffusion of voltage from high to low voltage

regions (vxx(x, t)); (ii) decay of voltage due to ions leaking out of the cell (v(x, t)); and (iii) a

source term due to externally applied current (Jext(x, t)).

(a) Stationary solutions. To understand the behavior of the cable equation for a passive membrane,

Eq. (8), let’s begin by looking at stationary solutions where vt(x, t) = 0. This will help us understand

the spatial effects of the external input Jext(x, t) before tackling the full dynamics (both spatial and

temporal). Start by solving for the homogeneous solutions (where Jext(x, t) ≡ 0). The stationary

solutions will only be functions of x: v(x, t) := V0(x). Consider an infinite-length cable so x ∈

(∞,∞), a modeling approximation which makes use of the fact that neural axons and dendrites

tend to be much, much longer than they are wide.

(b) Impulse response: stationary solution. Now, for a stationary input current localized at x = 0:

Jext(x, t) = δ(x), and boundary conditions limx→±∞ v(x, t) = 0, solve for the stationary solution:

v(x, t) = VF (x) of Eq. (8). Using Fourier transforms may be the most straightforward way, but you

can also use a fundamental solution or Green’s function approach.

vanishing boundary conditions. Utilizing the principle of superposition for translated fundamental

solutions, derive the formula for a stationary solution v(x, t) = VA(x) given an arbitrary stationary

input Jext(x). Hint: It should be an integral. This is often called a Green’s function approach to

fundamental solutions.

(d) Green’s function. Now, return to the full time-dependent Eq. (8). Solve for the Green’s function

or fundamental solution. That is, given a δ-impulse at t = 0 and x = 0, what is the solution v(x, t)

Solve the following equation using Fourier transform methods:

v(x, t)

t =2

v(x, t)

v(x, t) + δ(t)δ(x),

and call the solution G∞(x, t), the Green’s function on the infinite cable. Create some plots of this

function using mathematical software, at t = 0, t = 1, t = 10, t = 100. Demonstrate that the voltage

is monotone decreasing at x = 0 on t ∈ (0, ∞), whereas it is nonmonotonic at x = 2 on t ∈ (0,∞).

Finally, compute the total voltage Gˉ∞(t) = R ∞

∞ G∞(x, t)dx and show it is monotone decreasing on

t ∈ (0,∞).

(e) Generalize solutions. Now, utilize the principle of superposition for translated fundamental solutions

to derive the formula for a generalized solution given arbitrary spatiotemporal input Jext(x, t).

Hint: It should be an integral over both x and t. Create some plots of this function for Jext(x, t) =

[δ(x 1) + δ(x + 1)] δ(t) as well as Jext(x, t) = 1 as well as Jext(x, t) = sin(x). You will have to

compute the integrals for each input function.

(f) Extra credit: Comparison to numerical simulations. Use MATLAB (or your favorite mathematical

computing software) to numerically simulate Eq. (8). You will want to use a finite different

approximation to simulate this PDE. There are many great resources online that explain how to do

this. Here are a few:

http://www.ehu.eus/aitor/irakas/fin/apuntes/pde.pdf

https://www.mathworks.com/moler/pdes.pdf

http://www.colorado.edu/geography/class_homepages/geog_4023_s07/labs/

html/PDE_lab.html

http://www4.ncsu.edu/?zhilin/TEACHING/MA587/fd_fem.pdf You do not actually

need to specify boundary conditions when you simulate the PDE. Start off by using a “near-delta”

function, such as

Jext(x, t) = 10e?25x

δ(t). (9)

Rather than simulating a temporal delta function δ(t), just make Eq. (9) the initial condition. Compare

this with the results of your analysis in part (e) by computing the supposed solution by hand.

Note, we did specify the spatial (and time) domain to be infinite, but you will clearly have to use a

finite spatial domain and finite time to stop. Take x = ?10 and x = 10 to start, and use t = 50 as

your stopping time. When you discretize your PDE using finite differences on the spatial operator

and forward Euler on the timestep, take ?x = 0.1 and ?t = 0.1. Does changing these values make a

big difference in the result?

2. Traveling wave solutions to the PDE (7). When modeling “biological” neurons, one typically incorporates

some of the active properties of cells into the cable equation (7). In particular, we want to

account for some of the nonlinear effects that arise in the dynamics of ion currents that control and

are controlled by the local voltage of the neural membrane. While there are many different ways of

modeling this via the term f(v(x, t)) in Eq. (7), one typical example is to choose a function that makes

the membrane bistable. A bistable system has two possible stable equilibria, a nice simple model of

the activation of a neural membrane. Without any input, the membrane exhibits a low baseline potential,

but given enough input, the membrane becomes “active” and ion channels along the membrane

reenforce the potential excursion of the current input. The two stable states we will consider are the

inactive state at v = 0 and the active state at v = 1. It turns out that when a function f(v) possessing

these equilibria is incorporated into Eq. (7), one can identify traveling wave solutions. You will now

derive traveling wave solutions to Eq. (7) by projecting the PDE down to a system of ODEs that you

can solve analytically. Why do this? Why study this particular type of solution? One reason is that

these are the types of solutions that are often observed when recording from a single neuron. An

impulse current to a neural membrane does not always just spread out, as in the case of a passive

membrane, it tends to propagate like a wave approaching a beach. You will also have the opportunity

to compare your analytical results with numerical simulations. To facilitate this, we will assume a

simple form for the nonlinearity f(v(x, t)). Consider the Heaviside step nonlinearity:

f(v) = ?v + H(v ? θ) (10)

in Eq. (7), where H(v ? θ) is the Heaviside step function with threshold θ.

(a) To identify traveling wave solutions, first look at homogeneous solutions, which are constant in

space and time: v(x, t) = ˉv. What are the two possible homogeneous solutions (call them vˉ

+ and vˉ)

of Eq. (7) given (10)? What condition on θ must be satisfied for both to exist?

(b) A reasonable simplification of neuronal axons is that, because they are long in comparison to their

diameter, we can model them as infinitely long cables so x ∈ (?∞,∞). Physically realistic solutions

must therefore obey the boundary conditions

limx→∞

v(x, t)

x = 0 and lim x→∞

v(x, t)

x = 0. (11)

We are specifically looking for traveling wave solutions v(x, t) = V (ξ), where ξ = x ct is the

traveling wave coordinate. Apply this change of variables to the derivatives in Eq. (7) as well as the

boundary condition Eq. (11) to convert to a second order ordinary differential equation.

ODE, with a couple caveats. First, what is the value of H(V (ξ) ? θ)? It can either be 0 or 1 at any

point ξ. We will assume a traveling front solution that decreases as ξ increases, and crosses through

zero at ξ = 0, so that

V (ξ) > θ when ξ < 0 (12)

V (ξ) < θ when ξ > 0. (13)

This means that you have to solve two different 2nd order ODEs: one on ξ ∈ (?∞, 0) and one on

ξ ∈ (0,∞). What will the boundary conditions of the each solution be? One is given by the ξ-version

of Eq. (11), another by the requirement that V (ξ) tends to a homogeneous solution at ξ → ±∞, and

another requires that the solutions match at ξ = 0. It turns out that not only should V (ξ) match at

ξ = 0, but also V

(ξ). Why do you think this is? List all the boundary conditions and the two ODEs

you have derived.

(d) Solve the two ODEs you derived in part (iii) using the method of undetermined coefficient (or another

method if you choose). All constants should be specified by applying the boundary conditions.

(e) Finally, specify the speed c of the traveling front by applying the threshold condition V (0) = θ

and solving for c. The threshold θ represents how “active” the axon is (lower θ means more active).

Does the relationship between c and θ make sense in light of this?

(f) Traveling front solutions constructed in (a)-(e) are asymptotically stable, which means that if initial

conditions start out near a traveling front solution they will converge to such a solution. How can

we prove this? One way to do so is using a linear stability analysis. What does this mean? Basically,

we add a perturbation to the traveling wave solution, derive an equation for its evolution, and check

if it decays in time. If it does decay in time, that means solutions in the vicinity of the wave will tend

to be attracted to it.

To start, we will plug the solution v(x, t) = V (ξ) + ψ(ξ, t) (0 < 1) into Eq. (7) and linearize.

That is, Taylor expand f(V + ψ) using the fact that is small, and eliminate all terms that are not

on the order of . Note, you will need to apply the chain rule to differentiate ψ(ξ, t), since it depends

on ξ = x ? ct. This will be a linear PDE, which can be solved using separation of variables. One

solution corresponds to ψ(ξ, t) = U

(ξ). What do you think this solution corresponds to? There will

also be other solutions of the form ψ(ξ, t) = eλtψˉ(ξ) such that λ < 0, but there should be no solutions

with λ > 0.

(g) Extra credit: Comparison to numerical simulations.Use MATLAB (or your favorite mathematical

computing software) to numerically simulate Eq. (??) given (10). Also, compute the speed of the

traveling wave numerically for a few different values of θ. Figure out a good way to do this. How do

you estimate the speed of a moving object? You will want to use a finite different approximation to

simulate this PDE. There are many great resources online that explain how to do this. Here are a few:

http://www.ehu.eus/aitor/irakas/fin/apuntes/pde.pdf

https://www.mathworks.com/moler/pdes.pdf

http://www.colorado.edu/geography/class_homepages/geog_4023_s07/labs/

html/PDE_lab.html

http://www4.ncsu.edu/?zhilin/TEACHING/MA587/fd_fem.pdf

You do not actually need to specify boundary conditions when you simulate the PDE. Produce plots

of traveling waves using three different threshold values θ. Note, we did specify the spatial (and time)

domain to be infinite, but you will clearly have to use a finite spatial domain and finite time to stop.

Take x = 10 and x = 10 to start, and use t = 50 as your stopping time. When you discretize your

PDE using finite differences on the spatial operator and forward Euler on the timestep, take x = 0.1

and t = 0.1. Does changing these values make a big difference in the result?

(vii) Extra credit. Numerically simulate the system with

f(v, x) = v + H(v θ(1 + 0.5 cos(x))), (14)

a periodically varying threshold. How does this alter the behavior of the traveling wave solution? As

you increase the constant in front of cos(x) above, what happens to the average speed of the front.

Does the front ever stop completely? Such a modification models the spatial heterogeneity in the

excitability of the axon.

References

[1] A. L. HODGKIN AND A. F. HUXLEY, A quantitative description of membrane current and its application

to conduction and excitation in nerve, The Journal of physiology, 117 (1952), p. 500.

[2] J. P. KEENER AND J. SNEYD, Mathematical physiology, vol. 1, Springer, 1998.