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

single_axon.cxx

//
// Code to control a simulation of two axons with the surrouding
// extracellular space.
//
// Written by: Gary Holt  14 May 1997.
//

#include <iostream.h>
#include "single_axon.h"

//
// Scan a sorted array looking for the element closest in value to the given
// value.  Normally, the value we're looking for is close to the last value
// we returned, so we do a linear search in the direction specified.
//
// Arguments:
// 1) The array to search.
// 2) The value to look for.
// 3) The previous value.
//
// Returns the largest index whose corresponding element is less than or equal
// to the given value.
//
static unsigned
scan_array_for_value(const array<Float> &arr, Float val, unsigned idx)
{
  if (val > arr[idx])         // Do we have to move forward?
  {
    while (idx+1 < arr.length() && val > arr[idx+1]) // Advance forward until
      ++idx;                  // we find the appropriate element.
  }
  else                        // Move backward:
  {
    while (idx > 0 && val < arr[idx]) // Advance backward until
      --idx;                  // we find the appropriate element.
  }

  return idx;
}

///////////////////////////////////////////////////////////////////////////////
//
// Members of the SingleAxon class:
//

Float SingleAxon::E_syn = 0;  // Reversal potential for the synapse.

//
// Constructor for the ExAxon class.  Arguments:
//
// 1) The z grid (the coordinates of the centers of each compartment).
// 2) The number of compartments.
// 3) The diameter of the axon (um).
// 4) The capacitance (uF/cm^2).
// 5) The intracellular resistivity (Ohm-cm).
//
SingleAxon::SingleAxon(const Float *arg_z_grid, unsigned n_compartments,
                   Float axon_diameter, Float capacitance,
                   Float arg_resistivity) :
  g_syn(0),             // No synaptic input.
  last_k(0),                  // Set this to something valid (doesn't matter
                        // particularly what).
  last_dt(0.05),        // Set this to something reasonable so we
                        // don't get a divide by 0 error.
  axon_radius(0.5*axon_diameter),// Store this value.
  membrane_capacitance(1.e3*1.e-8*capacitance), // Convert from uF/cm^2 to nF/um^2.
  resistivity_int(arg_resistivity * 1.e-6 * 1.e4), // Ohm-cm to MOhm-um.
  compartments(n_compartments),     // Allocate space for the comparments.
  subdiag(n_compartments+1),  // Allocate space for subdiagonal terms.
  z_grid(n_compartments)
{
  for (unsigned k = 0; k < n_compartments; ++k)
  {
    z_grid[k] = arg_z_grid[k];      // Make a copy of the z grid.
    AxonCompartment *c_ptr = &compartments[k];

    if (k == 0)
    {
      if (n_compartments == 1)      // Only one compartment?
      c_ptr->compartment_length = 2 * arg_z_grid[0];
      else
      c_ptr->compartment_length = arg_z_grid[1]-arg_z_grid[0];
    }
    else if (k == n_compartments-1)
      c_ptr->compartment_length = arg_z_grid[n_compartments-1] - arg_z_grid[n_compartments-2];
    else
      c_ptr->compartment_length = 0.5*(arg_z_grid[k+1]-arg_z_grid[k-1]);

    c_ptr->capacitance = membrane_capacitance * M_PI * axon_diameter * c_ptr->compartment_length;
                        // Compute the capacitance of this segment.
    c_ptr->Vm = c_ptr->last_Vm = -65;// Initialize to -65.
  }

  setup_subdiag();            // Setup the subdiagonal terms.
}

//
// Update all the conductances in one of the membranes, and adjust the
// coefficients in the matrix.  Arguments:
// 1) The timestep.
// 2) The extracellular voltages.
//
void
SingleAxon::advance(Float timestep, const Float *Ve)
{
//
// Compute the intracellular voltage all along this cable, holding the
// extracellular voltage and the Hodgkin-Huxley state variables constant.
// The equation to solve is:
//        2                    2
//           d V     dV             d V
//    1     m        m                 1     e
//    -- --- = C --- + g  (E  - V) - --   ----
//    r    2       dt     hh  hh   m   r    2
//     a dz                   a     dz
//
// which when discretized by the usual finite difference method gives
//
//       Vm       [                    c  L           ]       Vm
//    1    k-1    [ 1  ( 1    1    )    a  k          ]    1    k+1
//    -- -----  - [ -- ( -- + ---- ) + ----- + g   L  ]  + -- ----- =
//    r  h        [ r  ( h    h    )     dt     hh  k ]    r   h
//     a  k       [  a (  k    k+1 )                  ]     a   k+1
//
//
//     c  L                             Ve      Ve                      Ve
//      a  k                         1    k-1     k ( 1    1    )   1    k+1
//   - ----- Vm(t-dt) - g  L  E    - -- ----- + --- ( -- + ---- ) - -- -----
//       dt              hh k  hh    r    h     r   ( h    h    )   r   h
//                                    a    k     a  (  k    k+1 )    a   k+1
//
// where h_k = z_k - z_{k-1}, c_a is the capacitance per unit length,
// r_a is the axial resistance per unit length, and L_k is the length of the
// k'th segment (L_k = 0.5*(h_k + h_{k+1})).
//
// The boundary condition is dVm/dz = -dVe/dz (i.e., dVi/dt = 0).
// Ve = 0 at the boundaries.  We implement the Neumann conditions by
// inventing a fictitious point an equal distance from the boundary on
// the other side (i.e., h_0 = h_1, k_{K+1} = h_K):
//
//    (Vm_1 - Vm_0)/h_1 = -(Ve_1 - 0)/h_1
//    (Vm_K - Vm_{K-1})/h_K = -(0 - Ve_K)/h_K
//
// Substituting this into the fully discretized equation causes the
// coefficients of Vm_1 and Vm_K to be adjusted.
// (Actually, for convenience we've used a Neumann condition on Ve.  I don't
// think this really matters.)
//
// This is a tridiagonal system, so we can solve it in O(N) time.  
// We store the right hand side in the Vm values themselves, since they
// aren't needed during the computation.  The subdiagonal/superdiagonal
// have already been precomputed by set_grid.  Note that since r_a is
// constant, the subdiagonal element of each row is equal to the
// superdiagonal element of the previous row.
//
  Float circum = 2 * M_PI * axon_radius;// Circumference of the axon.
  Float neg_diag[compartments.length()];// Where we store the diagonal elements.

  for (unsigned k = 0; k < compartments.length(); ++k)
  {
    AxonCompartment *c_ptr = &compartments[k];
                        // Point to this compartment.
    c_ptr->advance(timestep, c_ptr->Vm);// Update the HH state variables.
    Float g_hh, g_hh_E_hh;
    c_ptr->net_g_and_E(g_hh, g_hh_E_hh, c_ptr->compartment_length * circum);
                        // Get the total conductance.
    Float C_Lk_over_dt = c_ptr->capacitance/timestep;
                        // Compute the capacitive contribution.

    neg_diag[k] = subdiag[k] + subdiag[k+1] + C_Lk_over_dt + g_hh;
                        // Store the diagonal.  subdiag[k+1] is the
                        // superdiagonal element of this row.
    c_ptr->Vm = -C_Lk_over_dt * c_ptr->last_Vm - g_hh_E_hh 
      + (subdiag[k] + subdiag[k+1]) * Ve[k];

    if (k < compartments.length()-1) // Does Ve[k+1] exist?
      c_ptr->Vm -= subdiag[k+1] * Ve[k+1];

    if (k == 0)         // Is this the first compartment in the axon?
    {
      neg_diag[k] += g_syn;   // Add the synaptic conductance.
      c_ptr->Vm -= g_syn * E_syn;
    }
    else
    {                   // k >= 1:
      c_ptr->Vm -= subdiag[k] * Ve[k-1]; // Ve[k-1] exists, so add it in.
//
// At this point we've computed the previous row.  So now eliminate the
// subdiagonal term from this row:
//
      Float tmp_factor = subdiag[k]/neg_diag[k-1];
                        // What to multiply the previous row
                        // by before adding this row's leading
                        // coefficient.
      neg_diag[k] -= tmp_factor * subdiag[k];
                        // Adjust our diagonal term.
      c_ptr->Vm += c_ptr[-1].Vm * tmp_factor;// Adjust RHS.
    }
  }

//
// We've made one pass through the matrix, eliminating the subdiagonal
// elements.  Go through now and backsubstitute to solve:
 //
  Float last_vm = (compartments[compartments.length()-1].Vm /=
               -neg_diag[compartments.length()-1]);
                        // Compute the final value.
  for (int k = compartments.length()-2; k >= 0; --k)
  {
    AxonCompartment *c_ptr = &compartments[k];
    c_ptr->Vm = last_vm = (last_vm * subdiag[k+1]-c_ptr->Vm) / neg_diag[k];
                        // Compute the intracellular voltage.
                        // subdiag[k+1] = superdiagonal[k].

  }

  last_dt = timestep;         // Remember the timestep so we can compute
                        // what the capacitative currents are.
}

//
// Copy the current value of the Hodgkin-Huxley state variables into temporary
// storage so we can back up again.
//
void
SingleAxon::save_hh_state()
{
  for (unsigned k = 0; k < compartments.length(); ++k)
    compartments[k].last_hh = (Hodgkin_Huxley)compartments[k];
                        // Just copy the whole HH structure.
}

//
// Restore the last saved Hodgkin-Huxley state variables.
//
void
SingleAxon::restore_hh_state()
{
  for (unsigned k = 0; k < compartments.length(); ++k)
    *(Hodgkin_Huxley *)&compartments[k] = compartments[k].last_hh;
                        // Just copy the whole HH structure.
}

//
// Save the transmembrane potentials so we can compute time derivatives
// on the next time step.
//
void
SingleAxon::save_voltages()
{
  for (unsigned k = 0; k < compartments.length(); ++k)
    compartments[k].last_Vm = compartments[k].Vm;
}

//
// Initialize for a particular transmembrane potential.  The extracellular
// potential is assumed to be 0.
//
void
SingleAxon::initialize(Float new_vm)
{
//
// Now set the intracellular potentials to the given value, and initialize
// the HH variables:
//
  for (unsigned k = 0; k < compartments.length(); ++k)
  {
    AxonCompartment *c_ptr = &compartments[k];
    c_ptr->last_Vm = c_ptr->Vm = new_vm; // Set the voltage.
    c_ptr->initialize(new_vm);      // Set up the HH state variables.
    c_ptr->last_hh.initialize(new_vm); // Same for the previous HH state.
  }
}

//
// Return the membrane current density through a particular patch of membrane.
// Arguments:
// 1) The x position.
//
Float
SingleAxon::Jm(Float xpos) const
{
  ((SingleAxon *)this)->last_k = scan_array_for_value(z_grid, xpos, last_k);

  const AxonCompartment *c_ptr = &compartments[last_k]; // Get the compartments.
  Float compartment_area = c_ptr->compartment_length * 2 * M_PI * axon_radius;

  Float net_g, net_g_E;

  c_ptr->net_g_and_E(net_g, net_g_E, 1); // Get the HH conductance.
  if (last_k == 0)            // Is this at the tip?
  {                     // Add in the synapse.
    Float g = g_syn / compartment_area;   // Calculate conductance per area.
    net_g += g;               // Update the currents.
    net_g_E += g * E_syn;
  }

  return net_g * c_ptr->Vm - net_g_E +
    c_ptr->capacitance * (c_ptr->Vm - c_ptr->last_Vm) / (last_dt * compartment_area);
}

//
// Get all the transmembrane currents from each compartment.  This is the
// total current through the membrane of the compartment, not the current
// per unit area.
//
void
SingleAxon::Im_compartment(Float *Im) const
{
  Float circum = 2 * M_PI * axon_radius;

  for (unsigned k = 0; k < compartments.length(); ++k)
  {
    const AxonCompartment *c_ptr = &compartments[k]; // Get the compartments.
    Float net_g, net_g_E;

    c_ptr->net_g_and_E(net_g, net_g_E, c_ptr->compartment_length * circum);
                        // Get the HH conductance.
    if (k == 0)               // Is this at the tip?
    {                   // Add in the synapse:
      Float g = g_syn; // Calculate conductance per area.
      net_g += g;             // Update the currents.
      net_g_E += g * E_syn;
    }

    Im[k] = net_g * c_ptr->Vm - net_g_E +
      c_ptr->capacitance * (c_ptr->Vm - c_ptr->last_Vm) / last_dt;
  }
}

//
// Return pointers to the axon compartments:
//
void
SingleAxon::get_compartments(AxonCompartment **cvec) const
{
  for (unsigned k = 0; k < compartments.length(); ++k)
    cvec[k] = (AxonCompartment *)&compartments[k];
}

//
// Get the transmembrane voltage at a given position.  Arguments:
// 1) The x position (um).
//
Float
SingleAxon::vm(Float xpos) const
{
  ((SingleAxon *)this)->last_k = scan_array_for_value(z_grid, xpos, last_k);

  return compartments[last_k].Vm;
}

//
// Get the whole vector of voltages.  Arguments:
// 1) Where to return the vector.  Should have the correct dimensions.
//
void
SingleAxon::get_vm(Float *Vm_out) const
{
  for (unsigned k = 0; k < compartments.length(); ++k)
    Vm_out[k] = compartments[k].Vm;
}

//
// Set up the subdiagonal elements.  Arguments:
// 1) The intracellular resistivity (MOhm/um).  Note that this is resistance
//    per unit length, not bulk resistivity.
//
void
SingleAxon::setup_subdiag()
{
//
// The subdiagonal terms in our tridiagonal system are of the form
//
//   -1        1      while the superdiagonals    -1     1
//   --     ---------   have the form                 -- ---------
//   r      z  - z                              r  z    - z
//    a      k    k-1                      a  k+1        k
//
// (see the comments in ActiveCable::advance for more details).  This means
// that the value of the superdiagonal is just the value of the subdiagonal
// on the next row.
//
// Because of our Neumann boundary conditions, the first row has a
// subdiagonal of 0 and the last row has a superdiagonal of 0.
//
  Float ra = resistivity_int / (M_PI * axon_radius * axon_radius);
                        // Calculate resistance per unit length.

  if (subdiag.length() != z_grid.length()+1)
    subdiag.reallocate(z_grid.length()+1);
                        // Make sure it's the right size.

  subdiag[0] = 0;       // Fill in the first term.
  for (unsigned k = 1; k < z_grid.length(); ++k)
    subdiag[k] = 1./(ra * (z_grid[k] - z_grid[k-1]));
                        // Fill in the intermediate terms.
  subdiag[z_grid.length()] = 0;     // Fill in the last superdiagonal.
}

template class array<Float>;  // Instantiate our templates.
template class array<AxonCompartment>;

Generated by  Doxygen 1.6.0   Back to index