kschan.h

Go to the documentation of this file.

#ifndef kschan_h
#define kschan_h

#include <math.h>
#include <OS/string.h>
#include "nrnoc2iv.h"
#include "ivocvect.h"

extern "C" {
#include "spmatrix.h"
}

extern "C" {
   //extern double dt;
   extern double celsius;
}

class KSState;
class KSChan;
class KSSingle;
struct NrnThread;

class KSChanFunction {
public:
   KSChanFunction();
   virtual ~KSChanFunction();
   virtual int type() { return 0; }
   virtual double f(double v) { return 1.; }
   void f(int cnt, double* v, double* val) {
      int i; for (i=0; i < cnt; ++i) { val[i] = f(v[i]); }
   }
   static KSChanFunction* new_function(int type, Vect*, double , double);
   Vect* gp_; // gate parameters
   double c(int i) { return gp_->elem(i); }
   static double Exp(double x) { if (x > 700.) {return exp(700.);
            } else if (x < -700.) { return exp(-700.);
            } else { return exp(x); } }
};

class KSChanConst : public KSChanFunction {
public:
   virtual int type() { return 1; }
   virtual double f(double v) { return c(0); }
};

class KSChanExp : public KSChanFunction {
public:
   virtual int type() { return 2; }
   virtual double f(double v) { return c(0)*Exp(c(1) * (v - c(2))); }
};

class KSChanLinoid : public KSChanFunction {
public:
   virtual int type() { return 3; }
   virtual double f(double v) {
      double x = c(1)*(v - c(2));
      if (fabs(x) > 1e-6) {
         return c(0)*x/(1 -  Exp(-x));
      }else{
         return c(0)*(1 + x/2.);
      }
   }
};

class KSChanSigmoid : public KSChanFunction {
public:
   virtual int type() { return 4; }
   virtual double f(double v) {
      return c(0) / (1. + Exp(c(1) * (v - c(2))));
   }
};


// e/(kT) e/k=11.604589 from hoc's FARADAY and R values
#define ebykt (11.604589/(273.15 + celsius))

// from MODELING NEURONAL BIOPHYSICS Lyle J Graham
//A Chapter in the Handbook of Brain Theory and Neural Networks, Volume 2
//a' = K*exp(z*gam*(v-vhalf)*F/RT)
//b' = K*exp(-z*(1-gam)*(v-vhalf)*F/RT)
//but tau = 1/(a' + b') + tau0
//and inf = a'/(a' + b')   
// so there is no fast way to get either a,b or inf,tau individually

class KSChanBGinf : public KSChanFunction {
public:
   virtual int type() { return 5; }
   virtual double f(double v) {
      double x = ebykt*c(2)*(v - c(1));
      double ap = c(0) * Exp(c(3)*x);
      double bp = c(0) * Exp((c(3) - 1.)*x);
      tau = 1/(ap + bp);
      double inf = ap*tau; tau += c(4);
      return inf;
   }
   double tau; // may avoid duplicate call to KSChanBGtau
};

class KSChanBGtau : public KSChanFunction {
public:
   virtual int type() { return 6; }
   virtual double f(double v) {
      double x = ebykt*c(2)*(v - c(1));
      double ap = c(0) * Exp(c(3)*x);
      double bp = c(0) * Exp((c(3) - 1.)*x);
      double tau = 1/(ap + bp);
      inf = ap*tau; tau += c(4);
      return tau;
   }
   double inf; // may avoid duplicate call to KSChanBGinf
};

class KSChanTable : public KSChanFunction {
public:
   KSChanTable(Vect*, double vmin, double vmax);
   virtual int type() { return 7; }
   virtual double f(double v);
   double vmin_, vmax_;
private:
   double dvinv_;
};

class KSTransition {
public:
   KSTransition();
   virtual ~KSTransition();
   // vmin, vmax only for type KSChanTable
   void setf(int direction, int type, Vect* vec, double vmin, double vmax);
   // only voltage gated for now
   double alpha(double v) { return type_ == 0 ? f0->f(v) : f0->f(v)/f1->f(v); }
   double beta(double v) { return type_ == 0 ? f1->f(v) : (1. - f0->f(v))/f1->f(v); }
   double inf(double v) { return type_ == 1 ? f0->f(v) : f0->f(v)/(f0->f(v) + f1->f(v)); }
   double tau(double v) { return type_ == 1 ? f1->f(v) : 1./(f0->f(v) + f1->f(v)); }
   void ab(double v, double& a, double& b);
   void ab(Vect* v, Vect* a, Vect* b);
   void inftau(double v, double& inf, double& tau);
   void inftau(Vect* v, Vect* inf, Vect* tau);
   // hh tables
   // easily out of date!!
   // in anything about f0 or f1 changes then must call hh_table_make;
   void hh_table_make(double dt, int size=200, double vmin=-100., double vmax=50.);
   boolean usehhtable() { return (size1_ > 0); }
   void inftau_hh_table(int i, double& inf, double& tau) {
      inf = inftab_[i];
      tau = tautab_[i];
   }
   void inftau_hh_table(int i, double x, double& inf, double& tau) {
      inf = inftab_[i] + (inftab_[i+1] - inftab_[i])*x;
      tau = tautab_[i] + (tautab_[i+1] - tautab_[i])*x;
   }
   // the agent style
   virtual double alpha(Datum*);
   virtual double beta();
   void lig2pd(int pdoff);
public:
   Object* obj_;
   int index_; // into trans_ array
   int src_;
   int target_;
   KSChan* ks_;
   KSChanFunction* f0;
   KSChanFunction* f1;
   int type_; // 0-ab,voltage gated; 1-inftau voltage gated
         // 2-ligand outside;  3-ligand inside
   int ligand_index_;
   int pd_index_;
   int stoichiom_;
private:
   // for hh tables.
   double* inftab_;
   double* tautab_;
   int size1_;
};

class KSGateComplex {
public:
   KSGateComplex();
   virtual ~KSGateComplex();
   double conductance(double* state, KSState* st);
public:
   Object* obj_;
   KSChan* ks_;
   int index_; // into gc_ array
   int sindex_; // starting state index for this complex
   int nstate_; // number of states
   int power_; // eg. n^4, or m^3   
};

class KSIv {
public:
   // this one for ionic ohmic and nernst.
   virtual double cur(double g, double* p, Datum* pd, double v);
   virtual double jacob(double* p, Datum* pd, double v);
};
class KSIvghk : public KSIv {
public:
   // this one for ionic Goldman-Hodgkin-Katz
   virtual double cur(double g, double* p, Datum* pd, double v);
   virtual double jacob(double* p, Datum* pd, double v);
   double z;
};
class KSIvNonSpec : public KSIv{
   // this one for non-specific ohmic. There will be a PARAMETER e_suffix at p[1]
   virtual double cur(double g, double* p, Datum* pd, double v);
   virtual double jacob(double* p, Datum* pd, double v);
};

class KSPPIv : public KSIv {
public:
   // this one for POINT_PROCESS ionic ohmic and nernst.
   virtual double cur(double g, double* p, Datum* pd, double v);
   virtual double jacob(double* p, Datum* pd, double v);
   int ppoff_;
};
class KSPPIvghk : public KSPPIv {
public:
   // this one for POINT_PROCESS ionic Goldman-Hodgkin-Katz
   virtual double cur(double g, double* p, Datum* pd, double v);
   virtual double jacob(double* p, Datum* pd, double v);
   double z;
};
class KSPPIvNonSpec : public KSPPIv{
   // this one for POINT_PROCESS non-specific ohmic. There will be a PARAMETER e_suffix at p[1]
   virtual double cur(double g, double* p, Datum* pd, double v);
   virtual double jacob(double* p, Datum* pd, double v);
};

class KSState {
public:
   KSState();
   virtual ~KSState();
   const char* string() { return name_.string(); }
   double f_; // normalized conductance
   CopyString name_; 
   int index_; // into state_ array
   KSChan* ks_;
   Object* obj_;
};

class KSChan {
public:
   KSChan(Object*, boolean is_point = false);
   virtual ~KSChan();
   virtual void alloc(Prop*);
   virtual void init(int, Node**, double**, Datum**, NrnThread*);
   virtual void cur(int, Node**, double**, Datum**);
   virtual void jacob(int, Node**, double**, Datum**);
   virtual void state(int, Node**, double**, Datum**, NrnThread*);
#if CACHEVEC != 0
   virtual void cur(int, int *, double**, Datum**, NrnThread*);
   virtual void jacob(int, int *, double**, Datum**, NrnThread*);
   virtual void state(int, int *, Node**, double**, Datum**, NrnThread*);
#endif /* CACHEVEC */
   void add_channel(char**);
   //for cvode
   virtual int count();
   virtual void map(int, double**, double**, double*, Datum*, double*);
   virtual void spec(int, Node**, double**, Datum**);
   virtual void matsol(int, Node**, double**, Datum**, NrnThread*);
   virtual void cv_sc_update(int, Node**, double**, Datum**, NrnThread*);
   double conductance(double gmax, double* state);
public:
   // hoc accessibilty
   int state(const char* name);
   const char* state(int index);
   int trans_index(const char* src, const char* target); // index of the transition
   int trans_index(int src, int target); // index of the transition
   int gate_index(int state_index); // index of the gate
   boolean is_point() { return is_point_; }
   boolean is_single() { return is_single_; }
   void set_single(boolean, boolean update = true);
   void destroy_pnt(Point_process*);
   int nsingle(Point_process*);
   void nsingle(Point_process*, int);
   double alpha(double v, int, int) {return 0.;}
   double beta(double v, int, int) {return 0.;}
   void setstructure(Vect*);
   void setname(const char*);
   void setsname(int, const char*);
   void setion(const char*);
   void setligand(int i, const char*);
   void settype(KSTransition*, int type, const char*);
   void setivrelation(int);
   // hoc incremental management
   KSState* add_hhstate(const char*);
   KSState* add_ksstate(int igate, const char*);
   void remove_state(int);
   // these are only for kinetic scheme transitions since an hh
   // always has one and only one transition.
   KSTransition* add_transition(int src, int target, const char* ligand);
   void remove_transition(int);
   void setcond();
   void power(KSGateComplex*, int);
   void usetable(boolean, int size, double vmin, double vmax);
   void usetable(boolean);
   int usetable(double* vmin, double* vmax);// get info
   boolean usetable() { return usetable_; }
   void check_table_thread(NrnThread*);
private:
   void free1();
   void build();
   void setupmat();
   void fillmat(double v, Datum* pd);
   void mat_dt(double dt, double* p);
   void solvemat(double*);
   void mulmat(double*, double*);
   void ion_consist();
   void ligand_consist(int, int, Prop*, Node*);
   Prop* needion(Symbol*, Node*, Prop*);
   void state_consist(int shift = 0);
   void sname_install();
   Symbol* looksym(const char*, Symbol* tmplt = nil);
   Symbol* installsym(const char*, int, Symbol* tmplt = nil);
   void freesym(Symbol*, Symbol* tmplt = nil);
   Symbol** newppsym(int);
   void delete_schan_node_data();
   void alloc_schan_node_data();
   void update_prop(); // can add and remove Nsingle and SingleNodeData

   KSState* state_insert(int i, const char* name, double frac);
   void state_remove(int i);
   KSGateComplex* gate_insert(int ig, int is, int power);
   void gate_remove(int i);
   KSTransition* trans_insert(int i, int src, int target);
   void trans_remove(int i);
   void check_struct();
   int state_size_;
   int gate_size_;
   int trans_size_;
   boolean is_point_;
   boolean is_single_;
   // for point process
   int pointtype_;
   int mechtype_;
public:
   CopyString name_; // name of channel
   CopyString ion_; // name of ion , "" means non-specific
   double gmax_deflt_;
   double erev_deflt_;
   int cond_model_;
   KSIv* iv_relation_;
   int ngate_; // number of gating complexes
   int ntrans_; // total number of transitions
   int ivkstrans_; // index of beginning of voltage sensitive ks transitions
   int iligtrans_; // index of beginning of ligand sensitive ks transitions
         // 0 to ivkstrans_ - 1 are the hh transitions
   int nhhstate_; // total number of hh states, does not include ks states
   int nksstate_; // total number of kinetic scheme states.
   int nstate_;   // total number of all states, nhhstate_ to nstate_ - 1
         // are kinetic scheme states. The nhhstate_ are first
         // and the nhhstate_ = ivkstrans_
   KSState* state_; // the state names
   KSGateComplex* gc_;
   KSTransition* trans_; // array of transitions
   Symbol* ion_sym_; // if nil then non-specific and e_suffix is a parameter
   int  nligand_;
   Symbol** ligands_;
   Object* obj_;
   KSSingle* single_;
private:
   int cvode_ieq_;
   Symbol* mechsym_; // the top level symbol (insert sym or new sym)
   Symbol* rlsym_; // symbol with the range list (= mechsym_ when  density)
   char* mat_;
   double** elms_;
   double** diag_;
   int dsize_; // size of prop->dparam
   int psize_; // size of prop->param
   int soffset_; //STATE begins here in the p array.
   int gmaxoffset_; // gmax is here in the p array, normally 0 but if
      // there is an Nsingle then it is 1.
   int ppoff_; //2 or 3 for point process since area and Point_process* are
         // first two elements and there may be a KSSingleNodeData*
   // for hh rate tables
   double vmin_, vmax_, dvinv_, dtsav_;
   int hh_tab_size_;
   boolean usetable_;
};

#endif