ode
![[Picture of an egg]](egg.jpg)
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.
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 = awhere 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' = 1cannot 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.
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:
The returned dispatcher procedure can take one of the following arguments:
| 'add-guard! |
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:
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
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. |
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:
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 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
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. |
| Prints the values of the quantities specified by the set-print!procedure. |
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 |
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.
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 ...) |
An extension for rate equations in ODE systems.
The recognized syntax is:
(rate ((=> SRC DEST RATE-EXPR) | (<= DEST SRC RATE-EXPR) ...)
An extension for specifying Hodgkin-Huxley-type dynamics in ODE systems.
An extension for creating Scheme procedures in place right-hand side expressions in ODE systems.
;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; 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))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ;; 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)))
;;
;; 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)));;
;; 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 it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A full copy of the GPL license can be found at <http://www.gnu.org/licenses/>.