- 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

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

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
| ||||||||||

'defun! |
Returns a procedure of the form Argument This procedure modifies the given equation environment. | ||||||||||

'depgraph |
Returns a procedure of the form | ||||||||||

'env-extend! |
Returns a procedure of the form Argument
Argument Argument This procedure modifies the given equation environment. | ||||||||||

eqdef! |
Returns a procedure of the form Argument Optional argument This procedure modifies the given equation environment. | ||||||||||

eval-const |
Returns a procedure of the form 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 | ||||||||||

exam |
Returns a procedure of the form | ||||||||||

'set-indep! |
Returns a procedure of the form This procedure modifies the given equation environment. | ||||||||||

set-print! |
Returns a procedure of the form Argument The optional argument The optional argument This procedure modifies the given equation environment. | ||||||||||

'set-timestep! |
Returns a procedure of the form This procedure modifies the given equation environment. | ||||||||||

step! |
Returns a procedure of the form Argument Optional argument This procedure modifies the given equation environment. | ||||||||||

'system |
Returns a procedure of the from | ||||||||||

'sysname |
Returns a procedure of the from |

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 Argument Argument |

'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

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/>.