Logo Search packages:      
Sourcecode: matwrap version File versions  Download package

hodgkin_huxley.cxx

//
// This module contains the instantiation of various members of the
// Hodgkin_Huxley class.  We use an implicit Euler method to update the
// m, n, and h variables to compute the conductance.
//
#include "hodgkin_huxley.h"
#include <math.h>
#include <iostream.h>

//
// The various activation and inactivation functions.  These are tabulated.
// Note that we only need alpha+beta, not beta independently.  The code
// for these functions does not need to be very efficient, since they are
// only evaluated at the beginning and the results are tabulated.
//
// Alpha and beta are defined as:
//
// dm
// -- = (1-m) alpha(V) - m beta(V)
// dt
//
// and similarly for the other state variables.
// The constants here are those specified in the original Hodgkin-Huxley
// papers except that the sign of V has been reversed and the resting state
// is V=-65 mV instead of V=0.  Currently, you can't change the constants.
//
static Float m_alpha(Float V) {
  return fabs(V+40) < 0.00001 ? 1. : -0.1*(V+40.)/(exp(-0.1*(V+40.))-1.);
}
static Float m_alpha_plus_beta(Float V) {
  return m_alpha(V) + 4*exp(-(V+65.)/18.);
}

static Float h_alpha(Float V) {
  return 0.07*exp(-(V+65.)/20.);
}
static Float h_alpha_plus_beta(Float V) {
  return h_alpha(V) + 1/(exp(-(V+35.)/10.)+1);
}

static Float n_alpha(Float V) {
  return fabs(V+55.) < 0.00001 ? 0.1 : -0.01*(V+55.)/(exp(-0.1*(V+55.))-1);
}
static Float n_alpha_plus_beta(Float V) {
  return n_alpha(V) + 0.125*exp(-(V+65.)/80.);
}

//
// Tabulated versions of the above:
//

static const Float table_max = 50.; // Tables go up to 50 mV.
static const Float table_min = -100.; // and down to -100.
static const Float table_step = 5.; // From plotting, 5 mV gives very good
                        // interpolation.  We want this to be as big
                        // as possible so the tables are as small as
                        // possible so we don't flush cache.

Float m_alpha_table[int(1.00001*(table_max-table_min)/table_step)+1];
Float m_alpha_plus_beta_table[int(1.00001*(table_max-table_min)/table_step)+1];

Float h_alpha_table[int(1.00001*(table_max-table_min)/table_step)+1];
Float h_alpha_plus_beta_table[int(1.00001*(table_max-table_min)/table_step)+1];

Float n_alpha_table[int(1.00001*(table_max-table_min)/table_step)+1];
Float n_alpha_plus_beta_table[int(1.00001*(table_max-table_min)/table_step)+1];

static bool tables_initialized = false;
                        // True if we have filled out the tables.

///////////////////////////////////////////////////////////////////////////////
//
// Members of the Hodgkin_Huxley class:
//
Float Hodgkin_Huxley::g_Na_bar = 0.12*1.e-2; // Peak Na conductance (uS/um^2).
Float Hodgkin_Huxley::g_K_bar = 0.036*1.e-2; // Peak K conductance (uS/um^2).
Float Hodgkin_Huxley::g_leak = 0.0003*1.e-2; // Leak conductance (uS/um^2).

Float Hodgkin_Huxley::E_Na = 50; // Na reversal potential (mV).
Float Hodgkin_Huxley::E_K = -77.5; // K reversal potential (mV).
Float Hodgkin_Huxley::E_leak = -54.3; // Leak reversal potential (mV).

//
// Code to advance to the next timestep.  The algorithm is implicit Euler,
// like this:
//
//   dm
//   -- = (1-m) alpha(V) - m beta(V) = alpha - (alpha+beta) m
//   dt
//
//   m(t+dt)-m(t)
//   ------------ = alpha - (alpha(V)+beta(V)) m(t+dt)
//     dt
//
//             alpha + m(t)/dt       alpha dt + m(t)
//   m(t+dt) = --------------- = ------------------
//           1/dt+alpha+beta   1 + dt(alpha+beta)
//
// and similarly for the other variables.  Note that we use V as given, which
// means that this is not a fully implicit algorithm.  Note also that we don't
// care if the timestep is changed, since we don't have to precompute any
// constants that depend on the timestep.
//
// Arguments:
// 1) The time step.
// 2) The current value of the voltage.
//
void
Hodgkin_Huxley::advance(Float dt, Float V)
{
//
// Prepare to get the alpha and beta values by table lookup.  All the tables
// have the same size, so we compute the offset into the table just once.
//
  Float norm_v = (V - table_min)/table_step; // Compute the index into the
                        // table.
  unsigned idx = unsigned(norm_v);// Get an integer index.
#ifndef NDEBUG
  if (idx >= sizeof (m_alpha_table)/sizeof (Float) - 1)
  {
    cerr << "Hodgkin_Huxley::advance: voltage (" << V << ") too high\n";
    idx = sizeof (m_alpha_table)/sizeof (Float) - 2;
    norm_v = (Float)idx;
  }
#endif
  norm_v -= idx;        // Compute an offset from that index.
  Float one_minus_norm_v = 1-norm_v; // Precompute this, since we use it
                        // three times.

//
// Update m:
//
  Float alpha = m_alpha_table[idx]*one_minus_norm_v + m_alpha_table[idx+1]*norm_v;
  Float alpha_plus_beta =
    m_alpha_plus_beta_table[idx]*one_minus_norm_v + m_alpha_plus_beta_table[idx+1]*norm_v;
  m = (alpha*dt + m)/(1 + dt*alpha_plus_beta);

//
// Update h:
//
  alpha = h_alpha_table[idx]*one_minus_norm_v + h_alpha_table[idx+1]*norm_v;
  alpha_plus_beta =
    h_alpha_plus_beta_table[idx]*one_minus_norm_v + h_alpha_plus_beta_table[idx+1]*norm_v;
  h = (alpha*dt + h)/(1 + dt*alpha_plus_beta);

//
// Update n:
//
  alpha = n_alpha_table[idx]*one_minus_norm_v + n_alpha_table[idx+1]*norm_v;
  alpha_plus_beta =
    n_alpha_plus_beta_table[idx]*one_minus_norm_v + n_alpha_plus_beta_table[idx+1]*norm_v;
  n = (alpha*dt + n)/(1 + dt*alpha_plus_beta);
}

//
// Initialize the variables to their steady state value for some given voltage.
// Also fills in the tables for alpha and beta if they are not already set up.
// Arguments:
// 1) The voltage to initialize the variables to.
//
void
Hodgkin_Huxley::initialize(Float ini_v)
{
  if (!tables_initialized)    // Have not initialized tables of alpha and
  {                     // beta?
    for (unsigned idx = 0; idx < sizeof (m_alpha_table)/sizeof (Float); ++idx)
    {
      m_alpha_table[idx] = m_alpha(table_min + table_step *idx);
      m_alpha_plus_beta_table[idx] = m_alpha_plus_beta(table_min + table_step *idx);

      h_alpha_table[idx] = h_alpha(table_min + table_step *idx);
      h_alpha_plus_beta_table[idx] = h_alpha_plus_beta(table_min + table_step *idx);

      n_alpha_table[idx] = n_alpha(table_min + table_step *idx);
      n_alpha_plus_beta_table[idx] = n_alpha_plus_beta(table_min + table_step *idx);
    }

    tables_initialized = true;      // Don't reinitialize them.
  }

//
// Now initialize the state variables to their steady state value at this
// voltage.  The steady state value is
//             alpha
//   m          =  ----------
//    infty    alpha+beta
//
//
// Get the alpha and beta values by table lookup.  All the tables
// have the same size, so we compute the offset into the table just once.
// Speed probably isn't important here, but I already wrote the code for
// advance() and it's easy to duplicate it.
//
  Float norm_v = (ini_v - table_min)/table_step; // Compute the index into the
                        // table.
  unsigned idx = unsigned(norm_v);// Get an integer index.
#ifndef NDEBUG
  if (idx >= sizeof (m_alpha_table)/sizeof (Float) - 1)
  {
    cerr << "Hodgkin_Huxley::initialize: voltage (" << ini_v << ") too high\n";
    idx = sizeof (m_alpha_table)/sizeof (Float) - 2;
    norm_v = (Float)idx;
  }
#endif
  norm_v -= idx;        // Compute an offset from that index.
  Float one_minus_norm_v = 1-norm_v; // Precompute this, since we use it
                        // three times.

  Float alpha = m_alpha_table[idx]*one_minus_norm_v + m_alpha_table[idx+1]*norm_v;
  Float alpha_plus_beta =
    m_alpha_plus_beta_table[idx]*one_minus_norm_v + m_alpha_plus_beta_table[idx+1]*norm_v;
  m = alpha/alpha_plus_beta;

  alpha = h_alpha_table[idx]*one_minus_norm_v + h_alpha_table[idx+1]*norm_v;
  alpha_plus_beta =
    h_alpha_plus_beta_table[idx]*one_minus_norm_v + h_alpha_plus_beta_table[idx+1]*norm_v;
  h = alpha/alpha_plus_beta;

  alpha = n_alpha_table[idx]*one_minus_norm_v + n_alpha_table[idx+1]*norm_v;
  alpha_plus_beta =
    n_alpha_plus_beta_table[idx]*one_minus_norm_v + n_alpha_plus_beta_table[idx+1]*norm_v;
  n = alpha/alpha_plus_beta;
}

Generated by  Doxygen 1.6.0   Back to index