source: project/ode/trunk/extensions/ode-hhsm.scm @ 8001

Last change on this file since 8001 was 8001, checked in by Ivan Raikov, 13 years ago

Bug fixes and changes to ode-hhsm and the solvers.

File size: 15.5 KB
Line 
1;;
2;;
3;; An extension for specifying Hodgkin-Huxley type dynamics as
4;; aggregated Markov chains in ODE systems.
5;;
6;; Copyright 2007 Ivan Raikov and the Okinawa Institute of Science and Technology
7;;
8;; This program is free software: you can redistribute it and/or
9;; modify it under the terms of the GNU General Public License as
10;; published by the Free Software Foundation, either version 3 of the
11;; License, or (at your option) any later version.
12;;
13;; This program is distributed in the hope that it will be useful, but
14;; WITHOUT ANY WARRANTY; without even the implied warranty of
15;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
16;; General Public License for more details.
17;;
18;; A full copy of the GPL license can be found at
19;; <http://www.gnu.org/licenses/>.
20;;
21
22
23(require-extension ode)
24(require-extension srfi-1)
25(require-extension srfi-13)
26(require-extension environments)
27
28
29(define s+  string-append)
30(define s\  string-concatenate)
31(define c+  concatenate)
32
33(define ($ p n) (string->symbol (s+ (->string p) "_" (->string n))))
34
35(define (lookup-car k lst . rest)
36  (let-optionals rest ((default #f))
37   (let ((v (alist-ref k lst)))
38     (if v (first v) default))))
39
40
41(define (lookup-cdr k lst . rest)
42  (let-optionals rest ((default #f))
43   (let ((v (alist-ref k lst)))
44     (or v default))))
45
46                                                           
47(define (rrhs1 suf expr alst vv)
48  (match expr 
49         ((s . es)   (or (and (symbol? s) (cons s (map (lambda (x) (rrhs1 suf x alst vv)) es))) expr))
50         (id  (or (and (symbol? id) 
51                       (let ((id1 (alist-ref id alst)))
52                         (or (and id1 (first id1))
53                             (if (member id vv) (ode:error 'rrhs "duplicate variable name during renaming: " id)
54                                 id))))  id))))
55                                                           
56
57
58(define (rrhs suf expr)
59  (let* ((ss '(g gbar m h gamma delta eps minf hinf sing taum tauh taus 
60                 alpham betam alphah betah alphas betas area phi density total))
61         (vv (map (lambda (x) ($ x suf)) ss)))
62    (let ((expr1 (rrhs1 suf expr (zip ss vv) vv)))
63      expr1)))
64
65(define (check-names ion env . names)
66  (for-each (lambda (name)
67              (if (environment-includes? env name)
68                  (ode:error 'ode:hhsp-transformer "quantity " name " in ionic conductance declaration " ion
69                             "is already declared elsewhere")))
70            names))
71
72(define (check-decls ion names alst)
73  (for-each (lambda (name)
74              (if (not (alist-ref name alst))
75                  (ode:error 'ode:hhsp-transformer "required quantity " name 
76                             " is not present in ionic conductance declaration " ion)))
77            names))
78
79(define (mhssym suf)
80  (lambda (m h s) ($ (s+ "m" (number->string m) "h" (number->string h) "s" (number->string s))
81                     suf)))
82(define (mhsym suf)
83  (lambda (m h) ($ (s+ "m" (number->string m) "h" (number->string h)) suf)))
84(define (msym suf)
85  (lambda (m) ($ (s+ "m" (number->string m)) suf)))
86
87
88(define (make-state-names suf gamma . rest)
89  (define tab list-tabulate)
90  (let ((mhssym (mhssym suf))
91        (mhsym (mhsym suf))
92        (msym  (msym suf)))
93  (let-optionals rest ((delta 0) (eps 0))
94   (let ((suf (->string suf)) 
95         (m (inexact->exact (+ 1 gamma)))
96         (h (inexact->exact (+ 1 delta)))
97         (s (inexact->exact (+ 1 eps))))
98     (cond ((and (positive? eps) (positive? delta))
99            (concatenate 
100             (tab m (lambda (m) (concatenate (tab h (lambda (h) (tab s (lambda (s) (mhssym m h s))))))))))
101           ((positive? delta)
102            (concatenate (tab m (lambda (m) (tab h (lambda (h) (mhsym m h)))))))
103           (else (tab m (lambda (m) (msym m)))))))))
104
105
106(define (make-coupling chain-names state in-rate out-rate)
107  (let* ((si (map second state))
108         (stname (cond ((= 3 (length si)) mhssym)
109                       ((= 2 (length si)) mhsym)
110                       ((= 1 (length si)) msym))))
111    (let loop ((sts (map (lambda (cn) (apply (stname cn) si)) chain-names)) (lst (list)))
112      (match sts
113             ((a b . rest)
114              (loop (cons b rest) (cons* `(-> ,a ,b ,in-rate) `(-> ,b ,a ,out-rate) lst)))
115             (else  lst)))))
116     
117 
118(define (make-rates suf gamma . rest)
119  (let ((mhssym (mhssym suf))
120        (mhsym  (mhsym suf))
121        (msym   (msym suf)))
122  (define (tab i lam)  (concatenate (list-tabulate i lam)))
123  (define (mhsrate m h s mi hi si)
124    (let ((bmrate (+ 1 mi))
125          (amrate (- m mi))
126          (bhrate (+ 1 hi))
127          (ahrate (- h hi))
128          (bsrate (+ 1 si))
129          (asrate (- s si)))
130      (let ((rs (append (if (< mi m)
131                            `((-> ,(mhssym mi hi si) ,(mhssym (+ 1 mi) hi si)
132                                  (* ,amrate ,($ "alpham" suf)))
133                              (-> ,(mhssym (+ 1 mi) hi si) ,(mhssym mi hi si) 
134                                  (* ,bmrate ,($ "betam" suf))))
135                            (list))
136                        (if (< hi h)
137                            `((-> ,(mhssym mi hi si) ,(mhssym mi (+ 1 hi) si) 
138                                  (* ,ahrate ,($ "alphah" suf)))
139                              (-> ,(mhssym mi (+ 1 hi) si) ,(mhssym mi hi si) 
140                                  (* ,bhrate ,($ "betah" suf))))
141                            (list))
142                        (if (< si s)
143                            `((-> ,(mhssym mi hi si) ,(mhssym mi hi (+ 1 si) )
144                                  (* ,asrate ,($ "alphas" suf)))
145                              (-> ,(mhssym mi hi (+ 1 si)) ,(mhssym mi hi si) 
146                                  (* ,bsrate ,($ "betas" suf))))
147                            (list)))))
148        rs)))
149  (define (mhrate m h mi hi)
150    (let ((bmrate (+ 1 mi))
151          (amrate (- m mi))
152          (bhrate (+ 1 hi))
153          (ahrate (- h hi)))
154      (let ((rs (append (if (< mi m)
155                            `((-> ,(mhsym mi hi) ,(mhsym (+ 1 mi) hi) 
156                                  (* ,amrate ,($ "alpham" suf)))
157                              (-> ,(mhsym (+ 1 mi) hi) ,(mhsym mi hi) 
158                                  (* ,bmrate ,($ "betam" suf))))
159                              (list))
160                        (if (< hi h)
161                            `((-> ,(mhsym mi hi) ,(mhsym mi (+ 1 hi)) 
162                                  (* ,ahrate ,($ "alphah" suf)))
163                              (-> ,(mhsym mi (+ 1 hi)) ,(mhsym mi hi) 
164                                  (* ,bhrate ,($ "betah" suf))))
165                            (list)))))
166        rs)))
167
168    (let-optionals rest ((delta 0) (eps 0))
169      (let ((suf (->string suf)) 
170            (m (inexact->exact gamma)) 
171            (h (inexact->exact delta))
172            (s (inexact->exact eps)))
173        (cond ((and (positive? delta) (positive? eps))
174                 (concatenate 
175                  (tab (+ 1 m) (lambda (mi) 
176                                 (tab (+ 1 h) (lambda (hi) 
177                                                (list-tabulate (+ 1 s) (lambda (si) (mhsrate m h s mi hi si)))))))))
178             
179              ((positive? delta)
180               (concatenate
181                (tab (+ 1 m) (lambda (mi) (list-tabulate (+ 1 h) (lambda (hi) (mhrate m h mi hi)))))))
182
183              (else (tab  m (lambda (mi) 
184                              (let ((bmrate (+ 1 mi))
185                                    (amrate (- m mi)))
186                                `((-> ,(msym mi) ,(msym (+ 1 mi))
187                                      (* ,amrate ,($ "alpham" suf)))
188                                  (-> ,(msym (+ 1 mi)) ,(msym mi) 
189                                      (* ,bmrate ,($ "betam" suf)))))))))))))
190
191
192(define (ion-transformer new-sys env-extend! eqdef! eval-const ion alst)
193  (check-decls ion '(gamma delta minf taum density area phi) alst)
194  (let ((suffix (->string ion)))
195    (let ((g        ($ "g" suffix))
196          (gamma    ($ "gamma" suffix))
197          (delta    ($ "delta" suffix))
198          (eps      ($ "eps" suffix))
199          (minf     ($ "minf" suffix))
200          (hinf     ($ "hinf" suffix))
201          (sinf     ($ "sinf" suffix))
202          (taum     ($ "taum" suffix))
203          (tauh     ($ "tauh" suffix))
204          (taus     ($ "taus" suffix))
205          (alpham   ($ "alpham" suffix))
206          (betam    ($ "betam" suffix))
207          (alphah   ($ "alphah" suffix))
208          (betah    ($ "betah" suffix))
209          (alphas   ($ "alphas" suffix))
210          (betas    ($ "betas" suffix))
211          (density  ($ "density" suffix))
212          (total    ($ "total" suffix))
213          (area     ($ "area" suffix))
214          (phi      ($ "phi" suffix)))
215      (check-names ion new-sys g gamma delta eps
216                   minf hinf sinf taum tauh taus 
217                   alpham betam alphah betah alphas betas
218                   density total area phi)
219      (let ((seed-val  (lookup-car 'rng-seed alst))
220            (gamma-val (eval-const new-sys (lookup-car 'gamma alst)))
221            (delta-val (eval-const new-sys (lookup-car 'delta alst 0)))
222            (eps-val   (eval-const new-sys (lookup-car 'eps alst 0)))
223            (density-val (eval-const new-sys (lookup-car 'density alst)))
224            (area-val  (eval-const new-sys (lookup-car 'area alst)))
225            (phi-val   (eval-const new-sys (lookup-car 'phi alst))))
226        (if (positive? delta-val) (check-decls ion '(hinf tauh) alst))
227        (if (positive? eps-val) (check-decls ion '(sinf taus) alst))
228        (if (not (and (integer? gamma-val) (positive? gamma-val)))
229            (ode:error 'ode:hhsp-transformer 
230                       "gamma value in ionic conductance declaration " ion
231                       " must be a positive integer"))
232        (if (not (and (integer? delta-val) (or (positive? delta-val) (zero? delta-val))))
233            (ode:error 'ode:hhsp-transformer 
234                       "delta value in ionic conductance declaration " ion
235                       " must be a positive integer or zero"))
236        (if (not (and (integer? eps-val) (or (positive? eps-val) (zero? eps-val))))
237            (ode:error 'ode:hhsp-transformer 
238                       "eps value in ionic conductance declaration " ion
239                       " must be a positive integer or zero"))
240        (let* ((states        (make-state-names ion gamma-val delta-val eps-val))
241               (open-state    (last states)))
242          (apply check-names (cons ion (cons new-sys states)))
243          (let ((total-val   (* density-val area-val))
244                (nstates     (length states))
245                (g-rhs      `(/ (* ,phi ,open-state) ,area))
246                (alpham-rhs  `(/ ,minf ,taum))
247                (betam-rhs   `(/ (- 1 ,minf) ,taum))
248                (alphah-rhs  (and (positive? delta-val) `(/ ,hinf ,tauh)))
249                (betah-rhs   (and (positive? delta-val) `(/ (- 1 ,hinf) ,tauh)))
250                (alphas-rhs  (and (positive? eps-val) `(/ ,sinf ,taus)))
251                (betas-rhs   (and (positive? eps-val) `(/ (- 1 ,sinf) ,taus)))
252                (minf-rhs    (rrhs suffix (lookup-car 'minf alst)))
253                (hinf-rhs    (and (positive? delta-val) (rrhs suffix (lookup-car  'hinf alst))))
254                (sinf-rhs    (and (positive? eps-val) (rrhs suffix (lookup-car  'sinf alst))))
255                (taum-rhs    (rrhs suffix (lookup-car 'taum alst)))
256                (tauh-rhs    (and (positive? delta-val) (rrhs suffix (lookup-car  'tauh alst))))
257                (taus-rhs    (and (positive? eps-val) (rrhs suffix (lookup-car  'taus alst))))
258                (initial-m   (eval-const new-sys (lookup-car 'initial-m alst)))
259                (initial-h   (and (positive? delta-val) 
260                                  (eval-const new-sys (lookup-car 'initial-h alst))))
261                (initial-s   (and (positive? eps-val) 
262                                  (eval-const new-sys (lookup-car 'initial-s alst)))))
263            (env-extend! gamma   '(const)  gamma-val)
264            (env-extend! delta   '(const)  delta-val)
265            (env-extend! eps     '(const)  eps-val)
266            (env-extend! total   '(const)  total-val)
267            (env-extend! density '(const)  density-val)
268            (env-extend! area    '(const)  area-val)
269            (env-extend! phi     '(const)  phi-val)
270            (env-extend! g       '(asgn)   'none g-rhs)
271            (env-extend! alpham  '(asgn)   'none alpham-rhs)
272            (env-extend! betam   '(asgn)   'none betam-rhs)
273            (env-extend! minf    '(asgn)   'none minf-rhs)
274            (env-extend! taum    '(asgn)   'none taum-rhs)
275            (if (positive? delta-val) 
276                (begin
277                  (env-extend! alphah '(asgn) 'none alphah-rhs)
278                  (env-extend! betah  '(asgn) 'none betah-rhs)
279                  (env-extend! hinf   '(asgn) 'none hinf-rhs)
280                  (env-extend! tauh   '(asgn) 'none tauh-rhs)))
281            (if (positive? eps-val) 
282                (begin
283                  (env-extend! alphas '(asgn) 'none alphas-rhs)
284                  (env-extend! betas  '(asgn) 'none betas-rhs)
285                  (env-extend! sinf   '(asgn) 'none sinf-rhs)
286                  (env-extend! taus   '(asgn) 'none taus-rhs)))
287            (let ((open-pr (* (expt initial-m gamma-val)
288                              (or (and initial-h (expt initial-h delta-val)) 1.0) 
289                              (or (and initial-s (expt initial-s eps-val)) 1.0))))
290              (env-extend! open-state '(state) (eval-const new-sys `(* ,open-pr ,total-val)))
291              (for-each (lambda (x) 
292                          (let ((ival  `(/ (* ,(- 1 open-pr) ,total-val) (- ,nstates 1))))
293                            (env-extend! x  '(state) (eval-const new-sys ival))))
294                        (take states (- nstates 1)))
295              (let ((rates (make-rates ion gamma-val delta-val eps-val)))
296                (if seed-val (cons `(seed ,seed-val) rates) rates)
297                ))))))))
298
299
300(define (fold-couplings couplings chains)
301  (fold (lambda (x ax)
302          (match x  ((chain-names state in out)
303                     (let ((coupling (make-coupling chain-names state in out)))
304                       (cons coupling ax)))))
305        chains couplings))
306
307
308(define (ode:hhsm-transformer pr-transformer sys . rest)
309 (let-optionals rest ((debug #f))
310   (let ((ode (match (environment-ref sys (ode-intern 'odecore))
311                     (($ ode:quantity 'ODECORE value)  value)))
312        (new-sys (ode:env-copy sys #t)))
313    (let ((env-extend! ((ode 'env-extend!) new-sys))
314          (eqdef!      ((ode 'eqdef!) new-sys))
315          (eval-const  (ode 'eval-const)))
316      (environment-for-each sys 
317       (lambda (sym en)
318         (match en
319            ((or (('ionic 'conductance) ('name ion) . alst)
320                 (('ionic-conductance)  ('name ion) . alst))
321             (let ((rates (ion-transformer new-sys env-extend! eqdef! eval-const ion alst)))
322               (if debug (print "ode-hhsm:transformer: rates = " rates))
323               (env-extend! (gensym 'rate) '(rate) rates)))
324           
325            ((or (('coupled 'ionic 'conductance) ('name ion) . lst)
326                 (('coupled-ionic-conductance)   ('name ion) . lst))
327
328             (let* ((suffix (->string ion))
329                    (g      ($ "g" suffix))
330                    (chains
331                     (filter-map
332                      (lambda (x) 
333                        (match x (('species (species . alst))
334                                  (let ((ion-species ($ ion species)))
335                                    (cons ion-species
336                                          (ion-transformer new-sys env-extend! eqdef! eval-const ion-species alst))))
337                               (else #f)))
338                      lst))
339                    (chain-names (map car chains))
340                    (coupling-factor (lookup-car 'coupling-factor lst))
341                    (coupling-scheme (lookup-cdr 'coupling-scheme lst)))
342               
343               (cond ((and coupling-factor coupling-scheme)
344                      (ode:error 'ode:hhsm-transformer "coupled conductance declaration " ion
345                                 " has both coupling factor and coupling scheme specified"))
346                     
347                     (coupling-factor
348                      (let ((prob  (lookup-car 'prob lst))
349                            (in     (lookup-car 'in  lst))
350                            (out    (lookup-car 'out lst)))
351                        (let ((prob-val (eval-const new-sys prob)))
352                          (if (not (number? prob-val))
353                              (ode:error 'ode:hhsm-transformer "coupling probability in "
354                                         "coupled conductance declaration " lst
355                                         " must be a numeric constant"))
356                          (let ((k     (min 100 (abs (round (* 100 prob-val)))))
357                                (inq   ($ "in"  suffix))
358                                (outq  ($ "out" suffix)))
359                               
360                            (randomize)
361                            (let* ((couplings
362                                    (filter-map (lambda (x) 
363                                                  (and (<= (random 100) k) (list chain-names (car x) inq outq)))
364                                                chains))
365                                   (chains1 (fold-couplings couplings (map cdr chains))))
366                              (if debug (print "ode-hhsm:transformer: chains1 = " chains1))
367                              (env-extend! inq  '(asgn)   'none in)
368                              (env-extend! outq '(asgn)   'none out)
369                             
370                              (env-extend! (gensym 'rate) '(rate) (c+ chains1)))))))
371
372                       
373                     
374                     (coupling-scheme
375                      (let-values
376                       (((couplings coupling-rates)
377                         (let loop ((coupling-scheme coupling-scheme) (couplings (list)) 
378                                    (cpl-rates (list)) (i 0))
379                           (if (null? coupling-scheme) (values couplings cpl-rates)
380                               (let ((x (car coupling-scheme)))
381                                 (match x
382                                        (('coupling . alst)
383                                         (let ((inq    ($ (s+ "inq" (->string i)) suffix))
384                                               (outq   ($ (s+ "outq" (->string i)) suffix))
385                                               (state  (lookup-cdr 'state alst))
386                                               (in     (lookup-car 'in alst))
387                                               (out    (lookup-car 'out alst)))
388                                           (if (not state)
389                                               (ode:error 'ode:hhsm-transformer
390                                                          "coupling declaration " x " lacks state name"))
391                                           
392                                           (let ((coupling (list chain-names state in out)))
393                                             (loop (cdr coupling-scheme) (cons coupling couplings) 
394                                                   (cons* `(,inq ,in) `(,outq ,out) cpl-rates)
395                                                   (+ i 1)))))
396                                       
397                                        (else (ode:error 'ode:hhsm-transformer
398                                                         "unknown coupling declaration " x))))))))
399                       (let ((chains1 (fold-couplings couplings (map cdr chains))))
400                         (if debug (print "ode-hhsm:transformer: chains1 = " chains1))
401                         (for-each (lambda (x) (env-extend! (first x)  '(asgn)   'none (second x)))
402                                   coupling-rates)
403                         (env-extend! (gensym 'rate) '(rate) (c+ chains1)))))
404
405                     (else
406                      (ode:error 'ode:hhsm-transformer "coupled conductance declaration " ion
407                                 " has neither coupling factor nor coupling scheme specified")))
408               
409               (let ((g-rhs   `(+ . ,(map (lambda (x) ($ "g" (car x))) chains))))
410                 (env-extend! g '(asgn)   'none g-rhs))))
411           
412           
413            (else (void)))))
414      (pr-transformer new-sys)))))
415
416
Note: See TracBrowser for help on using the repository browser.