source: project/ode/tags/2.7a/extensions/ode-hhsm.scm @ 7359

Last change on this file since 7359 was 7359, checked in by Ivan Raikov, 12 years ago

Created releases with new author URL.

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