Description

Numerical solver framework for systems of first-order differential equations (ODE).

Ivan Raikov

Version

• 4.0 Introduced ode-stalin extension
• 3.3 Bug fix in the pattern matcher of ode-rate and ode-bpr
• 3.2 Bug fix in the coupled chains code of ode-hhsm
• 3.1 Changes to the solver interface for compatibility with other implementations of Scheme
• 3.0 Added support for configurable single/double floating point format
• 2.9 Some bug fixes to the bpr extension
• 2.8 Meta-file bug fix
• 2.7 Meta-file bug fix
• 2.6 Added support for user-specified error tolerances
• 2.5 Added the ode-ctranslator extension; some reorganization to support foreign-lambdas as rhs procedures
• 2.4 Added the ode-pr and ode-hhsp extensions
• 2.3 Restructuring of the interp1d extension into waveform extension and new interp1d
• 2.2 Bug fixes in the ode macros and in the meta file
• 2.1 Added support for vector state functions
• 2.0 Changes to the ode runtime interface
• 1.4 Bug fixes in the meta file
• 1.3 Added sysname method and ode-transform macro
• 1.2 Yet more solver bug fixes and adjustments
• 1.1 Some bug fixes and fine tuning of the numerical solvers
• 1.0 Initial release

Usage

(require-extension ode)

ode.egg

Documentation

The ode library provides a framework that can represent several different families of first order differential equations. It includes three distinct algorithms for the numerical solution of differential equations: Euler, Runge-Kutta-Fehlberg, and 4th order Adams-Moulton.

Definitions

A differential equation is an equation that involves an unknown function and its derivatives. A differential equation is ordinary if the unknown function depends on only one independent variable, often denoted t. The order of the differential equation is the order of the highest-order derivative in the equation.

A system of equations is a set of equations that contain multiple variables. If the equations are dependent on one another they are said to be coupled. A solution is any function satisfying the equations. An initial value problem (IVP) is a system of ordinary differential equations together with specified initial conditions on the unknown functions and their derivatives, all of which are given at the same value of the independent variable.

In practical problems, the solution of a differential equation is usually not known in terms of an explicit formula, and the solution must be numerically approximated.

A numerical method for a system of ordinary differential equations produces approximate solutions at particular points using only the operations of arithmetic and functional evaluations of the given equations. ode solves first order initial value problems of the form:

x' = f(t,x,y,...,z)
y' = g(t,x,y,...,z)
.
.
.
z' = h(t,x,y,...,z)
given the initial values for each dependent variable at the initial value of the independent variable,
x(a) = b
y(a) = c
.
.
.
z(a) = d
t = a
where a,b,c,...,d are constants.

In the equations above, f, g, ..., and h can be expressed with the operators +,-,*,/ and the mathematical functions described in the next section. In order to solve a system of equations with ode, it must be possible to solve explicitly for each derivative. ode can't solve an equation in which the first derivative is implicitly defined; for instance, the system

y' + y' = 1
cannot be solved with ode.

All numerical integration algorithms calculate an approximate solution at successive discrete values of the independent variable. The interval between any two successive points is called the time step (usually denoted dt) and it may be constant or variable. In general, as dt becomes smaller, the solution becomes more accurate.

ODE Interface

The ode interface is structured as an "object" -- a dispatcher procedure that takes a method name as a symbol, and returns the procedure that implements the respective operation. The ode object is created by procedure make-ode:

procedure: make-ode -> ODE

The returned dispatcher procedure can take one of the following arguments:

A guard is a boolean expression that always must hold true at the end of an integration cycle. It can be used to test the invariants of a system of equations being developed.

add-guard! returns a procedure of the form LAMBDA ODE-ENV -> (LAMBDA EXPR -> UNDEFINED) that when invoked with a particular equation environment, returns a procedure that takes a boolean expression and adds it to the list of guard expressions for that equation environment. This procedure modifies the given equation environment.

'defun!

Returns a procedure of the form LAMBDA ODE-ENV -> (LAMBDA NAME FORMALS BODY -> UNDEFINED) that when invoked with a particular equation environment, defines a new function within that environment. It is an error to define the same function more than once.

Argument NAME is the name of the function and it must be a symbol. Argument FORMALS is the list of formal arguments to the function, and it must be a list of symbols. Argument BODY is the body of the function, and it must be an arithmetic expression that only references the formal arguments.

This procedure modifies the given equation environment.

'depgraph

Returns a procedure of the form LAMBDA ODE-ENV -> DIGRAPH that when invoked with a particular equation environment, returns the dependency graph of the different quantities in the system. Please see the digraph egg for details about the graph API used.

'env-extend!

Returns a procedure of the form LAMBDA ODE-ENV -> (LAMBDA NAME TYPE INITIAL [RHS] -> UNDEFINED) that when invoked with a particular equation environment, defines a new quantity within that environment. It is an error to define the same quantity more than once.

Argument NAME is the name of the quantity and it must be a symbol. Argument TYPE is a symbol that has one of the following values:

 'const constant during integration 'state state variable 'asgn a quantity that is assigned a directly evaluated arithmetic expression during integration 'prim primitive (a Scheme value) other quantity type defined by an extension module

Argument INITIAL must be a numeric value in the case of constants and state variables, the symbol '(none) in the case of assigned quantities, or a Scheme value in the case of primitives.

Argument RHS is optional and allowed only in the case of states or assigned quantities and is an arithmetic expression that represents the right-hand side of the equation for that quantity.

This procedure modifies the given equation environment.

eqdef!

Returns a procedure of the form LAMBDA ODE-ENV -> (LAMBDA NAME RHS [DF] -> UNDEFINED) that when invoked with a particular equation environment, defines the right-hand side of the equation for the specified quantity. It is an error to specify the right-hand side of the same quantity more than once, or to refer to a quantity that has not been defined with env-extend!.

Argument RHS is required only in the case of states or assigned quantities and is an arithmetic expression that represents the right-hand side of the equation for that quantity.

Optional argument DF can specify if this is a difference equation, when the quantity being referred is a system state. DF can take the values of 'd (the default), to indicate a differential equation, or 'q, to indicate a difference equation.

This procedure modifies the given equation environment.

eval-const

Returns a procedure of the form LAMBDA ODE-ENV EXPR -> NUMBER that when invoked with an equation environment and an expression evaluates the expression and returns its numeric value.

The expression must refer only to constants and primitive procedures that have already been defined in the equation environment.

make-const-env

Returns a procedure of the form LAMBDA ODE-ENV -> ENV that when invoked with an equation environment returns an environment that contains only the constants and primitive functions from the given equation environment.

exam

Returns a procedure of the form LAMBDA ODE-ENV -> (LAMBDA NAME -> UNDEFINED) that when invoked with a particular equation environment, prints information about the specified quantity. It is an error to specify an undefined quantity.

'set-indep!

Returns a procedure of the form LAMBDA ODE-ENV -> (LAMBDA NAME -> UNDEFINED) that when invoked with a particular equation environment, sets the name of the independent variable (default: t).

This procedure modifies the given equation environment.

set-print!

Returns a procedure of the form LAMBDA ODE-ENV -> PRTLIST [EVERY [FROM]] -> UNDEFINED which specifies a list of quantities to be printed during the numerical integration process.

Argument PRTLIST is a list in which every element is of the form (TAG SYM). SYM is the name of a quantity in the system, and TAG is a symbol with one of the following values: VALUE, PRIME, ABSERR, RELERR

The optional argument EVERY specifies that the quantities must be printed only when the number of iterations is a multiple of EVERY.

The optional argument FROM specifies that the quantities must be printed only when the number of iterations exceeds FROM.

This procedure modifies the given equation environment.

'set-timestep!

Returns a procedure of the form LAMBDA ODE-ENV -> (LAMBDA NAME -> UNDEFINED) that when invoked with a particular equation environment, sets the name of the timestep variable (default: dt).

This procedure modifies the given equation environment.

step!

Returns a procedure of the form LAMBDA ODE-ENV -> (LAMBDA SOLVER START STOP [INITIAL-H] -> UNDEFINED) that when invoked, performs the numerical step from START to STOP.

Argument SOLVER is the solver procedure to use; please see the next section for information about the solver interface, and files euler.scm, rkf45.scm, abm4.scm in the egg distribution for examples of different solver implementations.

Optional argument INITIAL-H can be used to specify initial time step.

This procedure modifies the given equation environment.

'system

Returns a procedure of the from LAMBDA NAME -> ODE-ENV that when invoked, returns an environment that represents a system of equations.

'sysname

Returns a procedure of the from LAMBDA ODE-ENV -> NAME that when invoked with an ODE environment argument, returns the name for that system of equations.

procedure: ode:env-copy:: ODE-ENV -> ODE-ENV
This procedure makes a copy of the given equation environment.

Arithmetic Expressions

Each right-hand side expression is an s-expression in which the quantities of the system and the built-in arithmetic functions and constants are recognized.

The following built-in functions and constants are currently defined:

Built-in Arithmetic Functions
+ - * / pow neg abs atan asin acos sin cos ceiling floor exp ln sqrt tan cosh sinh tanh hypot gamma lgamma log10 log2 log1p ldexp cube round ceiling floor max min
Built-in Constants
E 1/E E^2 E^PI/4 LOG2E LOG10E LN2 LN3 LNPI LN10 1/LN2 1/LN10 PI PI/2 PI/4 1/PI 2/PI 2/SQRTPI SQRTPI PI^2 DEGREE SQRT2 1/SQRT2 SQRT3 SQRT5 SQRT10 CUBERT2 CUBERT3 4THRT2 GAMMA1/2 GAMMA1/3 GAMMA2/3 PHI LNPHI 1/LNPHI EULER E^EULER SIN1 COS1 ZETA3

In addition, a right hand side expression can include an if-statement of the form IF COND IFTRUE IFFALSE. The given condition must be a boolean expression; if it is true, then the IFTRUE expression is evaluated and its value returned; otherwise the IFFALSE expression is evaluated and its value returned.

Guard Expressions

Guard expressions are a sequence of boolean expressions that must hold true at the end of each integration cycle. Guard expressions can include all of the built-in arithmetic functions, system quantities, as well as the following boolean functions: > < <= >= = and or

Solver Interface

The step! procedure invokes a particular solver procedure with an argument ODE-RUNTIME which is a dispatcher procedure that recognizes the following messages:

 'order Returns the number of states in the system 'svec-value-idx The index at which the quantity value is stored in the numeric vector associated with each quantity during an integration run 'svec-deriv-idx The index at which the derivative value is stored in the numeric vector associated with each quantity during an integration run 'svec-abserr-idx) The index at which the absolute error value is stored in the numeric vector associated with each quantity during an integration run 'svec-relerr-idx The index at which the relative error value is stored in the numeric vector associated with each quantity during an integration run 'is-state? Given a symbol, returns true if the symbol is the name of a state in the given equation system 'is-asgn? Given a symbol, returns true if the symbol is the name of a state equation quantity in the given equation system 'eval-poset-foreach Given a procedure of the form LAMBDA ORDER LST evaluates the procedure once for each quantity in the equation system, grouped according to their evaluation order. Argument ORDER is the evaluation order index. Argument LST is the list of all quantities that have evaluation order ORDER. Each element of the list is of the form (SYM . RHS which are respectively the name and the equation right-hand side for that quantity. 'eval-poset-fold The same as eval-poset-foreach but performs a fold-like operation. The given procedure must take an additional argument for the fold accumulator value. 'eval-env-ref Returns the value of a quantity in the current evaluation environment. 'eval-env-set! Sets the value of a quantity in the current evaluation environment. 'solve-env-ref Returns the value of a quantity in the current solve environment. 'solve-env-set! Sets the value of a quantity in the current solve environment. 'solve-env-state-fold Performs a fold-like operation on the states in the solve environment 'solve-env-state-foreach Iterates over the states in the solve environment 'solve-env-asgn-fold Performs a fold-like operation on the assigned quantities in the solve environment 'eval-rhs Evaluates the given expression in the evaluation environment for the current system. 'relmax Returns the upper bound relative error for the current system. 'relmin Returns the lower bound relative error for the current system. 'absmax Returns the upper bound absolute error for the current system. 'absmin Returns the lower bound absolute error for the current system. 'maxerror Computes the current highest relative and absolute errors in the system. 'hierror? Given the return value of the maxerror procedure, returns true if the maximum relative or absolute error is exceeded. 'lowerror? Given the return value of the maxerror procedure, returns true if the error is less than the minimum relative and absolute errors. 'print Prints the values of the quantities specified by the set-print!procedure.

Macros

macro: (with-ode-system name declarations . body)

A helper macro for creating ODE systems and simulating them with the provided numerical solvers. NAME is the name of the variable that will be bound to the ODE system object when the statements in BODY are executed.

The declaration syntax is the following:

 ('prim id value) define Scheme values ('state id '= expr) define state quantity ('indep id) define independent variable (default: t) (id '= expr) define algebraic assignment ('defun id (arg1 arg2 ...) expr) define state function ('prim id scheme-expr) define primitive function or constant (Scheme value) ('d (id) '= expr) define differential equation ('q (id) '= expr) define difference equation ('guard expr) define guard expression ('guards (expr1 ...)) define several guard expressions at once ('print prtlist) define quantities to be printed during simulation run

macro: (ode-transform sys declarations)

This macro applies the given declarations to the given ODE system, and returns the modified equation environment. It recognizes the same syntax as with-ode-system.

ode-waveform

An extension for waveform generation in ODE systems.

The recognized syntax is :

(waveform (NAME (TYPE PROPERTIES ...)))
where TYPE and PROPERTIES are of the following form:
 interp1d ((step CONST) (src LIST-OR-STREAM) [(method INTERP-METHOD)]) sine ((amplitude CONST) (rate CONST) [(phase-shift CONST) (vertical-offset CONST)]) square ((amplitude CONST) (period CONST) [(duty-cycle CONST) (vertical-offset CONST)]) sawtooth ((amplitude CONST) (period CONST) [(vertical-offset CONST)]) compose (def1 ...) sum (def1 ...) multiply (def1 ...)

ode-rate

An extension for rate equations in ODE systems.

The recognized syntax is:

(rate ((=> SRC DEST RATE-EXPR) | (<= DEST SRC RATE-EXPR) ...)

ode-hhs

An extension for specifying Hodgkin-Huxley-type dynamics in ODE systems.

ode-lambda

An extension for creating Scheme procedures in place right-hand side expressions in ODE systems.

Examples

Simple Examples

;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Compute Euler's number
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;

(require-extension ode-macros)
(require-extension ode-abm4)

(with-ode-system euler
((state y = 1)
(d (y) = y)
(print ((value t) (value y))))

(((ode 'step!) euler) ode:abm4 0.0 1.0))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Example 203a in _Numerical Methods for Ordinary Differential
;; Equations_, by J.C. Butcher.
;;
;; The exact solution of this system is:
;;
;;  y1(x) = cos(x)
;;  y2(x) = sin(x)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(require-extension ode-macros)
(require-extension ode-rkf45)

(with-ode-system b203a
((indep x)
(state y1 = 1.0)
(state y2 = 0.0)
(d (y1) = (+ (+ (* -16.0 y1) (* 12.0 y2)) (+ (* 16.0 (cos x)) (* -13.0 (sin x)))))
(d (y2) = (+ (+ (* 12.0  y1) (* -9.0  y2)) (+ (* -11.0 (cos x)) (* 9.0  (sin x)))))
(print ((value x) (value y1) (value y2))))

(((ode 'step!) b203a) ode:rkf45 0.0 1.0))

Hodgkin-Huxley Model

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; This model is used to calculate a neuron's membrane potential
;; assuming some initial state. The calculation is based on sodium ion
;; flow, potassium ion flow and leakage ion flow. (Hodgkin, A. L. and
;; Huxley, A. F. (1952) "A Quantitative Description of Membrane
;; Current and its Application to Conduction and Excitation in Nerve"
;; Journal of Physiology 117: 500-544)
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(require-extension ode-macros)
(require-extension ode-rkf45)
(require-extension ode-waveform)

(with-ode-system hodgkin-huxley
(;; model states
(state v = -65)
(state m = 0.052)
(state h = 0.596)
(state n = 0.317)

;; input stimulus
(waveform (I_stim (sine (amplitude 10))))

;; model parameters
(const C_m      = 1)
(const E_Na     = 50)
(const E_K      = -77)
(const E_L      = -54.4)
(const gbar_Na  = 120)
(const gbar_K   = 36)
(const g_L      = 0.3)

;; rate functions
(defun amf (v)   (* 0.1    (/ (+ v 40)  (- 1.0 (exp (/ (neg (+ v 40)) 10))))))
(defun bmf (v)   (* 4.0    (exp (/ (neg (+ v 65)) 18))))
(defun ahf (v)   (* 0.07   (exp (/ (neg (+ v 65)) 20))))
(defun bhf (v)   (/ 1.0    (+ 1.0 (exp (/ (neg (+ v 35)) 10)))))
(defun anf (v)   (* 0.01   (/ (+ v 55) (- 1 (exp (/ (neg (+ v 55)) 10))))))
(defun bnf (v)   (* 0.125  (exp (/ (neg (+ v 65)) 80))))

;; transition rates at current step
(am = (amf v))
(an = (anf v))
(ah = (ahf v))
(bm = (bmf v))
(bn = (bnf v))
(bh = (bhf v))

(g_Na = (* gbar_Na  (* h (pow m 3))))
(g_K  = (* gbar_K   (pow n 4)))

;; currents
(I_Na   = (* (- v E_Na) g_Na))
(I_K    = (* (- v E_K)  g_K))
(I_L    = (* g_L  (- v E_L)))

;; state equations
(d (m)  = (- (* am (- 1 m))  (* bm m)))
(d (h)  = (- (* ah (- 1 h))  (* bh h)))
(d (n)  = (- (* an (- 1 n))  (* bn n)))
(d (v)  = (/ (- Istim I_L I_Na I_K) C_m))

(print ((value t)  (value v)
(value g_Na)  (value I_Na) (value g_K) (value I_K)
(value am) (value bm) (value ah) (value bh)
(value an) (value bn)
(value n))))

(let ((hh1 (ode:waveform-transformer hodgkin-huxley)))
(((ode 'step!) hh1) ode:rkf45 0.0 40.0)))

Hodgkin-Huxley Dynamics as Rate Equations

;;
;; This model is the deterministic variant of:
;;
;; _Ion Channel Stochasticity May Be Critical in Determining the
;; Reliability and Precision of Spike Timing_
;;
;; Elad Schneidman, Barry Freedman, and Idan Segev
;;
;; Neural Computation 10, 1679-1703 (1998)
;;

(require-extension ode-macros)
(require-extension ode-rkf45)
(require-extension ode-rate)
(require-extension format-graph)

(define format-graph (make-format-graph 'vcg))

(with-ode-system schneidman98
(;; model states
(state v = -65)

;; external inputs to the model
(const I_src = 10)

;; model parameters
(const C_m      = 1)    ;; Specific membrane capacitance, 1 uF/cm^2
(const T        = 6.3)  ;; Temperature, degC

(const E_Na     = 115)  ;; Sodium reversal potential, mV
(const E_K      = -12)  ;; Potassium reversal potential, mV
(const E_L      = 10.6) ;; Leakage reversal potential, mV

(const gamma_Na = 20)   ;; Sodium channel conductance, pS
(const gamma_K  = 20)   ;; Potassium channel conductance, pS

(const gbar_Na  = 120)   ;; Maximal sodium conductance, mS/cm^2
(const gbar_K   = 36)    ;; Maximal potassium conductance, mS/cm^2
(const g_L      = 0.3)   ;; Leakage conductance

(const D_Na     = 60)   ;; Sodium ion channel density, channels/um^2
(const D_K      = 18)   ;; Potassium ion channel density, channels/um^2
(const Area     = 200)  ;; Patch area, um^2

;; total number of K channels
(const n_K  = (* Area D_K))
;; total number of Na channels
(const n_Na = (* Area D_Na))

(const i_n0 = (* n_K 0.2))
(const i_n1 = (* n_K 0.2))
(const i_n2 = (* n_K 0.2))
(const i_n3 = (* n_K 0.2))
(const i_n4 = (* n_K 0.2))

(state n0 = i_n0)
(state n1 = i_n1)
(state n2 = i_n2)
(state n3 = i_n3)
(state n4 = i_n4)

(const i_m0h0 = (* n_Na 0.47))
(const i_m1h0 = (* n_Na 0.01))
(const i_m2h0 = (* n_Na 0.01))
(const i_m3h0 = (* n_Na 0.01))
(const i_m0h1 = (* n_Na 0.47))
(const i_m1h1 = (* n_Na 0.01))
(const i_m2h1 = (* n_Na 0.01))
(const i_m3h1 = (* n_Na 0.01))

(state m0h0 = i_m0h0)
(state m1h0 = i_m1h0)
(state m2h0 = i_m2h0)
(state m3h0 = i_m3h0)
(state m0h1 = i_m0h1)
(state m1h1 = i_m1h1)
(state m2h1 = i_m2h1)
(state m3h1 = i_m3h1)

;; rate functions
(defun anf (v)   (* 0.01   (/ (- 10 v)  (- (exp (/ (- 10 v) 10)) 1))))
(defun bnf (v)   (* 0.125  (exp (neg (/ v 80)))))
(defun amf (v)   (* 0.1    (/ (- 25 v)  (- (exp (/ (- 25 v) 10)) 1))))
(defun bmf (v)   (* 4.0    (exp (neg (/ v 18)))))
(defun ahf (v)   (* 0.07   (exp (/ (neg v) 20))))
(defun bhf (v)   (/ 1.0    (+ 1 (exp (/ (- 30 v) 10)))))

;; transition rates at current step
(am = (amf v))
(an = (anf v))
(ah = (ahf v))
(bm = (bmf v))
(bn = (bnf v))
(bh = (bhf v))

;; Transition rates for the K channels, according to
;; the following kinetic scheme:
;;
;; [n0] <-> [n1] <-> [n2] <-> [n3] <-> [n4]
;;
(rate (-> n0 n1 (* 4 an))
(-> n1 n2 (* 3 an))
(-> n2 n3 (* 2 an))
(-> n3 n4 (* 1 an))
(-> n1 n0 (* 1 bn))
(-> n2 n1 (* 2 bn))
(-> n3 n2 (* 3 bn))
(-> n4 n3 (* 4 bn)))

;; check that the total # of K channels remains
;; constant
(guard (= n_K (round (+ n0 n1 n2 n3 n4))))

;; check that we have a positive # of K channels in
;; each state
(guards  (> n0 0.0)
(> n1 0.0)
(> n2 0.0)
(> n3 0.0)
(> n4 0.0))

;; Transition rates for the Na channels,
;; according to the following kinetic scheme:
;;
;; [m0h1] <-> [m1h1] <-> [m2h1] <-> [m3h1]
;;   /\         /\         /\         /\
;;   ||         ||         ||         ||
;;   \/         \/         \/         \/
;; [m0h0] <-> [m1h0] <-> [m2h0] <-> [m3h0]
;;
(rate ;; transitions in row h1
(-> m0h1 m1h1 (* 3 am))
(-> m1h1 m2h1 (* 2 am))
(-> m2h1 m3h1 am)
(-> m1h1 m0h1 bm)
(-> m2h1 m1h1 (* 2 bm))
(-> m3h1 m2h1 (* 3 bm))
;; transitions in row h0
(-> m0h0 m1h0 (* 3 am))
(-> m1h0 m2h0 (* 2 am))
(-> m2h0 m3h0 am)
(-> m1h0 m0h0 bm)
(-> m2h0 m1h0 (* 2 bm))
(-> m3h0 m2h0 (* 3 bm))
;; transitions between rows
(-> m0h1 m0h0 bh)
(-> m1h1 m1h0 bh)
(-> m2h1 m2h0 bh)
(-> m3h1 m3h0 bh)
(-> m0h0 m0h1 ah)
(-> m1h0 m1h1 ah)
(-> m2h0 m2h1 ah)
(-> m3h0 m3h1 ah))

;; check that the total # of Na channels remains
;; constant
(guard (= n_Na (round (+ m0h0 m1h0 m2h0 m3h0 m0h1 m1h1 m2h1 m3h1))))

;; check that we have a positive # of Na channels in
;; each state
(guards  (> m0h0 0.0)
(> m1h0 0.0)
(> m2h0 0.0)
(> m3h0 0.0)
(> m0h1 0.0)
(> m1h1 0.0)
(> m2h1 0.0)
(> m3h1 0.0))

;; Potassium and sodium membrane conductances as
;; functions of V and t
(g_K   = (* gamma_K   (/ n4 Area)))
(g_Na  = (* gamma_Na  (/ m3h1 Area)))

(I_Na   = (* g_Na (- v E_Na)))
(I_K    = (* g_K  (- v E_K)))
(I_L    = (* g_L  (- v E_L)))
(I_stim = (* 1.0 I_src))

(d (v)  = (/ (- I_stim I_L I_Na I_K) C_m))

(print ((value t) (value v)
(value g_K) (value I_K) (value g_Na) (value I_Na)
(value an) (value bn) (value am) (value bm)
(value n0) (value n1) (value n2) (value n3) (value n4)
(value m0h0) (value m1h0) (value m2h0) (value m3h0)
(value m0h1) (value m1h1) (value m2h1) (value m3h1))) )

(with-output-to-file "schneidman98.vcg"
;; output the equation dependency graph
(lambda () (format-graph (current-output-port) ((ode 'depgraph) schneidman98))))
(let ((schneidman98-1 (ode:rate-transformer schneidman98)))
;; run the numerical integrator for the specified length of time
(((ode 'step!) schneidman98-1) ode:rkf45 0.0 250.0 0.01)))

Crustacean Stomatogastric Ganglion (STG) Model

;;
;; This model is the deterministic variant of:
;;
;; _ Whole Cell Stochastic Model Reproduces the Irregularities Found
;; in the Membrane Potential of Bursting Neurons_
;;
;; Pedro V. Carelli, Marcelo B. Reyes, Jose C. Sartorelli, and Reynaldo D. Pinto
;;
;; J Neurophysiol 94: 1169 1179, 2005.
;;

(require-extension ode-macros)
(require-extension ode-euler)
(require-extension ode-hhs)
(require-extension format-graph)

(define format-graph (make-format-graph 'vcg))

(with-ode-system carelli05
(;; membrane potential
(state v  = -65)

;; Ca2+ concentration, uM
(const Ca0 = 0.5)
(state Ca  = Ca0)
(const f   = 14.96) ;; uM/nA
(const tau_Ca = 200) ;; ms

;; external inputs to the model
(const I_src = 20)

;; model parameters
(const C_m      = 1)    ;; Specific membrane capacitance, uF/cm^2

(const Area = (* 0.628 1e-3)) ;; Membrane area, cm^2

;; Reversal potentials, mV
(const E_Na     = 50)
(const E_Kd     = -80)
(const E_KCa    = -80)
(const E_A      = -80)
(const E_L      = -50)
(const E_H      = -20)

;; Ionic conductances, mS/cm^2

(const g_L      = 0.01)

(ionic conductance (Na
(gbar  200)
(gamma 3)
(delta 1)
(initial-m 0.9)
(initial-h 0.1)
(minf  (/ 1 (+ 1 (exp (/ (+ v 25.5) -5.29)))))
(hinf  (/ 1 (+ 1 (exp (/ (+ v 48.9) 5.18)))))
(taum  (- 2.64 (/ 2.52 (+ 1 (exp (/ (+ v 120) -25))))))
(tauh  (* (/ 1.34 (+ 1 (exp (/ (+ v 62.9) -10))))
(+ 1.5 (/ 1 (+ 1 (exp (/ (+ v 34.9) 3.6)))))))))

(ionic conductance (CaT
(gbar 2.5)
(gamma 3)
(delta 1)
(initial-m 0.1)
(initial-h 0.1)
(minf  (/ 1 (+ 1 (exp (/ (+ v 27.1) -7.2)))))
(hinf  (/ 1 (+ 1 (exp (/ (+ v 32.1) 5.5)))))
(taum  (- 43.4 (/ 42.6 (+ 1 (exp (/ (+ v 68.1) -20.5))))))
(tauh  (- 210  (/ 179.6 (+ 1 (exp (/ (+ v 55) -16.9))))))))

(ionic conductance (CaS
(gbar  4)
(gamma 3)
(delta 1)
(initial-m 0.1)
(initial-h 0.1)
(minf  (/ 1 (+ 1 (exp (/ (+ v 33) -8.1)))))
(hinf  (/ 1 (+ 1 (exp (/ (+ v 60) 6.2)))))
(taum  (+ 2.8  (/ 14  (+ (exp (/ (+ v 27) 10)) (exp (/ (+ v 70) -13))))))
(tauh  (+ 120  (/ 300 (+ (exp (/ (+ v 55) 9)) (exp (/ (+ v 65) -16))))))))

(ionic conductance (A
(gbar 50)
(gamma 3)
(delta 1)
(initial-m 0.1)
(initial-h 0.1)
(minf  (/ 1 (+ 1 (exp (/ (+ v 27.2) -8.7)))))
(hinf  (/ 1 (+ 1 (exp (/ (+ v 56.9) 4.9)))))
(taum  (- 23.2 (/ 20.8 (+ 1 (exp (/ (+ v 32.9) -15.2))))))
(tauh  (- 77.2  (/ 58.4 (+ 1 (exp (/ (+ v 38.9) -26.5))))))))

;; In cases where delta=0, there is no inactivation,
;; and therefore the h expressions are omitted

(ionic conductance (KCa
(gbar  5)
(gamma 4)
(delta 0)
(initial-m 0.1)
(minf  (* (/ Ca (+ 3 Ca)) (/ 1 (+ 1 (exp (/ (+ v 28.3) -12.6))))))
(taum  (- 180.6 (/ 150.2 (+ 1 (exp (/ (+ v 46) -22.7))))))))

(ionic conductance (Kd
(gbar 100)
(gamma 4)
(delta 0)
(initial-m 0.1)
(minf  (/ 1 (+ 1 (exp (/ (+ v 12.3) -11.8)))))
(taum  (- 14.4  (/ 12.8  (exp (/ (+ v 28.3) -19.2)))))))

(ionic conductance (H
(gbar 0.01)
(gamma 1)
(delta 0)
(initial-m 0.1)
(minf  (/ 1 (+ 1 (exp (/ (+ v 75) 5.5)))))
(taum  (/ 2.2 (+ (exp (/ (+ v 169.7) -11.6)) (exp (/ (- v 26.7) 14.3)))))))

;; Assuming external Ca2+ concentration of 3000 uM
(E_Ca = (* 12.5  (ln (/ 3000 Ca))))

;; Ionic currents
(I_Na   = (* g_Na  (- v E_Na)))
(I_CaT  = (* g_CaT (- v E_Ca)))
(I_CaS  = (* g_CaS (- v E_Ca)))
(I_A    = (* g_A   (- v E_A)))
(I_KCa  = (* g_KCa (- v E_KCa)))
(I_Kd   = (* g_Kd  (- v E_Kd)))
(I_H    = (* g_H   (- v E_H)))
(I_L    = (* g_L   (- v E_L)))

;; stimulus current
(I_stim = I_src)

;; Calcium buffer
(d (Ca) = (* (neg (/ 1 tau_Ca)) (+ (* f (+ I_CaT I_CaS)) Ca (neg Ca0))))

(d (v)  = (/ (- I_stim I_L I_Na I_CaT I_CaS I_A I_KCa I_Kd I_H) C_m))

(print ((value t) (value v) (value I_stim))) )

(let ((carelli05-1 (ode:hhs-transformer carelli05)))
(with-output-to-file "carelli05.vcg"
;; output the equation dependency graph
(lambda () (format-graph (current-output-port) ((ode 'depgraph) carelli05-1))))
;; run the numerical integrator for the specified length of time
(((ode 'step!)  carelli05-1) ode:euler 0.0 500.0 5e-3)))

Copyright 2007 Ivan Raikov and the Okinawa Institute of Science and Technology

This program is free software: you can redistribute it and/or modify