spga: /home/kevin/Documents/PYLikeC/model_slow

00001 #define _GNU_SOURCE
00002 #include <stdio.h>
00003 #include <math.h>
00004 #include <gsl/gsl_errno.h>
00005 #include <gsl/gsl_odeiv.h>
00006 #include <time.h>
00007 #include "currents.h"
00008 #include "utility.h"
00009 
00010 #define KHH_SILLY_MODEL 65521
00011 
00012 int pylike(double t1, const double y[], double yp[], void * comm);
00013 
00014 size_t model_slow_wave(struct input_data *model_data_ptr, double ** wave_ptr)
00015 {
00016         int neq;
00017         double h, tnow, tend;
00018         double  tol;
00019         int i;
00020         double thres[NEQ], ynow[NEQ], ypnow[NEQ];
00021         clock_t start, now;
00022         double cpu_time_used;
00023         int dt_count = 0;
00024 
00025         double * wave;
00026         size_t sample_buffer_length;
00027         size_t samples_taken = 0;       
00028         const double sample_interval = 0.001;
00029         double last_sample_t;
00030         
00031         const gsl_odeiv_step_type * T 
00032                 = gsl_odeiv_step_rk8pd;
00033         
00034         gsl_odeiv_step * s 
00035                 = gsl_odeiv_step_alloc (T, NEQ);
00036         gsl_odeiv_control * c 
00037                 = gsl_odeiv_control_y_new (1e-6, 1e-6);
00038         gsl_odeiv_evolve * e 
00039                 = gsl_odeiv_evolve_alloc (NEQ);
00040         
00041         gsl_odeiv_system sys = {pylike, NULL, NEQ, model_data_ptr};
00042         
00043 #ifdef SINGLE_BUILD
00044         FILE *wavefile;
00045         
00046         wavefile = fopen("output_slow_wave.txt", "w");
00047 #endif
00048         
00049         
00050         neq = NEQ;
00051         
00052         /* INITIAL CONDITIONS */
00053         tnow     =-100;       /* s */
00054         ynow[0]  = -55;       /* v */
00055         ynow[1]  =   0.926;   /* na_h */
00056         ynow[2]  =   0.146;   /* kd_n */
00057         ynow[3]  =   0.237;   /* a_bfast */
00058         ynow[4]  =   0.309;   /* a_a */
00059         ynow[5]  =   0.237;   /* a_bslow */
00060         ynow[6]  =   0.105;   /* h_r */
00061         ynow[7]  =   0.146;   /* kv3_n */
00062         ynow[8]  =   0.2;     /* kv3_h */
00063         ynow[9]  =   0.05;    /* ca_i */
00064         ynow[10] =   0.183;   /* k_a */
00065         ynow[11] =   6.77e-7; /* k_ca_a */
00066         ynow[12] =   0.923;   /* k_ca_b */
00067         ynow[13] =   0.146;   /* ahp_a */
00068         ynow[14] =   0.00185; /* ca_f1_a */
00069         ynow[15] =   0.651;   /* ca_f1_b */
00070         ynow[16] =   1.65e-5; /* ca_f2_a */
00071         ynow[17] =   0.18;    /* ca_s_a */
00072         
00073         tend = 691;
00074         for (i = 0; i < neq; i++)
00075                 thres[i] = 1.0e-4;
00076         
00077         h = 1e-12;
00078         tol = 1.0e-3;
00079         
00080         last_sample_t = tnow;
00081         
00082         sample_buffer_length = 1024;
00083         wave = xmalloc(sample_buffer_length * sizeof(double));
00084         
00085         start = clock();
00086         while (tnow < tend)
00087         {
00088                 /* Apply a time step. */
00089                 int status = gsl_odeiv_evolve_apply (e, c, s,
00090                                 &sys,
00091                                 &tnow, tend,
00092                                 &h, ynow);
00093 
00094                 /* Check for NAN derivatives */
00095                 if (status == KHH_SILLY_MODEL)
00096                 {
00097                         tnow = tnow - h;
00098                         h = h / 10.0;
00099                 }
00100                 /* Check for some unknown failure */
00101                 else if (status != GSL_SUCCESS)
00102                 {
00103                         fprintf(stderr, "\n\nMODEL FAILURE %i at t = %f\n\n",  status, tnow);
00104                         samples_taken = 0;
00105                         tnow = tend + 1;
00106                 }
00107                 /* We apear to have made a good step. */
00108                 else
00109                 {
00110                         /* Is it time to take a sample? */
00111                         while ( tnow >= last_sample_t + sample_interval )
00112                         {
00113                                 last_sample_t += sample_interval;
00114 #ifdef SINGLE_BUILD
00115                                 fprintf(wavefile, "%a %a\n", last_sample_t, ynow[0]);
00116 #endif            
00117                                 wave[samples_taken] = ynow[0];
00118                                 samples_taken++;
00119                                 if ( samples_taken >= sample_buffer_length )
00120                                 {
00121                                         sample_buffer_length += 1024;
00122                                         wave = xrealloc(wave, sample_buffer_length * sizeof(double));
00123                                 }
00124                                 
00125                         }
00126                 }
00127                 
00128                 dt_count++;
00129                 if ((dt_count % 1000) == 0);
00130                 {
00131                         now = clock();
00132                         cpu_time_used = ((double) (now - start)) / CLOCKS_PER_SEC;
00133                         if (cpu_time_used > 300)
00134                         {
00135                                 fprintf(stderr, "\nslow death %f%%\n", tnow / tend * 100);
00136                                 for (i = 0; i < NUM_PARAMS; i++)
00137                                         fprintf(stderr, " %a", model_data_ptr->params[i]);
00138                                 fprintf(stderr, "\n");
00139                                 samples_taken = 0;
00140                                 tnow = tend + 1;
00141                         }
00142                 }  
00143         }
00144         
00145         
00146         gsl_odeiv_evolve_free(e);
00147         gsl_odeiv_control_free(c);
00148         gsl_odeiv_step_free(s);
00149         
00150 #ifdef SINGLE_BUILD
00151         fclose(wavefile);
00152 #endif
00153         
00154         *wave_ptr = wave;
00155         
00156         return samples_taken;
00157 }
00158 int pylike(double t1, const double y[], double yp[], void *comm)
00159 {
00160    /* PARAMETERS */
00161    double vr_leak = -50;     /* mv */
00162    double vr_na = 50;        /* mv */
00163    double vr_k = -80;        /* mv */
00164    double vr_h = -10;        /* mv */
00165    double vr_ca;             /* mv */ /* calculated later */
00166    double capacitance = 1.7; /* nF */
00167    double ms_per_s = 1000;   /* ms / s */
00168    double stim_amp = -4;     /* nA */ 
00169    
00170    /*  MAXIMUM CONDUCTANCES */
00171    double g_leak;  /* =    0.1;   */ /* muS */
00172    double g_na;    /* = 2300;     */ /* muS */
00173    double g_kd;    /* =    0.59;  */ /* muS */
00174    double g_afast; /* =    0.21;  */ /* muS */
00175    double g_aslow; /* =    0.047; */ /* muS */
00176    double g_h;     /* =    0.037; */ /* muS */
00177    double g_kv3;   /* =    0.08;  */ /* muS */
00178    double g_k;     /* =    1.2;   */ /* muS */
00179    double g_k_ca;  /* =    3.2;   */ /* muS */
00180    double g_ahp;   /* =    0.001; */ /* muS */
00181    double g_ca_f1; /* =    0.1;   */ /* muS */
00182    double g_ca_f2; /* =    0.001; */ /* muS */
00183    double g_ca_s;  /* =    0.008; */ /* muS */
00184    
00185    
00186    /* VARIABLES */
00187    double v;                     /* mv */
00188    double na_m, na_h;            /* unitless */ /* na_m will be instantaneous */
00189    double kd_n;                  /* unitless */
00190    double a_bfast, a_a, a_bslow; /* unitless */
00191    double h_r;                   /* unitless */
00192    double kv3_n, kv3_h;          /* unitless */
00193    double ca_i;                  /* muM */
00194    double k_a;                   /* unitless */
00195    double k_ca_a, k_ca_b;        /* unitless */
00196    double ahp_a;                 /* unitless */
00197    double ca_f1_a;               /* unitless */
00198    double ca_f1_b;               /* unitless */
00199    double ca_f2_a;               /* unitless */
00200    double ca_s_a;                /* unitless */
00201    
00202    /* CURRENTS */
00203    double i_leak; /* nA */
00204    double i_na;   /* nA */
00205    double i_kd;   /* nA */
00206    double i_a;    /* nA */
00207    double i_h;    /* nA */
00208    double i_kv3;  /* nA */
00209    double i_k;    /* nA */
00210    double i_k_ca; /* nA */
00211    double i_ahp;  /* nA */
00212    double i_ca_f; /* nA */
00213    double i_ca_s; /* nA */
00214    double i_ext;  /* nA */
00215    
00216    
00217    struct input_data *model_data_ptr; 
00218    
00219    /* model_data_ptr is a pointer to a structure of type input_data*/
00220    
00221    int i;
00222    int status = GSL_SUCCESS;
00223    
00224    model_data_ptr = (struct input_data *)comm;
00225    
00226    /* pointy structures are cast at one another*/
00227    
00228    v       = y[0];        /* mv */
00229    na_m    = na_m_inf(v); /* unitless */
00230    na_h    = y[1];        /* unitless */
00231    kd_n    = y[2];        /* unitless */
00232    a_bfast = y[3];        /* unitless */
00233    a_a     = y[4];        /* unitless */
00234    a_bslow = y[5];        /* unitless */
00235    h_r     = y[6];        /* unitless */
00236    kv3_n   = y[7];        /* unitless */
00237    kv3_h   = y[8];        /* unitless */
00238    ca_i    = y[9];        /* muM */
00239    k_a     = y[10];       /* unitless */
00240    k_ca_a  = y[11];       /* unitless */
00241    k_ca_b  = y[12];       /* unitless */
00242    ahp_a   = y[13];       /* unitless */
00243    ca_f1_a = y[14];       /* unitless */
00244    ca_f1_b = y[15];       /* unitless */
00245    ca_f2_a = y[16];       /* unitless */
00246    ca_s_a  = y[17];       /* unitless */
00247    
00248    g_leak  = model_data_ptr->params[0];
00249    g_na    = model_data_ptr->params[1];
00250    g_kd    = model_data_ptr->params[2];
00251    g_afast = model_data_ptr->params[3];
00252    g_aslow = model_data_ptr->params[4];
00253    g_h     = model_data_ptr->params[5];
00254    g_kv3   = model_data_ptr->params[6];
00255    g_k     = model_data_ptr->params[7];
00256    g_k_ca  = model_data_ptr->params[8];
00257    g_ahp   = model_data_ptr->params[9];
00258    g_ca_f1 = model_data_ptr->params[10];
00259    g_ca_f2 = model_data_ptr->params[11];
00260    g_ca_s  = model_data_ptr->params[12];
00261    
00262    vr_ca  = ca_reverse(ca_i);                                                   /* mv */
00263    i_leak = (v - vr_leak) *  g_leak;                                            /* nA */
00264    i_na   = (v - vr_na)   *  g_na    * na_m    * na_m  * na_m  * na_h;          /* nA */
00265    i_kd   = (v - vr_k)    *  g_kd    * kd_n    * kd_n  * kd_n  * kd_n;          /* nA */
00266    i_a    = (v - vr_k)    * (g_afast * a_a     * a_a   * a_a   * a_bfast + 
00267                              g_aslow * a_a     * a_a   * a_a   * a_bslow);      /* nA */
00268    i_h    = (v - vr_h)    *  g_h     * h_r;                                     /* nA */
00269    i_kv3  = (v - vr_k)    *  g_kv3   * kv3_n   * kv3_n * kv3_n * kv3_n * kv3_h; /* nA */
00270    i_k    = (v - vr_k)    *  g_k     * k_a     * k_a;                           /* nA */
00271    i_k_ca = (v - vr_k)    *  g_k_ca  * k_ca_a  * k_ca_b;                        /* nA */
00272    i_ahp  = (v - vr_k)    *  g_ahp   * ahp_a   * ahp_a;                         /* nA */
00273    i_ca_f = (v - vr_ca)   * (g_ca_f1 * ca_f1_a * ca_f1_b +
00274                              g_ca_f2 * ca_f2_a);                                /* nA */
00275    i_ca_s = (v - vr_ca)   *  g_ca_s  * ca_s_a;                                  /* nA */
00276    i_ext  = stim_amp * stim(t1, model_data_ptr);                                /* nA    stimulation.c */
00277       
00278    /* DERIVATIVES */
00279    yp[0]  = -(i_leak + i_na + i_kd + i_a + i_h + i_kv3 + i_k + 
00280               i_k_ca + i_ahp + i_ca_f + i_ca_s - i_ext) * ms_per_s / capacitance; /* mV / s  */
00281    if (g_na > 0)
00282      {
00283         yp[1]  = na_h_deriv(v, na_h);               /* 1 / s    fast_sodium.c */
00284      }
00285    else
00286      {
00287         yp[1]  = 0;
00288      }
00289    if (g_kd > 0)
00290      {
00291         yp[2]  = kd_n_deriv(v, kd_n);               /* 1 / s    delayed_rectifier.c */
00292      }
00293    else
00294      {
00295         yp[2]  = 0;
00296      }
00297    if (g_afast > 0)
00298      {
00299         yp[3]  = a_bfast_deriv(v, a_bfast);         /* 1 / s    depol_act_potasium.c */
00300      }
00301    else 
00302      {
00303         yp[3]  = 0;
00304      }
00305    if ((g_afast > 0) || (g_aslow > 0))
00306      {
00307         yp[4]  = a_a_deriv(v, a_a);                 /* 1 / s    depol_act_potasium.c */
00308      }
00309    else 
00310      {
00311         yp[4]  = 0;
00312      }
00313    if (g_aslow > 0)
00314      {
00315         yp[5]  = a_bslow_deriv(v, a_bslow);         /* 1 / s    depol_act_potasium.c */
00316      }
00317    else 
00318      {
00319         yp[5]  = 0;
00320      }
00321    if (g_h > 0)
00322      {
00323         yp[6]  = h_r_deriv(v, h_r);                 /* 1 / s    hyperpol_act_cation.c */
00324      }
00325    else 
00326      {
00327         yp[6]  = 0;
00328      }
00329    if (g_kv3 > 0)
00330      {
00331         yp[7]  = kv3_n_deriv(v, kv3_n);             /* 1 / s    abbot_kv3.c */
00332         yp[8]  = kv3_h_deriv(kv3_h, kv3_n);         /* 1 / s    abbot_kv3.c */
00333      }
00334    else
00335      {
00336         yp[7]  = 0;
00337         yp[8]  = 0;
00338      }
00339    yp[9]  = ca_i_deriv(i_ca_f + i_ca_s, ca_i); /* muM / s  calcium_regulation.c */
00340    if (g_k > 0)
00341      {
00342         yp[10] = k_a_deriv(v, k_a, ca_i);           /* 1 / s    calcium_and_v_dep_k-K.c */
00343      }
00344    else 
00345      {
00346         yp[10] = 0;
00347      }
00348    if (g_k_ca > 0)
00349      {
00350         yp[11] = k_ca_a_deriv(v, ca_i,  k_ca_a);    /* 1 / s    calcium_and_v_dep_k-K_Ca.c */
00351         yp[12] = k_ca_b_deriv(ca_i,  k_ca_b);       /* 1 / s    calcium_and_v_dep_k-K_Ca.c */
00352      }
00353    else 
00354      {
00355         yp[11] = 0;
00356         yp[12] = 0;
00357      }
00358    if (g_ahp > 0)
00359      {
00360         yp[13] = ahp_a_deriv(ca_i, ahp_a);          /* 1 / s    after_hyperpolarization.c */
00361      }
00362    else
00363      {
00364         yp[13] = 0;
00365      }
00366    if (g_ca_f1 > 0)
00367      {
00368         yp[14] = ca_f1_a_deriv(v, ca_f1_a);         /* 1 / s    fast_calcium.c */
00369         yp[15] = ca_f1_b_deriv(v, ca_f1_b);         /* 1 / s    fast_calcium.c */
00370      }
00371    else
00372      {
00373         yp[14] = 0;
00374         yp[15] = 0;
00375      }
00376    if (g_ca_f2 > 0)
00377      {
00378         yp[16] = ca_f2_a_deriv(v, ca_f2_a);         /* 1 / s    fast_calcium.c */
00379      }
00380    else
00381      {
00382         yp[16] = 0;
00383      }
00384    if (g_ca_s > 0)
00385      {
00386         yp[17] = ca_s_a_deriv(v, ca_s_a);           /* 1 / s    slow_calcium.c */
00387      }
00388    else
00389      {
00390         yp[17] = 0;
00391      }
00392 
00393   /*
00394    * fprintf(stderr, "%f ", t1);
00395    * for (i = 0; i < NEQ; i++)
00396    *   fprintf(stderr, "%f %f ", y[i], yp[i]);
00397    * fprintf(stderr, "\n");
00398    */
00399    
00400    for (i = 0; i < NEQ; i++)
00401      if (gsl_finite(yp[i]) == 0)
00402           status = KHH_SILLY_MODEL;
00403    
00404    if (status == KHH_SILLY_MODEL)
00405      for (i = 0; i < NEQ; i++)
00406        yp[i] = 0;
00407    
00408    return status;
00409 }
/home/kevin/Documents/PYLikeC/model_slow_wave.c
