API tools faq
paste
SpaceQuester

Untitled

SpaceQuester
Aug 3rd, 2020
288
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 45.17 KB | None | 0 0
raw download clone embed print report
  1. #define _CRT_SECURE_NO_WARNINGS
  2.  
  3. #define _USE_MATH_DEFINES
  4. #include "math.h"
  5. #include <stdlib.h>
  6. #include <stdio.h>
  7. #include <locale.h>
  8. #include <time.h>
  9. #include <stdbool.h>
  10. #include <list>
  11. #include <omp.h>
  12.  
  13. using namespace std;
  14.  
  15. #define Node_count 200
  16.  
  17. #define Equations_per_node 12 // !!! Don't change !!!
  18. #define Equations_count Node_count * Equations_per_node
  19.  
  20. int MaxDeep = 50;
  21.  
  22. double* f;
  23. double* f_diff;
  24.  
  25. double** k;
  26. double* phi_k1;
  27. double* phi_k2;
  28. double* phi_k3;
  29.  
  30. bool enable_I_syn_out = false;
  31.  
  32. double c_0 = 2; // uM
  33. double c_1 = 0.185;
  34. double v_1 = 6; // s^-1
  35. double v_2 = 0.11; // s^-1
  36. double v_3 = 2.2; // uM/s
  37. double* v_4; // uM/s - Controling parameter // 0.5 //double v_4[Node_count]; // uM/s - Controling parameter //0.495
  38. double v_5 = 0.025; // uM/s
  39. double v_6 = 0.2; // uM/s
  40. double k_1 = 0.5; // s^-1
  41. double k_2 = 1; // uM
  42. double k_3 = 0.1;
  43. double k_4 = 1.1; // uM/s
  44. double a_2 = 0.14; // uM/s
  45. double d_1 = 0.13; // uM
  46. double d_2 = 1.049; // uM
  47. double d_3 = 0.9434; // uM
  48. double d_5 = 0.082; // uM
  49. double alpha = 0.8;
  50. double tau_IP3 = 7.143; // s
  51. double IP3_star = 0.16; // uM
  52. double d_Ca = 0.001; // 0.001
  53. double d_IP3 = 0.2; // 0.12
  54. double alpha_Glu = 2; // 2
  55. double g_astro; // 3
  56.  
  57. // https://neuronaldynamics.epfl.ch/online/Ch2.S2.html
  58. double C_m = 1; // muF/cm^2
  59. double g_K = 35; // mS/cm^2
  60. double g_Na = 40; // mS/cm^2
  61. double g_L = 0.3; // mS/cm^2
  62. double E_K = -77; // mV
  63. double E_Na = 55; // mV
  64. double E_L = -65; // mV
  65.  
  66. double C_m_P = 1; // muF/cm^2
  67. double g_K_P = 35; // mS/cm^2
  68. double g_Na_P = 40; // mS/cm^2
  69. double g_L_P = 0.3; // mS/cm^2
  70. double E_K_P = -77; // mV
  71. double E_Na_P = 55; // mV
  72. double E_L_P = -65; // mV
  73.  
  74. double p_rewir;
  75. double p_inhib;
  76.  
  77. double I_app_min;
  78. double I_app_max;
  79.  
  80. double g_syn;
  81. double k_syn = 0.2;
  82. double* E_syn;
  83.  
  84. double g_syn_P;
  85. double k_syn_P = 0.2;
  86. double* E_syn_P;
  87.  
  88. double alpha_G_P = 25; //s^-1
  89. double beta_G_P = 500; //s^-1
  90.  
  91. double* I_app;
  92. double* I_app_P;
  93.  
  94. double** A_A;
  95. double** B_A;
  96. double* C_A;
  97.  
  98. double** A_N;
  99. double** B_N;
  100. double* C_N;
  101.  
  102. double** A_N_P;
  103. double** B_N_P;
  104. double* C_N_P;
  105.  
  106. list<double>* V_spikes;
  107. list<double>* V_spikes_Freq;
  108.  
  109. FILE* fp_I_syn;
  110.  
  111. double** tau;
  112.  
  113. #define tau_min 2 // ms
  114. #define tau_max 12 // ms
  115.  
  116. #define ms_to_step 40 // (0.001 / dt) !!! Don't forget !!!
  117.  
  118. #define Max_delay tau_max * ms_to_step
  119. double** V_old_array;
  120.  
  121. const double Poisson_Freq = 200; // Hz
  122. //const double Min_magintude = -0.20; // muA/cm^2
  123. double Max_magnitude; // muA/cm^2 // 0.20
  124. const double Duration = 0.002; // sec
  125.  
  126. double* Meander_start_from_zero;
  127. double* Meander_width;
  128. double* Meander_height;
  129. double* Meander_interval;
  130. double* last_meander_end;
  131.  
  132. //bool thread_count_printed = false;
  133.  
  134. double I_stim(int i, double t)
  135. {
  136. if (t < Meander_start_from_zero[i])
  137. return 0;
  138.  
  139. t -= Meander_start_from_zero[i];
  140. t = fmod(t, Meander_width[i] + Meander_interval[i]);
  141.  
  142. return t < Meander_width[i] ? Meander_height[i] : 0;
  143. }
  144.  
  145. double Ca(int i)
  146. {
  147. return f[i * Equations_per_node];
  148. }
  149.  
  150. void SetCa(int i, double value)
  151. {
  152. f[i * Equations_per_node] = value;
  153. }
  154.  
  155. double IP3(int i)
  156. {
  157. return f[i * Equations_per_node + 1];
  158. }
  159.  
  160. void SetIP3(int i, double value)
  161. {
  162. f[i * Equations_per_node + 1] = value;
  163. }
  164.  
  165. double z(int i)
  166. {
  167. return f[i * Equations_per_node + 2];
  168. }
  169.  
  170. void Setz(int i, double value)
  171. {
  172. f[i * Equations_per_node + 2] = value;
  173. }
  174.  
  175. double G_P(int i)
  176. {
  177. return f[i * Equations_per_node + 3];
  178. }
  179.  
  180. void SetG_P(int i, double value)
  181. {
  182. f[i * Equations_per_node + 3] = value;
  183. }
  184.  
  185. double V(int i)
  186. {
  187. return f[i * Equations_per_node + 4];
  188. }
  189.  
  190. void SetV(int i, double value)
  191. {
  192. f[i * Equations_per_node + 4] = value;
  193. }
  194.  
  195. double m(int i)
  196. {
  197. return f[i * Equations_per_node + 5];
  198. }
  199.  
  200. void Setm(int i, double value)
  201. {
  202. f[i * Equations_per_node + 5] = value;
  203. }
  204.  
  205. double n(int i)
  206. {
  207. return f[i * Equations_per_node + 6];
  208. }
  209.  
  210. void Setn(int i, double value)
  211. {
  212. f[i * Equations_per_node + 6] = value;
  213. }
  214.  
  215. double h(int i)
  216. {
  217. return f[i * Equations_per_node + 7];
  218. }
  219.  
  220. void Seth(int i, double value)
  221. {
  222. f[i * Equations_per_node + 7] = value;
  223. }
  224.  
  225. ///
  226.  
  227. double V_P(int i)
  228. {
  229. return f[i * Equations_per_node + 8];
  230. }
  231.  
  232. void SetV_P(int i, double value)
  233. {
  234. f[i * Equations_per_node + 8] = value;
  235. }
  236.  
  237. double m_P(int i)
  238. {
  239. return f[i * Equations_per_node + 9];
  240. }
  241.  
  242. void Setm_P(int i, double value)
  243. {
  244. f[i * Equations_per_node + 9] = value;
  245. }
  246.  
  247. double n_P(int i)
  248. {
  249. return f[i * Equations_per_node + 10];
  250. }
  251.  
  252. void Setn_P(int i, double value)
  253. {
  254. f[i * Equations_per_node + 10] = value;
  255. }
  256.  
  257. double h_P(int i)
  258. {
  259. return f[i * Equations_per_node + 11];
  260. }
  261.  
  262. void Seth_P(int i, double value)
  263. {
  264. f[i * Equations_per_node + 11] = value;
  265. }
  266.  
  267. double V_old(int i, int delay)
  268. {
  269. return V_old_array[i][Max_delay - 1 - delay];
  270. }
  271.  
  272. int RandomI(int min, int max)
  273. {
  274. return ((double)rand() / (RAND_MAX - 1)) * (max - min) + min;
  275. }
  276.  
  277. double RandomD(double min, double max)
  278. {
  279. return ((double)rand() / RAND_MAX) * (max - min) + min;
  280. }
  281.  
  282. double J_channel(double* f, int i)
  283. {
  284. return c_1 * v_1 * pow(IP3(i), 3) * pow(Ca(i), 3) * pow(z(i), 3) * (c_0 / c_1 - (1 + 1 / c_1) * Ca(i)) / pow((IP3(i) + d_1) * (Ca(i) + d_5), 3);
  285. }
  286.  
  287. double J_PLC(double* f, int i)
  288. {
  289. return v_4[i] * (Ca(i) + (1 - alpha) * k_4) / (Ca(i) + k_4);
  290. }
  291.  
  292. double J_leak(double* f, int i)
  293. {
  294. return c_1 * v_2 * (c_0 / c_1 - (1 + 1 / c_1) * Ca(i));
  295. }
  296.  
  297. double J_pump(double* f, int i)
  298. {
  299. return v_3 * pow(Ca(i), 2) / (pow(k_3, 2) + pow(Ca(i), 2));
  300. }
  301.  
  302. double J_in(double* f, int i)
  303. {
  304. return v_5 + v_6 * pow(IP3(i), 2) / (pow(k_2, 2) + pow(IP3(i), 2));
  305. }
  306.  
  307. double J_out(double* f, int i)
  308. {
  309. return k_1 * Ca(i);
  310. }
  311.  
  312. double J_Glu(double* f, int i)
  313. {
  314. /*double J = 0;
  315. if (E_syn_P[i] == 0)
  316. {
  317. J += alpha_Glu / (1 + exp(-(G_P(i) - 0.25) / 0.01));
  318. }
  319. return J;*/
  320.  
  321. return alpha_Glu / (1 + exp(-(G_P(i) - 0.25) / 0.01));
  322. }
  323.  
  324. double alpha_m(double* f, int i)
  325. {
  326. return 0.182 * (V(i) + 35) / (1 - exp(-(V(i) + 35) / 9));
  327. }
  328.  
  329. double beta_m(double* f, int i)
  330. {
  331. return -0.124 * (V(i) + 35) / (1 - exp((V(i) + 35) / 9));
  332. }
  333.  
  334. double alpha_n(double* f, int i)
  335. {
  336. return 0.02 * (V(i) - 25) / (1 - exp(-(V(i) - 25) / 9));
  337. }
  338.  
  339. double beta_n(double* f, int i)
  340. {
  341. return -0.002 * (V(i) - 25) / (1 - exp((V(i) - 25) / 9));
  342. }
  343.  
  344. double alpha_h(double* f, int i)
  345. {
  346. return 0.25 * exp(-(V(i) + 90) / 12);
  347. }
  348.  
  349. double beta_h(double* f, int i)
  350. {
  351. return 0.25 * exp((V(i) + 62) / 6) / exp((V(i) + 90) / 12);
  352. }
  353.  
  354. //
  355.  
  356. double alpha_m_P(double* f, int i)
  357. {
  358. return 0.182 * (V_P(i) + 35) / (1 - exp(-(V_P(i) + 35) / 9));
  359. }
  360.  
  361. double beta_m_P(double* f, int i)
  362. {
  363. return -0.124 * (V_P(i) + 35) / (1 - exp((V_P(i) + 35) / 9));
  364. }
  365.  
  366. double alpha_n_P(double* f, int i)
  367. {
  368. return 0.02 * (V_P(i) - 25) / (1 - exp(-(V_P(i) - 25) / 9));
  369. }
  370.  
  371. double beta_n_P(double* f, int i)
  372. {
  373. return -0.002 * (V_P(i) - 25) / (1 - exp((V_P(i) - 25) / 9));
  374. }
  375.  
  376. double alpha_h_P(double* f, int i)
  377. {
  378. return 0.25 * exp(-(V_P(i) + 90) / 12);
  379. }
  380.  
  381. double beta_h_P(double* f, int i)
  382. {
  383. return 0.25 * exp((V_P(i) + 62) / 6) / exp((V_P(i) + 90) / 12);
  384. }
  385.  
  386. double UllahJung_HodgkinHuxley(int i, double* f, double t)
  387. {
  388. int in = i / Equations_per_node;
  389. int il = i % Equations_per_node;
  390.  
  391. switch (il)
  392. {
  393. case 0: // Ca
  394. {
  395. double sum_1 = 0;
  396.  
  397. /*for (int j = 0; j < Node_count; j++)
  398. {
  399. sum_1 += d_Ca * (Ca(j) - Ca(in));
  400. }*/
  401.  
  402. for (int j = 0; j < C_A[in]; j++)
  403. {
  404. sum_1 += d_Ca * (Ca((int)B_A[in][j]) - Ca(in));
  405. }
  406.  
  407. return J_channel(f, in) - J_pump(f, in) + J_leak(f, in) + J_in(f, in) - J_out(f, in) + sum_1;
  408. }
  409.  
  410. case 1: // IP3
  411. {
  412. double sum_2 = 0;
  413.  
  414. /*for (int j = 0; j < Node_count; j++)
  415. {
  416. sum_2 += d_IP3 * (IP3(j) - IP3(in));
  417. }*/
  418.  
  419. for (int j = 0; j < C_A[in]; j++)
  420. {
  421. sum_2 += d_IP3 * (IP3((int)B_A[in][j]) - IP3(in));
  422. }
  423.  
  424. return (IP3_star - IP3(in)) / tau_IP3 + J_PLC(f, in) + sum_2 + J_Glu(f, in);
  425. }
  426.  
  427. case 2: // z
  428. {
  429. return a_2 * (d_2 * (IP3(in) + d_1) / (IP3(in) + d_3) * (1 - z(in)) - Ca(in) * z(in));
  430. }
  431.  
  432. case 3: // G_P
  433. {
  434. return -alpha_G_P * G_P(in) + beta_G_P * (1 / (1 + exp(-V_P(in) / 0.5)));
  435. }
  436.  
  437. case 4: // V
  438. {
  439. double I_syn = 0;
  440.  
  441. /*for (int j = 0; j < Node_count; j++)
  442. {
  443. //sum += A[in][j] * g_syn * (V(in) - V_old(j, tau[in][j]));
  444. //sum += A[in][j] * g_syn * (V(j) - V(in));
  445. //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old(j, tau[in][j]) / k_syn));
  446. //sum += 1 / (0.2 * Node_count_half) * A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V(j) / k_syn)); // i up, j down
  447. //sum += 1 / (0.2 * Node_count_half) * A[in][j] * g_syn * (V(j) - E_syn[j]) / (1 + exp(-V(in) / k_syn)); // j up, i down
  448. I_syn += A_N[in][j] * g_syn * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn)));
  449. //printf("in = %d\t Node_count = %d\t A[in][j] = %f\t V(in) = %f\t E_syn[in] = %f\t sum = %f\n", in, j, A[in][j], V(in), E_syn[in], sum);
  450. }*/
  451.  
  452. /*for (int j = 0; j < C[in]; j++)
  453. {
  454. sum += sigma[in][(int)B[in][j]] * (V((int)B[in][j]) - V(in));
  455. }*/
  456.  
  457. for (int j = 0; j < C_N[in]; j++)
  458. {
  459. if ((1 + g_astro * Ca(in)) > 0)
  460. {
  461. if (Ca(in) >= 0.3)
  462. {
  463. I_syn += g_syn * (1 + g_astro * Ca(in)) * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V((int)B_N[in][j]) / k_syn))); // версия с V (без V_old)
  464. I_syn += g_syn_P * (1 + g_astro * Ca(in)) * (E_syn_P[in] - V_P(in)) / (1 + exp(-(V_P(in) / k_syn_P))); // версия с V (без V_old)
  465. }
  466. else
  467. {
  468. I_syn += g_syn * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V((int)B_N[in][j]) / k_syn))); // версия с V (без V_old)
  469. I_syn += g_syn_P * (E_syn_P[in] - V_P(in)) / (1 + exp(-(V_P(in) / k_syn_P))); // версия с V (без V_old)
  470. }
  471. }
  472.  
  473. //I_syn += g_syn * (1 + g_astro * Ca(in)) * (E_syn[in] - V(in)) / (1 + exp(-(V_old((int)B_N[in][j], tau[in][(int)B_N[in][j]]) / k_syn)); // версия с V_old
  474. //I_syn += g_syn * (1 + g_astro * Ca(in)) * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия без с V (без V_old)
  475. //I_syn += g_syn * (E_syn[in] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия без с V (без V_old), упрощенная версия
  476. // sum_3 += g_syn * (1 + g_astro * Ca(in)) * (E_syn[i] - V(i)) / (1 + exp(-(V(j) / k_syn))); // образец из старой версии
  477. //I_syn += g_syn * (E_syn[in] - V(in)) / (1 + exp(-(V_old((int)B_N[in][j], tau[in][(int)B_N[in][j]])) / k_syn));*/
  478. //I_syn += g_syn * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V(j) / k_syn))); // версия с V (без V_old)
  479. //I_syn += g_syn * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V((int)B_N[in][j]) / k_syn))); // версия с V (без V_old), упрощенная версия !!!
  480. //I_syn += g_syn * (E_syn[(int)B_N[in][j]] - V(in)) / (1 + exp(-(V_old((int)B_N[in][j], tau[in][(int)B_N[in][j]]) / k_syn))); // версия с V (без V_old), упрощенная версия !!!
  481. //sum += g_syn * (V((int)B[in][j]) - V(in)); // устаревшая часть, нужна для проверки разностной схемы
  482. //sum += A[in][j] * g_syn * (V(in) - E_syn[in]) / (1 + exp(-V_old((int)B[in][j]) / k_syn)); // устаревшая часть
  483. //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));
  484. //sum += A[in][(int)B[in][j]] * g_syn * (V((int)B[in][j]) - E_syn[(int)B[in][j]]) / (1 + exp(-V_old(in, tau[in][(int)B[in][j]]) / k_syn));
  485. //sum += 1 / (0.2 * Node_count_half) * A[in][(int)B[in][j]] * g_syn * (V((int)B[in][j]) - E_syn[(int)B[in][j]]) / (1 + exp(-V(in) / k_syn)); // j up, i down
  486. //sum += 1 / (0.2 * Node_count) * /*(int)A_N[in][(int)B_N[in][j]] * */ g_syn * (1 + g_astro * Ca(in)) * (V(in) - E_syn[in]) / (1 + exp(-V((int)B_N[in][j]) / k_syn)); // i up, j down
  487. // i up, j down
  488. /*printf("i = %d\t j = %d\t A[i, j] = %d\n", in, (int)B_N[in][j], (int)A_N[in][(int)B_N[in][j]]);*/
  489. /*printf("i = %d\t V_old = %f\t exp = %f\n", in, Vold, ee);*/
  490. }
  491. /*printf("i = %d\t sum = %f\n", in, sum);*/
  492.  
  493. if (enable_I_syn_out)
  494. fprintf(fp_I_syn, i == Equations_count - 1 ? "%f" : "%f\t", I_syn);
  495.  
  496. 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_syn) / C_m); // V
  497. }
  498.  
  499. case 5: // m
  500. {
  501. return 1000 * (alpha_m(f, in) * (1 - m(in)) - beta_m(f, in) * m(in)); // m
  502. }
  503.  
  504. case 6: // n
  505. {
  506. return 1000 * (alpha_n(f, in) * (1 - n(in)) - beta_n(f, in) * n(in)); // n
  507. }
  508.  
  509. case 7: // h
  510. {
  511. return 1000 * (alpha_h(f, in) * (1 - h(in)) - beta_h(f, in) * h(in)); // h
  512. }
  513.  
  514. case 8: // V_P
  515. {
  516. return 1000 * ((g_Na_P * pow(m_P(in), 3) * h_P(in) * (E_Na_P - V_P(in)) + g_K_P * n_P(in) * (E_K_P - V_P(in)) + g_L_P * (E_L_P - V_P(in)) + I_app_P[in] + I_stim(in, t)) / C_m_P); // V_P
  517. }
  518.  
  519. case 9: // m_P
  520. {
  521. return 1000 * (alpha_m_P(f, in) * (1 - m_P(in)) - beta_m_P(f, in) * m_P(in)); // m_P
  522. }
  523.  
  524. case 10: // n_P
  525. {
  526. return 1000 * (alpha_n_P(f, in) * (1 - n_P(in)) - beta_n_P(f, in) * n_P(in)); // n_P
  527. }
  528.  
  529. case 11: // h_P
  530. {
  531. return 1000 * (alpha_h_P(f, in) * (1 - h_P(in)) - beta_h_P(f, in) * h_P(in)); // h_P
  532. }
  533. }
  534.  
  535. return 0;
  536. }
  537.  
  538. void RungeKutta(double t, double dt, double* f, double* f_next)
  539. {
  540. // k1
  541. #pragma omp parallel for
  542. for (int i = 0; i < Equations_count; i++)
  543. {
  544. //if (!thread_count_printed)
  545. //{
  546. // thread_count_printed = true;
  547. // printf("Threads = %d\n", omp_get_num_threads());
  548. //}
  549.  
  550. k[i][0] = UllahJung_HodgkinHuxley(i, f, t) * dt;
  551. phi_k1[i] = f[i] + k[i][0] / 2;
  552. k[i][1] = UllahJung_HodgkinHuxley(i, phi_k1, t) * dt;
  553. phi_k2[i] = f[i] + k[i][1] / 2;
  554. k[i][2] = UllahJung_HodgkinHuxley(i, phi_k2, t) * dt;
  555. phi_k3[i] = f[i] + k[i][2] / 2;
  556. k[i][3] = UllahJung_HodgkinHuxley(i, phi_k3, t) * dt;
  557. f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
  558. }
  559.  
  560. //for (int i = 0; i < Equations_count; i++)
  561. // phi_k1[i] = f[i] + k[i][0] / 2;
  562.  
  563. // k2
  564. //for (int i = 0; i < Equations_count; i++)
  565. // k[i][1] = UllahJung_HodgkinHuxley(i, phi_k1, t) * dt;
  566.  
  567.  
  568. //for (int i = 0; i < Equations_count; i++)
  569. // phi_k2[i] = f[i] + k[i][1] / 2;
  570.  
  571. // k3
  572. //for (int i = 0; i < Equations_count; i++)
  573. // k[i][2] = UllahJung_HodgkinHuxley(i, phi_k2, t) * dt;
  574.  
  575.  
  576. //for (int i = 0; i < Equations_count; i++)
  577. // phi_k3[i] = f[i] + k[i][2] / 2;
  578.  
  579. //enable_I_syn_out = true;
  580.  
  581. // k4
  582. //for (int i = 0; i < Equations_count; i++)
  583. // k[i][3] = UllahJung_HodgkinHuxley(i, phi_k3, t) * dt;
  584.  
  585. //enable_I_syn_out = false;
  586.  
  587. //for (int i = 0; i < Equations_count; i++)
  588. // f_next[i] = f[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
  589. }
  590.  
  591. void CopyArray(double* source, double* target, int N)
  592. {
  593. for (int i = 0; i < N; i++)
  594. target[i] = source[i];
  595. }
  596.  
  597. bool Approximately(double a, double b)
  598. {
  599. if (a < 0)
  600. a = -a;
  601.  
  602. if (b < 0)
  603. b = -b;
  604.  
  605. return a - b <= 0.000001;
  606. }
  607.  
  608. //bool CheckSameLine(int i, int j)
  609. //{
  610. // return i / Node_wire_width == j / Node_wire_width;
  611. //}
  612. //
  613. //bool IsWireNeighbors(int i, int j)
  614. //{
  615. // if (CheckSameLine(i, j) && (i == j - 1 || i == j + 1))
  616. // return true;
  617. //
  618. // if (i == j - Node_wire_width || i == j + Node_wire_width)
  619. // return true;
  620. //
  621. // return false;
  622. //}
  623.  
  624. // http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/
  625. double nextTime(double rateParameter)
  626. {
  627. return -log(1.0 - (double)rand() / (RAND_MAX)) / rateParameter;
  628. }
  629.  
  630. void GenerateRandomMeander(int i, double min_start_time)
  631. {
  632. double offset = nextTime(Poisson_Freq);
  633.  
  634. if (offset < 0)
  635. {
  636. int a = 0;
  637. }
  638.  
  639. Meander_start_from_zero[i] = min_start_time + offset;
  640. Meander_width[i] = Duration;
  641. //Meander_height[i] = RandomD(-Max_magnitude, Max_magnitude);
  642. //Meander_height[i] = RandomD(0, Max_magnitude);
  643. Meander_height[i] = Max_magnitude;
  644. }
  645.  
  646. void FillAMatrixZero()
  647. {
  648. for (int i = 0; i < Node_count; i++)
  649. {
  650. for (int j = 0; j < Node_count; j++)
  651. {
  652. A_A[i][j] = 0;
  653. A_N[i][j] = 0;
  654. A_N_P[i][j] = 0;
  655. }
  656. }
  657. }
  658.  
  659. void FillAstrociteMatrix()
  660. {
  661. for (int i = 0; i < Node_count; i++)
  662. {
  663. for (int j = 0; j < Node_count; j++)
  664. {
  665. if (i == j)
  666. {
  667. A_A[i][j] = 0;
  668. continue;
  669. }
  670.  
  671. if (i > j)
  672. {
  673. A_A[i][j] = A_A[j][i];
  674. continue;
  675. }
  676.  
  677. if (i == 0 && j == Node_count - 1)
  678. {
  679. A_A[i][j] = 1;
  680. continue;
  681. }
  682.  
  683. if (i == Node_count - 1 && j == 0)
  684. {
  685. A_A[i][j] = 1;
  686. continue;
  687. }
  688.  
  689. if (i == j - 1 || i == j + 1)
  690. {
  691. A_A[i][j] = 1;
  692. continue;
  693. }
  694. }
  695. }
  696. //A_A[0][1] = 0; // only for debug. diffusion Ca test
  697. //A_A[1][0] = 0; // only for debug. diffusion Ca test
  698. //A_A[1][3] = 0;
  699. //A_A[3][1] = 0;
  700. //A_A[0][2] = 0;
  701. //A_A[2][0] = 0;
  702. }
  703.  
  704. void FillBCMatrix_A()
  705. {
  706. for (int i = 0; i < Node_count; i++)
  707. {
  708. int bIndex = 0;
  709. C_A[i] = 0;
  710. for (int j = 0; j < Node_count; j++)
  711. {
  712. if (A_A[i][j] == 1)
  713. {
  714. B_A[i][bIndex] = j;
  715. bIndex++;
  716. C_A[i]++;
  717. }
  718. }
  719. }
  720. }
  721.  
  722. void FillBCMatrix_N()
  723. {
  724. for (int i = 0; i < Node_count; i++)
  725. {
  726. int bIndex = 0;
  727. C_N[i] = 0;
  728. for (int j = 0; j < Node_count; j++)
  729. {
  730. if (A_N[i][j] == 1)
  731. {
  732. B_N[i][bIndex] = j;
  733. bIndex++;
  734. C_N[i]++;
  735. }
  736. }
  737. }
  738. }
  739.  
  740. void FillBCMatrix_N_P()
  741. {
  742. for (int i = 0; i < Node_count; i++)
  743. {
  744. int bIndex = 0;
  745. C_N_P[i] = 0;
  746. for (int j = 0; j < Node_count; j++)
  747. {
  748. if (A_N_P[i][j] == 1)
  749. {
  750. B_N_P[i][bIndex] = j;
  751. bIndex++;
  752. C_N_P[i]++;
  753. }
  754. }
  755. }
  756. }
  757.  
  758. bool IsWireNeighbors(int i, int j, int deep)
  759. {
  760. int j_border_left = j - deep < 0 ? j + Node_count : j;
  761. int j_border_right = j + deep >= Node_count ? j - Node_count : j;
  762.  
  763. if (i == j_border_left - deep || i == j_border_right + deep)
  764. {
  765. return true;
  766. }
  767.  
  768. return false;
  769. }
  770.  
  771. void FillNeuronMatrix()
  772. {
  773. for (int i = 0; i < Node_count; i++)
  774. {
  775. for (int j = 0; j < Node_count; j++)
  776. {
  777. if (i == j)
  778. {
  779. A_N[i][j] = 0;
  780. continue;
  781. }
  782.  
  783. if (i > j)
  784. {
  785. A_N[i][j] = A_N[j][i];
  786. continue;
  787. }
  788.  
  789. for (int deep = 1; deep <= MaxDeep; deep++)
  790. {
  791. if (IsWireNeighbors(i, j, deep))
  792. A_N[i][j] = 1;
  793. }
  794. }
  795. }
  796. }
  797.  
  798. void RandomizeNeuronMatrix()
  799. {
  800. srand(time(NULL));
  801.  
  802. for (int i = 0; i < Node_count; i++)
  803. {
  804. //if (i == Node_count / 2)
  805. // srand(time(NULL));
  806.  
  807. for (int link = 0; link < MaxDeep * 2; link++)
  808. {
  809. double x = RandomD(0, 1);
  810.  
  811. if (x > p_rewir)
  812. continue;
  813.  
  814. int rndJ;
  815.  
  816. do
  817. {
  818. rndJ = RandomI(0, Node_count);
  819. } while (i == rndJ || A_N[i][rndJ] == 1);
  820.  
  821. int rndJ_last;
  822.  
  823. do
  824. {
  825. rndJ_last = RandomI(i - MaxDeep - 1, i + MaxDeep + 1);
  826.  
  827. if (rndJ_last < 0)
  828. rndJ_last += Node_count;
  829. else if (rndJ_last >= Node_count)
  830. rndJ_last -= Node_count;
  831.  
  832. } while (i == rndJ_last || A_N[i][rndJ_last] == 0);
  833.  
  834. A_N[i][rndJ_last] = 0;
  835.  
  836. A_N[i][rndJ] = 1;
  837. }
  838. }
  839. }
  840.  
  841. void FillNeuronPoissonMatrix()
  842. {
  843. for (int i = 0; i < Node_count; i++)
  844. {
  845. for (int j = 0; j < Node_count; j++)
  846. {
  847. A_N_P[i][j] = 0;
  848. }
  849. }
  850. }
  851.  
  852. void FillVOldFromCurrent()
  853. {
  854. for (int i = 0; i < Node_count; i++)
  855. for (int j = 0; j < Max_delay; j++)
  856. V_old_array[i][j] = V(i);
  857. }
  858.  
  859. void UpdateVOld()
  860. {
  861. for (int i = 0; i < Node_count; i++)
  862. {
  863. for (int j = 1; j < Max_delay; j++)
  864. V_old_array[i][j - 1] = V_old_array[i][j];
  865.  
  866. V_old_array[i][Max_delay - 1] = V(i);
  867. }
  868. }
  869.  
  870. //void FillFullTauMatrix()
  871. //{
  872. // for (int i = 0; i < Node_count; i++)
  873. // {
  874. // for (int j = 0; j < Node_count; j++)
  875. // {
  876. // if (i < Node_count || j < Node_count)
  877. // {
  878. // tau[i][j] = 0;
  879. // continue;
  880. // }
  881. //
  882. // int i_neuron = i - Node_count;
  883. // int j_neuron = j - Node_count;
  884. //
  885. // int i_wire_x = i_neuron / Node_wire_width;
  886. // int i_wire_y = i_neuron % Node_wire_width;
  887. //
  888. // int j_wire_x = j_neuron / Node_wire_width;
  889. // int j_wire_y = j_neuron % Node_wire_width;
  890. //
  891. // double distance_max = sqrt(2.) * (Node_wire_width - 1);
  892. // 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));
  893. //
  894. // tau[i][j] = (tau_min + distance / (distance_max) * (tau_max - tau_min)) * ms_to_step;
  895. // }
  896. // }
  897. //}
  898. //
  899. //void FillTauMatrix()
  900. //{
  901. // for (int i = 0; i < Node_count; i++)
  902. // {
  903. // for (int j = 0; j < Node_count; j++)
  904. // {
  905. // if (i == j || A_N[i][j] == 0)
  906. // {
  907. // tau[i][j] = 0;
  908. // continue;
  909. // }
  910. //
  911. // int i_neuron = i;
  912. // int j_neuron = j;
  913. //
  914. // int i_wire_x = i_neuron / Node_wire_width;
  915. // int i_wire_y = i_neuron % Node_wire_width;
  916. //
  917. // int j_wire_x = j_neuron / Node_wire_width;
  918. // int j_wire_y = j_neuron % Node_wire_width;
  919. //
  920. // double distance_max = sqrt(2.) * (Node_wire_width - 1);
  921. // 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));
  922. //
  923. // double t = (distance - 1) / (distance_max - 1);
  924. // tau[i][j] = (tau_min + t * (tau_max - tau_min)) * ms_to_step;
  925. // }
  926. // }
  927. //}
  928.  
  929. int main(int argc, char* argv[])
  930. {
  931. // run like: UJ_HH_Ring_acc.out 250 0.2 0.4 0.05 3.0 1.05 1.50 // (1) p_rewir (2) g_syn (3) g_astro (4) Max_magnitude (5) g_syn_P
  932. /*sscanf(argv[1], "%d", &Node_count);
  933. FILE* fp_Node_count;
  934. fp_Node_count = fopen("Node_count.txt", "w");
  935. fprintf(fp_Node_count, "%d\t", Node_count);
  936. fclose(fp_Node_count);
  937. printf("Node_count = %d\n", Node_count);*/
  938.  
  939. sscanf(argv[1], "%lf", &p_rewir);
  940. FILE* fp_p_rewir;
  941. fp_p_rewir = fopen("p_rewir.txt", "w");
  942. fprintf(fp_p_rewir, "%f\t", p_rewir);
  943. fclose(fp_p_rewir);
  944. printf("p_rewir = %f\n", p_rewir);
  945.  
  946. /*sscanf(argv[3], "%lf", &p_inhib);
  947. FILE* fp_p_inhib;
  948. fp_p_inhib = fopen("p_inhib.txt", "w");
  949. fprintf(fp_p_inhib, "%f\t", p_inhib);
  950. fclose(fp_p_inhib);
  951. printf("p_inhib = %f\n", p_inhib);*/
  952.  
  953. sscanf(argv[2], "%lf", &g_syn);
  954. FILE* fp_g_syn;
  955. fp_g_syn = fopen("g_syn.txt", "w");
  956. fprintf(fp_g_syn, "%f\t", g_syn);
  957. fclose(fp_g_syn);
  958. printf("g_syn = %f\n", g_syn);
  959.  
  960. sscanf(argv[3], "%lf", &g_astro);
  961. FILE* fp_g_astro;
  962. fp_g_astro = fopen("g_astro.txt", "w");
  963. fprintf(fp_g_astro, "%f\t", g_astro);
  964. fclose(fp_g_astro);
  965. printf("g_astro = %f\n", g_astro);
  966.  
  967. sscanf(argv[4], "%lf", &Max_magnitude);
  968. FILE* fp_Max_magnitude;
  969. fp_Max_magnitude = fopen("Max_magnitude.txt", "w");
  970. fprintf(fp_Max_magnitude, "%f\t", Max_magnitude);
  971. fclose(fp_Max_magnitude);
  972. printf("I_stim magnitude = %f\n", Max_magnitude);
  973.  
  974. sscanf(argv[5], "%lf", &g_syn_P);
  975. FILE* fp_g_syn_P;
  976. fp_g_syn_P = fopen("g_syn_P.txt", "w");
  977. fprintf(fp_g_syn_P, "%f\t", g_syn_P);
  978. fclose(fp_g_syn_P);
  979. printf("g_syn_P = %f\n", g_syn_P);
  980.  
  981. //sscanf(argv[6], "%lf", &I_app_min);
  982. //sscanf(argv[7], "%lf", &I_app_max);
  983.  
  984. f = new double[Equations_count];
  985. f_diff = new double[Equations_count];
  986. v_4 = new double[Node_count];
  987. E_syn = new double[Node_count];
  988. I_app = new double[Node_count];
  989. E_syn_P = new double[Node_count];
  990. I_app_P = new double[Node_count];
  991. V_spikes = new list<double>[Node_count];
  992. V_spikes_Freq = new list<double>[Node_count];
  993.  
  994. Meander_start_from_zero = new double[Node_count];
  995. Meander_width = new double[Node_count];
  996. Meander_height = new double[Node_count];
  997. Meander_interval = new double[Node_count];
  998. last_meander_end = new double[Node_count];
  999.  
  1000. tau = new double* [Node_count];
  1001. for (int i = 0; i < Node_count; i++)
  1002. tau[i] = new double[Node_count];
  1003.  
  1004. V_old_array = new double* [Node_count];
  1005. for (int i = 0; i < Node_count; i++)
  1006. V_old_array[i] = new double[Max_delay];
  1007.  
  1008. FILE* fp0;
  1009. FILE* fp_I_stim;
  1010. FILE* fp_Ca;
  1011. FILE* fp_IP3;
  1012. //FILE *fp_z;
  1013. //FILE* fp_G_P;
  1014. FILE* fp_V;
  1015. FILE* fp_V_P;
  1016. //FILE *fp_m;
  1017. //FILE *fp_n;
  1018. //FILE *fp_h;
  1019. FILE* fp_V_spikes;
  1020. FILE* fp_Esyn;
  1021. FILE* fp_Esyn_P;
  1022.  
  1023. //FILE* fp_res;
  1024. srand(time(NULL));
  1025.  
  1026. //for (int i = 0; i < Node_count; i++)
  1027. // V_old_length[i] = 0;
  1028.  
  1029. A_A = new double* [Node_count];
  1030. for (int i = 0; i < Node_count; i++)
  1031. A_A[i] = new double[Node_count];
  1032.  
  1033. B_A = new double* [Node_count];
  1034. for (int i = 0; i < Node_count; i++)
  1035. B_A[i] = new double[Node_count];
  1036.  
  1037. C_A = new double[Node_count];
  1038.  
  1039. A_N = new double* [Node_count];
  1040. for (int i = 0; i < Node_count; i++)
  1041. A_N[i] = new double[Node_count];
  1042.  
  1043. B_N = new double* [Node_count];
  1044. for (int i = 0; i < Node_count; i++)
  1045. B_N[i] = new double[Node_count];
  1046.  
  1047. C_N = new double[Node_count];
  1048.  
  1049. A_N_P = new double* [Node_count];
  1050. for (int i = 0; i < Node_count; i++)
  1051. A_N_P[i] = new double[Node_count];
  1052.  
  1053. B_N_P = new double* [Node_count];
  1054. for (int i = 0; i < Node_count; i++)
  1055. B_N_P[i] = new double[Node_count];
  1056.  
  1057. C_N_P = new double[Node_count];
  1058.  
  1059. FillAMatrixZero();
  1060. FillAstrociteMatrix();
  1061. FillNeuronMatrix();
  1062. RandomizeNeuronMatrix();
  1063. FillNeuronPoissonMatrix();
  1064. FillBCMatrix_A();
  1065. FillBCMatrix_N();
  1066. FillBCMatrix_N_P();
  1067. //FillTauMatrix();
  1068.  
  1069. fp0 = fopen("A_A.txt", "w+");
  1070. for (int i = 0; i < Node_count; i++)
  1071. {
  1072. for (int j = 0; j < Node_count; j++)
  1073. {
  1074. fprintf(fp0, "%d\t", (int)A_A[i][j]);
  1075. }
  1076. fprintf(fp0, "\n");
  1077. }
  1078. fclose(fp0);
  1079.  
  1080. fp0 = fopen("A_N.txt", "w+");
  1081. for (int i = 0; i < Node_count; i++)
  1082. {
  1083. for (int j = 0; j < Node_count; j++)
  1084. {
  1085. fprintf(fp0, "%d\t", (int)A_N[i][j]);
  1086. }
  1087. fprintf(fp0, "\n");
  1088. }
  1089. fclose(fp0);
  1090.  
  1091. fp0 = fopen("tau.txt", "w+");
  1092. for (int i = 0; i < Node_count; i++)
  1093. {
  1094. for (int j = 0; j < Node_count; j++)
  1095. {
  1096. fprintf(fp0, "%f\t", tau[i][j] / ms_to_step);
  1097. }
  1098. fprintf(fp0, "\n");
  1099. }
  1100. fclose(fp0);
  1101.  
  1102. // Write to file number of links for each neuron
  1103. fp0 = fopen("links.txt", "w+");
  1104. for (int i = 0; i < Node_count; i++)
  1105. {
  1106. int links_count = 0;
  1107. for (int j = 0; j < Node_count; j++)
  1108. {
  1109. if (A_N[i][j] == 1)
  1110. {
  1111. links_count++;
  1112. }
  1113. }
  1114. fprintf(fp0, "%d\n", (int)links_count);
  1115. }
  1116. fclose(fp0);
  1117.  
  1118. //setlocale(LC_NUMERIC, "French_Canada.1252");
  1119. fp0 = fopen("test_Poisson.txt", "w+");
  1120. for (int i = 0; i < 1000; i++)
  1121. fprintf(fp0, "%f\n", nextTime(Poisson_Freq));
  1122. fclose(fp0);
  1123.  
  1124. fp0 = fopen("B_A.txt", "w+");
  1125. for (int i = 0; i < Node_count; i++)
  1126. {
  1127. for (int j = 0; j < C_A[i]; j++)
  1128. {
  1129. fprintf(fp0, "%d\t", (int)B_A[i][j]);
  1130. }
  1131. fprintf(fp0, "\n");
  1132. }
  1133. fclose(fp0);
  1134.  
  1135. fp0 = fopen("B_N.txt", "w+");
  1136. for (int i = 0; i < Node_count; i++)
  1137. {
  1138. for (int j = 0; j < C_N[i]; j++)
  1139. {
  1140. fprintf(fp0, "%d\t", (int)B_N[i][j]);
  1141. }
  1142. fprintf(fp0, "\n");
  1143. }
  1144. fclose(fp0);
  1145.  
  1146. fp0 = fopen("C_A.txt", "w+");
  1147. for (int i = 0; i < Node_count; i++)
  1148. {
  1149. fprintf(fp0, "%d\n", (int)C_A[i]);
  1150. }
  1151. fclose(fp0);
  1152.  
  1153. fp0 = fopen("C_N.txt", "w+");
  1154. for (int i = 0; i < Node_count; i++)
  1155. {
  1156. fprintf(fp0, "%d\n", (int)C_N[i]);
  1157. }
  1158. fclose(fp0);
  1159.  
  1160. /*for (int i = 0; i < 6; i++)
  1161. {
  1162. v_4[i] = 0.6;
  1163. }*/
  1164. for (int i = 0; i < Node_count; i++)
  1165. {
  1166. v_4[i] = 0.4; // 0.4
  1167. }
  1168.  
  1169. // Initial values
  1170. /*for (int i = 0; i < Equations_count; i++)
  1171. {
  1172. f[i] = 0;
  1173. }*/
  1174.  
  1175. //I_app_min = 1.1;
  1176. //I_app_max = 1.5;
  1177.  
  1178. FILE* fp_I_app;
  1179. fp_I_app = fopen("I_app.txt", "w");
  1180. for (int i = 0; i < Node_count; i++)
  1181. {
  1182. I_app[i] = 0.8; //RandomD(I_app_min, I_app_max);
  1183. fprintf(fp_I_app, "%f\n", I_app[i]);
  1184. }
  1185. fclose(fp_I_app);
  1186.  
  1187. FILE* fp_I_app_P;
  1188. fp_I_app_P = fopen("I_app_P.txt", "w");
  1189. for (int i = 0; i < Node_count; i++)
  1190. {
  1191. I_app_P[i] = 0.7;
  1192. fprintf(fp_I_app_P, "%f\n", I_app_P[i]);
  1193. }
  1194. fclose(fp_I_app_P);
  1195.  
  1196. for (int i = 0; i < Node_count; i++) // init array for all nodes
  1197. {
  1198. SetG_P(i, 0); // G_P
  1199. }
  1200.  
  1201. double percent_stable_state = 0.50; // 0.40
  1202. double eps_persent = 0.05; //0.05
  1203.  
  1204. double Ca0 = 0.07;
  1205. double IP30 = 0.16;
  1206. double z0 = 0.67;
  1207.  
  1208. for (int i = 0; i < Node_count; i++)
  1209. {
  1210. /*SetCa(i, Ca0 + RandomD(-Ca0 * eps_persent, Ca0 * eps_persent)); // Ca
  1211. SetIP3(i, IP30 + RandomD(-IP30 * eps_persent, IP30 * eps_persent)); // IP3
  1212. Setz(i, z0 + RandomD(-z0 * eps_persent, z0 * eps_persent)); // z */
  1213. SetCa(i, Ca0); // Ca
  1214. SetIP3(i, IP30); // IP3
  1215. Setz(i, z0); // z
  1216. }
  1217.  
  1218. double V0 = -58.7085;
  1219. double m0 = 0.0953;
  1220. double n0 = 0.000913;
  1221. double h0 = 0.3662;
  1222.  
  1223. double V1 = 14.8409;
  1224. double m1 = 0.9174;
  1225. double n1 = 0.0140;
  1226. double h1 = 0.0539;
  1227.  
  1228. /*for (int i = 0; i < Node_count; i++) // init only for neurons
  1229. {
  1230. double random = RandomD(0, 1);
  1231.  
  1232. SetV(i, random < percent_stable_state ? V0 + RandomD(-V0 * eps_persent, V0 * eps_persent) : V1 + RandomD(-V1 * eps_persent, V1 * eps_persent)); // V
  1233. Setm(i, random < percent_stable_state ? m0 + RandomD(-m0 * eps_persent, m0 * eps_persent) : m1 + RandomD(-m1 * eps_persent, m1 * eps_persent)); // m
  1234. Setn(i, random < percent_stable_state ? n0 + RandomD(-n0 * eps_persent, n0 * eps_persent) : n1 + RandomD(-n1 * eps_persent, n1 * eps_persent)); // n
  1235. Seth(i, random < percent_stable_state ? h0 + RandomD(-h0 * eps_persent, h0 * eps_persent) : h1 + RandomD(-h1 * eps_persent, h1 * eps_persent)); // h
  1236. }*/
  1237.  
  1238. for (int i = 0; i < Node_count; i++) // init only for neurons
  1239. {
  1240. double random = RandomD(0, 1);
  1241.  
  1242. /*SetV(i, random < percent_stable_state ? V0 : V1); // V
  1243. Setm(i, random < percent_stable_state ? m0 : m1); // m
  1244. Setn(i, random < percent_stable_state ? n0 : n1); // n
  1245. Seth(i, random < percent_stable_state ? h0 : h1); // h*/
  1246.  
  1247. SetV(i, V0); // V
  1248. Setm(i, m0); // m
  1249. Setn(i, n0); // n
  1250. Seth(i, h0); // h
  1251.  
  1252. /*SetV_P(i, random < percent_stable_state ? V0 : V1); // V
  1253. Setm_P(i, random < percent_stable_state ? m0 : m1); // m
  1254. Setn_P(i, random < percent_stable_state ? n0 : n1); // n
  1255. Seth_P(i, random < percent_stable_state ? h0 : h1); // h*/
  1256.  
  1257. SetV_P(i, V0); // V
  1258. Setm_P(i, m0); // m
  1259. Setn_P(i, n0); // n
  1260. Seth_P(i, h0); // h
  1261. }
  1262.  
  1263. /*for (int i = 0; i < Node_count; i++) // init only for neurons
  1264. {*/
  1265. /*SetV(i, RandomD(-80, 20)); // V
  1266. Setm(i, RandomD(0, 1)); // m
  1267. Setn(i, RandomD(0, 1)); // n
  1268. Seth(i, RandomD(0, 1)); // h*/
  1269. /*SetV(i, V1); // V
  1270. Setm(i, m1); // m
  1271. Setn(i, n1); // n
  1272. Seth(i, h1); // h
  1273. }*/
  1274.  
  1275. double E_syn0 = 0; // Excitatory neuron
  1276. double E_syn1 = -90; // Inhibitory neuron
  1277.  
  1278. fp_Esyn = fopen("results_E_syn.txt", "w+");
  1279. for (int i = 0; i < Node_count; i++)
  1280. {
  1281. /*E_syn[i] = E_syn0;
  1282.  
  1283. double x = RandomD(0, 1);
  1284.  
  1285. if (x > p_inhib)
  1286. {
  1287. fprintf(fp_Esyn, "%f\n", E_syn[i]);
  1288. continue;
  1289. }*/
  1290.  
  1291. E_syn[i] = E_syn1;
  1292.  
  1293. fprintf(fp_Esyn, "%f\n", E_syn[i]);
  1294.  
  1295. E_syn_P[i] = E_syn0;
  1296. }
  1297. fclose(fp_Esyn);
  1298.  
  1299. for (int i = 0; i < Node_count; i++)
  1300. {
  1301. GenerateRandomMeander(i, 0);
  1302. last_meander_end[i] = Meander_start_from_zero[i] + Duration;
  1303. }
  1304.  
  1305. const double t_start = 0;
  1306. const double t_max = 40; // 100 msec = 0.1 sec // 240 // 270
  1307. const double dt = 0.00005; // 0.01 msec = 0.00001 sec; 0.1 msec = 0.0001 sec; 1 msec = 0.001 sec // 0.000025
  1308.  
  1309. double t = t_start;
  1310.  
  1311. //fp_res = fopen("matlab_res.txt", "a+");
  1312. //fprintf(fp_res, "%f\t%f\t", Max_magintude, g_syn);
  1313.  
  1314. //fp0 = fopen("results.txt", "w+");
  1315. //setlocale(LC_NUMERIC, "French_Canada.1252");
  1316.  
  1317. //double start_rk4, end_rk4;
  1318. clock_t start_rk4, end_rk4;
  1319. //start_rk4 = omp_get_wtime();
  1320. start_rk4 = clock();
  1321. int lastPercent = -1;
  1322.  
  1323. //FillVOldFromCurrent();
  1324.  
  1325. k = new double* [Equations_count];
  1326. for (int i = 0; i < Equations_count; i++)
  1327. k[i] = new double[4];
  1328.  
  1329. phi_k1 = new double[Equations_count];
  1330. phi_k2 = new double[Equations_count];
  1331. phi_k3 = new double[Equations_count];
  1332.  
  1333. fp_I_stim = fopen("results_I_stim.txt", "w+");
  1334. //fp_I_syn = fopen("results_I_syn.txt", "w+");
  1335. fp_Ca = fopen("results_Ca.txt", "w+");
  1336. //fp_IP3 = fopen("results_IP3.txt", "w+");
  1337. //fp_z = fopen("results_z.txt", "w+");
  1338. //fp_G_P = fopen("results_G_P.txt", "w+");
  1339. fp_V = fopen("results_V.txt", "w+");
  1340. fp_V_P = fopen("results_V_P.txt", "w+");
  1341. //fp_m = fopen("results_m.txt", "w+");
  1342. //fp_n = fopen("results_n.txt", "w+");
  1343. //fp_h = fopen("results_h.txt", "w+");
  1344. fp_V_spikes = fopen("results_V_spikes.txt", "w+");
  1345. //
  1346.  
  1347. double* f_next = new double[Equations_count];
  1348.  
  1349. while (t < t_max || Approximately(t, t_max))
  1350. {
  1351. fprintf(fp_I_stim, "%f\t", t);
  1352. fprintf(fp_Ca, "%f\t", t);
  1353. //fprintf(fp_IP3, "%f\t", t);
  1354. //fprintf(fp_z, "%f\t", t);
  1355. //fprintf(fp_G, "%f\t", t);
  1356. fprintf(fp_V, "%f\t", t);
  1357. fprintf(fp_V_P, "%f\t", t);
  1358. //fprintf(fp_m, "%f\t", t);
  1359. //fprintf(fp_n, "%f\t", t);
  1360. //fprintf(fp_h, "%f\t", t);
  1361. fprintf(fp_V_spikes, "%f\t", t);
  1362.  
  1363. for (int i = 0; i < Node_count; i++)
  1364. {
  1365. if (t > last_meander_end[i])
  1366. {
  1367. GenerateRandomMeander(i, t);
  1368. last_meander_end[i] = Meander_start_from_zero[i] + Duration;
  1369. }
  1370.  
  1371. fprintf(fp_I_stim, "%f\t", I_stim(i, t));
  1372. }
  1373. fprintf(fp_I_stim, "\n");
  1374.  
  1375. for (int i = 0; i < Equations_count; i += Equations_per_node)
  1376. fprintf(fp_Ca, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // Ca
  1377.  
  1378. //for (int i = 1; i < Equations_count; i += Equations_per_node)
  1379. // fprintf(fp_IP3, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // IP3
  1380.  
  1381. //for (int i = 2; i < Equations_count; i += Equations_per_node)
  1382. // fprintf(fp_z, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // z
  1383.  
  1384. //for (int i = 3; i < Equations_count; i += Equations_per_node)
  1385. // fprintf(fp_G_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // G
  1386.  
  1387. for (int i = 4; i < Equations_count; i += Equations_per_node)
  1388. fprintf(fp_V, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V
  1389.  
  1390. //for (int i = 5; i < Equations_count; i += Equations_per_node)
  1391. // fprintf(fp_m, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // m
  1392.  
  1393. //for (int i = 6; i < Equations_count; i += Equations_per_node)
  1394. // fprintf(fp_n, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // n
  1395.  
  1396. //for (int i = 7; i < Equations_count; i += Equations_per_node)
  1397. // fprintf(fp_h, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // h
  1398.  
  1399. for (int i = 8; i < Equations_count; i += Equations_per_node)
  1400. fprintf(fp_V_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // V_P
  1401.  
  1402. //for (int i = 9; i < Equations_count; i += Equations_per_node)
  1403. // fprintf(fp_m_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // m_P
  1404.  
  1405. //for (int i = 10; i < Equations_count; i += Equations_per_node)
  1406. // fprintf(fp_n_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // n_P
  1407.  
  1408. //for (int i = 11; i < Equations_count; i += Equations_per_node)
  1409. // fprintf(fp_h_P, i == Equations_count - 1 ? "%f" : "%f\t", f[i]); // h_P
  1410.  
  1411. fprintf(fp_Ca, "\n");
  1412. //fprintf(fp_IP3, "\n");
  1413. //fprintf(fp_z, "\n");
  1414. //fprintf(fp_G, "\n");
  1415. fprintf(fp_V, "\n");
  1416. //fprintf(fp_m, "\n");
  1417. //fprintf(fp_n, "\n");
  1418. //fprintf(fp_h, "\n");
  1419.  
  1420. //fprintf(fp_G_P, "\n");
  1421. fprintf(fp_V_P, "\n");
  1422. //fprintf(fp_m_P, "\n");
  1423. //fprintf(fp_n_P, "\n");
  1424. //fprintf(fp_h_P, "\n");
  1425.  
  1426. RungeKutta(t, dt, f, f_next);
  1427.  
  1428. #pragma omp parallel for
  1429. for (int i = 0; i < Node_count; i++)
  1430. {
  1431. int index = Equations_per_node * i + 4;
  1432. double diff = f_next[index] - f[index];
  1433.  
  1434. fprintf(fp_V_spikes, i == Equations_count - 1 ? "%d" : "%d\t", diff < 0 && f_diff[index] > 0 && f[index] > -10 && (V_spikes[i].size() == 0 || t - V_spikes[i].back() > 0.001) ? 1 : 0);
  1435.  
  1436. if (diff < 0 && f_diff[index] > 0 && f[index] > -10 && (V_spikes[i].size() == 0 || t - V_spikes[i].back() > 0.001))
  1437. {
  1438. V_spikes[i].push_back(t);
  1439. }
  1440.  
  1441. f_diff[index] = diff;
  1442. }
  1443.  
  1444. fprintf(fp_V_spikes, "\n");
  1445.  
  1446. CopyArray(f_next, f, Equations_count);
  1447.  
  1448. t += dt;
  1449.  
  1450. int percent = (int)(100 * (t - t_start) / (t_max - t_start));
  1451. if (percent != lastPercent)
  1452. {
  1453. printf("Progress: %d%%\n", percent);
  1454. lastPercent = percent;
  1455. }
  1456.  
  1457. //printf("V(24) = %f\t V_old(24) = %f\n", f[24*4], V_old(24));
  1458. //UpdateVOld();
  1459.  
  1460. //fprintf(fp_I_syn, "\n");
  1461. }
  1462.  
  1463. delete[] f_next;
  1464.  
  1465. double* V_mean_freqs = new double[Node_count];
  1466. double* V_STD_freqs = new double[Node_count];
  1467. //list<double> V_mean_freqs;
  1468.  
  1469. #pragma omp parallel for
  1470. for (int i = 0; i < Node_count; i++)
  1471. {
  1472. list<double>::iterator it_V_spikes = V_spikes[i].begin();
  1473.  
  1474. while (it_V_spikes != V_spikes[i].end() && *it_V_spikes < 0.25 * t_max)
  1475. {
  1476. V_spikes[i].pop_front();
  1477. it_V_spikes = V_spikes[i].begin();
  1478. }
  1479.  
  1480. list<double> V_freqs;
  1481.  
  1482. it_V_spikes = V_spikes[i].begin();
  1483.  
  1484. V_mean_freqs[i] = 0;
  1485. V_STD_freqs[i] = 0;
  1486.  
  1487. if (V_spikes[i].size() <= 1)
  1488. {
  1489. continue;
  1490. }
  1491. else
  1492. {
  1493. for (int j = 1; j < V_spikes[i].size(); j++)
  1494. {
  1495. double first = *it_V_spikes;
  1496. advance(it_V_spikes, 1);
  1497. double next = *it_V_spikes;
  1498.  
  1499. double T = next - first;
  1500. V_freqs.push_back(1 / T);
  1501. }
  1502. }
  1503.  
  1504. list<double>::iterator it_Freq = V_freqs.begin();
  1505.  
  1506. for (int j = 0; j < V_freqs.size(); j++)
  1507. {
  1508. V_mean_freqs[i] += *it_Freq;
  1509. advance(it_Freq, 1);
  1510. }
  1511.  
  1512. V_mean_freqs[i] /= V_freqs.size();
  1513.  
  1514. it_Freq = V_freqs.begin();
  1515.  
  1516. for (int j = 0; j < V_freqs.size(); j++)
  1517. {
  1518. V_STD_freqs[i] += pow(*it_Freq - V_mean_freqs[i], 2);
  1519. advance(it_Freq, 1);
  1520. }
  1521.  
  1522. V_STD_freqs[i] /= V_freqs.size();
  1523. V_STD_freqs[i] = sqrt(V_STD_freqs[i]);
  1524. }
  1525.  
  1526. fp0 = fopen("V_STD_freqs.txt", "w+");
  1527. for (int i = 0; i < Node_count; i++)
  1528. {
  1529. fprintf(fp0, "%f\t", V_STD_freqs[i]);
  1530. }
  1531. fclose(fp0);
  1532.  
  1533. double V_STD_mean_Freq = 0;
  1534. double V_STD_mean_Freq_not_zero = 0;
  1535. double V_STD_mean_Freq_not_zero_count = 0;
  1536.  
  1537. for (int j = 0; j < Node_count; j++)
  1538. {
  1539. if (V_STD_freqs[j] != 0)
  1540. {
  1541. V_STD_mean_Freq_not_zero_count++;
  1542. V_STD_mean_Freq += V_STD_freqs[j];
  1543. }
  1544. }
  1545.  
  1546. V_STD_mean_Freq_not_zero = V_STD_mean_Freq / V_STD_mean_Freq_not_zero_count;
  1547. V_STD_mean_Freq /= Node_count;
  1548.  
  1549. fp0 = fopen("V_STD_mean_Freq.txt", "w+");
  1550. fprintf(fp0, "%f\n", V_STD_mean_Freq);
  1551. //fprintf(fp_res, "%f\t", V_STD_mean_Freq);
  1552. fclose(fp0);
  1553.  
  1554. fp0 = fopen("V_STD_mean_Freq_not_zero.txt", "w+");
  1555. fprintf(fp0, "%f\n", V_STD_mean_Freq_not_zero);
  1556. //fprintf(fp_res, "%f\t", V_STD_mean_Freq_not_zero);
  1557. fclose(fp0);
  1558.  
  1559. double V_mean_mean_Freq = 0;
  1560. double V_mean_mean_Freq_not_zero = 0;
  1561. double V_mean_mean_Freq_not_zero_count = 0;
  1562.  
  1563. for (int j = 0; j < Node_count; j++)
  1564. {
  1565. if (V_mean_freqs[j] != 0)
  1566. {
  1567. V_mean_mean_Freq_not_zero_count++;
  1568. V_mean_mean_Freq += V_mean_freqs[j];
  1569. }
  1570. }
  1571.  
  1572. delete[] V_mean_freqs;
  1573. delete[] V_STD_freqs;
  1574.  
  1575. V_mean_mean_Freq_not_zero = V_mean_mean_Freq / V_mean_mean_Freq_not_zero_count;
  1576. V_mean_mean_Freq /= Node_count;
  1577.  
  1578. fp0 = fopen("V_mean_mean_Freq.txt", "w+");
  1579. fprintf(fp0, "%f\n", V_mean_mean_Freq);
  1580. fclose(fp0);
  1581.  
  1582. fp0 = fopen("V_mean_mean_Freq_not_zero.txt", "w+");
  1583. fprintf(fp0, "%f\n", V_mean_mean_Freq_not_zero);
  1584. fclose(fp0);
  1585.  
  1586. // CHUNKS
  1587. double chunk_t_start = 0.25 * t_max;
  1588. double chunk_t_step = 0.5;
  1589. int chunk_step_count = (0.75 * t_max / chunk_t_step);
  1590.  
  1591. double dt_chunk = 0.1 / V_mean_mean_Freq_not_zero; // 0.25
  1592. int chunks_count = chunk_t_step / dt_chunk;
  1593.  
  1594. double corr_aver_mean = 0, corr_aver_mean_not_zero = 0;
  1595.  
  1596. FILE* fp_corr_not_zero = fopen("corr_not_zero.txt", "w+");
  1597. FILE* fp_V_spikes_chunks = fopen("V_spikes_chunks.txt", "w+");
  1598.  
  1599. for (int s = 0; s < chunk_step_count; s++)
  1600. {
  1601. int** V_spikes_chunks = new int* [Node_count];
  1602.  
  1603. //#pragma omp parallel for
  1604. for (int i = 0; i < Node_count; i++)
  1605. {
  1606. V_spikes_chunks[i] = new int[chunks_count];
  1607.  
  1608. if (V_spikes[i].size() == 0)
  1609. for (int ch = 0; ch < chunks_count; ch++)
  1610. V_spikes_chunks[i][ch] = 0;
  1611.  
  1612. double ch_start = chunk_t_start + chunk_t_step * s;
  1613. double ch_end = ch_start + dt_chunk;
  1614. list<double>::iterator currentSpike = V_spikes[i].begin();
  1615.  
  1616. for (int ch = 0; ch < chunks_count; ch++)
  1617. {
  1618. while (currentSpike != V_spikes[i].end() && *currentSpike < ch_start)
  1619. currentSpike++;
  1620.  
  1621. if (currentSpike == V_spikes[i].end())
  1622. V_spikes_chunks[i][ch] = 0;
  1623. else
  1624. V_spikes_chunks[i][ch] = *currentSpike >= ch_start && *currentSpike <= ch_end;
  1625.  
  1626. ch_start += dt_chunk;
  1627. ch_end += dt_chunk;
  1628. }
  1629. }
  1630.  
  1631. char buffer[50];
  1632. //sprintf(buffer, "V_spikes_chunks_%d.txt", s);
  1633. //fp0 = fopen(buffer, "w+");
  1634.  
  1635. for (int ch = 0; ch < chunks_count; ch++)
  1636. {
  1637. for (int i = 0; i < Node_count; i++)
  1638. fprintf(fp_V_spikes_chunks, "%d\t", V_spikes_chunks[i][ch]);
  1639.  
  1640. fprintf(fp_V_spikes_chunks, "\n");
  1641. }
  1642.  
  1643. //fclose(fp0);
  1644.  
  1645. double corr = 0;
  1646. double corr_not_zero = 0;
  1647. double counter = 0;
  1648. double counter_not_zero = 0;
  1649.  
  1650. for (int i = 0; i < Node_count; i++)
  1651. {
  1652. for (int j = 0; j < Node_count; j++)
  1653. {
  1654. if (i == j)
  1655. continue;
  1656.  
  1657. int k1 = 0;
  1658. int k2 = 0;
  1659. int k3 = 0;
  1660.  
  1661. for (int l = 0; l < chunks_count; l++)
  1662. {
  1663. if (V_spikes_chunks[i][l] == 1 && V_spikes_chunks[j][l] == 1)
  1664. k1++;
  1665.  
  1666. k2 += V_spikes_chunks[i][l];
  1667. k3 += V_spikes_chunks[j][l];
  1668. }
  1669.  
  1670. if (k2 != 0 && k3 != 0)
  1671. {
  1672. corr += (double)k1 / sqrt((double)k2 * (double)k3);
  1673. counter_not_zero++;
  1674. }
  1675.  
  1676. counter++;
  1677. }
  1678. }
  1679.  
  1680. double corr_aver = corr / counter;
  1681. double corr_aver_not_zero = corr / counter_not_zero;
  1682.  
  1683. corr_aver_mean += corr_aver;
  1684. corr_aver_mean_not_zero += corr_aver_not_zero;
  1685.  
  1686. //sprintf(buffer, "corr_aver_%d.txt", s);
  1687.  
  1688. //fp0 = fopen(buffer, "w+");
  1689. //fprintf(fp0, "%f\n", corr_aver);
  1690. //fclose(fp0);
  1691.  
  1692. double ch_start = chunk_t_start + chunk_t_step * s;
  1693. fprintf(fp_corr_not_zero, "%f\t%f\n", ch_start, corr_aver_not_zero);
  1694.  
  1695. for (int i = 0; i < Node_count; i++)
  1696. delete[] V_spikes_chunks[i];
  1697.  
  1698. delete[] V_spikes_chunks;
  1699. }
  1700.  
  1701. fclose(fp_corr_not_zero);
  1702. fclose(fp_V_spikes_chunks);
  1703.  
  1704. corr_aver_mean /= chunk_step_count;
  1705. corr_aver_mean_not_zero /= chunk_step_count;
  1706.  
  1707. //fp0 = fopen("corr_aver_mean.txt", "w+");
  1708. //fprintf(fp0, "%f\n", corr_aver_mean);
  1709. //fclose(fp0);
  1710.  
  1711. //fp0 = fopen("corr_aver_mean_not_zero.txt", "w+");
  1712. //fprintf(fp0, "%f\n", corr_aver_mean_not_zero);
  1713. //fclose(fp0);
  1714. ////// CHUNKS END
  1715.  
  1716. //2
  1717. double* V_freq_sync_time = new double[Node_count];
  1718. double* V_freq_sync_time_relative = new double[Node_count];
  1719.  
  1720. #pragma omp parallel for
  1721. for (int i = 0; i < Node_count; i++)
  1722. {
  1723. list<double>::iterator it_V_spikes = V_spikes[i].begin();
  1724.  
  1725. while (it_V_spikes != V_spikes[i].end() && *it_V_spikes < 0.25 * t_max)
  1726. {
  1727. V_spikes[i].pop_front();
  1728. it_V_spikes = V_spikes[i].begin();
  1729. }
  1730.  
  1731. list<double> V_freqs_normalized;
  1732. list<double> V_freqs_time;
  1733.  
  1734. it_V_spikes = V_spikes[i].begin();
  1735.  
  1736. if (V_spikes[i].size() <= 1)
  1737. {
  1738. continue;
  1739. }
  1740. else
  1741. {
  1742. for (int j = 1; j < V_spikes[i].size(); j++)
  1743. {
  1744. double first = *it_V_spikes;
  1745. advance(it_V_spikes, 1);
  1746. double next = *it_V_spikes;
  1747.  
  1748. double T = next - first;
  1749. V_freqs_normalized.push_back(1 / T - V_mean_mean_Freq);
  1750. V_freqs_time.push_back(next);
  1751. }
  1752. }
  1753.  
  1754. list<double>::iterator it_Freq_normalized = V_freqs_normalized.begin();
  1755. list<double>::iterator it_Freq_time = V_freqs_time.begin();
  1756.  
  1757. V_freq_sync_time[i] = 0;
  1758.  
  1759. for (int j = 1; j < V_freqs_normalized.size(); j++)
  1760. {
  1761. double Freq_normalized_last = *it_Freq_normalized;
  1762. advance(it_Freq_normalized, 1);
  1763. double Freq_normalized_next = *it_Freq_normalized;
  1764.  
  1765. double Freq_time_last = *it_Freq_time;
  1766. advance(it_Freq_time, 1);
  1767. double Freq_time_next = *it_Freq_time;
  1768.  
  1769. if (abs(Freq_normalized_last) <= 0.5 && abs(Freq_normalized_next) <= 0.5 && (Freq_time_next - Freq_time_last) <= 0.035)
  1770. V_freq_sync_time[i] += Freq_time_next - Freq_time_last;
  1771. }
  1772.  
  1773. V_freq_sync_time_relative[i] = V_freq_sync_time[i] / (t_max - (0.25 * t_max));
  1774. }
  1775.  
  1776. double V_mean_freq_sync_time_relative = 0;
  1777.  
  1778. for (int j = 0; j < Node_count; j++)
  1779. {
  1780. V_mean_freq_sync_time_relative += V_freq_sync_time_relative[j];
  1781. }
  1782.  
  1783. V_mean_freq_sync_time_relative /= Node_count;
  1784.  
  1785. fp0 = fopen("V_freq_sync_time_relative.txt", "w+");
  1786. for (int i = 0; i < Node_count; i++)
  1787. {
  1788. fprintf(fp0, "%f\t", V_freq_sync_time_relative[i]);
  1789. }
  1790. fclose(fp0);
  1791.  
  1792. fp0 = fopen("V_mean_freq_sync_time_relative.txt", "w+");
  1793. fprintf(fp0, "%f\n", V_mean_freq_sync_time_relative);
  1794. //fprintf(fp_res, "%f\n", V_mean_freq_sync_time_relative);
  1795. fclose(fp0);
  1796.  
  1797. //fclose(fp_res);
  1798. fclose(fp_Max_magnitude);
  1799. fclose(fp_I_stim);
  1800. ///fclose(fp_I_syn);
  1801. fclose(fp_Ca);
  1802. //fclose(fp_IP3);
  1803. //fclose(fp_z);
  1804. fclose(fp_V);
  1805. //fclose(fp_m);
  1806. //fclose(fp_n);
  1807. //fclose(fp_h);
  1808. //fclose(fp_G_P);
  1809. fclose(fp_V_P);
  1810. //fclose(fp_m_P);
  1811. //fclose(fp_n_P);
  1812. //fclose(fp_h_P);
  1813.  
  1814. fclose(fp_V_spikes);
  1815.  
  1816. //end_rk4 = omp_get_wtime();
  1817. end_rk4 = clock();
  1818. double extime_rk4 = (double)(end_rk4 - start_rk4);// / CLOCKS_PER_SEC;
  1819. int minutes = (int)extime_rk4 / 60;
  1820. int seconds = (int)extime_rk4 % 60;
  1821. printf("\nExecution time is: %d minutes %d seconds\n ", minutes, seconds);
  1822.  
  1823. /*int nth;
  1824. #pragma omp parallel
  1825. {
  1826. #pragma omp master
  1827. nth = omp_get_num_threads();
  1828. }*/
  1829.  
  1830. fp0 = fopen("time_exec.txt", "a");
  1831. //fprintf(fp0, "%d %lf\n", nth, extime_rk4);
  1832. fprintf(fp0, "%lf\n", extime_rk4);
  1833. fclose(fp0);
  1834.  
  1835. for (int i = 0; i < Node_count; i++)
  1836. delete[] A_A[i];
  1837.  
  1838. delete[] A_A;
  1839.  
  1840. for (int i = 0; i < Node_count; i++)
  1841. delete[] B_A[i];
  1842.  
  1843. delete[] B_A;
  1844.  
  1845. delete[] C_A;
  1846.  
  1847. for (int i = 0; i < Node_count; i++)
  1848. delete[] A_N[i];
  1849.  
  1850. delete[] A_N;
  1851.  
  1852. for (int i = 0; i < Node_count; i++)
  1853. delete[] B_N[i];
  1854.  
  1855. delete[] B_N;
  1856.  
  1857. delete[] C_N;
  1858.  
  1859. delete[] A_N_P;
  1860.  
  1861. for (int i = 0; i < Node_count; i++)
  1862. delete[] B_N_P[i];
  1863.  
  1864. delete[] B_N_P;
  1865.  
  1866. delete[] C_N_P;
  1867.  
  1868. delete[] f;
  1869. delete[] f_diff;
  1870. delete[] v_4;
  1871. delete[] E_syn;
  1872. delete[] I_app;
  1873. delete[] V_spikes;
  1874. delete[] V_spikes_Freq;
  1875.  
  1876. delete[] Meander_start_from_zero;
  1877. delete[] Meander_width;
  1878. delete[] Meander_height;
  1879. delete[] Meander_interval;
  1880. delete[] last_meander_end;
  1881. delete[] V_freq_sync_time;
  1882.  
  1883. for (int i = 0; i < Node_count; i++)
  1884. delete[] tau[i];
  1885.  
  1886. delete[] tau;
  1887.  
  1888. for (int i = 0; i < Node_count; i++)
  1889. delete[] V_old_array[i];
  1890.  
  1891. delete[] V_old_array;
  1892.  
  1893. for (int i = 0; i < Equations_count; i++)
  1894. delete[] k[i];
  1895.  
  1896. delete[] k;
  1897.  
  1898. delete[] phi_k1;
  1899. delete[] phi_k2;
  1900. delete[] phi_k3;
  1901. }
  1902.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!