source: project/release/4/sundials/trunk/tests/hh.scm @ 30930

Last change on this file since 30930 was 30930, checked in by Ivan Raikov, 6 years ago

sundials: small tweaks to example model parameters

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