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

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

Renamed the hhsbp extension to hhsm, and changed the API to allow for any rate transformer.

File size: 11.5 KB
Line 
1;;
2;;
3;; An extension for specifying Hodgkin-Huxley type dynamics as Markov
4;; 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 environments)
26
27(define (sstr s)
28  (cond ((symbol? s) (symbol->string s))
29        (else s)))
30 
31(define ($ p n) (string->symbol (string-append (sstr p) "_" (sstr n))))
32
33(define (lookup-field k lst . rest)
34  (let-optionals rest ((default #f))
35   (let ((v (alist-ref k lst)))
36     (if v (first v) default))))
37
38                                                           
39(define (rrhs1 suf expr alst vv)
40  (match expr 
41         ((s . es)   (or (and (symbol? s) (cons s (map (lambda (x) (rrhs1 suf x alst vv)) es))) expr))
42         (id  (or (and (symbol? id) 
43                       (let ((id1 (alist-ref id alst)))
44                         (or (and id1 (first id1))
45                             (if (member id vv) (ode:error 'rrhs "duplicate variable name during renaming: " id)
46                                 id))))  id))))
47                                                           
48
49
50(define (rrhs suf expr)
51  (let* ((ss '(g gbar m h gamma delta eps minf hinf sing taum tauh taus 
52                 alpham betam alphah betah alphas betas area phi density total))
53         (vv (map (lambda (x) ($ x suf)) ss)))
54    (let ((expr1 (rrhs1 suf expr (zip ss vv) vv)))
55      expr1)))
56
57(define (check-names ion env . names)
58  (for-each (lambda (name)
59              (if (environment-includes? env name)
60                  (ode:error 'ode:hhsp-transformer "quantity " name " in ionic conductance declaration " ion
61                             "is already declared elsewhere")))
62            names))
63
64(define (check-decls ion names alst)
65  (for-each (lambda (name)
66              (if (not (alist-ref name alst))
67                  (ode:error 'ode:hhsp-transformer "required quantity " name 
68                             " is not present in ionic conductance declaration " ion)))
69            names))
70
71(define (mhssym suf)
72  (lambda (m h s) ($ (string-append "m" (number->string m) 
73                                    "h" (number->string h) 
74                                    "s" (number->string s)) 
75                     suf)))
76(define (mhsym suf)
77  (lambda (m h) ($ (string-append "m" (number->string m) "h" (number->string h)) suf)))
78(define (msym suf)
79  (lambda (m) ($ (string-append "m" (number->string m)) suf)))
80
81
82(define (make-state-names suf gamma . rest)
83  (define tab list-tabulate)
84  (let ((mhssym (mhssym suf))
85        (mhsym (mhsym suf))
86        (msym  (msym suf)))
87  (let-optionals rest ((delta 0) (eps 0))
88   (let ((suf (sstr suf)) 
89         (m (inexact->exact (+ 1 gamma)))
90         (h (inexact->exact (+ 1 delta)))
91         (s (inexact->exact (+ 1 eps))))
92     (cond ((and (positive? eps) (positive? delta))
93            (concatenate 
94             (tab m (lambda (m) (concatenate (tab h (lambda (h) (tab s (lambda (s) (mhssym m h s))))))))))
95           ((positive? delta)
96            (concatenate (tab m (lambda (m) (tab h (lambda (h) (mhsym m h)))))))
97           (else (tab m (lambda (m) (msym m)))))))))
98
99
100 
101(define (make-rates suf gamma . rest)
102  (let ((mhssym (mhssym suf))
103        (mhsym  (mhsym suf))
104        (msym   (msym suf)))
105  (define (tab i lam)
106    (concatenate (list-tabulate i lam)))
107  (define (mhsrate m h s mi hi si)
108    (let ((bmrate (+ 1 mi))
109          (amrate (- m mi))
110          (bhrate (+ 1 hi))
111          (ahrate (- h hi))
112          (bsrate (+ 1 si))
113          (asrate (- s si)))
114      (let ((rs (append (if (< mi m)
115                            `((-> ,(mhssym mi hi si) ,(mhssym (+ 1 mi) hi si)
116                                  (* ,amrate ,($ "alpham" suf)))
117                              (-> ,(mhssym (+ 1 mi) hi si) ,(mhssym mi hi si) 
118                                  (* ,bmrate ,($ "betam" suf))))
119                            (list))
120                        (if (< hi h)
121                            `((-> ,(mhssym mi hi si) ,(mhssym mi (+ 1 hi) si) 
122                                  (* ,ahrate ,($ "alphah" suf)))
123                              (-> ,(mhssym mi (+ 1 hi) si) ,(mhssym mi hi si) 
124                                  (* ,bhrate ,($ "betah" suf))))
125                            (list))
126                        (if (< si s)
127                            `((-> ,(mhssym mi hi si) ,(mhssym mi hi (+ 1 si) )
128                                  (* ,asrate ,($ "alphas" suf)))
129                              (-> ,(mhssym mi hi (+ 1 si)) ,(mhssym mi hi si) 
130                                  (* ,bsrate ,($ "betas" suf))))
131                            (list)))))
132        rs)))
133  (define (mhrate m h mi hi)
134    (let ((bmrate (+ 1 mi))
135          (amrate (- m mi))
136          (bhrate (+ 1 hi))
137          (ahrate (- h hi)))
138      (let ((rs (append (if (< mi m)
139                            `((-> ,(mhsym mi hi) ,(mhsym (+ 1 mi) hi) 
140                                  (* ,amrate ,($ "alpham" suf)))
141                              (-> ,(mhsym (+ 1 mi) hi) ,(mhsym mi hi) 
142                                  (* ,bmrate ,($ "betam" suf))))
143                              (list))
144                        (if (< hi h)
145                            `((-> ,(mhsym mi hi) ,(mhsym mi (+ 1 hi)) 
146                                  (* ,ahrate ,($ "alphah" suf)))
147                              (-> ,(mhsym mi (+ 1 hi)) ,(mhsym mi hi) 
148                                  (* ,bhrate ,($ "betah" suf))))
149                            (list)))))
150        rs)))
151
152    (let-optionals rest ((delta 0) (eps 0))
153      (let ((suf (sstr suf)) 
154            (m (inexact->exact gamma)) 
155            (h (inexact->exact delta))
156            (s (inexact->exact eps)))
157        (cond ((and (positive? delta) (positive? eps))
158                 (concatenate 
159                  (tab (+ 1 m) (lambda (mi) 
160                                 (tab (+ 1 h) (lambda (hi) 
161                                                (list-tabulate (+ 1 s) (lambda (si) (mhsrate m h s mi hi si)))))))))
162             
163              ((positive? delta)
164               (concatenate
165                (tab (+ 1 m) (lambda (mi) (list-tabulate (+ 1 h) (lambda (hi) (mhrate m h mi hi)))))))
166
167              (else (tab  m (lambda (mi) 
168                              (let ((bmrate (+ 1 mi))
169                                    (amrate (- m mi)))
170                                `((-> ,(msym mi) ,(msym (+ 1 mi))
171                                      (* ,amrate ,($ "alpham" suf)))
172                                  (-> ,(msym (+ 1 mi)) ,(msym mi) 
173                                      (* ,bmrate ,($ "betam" suf)))))))))))))
174
175
176(define (ode:hhsm-transformer pr-transformer sys . rest)
177 (let-optionals rest ((debug #f))
178  (let ((ode        (match (environment-ref sys (ode-intern 'odecore))
179                           (($ ode:quantity 'ODECORE value)  value)))
180        (new-env    (ode:env-copy sys #t)))
181    (let ((env-extend! ((ode 'env-extend!) new-env))
182          (eqdef!      ((ode 'eqdef!) new-env))
183          (eval-const  (ode 'eval-const)))
184      (environment-for-each sys 
185       (lambda (sym en)
186         (match en
187                ((('ionic 'conductance) ('name ion) . alst)
188                 (check-decls ion '(gamma delta minf taum density area phi) alst)
189                 (let ((suffix (sstr ion)))
190                   (let ((g        ($ "g" suffix))
191                         (gamma    ($ "gamma" suffix))
192                         (delta    ($ "delta" suffix))
193                         (eps      ($ "eps" suffix))
194                         (minf     ($ "minf" suffix))
195                         (hinf     ($ "hinf" suffix))
196                         (sinf     ($ "sinf" suffix))
197                         (taum     ($ "taum" suffix))
198                         (tauh     ($ "tauh" suffix))
199                         (taus     ($ "taus" suffix))
200                         (alpham   ($ "alpham" suffix))
201                         (betam    ($ "betam" suffix))
202                         (alphah   ($ "alphah" suffix))
203                         (betah    ($ "betah" suffix))
204                         (alphas   ($ "alphas" suffix))
205                         (betas    ($ "betas" suffix))
206                         (density  ($ "density" suffix))
207                         (total    ($ "total" suffix))
208                         (area     ($ "area" suffix))
209                         (phi      ($ "phi" suffix)))
210                     (check-names ion new-env g gamma delta eps
211                                  minf hinf sinf taum tauh taus 
212                                  alpham betam alphah betah alphas betas
213                                  density total area phi)
214                     (let ((seed-val  (lookup-field 'rng-seed alst))
215                           (gamma-val (eval-const new-env (lookup-field 'gamma alst)))
216                           (delta-val (eval-const new-env (lookup-field 'delta alst 0)))
217                           (eps-val   (eval-const new-env (lookup-field 'eps alst 0)))
218                           (density-val (eval-const new-env (lookup-field 'density alst)))
219                           (area-val  (eval-const new-env (lookup-field 'area alst)))
220                           (phi-val   (eval-const new-env (lookup-field 'phi alst))))
221                       (if (positive? delta-val) (check-decls ion '(hinf tauh) alst))
222                       (if (positive? eps-val) (check-decls ion '(sinf taus) alst))
223                       (if (not (and (integer? gamma-val) (positive? gamma-val)))
224                           (ode:error 'ode:hhsp-transformer 
225                                      "gamma value in ionic conductance declaration " ion
226                                      " must be a positive integer"))
227                       (if (not (and (integer? delta-val) (or (positive? delta-val) (zero? delta-val))))
228                           (ode:error 'ode:hhsp-transformer 
229                                      "delta value in ionic conductance declaration " ion
230                                      " must be a positive integer or zero"))
231                       (if (not (and (integer? eps-val) (or (positive? eps-val) (zero? eps-val))))
232                           (ode:error 'ode:hhsp-transformer 
233                                      "eps value in ionic conductance declaration " ion
234                                      " must be a positive integer or zero"))
235                       (let* ((states        (make-state-names ion gamma-val delta-val eps-val))
236                              (open-state    (last states)))
237                         (apply check-names (cons ion (cons new-env states)))
238                         (let ((total-val   (* density-val area-val))
239                               (nstates     (length states))
240                               (g-rhs      `(/ (* ,phi ,open-state) ,area))
241                               (alpham-rhs  `(/ ,minf ,taum))
242                               (betam-rhs   `(/ (- 1 ,minf) ,taum))
243                               (alphah-rhs  (and (positive? delta-val) `(/ ,hinf ,tauh)))
244                               (betah-rhs   (and (positive? delta-val) `(/ (- 1 ,hinf) ,tauh)))
245                               (alphas-rhs  (and (positive? eps-val) `(/ ,sinf ,taus)))
246                               (betas-rhs   (and (positive? eps-val) `(/ (- 1 ,sinf) ,taus)))
247                               (minf-rhs    (rrhs suffix (lookup-field 'minf alst)))
248                               (hinf-rhs    (and (positive? delta-val) (rrhs suffix (lookup-field  'hinf alst))))
249                               (sinf-rhs    (and (positive? eps-val) (rrhs suffix (lookup-field  'sinf alst))))
250                               (taum-rhs    (rrhs suffix (lookup-field 'taum alst)))
251                               (tauh-rhs    (and (positive? delta-val) (rrhs suffix (lookup-field  'tauh alst))))
252                               (taus-rhs    (and (positive? eps-val) (rrhs suffix (lookup-field  'taus alst))))
253                               (initial-m   (eval-const new-env (lookup-field 'initial-m alst)))
254                               (initial-h   (and (positive? delta-val) 
255                                                 (eval-const new-env (lookup-field 'initial-h alst))))
256                               (initial-s   (and (positive? eps-val) 
257                                                 (eval-const new-env (lookup-field 'initial-s alst)))))
258                           (env-extend! gamma   '(const)  gamma-val)
259                           (env-extend! delta   '(const)  delta-val)
260                           (env-extend! eps     '(const)  eps-val)
261                           (env-extend! total   '(const)  total-val)
262                           (env-extend! density '(const)  density-val)
263                           (env-extend! area    '(const)  area-val)
264                           (env-extend! phi     '(const)  phi-val)
265                           (env-extend! g       '(asgn)   'none g-rhs)
266                           (env-extend! alpham  '(asgn)   'none alpham-rhs)
267                           (env-extend! betam   '(asgn)   'none betam-rhs)
268                           (env-extend! minf    '(asgn)   'none minf-rhs)
269                           (env-extend! taum    '(asgn)   'none taum-rhs)
270                           (if (positive? delta-val) 
271                               (begin
272                                 (env-extend! alphah '(asgn) 'none alphah-rhs)
273                                 (env-extend! betah  '(asgn) 'none betah-rhs)
274                                 (env-extend! hinf   '(asgn) 'none hinf-rhs)
275                                 (env-extend! tauh   '(asgn) 'none tauh-rhs)))
276                           (if (positive? eps-val) 
277                               (begin
278                                 (env-extend! alphas '(asgn) 'none alphas-rhs)
279                                 (env-extend! betas  '(asgn) 'none betas-rhs)
280                                 (env-extend! sinf   '(asgn) 'none sinf-rhs)
281                                 (env-extend! taus   '(asgn) 'none taus-rhs)))
282                           (let ((open-pr (* (expt initial-m gamma-val)
283                                             (or (and initial-h (expt initial-h delta-val)) 1.0) 
284                                             (or (and initial-s (expt initial-s eps-val)) 1.0))))
285                             (env-extend! open-state '(state) (eval-const new-env `(* ,open-pr ,total-val)))
286                             (for-each (lambda (x) 
287                                         (let ((ival  `(/ (* ,(- 1 open-pr) ,total-val) (- ,nstates 1))))
288                                           (env-extend! x  '(state) (eval-const new-env ival))))
289                                       (take states (- nstates 1)))
290                           (let ((rates (make-rates ion gamma-val delta-val eps-val)))
291                             (if debug (print "ode-hhsm:transformer: rates = " rates))
292                             (env-extend! (gensym 'rate) '(rate) 
293                                          (if seed-val (cons `(seed ,seed-val) rates) rates))))))))))
294                           
295                (else (void)))))
296     (pr-transformer new-env)))))
297
Note: See TracBrowser for help on using the repository browser.