API tools faq
paste
SpaceQuester

Untitled

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