/* MemHH? Module */
/* */
/* patch file for the Hodgkin-Huxley membrane */
/* model, Vm goes 0 to 110 mV with AP */
/* */
/* Vrest is in mV */
/* Ek, Ena, and El are in mV */
/* Gk, Gna, and Gl are in mS (?) */
/* original paper: AL Hodgkin and AF Huxley, "A quantitative description */
/* of membrane current and its application to conductance and excitation */
/* in nerve," J Physiol, v.117, p.500, 1952 */
/* author of this computer code: John B. Pormann */
#include "CardioWave.h"
void SetPatches?_HH( domain_t tp, vector Vm, vector Q );
int GetF_HH( vector Vm, vector Qv, vector Fv, vector Fq, vector Av );
/* Hodgkin-Huxley patch */
typedef struct {
real m,n,h;
} HH_patch;
/* Auxliary variables - for ionic currents, etc. */
typedef struct {
real Ina, Ik, Il;
} HH_auxvar;
/* pseudo-globals */
int HH_PatchSize = sizeof(HH_patch)/sizeof(real);
int HH_AuxiliarySize = sizeof(HH_auxvar)/sizeof(real);
int HH_MemParamSize = 0;
byte HH_NodeType = ByteError?;
static char* RCSID = "$Id: MemHH.c,v 1.4 2005/02/17 20:57:04 john Exp $";
/* constants */
static real HH_Gna = 120.0;
static real HH_Gk = 36.0;
static real HH_Gl = 0.30;
static real HH_Ena = 115.0;
static real HH_Ek = -12.0;
static real HH_El = 10.6130;
static real HH_RestVoltage = 0.0;
/* patch OldRestPatch? = { 0.05293249, 0.31767914, 0.59612075 }; */
static HH_patch HH_RestPatch = { 0.0529380391521187335102105,
0.3177029287926248501960913,
0.5952926477586697462385246 };
static rword resources[] = {
{ "HH_Ek", 1001 },
{ "HH_El", 1002 },
{ "HH_Ena", 1003 },
{ "HH_Gk", 1004 },
{ "HH_Gl", 1005 },
{ "HH_Gna", 1006 },
{ "HH_IV", 1007 },
{ "HH_Node", 1100 },
{ "HH_Nodetype",1100 },
{ "HH_Patch", 1007 },
{ "HH_Type", 1100 },
{ "HH_Vr", 1008 },
{ "HH_Vrest", 1008 },
{ "showofs", 1200 },
{ "showoffset", 1200 },
{ "showoffsets",1200 },
{ "show_ofs", 1200 },
{ "show_offset",1200 },
{ "show_offsets",1200 },
{ NULL, 0 }
};
int InitMembrane?_HH( char** res ) {
int i,j,c;
int cmd;
real* iv;
real* p;
logical show_ofs = false;
i = 0;
while( res[i] != NULL ) {
cmd = FindCommand( resources, res[i] );
switch( cmd ) {
case 1001:
HH_Ek = GetRealValue( res[i] );
break;
case 1002:
HH_El = GetRealValue( res[i] );
break;
case 1003:
HH_Ena = GetRealValue( res[i] );
break;
case 1004:
HH_Gk = GetRealValue( res[i] );
break;
case 1005:
HH_Gl = GetRealValue( res[i] );
break;
case 1006:
HH_Gna = GetRealValue( res[i] );
break;
case 1007:
iv = GetRealArray( res[i] );
p = (real*)(&HH_RestPatch);
c = GetNumValues( res[i] );
if( c > HH_PatchSize ) {
c = HH_PatchSize;
}
for(j=0;j<c;j++) {
p[j] = iv[j];
}
break;
case 1008:
HH_RestVoltage = GetRealValue( res[i] );
break;
case 1100:
HH_NodeType = GetByteValue( res[i] );
break;
case 1200:
show_ofs = GetTFValue( res[i] );
break;
}
i++;
}
/* required!! */
if( HH_NodeType == ByteError? ) {
return( -1 );
}
PatchSize[HH_NodeType] = HH_PatchSize;
if( UseAuxvars ) {
AuxiliarySize[HH_NodeType] = HH_AuxiliarySize;
} else {
AuxiliarySize[HH_NodeType] = 0;
}
MemParamSize[HH_NodeType] = HH_MemParamSize;
FunctionTable[HH_NodeType] = GetF_HH;
SetpatchTable[HH_NodeType] = SetPatches?_HH;
if( (DebugLevel>0) and (SelfPE==0) ) {
printf("MemHH?: patchsize = 0.2lf %0.2lf %0.2lf %0.2lf %0.2lf %0.2lf\n",
HH_PatchSize,HH_Gna,HH_Gk,HH_Gl,
HH_Ena,HH_Ek,HH_El);
if( DebugLevel > 1 ) {
printf("MemHH?: RCSID: %s\n",RCSID);
}
}
if( show_ofs and (SelfPE==0) ) {
printf("HH patch offsets:\n");
printf(" M: %i\n", OffsetOf?(HH_patch,m) );
printf(" N: %i\n", OffsetOf?(HH_patch,n) );
printf(" H: %i\n", OffsetOf?(HH_patch,h) );
printf("HH auxvar offsets:\n");
printf(" Ina: %i\n", OffsetOf?(HH_auxvar,Ina) );
printf(" Ik: %i\n", OffsetOf?(HH_auxvar,Ik) );
printf(" Il: %i\n", OffsetOf?(HH_auxvar,Il) );
}
return( 0 );
}
void ExitMembrane?_HH( void ) {
return;
}
void SetPatches?_HH( domain_t tp, vector Vm, vector Q ) {
int i,j;
real* vm = Vm.data;
real* qp = Q.data;
real* rp = (real*)(&HH_RestPatch);
byte* nt = NodeType.data;
if( tp == Tissue ) {
for(i=0;i<TissueLocalSize;i++) {
if( nt[i] == HH_NodeType ) {
vm[i] = HH_RestVoltage;
for(j=0;j<HH_PatchSize;j++) {
qp[j] = rp[j];
}
}
qp += PatchSize[ nt[i] ];
}
}
return;
}
/* GetF() should fill in vectors Fv and Fq with the new currents/updates */
/* NOTE: this should be from the form: dV/dt + Fv = 0, dq/dt + Fq = 0 */
int GetF_HH( vector Vm, vector Qv, vector Fv, vector Fq, vector Av ) {
int i;
real* vmptr = Vm.data;
real* fv = Fv.data;
real* qm = Qv.data;
real* fq = Fq.data;
real* av = Av.data;
HH_patch* qp;
HH_patch* fp;
HH_auxvar* ap;
byte* nt = NodeType.data;
real vm,m,n,h;
real an,bn,am,bm,ah,bh;
real Ina,Ik,Il;
for(i=0;i<TissueLocalSize;i++) {
if( nt[i] == HH_NodeType ) {
/* map pointers to patch structures */
qp = (HH_patch*)qm;
fp = (HH_patch*)fq;
ap = (HH_auxvar*)av;
/* get local variables */
vm = vmptr[i];
m = qp->m;
n = qp->n;
h = qp->h;
/* calculate the alphas and betas */
am = (0.1*(25.0-vm))/(exp((25.0-vm)/10.0)-1.0);
bm = 4.0*exp(-vm/18.0);
an = (0.01*(10.0-vm))/(exp((10.0-vm)/10.0)-1.0);
bn = 0.125*exp(-vm/80.0);
ah = 0.07*exp(-vm/20.0);
bh = (1.0)/(exp((30.0-vm)/10.0)+1.0);
Ik = HH_Gk*(n*n*n*n)*( HH_Ek - vm );
Ina = HH_Gna*(m*m*m*h)*( HH_Ena - vm );
Il = HH_Gl*( HH_El - vm );
fv[i] += Ina + Ik + Il;
fp->m = am*(1.0-m) - bm*m;
fp->n = an*(1.0-n) - bn*n;
fp->h = ah*(1.0-h) - bh*h;
if( UseAuxvars ) {
ap->Ina = Ina;
ap->Ik = Ik;
ap->Il = Il;
}
} /* end-if */
/* go to next patch */
qm += PatchSize[ nt[i] ];
fq += PatchSize[ nt[i] ];
av += AuxiliarySize[ nt[i] ];
}
return( 0 );
}
