source: project/release/4/sundials/tags/1.4/tests/hh.scm @ 25541

Last change on this file since 25541 was 25541, checked in by Ivan Raikov, 9 years ago

sundials release 1.4

File size: 5.1 KB
Line 
1;;
2;; Hodgkin-Huxley model
3;;
4
5(use mathh sundials srfi-4)
6
7(define neg -)
8(define pow expt)
9
10(define TEND  500.0)
11
12                           
13;; Model parameters
14
15;; model parameters
16(define (I_stim t) 10)
17(define C_m       1)
18(define E_Na      50)
19(define E_K       -77)
20(define E_L       -54.4)
21(define gbar_Na   120)
22(define gbar_K    36)
23(define g_L       0.3)
24
25;; Rate functions
26
27(define (amf v)   (* 0.1    (/ (+ v 40)  (- 1.0 (exp (/ (neg (+ v 40)) 10))))))
28(define (bmf v)   (* 4.0    (exp (/ (neg (+ v 65)) 18))))
29(define (ahf v)   (* 0.07   (exp (/ (neg (+ v 65)) 20))))
30(define (bhf v)   (/ 1.0    (+ 1.0 (exp (/ (neg (+ v 35)) 10)))))
31(define (anf v)   (* 0.01   (/ (+ v 55) (- 1 (exp (/ (neg (+ v 55)) 10))))))
32(define (bnf v)   (* 0.125  (exp (/ (neg (+ v 65)) 80))))
33
34;; State functions
35
36(define (minf v) (* 0.5 (+ 1 (tanh (/ (- v v1) v2)))))
37(define (winf v) (* 0.5 (+ 1 (tanh (/ (- v v3) v4)))))
38(define (lamw v) (* phi (cosh (/ (- v v3) (* 2 v4)))))
39                           
40;; Model equations
41
42(define (rhs t yy)
43
44  (let ((v (f64vector-ref yy 0))
45        (m (f64vector-ref yy 1))
46        (h (f64vector-ref yy 2))
47        (n (f64vector-ref yy 3)))
48
49    ;; transition rates at current step
50    (let ((am  (amf v))
51          (an  (anf v))
52          (ah  (ahf v))
53          (bm  (bmf v))
54          (bn  (bnf v))
55          (bh  (bhf v))
56
57          (g_Na (* gbar_Na  (* h (pow m 3))))
58          (g_K  (* gbar_K   (pow n 4))))
59     
60      (let (
61
62            ;; currents
63            (I_Na   (* (- v E_Na) g_Na))
64            (I_K    (* (- v E_K)  g_K))
65            (I_L    (* g_L  (- v E_L))))
66                 
67        (let (
68              ;; state equations
69              (dm (- (* am (- 1 m))  (* bm m)))
70              (dh (- (* ah (- 1 h))  (* bh h)))
71              (dn (- (* an (- 1 n))  (* bn n)))
72              (dv (/ (- (I_stim t) I_L I_Na I_K) C_m))
73              )
74   
75          (f64vector dv dm dh dn)
76         
77          )))
78    ))
79
80
81(define (rhs/unsafe t yy yp _)
82  (let ((v (pointer-f64-ref yy))
83        (m (pointer-f64-ref (pointer+-f64 yy 1)))
84        (h (pointer-f64-ref (pointer+-f64 yy 2)))
85        (n (pointer-f64-ref (pointer+-f64 yy 3))))
86
87
88    ;; transition rates at current step
89    (let ((am  (amf v))
90          (an  (anf v))
91          (ah  (ahf v))
92          (bm  (bmf v))
93          (bn  (bnf v))
94          (bh  (bhf v))
95
96          (g_Na (* gbar_Na  (* h (pow m 3))))
97          (g_K  (* gbar_K   (pow n 4))))
98     
99      (let (
100
101            ;; currents
102            (I_Na   (* (- v E_Na) g_Na))
103            (I_K    (* (- v E_K)  g_K))
104            (I_L    (* g_L  (- v E_L))))
105                 
106        (let (
107              ;; state equations
108              (dm (- (* am (- 1 m))  (* bm m)))
109              (dh (- (* ah (- 1 h))  (* bh h)))
110              (dn (- (* an (- 1 n))  (* bn n)))
111              (dv (/ (- (I_stim t) I_L I_Na I_K) C_m))
112              )
113   
114          (pointer-f64-set! yp dv)
115          (pointer-f64-set! (pointer+-f64 yp 1) dm)
116          (pointer-f64-set! (pointer+-f64 yp 2) dh)
117          (pointer-f64-set! (pointer+-f64 yp 3) dn)
118         
119          )))
120    ))
121
122
123
124(define (main)
125 
126  (let ((yy (f64vector -65  0.052 0.596 0.317)) ;; v m h n
127
128        ;; Integration limits
129        (t0  0.0)
130        (tf  TEND)
131        (dt  1e-2))
132   
133    ;; CVODE initialization
134    (let ((solver (cvode-create-solver
135                   t0 yy rhs 
136                   tstop: tf
137                   abstol: 1e-4
138                   reltol: 1e-4)))
139
140      ;; In loop, call CVodeSolve, print results, and test for error.
141     
142      (let recur ((tnext (+ t0 dt)) (iout 1))
143
144        (let ((flag  (cvode-solve solver tnext)))
145          (if (negative? flag) (error 'main "CVODE solver error" flag))
146
147          ;;(print-results solver tnext)
148
149          (if (< tnext tf)
150              (recur (+ tnext dt) (+ 1 iout)))
151          ))
152     
153
154      (let ((yy (cvode-yy solver)))
155        (let ((v (f64vector-ref yy 0))
156              (m (f64vector-ref yy 1))
157              (h (f64vector-ref yy 2))
158              (n (f64vector-ref yy 3)))
159          (assert (< (abs (- v 13.1818250584672)) 1e-12))
160          (assert (< (abs (- m 0.988216867260496)) 1e-12))
161          (assert (< (abs (- h 0.152871810132604)) 1e-12))
162          (assert (< (abs (- n 0.680166108403409)) 1e-12))
163          ))
164     
165      (cvode-destroy-solver solver)
166     
167      )))
168
169(define (main/unsafe)
170 
171  (let ((yy (f64vector -65  0.052 0.596 0.317)) ;; v m h n
172
173        ;; Integration limits
174        (t0  0.0)
175        (tf  TEND)
176        (dt  1e-2))
177   
178    ;; CVODE initialization
179    (let ((solver (cvode-create-solver/unsafe
180                   t0 yy rhs/unsafe 
181                   tstop: tf
182                   abstol: 1e-4
183                   reltol: 1e-4)))
184
185      ;; In loop, call CVodeSolve, print results, and test for error.
186     
187      (let recur ((tnext (+ t0 dt)) (iout 1))
188
189        (let ((flag  (cvode-solve solver tnext)))
190          (if (negative? flag) (error 'main "CVODE solver error" flag))
191
192          (if (< tnext tf)
193              (recur (+ tnext dt) (+ 1 iout)))
194          ))
195     
196
197      (let ((yy (cvode-yy solver)))
198        (let ((v (f64vector-ref yy 0))
199              (m (f64vector-ref yy 1))
200              (h (f64vector-ref yy 2))
201              (n (f64vector-ref yy 3)))
202          (assert (< (abs (- v 13.1818250584672)) 1e-12))
203          (assert (< (abs (- m 0.988216867260496)) 1e-12))
204          (assert (< (abs (- h 0.152871810132604)) 1e-12))
205          (assert (< (abs (- n 0.680166108403409)) 1e-12))
206          ))
207     
208      (cvode-destroy-solver solver)
209     
210      )))
211
212
213(define (print-results solver t)
214  (let ((yy (cvode-yy solver)))
215    (printf "~A ~A ~A ~A ~A ~A ~A ~A~%" 
216            t 
217            (f64vector-ref yy 0)
218            (f64vector-ref yy 1)
219            (f64vector-ref yy 2)
220            (f64vector-ref yy 3)
221            (cvode-get-last-order solver)
222            (cvode-get-num-steps solver)
223            (cvode-get-last-step solver)
224            )))
225     
226     
227(main)
228(main/unsafe)
Note: See TracBrowser for help on using the repository browser.