Transcript Slide 1

Mathematical Models in
Neural and Neuroendocrine
Systems
Richard Bertram
Department of Mathematics and
Programs in Neuroscience and
Molecular Biophysics
Florida State University
What’s the Point?





Integrate biological data
Writing down equations often identifies holes in
biological knowledge
Use of math models helps to determine if
data are logically consistent
Making testable predictions
Models force simplification, which can help
in the identification of essential elements
How Much Detail Should Go Into a
Model?
Some reasons that simpler is better:
 More detail means more undetermined
parameters
 More equations decrease our ability to
comprehend the behavior of the model
 Added detail can be misleading. It can give the
false impression that all the components are
necessary for the behavior
How Should One Calibrate a Biological
Model?
In engineering, one does experiments to
determine the values of each parameter
In biological models it is rare to be able to
measure all or most of the parameters…
Calibration Method for Biological
Models
Make model predictions
Test them in the lab
Modify model
Single Cell Models
Example: Neural Excitability
Hodgkin-Huxley Model (1952)

dV
Conservation
Cm
  I Na  I K  I leak  I ap
of charge
dt
Gating
variables

dn
 n (V )  n  /  n (V )
dt
dm
 m (V )  m  /  m (V )
dt
dh
 h (V )  h  /  h (V )
dt
Hodgkin-Huxley Model (1952)
Strengths
dV
Cm
 I Na  I K  I leak  I ap 
dt
Accurate
Essential
elements
Nothing
extra
dn
 n (V )  n  /  n (V )
dt
dm
 m (V )  m  /  m (V )
dt
dh
 h (V )  h  /  h (V )
dt
Hodgkin-Huxley Model (1952)
Strengths
dV
Cm
 I Na  I K  I leak  I ap 
dt
Accurate
Essential
elements
Nothing
extra
dn
 n (V )  n  /  n (V )
dt
dm
 m (V )  m  /  m (V )
dt
dh
 h (V )  h  /  h (V )
dt
Difficulties
Nonlinear
High
dimensionality
Morris-Lecar Model (1981)
dV
Cm
 ( I Ca  I K  I leak  I ap )
dt
dw
 w (V )  w /  w (V )
dt
The m ODE is removed
by assuming that the
Ca2+ channels activate
instantaneously.
The h ODE is removed
since the Ca2+ channels
don’t inactivate.
The resulting system is two-dimensional or planar.
This greatly simplifies the analysis.
How do we analyze this?
Phase plane: Study the dynamics in the Vw-plane
rather than V or w versus time
Nullclines: Determine the curves along which one of
the time derivatives is 0
Steady states: At the intersections of the two nullclines
both derivatives are 0, so the system is at rest
Direction arrows: The nullclines divide up the plane,
and the direction of flow in each region can be
determined
Morris-Lecar with Iap=0
Morris-Lecar with Iap=0
Morris-Lecar with Iap=100 pA
Increasing the applied current translates the V-nullcline
upward, so that the nullcline intersection is on the
middle branch of the V-nullcline. The steady state is
unstable. The stable solution is a limit cycle.
Morris-Lecar with Iap=100 pA
Motion along the limit cycle is periodic.
Hopf Bifurcation
This transition from a stable steady state to a
stable limit cycle through variation of a parameter
is called a Hopf bifurcation. It is one way in which
periodic motion can arise from a previously
stationary system, or vice versa.
Morris-Lecar with Iap=250 pA
Increasing the applied current more puts the intersection
on the right branch of the V-nullcline. Here the steady
state is again stable. The system has gone through a
second Hopf bifurcation.
Bifurcation Diagrams
These different behaviors and how they
change with variation of a parameter can
be summarized with a Bifurcation Diagram.
Stationary Branches
For each value of the bifurcation parameter Iap
plot the V value of the steady state solution. If
stable make the curve solid; if unstable make the
curve dashed. Stability changes at bifurcation points.
Periodic Branches
Next plot the maximum and minimum V values of
periodic solutions. Stability changes at two
Saddle Node of Periodics (SNP) bifurcations.
Bistability
For some values of Iap the system is bistable. For
a single value of the parameter, a stable steady
state and a stable limit cycle coexist.
Bistability in the Phase Plane
The unstable limit cycle (dashed) separates the
basin of attraction of the steady state from that
of the stable limit cycle. It is the separatrix.
Time Courses of Bistable System
Example 2: Calcium Dynamics
Because Ca2+ plays a large role in cells, modelers
have used differential equations to describe the
changes over time of the concentration of free Ca2+
in the cytosol, endoplasmic reticulum, and in some
cases in the mitochondria.
We will derive equations for the free concentration
in the cytosol and the ER, starting with the cytosol.
The Ins and Outs of Calcium
Jin =influx of calcium
Jout =efflux of calcium
Cac = free calcium
concentration
d Cac
 f c  J in  J out 
dt
fc =fraction of calcium ions not bound to buffer proteins
Calcium Flux Terms
Influx of calcium is through ion channels, with current
ICa:
J in   I Ca
 converts current to flux
Efflux of calcium is through calcium pumps in the
plasma membrane. We can approximate this pumping
as proportional to the amount of free calcium in the
cytosol:
J out  k pmcaCac
kpmca is the pump rate
Altogether Now…
d Cac
  f c I Ca  k pmca Cac 
dt
This is just an expression of conservation of mass
Calcium Equations with an ER
d Cac
 f c J in  J PMCA  J leak  J IP 3  J SERCA 
dt
d Ca ER
 f ER  J SERCA  J leak  J IP 3 
dt
  Vc V
ER
ER Fluxes
Jleak is calcium leakage from the ER into the cytosol.
This is driven by the calcium concentration difference,
so:
J leak  kleak (CaER  Cac )
JSERCA is calcium pumping from the cytosol into the
ER through SERCA pumps. For simplicity, we assume
that the pump flux is proportional to the level of
cytosolic calcium:
J SERCA  kSERCACac
JIP3 is flux through IP3 channels, when activated.
Typical Time Courses
Mean Field Models
What is Mean Field?
These models are often used to describe the
time dynamics of a large population of cells,
where single-cell models would be
inappropriate. With a mean field model a single
variable would be used to describe, for example,
the mean firing rate of a population of neurons.
Using single-cell models here would result in
thousands of variables and equations.
Example: Circadian Prolactin Rhythm
Daily PRL rhythm occurs for first 10 days of
pregnancy in rats.
What is the Mechanism for This?
Our hypothesis is that the negative feedback of
DA from the arcuate nucleus onto lactotrophs in
the pituitary together with the delayed positive
feedback of PRL onto DA neurons serves as a
rhythm generator. This feedback loop provides
delayed negative feedback to lactotrophs.
The Dopamine Equation
d DA
 Td  k p PRL2  q DA
dt
DA is mean firing rate of DA neurons
PRL is the mean PRL secretory activity
Td is tonic stimulatory drive
kpPRL2 is the delayed positive feedback of
prolactin
-qDA is first-order recovery
The Prolactin Equation
rp
d PRL

 q PRL
2
dt
1  DA
The first term has DA in the denominator,
reflecting the inhibitory influence of DA on
PRL secretion.
The PRL-Oxytocin Feedback Loop
We and others have shown that Oxytocin (OT) stimulates
lactotrophs, and that PRL inhibits OT neurons. Unlike the
stimulatory action of PRL on DA neurons, the inhibitory
action of PRL on OT neurons is rapid.
Modified Prolactin Equation
d PRL rp  oOT

 q PRL
2
dt
1  DA
The first term now has OT in the numerator,
reflecting its stimulatory effect on PRL.
The Oxytocin Equation
To
d OT


q
OT
2
dt
1  ko PRL
The first term has PRL in the denominator,
since PRL inhibits OT neurons. There is no
time delay.
Daily Pulse of VIP
VIPergic neurons from the suprachiasmatic nucleus
(SCN) innervate DA neurons of the arcuate nucleus
and provide inhibition. The VIP activity is high in the
morning, and low other times throughout the day.
Modified Dopamine Equation
d DA
 Td  k p PRL2  q DA  rvVIP  DA
dt
The last term, which is negative, reflects
the inhibitory action of VIP.
The VIP is non-zero for 3 hours in the early
morning.
VIP
Model Reproduces the PRL Rhythm
First 6 days
Day 3
Uses of the Model
We have used the PRL rhythm model to help
in the design of experiments. Several predictions
have been made and tested, and based on the
results we have modified the model. It is a work
in progress!
Our goal with this mean field model is not to
make quantitative predictions, but to make
qualitative predictions. This is all that one can
expect when thousands of cells are represented
by a single variable!
Thanks to …
Marc Freeman
Marcel Egli
The NIH and NSF for financial support