Hypocretin and Locus Coeruleus model neurons (Carter et al 2012)

 Download zip file 
Help downloading and running models
Accession:145162
Conductance based model of the hypocretin neurons (HCRT) and another one of the Locus Coeruleus one (LC). The HCRT drive the LCs via the HCRT receptor on the LCs. The LCs lead to the awakening of the mice if the number of spikes raises over 10 spikes in 10 seconds window.
Reference:
1 . Carter ME, Brill J, Bonnavion P, Huguenard JR, Huerta R, de Lecea L (2012) Mechanism for Hypocretin-mediated sleep-to-wake transitions PNAS 109(39):E2635-44 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Axon; Channel/Receptor;
Brain Region(s)/Organism:
Cell Type(s): Locus Coeruleus neuron;
Channel(s): I Na,t; I T low threshold; I A; I K; I K,Ca; I Calcium;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s):
Implementer(s): Huerta, Ramon [rhuerta at ucsd.edu];
Search NeuronDB for information about:  I Na,t; I T low threshold; I A; I K; I K,Ca; I Calcium;
/*
 This program solves the problem with the BDF method,
 * Newton iteration with the CVDENSE dense linear solver, 
 * -----------------------------------------------------------------
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Header files with a description of contents used here */

#include <cvode/cvode.h>             /* prototypes for CVODE fcts. and consts. */
#include <cvode/cvode_dense.h>       /* prototype for CVDense */
#include <nvector/nvector_serial.h>  /* serial N_Vector types, functions, and macros */

/* User-defined vector and matrix accessor macros: Ith, IJth */

/* These macros are defined in order to write code which exactly matches
   the mathematical problem description given above.

   Ith(v,i) references the ith component of the vector v, where i is in
   the range [1..NEQ] and NEQ is defined below. The Ith macro is defined
   using the N_VIth macro in nvector.h. N_VIth numbers the components of
   a vector starting from 0.*/

#define Ith(v,i)    NV_Ith_S(v,i-1)       /* Ith numbers components 1..NEQ */
//#define IJth(A,i,j) DENSE_ELEM(A,i-1,j-1) /* IJth numbers rows,cols 1..NEQ */


/* Problem Constants */

#define NEQ   7                /* number of equations  */
#define ATOL  RCONST(1.0e-14)   /* vector absolute tolerance components */
#define T0    RCONST(0.0)      /* initial time           */
#define T1    RCONST(0.4)      /* first output time      */
#define TMULT RCONST(0.01)     /* output time factor     */


/* Functions Called by the Solver */

static int f(double t, N_Vector y, N_Vector ydot, void *user_data);

/* Private functions to output results */

static void PrintOutput(double t, double y1, double y2);

/* Private function to print final statistics */

static void PrintFinalStats(void *cvode_mem);

/* Private function to check function return values */

static int check_flag(void *flagvalue, char *funcname, int opt);

void set_initial_conditions(N_Vector );

/*
 *-------------------------------
 * Main Program
 *-------------------------------
 */
#define PARAMETERS 4
#define TIMEMAX   1
#define TOLERANCE 2
#define IDC       3

double Idc;

int main(int argc, char **argv)
{
  double t, tout;
  N_Vector y=NULL;
  void *cvode_mem=NULL;
  int flag;
  double tmax;

  if (argc!=PARAMETERS) {
    puts("Max time in seconds");
    puts("Relative tolerance");
    puts("Current");
    return -1;
  }

  Idc=atof(argv[IDC]);
  tmax=atof(argv[TIMEMAX])*1000;

  /* Create serial vector of length NEQ for I.C. */
  y = N_VNew_Serial(NEQ);
  if (check_flag((void *)y, "N_VNew_Serial", 0)) return(1);

  /* Initialize y */
  set_initial_conditions(y);

  /* Call CVodeCreate to create the solver memory and specify the 
   * Backward Differentiation Formula and the use of a Newton iteration */
  cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
  if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
  
  /* Call CVodeInit to initialize the integrator memory and specify the
   * user's right hand side function in y'=f(t,y), the inital time T0, and
   * the initial dependent variable vector y. */
  flag = CVodeInit(cvode_mem, f, T0, y);
  if (check_flag(&flag, "CVodeInit", 1)) return(1);

  flag = CVodeSStolerances(cvode_mem, atof(argv[TOLERANCE]), ATOL);

  /* Call CVDense to specify the CVDENSE dense linear solver */
  flag = CVDense(cvode_mem, NEQ);
  if (check_flag(&flag, "CVDense", 1)) return(1);
  /* Set the Jacobian routine to internal estimation */
  flag=CVDlsSetDenseJacFn(cvode_mem, NULL);
  
  // flag = CVDlsSetDenseJacFn(cvode_mem, Jac);
  // if (check_flag(&flag, "CVDlsSetDenseJacFn", 1)) return(1);

  /* In loop, call CVode, print results, and test for error.
     Break out of loop when NOUT preset output times have been reached.  */

  tout = T1;
  while(1) {
    flag = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
    PrintOutput(t, Ith(y,1), Ith(y,7));

    if (check_flag(&flag, "CVode", 1)) break;
    if (flag == CV_SUCCESS) {
      tout += TMULT;
      //printf("%lg %lg\n",tout,TMULT);
    }

    if (tout>tmax) break;
  }

  /* Free y vector */
  N_VDestroy_Serial(y);

  /* Free integrator memory */
  CVodeFree(&cvode_mem);

  return(0);
}


/*
 * f routine. Compute function f(t,y). 
 */
#define CA  10.0 /* pico F */
#define CS  10.0 /* pico F */
#define gl  1.6 /* nano S */
#define gAS 65.0  /* nano S */ 
#define gNa 260.0   /* nano S */
#define gKd 80.0   /* nano S */
#define gCa 0.88   /* nano S */
#define MU  1.5    /* Calcium dynamics dissipation*/

#define ABS( X ) (((X)>0.0)?(X):-(X))
#define Cagam(X,Y,Z) (1.0/(1.0+exp(((X)-(Y))/(Z)))) 

static int f(double t, N_Vector y, N_Vector ydot, void *user_data)
{
  double Va=Ith(y,1);
  double Vs=Ith(y,2);
  double Vt=-50.14905;
  double v;
  double a;
  double INa, IKd, ICa;
  
  /* INa current [Traub & Miles, 1991] */
  INa=gNa*Ith(y,3)*Ith(y,3)*Ith(y,4)*(Va-50);

  /* INa activation Ith(y,3)*/
  v=18.0+Vt-Va;
  if (ABS(v)<0.0001) {
    a=1.28-0.16*v;
  } else 
    a=0.32*v/(exp(v/4.0)-1);    
  
  Ith(ydot,3)=a*(1-Ith(y,3));

  v=-40.0-Vt+Va;
  if (ABS(v)<0.0001) {
    a=1.40-0.14*v;
  } else 
    a=0.28*v/(exp(v/5.0)-1.0);
  Ith(ydot,3)-=a*(Ith(y,3));

  /* INa inactivation Ith(y,4) */
  v=17+Vt-Va;
  Ith(ydot,4)=0.128*exp(v/18.0)*(1-Ith(y,4));
  v=40.0+Vt-Va;
  Ith(ydot,4)-=4.0*Ith(y,4)/(1.0+exp(v/5.0));

  /* IKd delayed rectifier current */
  IKd=gKd*Ith(y,5)*(Va+60);
  
  /* IKd activation Ith(y,5) */
  v=35.0+Vt-Va;
  if (ABS(v)<0.0001) {
    a=0.08-0.008*v;
  } else 
    a=0.016*v/(exp(v/5.0)-1.0);
  v=20.0+Vt-Va;
  Ith(ydot,5)=a*(1-Ith(y,5))-0.25*exp(v/40)*Ith(y,5);
  
  /* ICa current  */
  if (ABS(Vs)<0.001) {
    a=0.5*(- 24.42002442+Vs);
  } else {
    a=Vs/(1-exp(2*Vs/ 24.42002442));
  }
  ICa=gCa*Ith(y,6)*Ith(y,6)*Ith(y,6)*a;
  /* Calcium activation Ith(y,6) */
  Ith(ydot,6)=0.1*(Cagam(-Vs,39.1,2)-Ith(y,6));
  //printf("%lg ydot %lf\n",ICa,Ith(ydot,6));exit(0);
  /* Axon membrane potential */

  /* [Ca] Ith(y,7)*/
  Ith(ydot,7)=0.001*(-0.35*ICa-MU*MU*Ith(y,7)+0.04*MU*MU);

  Ith(ydot,1) = (1.0/CA)*(-gl*(Va+45.0 /*resting potential*/)-gAS*(Va-Vs)
			  -INa-IKd);
  
  //printf("ydo1=%lf \n",Ith(ydot,1));
  //exit(0);
  /* Soma membrane pontential */
  Ith(ydot,2) = (1.0/CS)*(-gl*(Vs+45.0 /*resting potential*/)-gAS*(Vs-Va)
			  +ICa
			 +Idc);
  //printf("ydo1=%lg (%lf,%lf,%lf,%lf)\n",Ith(ydot,1),gl*(Va+45.0 /*resting potential*/),gAS*(Va-Vs),INa,IKd);
  //exit(0);
  return(0);      
}


/*
 *-------------------------------
 * Private helper functions
 *-------------------------------
 */

static void PrintOutput(double t, double y1, double y2)
{

  printf("%0.4e %14.6e  %14.6e\n", t, y1, y2*100);

  return;
}

void set_initial_conditions(N_Vector y) {
  int i;
  for(i=3;i<=NEQ;++i)   Ith(y,i)=0.1;
  
  Ith(y,1) = -45.0;
  Ith(y,2) = -46.0;
  Ith(y,7) = 0.04;

}

/* 
 * Get and print some final statistics
 */

static void PrintFinalStats(void *cvode_mem)
{
  long int nst, nfe, nsetups, nje, nfeLS, nni, ncfn, netf, nge;
  int flag;

  flag = CVodeGetNumSteps(cvode_mem, &nst);
  check_flag(&flag, "CVodeGetNumSteps", 1);
  flag = CVodeGetNumRhsEvals(cvode_mem, &nfe);
  check_flag(&flag, "CVodeGetNumRhsEvals", 1);
  flag = CVodeGetNumLinSolvSetups(cvode_mem, &nsetups);
  check_flag(&flag, "CVodeGetNumLinSolvSetups", 1);
  flag = CVodeGetNumErrTestFails(cvode_mem, &netf);
  check_flag(&flag, "CVodeGetNumErrTestFails", 1);
  flag = CVodeGetNumNonlinSolvIters(cvode_mem, &nni);
  check_flag(&flag, "CVodeGetNumNonlinSolvIters", 1);
  flag = CVodeGetNumNonlinSolvConvFails(cvode_mem, &ncfn);
  check_flag(&flag, "CVodeGetNumNonlinSolvConvFails", 1);

  flag = CVDlsGetNumJacEvals(cvode_mem, &nje);
  check_flag(&flag, "CVDlsGetNumJacEvals", 1);
  flag = CVDlsGetNumRhsEvals(cvode_mem, &nfeLS);
  check_flag(&flag, "CVDlsGetNumRhsEvals", 1);

  flag = CVodeGetNumGEvals(cvode_mem, &nge);
  check_flag(&flag, "CVodeGetNumGEvals", 1);

  printf("\nFinal Statistics:\n");
  printf("nst = %-6ld nfe  = %-6ld nsetups = %-6ld nfeLS = %-6ld nje = %ld\n",
	 nst, nfe, nsetups, nfeLS, nje);
  printf("nni = %-6ld ncfn = %-6ld netf = %-6ld nge = %ld\n \n",
	 nni, ncfn, netf, nge);
}

/*
 * Check function return value...
 *   opt == 0 means SUNDIALS function allocates memory so check if
 *            returned NULL pointer
 *   opt == 1 means SUNDIALS function returns a flag so check if
 *            flag >= 0
 *   opt == 2 means function allocates memory so check if returned
 *            NULL pointer 
 */

static int check_flag(void *flagvalue, char *funcname, int opt)
{
  int *errflag;

  /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
  if (opt == 0 && flagvalue == NULL) {
    fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n",
	    funcname);
    return(1); }

  /* Check if flag < 0 */
  else if (opt == 1) {
    errflag = (int *) flagvalue;
    if (*errflag < 0) {
      fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n",
	      funcname, *errflag);
      return(1); }}

  /* Check if function returned NULL pointer - no memory allocated */
  else if (opt == 2 && flagvalue == NULL) {
    fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n",
	    funcname);
    return(1); }

  return(0);
}

Loading data, please wait...