Transcript Lecture 3

Introduction to Mathematical
Methods in Neurobiology:
Dynamical Systems
Oren Shriki
2009
Neural Excitability
References about neurons as
electrical circuits:
• Koch, C. Biophysics of Computation, Oxford
Univ. Press, 1998.
• Tuckwell, HC. Introduction to Theoretical
Neurobiology, I&II, Cambridge UP, 1988.
References about neural excitability:
• Rinzel & Ermentrout. Analysis of neural excitability and
oscillations. In, Koch & Segev (eds): Methods in Neuronal
Modeling, MIT Press, 1998.
Also as “Meth3” on www.pitt.edu/~phase/
• Borisyuk A & Rinzel J. Understanding neuronal dynamics by
geometrical dissection of minimal models. In, Chow et al, eds:
Models and Methods in Neurophysics (Les Houches Summer
School 2003), Elsevier, 2005: 19-72.
• Izhikevich, EM. Dynamical Systems in Neuroscience: The
Geometry of Excitability and Bursting. MIT Press, (soon).
Neural Excitability
• One of the main properties of neurons is
that they can generate action potentials.
• In the language of non-linear dynamics we
call such systems “excitable systems”.
• The first mathematical model of a neuron
that can generate action potentials was
proposed by Hodgkin and Huxley in 1952.
• Today there are many models for neurons
of diverse types, but they all have common
properties.
A Typical Neuron
Intracellular Recording
Electrical Activity of a Pyramidal Neuron
The Neuron as an Electric Circuit
• The starting point in models of neuronal
excitability is describing the neuron as an
electrical circuit.
• The model consists of a dynamical equation for
the voltage, which is coupled to equations
describing the ionic conductances.
• The non-linear behavior of the conductances is
responsible for the generation of action potentials
and determines the properties of the neuron.
Hodgkin & Huxley Model )1952(
• The model considers an isopotential
membrane patch )or a “point neuron”(.
• It contains two voltage-gated ionic
currents, K and Na, and a constantconductance leak current.
The Equivalent Circuit
External solution
Internal solution
K and Na Currents During an
Action Potential
Positive and
Negative
Feedback
Cycles During
an Action
Potential
The Hodgkin-Huxley Equations
dV
C
  I ion V , w  I (t )
dt
I ion V , m, h, n  
g Na m h(V  VNa )  g K n (V  VK )  g L (V  VL )
3
4
dm /dt m(V)-m/τ m
dh/dt  h(V)-h/τ h
dn/dt  n(V)-n/τ n
The Hodgkin & Huxley Framework
dV
C
  I ion V , W1 , , WN   I (t )
dt
Each gating variable obeys the following dynamics:

Wi , V   Wi 
dWi

dt
 i V 
i
- Time constant

- Represents the effect of temperature
The Hodgkin & Huxley Framework
The current through each channel has the form:
I j  g j j V ,W1,,WN (V  V j )
gj
j
- Maximal conductance (when all channels are open)
- Fraction of open channels
(can depend on several W variables).
The Temperature Parameter Φ
• Allows for taking into account different
temperatures.
• Increasing the temperature accelerates the
kinetics of the underlying processes.
• However, increasing the temperature does not
necessarily increase the excitability. Both
increasing and decreasing the temperature can
cause the neuron to stop firing.
• A phenomenological model for Φ is:
Temp6.3/10
3
Simplified Versions of the HH Model
• Models that generate action potentials can
be constructed with fewer dynamic
variables.
• These models are more amenable for
analysis and are useful for learning the
basic principles of neuronal excitability.
• We will focus on the model developed by
Morris and Lecar.
The Morris-Lecar Model (1981)
• Developed for studying the barnacle muscle.
• Model equations:
dV
C
  I ion V , w  I app (t )
dt

dw
w V   w

dt
 w V 
I ion V , w  gCa m (V )(V  VCa )  gK w(V  VK )  gL (V  VL )
Morris-Lecar Model
• The model contains K and Ca currents.
• The variable w represents the fraction of open K
channels.
• The Ca conductance is assumed to behave in an
instantaneous manner.
m (V )  0.51  tanhV  V1  / V2 
w (V )  0.51  tanhV  V3  / V4 
 w V   1/ coshV  V3  /2V4 
Morris-Lecar Model
• A set of parameters for example:
V1  1.2
g Ca  4
VCa  120
V2  18
gK  8
VK  84
V3  2
gL  2
VL  60
V4  30
μF
C  20 2
cm
  0.04
Morris-Lecar Model
• Voltage dependence
of the various
parameters (at long
times):
Nullclines
• We first analyze the nullclines of the model (the
curves on which the time derivative is zero).
• The voltage equation gives:
 Iion V , w  I app  0
• To see what this curve looks like, let us assume that
Iapp=0 and fix w on a certain value, w0.
• This is the instantaneous I-V relation of the system:
I ion V , w  I inst V ; w0 
 g Ca m (V )(V  VCl )  g K w0 (V  VK )  g L (V  VL )
Nullclines
• This curve is N-shaped and increasing w
approximately moves the curve upward:
• We are interested in the zeros of this curve for each
value of w.
Nullclines
• The nullcline of w is simply the activation curve:
Nullclines
• The nullclines in the phase plane:
Nullclines
• The steady state solution of the system is the point
where the nullclines intersect.
• In the presence of an applied current,
the steady-state satisfies:
 Iion V , w (V )  I app  0
Iion V , w (V )  I app
I ss V 
Steady-State Current I ss V   Iion V , w (V )
1500
I ss (V,w  )
1000
500
0
-500
-80 -60 -40 -20
0
V
20 40 60 80
General Properties of Nullclines
• The nullclines segregate the phase plane into
regions with different directions for a trajectory’s
vector flow.
• A solution trajectory that crosses a nullcline does so
either vertically or horizontally.
• The nullclines’ intersections are the steady states or
fixed points of the system.
The Excitable Regime
• We say that the system is in the excitable
regime when the following conditions are
met:
– It has just one stable steady-state in the
absence of an external current.
– An action potential is generated after a large
enough brief current pulse.
Effect of Brief Pulses on the ML Neuron
• A brief current pulse is equivalent to setting the initial
voltage at some value above the rest state.
• Below are responses for the following initial voltages:
-20, -12, -10, 0 mV
Order of Events
• If the pulse is large enough, autocatalysis starts,
the solution’s trajectory heads rightward toward the
V-nullcline.
• After the trajectory crosses the V -nullcline
(vertically, of course) it follows upward along the
nullcline’s right branch.
• After passing above the knee it heads rapidly
leftward (downstroke).
Order of Events
• The number of K+ channels open reaches a
maximum during the downstroke, as the w-nullcline
is crossed.
• Then the trajectory crosses the V -nullcline’s left
branch (minimum of V) and heads downward,
returning to rest (recovery).
Dependence of Action Potential Amplitude
on the Initial Conditions
• The amplitude of the action potential depends on the
initial conditions.
• This contradicts the traditional view that the action
potential is an all-or-none event with a fixed amplitude.
• If we plot the peak V vs. the
size of the pulse or initial
condition V0, we get a
continuous curve.
• The curve becomes steep
only at low temperature (Φ).
Dependence of Action Potential Amplitude
on the Initial Conditions
• For small Φ, w is very slow and the flow in the phase
plane is close to horizontal, since V is relatively much
faster than w:
dw  dW   dV 

 /
  0 ( )
dV  dt   dt 
• In this case, it takes fine tuning of the initial conditions
to evoke the graded responses.
Dependence of Action Potential
Amplitude on the Initial Conditions
• The analysis above suggests a surprising
prediction:
• If the experimentally observed action potentials look
like all-or-none events, they may become graded if
recordings are made at higher temperatures.
• This experiment was suggested by FitzHugh and
carried out by Cole et al. in 1970.
• It was found that if recordings in squid giant axon
are made at 38◦C instead of, say, 15◦C then action
potentials do not behave in an all-or none manner.
Repetitive Firing at the ML Model
• In response to prolonged stimulation the cell fires
spikes in a repetitive manner.
• From a dynamical point of view, the rest state is
no longer stable and the system starts to oscillate.
• Typically, when the applied current is increased,
the firing rate increases as well.
• We will examine to cases:
– Generation of oscillations with a non-zero frequency
– Generation of oscillations with a zero frequency
Linearization of the Morris-Lecar Model
• At what current will the steady state loose stability?
• To answer this, we will use linear stability theory.
• The model equations are:
dV
C
  I ion V , w  I (t )
dt

dw
w V   w

dt
 w V 
I ion V , w  gCa m (V )(V  VCa )  gK w(V  VK )  gL (V  VL )
Linearization of the Morris-Lecar Model
• Denote the steady state by:
(V , w )
• Consider the effects of a small perturbation:
V V  x
ww y
• Using linearization we obtain:
dx
 ax  by
dt
dy
 cx  dy
dt
Linearization of the Morris-Lecar Model
• The Jacobian matrix is:
1 I ion (V , w)

a b   C V
 c d     w

   w V
1 I ion (V , w)
 C
w


w


 V ,w 
• The eigenvalues are given by:
  (a  d )  (ad  bc)       0
2
1, 2
2

1
2
     4
2

Linearization of the Morris-Lecar Model
  ad  bc 
 ad 


I ion (V , w)
 wC
V

1 I ion (V , w)
C
V
V , w 
I ion (V , w) w
w
V
 w

V ,w 
Linearization of the Morris-Lecar Model
I ion (V , w)
m (V )
 g Ca
(V  VCa )  g Ca m (V )  g K w  g L
V
V
 V  V1 
m (V ) 1 


tanh
V
2 V
 V2 

2
2
sinh
x
cosh
x

sinh
x

 
2
tanh x   

1

tanh
x
 
2
cosh x
 cosh x 

m (V )
1 
2  V  V1 

 

1

tanh

V
2V2 
V
 2 
Linearization of the Morris-Lecar Model
I ion (V , w)
 g K (V  VK )
w

w (V )
1 
2  V  V3 

 

1

tanh

V
2V4 
V
 4 
Linearization of the Morris-Lecar Model

 
1 
2  V  V1 

 (V  VCa )  g Ca m (V )

1

tanh
 g Ca

 wC 
2V2 
V
 2 

1 
2  V  V3 

 
 g K w  g L  g K (V  VK )
1

tanh

2V4 
V
 4  


1
1 
2  V  V1 

 (V  VCa )  g Ca m (V )  g K w  g L 
    g Ca
1

tanh

C 
2V2 

 V2  


w
Onset of Oscillations
• For the parameter values we are using the eigenvalues
have a negative real part.
• There are two options for loosing stability:
• Case 1: One of the eigenvalues becomes zero. In this
case detJ=0 at the bifurcation.(Saddle node bifurcation).
• Case 2: The real part of both eigenvalues becomes
zero at the bifurcation. After the bifurcation there are
two complex eigenvalues with positive real part. In this
case, TrJ=0 at the bifurcation.(Hopf bifurcation).
Onset of Oscillations with Finite Frequency
• The determinant is given by:



I ion (V , w)
 wC
V

I ion (V , w) w
w
V

V ,w 
• When substituting the steady-state values we obtain:


 
dI ss
 wC dV
• For these parameters values the steady-state I-V curve
is monotonic and the loss of stability can happen only
through Hopf bifurcation.
Hopf Bifurcation
• After the bifurcation the trace is positive:

1 I ion (V , w )
 C

V
w
1 I ion (V , w )
 C
V

0

w
• This condition can be interpreted as the
autocatalysis rate being faster than the negative
feedback’s rate.
• When this condition is met we obtain a limit cycle
with a finite frequency.
• The frequency of the repetitive firing is proportional
to the imaginary part of the eigenvalues.
Hopf Bifurcation
• We can estimate the current at which repetitive
firing will start by plotting TrJ as a function of the
applied current and estimating when it is zero.
• We will do it numerically.
• For each value of the applied current we draw the
nullclines and estimate the steady-state values of V
and w (the intersection of the nullclines):
Nullclines for Different Values of Iapp
I = 10:10:140
1
0.8
w
0.6
0.4
I
0.2
0
-0.2
-50
0
V
50
Hopf Bifurcation
• We then estimate TrJ and plot it as a function of the
applied current:
0.2
trace
0.1
0
-0.1
-0.2
0
100
50
I
• The critical current is around I=100.
150
Hopf Bifurcation
• We can also estimate the frequency:
0.1
Im( )
0.08
0.06
0.04
0.02
0
0
50
100
150
I
• Around I=100 the imaginary part is close to 0.084.
• The frequency should be around:
0.084/(2*pi)*1000 =13.4 Hz
f-I Curve
• A plot of the f-I curve, obtained using numerical
simulations of the dynamics, confirms our analysis:
Bifurcation Diagram
• The full bifurcation diagram is:
• The system experiences a subcritical Hopf
bifurcation.
The Effect of the Applied Current on the
Nullclines and the State of the Neuron
• At low currents the neuron is excitable.
• For intermediate currents it presents oscillations.
• At high currents the neuron experiences nerve block –
there is no firing but the membrane is depolarized. The
K and Ca currents balance each other.
Three Steady-States
• We next examine the model with the same parameters
as before, except:
V3  12 mV
V4  17 mV
• In this case, the steady-state I-V curve becomes Nshaped.
• Therefore, for some values of the applied current there
are three steady-states.
• The number and location of the steady-states do not
depend on the temperature (Φ) .
• However, the stability and the onset of repetitive firing
do depend on Φ.
Steady-State Current
I ss V   I ion V , w (V )
1500
I ss (V,w  )
1000
500
0
-500
-80 -60 -40 -20
0
V
20 40 60 80
Large Φ - Bistability of Steady-States
• At large temperature (Φ) the rate of the w is not slow
enough and there are no oscillations.
• The variable w can be considered instantaneous,
w=w(V).
• The model is then reduced to one dynamic variable, V,
which simplifies the analysis.
• We will examine the model with Φ =1.
Large Φ - Bistability of Steady-States
• The middle steady-state is unstable and the upper and
lower ones are stable.
• The neuron presents bistability between two steadystates.
• The voltage can be switched from one constant level to
the other by brief pulses )“plateau behavior”(.
Small Φ - Onset of Repetitive Firing at Zero
Frequency
• We now consider the small temperature case, and
examine the model with Φ = 0.06666667.
• The phase-plane picture:
• The two unstable manifold
branches form (topologically)
a circle that has two fixedpoints on it.
• As the applied current is
increased the two points
coalesce and disappear.
• At this point a, limit cycle
with infinite period is formed.
Saddle-Node Bifurcation
• The system experiences a saddle-node bifurcation.
• The bifurcation can be illustrated schematically in the
following manner:
• After the bifurcation, the
trajectory slows down near
the “ghost” of the stable
steady-state.
• This way, an arbitrarily slow
frequency can be obtained.
I app
Saddle-Node Bifurcation
• Near the bifurcation the frequency is proportional to
the square root of the deviation from the critical
current:
f  I app  I c
• The f-I curve:
Saddle-Node Bifurcation
Voltage vs. time close to the bifurcation (Iapp=40):
Bifurcation Diagram
Intermediate Φ – Bistability of Rest State and
a Depolarized Oscillation
• At intermediate values of Φ the upper steady state is
surrounded by a stable periodic orbit.
• The phase-plane picture for Φ = 0.25 :
• In this case, there is
bistability between the lowvoltage stationary point and
depolarized repetitive firing.
Intermediate Φ – Bistability of Rest State and
a Depolarized Oscillation
• Brief current pulses can be used to switch between the
states:
Repetitive Firing – Types of Excitability
• We can now summarize the characteristic properties
of two generic types of transition from the excitable
to oscillatory mode in neuronal models:
Type I
Type II
Iss
N-shaped
Monotonic
Subthreshold
oscillations
No
Yes
Threshold
All-or-none
No distinct threshold.
excitability with strict Potentials can be
threshold
graded
Frequency at
bifurcation
Zero
Finite
Repetitive Firing – Types of Excitability
• The two cases have different f-I relations.
• In type II the frequency jumps to finite values at the
bifurcation.
• In type I the f-I curve is continuous and near the
bifurcation there is a universal square-root dependence.
• However, if there is noise in the system, the f-I curves for
both types will become qualitatively similar.
• The discontinuity of type II will disappear with noise
since the probability of firing becomes zero where it was
zero before.
• Type I’s f-I curve will also inherit a smooth foot rather
than the infinite slope at zero frequency.
Types of Excitability in the ML Model
• The first set of parameters we studied
exhibited type II excitability.
• The second set of parameters we studied
exhibited type I excitability.
Neuronal Excitability – Characteristics
of the Non-Linear Dynamics
• Excitability – the system responds to a
large enough brief pulse with an action
potential.
• Repetitive firing – the system exhibits
oscillations in response to prolonged
stimulation.
Neuronal Excitability – Main Properties
of Models
• Autocatalysis –
e.g.
– The Na current in the HH model
– The Ca current in the ML model
• Slow (enough) negative feedback–
e.g.
– K current
– Inactivation of Na current in the HH model
Back to the Hodgkin-Huxley Model
The First
Intracellular
Recording from the
Squid Axon
The Equivalent Circuit
External solution
Internal solution
Hodgkin & Huxley Model
dV
C
  I ion V , w  I (t )
dt
I ion V , m, h, n  
g Na m h(V  VNa )  g K n (V  VK )  g L (V  VL )
3
4
dm /dt m(V)-m/τ m
dh/dt  h(V)-h/τ h
dn/dt  n(V)-n/τ n
Hodgkin & Huxley Model
dm /dt m(V)-m/τ m
dh/dt  h(V)-h/τ h
dn/dt  n(V)-n/τ n
Ionic Conductances During an
Action Potential
Repetitive Firing through Hopf
Bifurcation in Hodgkin–Huxley Model
A: Voltage time courses in response to a step of constant
depolarizing current. from bottom to top: Iapp= 5, 15, 50, 100, 200
in μamp/cm2). Scale bar is 10 msec.
B: f-I curves for temperatures of 6.3,18.5, 26◦C, as marked.
Dotted curves show frequency of the unstable periodic orbits.
Repetitive Firing through Hopf
Bifurcation in Hodgkin–Huxley Model
• With the standard parameters, the HH model
exhibits type II excitability.
• The bifurcation is a Hopf bifurcation.
• In different parameter regimes the model can
also exhibit type I excitability.
Fast-Slow Dissection of the
Action Potential
• n and h are slow compared to m and V.
• Based on this observation, the system can be
dissected into two time-scales.
• This simplifies the analysis.
• For details see:
Borisyuk A & Rinzel J. Understanding neuronal
dynamics by geometrical dissection of minimal models.
In, Chow et al, eds: Models and Methods in
Neurophysics (Les Houches Summer School 2003),
Elsevier, 2005: 19-72.
Correlation between n and h
• During the action potential the variables
n and h vary together.
• Using this correlation one can construct
a reduced model.
• The first to observe this was Fitzhouh.
Conductance-Based
Models of Cortical Neurons
Conductance-Based Models of
Cortical Neurons
• Cortical neurons behave differently than the squid
axon that Hodgkin and Huxley investigated.
• Over the years, people developed several variations of
the HH model that are more appropriate for describing
cortical neurons.
• We will now see an example of a simple model which
will later be useful in network simulations.
• The model was developed by playing with the
parameters such that it has an f-I curve similar to that
of cortical neurons.
Frequency-Current Responses
of Cortical Neurons
Excitatory Neuron:
Ahmed et. al., Cerebral
Cortex 8, 462-476, 1998
Inhibitory Neurons:
Azouz et. al., Cerebral
Cortex 7, 534-545, 1997
Frequency-Current Responses
of Cortical Neurons
The experimental findings show what f-I curves
of cortical neurons are:
• Continuous – starting from zero frequency.
• Semi-Linear – above the threshold current the
curve is linear on a wide range.
How can we reconstruct this behavior in a model?
An HH Neuron with a Linear f-I Curve
Shriki et al., Neural Computation 15, 1809–1841 (2003)
Linearization of the f-I Curve
• We start with an HH neuron that has a continuous f-I
curve (type I, saddle-node bifurcation).
• The linearization is made possible by the addition of a
certain K-current called A-current.
• The curve becomes linear only when the time constant
of the A-current is slow enough (~20 msec).
• There are other mechanisms for linearizing f-I curves.
Model
Equations:
Shriki et al.,
Neural
Computation
15, 1809–1841
(2003)
Bursting at the Cellular Level
• Bursting is an oscillation mode with relatively slow
rhythmic alternations between an active phase of rapid
spiking and a phase of quiescence.
• Since its early discovery bursting has been found as a
primary mode of behavior in many nerve and endocrine
cells.
• Bursting can be a single-cell property or a network
property.
• We will now focus on a simple example of a bursting
mechanism at the cellular level by extending the MorrisLecar model.
ML Model with Ca-dependent K-current
• The model is based on the ML model with the set of
parameters that gives bistability of a rest state with a
depolarized limit cycle (intermediate Φ ).
• The Ca-dependent K-current is described by:
I K Ca  g K -Ca (1  z )(V  VK )
Ca0
z
Ca  Ca0
C a    g Ca m (V )V  VCa   Ca 
g K Ca
mS
 0.25
,
2
cm
Ca0  10,   0.005,   0.2
ML Model with Ca-dependent K-current
• During a series of action potentials Ca concentration
slowly builds up and increases the Ca-dependent Kcurrent.
• When the Ca-dependent K-current is large enough, it
brings the voltage down, terminating the active phase.
• This specific mechanism
is termed “square wave”
bursting.
Fast-Slow Dissection of
Bursting Dynamics
• A general framework for analysis of bursting dynamics
is to separate all dynamic variables into two subgroup:
“fast” (X) and “slow (Y).
• In the context of bursting analysis a variable is
considered fast if it changes significantly during an
action potential, and slow if it only changes significantly
on the time scale of burst duration.
• The general equations have the form:
X  F  X , Y 
Y  G  X , Y ,
  1
Fast-Slow Dissection of
Bursting Dynamics
• The analysis of the system can be conducted in two steps:
• Step 1. First we think of the slow variables as parameters and
describe the spike generating fast subsystem for X as a
function of Y. This description involves finding steady states,
oscillatory orbits and their periods, and transitions between all
these solutions (bifurcations) as a function of Y.
• Step 2. To describe the full system we overlay the slow
dynamics on the fast system behavior. Y evolves slowly in
time according to its equations, while X is tracking its stable
states. Therefore, we must understand the direction of change
of Y at each part of the bifurcation diagram for X.