Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #define _USE_MATH_DEFINES
- #include "math.h"
- #include <stdlib.h>
- #include <stdio.h>
- #include <locale.h>
- #include <time.h>
- #include <stdbool.h>
- #define Node_count 5*5*2 // должно быть таким, что бы из этого числа, разделенного на два, извлекался корень
- #define Node_count_half Node_count / 2
- #define Node_wire_width (int)sqrt(Node_count_half)
- #define Equations_per_node 4
- #define Equations_count Node_count * Equations_per_node
- double f[Equations_count];
- double f_diff[Equations_count];
- double C_m = 1; // muF/cm^2
- double g_K = 35; // mS/cm^2
- double g_Na = 40; // mS/cm^2
- double g_L = 0.3; // mS/cm^2
- double E_K = -77; // mV
- double E_Na = 55; // mV
- double E_L = -65; // mV
- double g_syn = 0.08; // 0.1 // 0.04 // 0.2
- double k_syn = 0.2; // 0.2
- double E_syn[Node_count];
- double I_app[Node_count];
- /*const double sigma_G = 0.0;
- const double sigma_GN = 0.0;
- const double sigma_N = 1.0;*/
- double** A;
- double** B;
- double* C;
- double** sigma;
- double tau[Node_count][Node_count];
- #define tau_min 2 // ms
- #define tau_max 12 // ms
- #define ms_to_step 200 // (0.001 / dt) !!! Don't forget !!!
- #define Max_delay tau_max * ms_to_step
- double V_old_array[Node_count][Max_delay];
- const double Freq = 1300; // Hz
- const double Min_magintude = -0.00; // pA
- const double Max_magintude = 0.1; // pA
- const double Duration = 0.001; // sec
- double Meander_start_from_zero[Node_count];
- double Meander_width[Node_count];
- double Meander_height[Node_count];
- double Meander_interval[Node_count];
- double last_meander_end[Node_count];
- double I_stim(int i, double t)
- {
- if (t < Meander_start_from_zero[i])
- return 0;
- t -= Meander_start_from_zero[i];
- t = fmod(t, Meander_width[i] + Meander_interval[i]);
- return t < Meander_width[i] ? Meander_height[i] : 0;
- }
- double V(int i)
- {
- return f[i * 4];
- }
- void SetV(int i, double value)
- {
- f[i * 4] = value;
- }
- double m(int i)
- {
- return f[i * 4 + 1];
- }
- void Setm(int i, double value)
- {
- f[i * 4 + 1] = value;
- }
- double n(int i)
- {
- return f[i * 4 + 2];
- }
- void Setn(int i, double value)
- {
- f[i * 4 + 2] = value;
- }
- double h(int i)
- {
- return f[i * 4 + 3];
- }
- void Seth(int i, double value)
- {
- f[i * 4 + 3] = value;
- }
- double V_old(int i, int delay)
- {
- return V_old_array[i][Max_delay - 1 - delay];
- }
- int RandomI(int min, int max)
- {
- return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
- }
- double RandomD(double min, double max)
- {
- return ((double)rand() / RAND_MAX) * (max - min) + min;
- }
- double alpha_m(double* f, int i)
- {
- return 0.182 * (V(i) + 35) / (1 - exp(-(V(i) + 35) / 9));
- }
- double beta_m(double* f, int i)
- {
- return -0.124 * (V(i) + 35) / (1 - exp((V(i) + 35) / 9));
- }
- double alpha_n(double* f, int i)
- {
- return 0.02 * (V(i) - 25) / (1 - exp(-(V(i) - 25) / 9));
- }
- double beta_n(double* f, int i)
- {
- return -0.002 * (V(i) - 25) / (1 - exp((V(i) - 25) / 9));
- }
- double alpha_h(double* f, int i)
- {
- return 0.25 * exp(-(V(i) + 90) / 12);
- }
- double beta_h(double* f, int i)
- {
- return 0.25 * exp((V(i) + 62) / 6) / exp((V(i) + 90) / 12);
- }
- double HodgkinHuxley(int i, double* f, double t)
- {
- int in = i / 4;
- int il = i % 4;
- switch (il)
- {
- case 0:
- {
- double sum = 0;
- //double ee = 0;
- //double Vold = 0;
- /*for (int j = 0; j < Node_count; j++)
- {
- //sum += A[in][j] * g_syn * (V(in) - V(j));
- sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old(j, tau[in][j]) / k_syn));
- }*/
- /*for (int j = 0; j < C[in]; j++)
- {
- sum += sigma[in][(int)B[in][j]] * (V((int)B[in][j]) - V(in));
- }*/
- for (int j = 0; j < C[in]; j++)
- {
- //sum += g_syn * (V((int)B[in][j]) - V(in)); // устаревшая часть, нужна для проверки разностной схемы
- //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j]) / k_syn)); // устаревшая часть
- sum += A[in][(int)B[in][j]] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j], tau[in][(int)B[in][j]]) / k_syn));
- //Vold = V_old((int)B[in][j], tau[in][(int)B[in][j]]);
- //ee = exp(-V_old((int)B[in][j], tau[in][(int)B[in][j]]) / k_syn);
- //printf("i = %d\t V_old = %f\t exp = %f\n", in, Vold, ee);
- }
- //printf("i = %d\t sum = %f\n", in, sum);
- return 1000 * ((g_Na * pow(m(in), 3) * h(in) * (E_Na - V(in)) + g_K * n(in) * (E_K - V(in)) + g_L * (E_L - V(in)) + I_app[in] + I_stim(in, t) + sum ) / C_m); // V
- }
- case 1:
- {
- return 1000 * (alpha_m(f, in) * (1 - m(in)) - beta_m(f, in) * m(in)); // m
- }
- case 2:
- {
- return 1000 * (alpha_n(f, in) * (1 - n(in)) - beta_n(f, in) * n(in)); // n
- }
- case 3:
- {
- return 1000 * (alpha_h(f, in) * (1 - h(in)) - beta_h(f, in) * h(in)); // h
- }
- }
- return 0;
- }
- void RungeKutta(double t, double dt, double* f, double* f_next)
- {
- double k[Equations_count][4];
- // k1
- for (int i = 0; i < Equations_count; i++)
- k[i][0] = HodgkinHuxley(i, f, t) * dt;
- double phi_k1[Equations_count];
- for (int i = 0; i < Equations_count; i++)
- phi_k1[i] = f[i] + k[i][0] / 2;
- // k2
- for (int i = 0; i < Equations_count; i++)
- k[i][1] = HodgkinHuxley(i, phi_k1, t) * dt;
- double phi_k2[Equations_count];
- for (int i = 0; i < Equations_count; i++)
- phi_k2[i] = f[i] + k[i][1] / 2;
- // k3
- for (int i = 0; i < Equations_count; i++)
- k[i][2] = HodgkinHuxley(i, phi_k2, t) * dt;
- double phi_k3[Equations_count];
- for (int i = 0; i < Equations_count; i++)
- phi_k3[i] = f[i] + k[i][2] / 2;
- // k4
- for (int i = 0; i < Equations_count; i++)
- k[i][3] = HodgkinHuxley(i, phi_k3, t) * dt;
- for (int i = 0; i < Equations_count; i++)
- f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
- }
- void CopyArray(double* source, double* target, int N)
- {
- for (int i = 0; i < N; i++)
- target[i] = source[i];
- }
- bool Approximately(double a, double b)
- {
- if (a < 0)
- a = -a;
- if (b < 0)
- b = -b;
- return a - b <= 0.000001;
- }
- // http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/
- double nextTime(double rateParameter)
- {
- return -log(1.0 - (double)rand() / (RAND_MAX)) / rateParameter;
- }
- void GenerateRandomMeander(int i, double min_start_time)
- {
- double offset = nextTime(Freq);
- if (offset < 0)
- {
- int a = 0;
- }
- Meander_start_from_zero[i] = min_start_time + offset;
- Meander_width[i] = Duration;
- Meander_height[i] = RandomD(Min_magintude, Max_magintude);
- }
- void FillAMatrixZero()
- {
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 0; j < Node_count; j++)
- {
- A[i][j] = 0;
- }
- }
- }
- void FillBCMatrix()
- {
- for (int i = 0; i < Node_count; i++)
- {
- int bIndex = 0;
- C[i] = 0;
- for (int j = 0; j < Node_count; j++)
- {
- if (A[i][j] == 1)
- {
- B[i][bIndex] = j;
- bIndex++;
- C[i]++;
- }
- }
- }
- }
- /*void FillSigmaMatrix()
- {
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 0; j < Node_count; j++)
- {
- if ((i >= 0 && i < Node_count_half) && (j >= 0 || j < Node_count_half))
- {
- sigma[i][j] = sigma_G;
- }
- if ((((i >= Node_count_half && i < Node_count) && (j >= 0 && j < Node_count_half)) || ((i >= 0 && i < Node_count_half) && (j >= Node_count_half && j < Node_count))))
- {
- sigma[i][j] = sigma_GN;
- }
- if ((i >= Node_count_half && i < Node_count) && (j >= Node_count_half && j < Node_count))
- {
- sigma[i][j] = sigma_N;
- }
- }
- }
- }*/
- void FillRandomMatrix()
- {
- double p_links = 0.2; // 0.2
- //for (int k = 0; k < 2 * Node_count_half; k++)
- for (int k = 0; k < p_links * Node_count_half * (Node_count_half - 1); k++)
- {
- int i, j;
- do
- {
- i = RandomI(Node_count_half, Node_count);
- j = RandomI(Node_count_half, Node_count);
- } while ((i == j) || (A[i][j] == 1));
- A[i][j] = 1;
- //A[j][i] = 1;
- }
- }
- void FillVOldFromCurrent()
- {
- for (int i = 0; i < Node_count; i++)
- for (int j = 0; j < Max_delay; j++)
- V_old_array[i][j] = V(i);
- }
- void UpdateVOld()
- {
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 1; j < Max_delay; j++)
- V_old_array[i][j - 1] = V_old_array[i][j];
- V_old_array[i][Max_delay - 1] = V(i);
- }
- }
- /*void FillFullTauMatrix()
- {
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 0; j < Node_count; j++)
- {
- if (i < Node_count_half || j < Node_count_half)
- {
- tau[i][j] = 0;
- continue;
- }
- int i_neuron = i - Node_count_half;
- int j_neuron = j - Node_count_half;
- int i_wire_x = i_neuron / Node_wire_width;
- int i_wire_y = i_neuron % Node_wire_width;
- int j_wire_x = j_neuron / Node_wire_width;
- int j_wire_y = j_neuron % Node_wire_width;
- double distance_max = sqrt(2.) * (Node_wire_width - 1);
- double distance = sqrt((i_wire_x - j_wire_x) * (i_wire_x - j_wire_x) + (i_wire_y - j_wire_y) * (i_wire_y - j_wire_y));
- tau[i][j] = (tau_min + distance / (distance_max) * (tau_max - tau_min)) * ms_to_step;
- }
- }
- }*/
- void FillTauMatrix()
- {
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 0; j < Node_count; j++)
- {
- if (i < Node_count_half || j < Node_count_half || i == j || A[i][j] == 0)
- {
- tau[i][j] = 0;
- continue;
- }
- int i_neuron = i - Node_count_half;
- int j_neuron = j - Node_count_half;
- int i_wire_x = i_neuron / Node_wire_width;
- int i_wire_y = i_neuron % Node_wire_width;
- int j_wire_x = j_neuron / Node_wire_width;
- int j_wire_y = j_neuron % Node_wire_width;
- double distance_max = sqrt(2.) * (Node_wire_width - 1);
- double distance = sqrt((i_wire_x - j_wire_x) * (i_wire_x - j_wire_x) + (i_wire_y - j_wire_y) * (i_wire_y - j_wire_y));
- double t = (distance - 1) / (distance_max - 1);
- tau[i][j] = (tau_min + t * (tau_max - tau_min)) * ms_to_step;
- }
- }
- }
- int main(int argc, char *argv[])
- {
- FILE *fp0;
- //FILE *fp_Ca;
- //FILE *fp_IP3;
- //FILE *fp_z;
- //FILE *fp_G;
- FILE *fp_I_stim;
- FILE *fp_V;
- FILE *fp_m;
- FILE *fp_n;
- FILE *fp_h;
- FILE *fp_V_spikes;
- srand(time(NULL));
- //for (int i = 0; i < Node_count; i++)
- // V_old_length[i] = 0;
- A = malloc(Node_count * sizeof(double));
- for (int i = 0; i < Node_count; i++)
- A[i] = malloc(Node_count * sizeof(double));
- B = malloc(Node_count * sizeof(double));
- for (int i = 0; i < Node_count; i++)
- B[i] = malloc(Node_count * sizeof(double));
- C = malloc(Node_count * sizeof(double));
- sigma = malloc(Node_count * sizeof(double));
- for (int i = 0; i < Node_count; i++)
- sigma[i] = malloc(Node_count * sizeof(double));
- FillAMatrixZero();
- FillRandomMatrix();
- FillBCMatrix();
- //FillSigmaMatrix();
- FillTauMatrix();
- fp0 = fopen("A.txt", "w+");
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 0; j < Node_count; j++)
- {
- fprintf(fp0, "%d\t", (int)A[i][j]);
- }
- fprintf(fp0, "\n");
- }
- fclose(fp0);
- fp0 = fopen("tau.txt", "w+");
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 0; j < Node_count; j++)
- {
- fprintf(fp0, "%f\t", tau[i][j] / ms_to_step);
- }
- fprintf(fp0, "\n");
- }
- fclose(fp0);
- //Пишем в файл число связей у каждого осциллятора в четвертом квадранте
- fp0 = fopen("links.txt", "w+");
- for (int i = Node_count_half; i < Node_count; i++)
- {
- int links_count = 0;
- for (int j = Node_count_half; j < Node_count; j++)
- {
- if (A[i][j] == 1)
- {
- links_count++;
- }
- }
- fprintf(fp0, "%d\n", (int)links_count);
- }
- fclose(fp0);
- //setlocale(LC_NUMERIC, "French_Canada.1252");
- fp0 = fopen("test_Poisson.txt", "w+");
- for (int i = 0; i < 10000; i++)
- fprintf(fp0, "%f\n", nextTime(Freq));
- fclose(fp0);
- fp0 = fopen("B.txt", "w+");
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 0; j < C[i]; j++)
- {
- fprintf(fp0, "%d\t", (int)B[i][j]);
- }
- fprintf(fp0, "\n");
- }
- fclose(fp0);
- fp0 = fopen("C.txt", "w+");
- for (int i = 0; i < Node_count; i++)
- {
- fprintf(fp0, "%d\n", (int)C[i]);
- }
- fclose(fp0);
- /*fp0 = fopen("sigma.txt", "w+");
- for (int i = 0; i < Node_count; i++)
- {
- for (int j = 0; j < Node_count; j++)
- {
- fprintf(fp0, "%f\t", sigma[i][j]);
- }
- fprintf(fp0, "\n");
- }
- fclose(fp0);*/
- // Initial values
- /*for (int i = 0; i < Equations_count; i++)
- {
- f[i] = 0;
- }
- for (int i = 0; i < Equations_count; i++)
- {
- I_app[i] = RandomD(9, 40);
- }*/
- double percent_stable_state = 0.40; // 0.40
- double V0 = -58.7085;
- double m0 = 0.0953;
- double n0 = 0.000913;
- double h0 = 0.3662;
- double V1 = 14.8409;
- double m1 = 0.9174;
- double n1 = 0.0140;
- double h1 = 0.0539;
- for (int i = 0; i < Node_count; i++) // init array for all neurons
- {
- SetV(i, 0); // V
- Setm(i, 0); // m
- Setn(i, 0); // n
- Seth(i, 0); // h
- }
- for (int i = Node_count_half; i < Node_count; i++) // init only neuron nodes
- {
- double random = RandomD(0, 1);
- SetV(i, random < percent_stable_state ? V0 : V1); // V
- Setm(i, random < percent_stable_state ? m0 : m1); // m
- Setn(i, random < percent_stable_state ? n0 : n1); // n
- Seth(i, random < percent_stable_state ? h0 : h1); // h
- }
- for (int i = 0; i < Node_count_half; i++)
- {
- I_app[i] = 0; // init for 1st half array
- }
- for (int i = Node_count_half; i < Node_count; i++)
- {
- I_app[i] = 0.8; // init for neurons; Bifurcation point: I_app = 0.82; I_app_up = 1.04
- }
- double percent_excitable = 0.8; // 0.8
- double E_syn0 = 0;
- double E_syn1 = -90;
- for (int i = 0; i < Node_count; i++)
- {
- E_syn[i] = 0; // init for neurons all array
- }
- for (int i = Node_count_half; i < Node_count; i++)
- {
- double random = RandomD(0, 1);
- E_syn[i] = random < percent_excitable ? E_syn0 : E_syn1;
- //printf("i = %d\t E_syn = %f\n", i, E_syn[i]);
- }
- for (int i = 0; i < Node_count; i++)
- {
- GenerateRandomMeander(i, 0);
- last_meander_end[i] = Meander_start_from_zero[i] + Duration;
- }
- const double t_start = 0;
- const double t_max = 40.0; // 100 msec = 0.1 sec
- const double dt = 0.000005; // 0.01 msec = 0.00001 sec; 1 msec = 0.001 sec
- double t = t_start;
- //fp0 = fopen("results.txt", "w+");
- //setlocale(LC_NUMERIC, "French_Canada.1252");
- clock_t start_rk4, end_rk4;
- start_rk4 = clock();
- int lastPercent = -1;
- FillVOldFromCurrent();
- //fp_Ca = fopen("results_Ca.txt", "w+");
- //fp_IP3 = fopen("results_IP3.txt", "w+");
- //fp_z = fopen("results_z.txt", "w+");
- //fp_G = fopen("results_G.txt", "w+");
- fp_I_stim = fopen("results_I_stim.txt", "w+");
- fp_V = fopen("results_V.txt", "w+");
- fp_m = fopen("results_m.txt", "w+");
- fp_n = fopen("results_n.txt", "w+");
- fp_h = fopen("results_h.txt", "w+");
- fp_V_spikes = fopen("results_V_spikes.txt", "w+");
- while (t < t_max || Approximately(t, t_max))
- {
- //fprintf(fp_Ca, "%f\t", t);
- //fprintf(fp_IP3, "%f\t", t);
- //fprintf(fp_z, "%f\t", t);
- //fprintf(fp_G, "%f\t", t);
- fprintf(fp_I_stim, "%f\t", t);
- fprintf(fp_V, "%f\t", t);
- fprintf(fp_m, "%f\t", t);
- fprintf(fp_n, "%f\t", t);
- fprintf(fp_h, "%f\t", t);
- fprintf(fp_V_spikes, "%f\t", t);
- for (int i = 0; i < Node_count; i++)
- {
- if (t > last_meander_end[i])
- {
- GenerateRandomMeander(i, t);
- last_meander_end[i] = Meander_start_from_zero[i] + Duration;
- }
- fprintf(fp_I_stim, "%f\t", I_stim(i, t));
- }
- fprintf(fp_I_stim, "\n");
- // Зашито на будущее
- //for (int i = 0; i < Equations_count; i += 8)
- // fprintf(fp_Ca, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // Ca
- //for (int i = 1; i < Equations_count; i += 8)
- // fprintf(fp_IP3, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // IP3
- //for (int i = 2; i < Equations_count; i += 8)
- // fprintf(fp_z, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // z
- //for (int i = 3; i < Equations_count; i += 8)
- // fprintf(fp_G, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // G
- for (int i = 0; i < Equations_count; i += 4)
- fprintf(fp_V, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V
- for (int i = 1; i < Equations_count; i += 4)
- fprintf(fp_m, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // m
- for (int i = 2; i < Equations_count; i += 4)
- fprintf(fp_n, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // n
- for (int i = 3; i < Equations_count; i += 4)
- fprintf(fp_h, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // h
- //fprintf(fp_Ca, "\n");
- //fprintf(fp_IP3, "\n");
- //fprintf(fp_z, "\n");
- //fprintf(fp_G, "\n");
- fprintf(fp_V, "\n");
- fprintf(fp_m, "\n");
- fprintf(fp_n, "\n");
- fprintf(fp_h, "\n");
- double f_next[Equations_count];
- RungeKutta(t, dt, f, f_next);
- for (int i = 0; i < Equations_count; i += 4)
- {
- double diff = f_next[i] - f[i];
- fprintf(fp_V_spikes, i == Equations_count - 1 ? "%d" : "%d\t", diff < 0 && f_diff[i] > 0 && f[i] > 0 ? 1 : 0);
- f_diff[i] = diff;
- }
- fprintf(fp_V_spikes, "\n");
- CopyArray(f_next, f, Equations_count);
- t += dt;
- int percent = (int)(100 * (t - t_start) / (t_max - t_start));
- if (percent != lastPercent)
- {
- printf("Progress: %d%%\n", percent);
- lastPercent = percent;
- }
- //printf("V(24) = %f\t V_old(24) = %f\n", f[24*4], V_old(24));
- UpdateVOld();
- }
- //fclose(fp_Ca);
- //fclose(fp_IP3);
- //fclose(fp_z);
- //fclose(fp_G);
- fclose(fp_I_stim);
- fclose(fp_V);
- fclose(fp_m);
- fclose(fp_n);
- fclose(fp_h);
- fclose(fp_V_spikes);
- end_rk4 = clock();
- double extime_rk4 = (double)(end_rk4 - start_rk4) / CLOCKS_PER_SEC;
- int minutes = (int)extime_rk4 / 60;
- int seconds = (int)extime_rk4 % 60;
- printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);
- fp0 = fopen("time_exec.txt", "w+");
- fprintf(fp0, "%f\n", extime_rk4);
- fclose(fp0);
- for (int i = 0; i < Node_count; i++)
- free(A[i]);
- free(A);
- for (int i = 0; i < Node_count; i++)
- free(B[i]);
- free(B);
- for (int i = 0; i < Node_count; i++)
- free(sigma[i]);
- free(sigma);
- free(C);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement