implemented Hodgkin-Huxley model with K,Na,leak conductances

Makefile

@@ -1,2 +1,22 @@
-proj1 : 
-	gcc main.c -O3 -o proj1 
+CC = gcc
+CFLAGS = -O3
+GSLFLAGS = -lgsl -lgslcblas -lm
+
+all : proj1
+
+proj1: neuron.o general.o ode.o main.o
+	$(CC) $(CFLAGS) $(GSLFLAGS) neuron.o general.o ode.o main.o -o proj1
+
+neuron.o: neuron.c
+	$(CC) $(CFLAGS) -c neuron.c
+
+general.o: general.c
+	$(CC) $(CFLAGS) -c general.c
+
+ode.o: ode.c
+	$(CC) $(CFLAGS) $(GSLFLAGS) -c ode.c
+
+main.o: main.c
+	$(CC) $(CFLAGS) $(GSLFLAGS) -c main.c
+clean:
+	rm *.o proj1

README

@@ -1,7 +1,7 @@
 ======================
 USAGE
 ======================
-
+	user$>make
 	user$>./proj1 [config_file] [output_filename]
 
 If config_file is not specified, template.txt is used. If output_filename is given, the state of the network at the end of all calculations is output to output_filename as a "defined" type network.

definitions.h

@@ -1,3 +1,6 @@
+#ifndef DEFINITIONS
+#define DEFINITIONS
 #define MAX_LINE_LEN 80
 #define DEBUG 0
 #define DEFAULT_CONFIG_FILE "template.txt"
+#endif

general.c

@@ -1,3 +1,6 @@
+#include "definitions.h"
+#include "includes.h"
+
 char * get_input_filename(int argc, char** argv)
 {
   char * input_filename;

general.h

@@ -1 +1,6 @@
+#ifndef GENERAL_H
+#define GENERAL_H
+
 char * get_input_filename(int argc, char** argv);
+
+#endif

includes.h

@@ -1,8 +1,19 @@
+#ifndef INCLUDES
+#define INCLUDES
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <time.h>
-#include "general.h"
+#endif
+
+#ifndef NEURON_H
 #include "neuron.h"
-#include "neuron.c"
-#include "general.c"
+#endif
+
+#ifndef GENERAL_H
+#include "general.h"
+#endif
+
+#ifndef ODE_H
+#include "ode.h"
+#endif

main.c

@@ -1,5 +1,6 @@
 #include "definitions.h"
 #include "includes.h"
+#include "math_includes.h"
 
 int main(int argc, char** argv)
 {
@@ -13,6 +14,8 @@ int main(int argc, char** argv)
   assoc_network_params(network, params);
   link_neurons(network, input_filename);
 
+  ode_test(network->neurons[0]->state,params);
+
   if(argc > 2) output_state(network, init_neuron_state, params, argv[2]);
   else if(argc < 2) free(input_filename);
   cleanup(network, init_neuron_state, params);

math_includes.h

@@ -0,0 +1,7 @@
+#ifndef MATHINCLUDES
+#define MATHINCLUDES
+#include <math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_odeiv.h>
+#endif

neuron.c

@@ -1,3 +1,6 @@
+#include "definitions.h"
+#include "includes.h"
+
 struct network * create_network(char * filename)
 {
   long i;

neuron.h

@@ -1,3 +1,6 @@
+#ifndef NEURON_H
+#define NEURON_H
+
 struct neuron {
   int num_links;
   struct neuron_state * state;
@@ -17,6 +20,11 @@ struct neuron_params {
   float * values;
 };
 
+struct neuron_ode_params {
+  struct neuron_state * state;
+  struct neuron_params * params;
+};
+
 struct network {
   long size;
   struct neuron ** neurons;
@@ -61,3 +69,5 @@ void create_queued_links(struct network * network, struct link_queue * link_queu
 void output_state(struct network * network, struct neuron_state * state, struct neuron_params * params, char * filename);
 void write_to_file(FILE * fp, char * line);
 void init_nondefault_states(struct network * network, char * filename);
+
+#endif

ode.c

@@ -0,0 +1,107 @@
+#include "definitions.h"
+#include "includes.h"
+#include "math_includes.h"
+
+int hh_ode (double t, const double y[], double f[], void *params)
+{
+
+  double V_T, C_m, g_bar_Na, g_bar_K, g_bar_L, E_Na, E_K, E_L, A_alpha, B_alpha, A_beta, B_beta;
+  double alpha_n, beta_n, alpha_m, beta_m, alpha_h, beta_h;
+  double i_m, A, I_e;
+  double tau_n, tau_m, tau_h, n_inf, m_inf, h_inf;
+  struct neuron_ode_params * state = (struct neuron_ode_params *)params;
+  // f[0] - V', f[1] - n', f[2] - m', f[3] - h'
+
+  V_T = state->params->values[0];
+  C_m = state->params->values[1];
+  g_bar_Na = state->params->values[2];
+  g_bar_K = state->params->values[3];
+  g_bar_L = state->params->values[4];
+  E_Na = state->params->values[5];
+  E_K = state->params->values[6];
+  E_L = state->params->values[7];
+  A_alpha = state->params->values[8];
+  B_alpha = state->params->values[9];
+  A_beta = state->params->values[10];
+  B_beta = state->params->values[11];
+  A = state->params->values[12];
+  I_e = state->params->values[13];
+
+  //alpha_n = A_alpha*exp(-1.0*B_alpha*y[0]/V_T);
+  //beta_n = A_beta*exp(-1.0*B_beta*y[0]/V_T);
+  alpha_n = 0.01*(y[0] + 55)/(1 - exp(-0.1*(y[0]+55.0)));
+  beta_n = 0.125*exp(-0.0125*(y[0]+65.0));
+  n_inf = alpha_n / (alpha_n + beta_n);
+  tau_n = 1.0/(alpha_n + beta_n);
+
+  alpha_m = 0.1*(y[0] + 40.0) / (1.0 - exp(-0.1*(y[0]+40.0)));
+  beta_m = 4.0*exp(-0.0556*(y[0]+65.0));
+  m_inf = alpha_m / (alpha_m + beta_m);
+  tau_m = 1.0/(alpha_m + beta_m);
+
+  alpha_h = 0.07*exp(-0.05*(y[0]+65.0));
+  beta_h = 1.0 / (1.0 + exp(-0.1*(y[0]+35.0)));
+  h_inf = alpha_h / (alpha_h + beta_h);
+  tau_h = 1.0/(alpha_h + beta_h);
+
+  i_m = g_bar_L*(y[0] - E_L) + g_bar_K*pow(y[1],4.0)*(y[0] - E_K) + g_bar_Na*pow(y[2],3)*y[3]*(y[0] - E_Na);
+
+  //f[0] = (1.0/c_m) * (-1.0*i_m + I_e/A);
+  f[0] = I_e + -1.0*i_m/C_m;
+  f[1] = (n_inf - y[1])/(tau_n);
+  f[2] = (m_inf - y[2])/(tau_m);
+  f[3] = (h_inf - y[3])/(tau_h);
+  /*
+  f[1] = alpha_n*(1.0-y[1]) - beta_n*y[1];
+  f[2] = alpha_m*(1-y[2]) - beta_m*y[2];
+  f[3] = alpha_h*(1.0-y[3]) - beta_h*y[3];
+  */
+  //printf("%.5f %.5f %.5f %.5f\n",f[0],f[1],f[2],f[3]);
+
+  return GSL_SUCCESS;
+}
+
+int hh_ode_jac (double t, const double y[], double *dfdy, double dfdt[], void *params)
+{
+  return GSL_SUCCESS;
+}
+
+int ode_test(struct neuron_state * neuron_state, struct neuron_params * neuron_params)
+{
+  int i = 0;
+  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;
+
+  gsl_odeiv_step * s = gsl_odeiv_step_alloc(T, 4);
+  gsl_odeiv_control * c = gsl_odeiv_control_y_new(1e-6, 0.0);
+  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc(4);
+
+  struct neuron_ode_params * params = (struct neuron_ode_params *)malloc(sizeof(struct neuron_ode_params));
+  params->state = neuron_state;
+  params->params = neuron_params;
+
+  // params: function, jacobian, dimension, void * params
+  gsl_odeiv_system sys = {hh_ode, NULL, 4, params};
+
+  double t = 0.0, t1 = 100.0;
+  double h = 1e-6;
+  double y[4] = { neuron_state->values[0], neuron_state->values[1], neuron_state->values[2], neuron_state->values[3]};
+
+  while (t < t1)
+    {
+      int status = gsl_odeiv_evolve_apply (e, c, s,
+					   &sys, 
+					   &t, t1,
+					   &h, y);
+
+      if (status != GSL_SUCCESS)
+	break;
+      if(i++ % 1 == 0)
+	printf("%.5e, %.5e, %.5e, %.5e, %.5e\n", t, y[0], y[1], y[2], y[3]);
+    }
+
+  gsl_odeiv_evolve_free (e);
+  gsl_odeiv_control_free (c);
+  gsl_odeiv_step_free (s);
+  free(params);
+  return 0;
+}

ode.h

@@ -0,0 +1,18 @@
+#ifndef DEFINITIONS
+#define DEFINITIONS
+#include "definitions.h"
+#endif
+
+#ifndef INCLUDES
+#define INCLUDES
+#include "includes.h"
+#endif
+
+#ifndef ODE_H
+#define ODE_H
+
+int ode_test(struct neuron_state * neuron_state, struct neuron_params * neuron_params);
+int hh_ode_jac (double t, const double y[], double *dfdy, double dfdt[], void *params);
+int hh_ode (double t, const double y[], double f[], void *params);
+
+#endif

template.txt

@@ -2,25 +2,41 @@
 2
 
 @PARAMS
-7
-a-ltp 0
-a-ltd 0
-theta-minus -41
-theta-plus -38
-tau-x 16
-tau-minus 10
-tau-plus 7
+14
+V_T 25
+C_m 1.0
+g-bar_Na 120.0
+g-bar_K 36.0
+g-bar_L 0.3
+E_Na 50.0
+E_K -77.0
+E_L -54.402
+A_alpha 1.22
+B_alpha -1.0
+A_beta 0.056
+B_beta 0.3125
+A 0.1e-5
+I_e 60.0
+
+c_m = 
+A =  0.00001 to 0.0001
+V_t = 24 to 27, say 25
+
+C_m: 0.1 to 1 nF
+R_m: 10 to 100 M\omega
+n 0.417938
+m 0.109496
+h 0.596121
 
 @DEFAULT_STATE
-5
-u -60
-u-minus-bar -60
-u-plus-bar -60
-x 0
-x-bar 0
+4
+V -65.0
+n 0.3177
+m 0.0529
+h 0.5981
 
 @INIT_STATES
-1 1 2 3 4 5
+
 
 @LINKS
 defined