1: //
   2: // Heccer : a compartmental solver that implements efficient Crank-Nicolson
   3: // integration for neuronal models.
   4: //
   5: 
   6: //////////////////////////////////////////////////////////////////////////////
   7: //'
   8: //' Heccer : testbed C implementation
   9: //'
  10: //' Copyright (C) 2006-2008 Hugo Cornelis
  11: //'
  12: //' functional ideas .. Hugo Cornelis, hugo.cornelis@gmail.com
  13: //'
  14: //' coding ............ Hugo Cornelis, hugo.cornelis@gmail.com
  15: //'
  16: //////////////////////////////////////////////////////////////////////////////
  17: 
  18: 
  19: #include <stdlib.h>
  20: 
  21: #include "../../heccer/addressing.h"
  22: #include "../../heccer/compartment.h"
  23: #include "../../heccer/heccer.h"
  24: #include "../../heccer/output.h"
  25: 
  26: 
  27: #define HECCER_TEST_CONSTRUCT \
  28:     memcpy(&pheccer->inter, &inter, sizeof(inter));     \
  29:     pheccer->iStatus = HECCER_STATUS_PHASE_2
  30: #define HECCER_TEST_REPORTING_GRANULARITY 10000
  31: #define HECCER_TEST_STEPS 10000
  32: #define HECCER_TEST_TESTED_THINGS ( HECCER_DUMP_VM_COMPARTMENT_MATRIX \
  33:                                     | HECCER_DUMP_VM_COMPARTMENT_DATA \
  34:                                     | HECCER_DUMP_VM_COMPARTMENT_OPERATIONS \
  35:                                     | HECCER_DUMP_VM_MECHANISM_DATA \
  36:                                     | HECCER_DUMP_VM_MECHANISM_OPERATIONS \
  37:                                     | HECCER_DUMP_VM_SUMMARY \
  38:         )
  39: #define HECCER_TEST_TIME_STEP (1e-6)
  40: 
  41: 
  42: struct Compartment compSoma =
  43: {
  44:     //m administrative overhead
  45: 
  46:     {
  47:         //m type of structure
  48: 
  49:         MATH_TYPE_Compartment,
  50:     },
  51: 
  52:     //m index of parent compartment, -1 for none
  53: 
  54:     -1,
  55: 
  56: /*     //m first mechanism */
  57: 
  58: /*     NULL, */
  59: 
  60: /*     //m number of mechanisms */
  61: 
  62: /*     0, */
  63: 
  64:     //m descriptive values, alphabetical order
  65: 
  66: /*     double dCm; */
  67: 
  68:     4.57537e-11, // unscaled 0.0164,
  69: 
  70: /*     double dEm; */
  71: 
  72:     -0.08,
  73: 
  74: /*     double dInitVm; */
  75: 
  76:     -0.065,
  77: 
  78: /*     double dInject;           */
  79: 
  80:     2e-9,
  81: 
  82: /*     double dRa; */
  83: 
  84:     360502, // unscaled 2.5,
  85: 
  86: /*     double dRm; */
  87: 
  88:     3.58441e+08, // unscaled 1
  89: };
  90: 
  91: 
  92: //v a regular delayed rectifier, taken from the genesis website
  93: 
  94: struct ChannelAct caKdr =
  95: {
  96:     //m administrative overhead
  97: 
  98:     {
  99:         //m type of structure
 100: 
 101:         MATH_TYPE_ChannelAct,
 102:     },
 103: 
 104:     //m first set of descriptive values, alphabetical order
 105: 
 106:     //m initial reversal potential
 107: 
 108:     -0.085,
 109: 
 110:     //m get reversal potential from this intermediary, -1 for none
 111: 
 112:     -1,
 113: 
 114:     //m maximal conductance when all channels are permissive
 115: 
 116:     1.004349542e-06,
 117: 
 118:     //m contributes to this concentration pool, -1 for none, boolean indicator only.
 119: 
 120:     -1,
 121: 
 122:     //m activation description
 123: 
 124:     {
 125:         //m power, for a standard heccer, something between 1 and 4 or so.
 126: 
 127:         4,
 128: 
 129:         //m gate definition
 130: 
 131:         {
 132:             //m initial value, commonly forward over backward steady states
 133: 
 134:             0.3176769506,
 135: 
 136:             //m corresponding index in tables, set to -1 for initialization.
 137: 
 138:             -1,
 139: 
 140:             {
 141:                 //m forward kinetiks, commonly denoted with alpha or non-perm to perm rate
 142: 
 143:                 {
 144:                     //m multiplier
 145: 
 146:                     -550.0,
 147: 
 148:                     //m multiplier membrane dependence, 0.0 for no dependence
 149: 
 150:                     10000,
 151: 
 152:                     //m 2b: multiplier membrane dependence offset, 0.0 for no dependence
 153: 
 154:                     0.0,
 155: 
 156:                     //m choose between nominator or denominator, 1 means nominator, -1
 157:                     //m means denominator
 158: 
 159:                     -1.0,
 160: 
 161:                     //m nominator or denominator offset
 162: 
 163:                     -1.0,
 164: 
 165:                     //m membrane offset
 166: 
 167:                     0.055,
 168: 
 169:                     //m denormalized time constant
 170: 
 171:                     -0.01,
 172:                 },
 173: 
 174:                 //m backward kinetiks, commonly denoted with beta or perm to non-perm rate
 175: 
 176:                 {
 177:                     //m multiplier
 178: 
 179:                     125.0,
 180: 
 181:                     //m multiplier membrane dependence, 0.0 for no dependence
 182: 
 183:                     0.0,
 184: 
 185:                     //m 2b: multiplier membrane dependence offset, 0.0 for no dependence
 186: 
 187:                     0.0,
 188: 
 189:                     //m choose between nominator or denominator, 1 means nominator, -1
 190:                     //m means denominator
 191: 
 192:                     -1.0,
 193: 
 194:                     //m nominator or denominator offset
 195: 
 196:                     0.0,
 197: 
 198:                     //m membrane offset
 199: 
 200:                     0.065,
 201: 
 202:                     //m denormalized time constant
 203: 
 204:                     0.08,
 205:                 },
 206:             },
 207:         },
 208:     },
 209: 
 210: };
 211: 
 212: 
 213: //v a simple fast sodium channel
 214: 
 215: struct ChannelActInact caiNaF =
 216: {
 217:     //m administrative overhead
 218: 
 219:     {
 220:         //m type of structure
 221: 
 222:         MATH_TYPE_ChannelActInact,
 223:     },
 224: 
 225:     //m first set of descriptive values, alphabetical order
 226: 
 227:     //m initial reversal potential
 228: 
 229:     0.045,
 230: 
 231:     //m get reversal potential from this intermediary, -1 for none
 232: 
 233:     -1,
 234: 
 235:     //m maximal conductance when all channels are permissive
 236: 
 237:     0.00001,
 238: 
 239:     //m contributes to this concentration pool, -1 for none, boolean indicator only.
 240: 
 241:     -1,
 242: 
 243:     //m activation description
 244: 
 245:     {
 246:         //m power, for a standard heccer, something between 1 and 4 or so.
 247: 
 248:         3,
 249: 
 250:         //m gate definition
 251: 
 252:         {
 253:             //m initial value, commonly forward over backward steady states
 254: 
 255:             0.01224204035,
 256: 
 257:             //m corresponding index in tables, set to -1 for initialization.
 258: 
 259:             -1,
 260: 
 261:             {
 262:                 //m forward kinetiks, commonly denoted with alpha or non-perm to perm rate
 263: 
 264:                 {
 265:                     //m multiplier
 266: 
 267:                     35.0e3,
 268: 
 269:                     //m multiplier membrane dependence, 0.0 for no dependence
 270: 
 271:                     0.0,
 272: 
 273:                     //m 2b: multiplier membrane dependence offset, 0.0 for no dependence
 274: 
 275:                     0.0,
 276: 
 277:                     //m choose between nominator or denominator, 1 means nominator, -1
 278:                     //m means denominator
 279: 
 280:                     -1.0,
 281: 
 282:                     //m nominator or denominator offset
 283: 
 284:                     0.0,
 285: 
 286:                     //m membrane offset
 287: 
 288:                     5.0e-3,
 289: 
 290:                     //m denormalized time constant
 291: 
 292:                     -10.0e-3,
 293:                 },
 294: 
 295:                 //m backward kinetiks, commonly denoted with beta or perm to non-perm rate
 296: 
 297:                 {
 298:                     //m multiplier
 299: 
 300:                     7.0e3,
 301: 
 302:                     //m multiplier membrane dependence, 0.0 for no dependence
 303: 
 304:                     0.0,
 305: 
 306:                     //m 2b: multiplier membrane dependence offset, 0.0 for no dependence
 307: 
 308:                     0.0,
 309: 
 310:                     //m choose between nominator or denominator, 1 means nominator, -1
 311:                     //m means denominator
 312: 
 313:                     -1.0,
 314: 
 315:                     //m nominator or denominator offset
 316: 
 317:                     0.0,
 318: 
 319:                     //m membrane offset
 320: 
 321:                     65.0e-3,
 322: 
 323:                     //m denormalized time constant
 324: 
 325:                     20.0e-3,
 326:                 },
 327:             },
 328:         },
 329:     },
 330: 
 331:     //m inactivation description
 332: 
 333:     {
 334:         //m power, for a standard heccer, something between 1 and 4 or so.
 335: 
 336:         1,
 337: 
 338:         //m gate definition
 339: 
 340:         {
 341:             //m initial value, commonly forward over backward steady states
 342: 
 343:             0.1930685015,
 344: 
 345:             //m corresponding index in tables, set to -1 for initialization.
 346: 
 347:             -1,
 348: 
 349:             {
 350:                 //m forward kinetiks, commonly denoted with alpha or non-perm to perm rate
 351: 
 352:                 {
 353:                     //m multiplier
 354: 
 355:                     0.225e3,
 356: 
 357:                     //m multiplier membrane dependence, 0.0 for no dependence
 358: 
 359:                     0.0,
 360: 
 361:                     //m 2b: multiplier membrane dependence offset, 0.0 for no dependence
 362: 
 363:                     0.0,
 364: 
 365:                     //m choose between nominator or denominator, 1 means nominator, -1
 366:                     //m means denominator
 367: 
 368:                     -1.0,
 369: 
 370:                     //m nominator or denominator offset
 371: 
 372:                     1.0,
 373: 
 374:                     //m membrane offset
 375: 
 376:                     80.0e-3,
 377: 
 378:                     //m denormalized time constant
 379: 
 380:                     10.0e-3,
 381:                 },
 382: 
 383:                 //m backward kinetiks, commonly denoted with beta or perm to non-perm rate
 384: 
 385:                 {
 386:                     //m multiplier
 387: 
 388:                     7.5e3,
 389: 
 390:                     //m multiplier membrane dependence, 0.0 for no dependence
 391: 
 392:                     0.0,
 393: 
 394:                     //m 2b: multiplier membrane dependence offset, 0.0 for no dependence
 395: 
 396:                     0.0,
 397: 
 398:                     //m choose between nominator or denominator, 1 means nominator, -1
 399:                     //m means denominator
 400: 
 401:                     -1.0,
 402: 
 403:                     //m nominator or denominator offset
 404: 
 405:                     0.0,
 406: 
 407:                     //m membrane offset
 408: 
 409:                     -3.0e-3,
 410: 
 411:                     //m denormalized time constant
 412: 
 413:                     -18.0e-3,
 414:                 },
 415:             },
 416:         },
 417:     },
 418: };
 419: 
 420: 
 421: struct SpikeGenerator sg =
 422: {
 423:     //m administrative overhead
 424: 
 425:     {
 426:         //m type of structure
 427: 
 428:         MATH_TYPE_SpikeGenerator,
 429: 
 430: #ifdef HECCER_SOURCE_NEUROSPACES
 431: 
 432:         //m identification
 433: 
 434:         ADDRESSING_NEUROSPACES_2_HECCER(1000),
 435: 
 436: #endif
 437: 
 438:     },
 439: 
 440:     //m refractory time
 441: 
 442:     0.001,
 443: 
 444:     //m spiking threshold
 445: 
 446:     0.0,
 447: 
 448:     //m reset value, FLT_MAX for none
 449: 
 450:     FLT_MAX,
 451: 
 452:     //m table in event distributor with targets
 453: 
 454:     0,
 455: };
 456: 
 457: 
 458: int piC2m[] =
 459: {
 460:     3,
 461:     -1,
 462: };
 463: 
 464: 
 465: struct MathComponentArray mca =
 466: {
 467:     //m number of math components
 468: 
 469:     3,
 470: 
 471:     //m math component data
 472: 
 473:     NULL,
 474: 
 475:     //m math component index, initialize to NULL
 476: 
 477:     NULL,
 478: 
 479: };
 480: 
 481: 
 482: struct Intermediary inter =
 483: {
 484:     //m compartment array
 485: 
 486:     1,
 487: 
 488:     &compSoma,
 489: 
 490:     //m all other mathematical components
 491: 
 492:     &mca,
 493: 
 494:     //m compartment 2 first mechanism number
 495: 
 496:     piC2m,
 497: };
 498: 
 499: 
 500: #include "../../heccer/eventdistributor.h"
 501: 
 502: struct EventDistributorMatrix pedm[] =
 503: {
 504:     {
 505:         NULL,
 506:         -1,
 507:         NULL,
 508:         -1,
 509:     },
 510:     {
 511:         NULL,
 512:         -1,
 513:         NULL,
 514:         -1,
 515:     },
 516: };
 517: 
 518: struct EventDistributorData edd =
 519: {
 520:     //m number of connections in the matrix
 521: 
 522:     2,
 523: 
 524:     //m last used
 525: 
 526:     1,
 527: 
 528:     //m array of targets
 529: 
 530:     pedm,
 531: };
 532: 
 533: 
 534: struct EventDistributor ed =
 535: {
 536:     //m service specific data
 537: 
 538:     &edd,
 539: 
 540:     //m distribute an event over the targets
 541: 
 542:     EventDistributorSend,
 543: };
 544: 
 545: 
 546: struct OutputGenerator *pogSpike = NULL;
 547: struct OutputGenerator * pogVm = NULL;
 548: 
 549: /* double *pdVm = NULL; */
 550: 
 551: char pcStep[100] = "";
 552: 
 553: #include "main.h"
 554: 
 555: 
 556: int main(int argc, char *argv[])
 557: {
 558:     //- set default result : ok
 559: 
 560:     int iResult = EXIT_SUCCESS;
 561: 
 562:     //- determine intermediary size, and allocate
 563: 
 564:     struct MathComponentInfo *pmciKdr = MathComponentInfoLookup(caKdr.mc.iType);
 565: 
 566:     struct MathComponentInfo *pmciNaF = MathComponentInfoLookup(caiNaF.mc.iType);
 567: 
 568:     struct MathComponentInfo *pmciSpiker = MathComponentInfoLookup(sg.mc.iType);
 569: 
 570:     int iChars = pmciKdr->iChars + pmciNaF->iChars + pmciSpiker->iChars;
 571: 
 572:     void *pmc = calloc(sizeof(char), iChars);
 573: 
 574:     //- prepare the mechanism intermediary
 575: 
 576:     struct ChannelAct *pca = (struct ChannelAct *)pmc;
 577: 
 578:     *pca = caKdr;
 579: 
 580:     struct ChannelActInact *pcai = (struct ChannelActInact *)&((char *)pca)[pmciKdr->iChars];
 581: 
 582:     *pcai = caiNaF;
 583: 
 584:     struct SpikeGenerator *psg = (struct SpikeGenerator *)&((char *)pcai)[pmciNaF->iChars];
 585: 
 586:     *psg = sg;
 587: 
 588:     //- link the intermediary
 589: 
 590:     mca.pmc = pmc;
 591: 
 592:     //- link spiking element to output generator
 593: 
 594:     pogSpike = OutputGeneratorNew("/tmp/output_spike");
 595: 
 596:     pedm[0].pvObject = pogSpike;
 597:     pedm[0].pvFunction = OutputGeneratorTimedStep;
 598: 
 599:     //- create output elements
 600: 
 601:     pogVm = OutputGeneratorNew("/tmp/output_vm");
 602: 
 603: //d prepare output of membrane potential and spikes
 604: 
 605: #define HECCER_TEST_INITIATE \
 606:     double *pdVm = HeccerAddressCompartmentVariable(pheccer, 0, "Vm");  \
 607:     OutputGeneratorAddVariable(pogVm, "Vm", pdVm);                      \
 608:     double *pdSpike = HeccerAddressMechanismVariable(pheccer, 2, "spike"); \
 609:     OutputGeneratorAddVariable(pogSpike, "spike", pdSpike)
 610: 
 611: //d generate output of membrane potential each step
 612: 
 613: #define HECCER_TEST_OUTPUT \
 614:     OutputGeneratorAnnotatedStep(pogVm, sprintf(pcStep, "%i", i) ? pcStep : "sprintf() failed")
 615: 
 616:     //- allocate the heccer, for the event distributor service
 617: 
 618:     //t the event distributor should be added to the heccer instances
 619:     //t during output object construction, such that this is closer to
 620:     //t SSP behaviour, and closer to CBI architecture specs.
 621: 
 622:     //t so that is during HECCER_TEST_INITIATE ?
 623: 
 624:     pheccer = HeccerNew("spiker1", NULL, &ed, NULL);
 625: 
 626:     //- do the simulation
 627: 
 628:     simulate(argc,argv);
 629: 
 630:     //- finish the simulation output
 631: 
 632:     OutputGeneratorFinish(pogVm);
 633: 
 634:     OutputGeneratorFinish(pogSpike);
 635: 
 636:     //- add the simulation output to the program output
 637: 
 638:     WriteOutput("/tmp/output_vm");
 639: 
 640:     WriteOutput("/tmp/output_spike");
 641: 
 642:     //- return result
 643: 
 644:     return(iResult);
 645: }
 646: 
 647: 
 648: #define main(argc,argv) simulate(argc,argv)
 649: 
 650: //t this prototype can give warning and perhaps errors.
 651: 
 652: int main(int argc, char *argv[]);
 653: 
 654: 
 655: #include "main.c"
 656: 
 657: 
 658: