/* ghk.c Subba Shankar 8/25/91 * * This object calculates the current through a membrane, * the reversal potential, and the chord conductance * from the Goldman Hodgkin Katz equation (constant field equation). * * modified from nernst.c written by Matt Wilson * * The equation used is: * * Cin*exp(K * Vm) - Cout * Ik = p * F * K * Vm ------------------------ * exp(K * Vm) - 1 * * valency * F * where K = ------------- * R * (T + 273) * * F = Faraday's constant (internal) * * R = universal gas constant (internal) * * * Ek is the Nernst (or reversal) potential, and Gk is the chord * conductance, calculated from Ohm's Law. * * A full explanation of the constant field equn can be found in * 1) _From Neuron to Brain_, Kuffler, et al., Sinauer Assoc. * 2nd edition, (1984) p. 123, 130. * 2) _Electric Current Flow in Excitable Cells_, Jack, et al., * Oxford Press, (1983) p. 237. * Note that both of these references differ in their use of the * valence, and that the Kuffler equation would cause outward currents * to be negative. * * The current field can be used directly, but if the * object that you are passing this information to * does not understand current messages (such as compartment) * then pass the Ek and Gk values (e.g. in a CHANNEL mesg). * This will yield the equivalent current flow from Ohm's law. * * See below for current convention. * See ghk.doc for more thorough discussion. */ #include #include "ghk_ext.h" /* Define the current convention: * * 1 = outward current is positive * -1 = inward current is positive * * If you change this, please change the documentation * in ghk.doc and ghklib.g to reflect the change. */ #define CURRENT_CONV 1 /* define universal constants needed for * calculations. Taken from _Ionic Channels of * Excitable Membranes_ by Bertil Hille (1984) * Sinauer Associates, Sunderland MA, pg. 7. * Note that GAS_CONSTANT is usually called R. */ #define GAS_CONSTANT 8.314 /* (V * C)/(deg K * mol) */ #define FARADAY 9.648e4 /* C / mol */ #define ZERO_CELSIUS 273.16 /* deg */ #define R_OVER_F 8.6171458e-5 /* volt/deg */ #define F_OVER_R 1.1605e4 /* deg/volt */ #define F_TIMES_L_TO_M 9.648e7 /* FARADAY times conversion factor * from m^3 to liters */ /* define the message types that are accepted */ #define CIN 0 #define COUT 1 #define TEMP 2 #define VOLTAGE 3 #define PERMEABILITY 4 /* ** a single type of ion ** voltage_scale is in units of voltage conversion, ** should be 1 for input of volts ** conc_scale is in units of concentration conversion, ** should be 1 for input of molar (moles/liter) ** perm_scale should be 1 for input of m^3/sec ** current_scale is in units of current conversion, ** should be 1 for output of amps ** T has units of degrees Celsius ** ** In order to save time the common expression: ** [voltage_scale*valency*F]/[R*(temp in deg K)] ** is calculated and saved in ghk->constant */ void Ghk(ghk,action) register struct ghk_type *ghk; Action *action; { MsgIn *msg; /* following are used for common expression reduction */ double constant; double exponent, e_to_exponent; /* code: */ if(debug > 1){ ActionHeader("Ghk",ghk,action); } SELECT_ACTION(action){ case PROCESS: case RESET: /* first pick up all of our messages */ MSGLOOP(ghk,msg){ case CIN: /* intracellular ionic concentration */ ghk->Cin = MSGVALUE(msg,0); break; case COUT: /* extracellular ionic concentration */ ghk->Cout = MSGVALUE(msg,0); break; case TEMP: /* temperature in degrees Celsius */ ghk->T = MSGVALUE(msg,0); break; case VOLTAGE: /* transmembrane voltage */ ghk->Vm = MSGVALUE(msg,0); break; case PERMEABILITY: /* permeability, */ ghk->p = MSGVALUE(msg,0); break; } /* end of MSGLOOP */ /* Now calculate the output fields */ /* recalculate the constant */ constant = ghk->voltage_scale*F_OVER_R*ghk->valency/ (ghk->T + ZERO_CELSIUS); /* Calculate the reversal potential: ** Ek = ln([out]/[in])/constant ** constant = voltage_scale*(zF/RT) */ ghk->Ek = log(ghk->Cout/ghk->Cin) / constant; /* calculate the current: * 1) check for singularities near zero voltage * 2) calculate the common expressions * 3) calculate the current and the chord conductance */ exponent = constant*(ghk->Vm); /* unitless value, raised * to power e later */ e_to_exponent = exp( exponent ); if ( fabs(exponent) < .001 ) { /* exponent near zero, * calculate current * some other way */ /* take Taylor expansion of V'/[exp(V') - 1], where * V' = constant * Vm * First two terms should be enough this close to zero */ ghk->Ik = CURRENT_CONV * ghk->current_scale * ghk->perm_scale * ghk->p * F_TIMES_L_TO_M * ghk->conc_scale * ((ghk->Cin * e_to_exponent) - ghk->Cout) / (1 + 0.5 * constant * ghk->Vm); }/* end if fabs */ else { /* exponent far from zero, calculate directly */ /* note that F_TIMES_L_TO_M is Faradays constant and also * converts permeability from liters/sec to m^3/sec */ ghk->Ik = CURRENT_CONV * ghk->current_scale * ghk->perm_scale * ghk->p * F_TIMES_L_TO_M * constant * ghk->Vm * ghk->conc_scale * ((ghk->Cin * e_to_exponent) - ghk->Cout) / (e_to_exponent - 1.0); }/* end else */ /* Now calculate the chord conductance, but * check the denominator for a divide by zero. */ if ( ghk->voltage_scale * fabs(ghk->Vm - ghk->Ek) < 1e-12 ) { /* we are very close to the rest potential, so just set the * current and conductance to zero. */ ghk->Ik = ghk->Gk = 0.0; } /* if */ else { /* calculate in normal way */ ghk->Gk = ghk->Ik / (ghk->Vm - ghk->Ek); } /* else */ break; /* end of case PROCESS and RESET */ case CHECK: if(ghk->valency == 0){ ErrorMessage("Ghk", "Invalid ionic valency z.", ghk); } if(ghk->Cin < 0.0){ ErrorMessage("Ghk", "Invalid intracellular ionic concentration.", ghk); } if(ghk->Cout < 0.0){ ErrorMessage("Ghk", "Invalid extracellular ionic concentration.", ghk); } if(ghk->voltage_scale <= 0){ ErrorMessage("Ghk", "Invalid voltage scale parameter.", ghk); } if(ghk->current_scale <= 0){ ErrorMessage("Ghk", "Invalid current scale parameter.", ghk); } if(ghk->conc_scale <= 0){ ErrorMessage("Ghk", "Invalid concentration scale parameter.", ghk); } if(ghk->perm_scale <= 0){ ErrorMessage("Ghk", "Invalid permeability scale parameter.", ghk); } if(ghk->T + ZERO_CELSIUS <= 0){ ErrorMessage("Ghk", "Invalid temperature parameter.", ghk); } if(ghk->p < 0.0){ ErrorMessage("Ghk", "Invalid permeability parameter.", ghk); } break; /* end of case CHECK */ } /* end of SELECT_ACTION */ } /* end of function Ghk */