MATH564 Homework 5: Neurophysiology

Topics

Chemical intra-organism communication

Diffusion

Estimate whether diffusion is a viable intra-organism communication process. First, obtain diffusivity data

imTab622 = ;
Tab622 = Map[N[ToExpression[#]]&,TextRecognize[imTab622,"Word"]]

Length[Tab622]

DifData = ArrayReshape[Drop[Tab622,{3}],{6,10}]

αData = (10^(-6))DifData[[4]]

Estimate the time required to transmit a chemical signal through diffusion over a distance of (10^(-3)) cm

Map[(10^(-6))/(2 #) &, αData]

Estimate the time required to transmit a chemical signal through diffusion over a distance of (10^(-1)) cm

Map[(10^(-2))/(2 #) &, αData]

Diffusion “velocities” over one second

Map[Sqrt[2 #]&, αData]

Convection

Typical blood flow velocities in various circulatory vessels

Time required to transmit a signal over 1 cm varies from 0.1 to 0.01 seconds, much faster than diffusion.

Electrical intra-organism communication

Ion properties

Define ions active in neurotransmitter channel

ions={"(SuperPlus[Na])","(SuperPlus[K])","(SuperMinus[Cl])"}

Define intra- and extra-cellular concentrations

(Subscript[C, i])={15.,150., 7.5}
(Subscript[C, o])={145.,5.,110.}

Nernst equivalent voltages

R=8.314 ;(* perfect gas constant J/(mol K) *)
F=96485; (* Faraday's constant C *)
V[Ci_,Co_,T_]=R T/F Log[Ci/Co]

TableForm[MapThread[ 1000V[#1,#2,282.]&,{(Subscript[C, i]),(Subscript[C, o])}]]

Equilibrium potentials

EqPot={55.,-82.,-59.}

Hodgkin-Huxley model

General considerations

Salts are molecules with typical ionic bonds (NaCl, KCl). When in an aquaeous solution some salt molecules dissociate forming ions

m (v^2)/2= n k T /2=3 k T/2

v=(Sqrt[3 1.38 (10^(-23)) 310/18/1.66/(10^(-27))])

Ohm’s law states that V = I R

Here are schematics and circuit diagrams modeling current flow through the axon membrane

;

Observations on the squid axon shows a sigmoid (looks like an S which the Greek Sigma) polarization/depolarization voltage versus time dependence. Here is the prototypical sigmoid-shape function

sigmoid[x_] = (1/(1+Exp[-x]))

Plot[sigmoid[x],{x,-5,5}]

A typical decay is provided by a solution to a linear first-order differential equation

n[t_]= (Subscript[n, ∞])-((Subscript[n, ∞])-(Subscript[n, 0])) Exp[-t/(Subscript[τ, n])]

Simplify[n'[t]]

params={ (Subscript[n, ∞])->0,(Subscript[n, 0])->1, (Subscript[τ, n])->1}

Plot[Evaluate[n[t] /. params],{t,0,5}]

Simplify[DSolve[{y'[t] == (yinf-y[t])/(Subscript[τ, n]) ,y[0]==y0},y[t],t]]

Clear[n];
n[t_,α_,β_]= n[t] /. DSolve[{n'[t] == α (1-n[t]) - β n[t],n[0]==0},n[t],t][[1,1]]

Manipulate[Plot[n[t,α,β],{t,0,1},PlotRange->All,Frame->True,FrameLabel->{"t","n(t)"}],{α,-3,3},{β,-3,3}]

Potassium conductance model

Using the above observations, the Hodgkin-Huxley model posits the voltage to vary as some power of the exponential decay solution. Empirically for potassium a fourth-power fits the data.

n[t_,nInf_,n0_,τn_]=nInf-(nInf-n0) Exp[-t/τn]

Series[((n[t,nInf,0,τ])^4),{t,0,5}]

n4[ t_,nInf_,τ_]=Normal[(((Series[n[t,nInf,0,τ],{t,0,2}]))^4)]

Potassium conductance constant

gKbar=24.34 (* m mhos/(cm^2) *)

αn[V_,Vr_]=((0.01 (10-(V-Vr)))/(Exp[1-0.1(V-Vr)]-1))

βn[V_,Vr_]=0.125 Exp[-((V-Vr)/80)]

gK [t,tinf,τ]= gKbar n4[t,tinf,τ]

Sodium conductance model

Sodium conductance model constant

gNabar = 70.7 (* m mhos/(cm^2) *)

gNa [t,tinf,τ]= gNabar (((m[t,tinf,τ]))^3)

Sodium, Potassium slope constants

Vr=-70. (* mV *)

αn[V_]=((0.01 (10-(V-Vr)))/(Exp[1-0.1(V-Vr)]-1));
βn[V_]=0.125 Exp[-((V-Vr)/80)];
αm[V_]=((0.01 (25-(V-Vr)))/(Exp[2.5-0.1(V-Vr)]-1));
βm[V_]=4 Exp[-((V-Vr)/18)];
αh[V_]=0.07 Exp[-0.05 (V-Vr)];βh[V_]=(1/(Exp[3-0.1(V-Vr)]+1));

Model equations

Additional model constants

ENa=55.; (* mV *)
EK =-82.; (* mV *)
El=-59.; (* mV *)
Cm = 0.001; (* F/(cm^2) *)
gl = 0.3; (* m mhos/(cm^2) *)

pulse[t_,Vp_]=Vp HeavisidePi[((t-0.001)/0.001)];

Plot[pulse[t,1],{t,-.005,.005}]

Differential equations Vp=20

eqs = { m'[t] == αm[V[t]] (1-m[t]) - βm[V[t]] m[t],
          n'[t] == αn[V[t]] (1-n[t]) - βn[V[t]] n[t],
        h'[t] == αh[V[t]] (1-h[t]) - βh[V[t]] h[t],
       V'[t] == -(( 1)/Cm)( gNabar (((m[t]))^3) h[t] (V[t] - ENa) +
gKbar (((n[t]))^4) (V[t]-EK) + gl (V[t]-El)+pulse[t,20])}

TableForm[eqs]

Initial conditions

conds = {V[0]== Vr, n[0]==0.315, m[0]==0.042, h[0]==0.608}

HH[t_] ={m[t]/0.042,n[t]/0.315,h[t]/0.608,V[t]/Vr} /.NDSolve[Flatten[{eqs,conds}],{m[t],n[t],h[t],V[t]},{t,0,10}][[1]]

pltVp20=GraphicsGrid[{{Plot[HH[t][[1;;3]],{t,0,10},Frame->True,FrameLabel->{"t","Na,K,Cl"},GridLines->Automatic,PlotLabel->"Ion concentrations in axon in response to synapse"]},{Plot[HH[t][[4]],{t,0,.1},Frame->True,FrameLabel->{"t","V"},GridLines->Automatic,PlotRange->All,PlotLabel->"Voltage in axon in response to synapse"]}}]

Differential equations Vp=10

eqs = { m'[t] == αm[V[t]] (1-m[t]) - βm[V[t]] m[t],
          n'[t] == αn[V[t]] (1-n[t]) - βn[V[t]] n[t],
        h'[t] == αh[V[t]] (1-h[t]) - βh[V[t]] h[t],
       V'[t] == -(( 1)/Cm)( gNabar (((m[t]))^3) h[t] (V[t] - ENa) +
gKbar (((n[t]))^4) (V[t]-EK) + gl (V[t]-El)+pulse[t,10])}

TableForm[eqs]

Initial conditions

conds = {V[0]== Vr, n[0]==0.315, m[0]==0.042, h[0]==0.608}

HH[t_] ={m[t]/0.042,n[t]/0.315,h[t]/0.608,V[t]/Vr} /.NDSolve[Flatten[{eqs,conds}],{m[t],n[t],h[t],V[t]},{t,0,10}][[1]]

pltVp10=GraphicsGrid[{{Plot[HH[t][[1;;3]],{t,0,10},Frame->True,FrameLabel->{"t","Na,K,Cl"},GridLines->Automatic,PlotLabel->"Ion concentrations in axon in response to synapse"]},{Plot[HH[t][[4]],{t,0,.1},Frame->True,FrameLabel->{"t","V"},GridLines->Automatic,PlotRange->All,PlotLabel->"Voltage in axon in response to synapse"]}}]

Differential equations Vp=5

eqs = { m'[t] == αm[V[t]] (1-m[t]) - βm[V[t]] m[t],
          n'[t] == αn[V[t]] (1-n[t]) - βn[V[t]] n[t],
        h'[t] == αh[V[t]] (1-h[t]) - βh[V[t]] h[t],
       V'[t] == -(( 1)/Cm)( gNabar (((m[t]))^3) h[t] (V[t] - ENa) +
gKbar (((n[t]))^4) (V[t]-EK) + gl (V[t]-El)+pulse[t,5])}

TableForm[eqs]

Initial conditions

conds = {V[0]== Vr, n[0]==0.315, m[0]==0.042, h[0]==0.608}

HH[t_] ={m[t]/0.042,n[t]/0.315,h[t]/0.608,V[t]/Vr} /.NDSolve[Flatten[{eqs,conds}],{m[t],n[t],h[t],V[t]},{t,0,10}][[1]]

pltVp05=GraphicsGrid[{{Plot[HH[t][[1;;3]],{t,0,10},Frame->True,FrameLabel->{"t","Na,K,Cl"},GridLines->Automatic,PlotLabel->"Ion concentrations in axon in response to synapse"]},{Plot[HH[t][[4]],{t,0,.1},Frame->True,FrameLabel->{"t","V"},GridLines->Automatic,PlotRange->All,PlotLabel->"Voltage in axon in response to synapse"]}}]

Differential equations Vp=1

eqs = { m'[t] == αm[V[t]] (1-m[t]) - βm[V[t]] m[t],
          n'[t] == αn[V[t]] (1-n[t]) - βn[V[t]] n[t],
        h'[t] == αh[V[t]] (1-h[t]) - βh[V[t]] h[t],
       V'[t] == -(( 1)/Cm)( gNabar (((m[t]))^3) h[t] (V[t] - ENa) +
gKbar (((n[t]))^4) (V[t]-EK) + gl (V[t]-El)+pulse[t,1])}

TableForm[eqs]

Initial conditions

conds = {V[0]== Vr, n[0]==0.315, m[0]==0.042, h[0]==0.608}

HH[t_] ={m[t]/0.042,n[t]/0.315,h[t]/0.608,V[t]/Vr} /.NDSolve[Flatten[{eqs,conds}],{m[t],n[t],h[t],V[t]},{t,0,10}][[1]]

pltVp01=GraphicsGrid[{{Plot[HH[t][[1;;3]],{t,0,10},Frame->True,FrameLabel->{"t","Na,K,Cl"},GridLines->Automatic,PlotLabel->"Ion concentrations in axon in response to synapse"]},{Plot[HH[t][[4]],{t,0,.1},Frame->True,FrameLabel->{"t","V"},GridLines->Automatic,PlotRange->All,PlotLabel->"Voltage in axon in response to synapse"]}}]

Show[{pltVp01,pltVp05,pltVp10,pltVp20}]

Function encapsulating the Hodgkin-Huxley model

Instead of the above separate instructions, build a function to find the solution for various input voltages. Note that model parameters use the globally defined values.

HodgkinHuxley[Vp_] := Module[{HH,pltVp},
eqs = { m'[t] == αm[V[t]] (1-m[t]) - βm[V[t]] m[t],
          n'[t] == αn[V[t]] (1-n[t]) - βn[V[t]] n[t],
        h'[t] == αh[V[t]] (1-h[t]) - βh[V[t]] h[t],
       V'[t] == -(( 1)/Cm)( gNabar (((m[t]))^3) h[t] (V[t] - ENa) +
gKbar (((n[t]))^4) (V[t]-EK) + gl (V[t]-El)+Evaluate[pulse[t,Vp]])};
conds = {V[0]== Vr, n[0]==0.315, m[0]==0.042, h[0]==0.608};
HH[t_] ={m[t]/0.042,n[t]/0.315,h[t]/0.608,V[t]/Vr} /. NDSolve[Flatten[{eqs,conds}],{m[t],n[t],h[t],V[t]},{t,0,60}][[1]];
pltVp=GraphicsGrid[{{Plot[Evaluate[HH[t][[1;;3]]],{t,0,1},Frame->True,FrameLabel->{"t","Na,K,Cl"},GridLines->Automatic,PlotLabel->"Ion concentrations due to synapse",PlotLegends->{"Na","K","Cl"}],
Plot[HH[t][[4]],{t,0,0.1},Frame->True,FrameLabel->{"t","V"},GridLines->Automatic,PlotRange->All,PlotLabel->"Voltage response to synapse"]},
{Plot[Evaluate[HH[t][[1;;3]]],{t,0,10},Frame->True,FrameLabel->{"t","Na,K,Cl"},GridLines->Automatic,PlotLabel->"Ion concentrations due to synapse",PlotLegends->{"Na","K","Cl"}],Plot[HH[t][[4]],{t,0,10},Frame->True,FrameLabel->{"t","V"},GridLines->Automatic,PlotRange->All,PlotLabel->"Voltage response to synapse"]}},ImageSize->Large];
Return[pltVp];
]

HodgkinHuxley[20.]

Plot[pulse[t,10.],{t,0,0.01}]

Also define a function that returns the voltage response for some input pulse

HHV[P_]:= Module[{HH,m,n,h,V,pltVp,sol},
eqs = { m'[t] == αm[V[t]] (1-m[t]) - βm[V[t]] m[t],
          n'[t] == αn[V[t]] (1-n[t]) - βn[V[t]] n[t],
        h'[t] == αh[V[t]] (1-h[t]) - βh[V[t]] h[t],
       V'[t] == -(( 1)/Cm)( gNabar (((m[t]))^3) h[t] (V[t] - ENa) +
gKbar (((n[t]))^4) (V[t]-EK) + gl (V[t]-El)+Evaluate[P])};
conds = {V[0]== Vr, n[0]==0.315, m[0]==0.042, h[0]==0.608};
sol =NDSolve[Flatten[{eqs,conds}],{m[t],n[t],h[t],V[t]},{t,0,60}][[1]];
Return[ V[t]/Vr /. sol];
]

Plot[Evaluate[HHV[10. HeavisidePi[1000.` (-0.001`+t)]]],{t,0,60},PlotRange->All]

Plot[Evaluate[HHV[Sin[t]]],{t,0,60},PlotRange->All]

Neurons in series

Consider a first neuron that receives a step function pulse (synapse model)

V1[t_]=HHV[10. HeavisidePi[1000.` (-0.001`+t)]]

V2[t_] = HHV[5V1[t]]

Plot[{V1[t],V2[t]},{t,0,.5},Frame->True,FrameLabel->{"t","V1,V2"},PlotLabel->"Voltage in two neurons in series",GridLines->Automatic,PlotLegends->Automatic,PlotRange->All]

Space-time Hodgkin-Huxley

Define a wave equation operator

WaveEq[u_] := D[u,{t,2}] - (a^2)  D[u,{x,2}]

Verify the general form of the solution, a sum of a forward and backward propagating wave with speed a.

Simplify[WaveEq[f[x-a t] + g[x + a t]]]

Consider a succession of ion channel positions along an axon. First compute the response of the first portion

V1[t_]=HHV[10. HeavisidePi[1000.` (-0.001`+t)]]

Define values for transmission of action potential along an axon

L=0.2; a=21; J=10; Δ=L/J

HHV[P_,t0_]:= Module[{HH,m,n,h,V,pltVp,sol},
eqs = { m'[t] == αm[V[t]] (1-m[t]) - βm[V[t]] m[t],
          n'[t] == αn[V[t]] (1-n[t]) - βn[V[t]] n[t],
        h'[t] == αh[V[t]] (1-h[t]) - βh[V[t]] h[t],
       V'[t] == -(( 1)/Cm)( gNabar (((m[t]))^3) h[t] (V[t] - ENa) +
gKbar (((n[t]))^4) (V[t]-EK) + gl (V[t]-El)+Evaluate[P])};
conds = {V[t0]== Vr, n[t0]==0.315, m[t0]==0.042, h[t0]==0.608};
sol =NDSolve[Flatten[{eqs,conds}],{m[t],n[t],h[t],V[t]},{t,t0,60}][[1]];
Return[ V[t]/Vr /. sol];
]

V2[t_]=HHV[V1[t-Δ/a],Δ/a]

V3[t_]=HHV[V2[t-Δ/a],2Δ/a]

Plot[{V1[t],V2[t],V3[t]},{t,2Δ/a,0.5},PlotRange->All,Frame->True,GridLines->Automatic,PlotLegends->Automatic]

Problems

With this problem set, you will make the transition to full scale modeling within the life sciences. With the expertise built up in previous assignments on coding, mathematical background on differential equations and calculus, and exposure to the physiological aspects of neuroscience, you can now begin to leverage this skill set to answer realistic questions.

Problem 1: Action potential threshold

Exercise 1, p. 226

To determine the pulse strength threshold, first modify the HodgkinHuxley function to return just the voltage response, expressed in volts.

HodgkinHuxleyV[Vp_] := Module[{VHH,pltVp},
eqs = { m'[t] == αm[V[t]] (1-m[t]) - βm[V[t]] m[t],
          n'[t] == αn[V[t]] (1-n[t]) - βn[V[t]] n[t],
        h'[t] == αh[V[t]] (1-h[t]) - βh[V[t]] h[t],
       V'[t] == -(( 1)/Cm)( gNabar (((m[t]))^3) h[t] (V[t] - ENa) +
gKbar (((n[t]))^4) (V[t]-EK) + gl (V[t]-El)+Evaluate[pulse[t,Vp]])};
conds = {V[0]== Vr, n[0]==0.315, m[0]==0.042, h[0]==0.608};
VHH[t_] =V[t]/. NDSolve[Flatten[{eqs,conds}],{m[t],n[t],h[t],V[t]},{t,0,60}][[1]];
pltVp=Plot[VHH[t],{t,0,0.02},Frame->True,FrameLabel->{"t","V"},GridLines->Automatic,PlotRange->{All,{-80,40}},PlotLabel->"Voltage response to synapse"];
Return[pltVp];
]

Next, construct a table of plots for various pulse strengths and display it.

Manipulate[HodgkinHuxleyV[Vp],{Vp,2.,20.}]

plots = Table[HodgkinHuxleyV[Vp],{Vp,2.,18.,2.}];
GraphicsGrid[ArrayReshape[plots,{6,3}]]

HodgkinHuxleyV[20.]

Problem 2: Membrane capacitance

Exercise 2, p. 226

Problem 3: Leakage conductance

Exercise 3, p. 226

Problem 4: Metabolic sodium concentration

Exercise 4, p. 227

Problem 5: Neurons in series

Place two or more neurons in series and study the transmission of an action potential. Think about the role of multiple dendrites versus a single axon. Attempt to study a network of neurons, with multiple neuron axons linked to a the dendrites of another neuron. Reflect upon the implications for complex signal processing. This is an open-ended question for which up to 4 extra credit course points can be awarded (essentially a complete, thoughtful analysis to this question is equivalent to a regular full homework assignment).