API tools faq
paste
SpaceQuester

Untitled

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