file:/local_home/local_home/hugo/neurospaces_project/heccer/source/c/snapshots/0/mechanism.c (Sun Jul 13 20:35:17 2008
) HOME
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 <limits.h>
20: #include <math.h>
21: #include <stdarg.h>
22: #include <stdio.h>
23: #include <stdlib.h>
24: #include <string.h>
25:
26: #include "heccer/addressing.h"
27: #include "heccer/compartment.h"
28: #include "heccer/event.h"
29: #include "heccer/heccer.h"
30: #include "heccer/random.h"
31:
32:
33: static
34: int
35: HeccerCheckParameters
36: (struct Heccer *pheccer,
37: char *pcDescription,
38: ...);
39:
40: static
41: int
42: HeccerMechanismReadDoubleFile
43: (struct Heccer *pheccer, char *pcFilename, double **ppdValues);
44:
45:
46: /// **************************************************************************
47: ///
48: /// SHORT: HeccerCheckParameters()
49: ///
50: /// ARGS.:
51: ///
52: /// pheccer.......: a heccer.
53: /// pcDescription.: description of these parameters.
54: /// va_list.......: -1 terminated list of boolean values.
55: ///
56: /// RTN..: int
57: ///
58: /// success of operation.
59: ///
60: /// DESCR: Check parameters utility function.
61: ///
62: /// Just give a list of boolean expressions, telling if the
63: /// parameters comply or not, a description where they occur.
64: /// This function will call HeccerError() for any FALSE boolean in
65: /// the stdarg list.
66: ///
67: /// **************************************************************************
68:
69: static
70: int
71: HeccerCheckParameters
72: (struct Heccer *pheccer,
73: char *pcDescription,
74: ...)
75: {
76: //- set default result: ok
77:
78: int iResult = 1;
79:
80: //v stdargs list
81:
82: va_list vaList;
83:
84: //- get start of stdargs
85:
86: va_start(vaList, pcDescription);
87:
88: //- loop over all arguments
89:
90: int iOk = va_arg(vaList, int);
91:
92: while (iOk != -1)
93: {
94: //- if current argument not ok
95:
96: if (!iOk)
97: {
98: char pcMessage[1000];
99:
100: sprintf(pcMessage, "%s invalid", pcDescription);
101:
102: //t HeccerError(number, message, varargs);
103:
104: HeccerError
105: (pheccer,
106: NULL,
107: pcMessage);
108:
109: //- break loop
110:
111: break;
112: }
113:
114: //- go to next argument
115:
116: iOk = va_arg(vaList, int);
117: }
118:
119: //- end stdargs
120:
121: va_end(vaList);
122:
123: //- return result
124:
125: return(iResult);
126: }
127:
128:
129: /// **************************************************************************
130: ///
131: /// SHORT: HeccerMechanismCompile()
132: ///
133: /// ARGS.:
134: ///
135: /// pheccer...: a heccer.
136: ///
137: /// RTN..: int
138: ///
139: /// success of operation.
140: ///
141: /// DESCR: Compile the intermediary of the mechanisms to byte code.
142: ///
143: /// NOTE.:
144: ///
145: /// Compartment leak is a current, so it is considered to be a
146: /// mechanism.
147: ///
148: /// **************************************************************************
149:
150: int HeccerMechanismCompile(struct Heccer *pheccer)
151: {
152: //- set default result : ok
153:
154: int iResult = TRUE;
155:
156: //- first build the mechanism index
157:
158: if (!HeccerIntermediaryBuildIndex(pheccer))
159: {
160: return(FALSE);
161: }
162:
163: //v time step
164:
165: double dt;
166:
167: //v number of operators and operands
168:
169: int iMops;
170: int iMats;
171:
172: //- for backward euler integration
173:
174: if (pheccer->ho.iOptions & HECCER_OPTION_BACKWARD_EULER)
175: {
176: //- remember to do full time step
177:
178: dt = pheccer->dStep;
179: }
180:
181: //- else : crank-nicholson
182:
183: else
184: {
185: //- halve the time step
186:
187: dt = pheccer->dStep / 2.0;
188: }
189:
190: //- check parameters
191:
192: HeccerCheckParameters
193: (
194: pheccer,
195: "HeccerMechanismCompile(): time step",
196: (dt > 0 && dt < 1),
197: -1
198: );
199:
200: //- first count, then index, next compile the following block
201:
202: int iMopNumber;
203: int iMatNumber;
204:
205: int iCountIndexCompile;
206:
207: void *pvMops = NULL;
208: void *pvMats = NULL;
209: void **ppvCMatsIndex = NULL;
210: void **ppvMopsIndex = NULL;
211: void **ppvMatsIndex = NULL;
212: int *piMC2Mop = NULL;
213: uMC2Mat *piMC2Mat = NULL;
214:
215: for (iCountIndexCompile = 0 ; iCountIndexCompile < 3 ; iCountIndexCompile++)
216: {
217: //- counters always start at zero
218:
219: iMopNumber = 0;
220: iMatNumber = 0;
221:
222: iMops = 0;
223:
224: iMats = 0;
225:
226: //- loop over all compartments via their schedule number
227:
228: int iSchedule;
229:
230: for (iSchedule = 0 ; iSchedule < pheccer->inter.iCompartments ; iSchedule++)
231: {
232: //- fill in compartment operation
233:
234: SETMOP_COMPARTMENT(ppvMopsIndex, iMopNumber, pvMops, iMops);
235:
236: //! Em/Rm
237: //! injected current
238: //! dt/cm
239: //! diagonal
240: //!
241: //! injected current needs a separate entry for interfacing.
242:
243: //- get intermediary number for the current compartment
244:
245: int iIntermediary = pheccer->indexers.md.piBackward[iSchedule];
246:
247: //- retreive compartment constants
248:
249: double dCm = pheccer->inter.pcomp[iIntermediary].dCm;
250:
251: double dEm = pheccer->inter.pcomp[iIntermediary].dEm;
252:
253: //t perhaps better to do current injection with a
254: //t hardcoded injector callout ?
255:
256: double dInject = pheccer->inter.pcomp[iIntermediary].dInject;
257:
258: double dRm = pheccer->inter.pcomp[iIntermediary].dRm;
259:
260: //- check parameters
261:
262: char pcDescription[1000];
263:
264: sprintf(pcDescription, "HeccerMechanismCompile(): compartment %i parameters (dCm == %f, dEm == %f, dRm == %f)", iIntermediary, dCm, dEm, dRm);
265:
266: HeccerCheckParameters
267: (
268: pheccer,
269: pcDescription,
270: (dCm != 0),
271: (dEm >= -1 && dEm <= 1),
272: (dRm != 0),
273: -1
274: );
275:
276: //- fill in compartment constants
277:
278: //! note : pdDiagonals was computed with schedule numbers.
279:
280: //t these need an extra check, probably wrong.
281:
282: //t perhaps need to split SETMAT_COMPARTMENT in SETMAT_COMPARTMENT_START
283: //t and SETMAT_COMPARTMENT_FINISH
284: //t between those two, we compile in the mechanisms.
285:
286: SETMAT_COMPARTMENT(ppvCMatsIndex, iSchedule, ppvMatsIndex, iMatNumber, pvMats, iMats, dEm / dRm, dInject, dt / dCm, pheccer->vm.pdDiagonals[iSchedule]);
287:
288: //- loop over mechanisms for this compartment
289:
290: int iMathComponent;
291:
292: int iStart = iIntermediary == 0 ? 0 : pheccer->inter.piC2m[iIntermediary - 1];
293:
294: if (pheccer->inter.pmca && iStart > pheccer->inter.pmca->iMathComponents)
295: {
296: //t HeccerError(number, message, varargs);
297:
298: HeccerError
299: (pheccer,
300: NULL,
301: "trying to fetch math component number %i, which is out of range",
302: iStart);
303:
304: //- return error
305:
306: return(FALSE);
307: }
308:
309: //- lookup the start of the mechanisms for this compartment
310:
311: struct MathComponent *pmc = HeccerIntermediaryLookup(pheccer, iStart);
312:
313: for (iMathComponent = iStart ;
314: iMathComponent < pheccer->inter.piC2m[iIntermediary] ;
315: iMathComponent++)
316: {
317: //- look at mechanism type
318:
319: int iType = pmc->iType;
320:
321: switch (iType)
322: {
323: //- for a callout
324:
325: case MATH_TYPE_CallOut_conductance_current:
326: {
327: //- get type specific data
328:
329: struct Callout *pcall = (struct Callout *)pmc;
330:
331: pmc = MathComponentNext(&pcall->mc);
332:
333: SETMOP_CALLOUT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
334:
335: SETMAT_CALLOUT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcall);
336:
337: break;
338: }
339:
340: //- for a spring mass equation
341:
342: case MATH_TYPE_ChannelSpringMass:
343: {
344: //- get type specific data
345:
346: struct ChannelSpringMass * pcsm = (struct ChannelSpringMass *)pmc;
347:
348: pmc = MathComponentNext(&pcsm->mc);
349:
350: //- check parameters
351:
352: HeccerCheckParameters
353: (
354: pheccer,
355: "MATH_TYPE_ChannelSpringMass parameters",
356: (pcsm->dReversalPotential > -1 && pcsm->dReversalPotential < 1),
357: (pcsm->parameters.dTau1 != pcsm->parameters.dTau2),
358: -1
359: );
360:
361: //- compute conductance normalizer
362:
363: double dNormalizer;
364:
365: if (pcsm->parameters.dTau1 == pcsm->parameters.dTau2)
366: {
367: dNormalizer = pcsm->dMaximalConductance * M_E / pcsm->parameters.dTau1;
368: }
369: else
370: {
371: double dPeak
372: = (pcsm->parameters.dTau1
373: * pcsm->parameters.dTau2
374: * log(pcsm->parameters.dTau1 / pcsm->parameters.dTau2)
375: / (pcsm->parameters.dTau1 - pcsm->parameters.dTau2));
376:
377: dNormalizer
378: = (pcsm->dMaximalConductance
379: * (pcsm->parameters.dTau1 - pcsm->parameters.dTau2)
380: / (pcsm->parameters.dTau1
381: * pcsm->parameters.dTau2
382: * (exp(- dPeak / pcsm->parameters.dTau1)
383: - exp(- dPeak / pcsm->parameters.dTau2))));
384: }
385:
386: //- for a constant reversal potential
387:
388: if (pcsm->iReversalPotential == -1)
389: {
390: SETMOP_INITIALIZECHANNEL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, dNormalizer, pcsm->dReversalPotential);
391: }
392:
393: //- else a solved reversal potential
394:
395: else
396: {
397: //- get math component number
398:
399: int iMathComponentReversalPotential = pcsm->iReversalPotential;
400:
401: int iMatsReversalPotential = -1;
402:
403: if (iMathComponentReversalPotential != -1)
404: {
405: //- convert math component to mat number, convert mat number to mat addressable
406:
407: iMatsReversalPotential = piMC2Mat ? piMC2Mat[iMathComponentReversalPotential].iMat : -1;
408: }
409:
410: SETMOP_INITIALIZECHANNELEK(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, dNormalizer, iMatsReversalPotential);
411: }
412:
413: //- preprocess the event list if any
414:
415: if (pcsm->pcEventTimes
416: && pcsm->pdEventTimes)
417: {
418: //t HeccerError(number, message, varargs);
419:
420: HeccerError
421: (pheccer,
422: NULL,
423: "found a file specification as well as an array specification for event times of a springmass channel, "
424: "model container component number %i\n",
425: #ifdef HECCER_SOURCE_NEUROSPACES
426: pcsm->mc.iSerial,
427: #else
428: //t bhlkjwe, should have a char * here
429:
430: -1,
431: #endif
432: iStart);
433:
434: //- return error
435:
436: return(FALSE);
437: }
438:
439: if (pcsm->pcEventTimes)
440: {
441: //- read in the file
442:
443: int iDoubles = HeccerMechanismReadDoubleFile(pheccer, pcsm->pcEventTimes, &pcsm->pdEventTimes);
444:
445: if (iDoubles == -1
446: || !pcsm->pdEventTimes)
447: {
448: //t HeccerError(number, message, varargs);
449:
450: HeccerError
451: (pheccer,
452: NULL,
453: "could not read file %s, specified for a springmass channel, "
454: "model container component number %i\n",
455: pcsm->pcEventTimes,
456: #ifdef HECCER_SOURCE_NEUROSPACES
457: pcsm->mc.iSerial,
458: #else
459: //t bhlkjwe, should have a char * here
460:
461: -1,
462: #endif
463: iStart);
464:
465: //- return error
466:
467: return(FALSE);
468: }
469:
470: //! it would probably be a good idea to be able to go back to the original specification.
471:
472: //- mark end of array
473:
474: pcsm->pdEventTimes[iDoubles] = FLT_MAX;
475:
476: //- remove the reference to the file, for next loop over processing
477:
478: pcsm->pcEventTimes = NULL;
479: }
480:
481: //- tabulate the channel
482:
483: //- tabulate activation, Genesis X
484: //- create A table, alpha, create B table, alpha + beta
485:
486: int iTabulated = HeccerTabulateAny(pheccer, pcsm, MATH_TYPE_ChannelSpringMass);
487:
488: #ifdef HECCER_SOURCE_NEUROSPACES
489:
490: //- convert serial to index in the connection matrix
491:
492: int iDiscreteTarget = EventQueuerSerial2ConnectionIndex(pheccer->peq, ADDRESSING_HECCER_2_NEUROSPACES(pcsm->mc.iSerial));
493:
494: #else
495:
496: int iDiscreteTarget = -1;
497:
498: #endif
499:
500: SETMOP_SPRINGMASS(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcsm->pdEventTimes, iDiscreteTarget, pcsm->iTable, pheccer->dStep * pcsm->dFrequency);
501:
502: double dNextEvent = -1.0;
503:
504: SETMAT_SPRINGMASS(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcsm->dInitX, pcsm->dInitY, dNextEvent);
505:
506: SETMOP_UPDATECOMPARTMENTCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
507:
508: //t would it be useful to retabulate here ?
509:
510: //- register pool index
511:
512: //t for reasons of easy initialization, this should be a check for zero.
513: //t this means that I have to offset all mechanisms with 1
514: //t (mmm, the hines solver did the same, but for other reasons).
515:
516: if (pcsm->iPool != -1)
517: {
518: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
519:
520: SETMOP_FLUXPOOL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
521:
522: //! initial flux is assumed to be zero, always
523:
524: SETMAT_FLUXPOOL(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0);
525:
526: }
527:
528: //- compute individual channel contributions
529:
530: if ( (pheccer->ho.iOptions & HECCER_OPTION_ENABLE_INDIVIDUAL_CURRENTS) )
531: {
532: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
533:
534: SETMOP_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
535:
536: SETMAT_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0, 0.0);
537:
538: }
539:
540: //- compute aggregate current and conductance for this mathematical component type
541:
542: if ((pheccer->ho.iOptions & HECCER_OPTION_ENABLE_AGGREGATORS)
543: && pcsm->mc.iModelSourceType != -1)
544: {
545: int iAggregator = pcsm->mc.iModelSourceType;
546:
547: SETMOP_AGGREGATECURRENT(&pheccer->vm, iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, iAggregator);
548:
549: }
550:
551: //- register result from tabulation for outcome of this function
552:
553: iResult = iResult && iTabulated;
554:
555: if (!iResult)
556: {
557: HeccerError
558: (pheccer,
559: NULL,
560: "Compilation of ChannelSpringMass failed");
561: }
562:
563: break;
564: }
565:
566: //- for a nernst operation with internal variable concentration
567:
568: case MATH_TYPE_InternalNernst:
569: {
570: //- get type specific data
571:
572: struct InternalNernst * pin = (struct InternalNernst *)pmc;
573:
574: pmc = MathComponentNext(&pin->mc);
575:
576: //- check parameters
577:
578: HeccerCheckParameters
579: (
580: pheccer,
581: "MATH_TYPE_InternalNernst parameters",
582: -1
583: );
584:
585: //- get math component number
586:
587: int iMathComponentActivator = pin->iInternal;
588:
589: int iMatsActivator = -1;
590:
591: if (iMathComponentActivator != -1)
592: {
593: //- convert math component to mat number, convert mat number to mat addressable
594:
595: iMatsActivator = piMC2Mat ? piMC2Mat[iMathComponentActivator].iMat : -1;
596: }
597:
598: SETMOP_INTERNALNERNST(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pin->dConstant, pin->dExternal, iMatsActivator);
599:
600: SETMAT_INTERNALNERNST(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pin->dInitPotential);
601:
602: break;
603: }
604:
605: //- for a channel with only activation
606:
607: case MATH_TYPE_ChannelAct:
608: {
609: //- get type specific data
610:
611: struct ChannelAct *pca = (struct ChannelAct *)pmc;
612:
613: pmc = MathComponentNext(&pca->mc);
614:
615: //- check parameters
616:
617: HeccerCheckParameters
618: (
619: pheccer,
620: "MATH_TYPE_ChannelAct parameters",
621: -1
622: );
623:
624: //- for a constant reversal potential
625:
626: if (pca->iReversalPotential == -1)
627: {
628: SETMOP_INITIALIZECHANNEL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pca->dMaximalConductance, pca->dReversalPotential);
629: }
630:
631: //- else a solved reversal potential
632:
633: else
634: {
635: //- get math component number
636:
637: int iMathComponentReversalPotential = pca->iReversalPotential;
638:
639: int iMatsReversalPotential = -1;
640:
641: if (iMathComponentReversalPotential != -1)
642: {
643: //- convert math component to mat number, convert mat number to mat addressable
644:
645: iMatsReversalPotential = piMC2Mat ? piMC2Mat[iMathComponentReversalPotential].iMat : -1;
646: }
647:
648: SETMOP_INITIALIZECHANNELEK(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pca->dMaximalConductance, iMatsReversalPotential);
649: }
650:
651: //- tabulate the channel
652:
653: //- tabulate activation, Genesis X
654: //- create A table, alpha, create B table, alpha + beta
655:
656: int iTabulated
657: = HeccerTabulateAny(pheccer, &pca->pgc.gc, MATH_TYPE_GateConcept);
658:
659: SETMOP_LOADVOLTAGETABLE(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
660:
661: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pca->pgc.gc.iTable, pca->pgc.iPower,-1);
662:
663: //! at the beginning of a simulation, you would expect this to be the steady state value
664:
665: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pca->pgc.gc.dInitActivation);
666:
667: SETMOP_UPDATECOMPARTMENTCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
668:
669: //t retabulate cannot be done yet, do not know yet how many tables
670:
671: //- register pool index
672:
673: //t for reasons of easy initialization, this should be a check for zero.
674: //t this means that I have to offset all mechanisms with 1
675: //t (mmm, the hines solver did the same, but for other reasons).
676:
677: if (pca->iPool != -1)
678: {
679: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
680:
681: SETMOP_FLUXPOOL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
682:
683: //! initial flux is assumed to be zero, always
684:
685: SETMAT_FLUXPOOL(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0);
686:
687: }
688:
689: //- compute individual channel contributions
690:
691: if ( (pheccer->ho.iOptions & HECCER_OPTION_ENABLE_INDIVIDUAL_CURRENTS) )
692: {
693: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
694:
695: SETMOP_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
696:
697: SETMAT_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0, 0.0);
698:
699: }
700:
701: //- compute aggregate current and conductance for this mathematical component type
702:
703: if ((pheccer->ho.iOptions & HECCER_OPTION_ENABLE_AGGREGATORS)
704: && pca->mc.iModelSourceType != -1)
705: {
706: int iAggregator = pca->mc.iModelSourceType;
707:
708: SETMOP_AGGREGATECURRENT(&pheccer->vm, iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, iAggregator);
709:
710: }
711:
712: //- register result from tabulation for outcome of this function
713:
714: iResult = iResult && iTabulated;
715:
716: if (!iResult)
717: {
718: HeccerError
719: (pheccer,
720: NULL,
721: "Compilation of ChannelAct failed");
722: }
723:
724: break;
725: }
726:
727: //- for a regular channel with activation and inactivation
728:
729: case MATH_TYPE_ChannelActInact:
730: {
731: //- get type specific data
732:
733: struct ChannelActInact *pcai = (struct ChannelActInact *)pmc;
734:
735: pmc = MathComponentNext(&pcai->mc);
736:
737: //- check parameters
738:
739: HeccerCheckParameters
740: (
741: pheccer,
742: "MATH_TYPE_ChannelActInact parameters",
743: -1
744: );
745:
746: //- for a constant reversal potential
747:
748: if (pcai->iReversalPotential == -1)
749: {
750: SETMOP_INITIALIZECHANNEL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcai->dMaximalConductance, pcai->dReversalPotential);
751: }
752:
753: //- else a solved reversal potential
754:
755: else
756: {
757: //- get math component number
758:
759: int iMathComponentReversalPotential = pcai->iReversalPotential;
760:
761: int iMatsReversalPotential = -1;
762:
763: if (iMathComponentReversalPotential != -1)
764: {
765: //- convert math component to mat number, convert mat number to mat addressable
766:
767: iMatsReversalPotential = piMC2Mat ? piMC2Mat[iMathComponentReversalPotential].iMat : -1;
768: }
769:
770: SETMOP_INITIALIZECHANNELEK(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcai->dMaximalConductance, iMatsReversalPotential);
771: }
772:
773: //- tabulate the channel
774:
775: //- tabulate activation, Genesis X
776: //- create A table, alpha, create B table, alpha + beta
777:
778: int iTabulatedActivation
779: = HeccerTabulateAny(pheccer, &pcai->pgcActivation.gc, MATH_TYPE_GateConcept);
780:
781: SETMOP_LOADVOLTAGETABLE(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
782:
783: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcai->pgcActivation.gc.iTable, pcai->pgcActivation.iPower,-1);
784:
785: //! at the beginning of a simulation, you would expect this to be the steady state value
786:
787: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcai->pgcActivation.gc.dInitActivation);
788:
789: //- tabulate inactivation, Genesis Y
790: //- create A table, alpha, create B table, alpha + beta
791:
792: int iTabulatedInactivation
793: = HeccerTabulateAny(pheccer, &pcai->pgcInactivation.gc, MATH_TYPE_GateConcept);
794:
795: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcai->pgcInactivation.gc.iTable, pcai->pgcInactivation.iPower,-1);
796:
797: //! at the beginning of a simulation, you would expect this to be the steady state value
798:
799: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcai->pgcInactivation.gc.dInitActivation);
800:
801: SETMOP_UPDATECOMPARTMENTCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
802:
803: //t retabulate cannot be done yet, do not know yet how many tables
804:
805: //- register pool index
806:
807: //t for reasons of easy initialization, this should be a check for zero.
808: //t this means that I have to offset all mechanisms with 1
809: //t (mmm, the hines solver did the same, but for other reasons).
810:
811: if (pcai->iPool != -1)
812: {
813: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
814:
815: SETMOP_FLUXPOOL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
816:
817: //! initial flux is assumed to be zero, always
818:
819: SETMAT_FLUXPOOL(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0);
820:
821: }
822:
823: //- compute individual channel contributions
824:
825: if ( (pheccer->ho.iOptions & HECCER_OPTION_ENABLE_INDIVIDUAL_CURRENTS) )
826: {
827: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
828:
829: SETMOP_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
830:
831: SETMAT_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0, 0.0);
832:
833: }
834:
835: //- compute aggregate current and conductance for this mathematical component type
836:
837: if ((pheccer->ho.iOptions & HECCER_OPTION_ENABLE_AGGREGATORS)
838: && pcai->mc.iModelSourceType != -1)
839: {
840: int iAggregator = pcai->mc.iModelSourceType;
841:
842: SETMOP_AGGREGATECURRENT(&pheccer->vm, iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, iAggregator);
843:
844: }
845:
846: //- register result from tabulation for outcome of this function
847:
848: iResult = iResult && iTabulatedActivation && iTabulatedInactivation;
849:
850: if (!iResult)
851: {
852: HeccerError
853: (pheccer,
854: NULL,
855: "Compilation of ChannelActInact failed");
856: }
857:
858: break;
859: }
860:
861: //- for a channel with a potential and an external dependence
862:
863: case MATH_TYPE_ChannelActConc:
864: {
865: //- get type specific data
866:
867: struct ChannelActConc *pcac = (struct ChannelActConc *)pmc;
868:
869: pmc = MathComponentNext(&pcac->mc);
870:
871: //- check parameters
872:
873: HeccerCheckParameters
874: (
875: pheccer,
876: "MATH_TYPE_ChannelActConc parameters",
877: -1
878: );
879:
880: //- for a constant reversal potential
881:
882: if (pcac->iReversalPotential == -1)
883: {
884: SETMOP_INITIALIZECHANNEL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcac->dMaximalConductance, pcac->dReversalPotential);
885: }
886:
887: //- else a solved reversal potential
888:
889: else
890: {
891: //- get math component number
892:
893: int iMathComponentReversalPotential = pcac->iReversalPotential;
894:
895: int iMatsReversalPotential = -1;
896:
897: if (iMathComponentReversalPotential != -1)
898: {
899: //- convert math component to mat number, convert mat number to mat addressable
900:
901: iMatsReversalPotential = piMC2Mat ? piMC2Mat[iMathComponentReversalPotential].iMat : -1;
902: }
903:
904: SETMOP_INITIALIZECHANNELEK(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcac->dMaximalConductance, iMatsReversalPotential);
905: }
906:
907: //- tabulate the membrane dependence
908:
909: //- tabulate membrane dependence, Genesis X
910: //- create A table, alpha, create B table, alpha + beta
911:
912: int iTabulatedMembraneDependence
913: = HeccerTabulateAny(pheccer, &pcac->pgc.gc, MATH_TYPE_GateConcept);
914:
915: SETMOP_LOADVOLTAGETABLE(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
916:
917: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcac->pgc.gc.iTable, pcac->pgc.iPower,-1);
918:
919: //! at the beginning of a simulation, you would expect this to be the steady state value
920:
921: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcac->pgc.gc.dInitActivation);
922:
923: //- tabulate concentration dependence, Genesis Z
924: //- create A table, alpha, create B table, alpha + beta
925:
926: int iTabulatedBasalActivator
927: = HeccerTabulateAny(pheccer, &pcac->pac.ca, MATH_TYPE_Concentration);
928:
929: //! gate computations are just fetching things from tables, and
930: //! multiplying the conductances, so it is not relevant if these
931: //! computations are done for membrane potential dependent gates or
932: //! concentration dependent gates.
933:
934: //- get math component number
935:
936: int iMathComponentActivator = pcac->pac.ca.iActivator;
937:
938: int iMatsActivator = -1;
939:
940: if (iMathComponentActivator != -1)
941: {
942: //- convert math component to mat number, convert mat number to mat addressable
943:
944: iMatsActivator = piMC2Mat ? piMC2Mat[iMathComponentActivator].iMat : -1;
945: }
946:
947: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcac->pac.ca.iTable, pcac->pac.iPower, iMatsActivator);
948:
949: //! at the beginning of a simulation, you would expect this to be the steady state value
950:
951: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcac->pac.ca.dInitActivation);
952:
953: SETMOP_UPDATECOMPARTMENTCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
954:
955: //t retabulate cannot be done yet, do not know yet how many tables
956:
957: //- register pool index
958:
959: //t for reasons of easy initialization, this should be a check for zero.
960: //t this means that I have to offset all mechanisms with 1
961: //t (mmm, the hines solver did the same, but for other reasons).
962:
963: if (pcac->iPool != -1)
964: {
965: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
966:
967: SETMOP_FLUXPOOL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
968:
969: //! initial flux is assumed to be zero, always
970:
971: SETMAT_FLUXPOOL(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0);
972:
973: }
974:
975: //- compute individual channel contributions
976:
977: if ( (pheccer->ho.iOptions & HECCER_OPTION_ENABLE_INDIVIDUAL_CURRENTS) )
978: {
979: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
980:
981: SETMOP_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
982:
983: SETMAT_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0, 0.0);
984:
985: }
986:
987: //- compute aggregate current and conductance for this mathematical component type
988:
989: if ((pheccer->ho.iOptions & HECCER_OPTION_ENABLE_AGGREGATORS)
990: && pcac->mc.iModelSourceType != -1)
991: {
992: int iAggregator = pcac->mc.iModelSourceType;
993:
994: SETMOP_AGGREGATECURRENT(&pheccer->vm, iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, iAggregator);
995:
996: }
997:
998: //- register result from tabulation for outcome of this function
999:
1000: iResult = iResult && iTabulatedMembraneDependence && iTabulatedBasalActivator;
1001:
1002: if (!iResult)
1003: {
1004: HeccerError
1005: (pheccer,
1006: NULL,
1007: "Compilation of ChannelActConc failed");
1008: }
1009:
1010: break;
1011: }
1012:
1013: //- for an exponential decaying variable
1014:
1015: case MATH_TYPE_ExponentialDecay:
1016: {
1017: //- get type specific data
1018:
1019: struct ExponentialDecay *pexdec = (struct ExponentialDecay *)pmc;
1020:
1021: pmc = MathComponentNext(&pexdec->mc);
1022:
1023: //- check parameters
1024:
1025: HeccerCheckParameters
1026: (
1027: pheccer,
1028: "MATH_TYPE_ExponentialDecay parameters",
1029: -1
1030: );
1031:
1032: int piExternal[EXPONENTIALDECAY_CONTRIBUTORS];
1033:
1034: int i;
1035:
1036: for (i = 0 ; i < EXPONENTIALDECAY_CONTRIBUTORS ; i++)
1037: {
1038: //- get math component number
1039:
1040: int iContributor = pexdec->piExternal[i];
1041:
1042: int iMatsExternal = -1;
1043:
1044: if (iContributor != -1)
1045: {
1046: //- convert math component to mat number, convert mat number to mat addressable
1047:
1048: iMatsExternal = piMC2Mat ? piMC2Mat[iContributor].iMat : -1;
1049:
1050: piExternal[i] = iMatsExternal;
1051: }
1052: else
1053: {
1054: piExternal[i] = -1;
1055: }
1056: }
1057:
1058: SETMOP_EXPONENTIALDECAY(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pheccer->dStep * pexdec->dBeta, pexdec->dSteadyState, 1 + pheccer->dStep / (2 * pexdec->dTau), piExternal);
1059:
1060: SETMAT_EXPONENTIALDECAY(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pexdec->dInitValue);
1061:
1062: break;
1063: }
1064:
1065: //- for a channel specified as steady state and a stepped tau
1066:
1067: case MATH_TYPE_ChannelSteadyStateSteppedTau:
1068: {
1069: //- get type specific data
1070:
1071: struct ChannelSteadyStateSteppedTau *pcsst = (struct ChannelSteadyStateSteppedTau *)pmc;
1072:
1073: pmc = MathComponentNext(&pcsst->mc);
1074:
1075: //- check parameters
1076:
1077: HeccerCheckParameters
1078: (
1079: pheccer,
1080: "MATH_TYPE_ChannelSteadyStateSteppedTau parameters",
1081: -1
1082: );
1083:
1084: //- for a constant reversal potential
1085:
1086: if (pcsst->iReversalPotential == -1)
1087: {
1088: SETMOP_INITIALIZECHANNEL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcsst->dMaximalConductance, pcsst->dReversalPotential);
1089: }
1090:
1091: //- else a solved reversal potential
1092:
1093: else
1094: {
1095: //- get math component number
1096:
1097: int iMathComponentReversalPotential = pcsst->iReversalPotential;
1098:
1099: int iMatsReversalPotential = -1;
1100:
1101: if (iMathComponentReversalPotential != -1)
1102: {
1103: //- convert math component to mat number, convert mat number to mat addressable
1104:
1105: iMatsReversalPotential = piMC2Mat ? piMC2Mat[iMathComponentReversalPotential].iMat : -1;
1106: }
1107:
1108: SETMOP_INITIALIZECHANNELEK(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcsst->dMaximalConductance, iMatsReversalPotential);
1109: }
1110:
1111: //- tabulate the channel
1112:
1113: int iTabulated
1114: = HeccerTabulateAny(pheccer, pcsst, MATH_TYPE_ChannelSteadyStateSteppedTau);
1115:
1116: SETMOP_LOADVOLTAGETABLE(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1117:
1118: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcsst->iFirstTable, pcsst->iFirstPower, -1);
1119:
1120: //! at the beginning of a simulation, you would expect this to be the steady state value
1121:
1122: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcsst->dFirstInitActivation);
1123:
1124: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcsst->iSecondTable, pcsst->iSecondPower, -1);
1125:
1126: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcsst->dSecondInitActivation);
1127:
1128: SETMOP_UPDATECOMPARTMENTCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1129:
1130: //t retabulate cannot be done yet, do not know yet how many tables
1131:
1132: //- register pool index
1133:
1134: //t for reasons of easy initialization, this should be a check for zero.
1135: //t this means that I have to offset all mechanisms with 1
1136: //t (mmm, the hines solver did the same, but for other reasons).
1137:
1138: if (pcsst->iPool != -1)
1139: {
1140: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1141:
1142: SETMOP_FLUXPOOL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1143:
1144: //! initial flux is assumed to be zero, always
1145:
1146: SETMAT_FLUXPOOL(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0);
1147:
1148: }
1149:
1150: //- compute individual channel contributions
1151:
1152: if ( (pheccer->ho.iOptions & HECCER_OPTION_ENABLE_INDIVIDUAL_CURRENTS) )
1153: {
1154: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1155:
1156: SETMOP_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1157:
1158: SETMAT_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0, 0.0);
1159:
1160: }
1161:
1162: //- compute aggregate current and conductance for this mathematical component type
1163:
1164: if ((pheccer->ho.iOptions & HECCER_OPTION_ENABLE_AGGREGATORS)
1165: && pcsst->mc.iModelSourceType != -1)
1166: {
1167: int iAggregator = pcsst->mc.iModelSourceType;
1168:
1169: SETMOP_AGGREGATECURRENT(&pheccer->vm, iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, iAggregator);
1170:
1171: }
1172:
1173: //- register result from tabulation for outcome of this function
1174:
1175: iResult = iResult && iTabulated;
1176:
1177: if (!iResult)
1178: {
1179: HeccerError
1180: (pheccer,
1181: NULL,
1182: "Compilation of ChannelSteadyStateSteppedTau failed");
1183: }
1184:
1185: break;
1186: }
1187:
1188: //- for a persistent channel specified as steady state and two tau constants
1189:
1190: case MATH_TYPE_ChannelPersistentSteadyStateDualTau:
1191: {
1192: //- get type specific data
1193:
1194: struct ChannelPersistentSteadyStateDualTau *pcpsdt = (struct ChannelPersistentSteadyStateDualTau *)pmc;
1195:
1196: pmc = MathComponentNext(&pcpsdt->mc);
1197:
1198: //- check parameters
1199:
1200: HeccerCheckParameters
1201: (
1202: pheccer,
1203: "MATH_TYPE_ChannelPersistentSteadyStateDualTau parameters",
1204: -1
1205: );
1206:
1207: //- for a constant reversal potential
1208:
1209: if (pcpsdt->iReversalPotential == -1)
1210: {
1211: SETMOP_INITIALIZECHANNEL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcpsdt->dMaximalConductance, pcpsdt->dReversalPotential);
1212: }
1213:
1214: //- else a solved reversal potential
1215:
1216: else
1217: {
1218: //- get math component number
1219:
1220: int iMathComponentReversalPotential = pcpsdt->iReversalPotential;
1221:
1222: int iMatsReversalPotential = -1;
1223:
1224: if (iMathComponentReversalPotential != -1)
1225: {
1226: //- convert math component to mat number, convert mat number to mat addressable
1227:
1228: iMatsReversalPotential = piMC2Mat ? piMC2Mat[iMathComponentReversalPotential].iMat : -1;
1229: }
1230:
1231: SETMOP_INITIALIZECHANNELEK(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcpsdt->dMaximalConductance, iMatsReversalPotential);
1232: }
1233:
1234: //- tabulate the channel
1235:
1236: int iTabulated
1237: = HeccerTabulateAny(pheccer, pcpsdt, MATH_TYPE_ChannelPersistentSteadyStateDualTau);
1238:
1239: SETMOP_LOADVOLTAGETABLE(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1240:
1241: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcpsdt->iFirstTable, pcpsdt->iFirstPower, -1);
1242:
1243: //! at the beginning of a simulation, you would expect this to be the steady state value
1244:
1245: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcpsdt->dFirstInitActivation);
1246:
1247: SETMOP_UPDATECOMPARTMENTCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1248:
1249: SETMOP_INITIALIZECHANNEL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcpsdt->dMaximalConductance, pcpsdt->dReversalPotential);
1250:
1251: SETMOP_LOADVOLTAGETABLE(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1252:
1253: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcpsdt->iSecondTable, pcpsdt->iSecondPower, -1);
1254:
1255: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcpsdt->dSecondInitActivation);
1256:
1257: SETMOP_UPDATECOMPARTMENTCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1258:
1259: //t retabulate cannot be done yet, do not know yet how many tables
1260:
1261: //- register pool index
1262:
1263: //t for reasons of easy initialization, this should be a check for zero.
1264: //t this means that I have to offset all mechanisms with 1
1265: //t (mmm, the hines solver did the same, but for other reasons).
1266:
1267: if (pcpsdt->iPool != -1)
1268: {
1269: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1270:
1271: SETMOP_FLUXPOOL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1272:
1273: //! initial flux is assumed to be zero, always
1274:
1275: SETMAT_FLUXPOOL(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0);
1276:
1277: }
1278:
1279: //- compute individual channel contributions
1280:
1281: if ( (pheccer->ho.iOptions & HECCER_OPTION_ENABLE_INDIVIDUAL_CURRENTS) )
1282: {
1283: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1284:
1285: SETMOP_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1286:
1287: SETMAT_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0, 0.0);
1288:
1289: }
1290:
1291: //- compute aggregate current and conductance for this mathematical component type
1292:
1293: if ((pheccer->ho.iOptions & HECCER_OPTION_ENABLE_AGGREGATORS)
1294: && pcpsdt->mc.iModelSourceType != -1)
1295: {
1296: int iAggregator = pcpsdt->mc.iModelSourceType;
1297:
1298: SETMOP_AGGREGATECURRENT(&pheccer->vm, iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, iAggregator);
1299:
1300: }
1301:
1302: //- register result from tabulation for outcome of this function
1303:
1304: iResult = iResult && iTabulated;
1305:
1306: if (!iResult)
1307: {
1308: HeccerError
1309: (pheccer,
1310: NULL,
1311: "Compilation of ChannelPersistentSteadyStateDualTau failed");
1312: }
1313:
1314: break;
1315: }
1316:
1317: //- for a persistent channel specified as steady state and a tau
1318:
1319: case MATH_TYPE_ChannelPersistentSteadyStateTau:
1320: {
1321: //- get type specific data
1322:
1323: struct ChannelPersistentSteadyStateTau *pcpst = (struct ChannelPersistentSteadyStateTau *)pmc;
1324:
1325: pmc = MathComponentNext(&pcpst->mc);
1326:
1327: //- check parameters
1328:
1329: HeccerCheckParameters
1330: (
1331: pheccer,
1332: "MATH_TYPE_ChannelPersistentSteadyStateTau parameters",
1333: -1
1334: );
1335:
1336: //- for a constant reversal potential
1337:
1338: if (pcpst->iReversalPotential == -1)
1339: {
1340: SETMOP_INITIALIZECHANNEL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcpst->dMaximalConductance, pcpst->dReversalPotential);
1341: }
1342:
1343: //- else a solved reversal potential
1344:
1345: else
1346: {
1347: //- get math component number
1348:
1349: int iMathComponentReversalPotential = pcpst->iReversalPotential;
1350:
1351: int iMatsReversalPotential = -1;
1352:
1353: if (iMathComponentReversalPotential != -1)
1354: {
1355: //- convert math component to mat number, convert mat number to mat addressable
1356:
1357: iMatsReversalPotential = piMC2Mat ? piMC2Mat[iMathComponentReversalPotential].iMat : -1;
1358: }
1359:
1360: SETMOP_INITIALIZECHANNELEK(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcpst->dMaximalConductance, iMatsReversalPotential);
1361: }
1362:
1363: //- tabulate the channel
1364:
1365: int iTabulated
1366: = HeccerTabulateAny(pheccer, pcpst, MATH_TYPE_ChannelPersistentSteadyStateTau);
1367:
1368: SETMOP_LOADVOLTAGETABLE(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1369:
1370: SETMOP_POWEREDGATECONCEPT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, pcpst->iTable, pcpst->iPower, -1);
1371:
1372: //! at the beginning of a simulation, you would expect this to be the steady state value
1373:
1374: SETMAT_POWEREDGATECONCEPT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, pcpst->dInitActivation);
1375:
1376: SETMOP_UPDATECOMPARTMENTCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1377:
1378: //t retabulate cannot be done yet, do not know yet how many tables
1379:
1380: //- register pool index
1381:
1382: //t for reasons of easy initialization, this should be a check for zero.
1383: //t this means that I have to offset all mechanisms with 1
1384: //t (mmm, the hines solver did the same, but for other reasons).
1385:
1386: if (pcpst->iPool != -1)
1387: {
1388: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1389:
1390: SETMOP_FLUXPOOL(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1391:
1392: //! initial flux is assumed to be zero, always
1393:
1394: SETMAT_FLUXPOOL(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0);
1395:
1396: }
1397:
1398: //- compute individual channel contributions
1399:
1400: if ( (pheccer->ho.iOptions & HECCER_OPTION_ENABLE_INDIVIDUAL_CURRENTS) )
1401: {
1402: SETMOP_REGISTERCHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1403:
1404: SETMOP_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops);
1405:
1406: SETMAT_STORESINGLECHANNELCURRENT(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, 0.0, 0.0);
1407:
1408: }
1409:
1410: //- compute aggregate current and conductance for this mathematical component type
1411:
1412: if ((pheccer->ho.iOptions & HECCER_OPTION_ENABLE_AGGREGATORS)
1413: && pcpst->mc.iModelSourceType != -1)
1414: {
1415: int iAggregator = pcpst->mc.iModelSourceType;
1416:
1417: SETMOP_AGGREGATECURRENT(&pheccer->vm, iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, iAggregator);
1418:
1419: }
1420:
1421: //- register result from tabulation for outcome of this function
1422:
1423: iResult = iResult && iTabulated;
1424:
1425: if (!iResult)
1426: {
1427: HeccerError
1428: (pheccer,
1429: NULL,
1430: "Compilation of ChannelPersistentSteadyStateTau failed");
1431: }
1432:
1433: break;
1434: }
1435:
1436: case MATH_TYPE_SpikeGenerator:
1437: {
1438: //- get type specific data
1439:
1440: struct SpikeGenerator *psg = (struct SpikeGenerator *)pmc;
1441:
1442: pmc = MathComponentNext(&psg->mc);
1443:
1444: //- check parameters
1445:
1446: HeccerCheckParameters
1447: (
1448: pheccer,
1449: "MATH_TYPE_SpikeGenerator parameters",
1450: -1
1451: );
1452:
1453: //- set operators and operands
1454:
1455: //! INT_MAX means the membrane potential is the source
1456:
1457: //! for other things, fill in the matindex of the
1458: //! source, the linker will link the mechanisms
1459: //! together (untested).
1460:
1461: int iSource = INT_MAX;
1462:
1463: int iTable = psg->iTable;
1464:
1465: if (iTable == INT_MAX)
1466: {
1467: //t from ssp viewpoint, the distributor_service is already setup, and knows that it connects to this spikegen
1468:
1469: //t the distributor_service builds a table that converts model container serials to entries in the connection matrix
1470:
1471: //t here we use that table, to convert the serial to the index in the connection matrix
1472:
1473:
1474: #ifdef HECCER_SOURCE_NEUROSPACES
1475: int iSerial = psg->mc.iSerial;
1476: #else
1477: int iSerial = INT_MAX;
1478: #endif
1479:
1480: struct EventDistributor *ped = pheccer->ped;
1481:
1482: //t EventDistributor does not know about anything yet, this cannot work.
1483:
1484: //t we need to register (indirectly) the
1485: //t memory location of iTable, probably by
1486: //t storing the iMopNumber in an index intended
1487: //t for resolving connections after CompileP3.
1488:
1489: iTable = EventDistributorSerial2Index(ped, ADDRESSING_HECCER_2_NEUROSPACES(iSerial));
1490:
1491: //t so hardcoded solution that makes the spiker1 test case work
1492:
1493: iTable = 0;
1494: }
1495:
1496: SETMOP_EVENTGENERATE(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, iSource, psg->dThreshold, psg->dRefractory, iTable);
1497:
1498: //t we are not in the refractory period, check randomspike2 for initial refractory probability calculation.
1499:
1500: double dRefractoryTime = -1.0;
1501:
1502: double dSpike = 0.0;
1503:
1504: SETMAT_EVENTGENERATE(iMathComponent, piMC2Mat, ppvMatsIndex, iMatNumber, pvMats, iMats, dRefractoryTime, dSpike);
1505:
1506: //- if there is a reset value set
1507:
1508: if (psg->dReset != FLT_MAX)
1509: {
1510: //- generate op for resetting the membrane potential
1511:
1512: SETMOP_RESET(iMathComponent, piMC2Mop, ppvMopsIndex, iMopNumber, pvMops, iMops, psg->dReset);
1513:
1514: }
1515:
1516: break;
1517: }
1518:
1519: default:
1520: {
1521: //t HeccerError(number, message, varargs);
1522:
1523: HeccerError
1524: (pheccer,
1525: NULL,
1526: "unknown pmc->iType (%i)",
1527: iType);
1528:
1529: break;
1530: }
1531: }
1532: }
1533: }
1534:
1535: //- sanity : is next compartment's mechanism invalid ?
1536:
1537: //! so pheccer->inter.piC2m can be NULL if no compartments have been found.
1538:
1539: if (pheccer->inter.piC2m
1540: && pheccer->inter.piC2m[iSchedule] != -1)
1541: {
1542: //t HeccerError(number, message, varargs);
1543:
1544: HeccerError
1545: (pheccer,
1546: NULL,
1547: "mechanisms found after last compartment's mechanism");
1548:
1549: return(FALSE);
1550: }
1551:
1552: //- finish all operations
1553:
1554: SETMOP_FINISH(ppvMopsIndex, iMopNumber, pvMops, iMops);
1555:
1556: //- if we were counting during the previous loop
1557:
1558: if (iCountIndexCompile == 0)
1559: {
1560: pheccer->vm.iMopNumber = iMopNumber;
1561:
1562: pheccer->vm.ppvCMatsIndex = (void **)calloc(pheccer->inter.iCompartments + 1, sizeof(void *));
1563:
1564: ppvCMatsIndex = pheccer->vm.ppvCMatsIndex;
1565:
1566: pheccer->vm.ppvMopsIndex = (void **)calloc(iMopNumber + 1, sizeof(void *));
1567:
1568: ppvMopsIndex = pheccer->vm.ppvMopsIndex;
1569:
1570: pheccer->vm.iMatNumber = iMatNumber;
1571:
1572: pheccer->vm.ppvMatsIndex = (void **)calloc(iMatNumber + 1, sizeof(void *));
1573:
1574: ppvMatsIndex = pheccer->vm.ppvMatsIndex;
1575:
1576: if (pheccer->inter.pmca)
1577: {
1578: //! note that this one does not index compartments, only the mechanism math components.
1579:
1580: pheccer->vm.piMC2Mat = (uMC2Mat *)calloc(pheccer->inter.pmca->iMathComponents + 1, sizeof(uMC2Mat));
1581:
1582: piMC2Mat = pheccer->vm.piMC2Mat;
1583:
1584: pheccer->vm.piMC2Mop = (int *)calloc(pheccer->inter.pmca->iMathComponents + 1, sizeof(int));
1585:
1586: piMC2Mop = pheccer->vm.piMC2Mop;
1587: }
1588: }
1589:
1590: //- if we were indexing during the previous loop
1591:
1592: else if (iCountIndexCompile == 1)
1593: {
1594: //- prepare for compilation : allocate operators and addressables
1595:
1596: pheccer->vm.pvMops = (void *)calloc(iMops, 1);
1597:
1598: pvMops = pheccer->vm.pvMops;
1599:
1600: pheccer->vm.iMops = iMops;
1601:
1602: pheccer->vm.pvMats = (void *)calloc(iMats, 1);
1603:
1604: pvMats = pheccer->vm.pvMats;
1605:
1606: pheccer->vm.iMats = iMats;
1607: }
1608: }
1609:
1610: //- return result
1611:
1612: return(iResult);
1613: }
1614:
1615:
1616: /// **************************************************************************
1617: ///
1618: /// SHORT: HeccerMechanismLink()
1619: ///
1620: /// ARGS.:
1621: ///
1622: /// pheccer...: a heccer.
1623: ///
1624: /// RTN..: int
1625: ///
1626: /// success of operation.
1627: ///
1628: /// DESCR: Link mechanism operations.
1629: ///
1630: /// **************************************************************************
1631:
1632: int HeccerMechanismLink(struct Heccer *pheccer)
1633: {
1634: //- set default result : ok
1635:
1636: int iResult = TRUE;
1637:
1638: int *piMop = (int *)pheccer->vm.pvMops;
1639:
1640: void *pvMats = pheccer->vm.pvMats;
1641:
1642: //- loop over mechanism operators
1643:
1644: while (piMop[0] == HECCER_MOP_COMPARTMENT)
1645: {
1646: //- go to next operator
1647:
1648: piMop = &piMop[1];
1649:
1650: //- get compartment mechanism data
1651:
1652: struct MatsCompartment *pmatsc
1653: = (struct MatsCompartment *)pvMats;
1654:
1655: pvMats = (void *)&((struct MatsCompartment *)pvMats)[1];
1656:
1657: //- loop over mechanism operators
1658:
1659: while (piMop[0] > HECCER_MOP_COMPARTMENT_BARRIER)
1660: {
1661: //- look at next mechanism
1662:
1663: switch (piMop[0])
1664: {
1665: //- for aggregate current
1666:
1667: case HECCER_MOP_AGGREGATECURRENT:
1668: {
1669: //- go to next operator
1670:
1671: struct MopsAggregateCurrent *pmops = (struct MopsAggregateCurrent *)piMop;
1672:
1673: piMop = (int *)&pmops[1];
1674:
1675: break;
1676: }
1677:
1678: //- for a call out
1679:
1680: case HECCER_MOP_CALLOUT:
1681: {
1682: //- go to next operator
1683:
1684: piMop = &piMop[1];
1685:
1686: //- go to next type specific data
1687:
1688: struct MatsCallout * pcall = (struct MatsCallout *)pvMats;
1689:
1690: pvMats = (void *)&((struct MatsCallout *)pvMats)[1];
1691:
1692: break;
1693: }
1694:
1695: //- for a spring mass equation
1696:
1697: case HECCER_MOP_SPRINGMASS:
1698: {
1699: //- go to next operator
1700:
1701: struct MopsSpringMass *pmops = (struct MopsSpringMass *)piMop;
1702:
1703: piMop = (int *)&pmops[1];
1704:
1705: //- go to next type specific data
1706:
1707: struct MatsSpringMass *pmats = (struct MatsSpringMass *)pvMats;
1708:
1709: pvMats = (void *)&((struct MatsSpringMass *)pvMats)[1];
1710:
1711: //t iDiscreteSource must be linked here ?
1712:
1713: break;
1714: }
1715:
1716: //- for a nernst operation with internal variable concentration
1717:
1718: case HECCER_MOP_INTERNALNERNST:
1719: {
1720: //- go to next operator
1721:
1722: struct MopsInternalNernst *pmops = (struct MopsInternalNernst *)piMop;
1723:
1724: piMop = (int *)&pmops[1];
1725:
1726: //- go to next type specific data
1727:
1728: struct MatsInternalNernst * pnernst = (struct MatsInternalNernst *)pvMats;
1729:
1730: pvMats = (void *)&((struct MatsInternalNernst *)pvMats)[1];
1731:
1732: //- get index of internal concentration
1733:
1734: int iInternal = pmops->uInternal.iMat;
1735:
1736: if (iInternal != -1)
1737: {
1738: //- get solved dependency
1739:
1740: double *pdConcentration = (double *)pheccer->vm.ppvMatsIndex[iInternal];
1741:
1742: //- store solved internal concentration
1743:
1744: pmops->uInternal.pdValue = pdConcentration;
1745: }
1746: else
1747: {
1748: pmops->uInternal.pdValue = NULL;
1749: }
1750:
1751: break;
1752: }
1753:
1754: //- for single channel initialization
1755:
1756: case HECCER_MOP_INITIALIZECHANNEL:
1757: {
1758: //- go to next operator
1759:
1760: struct MopsInitializeChannel *pmops = (struct MopsInitializeChannel *)piMop;
1761:
1762: piMop = (int *)&pmops[1];
1763:
1764: break;
1765: }
1766:
1767: //- for single channel initialization with variable reversal potential
1768:
1769: case HECCER_MOP_INITIALIZECHANNELEREV:
1770: {
1771: //- go to next operator
1772:
1773: struct MopsInitializeChannelErev *pmops = (struct MopsInitializeChannelErev *)piMop;
1774:
1775: piMop = (int *)&pmops[1];
1776:
1777: //- get index of reversal potential
1778:
1779: int iReversalPotential = pmops->uReversalPotential.iMat;
1780:
1781: if (iReversalPotential != -1)
1782: {
1783: //- get solved dependency
1784:
1785: double *pdNernst = (double *)pheccer->vm.ppvMatsIndex[iReversalPotential];
1786:
1787: //- store solved nernst potential
1788:
1789: pmops->uReversalPotential.pdValue = pdNernst;
1790: }
1791: else
1792: {
1793: pmops->uReversalPotential.pdValue = NULL;
1794: }
1795:
1796: break;
1797: }
1798:
1799: //- to compute a channel's conductance
1800:
1801: case HECCER_MOP_STORECHANNELCONDUCTANCE:
1802: {
1803: //- go to next operator
1804:
1805: piMop = &piMop[1];
1806:
1807: //- store the current conductance
1808:
1809: struct MatsStoreChannelConductance * pmats = (struct MatsStoreChannelConductance *)pvMats;
1810:
1811: //- go to next type specific data
1812:
1813: pvMats = (void *)&pmats[1];
1814:
1815: break;
1816: }
1817:
1818: //- for a new membrane potential
1819:
1820: case HECCER_MOP_LOADVOLTAGETABLE:
1821: {
1822: //- go to next operator
1823:
1824: struct MopsVoltageTableDependency *pmops = (struct MopsVoltageTableDependency *)piMop;
1825:
1826: piMop = (int *)&pmops[1];
1827:
1828: break;
1829: }
1830:
1831: //- for a conceptual gate (HH alike, with powers)
1832:
1833: case HECCER_MOP_CONCEPTGATE:
1834: {
1835: //- go to next operator
1836:
1837: struct MopsSingleGateConcept *pmops = (struct MopsSingleGateConcept *)piMop;
1838:
1839: piMop = (int *)&pmops[1];
1840:
1841: //- go to next type specific data
1842:
1843: struct MatsSingleGateConcept * pmats = (struct MatsSingleGateConcept *)pvMats;
1844:
1845: pvMats = (void *)&pmats[1];
1846:
1847: //- get possibly solved dependence
1848:
1849: int iMatsActivator = pmops->uState.iMat;
1850:
1851: if (iMatsActivator != -1)
1852: {
1853: //- get solved dependency
1854:
1855: double *pdMatsActivator = (double *)pheccer->vm.ppvMatsIndex[iMatsActivator];
1856:
1857: //- store solved dependency
1858:
1859: pmops->uState.pdValue = pdMatsActivator;
1860: }
1861: else
1862: {
1863: pmops->uState.pdValue = NULL;
1864: }
1865:
1866: break;
1867: }
1868:
1869: //- register single channel contribution
1870:
1871: case HECCER_MOP_REGISTERCHANNELCURRENT:
1872: {
1873: //- go to next operator
1874:
1875: piMop = &piMop[1];
1876:
1877: break;
1878: }
1879:
1880: //- add channel contribution to compartmental currents
1881:
1882: case HECCER_MOP_UPDATECOMPARTMENTCURRENT:
1883: {
1884: //- go to next operator
1885:
1886: struct MopsUpdateCompartmentCurrent *pmops = (struct MopsUpdateCompartmentCurrent *)piMop;
1887:
1888: piMop = (int *)&pmops[1];
1889:
1890: break;
1891: }
1892:
1893: //- compute current for single channel
1894:
1895: case HECCER_MOP_STORESINGLECHANNELCURRENT:
1896: {
1897: //- go to next operator
1898:
1899: struct MopsStoreSingleChannelCurrent *pmops = (struct MopsStoreSingleChannelCurrent *)piMop;
1900:
1901: piMop = (int *)&pmops[1];
1902:
1903: //- go to next type specific data
1904:
1905: struct MatsStoreSingleChannelCurrent * pmats = (struct MatsStoreSingleChannelCurrent *)pvMats;
1906:
1907: pvMats = (void *)&pmats[1];
1908:
1909: break;
1910: }
1911:
1912: //- compute exponential decay, mostly an ion concentration
1913:
1914: case HECCER_MOP_EXPONENTIALDECAY:
1915: {
1916: //- go to next operator
1917:
1918: struct MopsExponentialDecay *pmops = (struct MopsExponentialDecay *)piMop;
1919:
1920: piMop = (int *)&pmops[1];
1921:
1922: //- go to next type specific data
1923:
1924: struct MatsExponentialDecay * pmats = (struct MatsExponentialDecay *)pvMats;
1925:
1926: pvMats = (void *)&pmats[1];
1927:
1928: //- get possibly solved external flux contributions
1929:
1930: int i;
1931:
1932: for (i = 0 ; i < EXPONENTIALDECAY_CONTRIBUTORS ; i++)
1933: {
1934: int iExternal = pmops->puExternal[i].iMat;
1935:
1936: if (iExternal != -1)
1937: {
1938: int iOffset = 0;
1939:
1940: //- if individual channel currents are enabled
1941:
1942: if (pheccer->ho.iOptions & HECCER_OPTION_ENABLE_INDIVIDUAL_CURRENTS)
1943: {
1944: //- there is a mat entry in between, subtract one from the index
1945:
1946: iOffset = -1;
1947: }
1948:
1949: //- get solved dependency
1950:
1951: double *pdFlux = (double *)pheccer->vm.ppvMatsIndex[iExternal + iOffset];
1952:
1953: //- store solved external flux contribution
1954:
1955: pmops->puExternal[i].pdValue = pdFlux;
1956: }
1957: else
1958: {
1959: pmops->puExternal[i].pdValue = NULL;
1960: }
1961: }
1962:
1963: break;
1964: }
1965:
1966: //- register current contribution to a pool
1967:
1968: case HECCER_MOP_FLUXPOOL:
1969: {
1970: //- go to next operator
1971:
1972: struct MopsFluxPool *pmops = (struct MopsFluxPool *)piMop;
1973:
1974: piMop = (int *)&pmops[1];
1975:
1976: //- go to next type specific data
1977:
1978: struct MatsFluxPool * pmats = (struct MatsFluxPool *)pvMats;
1979:
1980: pvMats = (void *)&pmats[1];
1981:
1982: break;
1983: }
1984:
1985: //- for an event generator
1986:
1987: case HECCER_MOP_EVENTGENERATE:
1988: {
1989: //- go to next operator
1990:
1991: struct MopsEventGenerate *pmops = (struct MopsEventGenerate *)piMop;
1992:
1993: piMop = (int *)&pmops[1];
1994:
1995: //- go to next type specific data
1996:
1997: struct MatsEventGenerate * pmats = (struct MatsEventGenerate *)pvMats;
1998:
1999: pvMats = (void *)&pmats[1];
2000:
2001: //- get source for comparison
2002:
2003: /* int i; */
2004:
2005: /* for (i = 0 ; i < EVENT_SOURCES ; i++) */
2006: /* { */
2007: int iSource = pmops->uSource.iMat;
2008:
2009: //- if should be current membrane potential
2010:
2011: if (iSource == INT_MAX)
2012: {
2013: //- set pointer value to a sentinel value
2014:
2015: //! this will generate a warning on some architectures
2016:
2017: pmops->uSource.pdValue = (double *)-1;
2018: }
2019:
2020: //- if linked to a mechanism value
2021:
2022: else if (iSource != -1)
2023: {
2024: //- get solved dependency
2025:
2026: double *pdSource = (double *)pheccer->vm.ppvMatsIndex[iSource];
2027:
2028: //- store solved external flux contribution
2029:
2030: pmops->uSource.pdValue = pdSource;
2031: }
2032: else
2033: {
2034: pmops->uSource.pdValue = NULL;
2035: }
2036: /* } */
2037:
2038: break;
2039: }
2040:
2041: //- for a reset operation
2042:
2043: case HECCER_MOP_RESET:
2044: {
2045: //- go to next operator
2046:
2047: struct MopsReset *pmops = (struct MopsReset *)piMop;
2048:
2049: piMop = (int *)&pmops[1];
2050:
2051: break;
2052: }
2053:
2054: default:
2055: {
2056: //t HeccerError(number, message, varargs);
2057:
2058: HeccerError
2059: (pheccer,
2060: NULL,
2061: "unknown mechanism operation (%i)",
2062: piMop[0]);
2063:
2064: //! the best we can do is advance the pointer with one
2065:
2066: piMop = &piMop[1];
2067:
2068: break;
2069: }
2070: }
2071: }
2072: }
2073:
2074: //- check sanity of operator array
2075:
2076: if (piMop[0] != HECCER_MOP_FINISH)
2077: {
2078: //t add something like HeccerError(number, message, varargs);
2079:
2080: HeccerError
2081: (pheccer,
2082: NULL,
2083: "piMop[0] is %i, should be %i",
2084: piMop[0],
2085: HECCER_MOP_FINISH);
2086:
2087: //t depending on options, bail out
2088:
2089: //t set status : illegal mop hecc.
2090: }
2091:
2092: //- return result
2093:
2094: return(iResult);
2095: }
2096:
2097:
2098: /// **************************************************************************
2099: ///
2100: /// SHORT: HeccerMechanismReadDoubleFile()
2101: ///
2102: /// ARGS.:
2103: ///
2104: /// pheccer....: a heccer.
2105: /// pcFilename.: filename to read.
2106: /// ppdValues..: values that have been read, NULL for failure.
2107: ///
2108: /// RTN..: int
2109: ///
2110: /// Number of doubles read, -1 for failure.
2111: ///
2112: /// ppdValues..: values that have been read, NULL for failure.
2113: ///
2114: /// DESCR: Read a file with doubles, -1 terminated.
2115: ///
2116: /// **************************************************************************
2117:
2118: static
2119: int
2120: HeccerMechanismReadDoubleFile
2121: (struct Heccer *pheccer, char *pcFilename, double **ppdValues)
2122: {
2123: //- set default result: failure
2124:
2125: int iResult = -1;
2126:
2127: //- open file
2128:
2129: FILE *pfile = fopen(pcFilename, "r");
2130:
2131: if (!pfile)
2132: {
2133: return(iResult);
2134: }
2135:
2136: //- allocate result
2137:
2138: int iAllocated = 100;
2139:
2140: double *pdValues = (double *)calloc(iAllocated, sizeof(double));
2141:
2142: if (!pdValues)
2143: {
2144: fclose(pfile);
2145:
2146: return(iResult);
2147: }
2148:
2149: //- go through the file
2150:
2151: int iEOF = 0;
2152:
2153: int iDoubles = 0;
2154:
2155: while (!iEOF)
2156: {
2157: #define ELEMENT_NAME_SIZE 100
2158:
2159: //- read a record
2160:
2161: char pc[ELEMENT_NAME_SIZE * 10];
2162:
2163: if (fgets(pc, ELEMENT_NAME_SIZE * 10, pfile))
2164: {
2165: //- if not an element of a yaml array
2166:
2167: //! hardcoded indentation, 4 spaces required
2168:
2169: if ((pc[0] == ' '
2170: && pc[1] == ' '
2171: && pc[2] == ' '
2172: && pc[3] == ' '
2173: && pc[4] == '-'
2174: && pc[5] == ' '
2175: && pc[6] != ' ')
2176: ||
2177: (pc[0] == ' '
2178: && pc[1] == ' '
2179: && pc[2] == '-'
2180: && pc[3] == ' '
2181: && pc[4] != ' '))
2182: {
2183: int iScanned = sscanf(pc, " - %lf\n", &pdValues[iDoubles]);
2184:
2185: //- if number of scanned values is ok
2186:
2187: if (iScanned == 1)
2188: {
2189: //- next record
2190:
2191: iDoubles++;
2192:
2193: //- check for reallocation need
2194:
2195: if (iDoubles >= 100)
2196: {
2197: iAllocated += 100;
2198:
2199: pdValues = (double *)realloc(pdValues, iAllocated * sizeof(double));
2200:
2201: if (!pdValues)
2202: {
2203: break;
2204: }
2205: }
2206: }
2207: else
2208: {
2209: HeccerError
2210: (pheccer,
2211: NULL,
2212: "parse failure for HeccerMechanismReadDoubleFile(), record %i, scanned %i items (instead of 1)\n",
2213: iDoubles,
2214: iScanned);
2215:
2216: break;
2217: }
2218: }
2219: }
2220:
2221: //- or
2222:
2223: else
2224: {
2225: //- end of file
2226:
2227: iEOF = 1;
2228: }
2229: }
2230:
2231: //- set result
2232:
2233: if (ppdValues)
2234: {
2235: *ppdValues = pdValues;
2236: }
2237:
2238: //- close
2239:
2240: fclose(pfile);
2241:
2242: //- return result
2243:
2244: iResult = iDoubles;
2245:
2246: return(iResult);
2247: }
2248:
2249:
2250: /// **************************************************************************
2251: ///
2252: /// SHORT: HeccerMechanismSolveCN()
2253: ///
2254: /// ARGS.:
2255: ///
2256: /// pheccer...: a heccer.
2257: ///
2258: /// RTN..: int
2259: ///
2260: /// success of operation.
2261: ///
2262: /// DESCR: Perform the mechanisms operators once.
2263: ///
2264: /// **************************************************************************
2265:
2266: int HeccerMechanismSolveCN(struct Heccer *pheccer)
2267: {
2268: //- set default result : ok
2269:
2270: int iResult = TRUE;
2271:
2272: int *piMop = (int *)pheccer->vm.pvMops;
2273:
2274: //- get access to compartment variables and results
2275:
2276: double *pdVm = &pheccer->vm.pdVms[0];
2277:
2278: void *pvMats = pheccer->vm.pvMats;
2279:
2280: double *pdResults = &pheccer->vm.pdResults[0];
2281:
2282: //- loop over mechanism operators
2283:
2284: while (piMop[0] == HECCER_MOP_COMPARTMENT)
2285: {
2286: //- go to next operator
2287:
2288: piMop = &piMop[1];
2289:
2290: //- get compartment mechanism data
2291:
2292: struct MatsCompartment *pmatsc
2293: = (struct MatsCompartment *)pvMats;
2294:
2295: pvMats = (void *)&((struct MatsCompartment *)pvMats)[1];
2296:
2297: //- get membrane potential
2298:
2299: double dVm = pdVm[0];
2300:
2301: //- initialize current with leak and injected
2302:
2303: double dCompartmentCurrents = pmatsc->dLeak + pmatsc->dInjected;
2304:
2305: //- initial total conductance is zero
2306:
2307: double dConductances = 0;
2308:
2309: //- single channel contribution is zero
2310:
2311: double dChannelConductance = 0.0;
2312:
2313: double dReversalPotential = 0.0;
2314:
2315: //- loop over mechanism operators
2316:
2317: while (piMop[0] > HECCER_MOP_COMPARTMENT_BARRIER)
2318: {
2319: //v single channel current (if applicable)
2320:
2321: double dSingleChannelCurrent;
2322:
2323: //- look at next mechanism
2324:
2325: switch (piMop[0])
2326: {
2327: //- for aggregate current
2328:
2329: case HECCER_MOP_AGGREGATECURRENT:
2330: {
2331: //- go to next operator
2332:
2333: struct MopsAggregateCurrent *pmops = (struct MopsAggregateCurrent *)piMop;
2334:
2335: piMop = (int *)&pmops[1];
2336:
2337: pheccer->vm.pdAggregators[pmops->iIndex] += dSingleChannelCurrent;
2338:
2339: break;
2340: }
2341:
2342: //- for a call out
2343:
2344: case HECCER_MOP_CALLOUT:
2345: {
2346: //- go to next operator
2347:
2348: piMop = &piMop[1];
2349:
2350: //- go to next type specific data
2351:
2352: struct MatsCallout * pcall = (struct MatsCallout *)pvMats;
2353:
2354: pvMats = (void *)&((struct MatsCallout *)pvMats)[1];
2355:
2356: //- get function pointer
2357:
2358: struct Callout * pco = (struct Callout *)pcall->pco;
2359:
2360: ExternalFunction * pef = pco->pef;
2361:
2362: //- fill in internal results
2363:
2364: struct InternalResults * pir = pco->pir;
2365:
2366: pir->dVm = dVm;
2367:
2368: //- call the function
2369:
2370: struct ExternalResults * per = pco->per;
2371:
2372: int iResult = (*pef)(pco, pheccer, pir, per);
2373:
2374: //- handle external results
2375:
2376: dConductances += per->dConductance;
2377:
2378: dCompartmentCurrents += per->dCurrent;
2379:
2380: //t check signaling
2381:
2382: break;
2383: }
2384:
2385: //- for a spring mass equation
2386:
2387: case HECCER_MOP_SPRINGMASS:
2388: {
2389: //- go to next operator
2390:
2391: struct MopsSpringMass *pmops = (struct MopsSpringMass *)piMop;
2392:
2393: piMop = (int *)&pmops[1];
2394:
2395: //- go to next type specific data
2396:
2397: struct MatsSpringMass *pmats = (struct MatsSpringMass *)pvMats;
2398:
2399: pvMats = (void *)&((struct MatsSpringMass *)pvMats)[1];
2400:
2401: //- get pointer to table
2402:
2403: int iTable = pmops->iTable;
2404:
2405: struct HeccerTabulatedSpringMass *phtsm = &pheccer->tsmt.phtsm[iTable];
2406:
2407: //- if the firing frequence has been set
2408:
2409: if (pmops->dFrequency > 0)
2410: {
2411: //- generate random number
2412:
2413: int iRandom = RANDOM;
2414:
2415: //- check generated number with firing frequency
2416:
2417: if (iRandom < RAND_MAX * pmops->dFrequency)
2418: {
2419: //- add one to the activation, and apply decay
2420:
2421: //! weight not normalized to the time step
2422:
2423: pmats->dX += phtsm->dX1;
2424:
2425: /* fprintf(stdout, "Firing %i, %g\n", (int)(pheccer->dTime / pheccer->dStep), pheccer->dTime); */
2426:
2427: }
2428: }
2429:
2430: //- if there is an event time table
2431:
2432: if (pmops->pdEvents)
2433: {
2434: //- while the next event fires
2435:
2436: while (pmops->pdEvents[pmops->iEvent] < pheccer->dTime)
2437: {
2438: //- add one to the activation, and apply decay
2439:
2440: //! fixed weight of 1, normalized to the time step
2441:
2442: pmats->dX += phtsm->dX1;
2443:
2444: //- go to the next event
2445:
2446: //t should be a variable ?
2447:
2448: pmops->iEvent++;
2449: }
2450: }
2451:
2452: //- if there is an incoming event from the event distributor
2453:
2454: if (pmats->dNextEvent != -1.0
2455: && pmats->dNextEvent < pheccer->dTime
2456:
2457: //! target index for this object
2458:
2459: && pmops->iDiscreteTarget != -1)
2460: {
2461: //- translate incoming events to their activation
2462:
2463: double dActivation = HeccerEventReceive(pheccer, pmops->iDiscreteTarget);
2464:
2465: if (dActivation != FLT_MAX)
2466: {
2467: //- add the activation of (possibly multiple) events to the channel activation
2468:
2469: pmats->dX += dActivation * phtsm->dX1;
2470: }
2471: else
2472: {
2473: //t HeccerError(number, message, varargs);
2474:
2475: HeccerError
2476: (pheccer,
2477: NULL,
2478: "event reception failed for (%i) at time %g\n",
2479: pmops->iDiscreteTarget,
2480: pheccer->dTime);
2481:
2482: return(0);
2483: }
2484:
2485: //- reset the incoming event time
2486:
2487: pmats->dNextEvent = -1.0;
2488: }
2489:
2490: //- compute channel activation
2491:
2492: pmats->dX = pmats->dX * phtsm->dX2;
2493:
2494: pmats->dY = pmats->dX * phtsm->dY1 + pmats->dY * phtsm->dY2;
2495:
2496: //- compute channel conductance
2497:
2498: dChannelConductance *= pmats->dY;
2499:
2500: break;
2501: }
2502:
2503: //- for a nernst operation with internal variable concentration
2504:
2505: case HECCER_MOP_INTERNALNERNST:
2506: {
2507: //- go to next operator
2508:
2509: struct MopsInternalNernst *pmops = (struct MopsInternalNernst *)piMop;
2510:
2511: piMop = (int *)&pmops[1];
2512:
2513: //- go to next type specific data
2514:
2515: struct MatsInternalNernst * pnernst = (struct MatsInternalNernst *)pvMats;
2516:
2517: pvMats = (void *)&((struct MatsInternalNernst *)pvMats)[1];
2518:
2519: //- fetch variables
2520:
2521: double dExternal = pmops->dExternal;
2522:
2523: double dConstant = pmops->dConstant;
2524:
2525: double dInternal = 0.0;
2526:
2527: if (pmops->uInternal.pdValue)
2528: {
2529: dInternal = *pmops->uInternal.pdValue;
2530: }
2531:
2532: //- computer nernst potential
2533:
2534: //t produces invalid things if pdInternal is not set.
2535:
2536: double dPotential = log(dExternal / dInternal) * dConstant;
2537:
2538: //- store result
2539:
2540: pnernst->dPotential = dPotential;
2541:
2542: break;
2543: }
2544:
2545: //- for single channel initialization
2546:
2547: case HECCER_MOP_INITIALIZECHANNEL:
2548: {
2549: //- go to next operator
2550:
2551: struct MopsInitializeChannel *pmops = (struct MopsInitializeChannel *)piMop;
2552:
2553: piMop = (int *)&pmops[1];
2554:
2555: //- load channel variables
2556:
2557: dChannelConductance = pmops->dMaximalConductance;
2558:
2559: dReversalPotential = pmops->dReversalPotential;
2560:
2561: break;
2562: }
2563:
2564: //- for single channel initialization with variable reversal potential
2565:
2566: case HECCER_MOP_INITIALIZECHANNELEREV:
2567: {
2568: //- go to next operator
2569:
2570: struct MopsInitializeChannelErev *pmops = (struct MopsInitializeChannelErev *)piMop;
2571:
2572: piMop = (int *)&pmops[1];
2573:
2574: //- load channel variables
2575:
2576: dChannelConductance = pmops->dMaximalConductance;
2577:
2578: dReversalPotential = *pmops->uReversalPotential.pdValue;
2579:
2580: break;
2581: }
2582:
2583: //- to compute a channel's conductance
2584:
2585: case HECCER_MOP_STORECHANNELCONDUCTANCE:
2586: {
2587: //- go to next operator
2588:
2589: piMop = &piMop[1];
2590:
2591: //- store the current conductance
2592:
2593: struct MatsStoreChannelConductance * pmats = (struct MatsStoreChannelConductance *)pvMats;
2594:
2595: pmats->dChannelConductance = dChannelConductance;
2596:
2597: //- go to next type specific data
2598:
2599: pvMats = (void *)&pmats[1];
2600:
2601: break;
2602: }
2603:
2604: //- for a new membrane potential
2605:
2606: case HECCER_MOP_LOADVOLTAGETABLE:
2607: {
2608: //- go to next operator
2609:
2610: struct MopsVoltageTableDependency *pmops = (struct MopsVoltageTableDependency *)piMop;
2611:
2612: piMop = (int *)&pmops[1];
2613:
2614: //t this is a nop for the moment, but when table
2615: //t rearrangements get in, this should load the
2616: //t table pointed to by the current membrane potential.
2617:
2618: break;
2619: }
2620:
2621: //- for a conceptual gate (HH alike, with powers)
2622:
2623: case HECCER_MOP_CONCEPTGATE:
2624: {
2625: //- go to next operator
2626:
2627: struct MopsSingleGateConcept *pmops = (struct MopsSingleGateConcept *)piMop;
2628:
2629: piMop = (int *)&pmops[1];
2630:
2631: //- go to next type specific data
2632:
2633: struct MatsSingleGateConcept * pmats = (struct MatsSingleGateConcept *)pvMats;
2634:
2635: pvMats = (void *)&pmats[1];
2636:
2637: //- get type specific constants and variables
2638:
2639: int iPower = pmops->iPower;
2640:
2641: int iTable = pmops->iTableIndex;
2642:
2643: double *pdState = pmops->uState.pdValue;
2644:
2645: double dActivation = pmats->dActivation;
2646:
2647: //- fetch tables
2648:
2649: //t table rearrangements
2650:
2651: struct HeccerTabulatedGate *phtg = &pheccer->tgt.phtg[iTable];
2652:
2653: double *pdA = phtg->pdA;
2654: double *pdB = phtg->pdB;
2655:
2656: double dState;
2657:
2658: //- for a concentration dependent gate
2659:
2660: if (pdState)
2661: {
2662: //- state is coming from a solved mechanism variable
2663:
2664: dState = *pdState;
2665: }
2666:
2667: //- else is a membrane potential dependent gate
2668:
2669: else
2670: {
2671: //- state is membrane potential
2672:
2673: //t move this to load membrane potential
2674: //t need LOADVOLTAGETABLE or LOADVOLTAGEINDEX
2675:
2676: dState = dVm;
2677: }
2678:
2679: //- discretize and offset the state
2680:
2681: int iIndex = (dState - phtg->hi.dStart) / phtg->hi.dStep;
2682:
2683: if (iIndex < 0)
2684: {
2685: iIndex = 0;
2686: }
2687: else if (iIndex >= phtg->iEntries)
2688: {
2689: iIndex = phtg->iEntries - 1;
2690: }
2691:
2692: //- fetch A and B gate rates
2693:
2694: double dA = pdA[iIndex];
2695: double dB = pdB[iIndex];
2696:
2697: //t move to channel rearrangements
2698:
2699: dB = 1.0 + pheccer->dStep / 2.0 * dB;
2700:
2701: dA = pheccer->dStep * dA;
2702:
2703: //- compute gate activation state
2704:
2705: dActivation = dA / dB + dActivation * 2.0 / dB - dActivation;
2706:
2707: //- and store it for the next cycle
2708:
2709: pmats->dActivation = dActivation;
2710:
2711: //- apply gate power, multiply with conductance (note:
2712: //- also changes units)
2713:
2714: if (iPower == 1)
2715: {
2716: dChannelConductance = dChannelConductance * dActivation;
2717: }
2718: else if (iPower == 2)
2719: {
2720: dChannelConductance *= dActivation * dActivation;
2721: }
2722: else if (iPower == 3)
2723: {
2724: dChannelConductance *= dActivation * dActivation * dActivation;
2725: }
2726: else if (iPower == 4)
2727: {
2728: dActivation *= dActivation;
2729: dChannelConductance *= dActivation * dActivation;
2730: }
2731: else if (iPower == 5)
2732: {
2733: dChannelConductance *= dActivation;
2734: dActivation *= dActivation;
2735: dActivation *= dActivation;
2736: dChannelConductance *= dActivation;
2737: }
2738: else if (iPower == 6)
2739: {
2740: dActivation *= dActivation;
2741: dChannelConductance *= dActivation;
2742: dActivation *= dActivation;
2743: dChannelConductance *= dActivation;
2744: }
2745: else
2746: {
2747: //t HeccerError(number, message, varargs);
2748:
2749: HeccerError
2750: (pheccer,
2751: NULL,
2752: "invalid gate power (%i)",
2753: iPower);
2754:
2755: *(int *)0 = 0;
2756: }
2757:
2758: break;
2759: }
2760:
2761: //- register single channel contribution
2762:
2763: case HECCER_MOP_REGISTERCHANNELCURRENT:
2764: {
2765: //- go to next operator
2766:
2767: piMop = &piMop[1];
2768:
2769: //- compute single channel current
2770:
2771: dSingleChannelCurrent = dChannelConductance * (dReversalPotential - dVm);
2772:
2773: break;
2774: }
2775:
2776: //- add channel contribution to compartmental currents
2777:
2778: case HECCER_MOP_UPDATECOMPARTMENTCURRENT:
2779: {
2780: //- go to next operator
2781:
2782: struct MopsUpdateCompartmentCurrent *pmops = (struct MopsUpdateCompartmentCurrent *)piMop;
2783:
2784: piMop = (int *)&pmops[1];
2785:
2786: //- compute conductance for matrix left side
2787:
2788: dConductances += dChannelConductance;
2789:
2790: //- compute current for matrix right side
2791:
2792: dCompartmentCurrents += dChannelConductance * dReversalPotential;
2793:
2794: break;
2795: }
2796:
2797: //- compute current for single channel
2798:
2799: case HECCER_MOP_STORESINGLECHANNELCURRENT:
2800: {
2801: //- go to next operator
2802:
2803: struct MopsStoreSingleChannelCurrent *pmops = (struct MopsStoreSingleChannelCurrent *)piMop;
2804:
2805: piMop = (int *)&pmops[1];
2806:
2807: //- go to next type specific data
2808:
2809: struct MatsStoreSingleChannelCurrent * pmats = (struct MatsStoreSingleChannelCurrent *)pvMats;
2810:
2811: pvMats = (void *)&pmats[1];
2812:
2813: //- store the current, use appropriate indexing
2814:
2815: pmats->dConductance = dChannelConductance;
2816: pmats->dCurrent = dSingleChannelCurrent;
2817:
2818: break;
2819: }
2820:
2821: //- compute exponential decay, mostly an ion concentration
2822:
2823: case HECCER_MOP_EXPONENTIALDECAY:
2824: {
2825: //- go to next operator
2826:
2827: struct MopsExponentialDecay *pmops = (struct MopsExponentialDecay *)piMop;
2828:
2829: piMop = (int *)&pmops[1];
2830:
2831: //- go to next type specific data
2832:
2833: struct MatsExponentialDecay * pmats = (struct MatsExponentialDecay *)pvMats;
2834:
2835: pvMats = (void *)&pmats[1];
2836:
2837: double dBeta = pmops->dBeta;
2838:
2839: double dSteadyState = pmops->dSteadyState;
2840:
2841: double dTau = pmops->dTau;
2842:
2843: double dState = pmats->dState;
2844:
2845: //- fetch external contribution
2846:
2847: //! so after the channel mat comes the flux mat. This
2848: //! is a hack ... which could be solved, by having an
2849: //! intermediary that linearizes and specifies
2850: //! _all_ the equations (so including fluxes).
2851:
2852: //t In other words, the test 'if (pcac->iPool != -1)'
2853: //t should be removed and the removal should be
2854: //t counterbalanced by the neccessary specifications
2855: //t for a flux in the intermediary. Then, these
2856: //t specifications are compiled to HECCER_MOP_FLUXPOOL
2857: //t byte code.
2858:
2859: //t needs a careful check again.
2860:
2861: double dExternal = 0.0;
2862:
2863: int i;
2864:
2865: for (i = 0 ; i < EXPONENTIALDECAY_CONTRIBUTORS ; i++)
2866: {
2867: double *pdExternal = pmops->puExternal[i].pdValue;
2868:
2869: if (pdExternal)
2870: {
2871: dExternal += *pdExternal;
2872: }
2873: }
2874:
2875: //- exponential decay with possibly external influx
2876:
2877: /* dState = dSteadyState + ((dState - dSteadyState) * (2.0 - dTau) + (dExternal * dBeta)) / dTau; */
2878:
2879: dState = (dSteadyState * dTau * 2.0 + dState * 2.0 + dExternal * dBeta - dSteadyState * 2.0) / dTau - dState;
2880:
2881: //- save state
2882:
2883: pmats->dState = dState;
2884:
2885: break;
2886: }
2887:
2888: //- register current contribution to a pool
2889:
2890: case HECCER_MOP_FLUXPOOL:
2891: {
2892: //- go to next operator
2893:
2894: struct MopsFluxPool *pmops = (struct MopsFluxPool *)piMop;
2895:
2896: piMop = (int *)&pmops[1];
2897:
2898: //- go to next type specific data
2899:
2900: struct MatsFluxPool * pmats = (struct MatsFluxPool *)pvMats;
2901:
2902: pvMats = (void *)&pmats[1];
2903:
2904: //- register contribution
2905:
2906: pmats->dFlux = dSingleChannelCurrent;
2907:
2908: break;
2909: }
2910:
2911: //- for an event generator
2912:
2913: case HECCER_MOP_EVENTGENERATE:
2914: {
2915: //- go to next operator
2916:
2917: struct MopsEventGenerate *pmops = (struct MopsEventGenerate *)piMop;
2918:
2919: piMop = (int *)&pmops[1];
2920:
2921: //- go to next type specific data
2922:
2923: struct MatsEventGenerate * pmats = (struct MatsEventGenerate *)pvMats;
2924:
2925: pvMats = (void *)&pmats[1];
2926:
2927: //- if spiking
2928:
2929: if (pmats->dRefractory == FLT_MAX)
2930: {
2931: //- initialize refractory period
2932:
2933: pmats->dRefractory = pmops->dRefractoryReset - pheccer->dStep;
2934:
2935: //- set not spiking
2936:
2937: pmats->dSpike = 0.0;
2938: }
2939:
2940: //- if in refractory period
2941:
2942: else if (pmats->dRefractory > 0.0)
2943: {
2944: //- maintain refractory period since last event
2945:
2946: pmats->dRefractory -= pheccer->dStep;
2947: }
2948:
2949: //- if not in refractory period
2950:
2951: else
2952: {
2953: //- get source value
2954:
2955: double *pdSource = pmops->uSource.pdValue;
2956:
2957: if (pdSource)
2958: {
2959: double dSource;
2960:
2961: //- if using membrane potential
2962:
2963: if (pdSource == (double *)-1)
2964: {
2965: dSource = dVm;
2966: }
2967:
2968: //- if using a solved mechanism value
2969:
2970: else
2971: {
2972: dSource = *pdSource;
2973: }
2974:
2975: //- if source greater than threshold
2976:
2977: if (dSource > pmops->dThreshold)
2978: {
2979: //- register spiking
2980:
2981: pmats->dRefractory = FLT_MAX;
2982:
2983: pmats->dSpike = 1.0;
2984:
2985: //- generate event for associated targets
2986:
2987: int iResult = HeccerEventGenerate(pheccer, pmops->iTable);
2988:
2989: if (!iResult)
2990: {
2991: //t HeccerError(number, message, varargs);
2992:
2993: HeccerError
2994: (pheccer,
2995: NULL,
2996: "spike generation failed at time (%g)",
2997: pheccer->dTime);
2998: }
2999: }
3000: }
3001: }
3002:
3003: break;
3004: }
3005:
3006: //- for a reset operation
3007:
3008: case HECCER_MOP_RESET:
3009: {
3010: //- go to next operator
3011:
3012: struct MopsReset *pmops = (struct MopsReset *)piMop;
3013:
3014: piMop = (int *)&pmops[1];
3015:
3016: //- get reset value
3017:
3018: double dVmNew = pmops->dReset;
3019:
3020: //- set the membrane potential
3021:
3022: //t what if the user mixes this with active currents ?
3023: //t this is ok for simple integrate and fire,
3024: //t but for more complicated models I am not sure.
3025:
3026: pdVm[0] = dVmNew;
3027:
3028: break;
3029: }
3030:
3031: default:
3032: {
3033: //t HeccerError(number, message, varargs);
3034:
3035: HeccerError
3036: (pheccer,
3037: NULL,
3038: "unknown mechanism operation (%i)",
3039: piMop[0]);
3040:
3041: //! the best we can do is advance the pointer with one
3042:
3043: piMop = &piMop[1];
3044:
3045: break;
3046: }
3047: }
3048: }
3049:
3050: //- for single compartment neurons
3051:
3052: if (pheccer->vm.iVms == 1)
3053: {
3054: //- compute the membrane potential right here
3055:
3056: //t differentiation needed between CN and BE ?
3057:
3058: double dResult = ((dVm + dCompartmentCurrents * pmatsc->dCapacity)
3059: / (dConductances * pmatsc->dCapacity + pmatsc->dDiagonal));
3060:
3061: pdVm[0] = dResult + dResult - pdVm[0];
3062: }
3063:
3064: //- for multiple compartment neurons
3065:
3066: else
3067: {
3068: //- compute results for the morphology matrix
3069:
3070: //- right side
3071:
3072: pdResults[0] = dVm + dCompartmentCurrents * pmatsc->dCapacity;
3073:
3074: //- left side
3075:
3076: pdResults[1] = dConductances * pmatsc->dCapacity + pmatsc->dDiagonal;
3077:
3078: pdResults += 2;
3079: }
3080:
3081: //- go to next compartment
3082:
3083: pdVm++;
3084: }
3085:
3086: //- check sanity of operator array
3087:
3088: if (piMop[0] != HECCER_MOP_FINISH)
3089: {
3090: //t add something like HeccerError(number, message, varargs);
3091:
3092: HeccerError
3093: (pheccer,
3094: NULL,
3095: "piMop[0] is %i, should be %i",
3096: piMop[0],
3097: HECCER_MOP_FINISH);
3098:
3099: //t depending on options, bail out
3100:
3101: //t set status : illegal mop hecc.
3102: }
3103:
3104: //- return result
3105:
3106: return(iResult);
3107: }
3108:
3109:
3110: /// **************************************************************************
3111: ///
3112: /// SHORT: HeccerMechanismSort()
3113: ///
3114: /// ARGS.:
3115: ///
3116: /// pheccer...: a heccer.
3117: ///
3118: /// RTN..: int
3119: ///
3120: /// success of operation.
3121: ///
3122: /// DESCR: Sort mechanisms according to different characteristics.
3123: ///
3124: /// Mechanisms are sorted using a regular quicksort. A full order
3125: /// is enforced by looking at the serial identification number in
3126: /// the intermediary (as last resort for unordered mechanisms).
3127: ///
3128: /// NOTE.:
3129: ///
3130: /// Compartment leak is a current, so it is considered to be a
3131: /// mechanism.
3132: ///
3133: /// **************************************************************************
3134:
3135: int HeccerMechanismSort(struct Heccer *pheccer)
3136: {
3137: //- set default result : ok
3138:
3139: int iResult = TRUE;
3140:
3141: //- loop over all compartments
3142:
3143: int iCompartmentSchedule;
3144:
3145: for (iCompartmentSchedule = 0 ; iCompartmentSchedule < pheccer->inter.iCompartments ; iCompartmentSchedule++)
3146: {
3147: //- get model number for this compartment
3148:
3149: int iCompartmentModel = pheccer->indexers.md.piBackward[iCompartmentSchedule];
3150:
3151: //- loop over all the mechanisms for this compartment
3152:
3153: int iMechanismModel;
3154:
3155: for (iMechanismModel = pheccer->inter.piC2m[iCompartmentModel] ;
3156: iMechanismModel < pheccer->inter.piC2m[iCompartmentModel + 1] ;
3157: iMechanismModel++)
3158: {
3159: //t build array with mechanism references into the intermediary
3160:
3161: //t because of the structure of the intermediary, we do
3162: //t not need this, do we ?
3163: }
3164: }
3165:
3166: //t qsort mechanisms per compartment
3167:
3168: //t sort on Erev, see adaptive time step paper and hsolve implementation.
3169:
3170: //- return result
3171:
3172: return(iResult);
3173: }
3174:
3175:
3176:
Generated by Xrefactory version 2.0.14 on Thu Jul 24 22:41:20 2008