User Guide
sympy_conductance.cpp
Go to the documentation of this file.
1 /*
2  * Copyright 2023 Blue Brain Project, EPFL.
3  * See the top-level LICENSE file for details.
4  *
5  * SPDX-License-Identifier: Apache-2.0
6  */
7 
8 #include <catch2/catch_test_macros.hpp>
9 
10 #include "ast/program.hpp"
11 #include "parser/nmodl_driver.hpp"
12 #include "test/unit/utils/test_utils.hpp"
13 #include "visitors/checkparent_visitor.hpp"
14 #include "visitors/constant_folder_visitor.hpp"
15 #include "visitors/inline_visitor.hpp"
16 #include "visitors/local_var_rename_visitor.hpp"
17 #include "visitors/sympy_conductance_visitor.hpp"
18 #include "visitors/symtab_visitor.hpp"
19 #include "visitors/visitor_utils.hpp"
20 
21 using namespace nmodl;
22 using namespace visitor;
23 using namespace test;
24 using namespace test_utils;
25 
26 using ast::AstNodeType;
27 using nmodl::parser::NmodlDriver;
28 
29 
30 //=============================================================================
31 // SympyConductance visitor tests
32 //=============================================================================
33 
34 std::string run_sympy_conductance_visitor(const std::string& text) {
35  // construct AST from text
36  NmodlDriver driver;
37  const auto& ast = driver.parse_string(text);
38 
39  // construct symbol table from AST
40  SymtabVisitor(false).visit_program(*ast);
41 
42  // run constant folding, inlining & local renaming first
43  ConstantFolderVisitor().visit_program(*ast);
44  InlineVisitor().visit_program(*ast);
45  LocalVarRenameVisitor().visit_program(*ast);
46  SymtabVisitor(true).visit_program(*ast);
47 
48  // run SympyConductance on AST
49  SympyConductanceVisitor().visit_program(*ast);
50 
51  // check that, after visitor rearrangement, parents are still up-to-date
52  CheckParentVisitor().check_ast(*ast);
53 
54  // run lookup visitor to extract results from AST
55  // return BREAKPOINT block as JSON string
56  return reindent_text(to_nmodl(collect_nodes(*ast, {AstNodeType::BREAKPOINT_BLOCK}).front()));
57 }
58 
59 std::string breakpoint_to_nmodl(const std::string& text) {
60  // construct AST from text
61  NmodlDriver driver;
62  const auto& ast = driver.parse_string(text);
63 
64  // construct symbol table from AST
65  SymtabVisitor().visit_program(*ast);
66 
67  // run lookup visitor to extract results from AST
68  // return BREAKPOINT block as JSON string
69  return reindent_text(to_nmodl(collect_nodes(*ast, {AstNodeType::BREAKPOINT_BLOCK}).front()));
70 }
71 
72 void run_sympy_conductance_passes(ast::Program& node) {
73  // construct symbol table from AST
74  SymtabVisitor v_symtab;
75  v_symtab.visit_program(node);
76 
77  // run SympySolver on AST several times
78  SympyConductanceVisitor v_sympy1;
79  v_sympy1.visit_program(node);
80  v_symtab.visit_program(node);
81  v_sympy1.visit_program(node);
82  v_symtab.visit_program(node);
83 
84  // also use a second instance of SympySolver
85  SympyConductanceVisitor v_sympy2;
86  v_sympy2.visit_program(node);
87  v_symtab.visit_program(node);
88  v_sympy1.visit_program(node);
89  v_symtab.visit_program(node);
90  v_sympy2.visit_program(node);
91  v_symtab.visit_program(node);
92 }
93 
94 SCENARIO("Addition of CONDUCTANCE using SympyConductance visitor", "[visitor][solver][sympy]") {
95  // First set of test mod files below all based on:
96  // nmodldb/models/db/bluebrain/CortexSimplified/mod/Ca.mod
97  GIVEN("breakpoint block containing VERBATIM statement") {
98  std::string nmodl_text = R"(
99  NEURON {
100  SUFFIX Ca
101  USEION ca READ eca WRITE ica
102  RANGE gCabar, gCa, ica
103  }
104  BREAKPOINT {
105  CONDUCTANCE gCa USEION ca
106  SOLVE states METHOD cnexp
107  VERBATIM
108  double z=0;
109  ENDVERBATIM
110  gCa = gCabar*m*m*h
111  ica = gCa*(v-eca)
112  }
113  )";
114  std::string breakpoint_text = R"(
115  BREAKPOINT {
116  CONDUCTANCE gCa USEION ca
117  SOLVE states METHOD cnexp
118  VERBATIM
119  double z=0;
120  ENDVERBATIM
121  gCa = gCabar*m*m*h
122  ica = gCa*(v-eca)
123  }
124  )";
125  THEN("Do nothing") {
126  auto result = run_sympy_conductance_visitor(nmodl_text);
127  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
128  }
129  }
130  GIVEN("breakpoint block containing IF/ELSE block") {
131  std::string nmodl_text = R"(
132  NEURON {
133  SUFFIX Ca
134  USEION ca READ eca WRITE ica
135  RANGE gCabar, gCa, ica
136  }
137  BREAKPOINT {
138  CONDUCTANCE gCa USEION ca
139  SOLVE states METHOD cnexp
140  IF(gCabar<1){
141  gCa = gCabar*m*m*h
142  ica = gCa*(v-eca)
143  } ELSE {
144  gCa = 0
145  ica = 0
146  }
147  }
148  )";
149  std::string breakpoint_text = R"(
150  BREAKPOINT {
151  CONDUCTANCE gCa USEION ca
152  SOLVE states METHOD cnexp
153  IF(gCabar<1){
154  gCa = gCabar*m*m*h
155  ica = gCa*(v-eca)
156  } ELSE {
157  gCa = 0
158  ica = 0
159  }
160  }
161  )";
162  THEN("Do nothing") {
163  auto result = run_sympy_conductance_visitor(nmodl_text);
164  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
165  }
166  }
167  GIVEN("ion current, existing CONDUCTANCE hint & var") {
168  std::string nmodl_text = R"(
169  NEURON {
170  SUFFIX Ca
171  USEION ca READ eca WRITE ica
172  RANGE gCabar, gCa, ica
173  }
174 
175  UNITS {
176  (S) = (siemens)
177  (mV) = (millivolt)
178  (mA) = (milliamp)
179  }
180 
181  PARAMETER {
182  gCabar = 0.00001 (S/cm2)
183  }
184 
185  ASSIGNED {
186  v (mV)
187  eca (mV)
188  ica (mA/cm2)
189  gCa (S/cm2)
190  mInf
191  mTau
192  mAlpha
193  mBeta
194  hInf
195  hTau
196  hAlpha
197  hBeta
198  }
199 
200  STATE {
201  m
202  h
203  }
204 
205  BREAKPOINT {
206  CONDUCTANCE gCa USEION ca
207  SOLVE states METHOD cnexp
208  gCa = gCabar*m*m*h
209  ica = gCa*(v-eca)
210  }
211 
212  DERIVATIVE states {
213  m' = (mInf-m)/mTau
214  h' = (hInf-h)/hTau
215  }
216 
217  INITIAL{
218  m = mInf
219  h = hInf
220  }
221  )";
222  std::string breakpoint_text = R"(
223  BREAKPOINT {
224  CONDUCTANCE gCa USEION ca
225  SOLVE states METHOD cnexp
226  gCa = gCabar*m*m*h
227  ica = gCa*(v-eca)
228  }
229  )";
230  THEN("Do nothing") {
231  auto result = run_sympy_conductance_visitor(nmodl_text);
232  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
233  }
234  }
235  GIVEN("ion current, no CONDUCTANCE hint, existing var") {
236  std::string nmodl_text = R"(
237  NEURON {
238  SUFFIX Ca
239  USEION ca READ eca WRITE ica
240  RANGE gCabar, gCa, ica
241  }
242 
243  UNITS {
244  (S) = (siemens)
245  (mV) = (millivolt)
246  (mA) = (milliamp)
247  }
248 
249  PARAMETER {
250  gCabar = 0.00001 (S/cm2)
251  }
252 
253  ASSIGNED {
254  v (mV)
255  eca (mV)
256  ica (mA/cm2)
257  gCa (S/cm2)
258  mInf
259  mTau
260  mAlpha
261  mBeta
262  hInf
263  hTau
264  hAlpha
265  hBeta
266  }
267 
268  STATE {
269  m
270  h
271  }
272 
273  BREAKPOINT {
274  SOLVE states METHOD cnexp
275  gCa = gCabar*m*m*h
276  ica = gCa*(v-eca)
277  }
278 
279  DERIVATIVE states {
280  m' = (mInf-m)/mTau
281  h' = (hInf-h)/hTau
282  }
283 
284  INITIAL{
285  m = mInf
286  h = hInf
287  }
288  )";
289  std::string breakpoint_text = R"(
290  BREAKPOINT {
291  CONDUCTANCE gCa USEION ca
292  SOLVE states METHOD cnexp
293  gCa = gCabar*m*m*h
294  ica = gCa*(v-eca)
295  }
296  )";
297  THEN("Add CONDUCTANCE hint using existing var") {
298  auto result = run_sympy_conductance_visitor(nmodl_text);
299  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
300  }
301  }
302  GIVEN("ion current, no CONDUCTANCE hint, no existing var") {
303  std::string nmodl_text = R"(
304  NEURON {
305  SUFFIX Ca
306  USEION ca READ eca WRITE ica
307  RANGE gCabar, ica
308  }
309 
310  UNITS {
311  (S) = (siemens)
312  (mV) = (millivolt)
313  (mA) = (milliamp)
314  }
315 
316  PARAMETER {
317  gCabar = 0.00001 (S/cm2)
318  }
319 
320  ASSIGNED {
321  v (mV)
322  eca (mV)
323  ica (mA/cm2)
324  mInf
325  mTau
326  mAlpha
327  mBeta
328  hInf
329  hTau
330  hAlpha
331  hBeta
332  }
333 
334  STATE {
335  m
336  h
337  }
338 
339  BREAKPOINT {
340  SOLVE states METHOD cnexp
341  ica = (gCabar*m*m*h)*(v-eca)
342  }
343 
344  DERIVATIVE states {
345  m' = (mInf-m)/mTau
346  h' = (hInf-h)/hTau
347  }
348 
349  INITIAL{
350  m = mInf
351  h = hInf
352  }
353  )";
354  std::string breakpoint_text = R"(
355  BREAKPOINT {
356  LOCAL g_ca_0
357  CONDUCTANCE g_ca_0 USEION ca
358  g_ca_0 = gCabar*h*pow(m, 2)
359  SOLVE states METHOD cnexp
360  ica = (gCabar*m*m*h)*(v-eca)
361  }
362  )";
363  THEN("Add CONDUCTANCE hint with new local var") {
364  auto result = run_sympy_conductance_visitor(nmodl_text);
365  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
366  }
367  }
368  GIVEN("2 ion currents, 1 CONDUCTANCE hint, 1 existing var") {
369  std::string nmodl_text = R"(
370  NEURON {
371  SUFFIX Ca
372  USEION ca READ eca WRITE ica
373  USEION na READ ena WRITE ina
374  RANGE gCabar, gNabar, ica, ina
375  }
376 
377  UNITS {
378  (S) = (siemens)
379  (mV) = (millivolt)
380  (mA) = (milliamp)
381  }
382 
383  PARAMETER {
384  gCabar = 0.00001 (S/cm2)
385  gNabar = 0.00005 (S/cm2)
386  }
387 
388  ASSIGNED {
389  v (mV)
390  eca (mV)
391  ica (mA/cm2)
392  ina (mA/cm2)
393  gCa (S/cm2)
394  mInf
395  mTau
396  mAlpha
397  mBeta
398  hInf
399  hTau
400  hAlpha
401  hBeta
402  }
403 
404  STATE {
405  m
406  h
407  }
408 
409  BREAKPOINT {
410  CONDUCTANCE gCa USEION ca
411  SOLVE states METHOD cnexp
412  gCa = gCabar*m*m*h
413  ica = gCa*(v-eca)
414  ina = (gNabar*m*h)*(v-eca)
415  }
416 
417  DERIVATIVE states {
418  m' = (mInf-m)/mTau
419  h' = (hInf-h)/hTau
420  }
421 
422  INITIAL{
423  m = mInf
424  h = hInf
425  }
426  )";
427  std::string breakpoint_text = R"(
428  BREAKPOINT {
429  LOCAL g_na_0
430  CONDUCTANCE g_na_0 USEION na
431  g_na_0 = gNabar*h*m
432  CONDUCTANCE gCa USEION ca
433  SOLVE states METHOD cnexp
434  gCa = gCabar*m*m*h
435  ica = gCa*(v-eca)
436  ina = (gNabar*m*h)*(v-eca)
437  }
438  )";
439  THEN("Add 1 CONDUCTANCE hint with new local var") {
440  auto result = run_sympy_conductance_visitor(nmodl_text);
441  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
442  }
443  }
444  GIVEN("2 ion currents, no CONDUCTANCE hints, 1 existing var") {
445  std::string nmodl_text = R"(
446  NEURON {
447  SUFFIX Ca
448  USEION ca READ eca WRITE ica
449  USEION na READ ena WRITE ina
450  RANGE gCabar, gNabar, ica, ina
451  }
452 
453  UNITS {
454  (S) = (siemens)
455  (mV) = (millivolt)
456  (mA) = (milliamp)
457  }
458 
459  PARAMETER {
460  gCabar = 0.00001 (S/cm2)
461  gNabar = 0.00005 (S/cm2)
462  }
463 
464  ASSIGNED {
465  v (mV)
466  eca (mV)
467  ica (mA/cm2)
468  ina (mA/cm2)
469  gCa (S/cm2)
470  mInf
471  mTau
472  mAlpha
473  mBeta
474  hInf
475  hTau
476  hAlpha
477  hBeta
478  }
479 
480  STATE {
481  m
482  h
483  }
484 
485  BREAKPOINT {
486  SOLVE states METHOD cnexp
487  gCa = gCabar*m*m*h
488  ica = gCa*(v-eca)
489  ina = (gNabar*m*h)*(v-eca)
490  }
491 
492  DERIVATIVE states {
493  m' = (mInf-m)/mTau
494  h' = (hInf-h)/hTau
495  }
496 
497  INITIAL{
498  m = mInf
499  h = hInf
500  }
501  )";
502  std::string breakpoint_text = R"(
503  BREAKPOINT {
504  LOCAL g_na_0
505  CONDUCTANCE g_na_0 USEION na
506  CONDUCTANCE gCa USEION ca
507  g_na_0 = gNabar*h*m
508  SOLVE states METHOD cnexp
509  gCa = gCabar*m*m*h
510  ica = gCa*(v-eca)
511  ina = (gNabar*m*h)*(v-eca)
512  }
513  )";
514  THEN("Add 2 CONDUCTANCE hints, 1 with existing var, 1 with new local var") {
515  auto result = run_sympy_conductance_visitor(nmodl_text);
516  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
517  }
518  }
519  GIVEN("2 ion currents, no CONDUCTANCE hints, no existing vars") {
520  std::string nmodl_text = R"(
521  NEURON {
522  SUFFIX Ca
523  USEION ca READ eca WRITE ica
524  USEION na READ ena WRITE ina
525  RANGE gCabar, gNabar, ica, ina
526  }
527 
528  UNITS {
529  (S) = (siemens)
530  (mV) = (millivolt)
531  (mA) = (milliamp)
532  }
533 
534  PARAMETER {
535  gCabar = 0.00001 (S/cm2)
536  gNabar = 0.00005 (S/cm2)
537  }
538 
539  ASSIGNED {
540  v (mV)
541  eca (mV)
542  ica (mA/cm2)
543  ina (mA/cm2)
544  gCa (S/cm2)
545  mInf
546  mTau
547  mAlpha
548  mBeta
549  hInf
550  hTau
551  hAlpha
552  hBeta
553  }
554 
555  STATE {
556  m
557  h
558  }
559 
560  BREAKPOINT {
561  SOLVE states METHOD cnexp
562  ica = (gCabar*m*m*h)*(v-eca)
563  ina = (gNabar*m*h)*(v-eca)
564  }
565 
566  DERIVATIVE states {
567  m' = (mInf-m)/mTau
568  h' = (hInf-h)/hTau
569  }
570 
571  INITIAL{
572  m = mInf
573  h = hInf
574  }
575  )";
576  std::string breakpoint_text = R"(
577  BREAKPOINT {
578  LOCAL g_ca_0, g_na_0
579  CONDUCTANCE g_na_0 USEION na
580  CONDUCTANCE g_ca_0 USEION ca
581  g_ca_0 = gCabar*h*pow(m, 2)
582  g_na_0 = gNabar*h*m
583  SOLVE states METHOD cnexp
584  ica = (gCabar*m*m*h)*(v-eca)
585  ina = (gNabar*m*h)*(v-eca)
586  }
587  )";
588  THEN("Add 2 CONDUCTANCE hints with 2 new local vars") {
589  auto result = run_sympy_conductance_visitor(nmodl_text);
590  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
591  }
592  }
593  GIVEN("1 ion current, 1 nonspecific current, no CONDUCTANCE hints, no existing vars") {
594  std::string nmodl_text = R"(
595  NEURON {
596  SUFFIX Ca
597  USEION ca READ eca WRITE ica
598  NONSPECIFIC_CURRENT ihcn
599  RANGE gCabar, ica
600  }
601 
602  UNITS {
603  (S) = (siemens)
604  (mV) = (millivolt)
605  (mA) = (milliamp)
606  }
607 
608  PARAMETER {
609  gCabar = 0.00001 (S/cm2)
610  }
611 
612  ASSIGNED {
613  v (mV)
614  eca (mV)
615  ica (mA/cm2)
616  ihcn (mA/cm2)
617  gCa (S/cm2)
618  mInf
619  mTau
620  mAlpha
621  mBeta
622  hInf
623  hTau
624  hAlpha
625  hBeta
626  }
627 
628  STATE {
629  m
630  h
631  }
632 
633  BREAKPOINT {
634  SOLVE states METHOD cnexp
635  ica = (gCabar*m*m*h)*(v-eca)
636  ihcn = (0.1235*m*h)*(v-eca)
637  }
638 
639  DERIVATIVE states {
640  m' = (mInf-m)/mTau
641  h' = (hInf-h)/hTau
642  }
643 
644  INITIAL{
645  m = mInf
646  h = hInf
647  }
648  )";
649  std::string breakpoint_text = R"(
650  BREAKPOINT {
651  LOCAL g_ca_0, g__0
652  CONDUCTANCE g__0
653  CONDUCTANCE g_ca_0 USEION ca
654  g_ca_0 = gCabar*h*pow(m, 2)
655  g__0 = 0.1235*h*m
656  SOLVE states METHOD cnexp
657  ica = (gCabar*m*m*h)*(v-eca)
658  ihcn = (0.1235*m*h)*(v-eca)
659  }
660  )";
661  THEN("Add 2 CONDUCTANCE hints with 2 new local vars") {
662  auto result = run_sympy_conductance_visitor(nmodl_text);
663  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
664  }
665  }
666  GIVEN("1 ion current, 1 nonspecific current, no CONDUCTANCE hints, 1 existing var") {
667  std::string nmodl_text = R"(
668  NEURON {
669  SUFFIX Ca
670  USEION ca READ eca WRITE ica
671  NONSPECIFIC_CURRENT ihcn
672  RANGE gCabar, ica, gihcn
673  }
674 
675  UNITS {
676  (S) = (siemens)
677  (mV) = (millivolt)
678  (mA) = (milliamp)
679  }
680 
681  PARAMETER {
682  gCabar = 0.00001 (S/cm2)
683  }
684 
685  ASSIGNED {
686  v (mV)
687  eca (mV)
688  ica (mA/cm2)
689  ihcn (mA/cm2)
690  gCa (S/cm2)
691  gihcn (S/cm2)
692  mInf
693  mTau
694  mAlpha
695  mBeta
696  hInf
697  hTau
698  hAlpha
699  hBeta
700  }
701 
702  STATE {
703  m
704  h
705  }
706 
707  BREAKPOINT {
708  SOLVE states METHOD cnexp
709  gihcn = 0.1235*m*h
710  ica = (gCabar*m*m*h)*(v-eca)
711  ihcn = gihcn*(v-eca)
712  }
713 
714  DERIVATIVE states {
715  m' = (mInf-m)/mTau
716  h' = (hInf-h)/hTau
717  }
718 
719  INITIAL{
720  m = mInf
721  h = hInf
722  }
723  )";
724  std::string breakpoint_text = R"(
725  BREAKPOINT {
726  LOCAL g_ca_0
727  CONDUCTANCE gihcn
728  CONDUCTANCE g_ca_0 USEION ca
729  g_ca_0 = gCabar*h*pow(m, 2)
730  SOLVE states METHOD cnexp
731  gihcn = 0.1235*m*h
732  ica = (gCabar*m*m*h)*(v-eca)
733  ihcn = gihcn*(v-eca)
734  }
735  )";
736  THEN("Add 2 CONDUCTANCE hints, 1 using existing var, 1 with new local var") {
737  auto result = run_sympy_conductance_visitor(nmodl_text);
738  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
739  }
740  }
741  // based on bluebrain/CortextPlastic/mod/ProbAMPANMDA.mod
742  GIVEN(
743  "2 ion currents, 1 nonspecific current, no CONDUCTANCE hints, indirect relation between "
744  "eqns") {
745  std::string nmodl_text = R"(
746  NEURON {
747  THREADSAFE
748  POINT_PROCESS ProbAMPANMDA
749  RANGE tau_r_AMPA, tau_d_AMPA, tau_r_NMDA, tau_d_NMDA
750  RANGE Use, u, Dep, Fac, u0, mg, NMDA_ratio
751  RANGE i, i_AMPA, i_NMDA, g_AMPA, g_NMDA, g, e
752  NONSPECIFIC_CURRENT i, i_AMPA,i_NMDA
753  POINTER rng
754  RANGE synapseID, verboseLevel
755  }
756  PARAMETER {
757  tau_r_AMPA = 0.2 (ms) : dual-exponential conductance profile
758  tau_d_AMPA = 1.7 (ms) : IMPORTANT: tau_r < tau_d
759  tau_r_NMDA = 0.29 (ms) : dual-exponential conductance profile
760  tau_d_NMDA = 43 (ms) : IMPORTANT: tau_r < tau_d
761  Use = 1.0 (1) : Utilization of synaptic efficacy (just initial values! Use, Dep and Fac are overwritten by BlueBuilder assigned values)
762  Dep = 100 (ms) : relaxation time constant from depression
763  Fac = 10 (ms) : relaxation time constant from facilitation
764  e = 0 (mV) : AMPA and NMDA reversal potential
765  mg = 1 (mM) : initial concentration of mg2+
766  mggate
767  gmax = .001 (uS) : weight conversion factor (from nS to uS)
768  u0 = 0 :initial value of u, which is the running value of Use
769  NMDA_ratio = 0.71 (1) : The ratio of NMDA to AMPA
770  synapseID = 0
771  verboseLevel = 0
772  }
773  ASSIGNED {
774  v (mV)
775  i (nA)
776  i_AMPA (nA)
777  i_NMDA (nA)
778  g_AMPA (uS)
779  g_NMDA (uS)
780  g (uS)
781  factor_AMPA
782  factor_NMDA
783  rng
784  }
785  STATE {
786  A_AMPA : AMPA state variable to construct the dual-exponential profile - decays with conductance tau_r_AMPA
787  B_AMPA : AMPA state variable to construct the dual-exponential profile - decays with conductance tau_d_AMPA
788  A_NMDA : NMDA state variable to construct the dual-exponential profile - decays with conductance tau_r_NMDA
789  B_NMDA : NMDA state variable to construct the dual-exponential profile - decays with conductance tau_d_NMDA
790  }
791  BREAKPOINT {
792  SOLVE state METHOD cnexp
793  mggate = 1.2
794  g_AMPA = gmax*(B_AMPA-A_AMPA) :compute time varying conductance as the difference of state variables B_AMPA and A_AMPA
795  g_NMDA = gmax*(B_NMDA-A_NMDA) * mggate :compute time varying conductance as the difference of state variables B_NMDA and A_NMDA and mggate kinetics
796  g = g_AMPA + g_NMDA
797  i_AMPA = g_AMPA*(v-e) :compute the AMPA driving force based on the time varying conductance, membrane potential, and AMPA reversal
798  i_NMDA = g_NMDA*(v-e) :compute the NMDA driving force based on the time varying conductance, membrane potential, and NMDA reversal
799  i = i_AMPA + i_NMDA
800  }
801  DERIVATIVE state{
802  A_AMPA' = -A_AMPA/tau_r_AMPA
803  B_AMPA' = -B_AMPA/tau_d_AMPA
804  A_NMDA' = -A_NMDA/tau_r_NMDA
805  B_NMDA' = -B_NMDA/tau_d_NMDA
806  }
807  )";
808  std::string breakpoint_text = R"(
809  BREAKPOINT {
810  CONDUCTANCE g
811  CONDUCTANCE g_NMDA
812  CONDUCTANCE g_AMPA
813  SOLVE state METHOD cnexp
814  mggate = 1.2
815  g_AMPA = gmax*(B_AMPA-A_AMPA) :compute time varying conductance as the difference of state variables B_AMPA and A_AMPA
816  g_NMDA = gmax*(B_NMDA-A_NMDA) * mggate :compute time varying conductance as the difference of state variables B_NMDA and A_NMDA and mggate kinetics
817  g = g_AMPA + g_NMDA
818  i_AMPA = g_AMPA*(v-e) :compute the AMPA driving force based on the time varying conductance, membrane potential, and AMPA reversal
819  i_NMDA = g_NMDA*(v-e) :compute the NMDA driving force based on the time varying conductance, membrane potential, and NMDA reversal
820  i = i_AMPA + i_NMDA
821  }
822  )";
823  THEN("Add 3 CONDUCTANCE hints, using existing vars") {
824  auto result = run_sympy_conductance_visitor(nmodl_text);
825  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
826  }
827  }
828  // based on neurodamus/bbp/lib/modlib/GluSynapse.mod
829  GIVEN("1 nonspecific current, no CONDUCTANCE hints, many eqs & a function involved") {
830  std::string nmodl_text = R"(
831  NEURON {
832  GLOBAL tau_r_AMPA, E_AMPA
833  RANGE tau_d_AMPA, gmax_AMPA
834  RANGE g_AMPA
835  GLOBAL tau_r_NMDA, tau_d_NMDA, E_NMDA
836  RANGE g_NMDA
837  RANGE Use, Dep, Fac, Nrrp, u
838  RANGE tsyn, unoccupied, occupied
839  RANGE ica_NMDA
840  RANGE volume_CR
841  GLOBAL ljp_VDCC, vhm_VDCC, km_VDCC, mtau_VDCC, vhh_VDCC, kh_VDCC, htau_VDCC
842  RANGE gca_bar_VDCC, ica_VDCC
843  GLOBAL gamma_ca_CR, tau_ca_CR, min_ca_CR, cao_CR
844  GLOBAL tau_GB, gamma_d_GB, gamma_p_GB, rho_star_GB, tau_Use_GB, tau_effca_GB
845  RANGE theta_d_GB, theta_p_GB
846  RANGE rho0_GB
847  RANGE enable_GB, depress_GB, potentiate_GB
848  RANGE Use_d_GB, Use_p_GB
849  GLOBAL p_gen_RW, p_elim0_RW, p_elim1_RW
850  RANGE enable_RW, synstate_RW
851  GLOBAL mg, scale_mg, slope_mg
852  RANGE vsyn, NMDA_ratio, synapseID, selected_for_report, verbose
853  NONSPECIFIC_CURRENT i
854  }
855  UNITS {
856  (nA) = (nanoamp)
857  (mV) = (millivolt)
858  (uS) = (microsiemens)
859  (nS) = (nanosiemens)
860  (pS) = (picosiemens)
861  (umho) = (micromho)
862  (um) = (micrometers)
863  (mM) = (milli/liter)
864  (uM) = (micro/liter)
865  FARADAY = (faraday) (coulomb)
866  PI = (pi) (1)
867  R = (k-mole) (joule/degC)
868  }
869  ASSIGNED {
870  g_AMPA (uS)
871  g_NMDA (uS)
872  ica_NMDA (nA)
873  ica_VDCC (nA)
874  depress_GB (1)
875  potentiate_GB (1)
876  v (mV)
877  vsyn (mV)
878  i (nA)
879  }
880  FUNCTION nernst(ci(mM), co(mM), z) (mV) {
881  nernst = (1000) * R * (celsius + 273.15) / (z*FARADAY) * log(co/ci)
882  }
883  BREAKPOINT {
884  LOCAL Eca_syn, mggate, i_AMPA, gmax_NMDA, i_NMDA, Pf_NMDA, gca_bar_abs_VDCC, gca_VDCC
885  g_AMPA = (1e-3)*gmax_AMPA*(B_AMPA-A_AMPA)
886  i_AMPA = g_AMPA*(v-E_AMPA)
887  gmax_NMDA = gmax_AMPA*NMDA_ratio
888  mggate = 1 / (1 + exp(slope_mg * -(v)) * (mg / scale_mg))
889  g_NMDA = (1e-3)*gmax_NMDA*mggate*(B_NMDA-A_NMDA)
890  i_NMDA = g_NMDA*(v-E_NMDA)
891  Pf_NMDA = (4*cao_CR) / (4*cao_CR + (1/1.38) * 120 (mM)) * 0.6
892  ica_NMDA = Pf_NMDA*g_NMDA*(v-40.0)
893  gca_bar_abs_VDCC = gca_bar_VDCC * 4(um2)*PI*(3(1/um3)/4*volume_CR*1/PI)^(2/3)
894  gca_VDCC = (1e-3) * gca_bar_abs_VDCC * m_VDCC * m_VDCC * h_VDCC
895  Eca_syn = FARADAY*nernst(cai_CR, cao_CR, 2)
896  ica_VDCC = gca_VDCC*(v-Eca_syn)
897  vsyn = v
898  i = i_AMPA + i_NMDA + ica_VDCC
899  }
900  )";
901  std::string breakpoint_text = R"(
902  BREAKPOINT {
903  LOCAL Eca_syn, mggate, i_AMPA, gmax_NMDA, i_NMDA, Pf_NMDA, gca_bar_abs_VDCC, gca_VDCC, nernst_in_0, g__0
904  CONDUCTANCE g__0
905  g__0 = (0.001*gmax_NMDA*mg*scale_mg*slope_mg*(A_NMDA-B_NMDA)*(E_NMDA-v)*exp(slope_mg*v)-0.001*gmax_NMDA*scale_mg*(A_NMDA-B_NMDA)*(mg+scale_mg*exp(slope_mg*v))*exp(slope_mg*v)+(g_AMPA+gca_VDCC)*pow(mg+scale_mg*exp(slope_mg*v), 2))/pow(mg+scale_mg*exp(slope_mg*v), 2)
906  g_AMPA = 1e-3*gmax_AMPA*(B_AMPA-A_AMPA)
907  i_AMPA = g_AMPA*(v-E_AMPA)
908  gmax_NMDA = gmax_AMPA*NMDA_ratio
909  mggate = 1/(1+exp(slope_mg*-v)*(mg/scale_mg))
910  g_NMDA = 1e-3*gmax_NMDA*mggate*(B_NMDA-A_NMDA)
911  i_NMDA = g_NMDA*(v-E_NMDA)
912  Pf_NMDA = (4*cao_CR)/(4*cao_CR+0.7246376811594204*120(mM))*0.6
913  ica_NMDA = Pf_NMDA*g_NMDA*(v-40.0)
914  gca_bar_abs_VDCC = gca_bar_VDCC*4(um2)*PI*(3(1/um3)/4*volume_CR*1/PI)^0.6666666666666666
915  gca_VDCC = 1e-3*gca_bar_abs_VDCC*m_VDCC*m_VDCC*h_VDCC
916  {
917  LOCAL ci_in_0, co_in_0, z_in_0
918  ci_in_0 = cai_CR
919  co_in_0 = cao_CR
920  z_in_0 = 2
921  nernst_in_0 = 1000*R*(celsius+273.15)/(z_in_0*FARADAY)*log(co_in_0/ci_in_0)
922  }
923  Eca_syn = FARADAY*nernst_in_0
924  ica_VDCC = gca_VDCC*(v-Eca_syn)
925  vsyn = v
926  i = i_AMPA+i_NMDA+ica_VDCC
927  }
928  )";
929  THEN("Add 1 CONDUCTANCE hint using new var") {
930  auto result = run_sympy_conductance_visitor(nmodl_text);
931  REQUIRE(result == breakpoint_to_nmodl(breakpoint_text));
932  }
933  }
934 }
test_utils.hpp
nmodl::parser::NmodlDriver
Class that binds all pieces together for parsing nmodl file.
Definition: nmodl_driver.hpp:67
nmodl::to_nmodl
std::string to_nmodl(const ast::Ast &node, const std::set< ast::AstNodeType > &exclude_types)
Given AST node, return the NMODL string representation.
Definition: visitor_utils.cpp:234
nmodl::test_utils::reindent_text
std::string reindent_text(const std::string &text, int indent_level)
Reindent nmodl text for text-to-text comparison.
Definition: test_utils.cpp:53
nmodl
encapsulates code generation backend implementations
Definition: ast_common.hpp:26
run_sympy_conductance_visitor
std::string run_sympy_conductance_visitor(const std::string &text)
Definition: sympy_conductance.cpp:34
constant_folder_visitor.hpp
Perform constant folding of integer/float/double expressions.
nmodl::ast::AstNodeType
AstNodeType
Enum type for every AST node type.
Definition: ast_decl.hpp:166
visitor_utils.hpp
Utility functions for visitors implementation.
sympy_conductance_visitor.hpp
Visitor for generating CONDUCTANCE statements for ions
program.hpp
Auto generated AST classes declaration.
driver
nmodl::parser::UnitDriver driver
Definition: parser.cpp:28
run_sympy_conductance_passes
void run_sympy_conductance_passes(ast::Program &node)
Definition: sympy_conductance.cpp:72
local_var_rename_visitor.hpp
Visitor to rename local variables conflicting with global scope
nmodl::parser::UnitDriver::parse_string
bool parse_string(const std::string &input)
parser Units provided as string (used for testing)
Definition: unit_driver.cpp:40
nmodl::collect_nodes
std::vector< std::shared_ptr< const ast::Ast > > collect_nodes(const ast::Ast &node, const std::vector< ast::AstNodeType > &types)
traverse node recursively and collect nodes of given types
Definition: visitor_utils.cpp:206
checkparent_visitor.hpp
Visitor for checking parents of ast nodes
inline_visitor.hpp
Visitor to inline local procedure and function calls
nmodl_driver.hpp
nmodl::ast::Program
Represents top level AST node for whole NMODL input.
Definition: program.hpp:39
SCENARIO
SCENARIO("Addition of CONDUCTANCE using SympyConductance visitor", "[visitor][solver][sympy]")
Definition: sympy_conductance.cpp:94
symtab_visitor.hpp
THIS FILE IS GENERATED AT BUILD TIME AND SHALL NOT BE EDITED.
breakpoint_to_nmodl
std::string breakpoint_to_nmodl(const std::string &text)
Definition: sympy_conductance.cpp:59